Data were obtained from 118 Lacaune dairy ewes reared on the Experimental Unit La Fage (INRAE UE 321 agreement A312031, Roquefort, France). The genetic structure of the INRAE La Fage flock includes independent divergent genetic lines of Lacaune dairy ewes: two selected for milk SCS, based on estimated breeding values (EBVs) for milk SCS, and the other two for PERS, based on EBVs for the coefficient of variation of milk production on the testing day. Two groups of ewes with extreme EBVs were created according to the log-transformed somatic cell count (SCC): a high-SCS line (SCS+) and a low-SCS line (SCS−). And also, two extreme groups of ewes were generated, one with high persistence (PERS+) and one with low persistence (PERS-) in the milk production curve (Table 1). Ewes were in second (n = 67) or third lactation (n = 51), ranging from 119 to 133 days in milk, and were milked twice a day (8 a.m. and 5 p.m.) with an average daily production of 1.53 ± 0.30 kg. Ewes were fed with a total of 7 kg of a mixture composed of 72% grass silage, 21% hay, and 7% barley (on a DM basis) after morning and evening milkings. They were also supplied with 100 g of barley and 100 g of a commercial protein-rich concentrate (Brebitanne®, RAGT, Albi) in the milking parlor. Ewes had free access to water and stayed indoors with no access to grazing.

Samples were collected on 2 days that were 1 week apart, following the same protocol. The animal was immobilized standing in a restraint cage, and a medical gastric tube introduced into the esophagus until reaching the rumen, the introduction depth being standardized by a graduation on the tube. The vacuum pump was turned on once the gastric tube in the rumen and turned off before removing it. The sample was placed in a cold storage box for transport to the lab, where three aliquots were collected: (1) the first aliquot of 2 ml for microbiota analysis; (2) the second aliquot of 5 ml of rumen fluid was mixed with 0.2 ml of sulfuric acid (25% v/v) for VFA analysis; and (3) the last aliquot of 40 ml for LCFA analysis. All three aliquots were frozen and stored at −20°C until analysis, except for those used for microbiota analysis, which were stored at −80°C. The gastric tube was cleaned with hot water after each collection.

Rumen samples were collected before milking from each animal. For both days of sampling, a first group of 60 animals were collected from 8:30 to 11:30 a.m. and a second group of 60 animals were collected from 1:30 to 4:30 p.m. The same ewes were always sampled in the morning or the evening. To avoid dilution of samples by water, the ewes were denied access to water 2 h before sampling. Animals did not have access to feed from the previous evening (10:00 p.m.) for the group sampled in the morning, and from the early morning (7:00 a.m.) for the group sampled in the afternoon. Consequently, the number of hours of fasting varied from 6.5 to 13.5 depending on the order and moment of sampling.

For each of the 2 weeks, milk yields were recorded on the previous evening and the morning of the sampling day, and milk samples were collected and preserved with bronopol (Agrolab, Aurillac, France) for milk composition analysis or by immediate freezing for LCFA analysis. Somatic cell count and fat and protein contents analyses were performed on both the evening and the morning milk samples. LCFA analysis was performed on the morning milk samples only.

Total DNA from 80 μl of rumen samples was extracted and purified using the QIAamp DNA Stool Mini Kit (Qiagen Ltd., West Sussex, United Kingdom) according to the manufacturer’s instructions, with a bead-beating step in a FastPrep instrument (MP Biomedicals, Illkirch, France). The 16S rRNA V3–V4 regions gene of the extracted DNA strands were amplified (first PCR: 30 cycles) with the primers F343 (5′-CTTT CCCTACACGACGCTCTTCCGATCTACGGRAGGCAGCAG-3′; Liu et al., 2007) and R784 (5′-GGAGTTCAGACGTGTGCTCTTCCGATC TTACCAGGGTATCTAATCCT-3; Andersson et al., 2008). As Illumina MiSeq technology enables 250 bp reads, the ends of each read were overlapped and stitched together to generate full-length reads of the entire V3 and V4 regions in a single run. Single multiplexing was performed using a 6 bp index, which was added to R784 during a second round of PCR with 12 cycles using the forward primer (AATGATACGGCGACCACCGAGATCTACACTCT TTCCCTACACGAC) and reverse primer (CAAGCAGAAGACGGCATACGAGATGTGACT GGAGTTCAGACGTGT). The PCR products were purified and loaded onto an Illumina MiSeq cartridge (Illumina, San Diego, CA, United States) at the Genomic and Transcriptomic Platform (INRAE, Toulouse, France) according to the manufacturer’s instructions.

