Characterizing rumen microbiota and CAZyme profile of Indian dromedary camel (Camelus dromedarius) in response to different roughages

In dromedary camels, which are pseudo-ruminants, rumen or C1 section of stomach is the main compartment involved in fiber degradation, as in true ruminants. However, as camels are adapted to the harsh and scarce grazing conditions of desert, their ruminal microbiota makes an interesting target of study. The present study was undertaken to generate the rumen microbial profile of Indian camel using 16S rRNA amplicon and shotgun metagenomics. The camels were fed three diets differing in the source of roughage. The comparative metagenomic analysis revealed greater proportions of significant differences between two fractions of rumen content followed by diet associated differences. Significant differences were also observed in the rumen microbiota collected at different time-points of the feeding trial. However, fraction related differences were more highlighted as compared to diet dependent changes in microbial profile from shotgun metagenomics data. Further, 16 genera were identified as part of the core rumen microbiome of Indian camels. Moreover, glycoside hydrolases were observed to be the most abundant among all Carbohydrate-Active enzymes and were dominated by GH2, GH3, GH13 and GH43. In all, this study describes the camel rumen microbiota under different dietary conditions with focus on taxonomic, functional, and Carbohydrate-Active enzymes profiles.

Highly significant community level differences were observed between Fraction (PERMANOVA p-value < 0.001), Collection (PERMANOVA p-value < 0.001) and feed (PERMANOVA p-value < 0.001) based on the Bray-Curtis distances while less significant changes were contributed by breed (PERMANOVA p-value = 0.037) ( Table 1). Structural similarities of these communities were visualized through NMDS (nonmetric multidimensional screening) plot based on Bray-Curtis distance (Fig. 2). A clear separation of groups forming Liquid-Solid fractions and 5 collections were observed. Further, feed and collection were the most significant factors within liquid and solid fraction samples. Similarly, fraction and collection were the most significant factors within samples of each feed.
Shotgun metagenomics. Total 118.34 GB of shotgun metagenomic sequencing data (Table S8) from samples of the last two collections (n = 48) was analyzed using SqueezeMeta co-assembly pipeline. Co-assembly resulted in 4,277,503 contigs (> 200 bp) with 2.27 Gbp size. Additionally, 47.91%-77.36% reads per sample were mapped back to the assembly with an average of 65.69%. A total of 5,145,814 genes were predicted from all the contigs and annotated for COG and CAZymes. SqueezeMeta pipeline expresses gene abundance in the form of TPM (Transcripts per million) 15 . TPM is similar to RPKM (reads per kilobase per million reads) and represents the number of times a gene is observed per randomly sampled million genes.

