Diverse methanogens, bacteria and tannase genes in the feces of the endangered volcano rabbit (Romerolagus diazi)

Sample collection

Four fecal samples were collected from Coajomulco, Morelos Mexico (14Q UTM 478960.07 m E, 2109691.32 m N) under a permit by Secretaría de Medio Ambiente y Recursos Naturales (SGPA/DGVS/006985/18). These pellets are ochre in color, with a smooth, shiny texture and their maximum diameter is one centimeter in adults. We kept more than 25 m between one latrine and another within a collection site, this distance is the estimated home range for these rabbits (Cervantes & Martínez-Vázquez, 1996), except for two samples (S1 and S2) that were collected in the same latrine with 0.1 m distance between them, but the pellets were of different sizes. Feces were stored individually in Eppendorf-type tubes in an icebox and later taken to the laboratory and stored in a refrigerator at 4 °C for 6 h.

DNA extraction, sequencing, and Illumina library generation

Total DNA of the four samples was extracted using the commercial DNA Isolation Kit from Roche Life Science. A total of 230 mg of feces were processed for each sample. DNA concentration was measured by Nanodrop (Thermo Scientific) and Qubit™ dsDNA HS Assay Kits and the quality was determined by visualization in a 1% agarose gel. The sample requirements for the preparation of the libraries were: DNA concentration >  200 ng/µL, total amount of DNA >  1 µg and DNA bands visualized in agarose had to be clear and of good quality. The genomic DNA was randomly cut into small fragments and subcloned into a “universal” cloning vector. The sub fragment library was randomly sampled, and several sequences reads (using a universal primer that directs sequencing from within the vector) were generated from both strands. Sequencing was obtained with Illumina HiSeq 4000 using paired-end 2 × 151 bp, at the sequencing unit Functional Genomics Laboratory at UC Berkeley, California.

Sequence filtering and metagenome assembly

Raw fastq reads were quality filtered with FASTP 0.20.0 (Chen et al., 2018). We considered good quality reads when the quality score was equal to or greater than 30. Clean reads were mapped against genomes such as human, rat, rabbit, yeast, mouse, pine, oak, fly, worm, adapters, and pika to know whether the reads had contamination using FastQ Screen 0.14.0 (Wingett & Andrews, 2018). Data is available in the SRA at NCBI, BioProject: PRJNA721235, accession number: SRR14209496, SRR14209495, SRR14209494, SRR14209493). Metagenome assembly was performed using metaSPADES 3.12.0 (Nurk et al., 2017) and MEGAHIT 1.1.3 (Li et al., 2015). Comparison of the assemblies was conducted to identify the best assembly with MetaQUAST 5.0.2 (Mikheenko, Saveliev & Gurevich, 2016). A good quality assembly considers the proportion of contigs <1000 nt, total length of the contigs, longest contig, N50 (length), L50 (position), fewer contigs and identifying single-copy genes (>90%) (Nurk et al., 2013). Clean reads (R1 and R2) were aligned to each corresponding assembly to obtain the percentage of reads that were assembled using Bowtie2 2.4.2 (Langmead & Salzberg, 2012). We confirmed that the samples were from volcano rabbits because we encountered DNA sequences in feces that matched those reported from the same rabbit species using the nucleotide sequences of the cytochrome b gene. We assessed taxonomic identity on fecal samples using Blast 2.10.0+ (Camacho et al., 2009) and a phylogenetic reconstruction using maximum likelihood (ML). Additional sequences of related organisms were obtained by Blast 2.10.0+ at NCBI. Sequences were aligned with Mafft 7.149b (Katoh & Standley, 2013). The best substitution model (TPM2+F+G4) and ML analysis were performed with IQ-TREE 1.6.12 (Kalyaanamoorthy et al., 2017; Trifinopoulos et al., 2016), using the UFBoot2 method (Hoang et al., 2018) with 10,000 replicates.

Taxonomic Assignment