Sequence reads were demultiplexed, and each paired-end read was assigned to its sample based on the previously integrated index, and processed with the FROGS 3.0 pipeline (Escudié et al., 2018). The procedure consisted of the following steps: (1) read pre-processing, removing sequences with primer mismatch, displaying unexpected length (<300 or >500 bp), or with ambiguous bases; (2) sequence clustering with denoising and one sequence difference between each of the three aggregation steps of clustering; (3) chimera removal; (4) cluster filtering with Bokulich filter (removing clusters with abundances <0.005%; Bokulich et al., 2013); and (5) taxonomy assignment to operational taxonomic units (OTUs) using the SILVA 138.16S pintail 100 database. From this process, an abundance table containing the number of sequences per OTU and rumen sample was obtained.

Two gas chromatography analyses were performed at the National Veterinary School of Toulouse (Toulouse, France) with rumen samples: one for VFAs and the other for LCFAs.

Six VFAs, acetic acid (C2:0), propionic acid (C3:0), butyric acid (C4:0), valeric acid (C5:0), isobutyric acid (iso-C4:0), and isovaleric acid (iso-C5:0), were quantified using automated gas separation, according to the method of Playne (1985) and modified as follows. The rumen samples were first centrifuged at 2,880 × g for 20 min to separate the liquid phase. For protein removal, 1 ml of supernatant was mixed with 200 μl of (25% v/v) metaphosphoric acid and further centrifuged at 20,000 × g for 15 min. Then, 100 μl of the supernatant was added to 75 μl (0.2% v/v) of 4-methylvaleric acid as an internal standard and 900 μl of ultrapure water. From this mixture, 1 μl was then injected into a gas chromatograph (Hewlett Packard, Model 7890A) equipped with a DB-FFAP column (30 m × 0.53 mm i.d., 1-μm film thickness, Agilent Technologies, Palo Alto, CA, United States) and an FID detector (Avondale, PA, United States). Chromatograms were integrated using Chromeleon software (Thermo Fisher Scientific, version 6.8, Waltham, MA, United States). The sum of the six VFA concentrations was defined as the total concentration and was used to obtain the molar proportions of each VFA.

The LCFAs of rumen content were extracted and methylated in situ using the procedure described by Park and Goins (1994), except that the solution of 14% boron trifluoride in methanol was replaced by a solution of methanol–acetylchloride (10:1). Nonadecanoic acid (C19:0) was used as the internal standard at a dose of 0.8 mg. The fatty acid methyl esters (FAMEs) were then quantified by gas chromatography (Agilent 6890N, Network GC System, equipped with a model 7,683 auto injector, Agilent Technologies, Palo Alto, CA, United States) using a fused silica capillary column (100 m × 0.25 mm i.d., 0.20 μm film thickness, CPSil 88, Varian, Middelburg, the Netherlands) as described by Zened et al. (2011). Peaks were identified and quantified by comparison with commercial standards (Sigma Co., St Louis, MO, United States), except for C18:1, C18:1 trans-9, C18:1 trans-11, and C18:1 cis-9, which were identified by the order of elution. Chromatograms were integrated using the Peak Simple software (Peak Simple Data System, version 2.83, SRI, Torrance, CA, United States). Results are expressed as the percentage of total FAME. The 29 measured fatty acids (FAs) were: C12:0, C13:0, anteiso-C13:0, iso-C13:0, C14:0, iso-C14:0, C15:0, anteiso-C15:0, iso-C15:0, C16:0, C17:0, anteiso-C17:0, iso-C17:0, C18:0, C18:1 cis-9, a mix of C18:1 cis-11 and C18:1 trans-15, C18:1 cis-12, C18:1 cis-15, a mix of C18:1 trans-6, trans-7, trans-8, C18:1 trans-9, C18:1 trans-10, C18:1 trans-11, C18:1 trans-12, C18:1 trans-16, C18:2, C18:2 cis-9, trans-11, C18:2 trans-11, cis-15, C18:3, and C20:1.