Discussion
Rumen ecosystem harbors a great diversity of microbes with wide variety of roles. It is speculated that a diverse set of microbes start fermenting incoming feed particles and then a different set of microbes starts acting on fermented feed 17,18 . The microbes adhered to feed particles are quite different from those present in the fluid. These differences are not only limited to the taxonomy of microbes but to their functions/metabolism as well 12,17,19 . In this study, significant differences were observed among liquid and solid fraction microbiota at both taxonomic and functional levels as reported in previous studies on camel rumen 12 . At a higher taxonomic level, the proportion of two major phyla Bacteroidota and Firmicutes differed between both the fractions. In accordance with the previous studies on rumen, Firmicutes (Firmicutes_A and Firmicutes_C) were comparatively higher in solid fraction as compared to liquid fraction and vice versa for Bacteroidota phylum 20,21 . Firmicutes phylum was split in multiple phyla according to the taxonomy of GTDB 22 . GTDB follows the standardized bacterial taxonomy based on genome phylogeny and hence, differs from the traditional taxonomy from NCBI or other databases. However, significant differences were observed only in Firmicutes_A and Firmicutes_C and not in Firmicutes phylum. The reason being that Firmicutes_A and Firmicutes_C phylum includes class Clostridia and Negativicutes, respectively which are commonly associated with fiber degradation and reported in higher abundances in rumen, justifying their higher abundance observed in Solid fraction 12,17,20 . The Firmicutes phylum includes class Bacilli having lesser abundance in anaerobic rumen environment and are not modulated by other parameters 17,20 .
At genus level, Prevotella was the most abundant genus from Bacteroidota phylum. Other abundant genera from Bacteroidota phylum included CAG-462, RC9, RF16, F0040, F082 and some other taxa annotated at higher level. All these taxa belong to Bacteroidales order and were reconstructed from metagenomes. Bacteroides, a major genus from this order along with Prevotella are some of the most commonly observed genera in rumen of both ruminants and pseudo-ruminants and were reported in previous studies on bovine, sheep, goat, camel, alpaca [23][24][25][26] . The RC9 and RF16 genera were observed in higher abundance especially in liquid fraction. The previous studies on camel and moose rumen have also reported the higher abundances of RC9 12,27 and RF16 genera 25,28 , respectively.
Fibrobacter is yet another important member of the rumen community, reported in several studies. Fibrobacter is mainly associated with cellulolytic-fiber degradation and hence is an integral part of the rumen community. www.nature.com/scientificreports/ However, unlike previous studies in camel 8,12 present study did not reported significant difference in its abundance between liquid and solid fractions but observed similar abundance levels. Previous studies on rumen have also observed similar abundance of Fibrobacter in cattle and buffalo rumen 20,25 but lower abundance in other members of camelids 25 . Another genus Succiniclasticum known to ferment only succinate to propionate and not any other carbohydrates or amino acids 29 , which explains its significantly higher abundance in solid fraction compared to liquid fraction. Several bacteria from families Opitutaceae and Muribacullaceae were also observed in comparatively higher abundance in liquid and solid fractions, respectively. Members of Opitutaceae family have been isolated from soil, terrestrial environment and gut of ants and wood-feeding termites [30][31][32] . They have the ability of degrade lignocellulosic biomass and explains their presence in environments with plant biomass such as camel rumen. Muribacullaceae family was described very recently and the members of the family have been reported from mouse gut (for which it was named) and chicken caecum 33,34 . The members of this family are also shown to possess lignocellulosic degradation capabilities 33,34 .
From the functional point of view, around 33% of the annotated genes were annotated as "Functions Unknown" (COG Class-S). While higher (> 5% average) annotated COG class includes Class L (Replication, recombination and repair), Class J (Translation, ribosomal structure and biogenesis), Class G (Carbohydrate transport and metabolism), Class M (Cell wall/membrane/envelope biogenesis) and Class E (Amino acid transport and metabolism). The findings are in line with earlier studies on rumen, wherein a major proportion of genes involved in Genetic Information processing (Class L and J) followed by genes involved in Carbohydrate and Amino acid metabolism (Class G and E) were reported 14,17 . Previous studies have also described a higher degree of functional differences between fractions 17,20 .
Comparatively, lesser functional differences were observed between Collections (seven classes), and Feed (five classes), while no differences were observed between the breeds. While diet is one of the important factors influencing the shape of the rumen microbiome, few studies have observed that the taxonomic changes due to diet are more evident than functional changes 17,20,35,36 . Based on the results obtained in the study, it is speculated that the changes in the diet leads to the change in the abundance of the microbiota which gradually becomes stable under the influence of the same diet. Previous studies have shown that a period of 4-6 weeks can stabilize these diet related changes in rumen microbiota 37,38 . Incidentally, we did observe the greater number of significantly differentiating (Kruskal-Wallis, p-value < 0.05; not as per BH adjusted values) genera in Collection-2 (50) followed by Collection-4 (37), Collection-3 (33), Collection-5 (18), and Collection-1 (8), further substantiating the fact that indeed the most variation observed was immediately after diet change and decreased with time, whereas 0th day had the least feed-dependent variations as expected. This could also be the reason why lesser changes were observed between collections in the functional profile as they were based only on the last two collections.
The feeds included in present study were selected on the basis of their lignocellulosic content and therefore, we expected differences among the diet groups. However, we found less or no significant variations in taxa and diversity among the diet groups. The observed changes were also comparatively less pronounced as compared to similar experiments across multiple ruminants 25,39 as well as compared to variations in roughage-concentrate proportion of same feed 40,41 . While the functions of rumen in camels and cattle are similar, there are some of the differences associated with the physiology of animal which might result in less pronounced differences among feed associated microbiota. It is also probable that camels being able to survive on a wide variety of plant-based diets available in scarce environments, change in the diet might have lesser impact on camel rumen microbiota as compared to true ruminants. We also speculate that including larger group of animals in further studies can provide more reliable findings confirming the effects of change in diet on rumen microbiota. Another point worth mentioning here is the probability of the introduction of sequencing biases due to the presence of reagent and laboratory contaminants affecting the analysis 42 . The results of this study are therefore to be interpreted with caution as no negative-controls (no-template controls) were included in the study using which such contaminants can be identified and removed 43 .
With respect to CAZYme profile, most of the previous studies have reported a high proportion of GH followed by GT and other classes of CAZymes similar to present study 14,28,44,45 . We also observed similar dominant organisms containing these CAZymes, i.e., members of Bacteroidetes, Firmicutes, Fibrobacter. However, unlike previous study on camel rumen, we didn't observe higher contributions from Spirochaetes (0.6% in our study compared to 4% in other study) 14 . In line with previous studies on rumen 14,28,45 , we observed comparable proportions of GH families acting as cellulases (GH5, GH9, GH88, GH95), hemicellulases (GH8, GH10, GH11, www.nature.com/scientificreports/ (arabino/xylosidases) and GH2 (β-Galactosidases) were the most dominant families observed as in other studies of rumen of camel 14 , cattle 18,45,46 , buffalo 47 and moose 28 . The members of Bacteroidetes and Firmicutes were the major contributors for these four GH families, especially GH2, while more than 1% ORFs in GH13 and GH43 were coded by Fibrobacteres; and by Proteobacteria in GH3 and GH13. The contributions from Eukaryotes were also observed in GH3 (0.9%; 0.3% by Neocallimastigaceae, rest unclassified), GH43 (0.03%; entirely by Neocallimastigaceae) and GH13 (2%; 0.03% by Neocallimastigaceae family and 0.3% by Eudiplodinium genus). Eudiplodinium genus is a group of rumen ciliates belonging to family ophryoscolecids and have been linked with their cellulolytic and amylolytic activities 48,49 .