The microbial taxonomy assignation of each assembly per sample was estimated using Kraken 2 2.0.8 (Wood, Lu & Langmead, 2019). The abundance of the microbial community by phylum and genus was estimated using Bracken 2.5.0 (Bayesian Abundance Re-estimation with KrakEN). Bracken could reassign sequences that kraken 2 could not classify within the genus or species levels with a reliable percentage above 98% (Lu et al., 2017). Both programs were run against a Minikraken2_v1 database containing bacterial, archaeal genomes and viral libraries. In addition, we searched for the 16S rRNA gene sequences in the metagenomes by Blast 2.10.0+. Redundant sequences were removed with CD-HIT 4.8.1 (Li & Godzik, 2006). Phylogenetic reconstruction was performed using the maximum likelihood (ML) method. Additional sequences of related organisms were obtained by Blast at NCBI and with the Refseq-RDP database. Sequences were aligned with Mafft 7.149b. Spurious bases, editing and trimming were performed with trimAl 1.4 (Capella-Gutierrez, Silla-Martinez & Gabaldon, 2009). TThe best substitution model (GTR+F+R5) and ML analysis were performed with IQ-TREE 1.6.12 (Kalyaanamoorthy et al., 2017; Trifinopoulos et al., 2016), using the UFBoot2 method (Hoang et al., 2018) with 10,000 replicates.

Alpha diversity and microbial composition

The microbial diversity of each sample at both phylum and genus level of Bacteria and Archaea was estimated with the Hill numbers in terms of effective numbers of elements under the same sample coverage (Chao & Jost, 2012). Hill numbers are a mathematically unified family of diversity indices differing only by an exponent q (^qD) (Chiu & Chao, 2016), where q = 0 is equivalent to species richness, q = 1 corresponds to the exponential of Shannon entropy (effective number of common elements), and q = 2 is equivalent to the inverse of Simpson index, interpreted as the effective number of dominant elements (Ma & Li, 2018; Ma, Li & Gotelli, 2019). To compare the microbial diversity between the four samples, we used the 95% confidence intervals (CI), where no overlap between CI values indicates significant differences (Cumming, Fidler & Vaux, 2007). The qD diversity, sample coverage, and their respective confidence intervals were obtained with the iNEXT R package 3.5.3 (Hsieh et al., 2016), using as endpoint the maximum number of contigs in each sample and 1,000 bootstraps for the construction of the rarefaction curves and CI.

The difference in composition of bacterial and archaeal communities was assessed by one-way ANOSIM based on Bray-Curtis similarity metrics and with 10,000 permutations. When the ANOSIM was significant, a pairwise comparison was conducted between samples using the Bonferroni p-value correction as implemented in Past 4.2 (Hammer, Harper & Ryan, 2001). This analysis was conducted considering only those genera with relative abundance equal to or greater than 0.5%, meaning that rarely occurring microorganisms were not considered.

Predicted functional annotation

Gene annotation and coding sequence identification in feces microbiomes were performed using Prokka 1.12 (Seemann, 2014). One step for cleaning the protein fasta file to eliminate redundant sequences was performed using CD-HIT. We used blastp 2.10.0+ on the sequences that Prokka did not score (hypothetical proteins) against the UniProtKB/Swiss-Prot database. Additionally, non-redundant protein fasta files of each sample were obtained. We performed a second annotation of protein sequences with the online program GhostKOALA 2.2 (Kanehisa, Sato & Morishima, 2016), using GHOSTX search (it uses suffix arrays to find matching sequences and runs 100 times faster than BLAST), against a non-redundant set of KEGG genes (Kyoto Encyclopedia of Genes and Genomes).

Sequence and metagenome assembly

