Microbial inoculum effects on the rumen epithelial transcriptome and rumen epimural metatranscriptome in calves

Manipulation of the rumen microbial ecosystem in early life may affect ruminal fermentation and enhance the productive performance of dairy cows. The objective of this experiment was to evaluate the effects of dosing three different types of microbial inoculum on the rumen epithelium tissue (RE) transcriptome and the rumen epimural metatranscriptome (REM) in dairy calves. For this objective, 15 Holstein bull calves were enrolled in the study at birth and assigned to three different intraruminal inoculum treatments dosed orally once weekly from three to six weeks of age. The inoculum treatments were prepared from rumen contents collected from rumen fistulated lactating cows and were either autoclaved (control; ARF), processed by differential centrifugation to create the bacterial-enriched inoculum (BE), or through gravimetric separation to create the protozoal-enriched inoculum (PE). Calves were fed 2.5 L/d pasteurized waste milk 3x/d from 0 to 7 weeks of age and texturized starter until euthanasia at 9 weeks of age, when the RE tissues were collected for transcriptome and microbial metatranscriptome analyses, from four randomly selected calves from each treatment. The different types of inoculum altered the RE transcriptome and REM. Compared to ARF, 9 genes were upregulated in the RE of BE and 92 in PE, whereas between BE and PE there were 13 genes upregulated in BE and 114 in PE. Gene ontology analysis identified enriched GO terms in biological process category between PE and ARF, with no enrichment between BE and ARF. The RE functional signature showed different KEGG pathways related to BE and ARF, and no specific KEGG pathway for PE. We observed a lower alpha diversity index for RE microbiome in ARF (observed genera and Chao1 (p < 0.05)). Five microbial genera showed a significant correlation with the changes in host gene expression: Roseburia (25 genes), Entamoeba (two genes); Anaerosinus, Lachnospira, and Succiniclasticum were each related to one gene. sPLS-DA analysis showed that RE microbial communities differ among the treatments, although the taxonomic and functional microbial profiles show different distributions. Co-expression Differential Network Analysis indicated that both BE and PE had an impact on the abundance of KEGG modules related to acyl-CoA synthesis, type VI secretion, and methanogenesis, while PE had a significant impact on KEGGs related to ectoine biosynthesis and D-xylose transport. Our study indicated that artificial dosing with different microbial inocula in early life alters not only the RE transcriptome, but also affects the REM and its functions.


