Characterizing the fecal bacteria and archaea community of heifers and lactating cows through 16S rRNA next-generation sequencing

The aim of this study was to describe the fecal bacteria and archaea composition of Holstein-Friesian and Simmental heifers and lactating cows, using 16S rRNA gene sequencing. Bacteria and archaea communities were characterized and compared between heifers and cows of the same breed. Two breeds from different farms were considered, just to speculate about the conservation of the microbiome differences between cows and heifers that undergo different management conditions. The two breeds were from two different herds. Firmicutes, Bacteroidetes, Actinobacteria, and Proteobacteria were the most abundant phyla in all experimental groups. Alpha- and beta-diversity metrics showed significant differences between heifers and cows within the same breed, supported by principal coordinate analysis. The analysis of Holstein-Friesian fecal microbiome composition revealed 3 different bacteria families, 2 genera, and 2 species that differed between heifers and cows; on the other hand, Simmental heifers and cows differed only for one bacteria family, one archaeal genus, and one bacteria species. Results of the present study suggest that fecal communities of heifers and cows are different, and that fecal microbiome is maintained across experimental groups.


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.

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.
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.

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). 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 twothirds of the total microbiome ( Fig. 1).

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 Grampositive 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 10 (11) 13 (6) 22 (14) 13 (10) Uncultured Eggerthellaceae 1 (2) 5 (6) 2 (3) 13 (8) Unclassified Actinobacteria 12 (7) 31 (16) 14 (20)    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  (Fig. 2).

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.

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. 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 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  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.

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.