We extracted DNA from four fecal samples. According to the size of the collected feces, S1 (0.5 cm) corresponded to a young individual, S2, S3 and S4 (1.0 cm) corresponded to adult rabbits (Velázquez, Romero & López-Paniagua, 1996). We confirmed that the feces samples S1, S3, and S4 were from volcano rabbits because we found cytochrome b sequences in feces matched those reported previously from the same rabbit species (Osuna et al., 2020). Although cytochrome b sequence was not detected from sample S2, we are confident that it corresponds to volcano rabbit as feces of other rabbits inhabiting the zone are clearly distinctive (Velázquez, Romero & López-Paniagua, 1996). Maximum likelihood analysis indicated that the cytochrome b gene sequences found in the metagenomes were phylogenetically placed with the volcano rabbit (Fig. S1).

We obtained a total of 332.8 million raw reads from the four samples. Once the cleaning was completed, we rescued 316.6 million reads and the average number of reads per sample was 79.15 million (ranging from 69.6 to 96.6) (Table S1). Extra filtering of the clean sequences was conducted to obtain external contamination. The results demonstrated that between 98.2% and 99.6% of the clean reads were not assigned to adapters, human, rat, rabbit, yeast, mouse, pine, oak, fly, worm, and pika genomes. Individual assemblies of each sample S1, S2, S3 and S4 were performed. Fewer and larger contigs were obtained with metaSPADES than with MEGAHIT. Moreover, 83%–85% of the clean reads were mapped against each assembly.

Taxonomic assignment

The percentage of contigs that could be classified ranged between 14% and 18% for the taxonomic levels of phylum and genus of the Bacteria and Archaea domains. This could be due to the limitations of the non-human database where rare species are present (Tamames, Cobo-Simon & Puente-Sanchez, 2019). Twenty-nine phyla were obtained from the Bacteria domain (Fig. 1A). The phylum Proteobacteria was the most abundant in sample two (S2) and sample from the young rabbit (S1) that has been found in other young animals and humans (Moon et al., 2018). The most abundant phylum in samples S3 and S4 was Firmicutes. Additionally, we found three phyla of the Archaea domain: Euryarchaeota, Crenarchaeota and Thaumarchaeota with the most abundant being Euryarchaeota (Fig. 1B). Nine families were the most abundant in all samples (Table 1). The most abundant genus of the domain Bacteria in the young rabbit feces was Acinetobacter (Proteobacteria) while in adult-rabbit feces was Lachnoclostridium (Firmicutes) (Fig. 2A). Candidatus tachikawaea, Leclercia, Candidatus riesia, and Obesumbacterium were only found in the gut metagenome from young rabbit feces (Fig. 2A). The number of archaeal genera we encountered here in the fecal microbiome of the volcano rabbit is remarkable. We identified 15 genera of the domain Archaea, Methanosarcina, Methanoculleus and Methanococcus were the most abundant in all samples. The genus Halorabdus was the most abundant only in sample two (Fig. 2B), which is the only sample for which there was no confirmation of host rabbit identity.

We recovered 169 16S rRNA sequences (1,273 ∼1,500 bp) from the assembled metagenomes. Maximum likelihood analysis indicated that the 16S rRNA gene sequences found in the metagenomes were phylogenetically placed in the following phyla: Firmicutes, Bacteroidetes, Proteobacteria, Synergistetes, Tenericutes, Verrucomicrobia, Cyanobacteria, Lentisphaerae and Spirochaetes (Fig. S2). Besides, we found some 16S rRNA gene sequences that did not cluster with any reported phyla. Many 16S rRNA gene sequences were related to the species Marvinbryantia formatexigens, which has been found in the human gut and has been reported to degrade plant oligosaccharides (stachyose and raffinose) (Rey et al., 2010).

Alpha diversity and microbial composition