The four milk samples from each ewe were analyzed at Agrolab (Aurillac, France). Fat and protein contents were obtained using mid-infrared spectrometry (Milk-Scan™ FT 6000 instrument, Foss, Nanterre, France) and somatic cell count was quantified using a Fossomatic cell counter (Foss, Nanterre, France), to which a log-transformation was applied to obtain the SCS (Ali and Shook, 1980). The data were translated into daily variables. Daily milk yield was obtained by summing the morning and evening milk yields. Daily protein content, fat content, and SCS were computed as the average of the morning and evening values weighted by the corresponding milk yields.

Long-chain fatty acid percentages in total FAME were measured in the morning milk for the 2 weeks by gas chromatography at the National Veterinary School of Toulouse (Toulouse, France), following the same method as for the rumen samples. The 38 measured FAs included the 29 measured in the rumen in addition to C4:0, C6:0, C7:0, C8:0, C9:0, C10:0, C11:0, C14:1, and C16:1.

All analyses were performed using R software (R Core Team, 2021)¹. Statistical significance was set at p < 0.05.

Microbiota and fatty acids data are compositional data (Gloor et al., 2017), meaning that the information is contained in the ratios between variables, not in the values themselves, due to the restrictions imposed by either the sequencing technology or the measurement unit (percentages). Thus, we applied the compositional data approach proposed by Aitchison (1986) for the composition of OTUs and FAs. It consisted of imputing zero values in the data set to transform the counts into log-ratios using the centered log-ratio (CLR). Then, the multivariate analysis described in a previous study by Martinez Boggio et al. (2021) was applied. Two approaches were applied to impute the zero values, considering their nature. In the bacterial abundance table, as zeros referred to a probability of count, geometric Bayesian multiplicative replacement was applied (cmultRepl function from the zCompositions package); in the LCFA datasets, the zero values reflected the detection limit of the chromatograph; thus, the expectation–maximization procedure was applied (lrEM function from the zCompositions package). Then, the CLR transformation was applied to OTUs, VFA, and LCFA data (clr function from the compositions package).

Table 2 presents the models used to correct the data for fixed effects and repeated measures (random animal effect) and the number of animals considered. Ewes were removed from OTUs models when one of their samples had a low sequencing quality (<500 OTUs). Ewes missing one of the two samples were removed from rumen and milk LCFA models. Ewes presenting values of VFAs under the detection limit of the gas chromatograph were removed. The fixed effects included in the models were: genetic lines (four levels: PERS−, PERS+, SCS-, and SCS+), week of sampling (two levels: “week 1” and “week 2”), parity (two levels), SOCS2 genotype (two levels: “TT/CT” or “CC”), time after feeding (six levels of equal size), and number of sequences per sample (four levels of equal size).

The significance of the fixed effects was assessed using the ANOVA function from the sasLM package. For milk yield, fat and protein contents, SCS, and the six VFAs, the fixed effects were included in the model when they were significant. For OTUs, rumen, and milk LCFAs, fixed effects were retained in the models when they were significant for at least 10% of the variables. Linear mixed models were defined for each trait (one OTU being one trait), using the lmer function from the lmerTest package, and the corresponding variances of random animal and residual effects were obtained (VarCorr function from the lme4 package), allowing to compute the repeatability, defined as the animal variance divided by the total variance. Then, Spearman correlations were performed between the OTU repeatability values and their percentage of zeros or average abundance.

To perform multivariate analyses assessing the effect of the week, residuals of OTUs and FAs obtained from ANOVA were used without correcting for the week effect (Table 2). Sparse partial least square discriminant analyses (sPLS-DA) and sparse partial least square analyses (sPLS) were performed, to assess the influence of the week on bacterial composition and identify relationships between the OTUs and the FAs, respectively (spls and splsda functions from the mixOmics package). The sparse procedure allowed to reduce the number of components, i.e., the dimensionality of the analyses, and select only the most relevant variables, i.e., the OTUs, for each component. The number of components was previously chosen to explain 90% of the variance, using a principal component analysis. The number of variables retained by component was determined using the CLR lasso penalty method considering the penalization obtained with a 10-fold cross-validation (cv.glmnet and glmnet functions from the glmnet package).

