Comparative analysis of the metabolically active microbial communities in the rumen of dromedary camels under different feeding systems using total rRNA sequencing

Rumen samples

Rumen samples were collected from 11 adult dromedary camels reared at three stations that use different feeding systems. Camels in group G1 (n = 3) were housed in the Maryout Research Station, Alexendria, Egypt and were fed on Egyptian clover hay (Trifolium alexandrinum), wheat straw and concentrates feed mixture. Camels in group G2 (n = 6) were housed at the commercial farm in the Kom Hammada and fed on fresh Egyptian clover (100% high-quality forage diet) then slaughtered in the Kom Hammada slaughterhouse, Elbehera, Egypt. Camels of group G3 (n = 2) were housed at the commercial farm in Cairo area and fed on wheat straw (100% low-quality forage diet) then were slaughtered in Pasateen slaughterhouse, Cairo, Egypt. Animals were kept on these diets for at least 1 month before the sampling time. The proximate analysis of feeds is illustrated in Table S1. Details regarding camel rumen samples in this study presented in Table S2. Rumen contents were strained immediately by two layers cheesecloth to separate the liquid and solid to form 22 samples, frozen using liquid nitrogen and stored at −80 °C before further processing (Elekwachi et al., 2017). The project was approved and all samples were collected according to the Institutional Animal Care and Use Committee, Faculty of Veterinary Medicine, University of Sadat City, Egypt (Approval number: VUSC00003).

RNA isolation, quality and quantity estimation and sequencing

The frozen rumen samples were ground using liquid nitrogen. About 0.5 g of frozen fine powder was used for total RNA isolation using Trizol-Reagent protocol (Invitrogen, Carlsbad, CA, USA), followed by RNA clean up using MEGA clear Kit (Invitrogen, Carlsbad, CA, USA). Total RNA quality and quantity were estimated using an Agilent 2100 bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and RNA 6000 Nano kit (Agilent Technologiess, Santa Clara, CA, USA, USA). One hundred nanogram of total RNA was reverse-transcribed into first strand cDNA and sequenced using Illumina rRNA MiSeq preparation kit (Illumina, San Diego, CA, USA) by Illumina MiSeq platform.

Bioinformatic data analysis

The generated RNA sequence reads were analyzed using pipeline developed by Elekwachi et al. (2017). Briefly, the sequence quality was checked using the FastQC program v. 0.11.4 (Andrews, 2010), then Trimmomatic program v. 0.35 (Bolger, Lohse & Usadel, 2014) was used to trim adaptors, barcodes, ambiguous and low quality reads. PEAR program v. 0.9.6 (Zhang et al., 2014) was used to merge read 1 and read 2 using default options. Then after, the hidden Markov models rRNA-HMM tool of the rapid analysis of multiple metagenomes with a clustering and annotation pipeline (RAMMCAP) (Li, 2009) was used to sort the reads into archaea and bacteria (16S, 23S), and eukaryote (18S, 23S) rRNA sequences. Merged sequence files were then sub-sampled as needed using MEME program v. 4.10.2 (Bailey et al., 2009). For each sample, 70,000 reads were run through the pipeline. For subsequent analysis steps, 20,000, 10,000, and 2,000 sequences were used for bacteria, eukaryote and archaea, respectively. Taxonomy binning for eukaryote and archaeal SSU rRNA sequences was performed using BLASTN. The sub-sampled query sequences were searched against the SILVA SSURef-111 database using an e-value of 1e⁻⁵. Bacterial SSU sequences were binned into operational taxonomic units (OTUs) using the “classify. seqs” command of Mothur v. 1.33.1 program (Schloss et al., 2009). The SSURef-108 gene and the SSURef-108b taxonomy databases were used. Principal co-ordinate analysis (PCoA) using Bray Curtis dissimilarity and alpha diversity indices (Chao1, Shannon and Inverse Simpson) were evaluated by Mothur (Schloss et al., 2009) based on sub-sampling of 70,000 reads per sample according the protocol “Community Structure Analysis Based on OTU Clustering” outlined in Elekwachi et al. (2017).

