Unveiling the Bovine Epimural Microbiota Composition and Putative Function

Numerous studies have used the 16S rRNA gene target in an attempt to characterize the structure and composition of the epimural microbiota in cattle. However, comparisons between studies are challenging, as the results show large variations associated with experimental protocols and bioinformatics methodologies. Here, we present a meta-analysis of the rumen epimural microbiota from 11 publicly available amplicon studies to assess key technical and biological sources of variation between experiments. Using the QIIME2 pipeline, 332 rumen epithelial microbiota samples were analyzed to investigate community structure, composition, and functional potential. Despite having a significant impact on microbial abundance, country of origin, farm, hypervariable region, primer set, animal variability, and biopsy location did not obscure the identification of a core microbiota. The bacterial genera Campylobacter, Christensenellaceae R-7 group, Defluviitaleaceae UCG-011, Lachnospiraceae UCG-010, Ruminococcaceae NK4A214 group, Ruminococcaceae UCG-010, Ruminococcaceae UCG-014, Succiniclasticum, Desulfobulbus, and Comamonas spp. were found in nearly all epithelium samples (>90%). Predictive analysis (PICRUSt) was used to assess the potential functions of the epithelial microbiota. Regularized canonical correlation analysis identified several pathways associated with the biosynthesis of precursor metabolites in Campylobacter, Comamonas, Desulfobulbus, and Ruminococcaceae NK4A214, highlighting key metabolic functions of these microbes within the epithelium.


Introduction
The rumen epimural bacteria refers to the microbiota that are firmly attached to the rumen wall [1,2]. Despite contributing to less than 1% of the overall rumen bacterial biomass [3], these bacteria play a critical role in host physiological development, tissue turnover, urea recycling, and competitive exclusion of transient pathogens [3][4][5][6][7][8]. Within the epimural population, microbes also confer protection to the epithelium from challenging ruminal conditions such as low pH, osmotic stress, and pathogen-blooms associated with ruminal dysbiosis [9,10]. However, unlike the feed and fluid-associated rumen microbiota, little is known about the impact of the host, age, breed, gender, and geographic location on the rumen epimural community, or whether there is a definitive core microbiota associated with the rumen epithelium. The current definition is that a core microbiota is not found at any particular taxonomic level, but instead is the core set of genes available to a community for performing necessary metabolic functions associated with its niche [11,12]. This is critical information for being able to define a healthy ecosystem, as well as to understand dysbiosis. Since modern feeding practices in cattle often create digestive disturbances, the stability and the metabolic functions of those microbiota at the interface between the digestive milieu and the host tissue are of great interest. Specifically, the rumen epithelium is densely populated with Proteobacteria, most of them phylogenetic neighbors to pathogenic species [13]. It is therefore plausible that a dysbiotic epithelial community would have the greatest potential for increased numbers of opportunistic pathogens and their derived products with high pro-inflammatory properties (i.e., endotoxins) that induce local and systemic inflammation and other health disturbances [14]. Consequently, it is necessary to define a core epithelial microbiota and describe its diversity and function to better understand and assess dysbiosis.
Compounding this lack of knowledge is a lack of standardization in the non-culturebased methods of analysis. The rapid advancement of the fields of sequencing and bioinformatics has created a large disparity in the methods of analysis, making it difficult to interpret data and advance the field of rumen microbiology. Technical differences are mainly associated with the selection of the amplified 16S rRNA hypervariable region, the use of different sets of PCR primers, DNA extraction protocols, sequencing platforms, and bioinformatics workflows [12,15,16]. It is essential that methodologies are standardized, and bioinformatics are closely assessed to ensure that accurate and comparable data is presented for interpretation and analysis. Therefore, a meta-analysis of all publicly available datasets encompassing the rumen bovine epimural community was performed to assess the weight that various environmental and analytical factors have on the data obtained from 16S rRNA sequencing.