The predictive ability of the sPLS-DA model was evaluated through the overall misclassification error rate after a 5-fold cross-validation repeated 10 times (perf function). For each sPLS, the 20 OTUs presenting the highest association with the FAs were selected and Pearson correlations were computed. Only the FAs common in both weeks and having at least one significant correlation with the 20 selected OTUs were presented. Fisher exact tests were performed to estimate over- or under-representation of specific phyla or highly repeatable OTUs among those selected for sPLS-DA and sPLS compared to all OTUs.

A total of 2,500,763 sequences were obtained, with an average of 10,736 ± 3,745 reads per sample, grouped into 2,079 OTUs, with an average of 1,406 ± 211 OTUs per sample. Ten phyla were identified; the three major phyla, based on the percentage of sequences assigned to each phylum out of the total number of DNA sequences, were Bacteroidetes (53.5%), Firmicutes (35.4%), and Proteobacteria (6.2%). The other phyla were Fibrobacterota (2.7%), Spirochaetota (1.2%), Actinobacteria (0.5%), Patescibacteria (0.3%), and Desulfabacterota, Elusmicrobiota, and Campylobacterota, each lower than 0.2%. We found that 81.2% of the DNA sequences enabled genus-level classification, resulting in the identification of 117 genera. The most important ones were Prevotella (30.78%), Lachnospiraceae_NK3A20_group (6.3%), Ruminococcus (5.3%), Christensenellaceae_R_7_group (4.7%), Sphingomonas (4.5%), and Rikenellaceae_RC9_gut_group (3.9%). The lowest taxonomic rank considered in this study was genus, as 92.7% of the reads were assigned to unknown species or were multi-affiliated species as expected with short 16S reads.

Operational taxonomic unit data were corrected by six effects that were significant for at least 10% of the OTUs: genetic lines (25.9%), week of sampling (16.2%), time after feeding (13.9%), number of sequences per sample (12.8%), parity (11.4%), and the interaction between the genetic line and parity (11.9%), with percentages in brackets indicating the percentage of OTUs for which the given factor was significant. The explained variance was estimated for all OTUs, resulting in an average of 0.12 ± 0.04, ranging from 0.03 to 0.50 according to OTUs.

Out of the 2,079 OTUs, 1,665 showed an estimated repeatability ranging from 0 to 0.93 with a median of 0.15 ± 0.14 (Figure 1). The remaining 414 OTUs had no estimate; linear models did not converge because of the small number of samples. Of the 1,665 OTUs, 66 presented a repeatability of over 0.50, corresponding to a frequency of 4%. Spearman correlations were performed between the repeatability values of the OTUs and their percentage of zeros (R = −0.01, p = 0.50) and between the repeatability of the OTUs and their average abundance (R = 0.24, p < 0.001).

The repeatabilities ranged from 0.21 to 0.57 for VFAs, from 0.05 to 0.50 for rumen LCFAs, and from 0.02 to 0.86 for milk LCFAs (Figure 1). The repeatability was 0.83 for milk yield, 0.74 for fat content, 0.77 for protein content, and 0.74 for SCS. The models did not converge for the LCFAs C18:1 trans-12 and C18:2 cis-9, trans-11 in the rumen and C20:1 in the milk.

The model of the sPLS-DA on bacterial composition over a week included 160 components and 57 variables per component. Based on OTU residual abundances, rumen samples were discriminated by the week of collection (Figure 2). The overall error rate of the model was 0.21. In the first component (Table 3), week 2 was characterized by a lower number of OTUs belonging to the phylum Firmicutes (p < 0.001) and a higher number of OTUs belonging to the phylum Bacteroidetes (p < 0.01) than when considering all OTUs. In the second component (Table 3), week 1 had a higher presence of OTUs from the phyla Firmicutes (p < 0.01) and a lower presence of OTUs from the phyla Bacteroidetes (p = 0.01) than when considering all OTUs.

In the first component, 49 OTUs among the 57 selected to discriminate week 1 and week 2 had a repeatability of less than 0.50 (Supplementary Table S1), three had a repeatability of over 0.50, and the remaining five OTUs did not have an estimated repeatability. In the second component (Supplementary Table S2), 44 OTUs of 57 had a repeatability of less than 0.50, one had a repeatability of over 0.50, and the remaining 12 OTUs had no estimated repeatability. The proportion of OTUs presenting repeatability values of less than 0.50 among the 57 selected for each component of the sPLS-DA was not significantly different from that considering all OTUs (p = 0.46 and p = 0.47 for components 1 and 2, respectively).