Statistical analyses

Data of relative abundance of bacterial phyla, protozoal genera, fungal genera and archaea genera and order Thermoplasmatales were tested for normality and homogeneity using Shapiro–Wilk test and variables that were deemed non-normal were then arcsine transformed. Linear Discriminate Analysis (LDA) and Bray Curtis Permutational Multivariate Analysis of Variance (PERMANOVA) tests depended on the relative abundance of bacterial phyla. All the protozoal, fungal and archaeal genera and the order Thermoplasmatales were used to show the differences in community structure and to compare the clustering of samples. Pearson correlation analysis was used to identify correlation within and between microbial communities and the correlation scores were visualized as a heatmap. The statistical analyses were performed using the SPSS v. 20.0 software package (SPSS, 1999) and PAST (Hammer, Harper & Ryan, 2001). Sequences were deposited to the sequence read archive (SRA) under the accession number: SRP107370.

The composition and diversity of active microbial community

Total rRNA sequencing in 22 solid and liquid rumen samples from 11 camels resulted in a total of 3,958,591 reads with average of 359,872 ± 85,366 (mean ± standard error (SE)) reads per animal in the solid fraction (SF) and 3,386,392 reads with an average of 307,854 ± 60,989 reads per animal in the liquid fraction (LF). The sequence reads of bacteria dominated the active microbial community, followed by protozoa, archaea and fungi (Table 1). Relative abundance of protozoa was higher in liquid fraction of G1 (LF-G1), while relative abundance of bacteria was higher in solid fraction of G1 (SF-G1). The highest population of archaea was observed in G2 camels. Additionally, G3 camels showed the highest relative abundance of fungi (Table 1; Fig. S1). Number of OTUs and Alpha-diversity indices, Chao1, Shannon and Inverse Simpson, were higher in the rumen of LF-G2 samples (Table 1).

Bacterial community

The composition of bacterial community varied little between groups and consisted of 12 phyla. The five most predominant phyla were Firmicutes, Bacteroidetes, Proteobacteria, Spirochaetes and Fibrobacteres, respectively (Table 2). Firmicutes dominated the bacterial community in all groups and was higher in G2 followed by G1 and G3 camels, respectively, and was also higher in SF compared to LF (Table 2). At the family level, Lachnospiraceae and Ruminococcuceae dominated the Firmicutes. In addition, six genera dominated this phylum, including Butyrivibrio, RFN8-YE57, Ruminococcus, vadinHA42, Acetitomaculum and Blautia (Fig. 1A; Table S3). The second largest phylum, Bacteroidetes, showed the highest relative abundance in G3 followed by G1 and G2 camels and was higher in LF than SF (Fig. 1A; Table S3). At the family level, Prevotellaceae, BS11_ gut_ group, and Rikenellaceae dominated the Bacteroidetes; and at the genus level, Prevotella, RC9_gut_group dominated the Bacteroidetes. Proteobacteria, phylum showed a higher relative abundance in LF-G1 samples and was dominated by Succinivibrionaceae family and Desulfovibrio genus (Table 2; Fig. 1A; Table S3). The Spirochaetes phylum was higher in the SF-G3 and it was classified into two families including Spirochaetaceae and PL-11B10 and was dominated by Treponema genus. Fibrobacteres phylum was higher in SF-G3 (Table 2; Fig. 1A; Table S3). Actinobacteria were higher in SF-G2 samples, Tenricutes phylum was higher in LF-G1 samples and Lentisphaerae phylum, was about 3-fold higher in LF as relative to SF and accounted for a large population in the camels of G3 (Table 2). Additionally, several minor bacterial phyla were also observed in the rumen of camels such as Verrucomicrobia, Elusimicrobia, Cyanobacteria and Chloroflexi (Table 2).

