Introduction

The fecal microbiome consists of a complex community of microorganisms and represents a central issue in relation to cattle welfare and feed efficiency. In particular, the associations between fecal microbiome and animal health have been shown in the intestinal microbiota of calves (Oikonomou et al. 2013). The main factor that influences fecal microbiome composition is animal diet. Callaway et al. (2010) carried out an evaluation of bacterial diversity of 6 cattle (3 Jersey cows and 3 Angus steers) through a comparison of 3 different diets in terms of amount of dried distillers grain; Shanks et al. (2011) analyzed the structure of fecal community in 30 adult beef cattle equally divided in 3 diet groups; and Rice et al. (2012) evaluated the influence of different types and amount of distillers grains on fecal microbial assemblages in 20 crossbreed cattle. The forage to concentrates ratio in the diet is the major factor affecting fecal microbiome composition in cattle (Kim et al. 2014). According to the meta-analysis of Kim and Wells (2016), the fecal cattle microbiome is composed of 10 phyla, 17 classes, 28 orders, 59 families, and 100 genera. Firmicutes is the most represented phylum, followed by Bacteroidetes and Proteobacteria. Within Firmicutes, Clostridia and Fecalibacterium are the largest class and genus, respectively. Within Bacteroidetes, Bacteroidia is the largest class and Prevotella the largest genus. Finally, Proteobacteria includes Gammaproteobacteria and Succinivibrio as the most abundant class and genus, respectively. To investigate the microbiome in cattle, most studies have used DNA-based methodologies such as Sanger sequencing technology, quantitative real-time polymerase chain reaction, and phylogenetic microarrays (Kim et al. 2017; Mende et al. 2012). Currently, the next-generation sequencing (NGS) is considered the most reliable approach to evaluate the diversity of bacteria, both in rumen and feces of cattle (Kim et al. 2017). The 16S rRNA is widely used as reference gene to determine the composition of bacterial community due to its phylogenetic variability (Tringe and Hugenholtz 2008); indeed, it includes 9 hypervariable regions and 10 conserved regions. The conserved regions (C1 to C10) are shared among bacterial and archaeal species, whereas 16S rRNA hypervariable regions (V1 to V9) are different. The latter can be targeted to identify individual bacterial or archaeal species using PCR with species-specific primers for the 16S rRNA gene. Data analysis assigns 16S rRNA sequences to operational taxonomic units (OTUs) that can be identified according to available database. The literature reports differences of fecal microbiome composition within beef cattle breeds, across dairy and beef cattle breeds, and within crossbreed cattle (Durso et al. 2012). Comparisons between heifers and cows of dairy and dual-purpose cattle breeds are currently lacking. The aim of this study was to characterize and analyze the difference of the fecal microbiome community of heifers and cows of dairy and dual-purpose cattle breeds, targeting the hypervariable regions of the bacterial 16S; moreover, we evaluated if the microbiome composition is conserved between the breeds that underwent to different management and diet composition. Gaining knowledge on these aspects is expected to be beneficial to investigate changes in methane emissions and variation of feed efficiency, as well as to develop non-invasive routine controls to evaluate animal welfare and health.

Materials and methods

Sample collection and DNA extraction

Fecal samples were collected through rectal picking during routine health monitoring of animals by authorized veterinarians of the Breeders Association of Veneto Region (Italy). Twenty individual samples (one sample per animal) from 2 single-breed herds (one rearing Holstein-Friesian and the other Simmental breed) were collected for microbiome analysis, considering two categories: cows and heifers. Animals were divided in 2 experimental groups: (1) Holstein-Friesian heifers (HFH, n = 5) and Holstein-Friesian cows (HFC, n = 5); (2) Simmental heifers (SIH, n = 5) and Simmental cows (SIC, n = 5). The sample size was chosen after a literature review (Callaway et al. 2010; Sandri et al. 2018). A description of rations used in the 2 farms is presented in Table 1.