Bacteroidetes, Firmicutes, and Proteobacteria were the three main phyla observed in the rumen of adult dairy ewes fed a mixed diet of forages and concentrates, in accordance with that observed in sheep (Belanche et al., 2019; Liu et al., 2020) and cows (Bainbridge et al., 2016; Plaizier et al., 2017). However, the rumen bacterial composition is affected by various factors, mainly biological ones. The most known is the diet composition, often represented by the forage/concentrate ratio that can influence family and genus abundances (Henderson et al., 2015). Another one is the feeding time, influencing diurnal variations in bacterial concentrations in sheep fed once daily (Warner, 1966), with the abundance of families belonging to the phylum Firmicutes increasing with time after feeding (de Assis Lage et al., 2020). In addition, both parity and lactation stage affect bacterial composition by modifying the abundance of phyla and genera (Pitta et al., 2014; Bainbridge et al., 2016; Xue et al., 2018). Although many authors have found some bacterial taxa to be common in all the animals in their study (Jami and Mizrahi, 2012; Xue et al., 2018; Huang et al., 2021), interindividual variations in bacterial composition were observed, and evidence of genetic determinism of abundance of some bacteria has been raised by Sasson et al. (2017) and Difford et al. (2018).

Technical and bioinformatics processing factors are also known to greatly influence the observed bacterial abundance (Pollock et al., 2018; Wang and LêCao, 2020). Diversity and relative abundances vary depending on the type of sampling (de Assis Lage et al., 2020) and DNA extraction methods (Henderson et al., 2013; Gerasimidis et al., 2016). Considering sequencing strategies, the choice of primers influences the estimation of bacterial abundances (Tremblay et al., 2015; Fouhy, 2016) while the sequencing platform influences the read length and error rate (D’Amore et al., 2016; Kchouk et al., 2017). The number of observed OTUs depends on the chosen pipeline, in particular, the clustering and filtering steps (Schloss, 2010). As the sequencing depth is highly variable among the samples, at random or depending on the sequencing run, normalization is usually applied to the data (McMurdie and Holmes, 2014; Weiss et al., 2017).

As the current objective was to determine if the bacterial composition was stable over time, the week was the factor of greatest interest in this study. To allow proper observation, breed, diet, and lactation stages were experimentally controlled. Moreover, the same DNA extraction method was used for the samples of both weeks, only one sequencing run was performed, and the bioinformatic analysis included all samples. The other effects, including genetic lines, SOCS2 genotype, parity, time after feeding, and sequencing depth were corrected through linear models following the methodology used by Martinez Boggio et al. (2021).

Assessing the stability of the bacterial composition over a short period of time will support minimizing repeated sampling, as they are invasive for the animals and time-consuming. It is essential to know if the results obtained at one sampling are reliable, whether the conclusions are the same 1 week later. In this study, two criteria were used to assess the stability of the bacterial composition. In the first criterion, bacterial composition was defined as stable when phyla, genera, and OTU proportions were conserved from 1 week to the next for the same animal. In the literature, various methods have been used for assessing bacterial composition stability, such as two one-sided tests (Skarlupka et al., 2019), principal coordinate analysis, and analysis of similarity (Huang et al., 2020; Mamun et al., 2020). In the current study, considering the compositional nature of the data, we used multivariate analyses (working in the Euclidian space), such as sPLS-DA and sPLS. In the second criterion, at the OTU level, the repeatability of the OTUs was used to determine their stability, ranging from 0 (completely unstable) to 1 (completely stable), with OTUs having a repeatability greater than 0.50 being defined as stable. Another method, the Lin’s concordance correlation coefficient, which is also conceptually similar, was used by Zhu et al. (2021) to estimate OTU stability.