Of the 74 genera observed, only seven were observed exclusively in libraries generated from a specific facility, including uncultured Marinilabiaceae (Bacteroidetes), Quinella (Firmicutes) and Streptococcus (Firmicutes) that were observed only in G2 and G3 camels. Ruminobacter (Proteobacteria) was observed only in G1 and G2 camels. On the other hand, Arcobacter and Succinivibrio within phylum Proteobacteria were observed only in G1 camels and Betaproteobacteria (Proteobacteria) was observed only in G3 camels. Moreover, many unclassified bacteria were observed across samples and accounted for 39% of total bacterial reads. Most of these unclassified bacterial reads were observed in phylum Firmicutes and Bacteroidetes.

Archaeal community

Reads that classified as archaea were further classified to three orders within the phylum Euryacheota: Thermoplasmatales, Methanobacteriale and Methanomicrobial. Thermoplasmatales dominated the archaeal community and showed the highest population in LF-G3 samples, this order was not classified out of order level (Table 3; Fig. 1B). Reads that classified in the Methanobacteriale were further classified to family Methanobacteriacea that includes three genera: Methanobrevibacter, Methanophera and Methanobacterium. Methanobrevibacter is the second largest contributor in archaeal population and was higher in SF-G1 samples. Methanosphaera exhibited higher relative abundance in SF-G2 samples. Methanobacterium was absent in G3 camels; however, a small proportion of this genus was found in the camels of G1 and G2. Methanomicrobium genus, which belongs to order Methanomicrobiales and family Methanomicrobiaceae was the least contributor in archaeal population and was more prevalent in LF-G3 samples (Table 3; Fig. 1B).

Protozoal community

Reads that classified as protozoa were further classified to two families: Ophryoscolecidae and Isotrichidae (Table 4). Reads that classified in the Ophryoscolecidae were further classified to seven genera, Diplodinium, Ophryoscolex, Entodinium, Polyplastron, Eudiplodinium, Epidinium and Trichostomatia. Reads that classified in the Isotrichidae were further classified to two genera, Dasytricha and Isotricha. The variation among the camels in protozoal population was clearly observed and seemed to be higher than other microbial communities; however, the protozoal community composition was similar among the camels (Table 4; Fig. 1C). The most dominant protozoal genera were Diplodinium, Ophryoscolex and Entodinium. Camels in G1 had the highest population of Entodinium and Epidinium. Camels in G2 had the greatest population of Eudiplodinium, Ophryoscolex, Isotricha and Dasytricha and camels in G3 had the greatest population of Diplodinium, Polyplastron and Trichostomatia. On the sample fraction level, solid fraction had a higher representation of Ophryoscolex, Polyplastron, Eudiplodinium, Epidinium and Diplodinium, while liquid fraction had a higher representation of Entodinium, Isotricha and Dasytricha (Table 4; Fig. 1C).

Anaerobic rumen fungal community

Reads that classified as rumen fungi were further classified to two phyla: Neocallimastigomycota and Chytridiomycota. Reads that classified in the Neocallimastigomycota were further classified to family Neocallimasticeceae that includes three genera, Neocallimastix, Piromyces and Cyllamyces. Neocallimastix dominated the fungal community, followed by Piromyces and Cyllamyces (Table 5; Fig. 1D). These anaerobic fungal genera represented >99.5% of the fungal population. In addition, reads that classified in the Chytridiomycota were further classified to family Spizellomycetaceae that includes genus Spizellomyces, which was noted in a very small proportion (<0.5%) (Table 5). Neocallimastix was more abundant in the SF-G1 samples while Piromyces and Cyllamyces were more abundant in LF-G2 and SF-G3 respectively (Table 5; Fig. 1D).

Effect of feeding system and facility on the composition of microbial communities