Table 1 Diet composition (kg/day) in each farm and cattle category
Full size table

Feces were stored at −80 °C within 1 h from sampling. DNA extraction was performed through AllPrep PowerFecal® DNA/RNA Kit (Qiagen, Hilden, Germany), and the quantity and quality of total DNA were checked through spectrophotometer assay (Multiscan Sky, Thermo Fisher Scientific, USA).

Next-generation sequencing

Total genomic DNA was amplified by using a standard protocol and modified primers (Takahashi et al. 2014). Amplicons were purified through magnetic beads Agencourt XP 0.8× (Beckman Coulter, Brea, CA, USA) and amplified through HiSeq by using Index Nextera XT kit (Illumina, San Diego, CA, USA). All amplified sequences were normalized by SequalPrep (Thermo Fisher, Waltham, MA, USA) and precipitated through magnetic beads Agencourt XP 0.8×. Libraries were loaded onto MiSeq (Illumina, San Diego, CA, USA) and sequenced following V3-300PE strategy.

Statistical analyses

The OTUs obtained from 16S rRNA sequencing results were filtered for 0.005% frequency, and organized in an OTU table. The taxonomic survey was obtained from a cross comparison between the QIIME2 software package (http://qiime.org/) and the two databases SILVA v.1.132 and Geengenes v.13.8 (the last as a comparison); clustered OTUs were matched against references from databases. Alpha-diversity analyses were conducted considering Observed OTUs, Shannon Index, Pielou’s Evenness, and Faith’s Phylogenetic Diversity Index using QIIME2 platform. All sequences were clustered with representative OTUs and cleaned considering 97% of identity as cutoff. The statistical significance of each index was analyzed by Kruskal-Wallis non-parametric test, comparing cows and heifers within the same breed. Beta-diversity was calculated through the Bray-Curtis Metric, Jaccard Metric, and the UniFrac Metric (weighted and unweighted) to evaluate the dissimilarity and distance between the animals of the same breed. Dissimilarities in fecal bacteria and archaea communities were visualized using principal coordinate analysis (PCoA) method. The permutational analysis of variance (PERMANOVA) (Anderson 2005) and analysis of similarities (ANOSIM) were performed in R-vegan package adonis (Oksanen et al. 2017). Finally, differential abundance test was performed using ANCOM packages of R software (Team R, Core Team 2015). Significance was determined through W statistic, which indicates the number of times the null hypothesis was rejected. Positive values of W statistic correspond to more abundant taxa in the comparison of HFH vs HFC, and SIH vs SIC. The test was performed accounting for the percentile abundance.

Results

Taxonomic identification

The total OTUs obtained (2,302) were clustered trough SILVA and Geengenes for the taxonomic analysis; this identified the presence of 2 kingdoms, 14 phyla, 22 classes, 34 orders, 74 families, 212 genera, and 350 species, while the remaining sequences were not assigned to known phyla (Table 2).

Table 2 Average number of sequences (standard deviation) per taxon obtained from fecal samples of Holstein-Friesian heifers (HFH), Holstein-Friesian cows (HFC), Simmental heifers (SIH), and Simmental cows (SIC)
Full size table

Archaea were represented by Euryarchaeota phylum which includes 5 genera equally distributed among the 2 experimental groups (Table 2). Within this phylum, we found several microorganisms that colonize the rumen and are involved in methane production (Holmes and Smith 2016). It is worth noting that the difference regarding Methanosphaera genus was larger for HFC (141 sequences) and SIC (157 sequences) than HFH (22 sequences) and SIH (74 sequences), likely due to the physiological status and diet composition (Hook et al. 2010). As expected, Bacteria was the largest domain in all experimental groups, representing about two-thirds of the total microbiome (Fig. 1).

Fig. 1
figure 1

Next-generation sequencing (NGS) relative abundance at the kingdom, phylum, and class levels. The figure represents 5 experimental replicates for each experimental group (x-axis). Annotation was done using SILVA database: kingdoms are shown on the left side of the figure, phyla on the central part, and classes on the right side

Full size image

Abundance of bacterial and archeal communities differs between heifers and lactating cows

As expected, within bacteria domain, Firmicutes represented the most abundant phylum (Table 2 and Fig. 1). In HFH, Paeniclostridium and Romboutsia were the largest genera, covering 1,543 and 1,119 sequences, respectively, followed by Clostridium, Eubacterium, and Turicibacter. The same order of abundance was maintained in SIH and SIC; however, in HFC, Anaerovibrio (97 sequences), Blautia (137 sequences), Marvinbryantia (157 sequences), Oscillibacter (148 sequences), Roseburia (202 sequences), and Lachnospiraceae AC2044group (274 sequences) showed high abundance. In SIH, the Candidatus class showed good number of sequences (138 sequences), while in SIC Syntrophococcus sequences were greater than the other groups. Moreover, Fourneriella genus was identified only in HFC (13 sequences) and SIC (12 sequences), being probably related to a diet rich in corn. Intestinimonas was represented only in HFH (3 sequences) and SIH (6 sequences).

Bacteroidetes represented the second most abundant phylum with 12,529 sequences in HFH; 10,308 in HFC; 9,997 in SIH; and 10,019 in SIC (Table 2). Prevotellaceae was the largest family and in all groups comprised more than 3,000 sequences, with a particular relative abundance of some genera in heifers, as Prevotellaceae UCG-004, or in cows, as Prevotellaceae UCG-003. Rikenellaceae RC9 gut group was the second largest genus and comprised 2,380 sequences in HFH; 2,015 in HFC; 2,238 in SIH; and 1,918 in SIC. The third largest genus identified was Bacteroides with high differences between heifers and cows; indeed, in HFH (1,715 sequences) and SIH (1,759 sequences), the abundance was lower compared with HFC (2,807 sequences) and SIC (2,350 sequences). Alistipes was another important genus (Table 2) that showed great variability among the experimental groups: 1,489 sequences in HFH; 866 in HFC; 1,108 in SIH; and 709 in SIC, and it is typical component of gut microbiome (Xin et al. 2019).

The third largest phylum was Actinobacteria; these Gram-positive bacteria have been studied mainly as a minor component of rumen microflora, representing about 3% of entire community, and this aspect reflect our results (Sulak et al. 2012). This group is involved in the amylase, caseins, gelatinase, lipase, chitins, and cellulose enzyme production (Borsanelli et al. 2018).

Minor phyla

The incidence of Proteobacteria phylum was variable among the 2 experimental groups. In particular, 251 sequences were detected in HFH, 760 in HFC, 566 in SIH, and 216 in SIC (Table 2). Succinivibrio was the largest genus of Proteobacteria especially in HFC (512 sequences); the high presence of genus Succinivibrio in cows is probably related to an abundant corn-based diet as described by Kim et al. (2014). Mailhella genus was the largest phylotype in heifers (178 sequences in HFH and 279 in SIH); however, their role is almost unknown in cattle and probably their different abundance could be related to animal diet. It is worth noting that the genera of Proteobacteria phylum reported in our study were significantly different from those reported in the study of Kim and Wells (2016). Chloroflexi, Cyanobacteria, Elusimicrobia, Epsilonbacteraeota, Fibrobacteres, and Verrucomicrobia comprised the classes that showed differences in number of sequences among the experimental groups (Table 2). Indeed, their abundance was higher in heifers than cows. Chloroflexi phylum was observed only in HFH (3 sequences) and SIH (34 sequences) due to contamination, since this phylum includes environmental photosynthetic bacteria (Borsanelli et al. 2018). Elusimicrobia and Epsilonbacteraeota are 2 phyla not commonly present in cattle feces (Kim et al. 2017), while Cyanobacteria have been reported in other species, but their role in cattle remains still unknown (Shepherd et al. 2012). Fibrobacteres phylum belongs to the group of bacteria that colonize the rumen; their function in cattle is related to fiber digestion, and thus, their presence is associated to a diet rich in forage (AlZahal et al. 2017). Verrucomicrobia phylum was the most variable among the experimental groups. Indeed, we found 481 sequences in HFH, 82 in HFC, 850 in SIH, and 379 in SIC (Table 2). These bacteria are usually not affected by the different diets or sample fractions and their role is still unknown (Deusch et al. 2017). Patescibacteria and Spirochaetes are commonly present in feces of cattle (Nyonyo et al. 2014); indeed, their level has been considered constant among the experimental groups (Table 2). Patescibacteria function is still unknown, while Spirochaetes phylum is associated with the cellulolytic activity (Nyonyo et al. 2014).

Alpha-diversity

Results of α-diversity test are shown in Table 3. The number of OTUs for all experimental groups ranged from 676 to 798; in particular, the average number of identified OTUs ± SD in HFH, HFC, SIH, and SIC was 798 ± 144, 676 ± 75, 733 ± 88, and 716 ± 75, respectively. The number of observed OTUs was not statistically different among samples of the same breed, whereas significant differences were observed for the other indexes. Pielou’s Evenness Index showed differences (p < 0.05) between cows and heifers of the same breed: HFH (0.915 ± 0.003) vs HFC (0.883 ± 0.016), and SIH (0.918 ± 0.008) vs SIC (0.899 ± 0.007), related to a different abundance of species, due to the diet composition and animal category. Through Shannon Index the observed number of species was different only between HFH (32.970 ± 2.997) and HFC (28.557 ± 2.228) (p < 0.05). Faith’s Phylogenetic Diversity Index underlined that just the comparison between HFH (8.799 ± 0.241) and HFC (8.299 ± 0.234) was statistically significant (p < 0.05), proving that within Holstein-Friesian breed, the microbial species seem to have different phylogenetic taxon. Finally, the OTU table was normalized according to the OTU abundance across samples, and the results were clustered as the heatmap to provide a better pattern across the experimental groups (Fig. 2).

Table 3 Summary of 16S rRNA OTU data and α-diversity values (Evenness Index, Shannon Index, and Faith’s Phylogenetic Diversity Index) in association with statistical test
Full size table
Fig. 2
figure 2

Hierarchical clustering of core OTUs in the 4 experimental groups, divided by quadrants. All core OTUs were clustered among the animals of the experimental groups. The blue background represents 0 counts, whereas red color indicates higher counts for that particular taxonomic unit in a specific sample

Full size image

Beta-diversity

Results of β-diversity were computed on the rarefied OTU table using different metrics (Table 4). The Bray-Curtis Metric showed high dissimilarity across experimental groups of the same breed, where the index was weakly high in the comparison HFH vs HFC (0.74), and SIH vs SIC (0.64). The PERMANOVA test underlined that the difference was statistically different between the samples of the same breed (HFH vs HFC and SIH vs SIC; p < 0.05). The Jaccard Metric provides the difference in the microbial composition among the experimental groups, HFH vs HFC (0.80) and SIH vs SIC (0.78), underlying a significant difference in fecal microbial composition (p < 0.05). The UniFrac Metric was considered unweighted and weighted (Table 4). The UniFrac Metric was low in the comparisons HFH vs HFC (unweighted 0.51 and weighted 0.23) and SIH vs SIC (unweighted 0.48 and weighted 0.21). Such results were supported by the statistical significance of the PERMANOVA test (p < 0.05) in inter-breed comparisons.

Table 4 β-Diversity values through Bray-Curtis, Jaccard, and unweighted and weighted UniFrac Metrics. P value for each comparison was obtained from PERMANOVA and considered significant at p < 0.05
Full size table

PCoA of the Bray-Curtis, Jaccard, and unweighted UniFrac Metrics revealed that microbiota in fecal samples were clustered in 4 different groups identified as our 2 experimental groups (Bray-Curtis—PC1: 38.12%, PC2: 12.97%, PC3: 9.78%; Jaccard—PC1: 30.65%, PC2: 12.72%, PC3: 10.44%; unweighted UniFrac—PC1: 42.75%, PC2: 12.23%, PC3: 8.31%); the repeatability of the results of the 5 animals within the same group was very high, particularly among the HFH and SIH groups (Fig. 3). PCoA of the weighted UniFrac distance revealed that microbiota of the animals of the same group had high variability (weighted UniFrac—PC1: 58.53%, PC2: 11.55%, PC3: 8.36%) (Fig. 3).

Fig. 3
figure 3

Principal coordinate analysis (PCoA) of 16S metagenomics data of microbial population in feces of Holstein-Friesian and Simmental heifers and cows. Comparison between community diversity based on different metric distances of microbial communities in cattle feces: a Bray-Curtis; b Jaccard; c unweighted UniFrac; d weighted UniFrac. Experimental groups are HFH (red), HFC (blue), SIH (orange), and SIC (green), each including 5 animals

Full size image

Abundance of species in comparative context: Holstein-Friesian heifers vs Holstein-Friesian cows and Simmental heifers vs Simmental cows

The analysis of composition of microbiomes (ANCOM) (Fig. 4) and volcano plot were performed in order to identify the relative abundance of species that was significantly different in the comparisons HFH vs HFC and SIH vs SIC (Kuczynski et al. 2012).

Fig. 4
figure 4

Volcano plot representation of ANCOM analyses. The horizontal axis is the centered log ratio transformation (clr) representative of the difference in abundance of a given taxonomical unit between Holstein-Friesian cows and Holstein-Friesian heifers at the family (red dots), genus (green dots), and species (blue dots) level. The W statistic indicates the value of the statistical test corresponding to the number of times the null hypothesis was rejected for each taxonomical class. Taxonomical classes with relative abundance significantly different (Holm-Bonferroni corrected p value < 0.05) are evidenced at the top of volcano plot and their taxonomy is reported. A darker shade of color indicates an overlap of OTUs

Full size image

The comparison between HFH and HFC associated with the ANCOM statistical analysis detected significant differences for 3 bacterial families according to Holm-Bonferroni post hoc test (p < 0.05); Anaeroplasmataceae (W 54) and Bifidobacteriaceae (W 52) were the most important families in HFC, and p-2534-18B5 (W 50) in HFH (Fig. 4). The statistical analysis of genera abundance in HFH vs HFC showed significant differences (p < 0.05) for Bifidobacterium (Bifidobacteriaceae) (W 103) and for one unclassified genus (p-2534-18B5). Results of the lowest taxonomic level revealed significant differences of species (p < 0.05) between HFH and HFC, suggesting a greater presence of Bifidobacterium pseudolongum (W 123) in HFC and p-25434-18B unclassified bacterium (W 116) in HFH.

In Simmental breed, Anaerolineales family (W 43) and 1 unclassified Methanobacteriaceae genus (W 96) and species (W 101) showed significant differences (p < 0.05) when comparing heifers and cows (Fig. 5). Anaerolinaceae family was more abundant in histidine-fed heifers (Klevenhusen et al. 2017). Indeed, SIH showed more abundant archaeal Methanocorpusculaceae and bacterial Elusimicrobiaceae family. As reported for Holstein-Friesian breed, Methanocorpuscolaceae role in rumen has not been studied yet and the presence in feces could be related to a shift outside from the rumen.

Fig. 5
figure 5

Volcano plot representation of ANCOM analyses. The horizontal axis is the centered log ratio transformation (clr) representative of the difference in abundance of a given taxonomical unit between Simmental cows and Simmental heifers at the family (red dots), genus (green dots), and species (blue dots) level. The W statistic indicates the value of the statistical test corresponding to the number of times the null hypothesis was rejected for each taxonomical class. The family, genus, and species with relative abundance significantly different (Holm-Bonferroni corrected p value < 0.05) are evidenced at the top of volcano plot and their taxonomy are reported. A darker shade of color indicates an overlap of OTUs

Full size image

Discussion

In the present study, we assessed the diversity of microbiome in cow and heifer feces within Holstein-Friesian and Simmental cattle breeds. The 16S rRNA sequencing approach was used to perform analysis of bacterial communities. The members of Firmicutes, Bacteroidetes, Proteobacteria, and Actinobacteria phyla were the most abundant and showed great variability among experimental groups (Kim et al. 2014). The presence of genus Fecalibacterium in the HFC and SIC might be related to a high corn percentage in the diet fed to cows compared with the diet fed to heifers. This hypothesis is corroborated by the greater abundance of Blautia and Roseburia which seem to be correlated with a corn-rich diet and are involved in the production of butyrate; the latter plays a role in the energy source for the mucosa (Kim et al. 2014), similar to Fecalibacterium which promotes the main energy sources for the gut epithelial cells (Pryde et al. 2002). The most representative genera, in all experimental groups, were correlated to digest functions, such as Clostridium, Eubacterium, Paeniclostridium, Romboutsia, and Turicibacter as related to diet composition and probably to physiological traits (Fan et al. 2017; Gerritsen et al. 2017; Xin et al. 2019).

Results of α-diversity analyses showed a significant difference between HFH and HFC, which is probably related to animal categories and diet composition. Finally, the OTU table was normalized according to the OTU abundance across samples, and the results were clustered as the heatmap to provide a better pattern across the experimental groups (Fig. 2). There was high difference of the clustered OTUs across the experimental groups; nevertheless, we cannot compare the breeds due to different farming conditions (diet and management conditions).

The PCoA showed the presence of two distinct groups (Fig. 3), which represent heifers and cows; indeed, the microbiome seems to be notably different between these two animal categories. We can indirectly justify the microbiome difference between heifers and cows, since diet (Callaway et al. 2010; Fan et al. 2017; Xin et al. 2019) and lactation (AlZahal et al. 2017) are important factors for the determination of bacteriological composition of fecal community. The variation of data could be explained as there is evidence in scientific literature indicating that the composition of fecal microbiome is completely different from the rumen and the differences between the taxonomical groups were less in the fecal samples compared with rumen (Azad et al. 2019). The differences between HFH and HFC and between SIH and SIC are strongly marked, especially with regard to some microbial populations, such as Tenericutes, Actinobacteria, and Bacteroidetes between HFC and HFH, and Chloroflexi and Euryarchaeota between SIC and SIH. This variability can be explained by diet composition in relation to the physiological status of the animals. However, the statistical analysis identified families, genus, and species present just in a specific category and breed. The variability of feces microbiome between the animals of the same experimental group was low, denoting a good reproducibility within samples. Interestingly the taxonomic class seems to cluster between heifers and cows of the 2 different breeds. Even if we could not perform a direct comparison between breeds, we can note that the trend in cows and heifers is conserved. Therefore, we may hypothesize that this trend is fixed, even when different interactions among breeds and diets are considered. Recently, a specific study reported a genetic link between breed and rumen microbiome community (Sandri et al. 2018); thus, as future perspective, we will investigate the microbiome composition of feces samples from animals of different breeds under the same farming and diet conditions. Finally, microbiome composition is crucial to evaluate feed efficiency, animal welfare, and methane emissions from cattle (Oikonomou et al. 2013; Holmes and Smith 2016; Sandri et al. 2018). Indeed, several alterations of the microbiome can compromise the physiological functions of the animals, altering their production capacity and health. Thus, our study lays the groundwork for future investigations aimed at developing minimally invasive and routine screening system to collect information on microbiome composition at population level. This would open new opportunities to select for more efficient and healthy animals.