In all samples, bacterial and archaeal phyla and genera showed a sample coverage of 100%, indicating that sampling is complete for all samples (Fig. 3). Then the diversity comparisons were made directly according to their confidence intervals (Chao & Jost, 2012). For bacterial communities, the richness of S1 was significantly higher than in the other samples, at both phylum (Fig. 4A) and genus level (Fig. 4C). For the phyla of Archaea domain, the richness of both S1 and S2 were not significantly different while a similar pattern was observed between S3 and S4 (Fig. 4B). Richness estimates at the genus level in both domains differed significantly between all samples (Fig. 4C; Fig. 4D). The diversity values of the effective number of common and dominant elements are shown in Table S2.

According to ANOSIM, community composition at Bacteria and Archaea domains was not significantly different between samples (p >  0.05). For genera of domain Bacteria, the composition was significantly different between S1 and the other three samples (p = 0.0002), while for the genera of domain Archaea the composition was significantly different between all samples (p = 0.0002).

Functional annotation

A similar number of sequences was annotated from all samples (101,832–138,018) in all categories (Table 2). According to the function of their specific catalytic action, genes encoded oxidoreductases, transferases, hydrolases, lyases, isomerases, ligases, and translocases (Chang et al., 2021). The most abundant genes in all samples were transferases (28,523–21,324 sequences) and hydrolases (17,885–24,675) followed by oxidoreductases (10,797–7,757), ligases (5,081-7,013), lyases (4,896–6,064) and isomerases (4,227-5,500), no translocase-type enzymes were found with the PROKKA annotation (Fig. 5). The highest number of annotations of all the enzymes (82,572) was in S1 while the sample with the lowest number of annotated enzymes was S2 (61,201) (Fig. 5). Additional enzyme search using blastp showed between 954–1,125 total in different samples. Genes encoding translocases were recovered, which ranged from 16 to 37 sequences in the different samples. We found 21 functional categories. Protein Families (17,065–11,711 contigs per sample), genetic information processing (17,065–12,536), carbohydrate metabolism (glycolysis/gluconeogenesis, citrate cycle, TCA cycle; 12,571–9,776), protein families: signaling and cellular processes (7,400–5,197) and amino acid metabolism (7,103–5,730) were the most abundant in all samples. Meanwhile, xenobiotics biodegradation (127–58) and biosynthesis of secondary metabolites (49–66) had the lowest number of annotated sequences (Fig. 6).

Tannase enzyme sequences from different taxonomic groups were searched and downloaded from UniProtKB/Swiss-Prot and NCBI databases. Blastp search was performed to obtain the possible gene sequences for tannase enzymes. 18 sequences were found in the four metagenomic samples of volcano rabbits that were affiliated to phylum Verrucomicrobia, family Ruminococcaceae, Victrivallis vadensis, Blautia sp. and Clostridium sp. Finally, enzymes that participate in the degradation of cellulose and hemicellulose of grasses were searched. The highest number of sequences was obtained from young rabbit feces (849), while the lowest number was obtained in S2 (526) (Fig. 7). The most abundant were endo 1,4 β-xylanases for samples S1 (363) and S4 (248) and arabinofuranosidases for samples S2 (160) and S3 (149). The enzymes found in lower proportion were those encoding β-glucosidases in sample S2 (94) and endoglucanases in samples S3 (142) (Fig. 7).

Gut microbiota

There are few studies on the microbiota of wild Lagomorphs. To our knowledge, this is the first one to characterize the fecal microbiome of the volcano rabbit (R. diazi) that is endemic to central highlands in Mexico. As an herbivore, the volcano rabbit would feed on plants that contain endophytes which may become part of the gut microbiota. This occurs in other herbivorous animals’ endophytes that have tannases and enzymes to degrade plant cell wall components (Martinez-Romero et al., 2021). A critical issue that should be considered is whether the microbiota identified from herbivorous animals represent bona fide gut bacteria or reflect endophytes that are still contained within the plant tissue fragments macerated during DNA extraction. Here a mild procedure not including physical maceration was performed and all plant debris was eliminated by centrifugation.