Multivariate analysis separated libraries by feeding system and housing facility distinctly (Figs. 2 and 3). Also, bacteria, dominated by Firmicutes, drove differences between animals (Fig. 3). Furthermore, Entodinium, Thermoplasmatales, Neocallimastix drove differences in protozoal, archaeal and fungal communities, respectively. PERMANOVA analysis revealed that the difference between camel groups was significant (P < 0.01) in all microbial groups (Table S4). Pairwise comparison between camel groups based on Bonferroni-corrected P-value demonstrated that the difference was significant (P < 0.05) between camels of G2 and G3 in bacterial and archaeal communities (Table S4). Moreover, the difference was significant between the three groups in the protozoal community (P < 0.05), whereas, in the fungal community, the difference was significant only between camels in group G1 and G2 (Table S4).

Pearson correlation between microbes in the rumen of dromedary camel

Pearson correlation analysis (Figs. 4A and 4B), revealed many significant positive and negative correlations (P < 0.05). For example, in active bacteria, Bacteroidetes correlated positively with Cyllamyces and negatively with Butyrivibrio, Methanosphaera and Trichostomatia. Prevotellaceae correlated positively with Neocallimastix and Entodinium and negatively with Ruminococcaceae, Methanosphaera and Diplodinium. Fibrobacteres correlated positively with Cyllamyces, Methanomicrobium, Thermoplasmatales and Diplodinium and negatively with Methanosphaera, Epidinium, Ruminococcaceae and Butyrivibrio. Firmicutes correlated positively with Methanosphaera and negatively with Piromyces, Thermoplasmatales and Methanomicrobium.

In active archaea, Thermoplasmatales correlated positively with Diplodinium and negatively with Methanobrevibacter and Methanosphaera. In active protozoa, there was a negative correlation between Polyplastron, Entodinium, Ophryoscolex and Epidinium. In active fungi, a negative correlation was observed between Cyllamyces, Neocallimastix and Piromyces and between Piromyces and Entodinium.

Bacterial community

Firmicutes phylum was more abundant than Bacteroidetes and both phyla comprised >75% of all bacterial reads (Table 2), which agrees with studies of camels (Samsudin et al., 2011), Surti Buffalo (Pandya et al., 2010) and muskoxen (Salgado-Flores et al., 2016). The majority of Firmicutes’ members have a potential role in fiber digestion, which might illustrate their higher population in G2 camels that were fed on high-quality forage and also in solid fraction. The high proportion of Ruminococcaceae and Lachnospiraceae supports this speculation (Pitta et al., 2014a; Nathani et al., 2015). Blautia and Acetitomaculum genera have a key role as reductive acetogens (Le Van et al., 1998; Yang et al., 2016) and varied among the camel groups in this study. This supports the observation that manipulation of diet can enhance reductive acetogenesis in rumen and minimize methanogenesis (Le Van et al., 1998).

Bacteroidetes were higher in samples collected from animals reared in the station that used low-quality feed (G3), which was similar to results on cattle (Pitta et al., 2014b). The phylum was dominated by family Prevotellaceae, which confirms Gharechahi et al. (2015). Members of Bacteroidetes possess diverse enzymes that can target cellulose, pectin and soluble polysaccharides released in the liquid phase (Mackenzie et al., 2015). Additionally, Prevotella genus produces propionate that is used for energy by the host (Nathani et al., 2015). We speculate that Bacteroidetes species contribute to the adaptation of camels to arid conditions.