Search Overview, Inclusion Criteria and Data Acquisition
All available studies related to cow epimural microbiota were systematically reviewed. The literature search included all studies in NCBI PubMed (http://www.ncbi. nlm.nih.gov/pubmed (accessed on 30 April 2020)) and Google Scholar (http://scholar. google.com (accessed on 30 April 2020)) that included (i) the terms "papillae", "epithelium, OR "epimural", (ii) published since 1998, (iii) in "cow" OR "bovine" OR "cattle" (iv) and using "amplicon" OR "16S rRNA" sequencing data. Additionally, NCBI Bio-Project (http://www.ncbi.nlm.nih.gov/bioproject (accessed on 30 April 2020)) was also screened for unpublished data. Sequence files from each study were downloaded from SRA (http://www.ncbi.nlm.nih.gov/sra (accessed on 30 April 2020)) or obtained directly from authors of previous studies. Our final dataset contained 16S rRNA gene sequences of rumen wall samples in cattle from 11 previously published datasets and 332 samples ( Table 1). The datasets were further merged into 8 projects when belonging to the same study.

Sequence Analysis
Raw sequences from each project were analyzed and checked for sequence quality and adapter contamination using FASTQC [25] and AdapterRemoval [26]. Microbial community analysis was performed with QIIME2 v2020.2 [27]. Paired-end sequence data was merged using vsearch join-pairs [28] with the option -p-maxee 2 and quality filtered using the q-score-joined plugin with a minimum acceptable PHRED score of 20 (-p-minquality 20). For some datasets (ADDA and Sugarhay), these criteria had to be adjusted in order to allow more input reads (Supplementary Table S1). Denoising was performed with deblur [29] using a-p-trim-length of 200 bp. Representative sequences and feature tables obtained for each project were merged and filtered in order to exclude all sub-operational taxonomic units (sOTU) classified as mitochondria or chloroplast sequences. Taxonomy was assigned to sOTU using a classify-sklearn naive Bayes taxonomy classifier against the SILVA 132 99% OTUs reference sequences [30], resulting in a total of 56,025 features. Alphaand beta-diversity metrics were estimated after samples were rarefied to 8000 sequences per sample to account for uneven sequencing depth, resulting in the exclusion of 124 samples and two entire studies [23,24]. Core microbiota was determined by extracting the core features at phylum, family, genus, and sOTU level in at least 90% of the samples using the compute_core_microbiome.py command in QIIME 1.9.1. [31]. The final sOTU table was used for functional metagenomic prediction using PICRUSt2 [32]. The pathway level predictions are given based on MetaCyc [33] pathways and IDs (Supplementary Table S2).
To determine pairwise associations between core microbiota and predicted functions, regularized canonical correlation analysis was performed using the R package "mixOmics" (version 6.12.2) in R studio (version 4.0.2) [34].

Statistical Evaluation of Data
Alpha diversity metrics were checked for normality using the Shapiro-Wilk test (p < 0.05) and compared among studies using a non-parametric Kruskal-Wallis test with a pairwise Wilcoxon rank sum test in the R software environment [34]. Significance was considered at a p-value of < 0.05 after correction with the Benjamin-Hochberg procedure [35]. Principal-coordinate analysis (PCoA) was used to visualize distances between beta-diversity matrices. Permutational multivariate analysis of variance (PERMANOVA) using adonis [36] with 999 permutations and analysis of similarities (ANOSIM) was used to analyze the unweighted and weighted UniFrac distances and the Bray-Curtis dissimilarities matrices regarding study, country of origin, farm, gender (cow vs. steer), individual variability, biopsy location, hypervariable region, and primer set used. Permutational analysis of multivariate dispersions (PERMDISP) was used to test the homogeneity of dispersion for each metadata category analyzed previously [36]. Linear discriminant analysis effect size (LEfSe) (https://huttenhower.sph.harvard.edu/galaxy/ (accessed on 30 April 2020), The Huttenhower Lab, Boston, MA, USA) was used to determine which genera were significantly enriched according to the hypervariable region sequenced, and which metabolic pathways were overestimated in the different rumen locations. Genera and pathways that were more abundant in a particular group were identified by LEfSe using the non-parametric factorial Kruskal-Wallis (KW) sum-rank test (p < 0.05) and the effect size of each of these genera was estimated using linear discriminant analysis [37]. A LDA score (log10) of 4.0 was used as the cut-off to identify differentially abundant microbial genera and a score of 3.0 was used to identify differentially abundant metabolic pathways in the different rumen locations.

Composition and Structure of the Core Epimural Microbiota
We were able to identify 11 studies, grouped in 8 datasets, that met the search criteria used (Tables 1 and 2) comprising of a total of 332 samples. Core microbiota was defined as the taxa present in 90% of the samples, allowing the identification of 6 phyla

Core Bacterial Phyla Are Shared across Studies
Only Firmicutes, Proteobacteria, and Epsilonbacteraeota were present in each single sample analyzed. The latter is a newly proposed phyla, that derives from the combination of Epsilonproteobacteria, formerly included in Proteobacteria, and of the order Desulfurellales (Deltaproteobacteria) [38]. The same authors recently proposed to replace the name of this new phylum with Campylobacterota [38,39], given that Campylobacter is one of its most widely known members. The six core phyla detected correspond to the most abundant phyla across all the studies, with the exception of Spirochaetes, which was not identified in the project ADDA [18,19]. While Neubauer et al. [19] identified this phylum only in the particle-associated microbiota, Petri et al. [18] found Spirochaetes to be present at low abundance on the ruminal epithelium. Since the methods were extremely similar for both studies, the presence or absence of this phylum could be due to the bioinformatics analyses performed. Proteobacteria and Firmicutes, in different proportions, were the most abundant phyla across all the studies analyzed. In accordance with our findings, the studies from De Mulder and Neubauer [19,23] found Bacteroidetes to be less present in the epimural niche compared to other ruminal microenvironments.

Lower Taxonomical Levels Allow the Identification of a Core
Campylobacteraceae was represented in 100% of the samples, confirming the importance of the recently proposed Campylobacterota phylum in the epimural fraction of the rumen microbiota. Interestingly, Campylobacter (12.1 ± 0.66%) was also the only genus to be detected in all the samples and was also the most commonly described genus across all studies included in this meta-analysis ( Figure 3). Out of 344 bacterial genera, Campylobacter (15.5%), Kingella (7.8%), Desulfobulbus (4.7%), and Brachymonas (4.2%) were previously found to be the most abundant epimural bacteria [9]. Project Inflacow described the most abundant epithelial OTUs to be Campylobacter hyointestinalis, followed by Kingella oralis [9,10]. Results of our meta-analysis partially confirm these findings, as the feature classified as genus Campylobacter is 100% similar to Campylobacter hyointestinalis (NCBI Accession number: NR_115710.1). As previously discussed, most of the studied members of the Campylobacteraceae family are pathogenic strains, such as C. jejuni [40,41]. However, the exact role of these microbes in the epithelial community remains unknown. No other studies report the presence of K. oralis, which might suggest that this is a bioinformatics artifact or a study-specific effect. Two features identified as an uncultured Neisseriaceae (2.79 ± 0.19%) are classified in the NCBI database as Neisseria oralis (NCBI Accession number: NR_118249.1). Neisseria has been previously found to be associated with oxidative stress response in a transcriptome study of the cattle epimural microbiota using the same samples as Inflacow [42]. Likewise, most of the literature concerning Neisseriaceae refers to pathogenic bacteria. Members of this family are found in many different animal species, and a few strains are associated with microbiological communities populating the oral cavity in ruminants [43]. Some species, such as N. sicca, have been suggested to have urea hydrolysis activity in the rumen [44], whereas N. oralis has been demonstrated to be able to ferment lactose [45].
Features belonging to Ruminococcaceae UCG-014 (1.55 ± 0.07%), Ruminococcaceae UCG-010 (0.69 ± 0.04%), and Ruminococcaceae NK4A214 (6.76 ± 0.41%) group were identified as core in this meta-analysis, and have been previously reported as present in the whole gastrointestinal (GI) tract of cattle [46]. The latter is a relatively new genus, which has been recently added to the SILVA database, as it was previously included in broader classification at the family level [47]. Therefore, this group was not described in any of the studies included in this meta-analysis, as it was before considered among the Ruminococcaceae. Although bacteria belonging to this family are well known for their fibrolytic activity, Ruminococcaceae NK4A214 group were identified and cultured for the first time in 2011 [48], and their exact function in the rumen is still unknown. Burkholderiaceae is a Gram-negative family with an incredible ecological and metabolic variability, with several species classified as pathogens for humans, animals, and plants. Their exact role in epimural ruminal microbial community requires further investigation, but their capacity to adapt to challenging conditions could explain their presence in the core epimural microbiota [49], making the genus Comamonas (4.53 ± 0.32%) a very interesting epimural microbe. Family Desulfovibrionaceae, also identified as core, being able to reduce sulfate and produce acetate through the metabolism of lactate and hydrogen [50,51]. Families Lachnospiraceae and Desulfobulbaceae, identified as highly prevalent in three studies, [21,23,24], were confirmed as part of the core microbiota in our study. Lachnospiraceae UCG-010 (0.46 ± 0.02%) was also identified as part of the core microbiota. Lachnospiraceae are Gram-positive bacteria responsible for pectin fermentation and formate production [52], according to some authors. In other studies, the known role as fibrolytic bacteria [12,53] of unclassified Lachnospiraceae together with unclassified Ruminococcaceae was confirmed by being positively correlated with acetate production and an increased acetate to propionate ratio [54]. Despite being frequently identified as part of the ruminal microbiota, family Desulfobulbaceae is rarely examined and discussed. Some of its components, in particular belonging to the genus Desulfobulbulbus (5.79 ± 0.35%), have been described as sulfur compounds reducing and oxygen scavenging bacteria [23,24]. In contrast with many other genera, the role of Succiniclasticum (2.50 ± 0.18%) in the rumen seems to be clear: bacteria belonging to this group convert succinate to propionate, which is the most important source for gluconeogenesis in cows [53,55]. As a confirmation of its important role in the core, genus Succiniclasticum was identified as key component of the epithelial microbiota in five of the sequencing projects that we analyzed. The importance of the reference database is also evident for other taxa. Clostridiales Family XIII Incertae Sedis and Cardiobacteriaceae were previously identified as key components of the epithelial microbiota in the three studies analyzing the results at the family level. However, our results do not confirm these findings, even if we lower the threshold for microbiota core detection to 75% of presence over all the samples. Other studies have identified Prevotellaceae and Prevotella as part of the epithelial microbiota as well as of the other ruminal microenvironments [21,23]. The relatively high prevalence of this taxon in our meta-analysis would suggest an important role played by these microbes in epimural community, as well as in the other niches of the rumen [46]. Bacteria belonging to the family Prevotellaceae are known to be able to process monosaccharides and polysaccharides, peptides, and proteins, and are usually enriched in high-grain diets [12,53,56]. Nevertheless, since Prevotellaceae were detected in studies with different feeding regimes, it is unlikely that diet was the factor responsible for the differences observed. It is possible that diverse species within this family prefer metabolic pathways that are also highly adaptable to environmental stressors such as low pH conditions [57].
Methanobrevibacter was identified in all the studies and in 83.7% of the samples. Methanogens are the dominant Archaea within the GI tract of animals [58] and, depending on the hypervariable region of the 16S rRNA gene, some primers also amplify archaeal gene sequences [59]. The majority of the studies included in this meta-analysis did not report the presence of archaeal features amongst their sequences. Methanogens have been previously reported to have preferably either a particle-adherent, protozoa-adherent, or a free-living lifestyle [23]. However, metatranscriptome sequencing also found Methanobrevibacter as a highly abundant methanogen in epimural microbiota, with a high number of transcripts belonging to glycan biosynthesis, methane metabolism, and glyoxylate and dicarboxylate metabolism [42]. A high prevalence in epithelial samples might indicate that their presence in the rumen and significance in this niche has been underestimated.

Canonical Correlation Analysis Discloses Predicted Functions of the Core Microbiota
Regularized canonical correlation analysis was conducted between the pathways found in all samples, and the core genus microbiota of the epithelium (Figure 4). Several metabolic pathways were identified associated with the Ruminococcaceae NK4A214 group, Campylobacter, Comamonas, and Desulfobulbus (|r| ≥ 0.7).
Campylobacter was mainly associated with pathways involved in cytosine-monophosphate- Interestingly, Campylobacter is negatively correlated with starch degradation (PWY-6737, r = −0.72) and nicotinamide adenine dinucleotide (NAD) biosynthesis (PYRIDNUCSAL-PWY, r = −0.73). The function of Campylobacter in the rumen epithelium has for long been subject to debate. Campylobacteraceae are Gram-negative microaerophilic bacteria, which are usually associated with foodborne diseases in humans and diarrhea in pigs. Their presence in the cattle GI tract has been studied for years, as their commensal role in livestock is considered to be a reservoir for pathogen shedding [40,41]. Glutamate dehydrogenase and glutamine synthetase have been shown to be highly expressed by Campylobacter within the epimural microbiota [42]. This matches previous studies that have implied that Campylobacteraceae could be involved in protein metabolism. Quinol and quinones play a fundamental role in energy generating processes, as these compounds mediate respiratory electron transport. A periplasmic nitrate reductase napA and a thioredoxin reductase were previously found to be highly abundant in this genus [42]. This would be fundamental for anaerobic growth or auxiliary electron transport, allowing Campylobacter to anaerobically and aerobically succeed as a colonizer of the bovine gut, indicating a very likely role of Campylobacter in oxygen scavenging. . Some studies suggest that Comamonas possibly intervenes in nitrate reducing and complex organic compound degrading processes [24] or that it may have an important role as an oxygen scavenger, similarly to Desulfobulbus [23]. The function of these families also appears to be fundamental for the processes of transformation of volatile fatty acids at the epimural level [60]. Ruminococ-caceae NK4A214 group and Comamonas were positively correlated with sulfur metabolism, assimilation, and reduction pathways. Ruminococcaceae NK4A214 group was also negatively correlated with the biosynthesis of methionine, a sulfur-containing amino acid [61]. Microorganisms in the rumen reduce sulphate to sulfur, that can be included in organic compounds, such as proteins and essential amino acids [62]. Sulfate reduction has been previously assumed to have a negative effect on the rumen epithelium, given the toxicity of H 2 S [42]. Comamonas, on the other hand, was positively correlated with the TCA cycle, in which 5 -methylthioadenosine (MTA) is produced. This by-product can be used as source of methionine, sulfur, or purines [63]. In addition to their role in sulfur metabolism, our meta-analysis confirmed the role of Comamonas in nitrogen metabolism in the rumen. It was in fact positively correlated with proteinogenic amino acid degradation, and with purine nucleotide degradation. The latter pathway was also positively correlated with the Ruminococcaceae NK4A214 group. These findings suggest a fundamental function of both of these families in amino acidic and protein metabolism [64]. Our analysis did not confirm the role of Desulfobulbus in the metabolism of sulfur compounds, but we interestingly found other families being strongly correlated with such pathways. The core microbes identified in our meta-analysis were involved in amino acid metabolism to some extent, as previously suggested [42,65], highlighting a possible role of epimural microbes in nitrogen recycling, particularly the Ruminococcaceae NK4A214 group.

Factors Affecting the Epithelial Microbiota Composition and Structure
Comparisons between studies are not straightforward due to the sources of technical and biological variation inherent to each study. Differences in the laboratory methodologies employed, such as the DNA extraction method, the hypervariable region, PCR primers, sequencing platform, and bioinformatics pipeline all contributed to this variation [16]. Beta-diversity analyses showed unequivocal patterns of clustering mainly by hypervariable region sequenced. The PCoAs of Bray-Curtis dissimilarities ( Figure 5) and both weighted ( Figure 6) and unweighted UniFrac (Figure 7) distances also reveal a separation between studies, indicating a substantial similarity of taxonomic composition within sequencing projects. In the unweighted UniFrac plot, it is evident how each study represents a cluster, indicating a considerable qualitative difference in the microbiota composition between studies.

Technical Factors: Study Characteristics and Biases
All metadata categories were tested using PERMANOVA, ANOSIM, and PERMDISP ( Table 3). The individual study had the largest effect in the inter-sample diversity, according to the PERMANOVA results on weighted (R 2 = 0.97, p = 0.001) and unweighted (R 2 = 0.70, p = 0.001) UniFrac distance matrices, whilst the individual variability seems to have the largest effect in the Bray-Curtis dissimilarity matrix (R 2 = 0.65, p = 0.001). ANOSIM using the weighted UniFrac distance matrix revealed that the most significant separation is given by hypervariable region (R = 1.00, p = 0.001), followed by study (R = 0.94, p = 0.001), country (R = 0.94, p = 0.001), farm (R = 0.90, p = 0.001) and primer set (R = 0.89, p = 0.001). In the unweighted UniFrac, however, primer set (R = 0.95, p = 0.001) is the most important factor after the hypervariable region, while in Bray-Curtis dissimilarities, country (R = 0.94, p = 0.001) shows up as the second most important factor. However, almost all of these results can be associated with a possible dispersion effect, as PERMDISP was significant (p < 0.05) in most of the cases. Exceptions were found in the weighted UniFrac distances for country (p = 0.15) and in the Bray-Curtis dissimilarities for the hypervariable region (p = 0.98). Comparisons of alpha diversity measures revealed significant differences between V3-5 and V4 (p < 0.001). In fact, samples obtained through V3-5 amplicon sequencing had overall fewer of observed features and lower diversity than samples sequenced using the V4 hypervariable region ( Figure 8A,B). Table 3. Metadata categories affecting the epimural microbiota community structure assessed using permutational multivariate analysis of variance (PERMANOVA), permutational analysis of multivariate dispersions (PERMDISP), and analysis of similarities (ANOSIM).

Influence of Hypervariable Region
In order to understand how the hypervariable region affected the microbiota composition, linear discriminant analysis was used to determine which genera were significantly over-and under-represented. Sixteen bacterial genera were identified as being affected by the bias introduced with the choice of the hypervariable region (Supplementary Figure S2). The methanogen Methanobrevibacter and the bacteria Rikenellaceae RC9 gut group and Candidatus Saccharimonas seem to be more prevalent in samples sequenced using the V4 or V3-4 than in samples using V3-5. The use of the region V3-5 led to an underrepresentation of sequences belonging to Prevotella 1, Fibrobacter, Eubacterium nodatum group, Mogibacterium, Butyrivibrio 2, and Succinivibrionaceae UCG-001, while incrementing the relative abundance of Campylobacter, the Christensenellaceae R7 group, the Ruminococcaceae NK4A214 group, Ruminococcus 1, Desulfobulbus, Comamonas, and an uncultured Neisseriaceae. The high relative abundance of Epsilonbacteraeota in the V3-5 studies is reflected at the family level, where the family Campylobacteraceae has a high relative abundance, especially in the projects Inflacow and Sugarhay. The detection of Prevotellaceae may be depending on the hypervariable region of choice: in fact, Prevotella is absent in all the Austrian studies, which employed the region V3-V5. This data show that future research looking at the identification of microbial shifts as biomarkers of disease should consider studies using the same hypervariable region to improve comparability.

Biological Factors: Animal and Biopsy
Several studies on the structure of the epithelium microbiome of cattle have focused on the influence of diet [9,10,18,20,24], or feed additives [18,19,22], but fewer have focused on the heterogeneity among ruminal locations [21]. The number of animals sampled varied greatly among studies. Some studies only sampled the same animal once, while other studies included several time-points (longitudinal studies) under different experimental conditions. In the weighted UniFrac distance matrix, individual variability is as significant as the study per se (R 2 = 0.97, p = 0.001). The second most important factor is the biopsy location according to ANOSIM (R = 0.49, p = 0.001), which is again consistent with the results obtained for the other matrices. The most sampled ruminal location was the ventral sac (84% of the samples), followed by the dorsal sac (9% of the samples). Only one project included samples from different ruminal locations (ventral, cranial, caudodorsal, and caudoventral blind sac). Comparisons of alpha diversity measures revealed significant differences between biopsy locations (p < 0.001). A pairwise Wilcoxon rank sum test on the number of observed features identified significant differences between caudodorsal and dorsal (p < 0.01), caudodorsal and ventral (p < 0.01), caudoventral and ventral (p < 0.01), cranial and ventral (p < 0.01), and dorsal and ventral (p < 0.01) ruminal regions. The only significant difference was found in community richness between caudodorsal and dorsal sac. The ventral fraction was distinct from all the other ruminal locations. No differences in diversity were found between ruminal locations, except for the ventral sac (p < 0.01). However, this bias seems to be again associated with the use of the V3-5 hypervariable region. Samples of the ventral sac sequenced using V4 seem to be comparable to the other ruminal locations ( Figure 8C,D).

Predicted Functional Potential Varies across Ruminal Locations
The ruminal epithelium consists of papillae, which increase the surface area for the absorption of SCFA [66] but considerable functional variation exists across the spatial distribution of this tissue layer. To further investigate which functions are associated with the epithelial microbiota, PICRUSt2 and LEfSE were used to identify which pathways were enriched among the samples relatively to biopsy location (Supplementary Figure S2). Pathways associated with the caudodorsal region are mainly involved in degradation/utilization/assimilation, biosynthesis, and generation of precursor metabolites and energy. These pathways are involved in D-galacturonate degradation (GALACTUROCAT_PWY), carbohydrate biosynthesis (COLAN-SYN_PWY), fermentation to SCFAs (P108_PWY), guanosine-5 -diphosphate sugar biosynthesis (PWY_7323), starch (PWY_6737) and sugar (PWY_6507, PWY_7242) degradation, pyrimidine nucleotide biosynthesis (PWY_7199), and glycogen degradation (GLYCOCAT_PWY). Interestingly, pathways enriched in the dorsal area of the rumen are mainly associated with quinol and quinone biosynthesis (PWY_5897, PWY_5898, PWY_5899, PWY_5840, PWY_5838), metabolic regulator biosynthesis (PPGPPMET_PWY), and lipopolysaccharide biosynthesis (KDO_NAGLIPASYN_PWY). The ventral side of the rumen had two overexpressed pathways mainly associated with the generation of precursor metabolites and energy, such as the TCA cycle (P105_PWY) and aerobic respiration/electron transfer chains (PWY_3781). The ventral epithelium maintains a higher blood flow and is more papillated than the dorsal epithelium, being thus more implicated in nutrient absorption. Ruminal contractions and stratification of rumen contents are likely to affect dorsal and medial locations more than ventral sites [3]. Therefore, differences in specific metabolic functions within the different ruminal locations are very plausible.

A Closer Look into Austrian Samples
Samples from studies conducted in Austria (Projects Inflacow, ADDA, and Sugarhay) always used the same hypervariable region (V3-5) and biopsy location (ventral sac). Furthermore, the animals were housed in the same farm, allowing us to control the bias induced both by the farm and country. These samples were tested for the effect of the individual study, primer set, individual cow and diet. PERMANOVA on the weighted UniFrac distance matrix identified study (R 2 = 0.37, p = 0.001), followed by primer (R 2 = 0.22, p = 0.001) diet (R 2 = 0.10, p = 0.001), and individual cow (R 2 = 0.02, p = 0.015) as significant factors impacting the structure of the epimural microbiota. This is consistent with the permutational tests performed on the unweighted UniFrac distance matrix and the Bray-Curtis dissimilarity matrix. ANOSIM on the weighted UniFrac distance matrix shows that individual cow is not statistically significant (R = 0.02, p = 0.10). The composition of the GI tract microbiome, when geography and experimental setup is similar, is mainly affected by diet, despite the recognized influence of age, breed, environment and host genetics [12,46,67]. In fact, if we remove the sources of technical variation (study and primer), diet (p = 0.03) has a stronger effect than individual variability (Figure 9).

Conclusions
The richness, diversity, and complexity of the rumen ecosystem makes it a wellstudied microbial community. However, despite having been clearly identified as a unique component nearly 40 years ago, the epimural microbiome remains the least understood of the rumen niches. The re-analysis of 11 studies showed that a core set of microbes are associated with the rumen epithelium despite differences in experimental factors and methodologies. Sampling sites within the rumen showed variation in microbial communities and their predicted metabolism. However, this variation was compounded by the methods used for sequencing and bioinformatics analysis. Specifically, the use of a consistent variable region will be critical for future evaluations of the epimural community, since this factor alone greatly impacts diversity measurements and the comparability of datasets. Therefore, bioinformatics standards need to be addressed to ensure the accurate analysis of future data. Only by advancing our understanding of the microbes, their metabolism, and the critical factors affecting these components, can we apply precision methods for stabilizing rumen and animal health.
Supplementary Materials: The following are available online at https://www.mdpi.com/2076-260 7/9/2/342/s1, Figure S1. Abundance histograms of features at the genus level identified by LEfSe as differentially abundant according to the target region of the 16S rRNA gene, Figure S2. Metabolic pathways identified by LEfSe as differentially abundant according to the biopsy location, Table  S1. Bioinformatics processing of the sequences included in the meta-analysis, Table S2 MetaCyc pathways and BioCyc IDs.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.