Although the week effect was the second largest in the ANOVA, it was significant for only 16% of OTUs. A clustering of the samples by week was observed in the sPLS-DA, mainly characterized by different proportions of Bacteroidetes and Firmicutes, but less than 6% of the variance was explained. Even if the results of these two analyses suggested bacterial composition instability, the high unexplained variance made them inconclusive. As explained previously, various factors influencing bacterial composition cannot be controlled or corrected for. In our study, such factors impacted the variance explained in the ANOVA and the sPLS-DA analyses. They could be unrecorded technical effects such as (1) contamination from the oral cavity or saliva (Lodge-Ivey et al., 2009; Terré et al., 2012), (2) multiple relocation of the gastric tube in the rumen (Geishauser and Gitzel, 1996), and (3) variable proportions of solid and liquid phases (Henderson et al., 2013; Vaidya et al., 2018); or biological effects such as the stress induced by handling before and during sampling (Yoshihara and Ogawa, 2021).

However, the repeatability of 96% of the OTUs was under 0.5, and it was this large proportion of non-stable OTUs (median repeatability of 0.15) that allowed us to state bacterial composition instability. Similar results were obtained by Jewell et al. (2015); Bainbridge et al. (2016); and Zhu et al. (2021) when working over long periods of 75–122 days, which was expected because of the numerous effects that could alter bacterial composition. The 96% of non-stable OTUs (96%) obtained is consistent with the 97% found by Zhu et al. (2021). Contrasting with our results, some authors obtained stability over short periods, from 3 days to 1 week using analyses comparable to ANOVA and sPLS-DA (Skarlupka et al., 2019; Huang et al., 2020; Mamun et al., 2020). However, they did not consider the compositional nature of the data neither computed the OTU repeatability that demonstrated that most of them were moderately repeatable (<0.50). Moreover, the number of animals used in this study was larger than in the previous ones, allowing for revealing small differences between weeks covered up by uncontrolled alterations of the bacterial composition. Those three specificities of our study allowed to clearly reveal bacterial composition instability over short periods.

The biological links between OTUs and FA phenotypes were compared between week 1 and week 2 to discuss functional ruminal bacterial stability over time. At the OTU level, only three or four OTUs were conserved from week 1 to week 2, with similar correlations to FA phenotypes (Figures 3–5). However, we found that OTUs belonging to a same family presented all correlations of the same sign for a given FA. This phenomenon was observed for some of the most predominant families in the rumen, namely Lachnospiraceae, Rikenellaceae, and Ruminococcaceae, and may suggest functional redundancy (Wohl et al., 2004; Weimer, 2015; Louca et al., 2018). That is to say, different OTUs belonging to the same family may have the same function, allowing them to maintain the organism function despite a variable microbial composition. In addition, OTUs from Prevotellaceae family presented both positive and negative correlations with the FAs. The main hypothesis to explain this observation is the large genetic and functional diversity of this family (Stewart et al., 1997; Matsui et al., 2000). Further investigation of these aspects is not possible with 16S rRNA gene sequencing, because as stated by Plummer and Twin (2015), it does not allow access to the species and strain classification necessary to investigate the function of the bacteria identified.

The OTUs generally presented low repeatability, with only 4% of them defined as stable. It is notable that stable OTUs were overrepresented among the OTUs highly correlated with FA phenotypes in both weeks, ranging from 3 to 11 among 20. If the stable OTUs are considered as those having high genetic determinism, it was not surprising to find them highly associated with FA phenotypes, for which a significant part of the variability is related to animal genetics. The correlations between bacteria and FA phenotypes from week 1 to week 2 were not conserved, but OTUs related to FA phenotypes appeared to be the most repeatable.

Week-to-week variation in the rumen bacterial composition of Lacaune dairy ewes was observed, with only 4% of the OTUs being stable, which seemed to alter the correlations between microbiota and FAs phenotypes, as only a few OTUs were conserved between 1 week and the next. Even if FA phenotypes were not linked to the same OTUs from 1 week to the next, it is noticeable that they were particularly associated with stable OTUs, possibly due to common genetic determinism of those traits and OTUs. In conclusion, we proved with a large dataset of 188 dairy ewes that bacterial composition and its phenotypic correlations with fatty acids are not transposable from 1 week to the next. Further work is necessary to confirm those results, as very few studies were performed on microbiota stability over short periods of time, and none using compositional data approach. We hope that our study will be a first step in identifying microbiota members with high repeatability and thus potential heritability, to then consider the possibility of genetic selection.