The RC9_gut_group found in this study belongs to uncultured genera and was found also in the Rhinoceros hindgut (Bian et al., 2013). Unclassified Bacteroidetes specialize in lignocellulose degradation (Mackenzie et al., 2015), which could support their high proportion in G3 camels. Fibrobacteres was higher (3.1%) in this study compared to the other findings on camels (Gharechahi et al., 2015); this phylum is the principal cellulolytic bacteria in the rumen (Ransom-Jones et al., 2012; Nathani et al., 2015), which might illustrate its higher relative abundance in solid fraction and in the rumen of G3 camels that fed on wheat straw (Table 2). The members of Proteobacteria were lower in G2 and G3 camels that were fed on diet rich in fiber contents. These findings highlighted this phylum’s function as a protein-degrading bacteria (Liu et al., 2017). The abundance of Treponema was higher in the solid fraction and in G3 camels (Fig. 1A). Treponema is the dominant genus in Spirochaetes phylum and it is fiber-associated bacteria, which could indicate to its cellulytic and xylanolytic activities (Ishaq & Wright, 2012).

The dominant bacterial genera in this study were Butyriovibrio, RFN8-YE57, Ruminococcus, Prevotella, Fibrobacter, Treponema and VadinHA. These genera were higher in the SF except RFN8-YE57 compared to the LF; this finding was consistent with a study on camels (Gharechahi et al., 2015), and confirms that solid-attached microbes could play a major role in ruminal fiber digestion (Jewell et al., 2015; Noel et al., 2017).

Archaeal community

Since some archaea produce CH₄ from H₂ and CO₂, this phyla may control methane emission from ruminants (Hook, Wright & McBride, 2010). Additionally, acetate produced in fiber breakdown provides a methyl group for methanogenesis; therefore, alteration of diet shifts the structure of methanogen populations (Hook, Wright & McBride, 2010; Tapio et al., 2017), which could demonstrate the variation in the relative abundance of archaea between camel groups. Camels of the second group (G2) that were fed fresh clover, showed the highest archaeal population (Table 2) and archaeal community was dominated by Thermoplasmatales, a methylotrophic methanogens order (Table 3), which was consistent with the results on cattle (Carberry et al., 2014) and camels (Gharechahi et al., 2015). Thermoplasmatales produce methane from methyl amine and supplementing of animal’s diet with rapeseed oil decreases the abundance of this order, making it a high potential target in future strategies to mitigate methane emissions (Poulsen et al., 2013). The Methanobrevibacter, Methanosphaera, Methanomicrobium and Methanobacterium (Table 4) dominated the reads classified as archaea in this study, which agrees with trends reported for beef cattle (Carberry et al., 2014). Methanobrevibacter dominated the methanogens in other ruminant (Henderson et al., 2015) and was associated with high methane emissions (Tapio et al., 2017). Moreover, Methanomicrobium was higher in the camels of G3 that were fed on poor quality forage, which was similar to results of buffalo (Franzolin & Wright, 2016), and in vitro (Wang et al., 2018). In rumen, Methanomicrobium converts H₂ and/or formate into CH₄ (Leahy et al., 2013). The abundance of Thermoplasmatales was also negatively correlated with Methanobrevibacter, which is consistent with previous results (Danielsson et al., 2017; McGovern et al., 2017).

Protozoal community

The majority of protozoal reads were classified as Diplodinium, Ophryoscolex, Entodinium, Polyplastron, Eudiplodinium and Epidinium (Table 4). Similar findings were observed on different ruminants (Baraka, 2012). Feed appeared to influence the relative abundance of protozoa, as reported previously for cattle (Hristov et al., 2001; Weimer, 2015); however, we cannot differentiate the effects of feed from facility. Diplodinium dominated protozoal community and was prevalent in G3 camels, which highlighted the cellulolytic activity of this genus (Coleman et al., 1976). Some species of genus Diplodinium, such as Diplodinium cameli, were discovered in, and are unique to, the rumen of Egyptian camel (Kubesy & Dehority, 2002). In addition, Entodinium was higher in G1 camels that were fed on concentrates feed mixture that increase the protozoa. Also, this genus predominates rumen of camels (Selim et al., 1999; Ghali, Scott & Jassim, 2005) and cattle (Carberry et al., 2012). Moreover, Kittelmann & Janssen (2011) showed that the Polyplastron was the dominant genus in cattle. On the function level, all the genus Eudiplidinum, Epidinum and Diplodinum have cellulolytic activity (Coleman et al., 1976), whereas, Polyplastrone and Epidinium have a xylanolytic activity (Devillard et al., 1999; Béra-Maillet et al., 2005).