Diet changes may have fast and important consequences in gut microbiota (Human Microbiome Project Consortium, 2012). The uncontrolled diets of R. diazi in their natural habitat may explain differences among adult samples. Otherwise, differences could be due to the distinct lapse from fecal excretion and collection for which we have no record. Peculiarities of the gut microbiota of babies, infants and young animals have been recurrently reported, and here we corroborated this with the only young rabbit that was sampled. There are reports from three-week-old Asian elephants (Loxodonta cyclotis) (Ilmberger et al., 2014), young bats (Leptonycteris yerbabuenae) (Gaona et al., 2019) and puppies and kittens (Moon et al., 2018) where Proteobacteria are abundant. Since the gut of newborns is rich in oxygen, the role of Proteobacteria may be involved in oxygen consumption, thus preparing the neonatal gut for colonization by strict anaerobic microorganisms necessary for a healthy gut (Shin, Whon & Bae, 2015).

Here we report that the fecal microbiome of the volcano rabbit was dominated by Proteobacteria, Firmicutes, Actinobacteria, for S1 and S2 and Firmicutes, Actinobacteria and Cyanobacteria for S3 and S4. We suggest that the similarity between S1 and S2 may be due to the individuals being related because young rabbits tend to be close to their mothers and these samples were collected at the same location. However, we remark that in zoos and in houses, rabbits would have an artificial diet. It is also worth noting that fecal and gut bacteria are not the same and more Proteobacteria were found in intestinal samples than in fecal samples in the European brown hare (Stalder et al., 2019). Additionally, Proteobacteria has been reported to be more abundant in rabbits with the symptoms of epizootic rabbit enteropathy (ERE) (Bauerl et al., 2014) and in other mammals, such as bats (Artibeus lituratus) (Ingala et al., 2018) and Asian elephants with a high fruit diet (Budd et al., 2020). On the other hand, a study in New Zealand white rabbits has reported a high abundance of Actinobacteria and Tenericutes (North, Dalle Zotte & Hoffman, 2019). Phylum Lentisphaerae was previously reported in hares (L. europaeus).

On the other hand, 16S rRNA reads from shotgun metagenomics provide sequences that are not subject to PCR primer bias and covers taxa that are not detected by primer sets (Yuan et al., 2015); therefore, we complemented the kraken 2 taxonomic assignment with a global 16S rRNA gene bacterial identification and phylogeny from metagenomes. We obtained only nine phyla in contrast to the 29 found with the kraken 2 analysis. Phylum Firmicutes was the most abundant in both. It has been reported that kraken 2 is suitable for the classification of gut microbiome with assemblies’ recruit of 85% or more of the original raw reads (Tamames, Cobo-Simon & Puente-Sanchez, 2019). Nevertheless, samples that have more diversity and complexity would not be correctly classified with kraken 2. Less studied habitat metagenomes such as non-human environments show fewer similarities with kraken 2 which is quite sensitive to the composition of databases and their performance decreases when rare species are present in the metagenome (Tamames, Cobo-Simon & Puente-Sanchez, 2019).