Gene expression profile and related pathway analysis
The RE tissue had a total of 270 differentially expressed genes (DEG) by pairwise comparison between the treatment groups.There were 127 DEGs between the BE and PE groups (adjusted-p < 0.05 and fold-change > 1.5), with 13 up-regulated genes (URGs) in the BE group and 114 URG in the PE group.There were 36 DEGs between ARF and BE, with 27 URGs in ARF and 9 URGs in BE.There were 107 DEGs between ARF and PE, with 15 URGs in ARF and 92 in PE (Supplementary Table S1).
Using the DEGs obtained by pairwise comparison between the treatment groups, we performed the GO analysis and used REVIGO to remove redundant GO terms.We found 21 enriched GO terms between ARF and PE, 33 enriched terms between BE and PE, and no enriched terms were found between the ARF and BE groups.Then, we used the DEGs to perform a biological process analysis of GO.The analysis of the GO terms related to Biological Process (BP) presented the same GO terms than the analysis of all the GO categories (molecular function (MF), cellular component (CC) and biological process (BP), with 21 enriched GO terms between ARF and PE, 33 enriched terms between BE and PE, and no enrichment between ARF and BE (Supplementary Table S1).

Rumen epimural microbiota community structure
The three types of inoculums resulted in no differences in Evenness, Shannon diversity index, and Simpson's index.However, the Observed genera and Chao1 estimates differ by inoculations (Table 1), whereas microbial Observed genera was significantly higher in PE inoculum (p < 0.05).

Taxonomic microbial profile and microbial signature by inoculum type
After using SILVA reference to classify the microbial data, taxonomic profiling revealed a total of 125 taxa at the phylum level and 2,341 taxa at the genus level in the rumen epithelium.A complete list of all taxa classified is provided in Supplementary Table S2.
The sPLS-DA multivariate analysis implemented in the mixOmics R package was used to identify microbial taxa that best characterize each treatment group with 95% of confidence.For this analysis, only microbial taxa with a relative abundance > 0.01% across all the samples were considered.After the centered log-ratio transformation procedures, it was observed a separation in microbial taxonomic profile differentiating the rumen microbiota in ARF, BE and PE groups (Fig. 2A).The BE and PE taxonomic profile overlapped, indicating that the structures of the microbial communities from both groups were partially similar.On the other hand, the functional profile of the BE group showed a clear separation from ARF and PE (Fig. 2B).In contrast to the taxonomic profile, the functional profile showed that ARF and PE groups overlap, while the BE group is clearly different (Fig. 2B).
Overall, 40% of the microbial signatures selected in component 1 of the sPLS-DA characterized the rumen microbiota of animals treated with PE inoculum, which included members of the taxa RBG−16−49−21, Spirochaeta, Butyrivibrio, and Fretibacterium.Alternatively, the microbial signature that characterized the BE group represented 30% of the microbial signature selected in component 1, and contain the taxa U29−B03, Lachnospiraceae UCG−008, and Synergistes.The microbial signature that characterized the ARF group represented 30% of the microbial signature and included the taxa Alistipes, Prevotella, and NK4A214 group (Fig. 3).
In accordance with the sPLS-DA analysis (Fig. 2), the heatmap showed a clear difference in the rumen microbiota related to the different types of microbial inoculations (Fig. 4), with specific genera strongly related to the ARF inoculum type, differentiating this treatment group from the BE and PE.

Functional microbial signature by inoculum type
A total of 304 KEGG modules and 202 metabolic pathways were identified using the HUMAnN2.Among these, 201 modules and 170 metabolic pathways showed the relative abundance of > 0.05% across all the samples.Analysis using microbial protein coding reads indicated that both BE and PE had a significant impact on the abundance of KEGG modules (p < 0.05) related to acyl-CoA synthesis (M00086), type VI secretion (M00334), and methanogenesis (M00347), while PE had a significant impact on KEGGs related to ectoine biosynthesis (M00033) and D-xylose transport (M00215) (Supplementary Table S3).

Rumen epimural microbiota and its association with host rumen epithelium expression changes
To investigate the relationship between host epithelial expression genes and the epimural microbiota, we performed correlation analysis between gene expression level and microbial taxa abundance using Spearman analysis in R. Microbial taxa were collapsed at the genus or last characterized level and filtered at 0.1% of relative abundance.We found the 30 most significant, unique gene-microbe correlations in the rumen epithelium www.nature.com/scientificreports/(p < 0.00001), where the magnitude of correlation (Spearman's Rho) ranged between 0.994 and 0.999 (Supplementary Table S4).

Discussion
This study aimed to elucidate the impact of different types of microbial inoculum on the RE transcriptome and REM in calves.We used an RNA-seq read-based characterization of the rumen epithelium of calves receiving three different types of artificial microbial inoculum.Although rumen microbiome analysis techniques have been mostly based on DNA sequences, RNA-seq analysis can be advantageous by elucidating accurately the active microbial functions 25 .
To our knowledge, this is the first study to use RNA-seq to explore the rumen epimural microbiota and its influence on the rumen epithelium transcriptome in inoculated, pre-weaned calves.Several genes in the  in the microbial taxonomic profile between BE and PE, which differed drastically from the taxonomic profile of the ARF group.Last, we found microbiome-host interactions that indicate the influence rumen microbiome composition on rumen epithelial gene expression.

Rumen epithelial transcriptome changes in calves treated with different types of microbial inoculum
For the top 5% most highly expressed genes in RE tissue across ARF, BE, and PE groups, 31 of them were uniquely expressed in the ARF group, 20 were uniquely expressed in the BE group, and 6 genes were uniquely expressed in the PE group.These results suggested that different sets of highly expressed genes responded to each inoculum type.The ARF group showed a greater number of uniquely expressed genes than the groups that received the inoculum with different microbial compositions.Compared to ARF, PE had a higher number of DEGs than BE did, suggesting that PE had a higher impact on the rumen epithelium transcriptome expression.PE is a protozoa enriched inoculum, and the functional role of protozoa in the cattle rumen is largely unexplored.Our findings indicated that further investigation on the molecular mechanisms underlying the crosstalk between rumen epithelial and its protozoan population may help yield meaning insights.Gene ontology (GO) is a widely used and well-documented annotation system that assigns molecular function, cellular component information, and biological process to gene products, allowing the comparison among subcellular structures 26 .In this study, we used the DEGs (Supplementary Table S1) for the Gene Ontology analysis in DAVID 27 .We found that the enriched GOs terms between ARF and PE groups were 100% in the category Cellular Component (Bonferroni ≤ 0.05) (Supplementary Table S1).Cellular components refer to the anatomic structures of the cell, where the gene product execute its function, and in our study the most enriched GOs were related to single-organism cellular process (GO:0044763, 62 genes, p-value < 0.01), developmental process (GO:0032502, 47 genes, p-value < 0.01), and anatomical structure development (GO:0048856, 45 genes, p-value < 0.01).Comparatively, animals with increased rumen absorption, as seen in those with high feed efficiency, have an increased expression of genes involved in metabolic processes, tissue morphogenesis, energy-generating pathways 28 ; cellular functions 21,28,29 ; and in short chain fatty acid (SCFA) absorption and metabolism 30,31 .On the other hand, for the PE group, no significant GO enrichment was detected.Using the unique genes for each group, we did a focused GO terms analysis in the category on biological process (BP) to further explore the effects of the treatments on the rumen epithelial function.We found enriched GO terms in ARF and BE, but not in the PE group.The BE group had several BP enriched terms related to the epithelium (GO:0030855-epithelial cell differentiation, 5 genes, p < 0.007; GO:0002064-epithelial cell development, 3 genes, p < 0.03; GO:0060429-epithelium development, 6 genes, p < 0.008), and related to proteins (GO:0008104-protein localization, 9 genes, p < 0.002; GO:1990778-protein localization to cell periphery, 3 genes, p < 0.05; GO:0072657protein localization to membrane, 6 genes, p < 0.0002; GO:0019538-protein metabolic process, 11 genes, p < 0.02; GO:1903076-regulation of protein localization to plasma membrane, 3 genes, p < 0.005; GO:0051246-regulation of protein metabolic process, 9 genes, p < 0.003).The ARF group had several BP enriched terms related to the cell localization (GO:0,051,641-cellular localization, 16 genes, p < 0.05; GO:0032879-regulation of localization, 9 genes, p < 0.05; GO:0040012-regulation of locomotion, 6 genes, p < 0.03), and related to mRNA processing (GO:0016071-mRNA metabolic process, 5 genes, p < 0.02; GO:0000398-mRNA splicing, via spliceosome, 3 genes, p < 0.02) (Supplementary Table S1).

Different inoculum types modified the taxonomic and functional composition of the rumen microbiota
Our results showed that the most abundant epimural microbial genera are prokaryotic.We found a high abundance of the genus Klebsiella, that are known for being present in mammary gland of cows and related to mastitis 32 , and has also been detected in samples collected in farms from water, soil, cattle feces, bedding, alleyway, and holding pen 32 .Klebsiella is not commonly found in a high abundance on rumen content samples, and the high abundance of this taxa found in our study might be explained by the fact that a different sample type, rumen epithelium, was analyzed in this study.When analyzing the rumen liquid metatranscriptome from the same animals, Park and co-authors 24 did not find this microbial genus, suggesting its presence is in the rumen epithelium only.Additionally, most of the studies exploring the rumen microbiome use DNA-based 16S rRNA amplicon sequencing methods, while we used whole transcriptome RNA sequencing approach.RNA-sequencing approach has two major advantages over 16S rRNA amplicon sequencing.RNA-sequencing is a random fragmentation based sequencing of the entire transcriptome, which is not limited by any pre-determined genomic regions like the 16S rRNA sequencing.Additionally, RNA-seq has the ability to identify actively transcribed microbial transcripts.For microbial taxa that have low abundance in a community, they might be difficult to identify using DNA-based, 16S rRNA methods.However, if these microbial taxa produce high number of transcripts, they can be identified using RNA-sequencing method.Thus, there can be "microbial blind spots" when different sequencing methods are used.Consistent with this, Park and co-authors 24 showed that RNA-sequencing and DNA-based 16S rRNA methods showed different microbiota profiles for the same rumen liquid samples, with a sizable number of microbial taxa being identified exclusively by RNA-seq method.
The overall diversity of the RE microbiome was lower in ARF when compared with BE and PE treatments, with a significant decrease in the number of observed genera and in the alpha diversity measured by Chao1.The greater number of observed genera with BE and PE could be due to the live microbial inocula promoting epimural microbial establishment.These results are consistent with the diversity measurements from rumen fluid using a 16S sequencing reported previously for these treatments 33 .However, these differences were not detected in the ruminal fluid from these same calves using RNA-seq 24 , further indicating that the microbial communities in the RE and rumen fluid are different.
The microbial taxa related to the BE treatment have several potential functions.The taxon U29_B03 was found to have a positive correlation with SCFA concentration on the rumen 34 .The UCG−008 is a butyrate producing bacteria belonging to the family Lachnospiraceae 35 .Butyrate is a potent stimulator of epithelial proliferation 36 , and an increase in epithelial proliferation results in an increase in rumen absorptivity and is related to high feed efficiency in cattle 30 .The genus Synergistes plays a role in methane production by cooperating with methanogens (e.g., Methanomicrobia), conducting interspecies hydrogen transfer 37 .Since the methane produced in the rumen is not absorbed by the animal, the increase in the methane production represents an energy loss to the animal 19 .
Among the genera related to the PE taxonomic signature, only three genera have previously reported functions in the rumen.Spirochaeta plays a role in plant biomass degradation through the secretion of glycoside hydrolases 38 .Butyrivibrio encodes a diverse spectrum of degradative carbohydrate-active enzymes (CAZymes), www.nature.com/scientificreports/which degrade polysaccharides to yield volatile fatty acids, that are used by the host for growth 39 .Additionally, Fretibacterium was found to be positively correlated with greater amounts of long-chain fatty acids, including alpha-linolenic acid, nervonic acid, and palmitic acid, that are related to the growth of specific microbial groups on the rumen 40 .
For the microbial genera related to the ARF taxonomic signature, two genera have known functions.Alistipes play a role in the degradation of plant-derived polysaccharides 41 .Members of the genus Prevotella have been associated with animals with high feed efficiency and animals with low feed efficiency, indicating that the species within this genus play different roles in rumen fermentation 19 .
According to the sPLS-DA analysis and the alpha-diversity indexes, the taxonomic profile of the RE microbial community is different in calves receiving BE and PE from the animals that received ARF.However, the functional profile shows that the BE and ARF groups are closer to compared to the PE group.The BE and PE groups shared several KEGG modules, and some modules were found exclusively in PE, and no specific KEGG module was found for the ARF group.Studies have shown that the functional profile of the rumen microbiome is more related to the health of the host and the milk production than the taxonomic profile 19,42 .Therefore, despite the perceived composition and functional differences and similarities among the treatment groups, it is still unclear if or how the inoculation of different microbial groups could improve the rumen fermentation and, consequently, improve animal health and production traits.

Microbiome-host interactions may influence rumen epithelial gene expression
The RE microbiome lies at the interface between the host and its gut environment, and the microbial activity in the RE directly influences the metabolism and physiology of the host animal 43 .The top 30 overall interactions within the rumen epithelial microbiome, which is represented by the abundance of major microbial genera and the DEG, were compared.A total of 38 nodes and 30 significant interactions between the RE microbiome and DEG were significant.The genera Roseburia, Anaerosinus, and Succiniclasticum were considered keystone taxa based on centrality measurements.Roseburia was strongly related to the 25 URGs, Entamoeba related to 2 URGs, and Anaerosinus, Lachnospira, and Succiniclasticum were related to one URG each (Fig. 5).
The genus Roseburia plays a role in starch fermentation, being a prominent butyrate producer 44 , and in sheep was found to have an increased abundance in animals with high feed efficiency 45 .In our study, the enriched GO terms from DEGs correlated to Roseburia were related the cell structure (GO:0030018-Z disc) and (GO:0042383sarcolemma) (Supplementary Table S1), while the DEGs correlated the other microbial genera presented no enriched GO terms.The genus Anaerosinus has the ability to ferment glycerol to propionate 46 , and the genus Succiniclasticum was reported to convert succinate to propionate 47 .Butyrate and propionate syntheses compete for the same substrate that archaeas use for methanogenesis in rumen, the hydrogen 48,49 .The genus Entamoeba comprises protozoan parasites hosted by vertebrates and invertebrates animals 50 .Entamoeba spp.have been detected from cattle feces in animals without any clinical symptoms 51,52 .It suggests that these protozoa can occur in cattle in a nonpathogenic form, and that their occurrence might be more common than previously thought.Our study has the first finding of Entamoeba spp. on the rumen of cattle.Our study used RNA-based methods to explore the rumen epithelium, and as cited before, it can explain differences in the microbial community composition when comparing with other studies about the theme.However, it is important to highlight that there are no rumen protozoa specific databases, and the bioinformatics analysis match the sequence against organisms described in other environment, which can cause potential misclassification of the microbes.
Our study suggests that different types of microbial inoculum alter the RE transcriptome and metatranscriptome.However short-and long-term implications of these results are not clear.Despite the ARF group having received autoclaved rumen fluid and have lower microbial diversity on RE, more DEGs were found in ARF when compared to PE, but not when compared to BE. Future studies should be conducted to elucidate the mechanisms of dosing on the epimural microbiome and their relation to the RE tissue function.For example, it is still not clear if the KEGG modules found in our study are related to the rumen establishment of the microbial taxa present from the inoculum or whether it is related to the interaction of these microorganisms with the microbes already present in the rumen (e.g., by predation, competition, etc.).Nonetheless, our findings suggest the possibility of manipulating the rumen epimural microbiome on calves in early life, and that the differences in the epimural microbial community influences the rumen epithelium tissue gene expression.However, these findings pointed out the potential manipulation of the epimural microbiome and RE in calves on early life, with the potential to alter host phenotype in dairy cattle.These results should be confirmed with further long-term studies to evaluate the effects on future productivity and feed efficiency in adulthood.
Our study suggested that different types of microbial inoculum alter the taxonomic and functional composition of the rumen epimural microbiota and transcriptome profile in the RE tissue.However, further studies are needed to better understand the functional impact of exogenous rumen content inoculation in young calves.Specifically, a better understanding is needed for how different types of rumen inoculum could improve rumen fermentation and, consequently, impact the host health, productivity, and efficiency.This work adds empirical evidence suggesting the feasibility of manipulating the RE microbiota through interventions in early life.More rumen developmental studies across different time points using RNA-seq are warranted to better understand how microbial populations and their functions influence host gene expression.

Conclusions
The rumen of dairy cattle contains microbiota from all domains of life that play a prominent role in digestion and may also affect animal health, production, and efficiency.To evaluate the effects of directed rumen microbial establishment on the gene expression of the rumen tissue and epimural microbiota, we dosed the rumen of pre-weaned dairy bull calves with either bacteria-or protozoa-enriched inoculum from adult cows.This is www.nature.com/scientificreports/95 °C for 2 min, then 40 cycles at 95 °C for 15 s and 60 °C for 60 s.The analyses were carried out in triplicate for each data point.The fold change in gene expression was obtained following normalization to two reference genes, ACTB and HMBS.These two reference genes were found to be very consistent in cattle 58 .The relative quantification of gene expression was determined using the 2 −ΔΔCt method 59 .

Mapping of RNA sequencing raw reads and accessing differential gene expression analysis
FastQC was used to check the quality of the raw reads (https:// www.bioin forma tics.babra ham.ac.uk/ proje cts/ fastqc/), and the raw reads were filtered to remove those shorter than 50 bp.For sequence alignment, the genome reference and the annotation file of Bos taurus were used (NCBI, ARS-UCD v.1.2) 60 and raw sequencing reads were aligned using STAR 61 .Cufflinks 62 was used to determine the expression level of mRNAs in each sample and to calculate the Fragments Per Kilobase of transcript per Million mapped reads (FPKM) for each gene.The total number of expressed genes was calculated using a FPKM cutoff value of 1. Differential expressed gene (DEG) analysis across the groups (ARF, BE, and PE) was performed using cuffdiff module in cufflinks v. 2.2 62 .DEGs analysis was done for the comparison between control group (ARF) and the groups BE and PE, and between BE and PE.Gene function annotation and pathway analysis were performed using DAVID 27 and stringDB v.11.5 63 .
The most highly expressed genes (top 5%) were first identified for each sample using FPKM values.Then, the most abundantly expressed genes across the three groups and the most highly expressed genes that are unique to each group were identified.To gain new insights into the underlying biological functions of DEGs, the Gene ontology (GO) pathway analysis was performed using stringDB v.11.5 63,64 and used to analyze the identified DEGs (p ≤ 0.05).

Figure 1 .
Figure 1.Taxa summary plot for microbial abundance across the ARF, BE, and PE groups.Considering the total abundance, the ten most abundant microbial genera were plotted.

Figure 2 .
Figure 2. Results of sPLS-DA for microbial the profile at genus level in the microbiome attached to the rumen epithelium of calves that received three types of inoculums.Individual score plot of the samples along the first two components, with a 95% confidence level.(A) taxonomic profile; (B) functional profile.

Figure 3 .
Figure 3. Discrimination of the taxa that best characterize each treatment.The loading plot displays the abundance of the microbial taxa on the treatment they are the most relevant.Ranked from the bottom (most important) to the top.Colors indicate the type of inoculum received in which the microbial taxa was most relevant.

Figure 4 .
Figure 4. Heatmap showing the most relevant microbial genera related to the different types of microbial inoculum.The ARF group was indicated by the blue bar, and the BE group was indicated by the orange bar, and the PE group was indicated by the gray bar.

Figure 5 .
Figure 5. Interactions between host tissue rumen epithelial genes associated with epithelial rumen microbiome.Network visualizing top 30 significant gene-microbe correlations.

Table 1 .
Microbial diversity measures from calves dosed with three different types of artificial dosing of rumen inoculum.Diversity measures are represented as the mean ± standard error (SE).