Materials and methods
Experimental design and sample collection. To access the dietary impact on the camel rumen microbiome, Kachchhi (K) and Bikaneri (B) breeds of camels were fed with three different diets, Bajra (B) (Pennisetum glaucum, pearl millet), Jowar (J) (Sorghum bicolor, sorghum) and Makai (M) (Zea mays, maize). The experimental animals were housed at the National Research Centre on Camel (NRCC), Bikaner, Rajasthan and provided ad libitum feed consumption and free access to drinking water. Twelve animals were divided into three groups (four animals in each group; two animals each of Bikaneri and Kachchhi breed) for 63 days ( Figure S1). Prior to the experiment, all the animals were maintained on the same diet based on Guar (cluster bean, Cyamopsis tetragonoloba), different from the experimental diets. Rumen liquor samples were collected using probang as mentioned earlier 50 under mild sedation. The samples were collected at 0 day before starting the feeding trial and subsequent collections were made on 10th, 21st, 42nd and 63rd days of experiment. We decided to collect the samples on every 21 days (21st, 42nd and 63rd) to cover the period of feed adaptation 37,38 and intermediary collection during initial week on 10th day. Collected rumen content was filtered through four-layered sterile muslin cloth to separate the liquid and solid fractions to be collected in 2 ml cryovials prefilled with Qiagen RNAprotect Bacteria reagent (Qiagen, Germany) at an approximate 1:1 ratio. Samples were immediately stored at − 20 °C in a portable freezer and transported to the laboratory where these samples were stored at − 80 °C until further processing.
Extraction of metagenomic DNA. Metagenomic DNA was isolated from liquid and solid fractions of rumen samples using QIAamp Fast DNA stool Mini Kit (Qiagen, Germany) following the manufacturer's instructions with minor modifications. Briefly, liquid samples were subjected to bead beating in Qiagen Tis-sueLyser for 30 s at 25 Hz and subsequently processed for lysis as per manufacturer's instructions. Rumen solid samples were vortexed for 20 min to completely dissociate bacteria attached with feed particles followed by centrifugation at 2600 g for 30 s to separate solid particles. Approximately, 600 μl of supernatant was processed from the previous step for DNA isolation as recommended by the kit manufacturer. Quantity and quality of metagenomic DNA was assessed using a Qubit 3.0 fluorometer (ThermoFisher scientific, MA) and agarose gel electrophoresis, respectively.
Library preparation and sequencing. V3-V4 hypervariable region of 16S rRNA gene was amplified using universal primer pair, 341F and 785R 51 and library was prepared according to Illumina 16S Metagenomics library preparation guide (Illumina, USA). The final library size and concentration was checked using Agilent Bioanalyzer DNA 1000 chip (Agilent, USA) and Qubit fluorometer (Invitrogen, USA), respectively. Four sequencing runs were carried out using prepared libraries on Illumina MiSeq sequencer employing 2 × 250 v2 chemistry. Shotgun metagenomic libraries were prepared from samples of collection 4 and 5 (n = 24). Libraries were prepared from 1 ng of metagenomic DNA with Nextera XT DNA Library Prep Kit (Illumina, USA) using the manufacturer's protocol. Prepared libraries were quantified using Qubit 3.0 and checked for size on Agilent Bioanalyzer 2100 using DNA HS kit. Five sequencing runs were performed on Illumina MiSeq using 2 × 250 v2 sequencing chemistry to sequence all the metagenomic libraries.
Data analysis. The raw data of amplicon sequencing was manually curated and quality filtered (average qual score < Q30 and trimming last 10 nucleotides from R2 reads) using Prinseq-lite Perl script 52 . The quality filtered data was then imported in the R v3.6.1 environment and analyzed with the DADA2 package v1.14.0 53 . As per DADA2 pipeline for 16S data (https:// benjj neb. github. io/ dada2/ tutor ial. html) and big data pipeline (https:// benjj neb. github. io/ dada2/ bigda ta. html), data from four runs was analyzed separately and then merged at a later stage. Briefly, the steps followed were quality check, trimming (primers were trimmed from both pairs) and filtering (no Ns and no PhiX), and sequence variants were inferred by estimating error rates and denoising. Sequence variants were merged across paired data and then data from all the runs were merged to construct the amplicon sequence variant (ASV) table followed by chimera/bimera removal and taxonomy assignment. GTDBr89 (Genome Taxonomy Database) database was used to assign taxonomy to the ASVs 22 59 . Observed ASVs and Shannon diversity were calculated and compared among groups. Further, between sample/groups comparison was done based on Bray-Curtis distance and visualized by plotting Non-metric multidimensional scaling (NMDS) plot followed by group level comparisons using PERMANOVA test. The Phylum and Genus level taxonomy was compared between groups to identify group specific differences. All the comparisons of diversity indices and taxa abundance across different groups were done using non-parametric Kruskal-Wallis (for multigroup comparisons) and Wilcoxon tests (for two-group comparison). The p-values were adjusted by Benjamini- www.nature.com/scientificreports/ Hochberg correction and have been mentioned accordingly throughout the manuscript. All the statistical testing between multiple groups were done using the R packages ggpubr and vegan. The raw reads obtained from shotgun metagenomics were curated using Prinseq-lite Perl script with following parameters: minimum length = 50, length trimmed to = 190 (to remove G-biased tails; one of the runs with very poor tail-quality was trimmed to 150 nucleotides), and minimum average quality = 30. The quality filtered reads were analyzed using SqueezeMeta employing a co-assembly pipeline 60 . Within the pipeline, assembly was done using MetaSpades 61 , ORF prediction using MetaProdigal 62 , taxonomy assignment using Diamond 63 against NCBI RefSeq database and functional prediction using Diamond/HMM against COG database 64 . Further, predicted ORFs were annotated for Carbohydrate Active Enzymes using HMMer based approach within dbCAN2 65,66 . Ethical permission. The work described in this article was carried out with prior ethical approval of the institutional animal ethics committee of the National Research Center on Camel, Bikaner, Rajasthan (NRCC/ PSME/6(141)2000-Tech/). All procedures performed in studies involving animals were in accordance with the ethical standards of the institution or practice at which the studies were conducted. The work included noninvasive sample collection and no animals were harmed during the experiment. The study was carried out in compliance with the ARRIVE guidelines.

Conclusion
In all, we report an extensive overview of camel rumen microbiota under influence of different diets. We observed the differences among three different feed roughages although the differences were not as much prominent as those reported in true ruminants. The study also tracked the microbiota diversity changes through time-points. We observed the highest number of significantly differentiating taxa in Collection-3 (21st day) with respect to Collection-1 (0 day). This points to the fact that on introduction of a new diet, microbiota starts changing slowly and more prominently during the third to sixth week and reaches a stable level thereafter. This was also observed in case of lesser functional differences between Collection-4 and Collection-5. However, the highest degree of variations were observed between two fractions of rumen content similar to that of previous studies of similar nature. We also observed a higher proportion of GH2, GH3, GH13 and GH43 CAZy families prominently involved in biomass degradation and reported in several rumen microbiota. Overall, this study presents important insights into camel rumen microbiome which can serve as critical information to increase feed digestibility in camels through selective enrichment of rumen microbiota.

Data availability
All the raw sequencing data is submitted in NCBI under BioProject PRJNA603266 and available from SRA under accessions SRR13178665 to SRR13178784 for 16S data and SRR13205818 to SRR13205865 for Shotgun data. The R script used for analysis is available from github.com/ankit4035/camelrumenproject (https:// doi. org/ 10. 5281/ zenodo. 43089 48) for reproduction of the entire work.