Domain Archaea has been reported in Caldes rabbit microbiome (Velasco-Galilea et al., 2018) with a single genus Methanobrevibacter (Euryarchaeota). Similarly, genus Methanobrevibacter was identified in a molecular profiling from a rabbit caecum (Kusar & Avgustin, 2010). These differences could be due to different strategies to characterize microbiota (16S rRNA amplicon sequencing vs shotgun). Here we found an unusual high number (15) of archaeal genera in the fecal microbiome of the volcano rabbit. Surprisingly, the most abundant genera were Methanosarcina, Methanoculleus, Methanococcus, Halorhabdus (Euryarchaeota) besides Methanobrevibacter. All of them have been found in the human gut microbiome (Borrel et al., 2020). Halorhabdus is usually associated with microaerophilic and halophilic environments. The axenic species described may use the monosaccharides glucose, fructose, and xylose for growth (Waino, Tindall & Ingvorsen, 2000). Xylose is produced mainly in plants (wood, fibers, fruits) and as a degradation product of hemicellulose. Except for Halorhabdus, the other archaea may perform the terminal step in the degradation of organic matter and produce energy and methane through methanogenesis pathways (Ferrer & Pérez, 2010; Ferry, 2020). These methanogenic archaea are generally associated with strict anaerobic environments, due to the high sensitivity to oxygen exposure of the enzymes involved in methanogenesis and to the lack of genes that code for proteins involved in protection against oxidative stress (Jasso-Chavez et al., 2015). Furthermore, these methanogenic archaea use compounds with one or two carbons to grow, for the generation of biosynthetic intermediates and polysaccharide storage (Santiago-Martinez et al., 2016). Methanoculleus, Methanococcus and Methanobrevibacter use carbon dioxide plus hydrogen (CO₂ + H₂) and formate (HCOO-), while Methanosarcina uses more diverse carbon sources, e.g., acetate (CH₃COO-), methanol (CH₃OH) and methylamines (CH₃NH₂), plus CO₂ + H₂ (Borrel et al., 2020; Buan & Robinson, 2018). These compounds are produced by the metabolism of eukaryotic and bacterial species. Specifically, there is evidence that the high intake of dietary fiber in animals, followed by fermentation mediated by the intestinal microbiota, changes the concentration and metabolism of methylamines (Li et al., 2017), which are a very good carbon source for methanogenic archaea. This suggests that the type of diet could increase the diversity and abundance of methanogenic species. The diversity and abundance of methanogenic archaea found in volcano rabbits may indicate a rich diversity of ingested vegetable and very complex interaction between bacteria and archaea. Seemingly there is no competition for food between these methanogenic and non-methanogenic archaea, due to the well-established difference in preference for carbon sources.

Prediction of functional annotation

In the present study, we found several functional categories suggesting that the fecal microbiome of R. diazi may have biosynthetic amino acid, carbohydrate metabolism, lipid metabolism, metabolism of terpenoids, polyketides and secondary metabolites and xenobiotic biodegradation capabilities. Enzymes found in the fecal microbiome of the volcano rabbit were mainly transferases, hydrolases, and oxidoreductases enzymes. Categories such as carbohydrate metabolism and amino acid metabolism had been reported in other herbivores (Muegge et al., 2011). Pathways such as bile secretion, mineral absorption, and xenobiotic biodegradation have been obtained in other rabbits (Zeng et al., 2015). Recently, metabolism of terpenoids –polyketides, amino acid metabolism, carbohydrate metabolism, nucleotide metabolism, lipid metabolism, energy metabolism, and metabolism of cofactors and vitamins pathways were reported in African elephants (Budd et al., 2020).

Enzymatic activity of tannases was identified in phyla Firmicutes and Proteobacteria: Enterobacter, Weissella and Lactobacillus (de las (de las Rivas et al., 2019), Bacillus, Streptococcus, Klebsiella, Enterococcus fecalis, Pantoea agglomerans, Staphylococcus lugdunensis, Lactobacillus plantarum, L. paraplantarum (Rodríguez-Durán et al., 2010) and Bacillus licheniformis (Palacio-Arango et al., 2018). In Koala feces, a Gram-negative, facultatively anaerobic and tannase-producing bacterium called Lonepinella koalarum was isolated (Osawa et al., 1995) which was not found in the microbiome of volcano rabbit feces pointing out the specialization or adaptation of some bacterial species to animal hosts.

Fiber degrading enzymes have been reported in bacteria such as Clostridium, Cellulomonas, Bacillus, Thermomonospora, Ruminococcus, Bacteroides, Erwinia and Streptomyces (Saratale, Saratale & Oh, 2012). Clostridium was abundant in S2, S3 and S4 whereas Erwinia in S1. In all samples, we also found phylum Actinobacteria, which is abundant in soil producing secondary metabolites, antibiotics and related with degrading plant cell wall components (Saratale, Saratale & Oh, 2012).