Anaerobic rumen fungal community

The highest fungal population was observed in the solid fraction and rumen of G3 camels (Table 1). These findings were in agreement with the results of studies stated that the fiber-based diets stimulated the fungal growth (Orpin, 1977; Roger et al., 1993; Kamra, 2005; Haitjema et al., 2014). This could explain the low fungal population in G1 camels in our study. Moreover, the longer retention time and neutral pH in camel’s forestomach (Russell & Wilson, 1996) make it more suitable for the survival of rumen fungi. Neocallimastix dominated the fungal community and was higher in the G1 camels, which was similar to other results on sheep and camels (Kittelmann et al., 2013, Rabee et al., 2019). This genus produces enzymes capable of hydrolyzing cellulose, xylan and starch (Pearce & Bauchop, 1985). Cyllamyces that was observed in small population, has the ability to degrade poor-quality feeds (Sridhar, Kumar & Anadan, 2014), which might explain its high population in solid fraction and G3 camels. Piromyces was the second dominant genus in the camel rumen of this study and produces cellulolytic and xylanolytic enzymes (Teunissen et al., 1992). Therefore, this genus was most abundant in rumen collected from the G2 group of camels. The genus Spizellomyces is closely related to Chytridiomctes (Bowman et al., 1992), and common in grassland and crop soil (Lozupone & Klein, 2002, Kittelmann et al., 2012). Thus, contamination of forages by soil could explain the presence of this fungus in camel rumen.

Correlation between rumen microbes

Interactions between rumen microbes drive feed degradation and methane formation in the rumen, which influence the animal production and the environment (Williams et al., 1994; Lee et al., 2012; Henderson et al., 2015). Positive and negative correlations were observed within and between microbial communities in this study (Fig. 4). Methanogens colonize protozoa and this relationship enhances methane formation (Newbold, Lassalas & Jouany, 1995). Additionally, fibrolytic bacteria produce hydrogen and methyl groups that methanogens use for growth (Johnson & Johnson, 1995), which demonstrated positive correlations found between Fibrobacteres and some methanogens. Also, positive correlation between methylotrophic Methanosphaera and Lachnospiraceae that has been implicated in pectin degradation and provides methanol as a substrate for the methylotrophs (Dehority, 1969). On the other hand, Prevotella is a hydrogen utilizer and produces propionate that impact the methanogenesis in the rumen negatively (Pitta et al., 2014a; Liu et al., 2017), which illustrates negative correlation between Prevotellaceae and archaea.

Since the rumen anaerobic fungi produce abundant H₂ through the fermentation of carbohydrate; they can interact positively with H₂ utilizers such as archaea, Prevotellaceae, Blautia and Acetitomaculum (Orpin & Joblin, 1997; Le Van et al., 1998; Yang et al., 2016; Liu et al., 2017). Additionally, anaerobic fungi penetrate plant tissue, providing an increased surface area for bacterial colonization (Orpin & Joblin, 1997), which could explain positive correlation between fungi and both Butyrivibrio and Fibrobacteres. However, some bacteria and protozoa prey on fungal zoospores (Morgavi et al., 1994), which demonstrated the negative correlation between both Neocallimastix and Piromyces with Diplodinium and Entodinium. Furthermore, Ruminococcus produces compounds that inhibit the growth of rumen fungi (Stewart, Duncan & Richardson, 1992), which supports the negative correlation between Neocallimastix and Ruminococcaceae. Polyplastron predates upon other protozoa like Epidinium, Eudiplodinium, Diplodinium, and Ostracodinium (Eadie, 1967).