Multiomics characterisation of the zoo-housed gorilla gut microbiome reveals bacterial community compositions shifts, fungal cellulose-degrading, and archaeal methanogenic activity

Abstract

We carried out a comparative analysis between the bacterial microbiota composition of zoo-housed western lowland gorillas and their wild counterparts through 16S rRNA gene amplicon sequencing. In addition, we characterised the carbohydrate-active and methanogenic potential of the zoo-housed gorilla (ZHG) microbiome through shotgun metagenomics and RNA sequencing. The ZHG microbiota showed increased alpha diversity in terms of bacterial species richness and a distinct composition from that of the wild gorilla microbiota, including a loss of abundant fibre-degrading and hydrogenic Chloroflexi. Metagenomic analysis of the CAZyome indicated predominant oligosaccharide-degrading activity, while RNA sequencing revealed diverse cellulase and hemi-cellulase activities in the ZHG gut, contributing to a total of 268 identified carbohydrate-active enzymes. Metatranscriptome analysis revealed a substantial contribution of 38% of the transcripts from anaerobic fungi and archaea to the gorilla microbiome. This activity originates from cellulose-degrading and hydrogenic fungal species belonging to the class Neocallimastigomycetes, as well as from methylotrophic and hydrogenotrophic methanogenic archaea belonging to the classes Thermoplasmata and Methanobacteria, respectively. Our study shows the added value of RNA sequencing in a multiomics approach and highlights the contribution of eukaryotic and archaeal activities to the gut microbiome of gorillas.

Introduction

Gut microbes are necessary for mammals to digest complex polysaccharides in plants (fibres), as the microbial genomes encode the metabolic potential to support the hydrolysis and fermentation of plant material (Gomez, 2014). The gut microbiome has even allowed the evolution of herbivory as a dietary niche, and many mammals now rely completely on plant material as a dietary source (Moeller and Sanders, 2020). After millions of years of co-evolution, the gut microbiome is embedded not only in host metabolism but also in the immune and endocrine systems, with hosts and microbiomes shaping each other’s genomes (Moeller et al., 2016). The host microbiome can rapidly adapt to composition and function in response to dietary and other lifestyle changes. Although all mammals may have experienced such changes, the most frequent and drastic lifestyle transitions have occurred in human history, with industrialisation being the most recent and fast-paced example. Human exposure to microbes has been drastically altered in industrialised societies owing to antibiotics, higher sanitation levels, non-emergency (planned) C-sections, formula infant feeding, and processed foods, which are low in dietary fibre (Sonnenburg and Sonnenburg, 2019a, 2019b). Comparisons between industrialised microbiomes and those of hunter-gatherer societies indicate that, while the hunter-gatherer microbiome has a rich diversity of species and carbohydrate-active enzymes, the industrialised microbiome has a lower taxonomic diversity and is enriched with mucin-degrading enzymes. An increased abundance of mucin-degrading genes is a sign of low nutrient availability for the gut microbiome, allowing host mucin-degrading microbes to flourish. While the microbiome may have adapted to industrialisation, concerns have been raised about host biology because of the “asymmetric plasticity” between the rapidly adapting microbiome and the relatively stable mammalian genome. Sonnenburg and Sonnenburg (2019b) postulated that because of rapid lifestyle transitions, the microbiome can become incompatible with its host, possibly leading to the emergence of noncommunicable chronic diseases in industrialised societies.

Interestingly, the effects of industrialisation on the microbiome also apply to zoo-born primates. Indeed, zoo-housed primates are also exposed to higher sanitation levels, antibiotics, sometimes formula infant feeding, and, most consistently, a diet lower in plant complex carbohydrates than their wild counterparts. Efforts have been made to understand the effects of zoo housing on the mammalian microbiome to gain a general understanding of lifestyle effects, improve the health of zoo-housed animals, and support conservation efforts. Clayton et al. (2016, 2018) found that the microbiota of zoo-housed primates (doucs, Pygathrix nemaeus, and mantled howler monkeys, Alouatta palliata) was reduced in diversity and enriched in genera that dominated the microbiota of industrialised humans, such as Bacteroides and Prevotella. According to observations in North American zoos, the gut microbiomes of chimpanzees, bonobos, and gorillas contain 30% fewer genera than those of their wild fellow species, which is similar to the loss of genera in industrialised humans compared to humans living a hunter-gatherer lifestyle (Nishida and Ochman, 2021). In contrast, Campbell et al. (2020) showed that zoo-housed gorillas (ZHGs) had almost two times higher amplicon sequencing variant (ASV) richness than wild gorillas (WGs), while Shannon diversity values were similar between the two groups. A recent study by Narat et al. (2020) confirmed a significantly higher microbial richness in ZHGs than in WGs, although Shannon diversity values also increased. Despite this, the microbiota of zoo-housed apes is still most similar to that of wild apes and more similar to that of non-industrialised humans than Westernised humans. This is attributed mainly to the presence of Treponema species, which are rarely present in industrialised human microbiomes but are common members of the hunter-gatherer microbiome and zoo-housed and wild great apes (Campbell et al., 2020). However, species that are considered characteristic of the wild great ape microbiome are often lost in the zoo-housed microbiomes (Clayton et al., 2016; Campbell et al., 2020). The effects of lifestyle and diet on microbiome diversity and composition are evident in both humans and non-human primates. In gorillas, this is visible to the naked eye. The abdomens of WGs are bloated, a result of fermentation of plant material by the microbiome leading to gas production by methanogens. The abdomens of ZHGs often appear flat, indicating low levels of fermentation of plant material (Masi, 2011). WGs are herbivores that rely strongly on the gut microbiome to degrade plant material. The microbiota of western lowland gorillas is characterised by a high relative abundance of the phylum Chloroflexi compared to other great apes and humans and a high abundance of the phylum Spirochaetes, which is only present in high abundance among other primates in chimpanzees (Hicks et al., 2018).

Most data on the non-human primate microbiome originate from 16S rRNA gene amplicon sequencing and are thus compositional in nature. In human microbiome studies, shotgun metagenomics is now widely used for taxonomic and functional microbiome characterisation (Zhang et al., 2019). For gorillas, such data are still limited; previous studies utilised shotgun metagenomics to study the prevalence of antibiotic-resistance genes in non-human primate microbiomes (Campbell et al., 2020), or to show seasonal fluctuations in the abundance of functional pathways in the great ape microbiome (Hicks et al., 2018). Whereas metagenomic studies can identify the potential of carbohydrate-active enzymes in the gut microbiome, metatranscriptomics can reveal which genes are expressed under certain conditions (Zhang et al., 2019), and thus reveal the functional activity of the microbiome. To the best of our knowledge, metatranscriptomics has never been used to study the microbiome of gorillas. Here, we studied the compositional differences between the gut microbiota of wild and ZHGs by combining the 16 s rRNA gene amplicon sequencing data of the gorilla population in ARTIS (the Amsterdam Royal Zoo) with previously published 16S rRNA gene amplicon sequencing data from WGs and ZHGs (Campbell et al., 2020; Narat et al., 2020). Furthermore, we present a metagenomic and metatranscriptomic analysis of the microbiome of the dominant male (“silverback”) of the ARTIS gorilla population. We report the composition, fibre-degrading potential, and activity of the ZHG microbiome.

Methods

Diet of ARTIS ZHGs

The ARTIS gorillas were daily fed gorilla pellets (Marion Leaf Eater Food, Plymouth, MA), 0.5–1% of their body weight. They received a wide array of vegetables that differed each day (rotated weekdays), including celeriac, endive, fennel, tomatoes, chicory, carrots, and a handful of nuts, seeds, rice, and muesli per animal. Tree branches of the willow tree and, when available, herbs from the ARTIS garden were added daily. An egg or tofu was offered once or twice a week. Hay was always available in the enclosure.

Comparative analysis of zoo-housed gorilla and WG microbiota

We utilised previously published 16S rRNA gene V4 and V3–V4 amplicon sequencing data (Illumina technology) of WG and ZHG faecal samples obtained by Campbell et al. (2020) and Narat et al. (2020) for a comparative analysis in terms of diversity and composition with the microbiota data of ZHG faecal samples obtained in this study. In the study by Campbell et al. (2020), faecal samples from 28 WGs were collected in Nouabalé-Ndoki National Park in the Republic of the Congo in February and March 2013. This period corresponded to the transition period between the dry and wet seasons and was a period of lowered dietary diversity. In addition, faecal samples from 15 ZHGs were collected in 2017 in two zoos in the United States of America: 4 in the St. Louis Zoo (St. Louis, MO) and 11 in the Lincoln Park Zoo (Chicago, IL). Narat et al. (2020) collected 15 faecal WG samples from a forest site in southeastern Cameroon (between the towns of Yokadouma and Moloundou) in January 2016. This period corresponds to the dry season, which is characterised by low fruit availability and consumption of high-fibre fallback food. Additionally, faecal samples from the six ZHGs were collected from a European zoo site in November 2017. A more detailed description of these datasets can be found in Campbell et al. (2020) and Narat et al. (2020). The 16S rRNA gene amplicon sequence data, generated with an Illumina MiSeq instrument using paired-end 2 × 250 bp and 2 × 270 bp protocols, respectively, were retrieved from the Sequence Read Archive (accession numbers: PRJNA539933 and PRJNA666756).

Collection of faecal samples

A total of 15 faecal samples for 16S rRNA gene V3–V4 amplicon sequencing were collected from the ARTIS gorilla population between August 2019 and the end of January 2020. In July and August 2021, additional faecal samples from the silverback were collected in duplicates at two time points. Samples were collected within 30 min of defecation and frozen at −60°C. These four samples were subjected to shotgun metagenomics and RNA-seq analyses.

DNA and RNA extraction, library preparation, and illumina sequencing

DNA was extracted from gorilla faeces using the ZymoBIOMICS 96 MagBead DNA Kit (Zymo Research). DNA QC quantitation was performed with Quant-iT dsDNA Broad-Range Assay Kit (Invitrogen) and agarose gel electrophoresis for DNA integrity. The DNA served as a template for PCR amplification of a part of the 16S rRNA gene and subsequent sequencing using an Illumina MiSeq instrument. In short, amplicons from the V3–V4 regions of the 16S rRNA gene were generated by PCR using primers 341F and 785R, which generated a 444 bp amplicon (Klindworth et al., 2013) complemented with standard Illumina adapters. Index Primers were added to the amplicons of each sample through a second PCR cycle. The PCR products were purified using Agencourt AMPure XP (Beckman Coulter), and the DNA concentration was measured by fluorometric analysis (Quant-iT, Invitrogen). Subsequently, PCR amplicons were pooled using equimolar quantities, followed by sequencing on an Illumina MiSeq instrument using a paired-end 2 × 300 bp protocol and indexing. For samples subjected to shotgun metagenomics and RNA-seq, the ZymoBIOMICS ZymoD4300 kit was used to extract DNA and the ZymoBIOMICS kit Zymo R2040 to extract RNA. Library preparation was performed using Illumina RNA and DNA library preparation kits (20037135, 20020595, 20022371, and FC-131-1096). RNA integrity was assessed using a Bioanalyzer system (Agilent RNA 6000 Nano kit art5067-1511). RNA quality was ensured through the RIN value in combination with manual evaluation of the traces. Sequencing was performed on an Illumina NovaSeq 6000 instrument (paired-end, 2 × 150 bp) using an SP and S2 flow cell. Reads containing phiX control signals were removed using an in-house filtering protocol (Baseclear BV, Leiden, the Netherlands). In addition, reads containing (partial) adapters were clipped to a minimum read length of 50 bp. The second quality assessment was based on the remaining reads, using the FASTQC quality control tool version 0.11.8 (Andrews, 2010).

Analysis of 16S rRNA gene amplicon sequencing data

MiSeq paired-end sequence reads from Campbell et al. (2020) and Narat et al. (2020), and our own study were subsampled to 10 MBp, corresponding to 20.000, 18.518/18.519, and 16.628–16.821 reads, respectively (rarefaction curve, Supplementary Figure S1). For Narat et al. (2020) and our own study, the 16S rRNA gene V3–V4 amplicons were trimmed to retain only the V4 region, allowing for a consistent comparison with the samples of the Campbell et al. (2020) study. Hereafter, the paired-end sequences were collapsed into pseudoreads (reads that have been merged based on their overlapping sequences) using sequence overlap with USEARCH version 9.2 (Edgar and Bateman, 2010). Classification of these pseudoreads was based on the results of alignment with SNAP version 1.0.23 against the RDP database version 11.5 for bacterial organisms (Cole et al., 2014). OTU clusters were generated at 97% similarity with the USEARCH function “cluster_otus” based on pseudoreads. Reads that could not be assigned to sequences originating from the Domain of Bacteria were excluded. To analyse beta diversity, sequence counts were transformed to relative abundances, and Bray–Curtis distances were calculated for Principal Coordinate Analysis (PCoA). Principal Coordinate Analysis significance was tested using a non-parametric multivariate analysis of variance (ADONIS). Alpha diversity metrics were calculated from untransformed genus-level counts. The significance of differences in alpha diversity metrics between groups was determined using the Wilcoxon rank-sum test LEfSE (Segata et al., 2011). The differential abundance analysis between wild and zoo-housed gorilla samples was performed using default CPM normalization on taxa counts. Genera with an LDA score of at least 4 and a p-value below 0.05 were considered significantly enriched. The above-mentioned analyses and visualisations were carried out with R studio (R version 4.1.3) using the following packages: phyloseq 1.38.0 (McMurdie and Holmes, 2013), vegan 2.6–2 (Dixon, 2003), metacoder 0.3.5 (Foster et al., 2017), and microbiomeMarker 1.2.1 (Cao et al., 2022).

Taxonomic assignment of DNA and RNA reads

Kraken2 (2.1.1)/Bracken (2.6.1) (Lu et al., 2017; Wood et al., 2019) was used to estimate the taxonomic composition of the metagenomes and metatranscriptomes using the Kraken/Braken database, which is based on the RefSeq of the NCBI (feb-2022). We extended this database with the addition of six Neocallimastigomycetes genomes: Anaeromyces sp. S4 v1.0 (GenBank: GCA_002104895.1), and Neocallimastix sp. G1 v1.0 (GenBank: GCA_002104975.1), ASM1694683v1 (GenBank: GCA_016946835.1), Orpinomyces sp. strain C1A (GenBank: GCA_000412615.1), Piromyces sp. finnis v3.0 (GenBank: GenBank: GCA_002104945.1), and Piromyces sp. E2 v1.0 (GenBank: GCA_002157105.1). Assignments to the order of primates were filtered out, as these were considered to be contamination, and total sum scaling was applied.

Functional profiling of metagenome and metatranscriptome

Shotgun metagenomic reads were filtered for polyG tails, a known artefact of Illumina NextSeq two-colour dye technology (fastp version 0.23.1 (Chen et al., 2018) using the default -g option (polyG tail length = 10)). In addition, we used fastp to remove low-quality reads using the default -q option (bases with a phred quality >15 were considered qualified) and -u option (filter reads with more than 40% unqualified bases). Metagenome and metatranscriptome reads were filtered for host cDNA using STAR aligner version 2.7.9a (Dobin et al., 2013) and Kamilah_GGO_v0 reference genome (Refseq GCF_008122165.1). The RNA-seq reads were subjected to the same preprocessing, but additionally filtered for ribosomal RNA with sortmerna version 4.3.4 (Kopylova et al., 2012) by aligning reads against all SILVA databases of release 138.1, which contain small (16S/18S, SSU) and large (23S/28S, LSU) ribosomal subunits of Archaea, Bacteria, and Eukarya (Quast et al., 2013). Preprocessed DNA and RNA reads were annotated using HUMAnN3.0, version 3.0.0 (Beghini et al., 2021), with UniRef90 (version 20191b_full) as a reference database. To increase mapping rates, reads were assembled with MEGAHIT (1.2.9) (Li et al., 2015), genes were clustered using CD-HIT (4.8.1) (Fu et al., 2012), annotated with Prokka (v1.14.5) (Seemann, 2014), and EggNOG v. 5.0 (Huerta-Cepas et al., 2019) using the EggNOG-mapper (v2.0.0) (Cantalapiedra et al., 2021). Subsequently, the original preprocessed reads were realigned against the deduplicated contigs with bowtie2 v2.4.2 (Langmead and Salzberg, 2012) and subjected to a second analysis with HUMAnN3.0.0, using the contig KEGG and COG annotation as custom mapping files, to obtain relative abundances of these functional groups. Normalisation was performed using HUMAnN3’s normalisation function. Analyses were executed with Nextflow version 21.10.5 (build 5658) (di Tommaso et al., 2017) for reproducibility.

Annotation of carbohydrate-active enzymes

To quantify the relative abundances of carbohydrate-active enzyme (CAZyme) families, we ran HUMAnN3.0, on the read-contig alignments with a custom id-mapping file containing the corresponding CAZyme family for each contig. To create this file, clustered contigs were aligned to the dbCAN2 version of the CAZy database (version 10: CAZyDB.09242021). fa) using run_dbcan (v2.0.11) (Zhang et al., 2018), employing the option for alignment with HMMER in metagenome mode. Normalisation was performed using the recommended normalisation function of HUMAnN3.0.

Detection of fungal genes and transcripts

After preprocessing, ShortBRED (v0.9.4) (Kaminski et al., 2015) was used for sensitive detection of proteins of interest. We created markers for several proteins originating from Neocallimastigomycetes species downloaded from UniProt (release 2022_01). Preprocessed DNA and RNA reads were searched against several markers for proteins of interest. A first set of markers contained all Neocallimastigomycetes sequences in UniProt annotated with the subcellular localisation “hydrogenosomal.” A second set of markers consisted of all malic enzyme sequences (malate dehydrogenase [decarboxylating] EC 1.1.1.39) assigned to Neocallimastigomycetes and the latest sets of GH48 and GH6 proteins. Clustid was set to 85% and the minimum marker length to eight amino acids, based on the authors of ShortBRED’s own parameter optimisation, except for the marker set of malic enzymes, for which the clustid option was set to 100% because of high sequence similarity. A workflow overview of the shotgun metagenomics and RNA-seq data analysis is shown in Supplementary Figure S2.

Results

Beta diversity of bacterial communities differs between zoo-housed gorillas and WGs

In this study, we compared the composition of ZHG microbiota to that of WGs. To accomplish this, we combined previously published amplicon sequencing data (Campbell et al., 2020; Narat et al., 2020) from wild and zoo-housed western lowland gorillas with our amplicon sequencing data from 15 faecal samples from a group of western lowland gorillas in ARTIS, the Amsterdam Royal Zoo. Principal Coordinate Analysis based on the Bray–Curtis distance between samples indicated that bacterial community compositions from different studies clustered according to whether they originated from zoo-housed gorillas or WGs (Figure 1). ZHG samples from the three different studies converged and were clearly separated from WG samples. Non-parametric multivariate analysis of variance (ADONIS) calculated for the Bray–Curtis distances between ZHG and WG samples revealed that intra-group variability was lower than inter-group variability (p = 0.001), indicating that the microbiota composition of ZHGs was significantly different from that of WGs.

Figure 1. Principal Coordinate Analysis based on Bray–Curtis dissimilarities. The first two principal coordinates showed that, combined, they explain 53.6% of the observed variation in the zoo-housed and wild gorilla gut microbiota based on 16S rRNA gene V4 amplicon sequences. Colours indicate the study groups from which the samples originated. To aid interpretation, shapes indicate whether samples originated from wild or zoo-housed gorillas, regardless of the study.

Microbiota alpha diversity is elevated in ZHGs

Alpha diversity indices revealed a significant increase in genus-level bacterial diversity in ZHG samples compared to WG samples with respect to estimated richness (Chao1 index, p = 0.0001), abundance-based coverage estimator (ACE) index (p < 0.0001), and Fisher diversity index (p < 0.0001). No significant difference was found in the Shannon index between the ZHG and WG samples (p = 0.13) (Figure 2). Chao1, ACE, and Fisher’s exact test are known to be sensitive to species richness, whereas Shannon’s diversity index also reflects species evenness. Therefore, we concluded that microbiota richness was higher in the ZHG samples than in the WG samples.

Figure 2. Microbiota alpha diversity of wild and zoo-housed gorillas. Alpha diversity metrics were calculated at the genus level for the zoo-housed and wild gorilla gut microbiota, based on 16S rRNA gene V4 amplicon sequences. Alpha diversity metrics for wild gorillas (red) and zoo-housed gorillas (blue–green). The diversity indices included (A) Chao1, (B) abundance-based coverage estimator index (ACE), (C) Shannon, and (D) Fisher. The indicators of significance levels ****p < 0.0001 and ***p < 0.001 were computed using Wilcoxon rank-sum tests.

Characteristic members of the WG bacterial microbiota vanish in ZHGs

At the phylum level, the WG samples were, in order of abundance, dominated by Firmicutes, Chloroflexi, Bacteroidetes, Actinobacteria, Proteobacteria, and Spirochaetes (Supplementary Table S2). The most striking difference was the reduction in Chloroflexi, from 15.1% in WGs to 0.13% in ZHGs (Supplementary Figure S4). In ZHGs from ARTIS (this study) and those described by Narat et al. (2020), Chloroflexi could not be detected at all, whereas they were detectable at low levels in ZHGs, as described in the study by Campbell et al. (2020) (Supplementary Figure S4). At the genus level, the dominant taxa differed greatly between the ZHGs and WGs. An overview of the top 10 genera in both the ZHGs and WG samples and their relative abundances is provided in Supplementary Table S1; only the genera Prevotella and Clostridium were present in the top 10 of both groups. The per-sample Z-transformed relative abundances are shown in the heatmap in Figure 3.

Figure 3. Heatmap depicting the top 10 dominant bacterial genera in WG and ZHG. Relative abundances have been Z-transformed and from left to right the top 10 genera of WG and the top 10 genera of ZHG samples have been selected for display. As two of the bacterial genera were in the top 10 of both groups, this combined heatmap shows the abundance per sample for 18 genera.

Linear discriminant analysis effect size (LEfSe) was used to identify specific enriched taxa in the gut microbiota of wild and ZHGs (Figure 4). WG samples were significantly enriched for nine bacterial genera, including Flexilinea (Sun et al., 2016b), in line with the above-mentioned lack of detection of species from the phylum Chloroflexi in ZHG samples. This reduction in Chloroflexi in ZHG samples is a striking observation, since Chloroflexi are considered a characteristic species of the native gorilla microbiota that degrade carbohydrates in high-fibre fallback foods consumed during dry seasons (Hicks et al., 2018). A second characteristic taxon of the native gorilla microbiota that was reduced in ZHG samples was Olsenella (phylum Actinobacteria). However, the characteristic genus Treponema was preserved in ZHG samples and was found to be enriched. The oligosaccharide-consuming genera Lactobacillus and Lentimicrobium (Sun et al., 2016a) were also enriched in the ZHG samples.

Figure 4. Genera enriched in the microbiota of wild and ZHGs. Log10 transformed LDA scores and the corresponding log10 transformed p-values resulting from LEfSE analysis at the genus level between WG and ZHG samples. Genera with an LDA score > 4.0 and p-value <0.05 are shown.

Altered composition and substantial fungal and archaeal activity in the ZHG gut microbiota

We proceeded with metagenomics (MG) and metatranscriptomics (MT) analysis of four new samples from two time points of the dominant male of the ARTIS gorillas. We used Kraken/Braken2 to determine the taxonomic origin of metagenomic and metatranscriptomic reads. Bacteria accounted for 98.9% of the gorilla gut metagenome. Eukarya (0.07%), Archaea (0.4%), and viruses (0.01%) together accounted for approximately 1% of the remaining genetic material in the gorilla gut. Interestingly, in the gut metatranscriptome only 60.1% of transcripts were of bacterial origin. Eukaryotic transcripts constituted 20.1% of the metatranscriptome, and 19.3% of the transcripts originated from Archaea. The remaining 0.03% of the transcripts were attributed to viruses in MT samples. Taxonomic assignment of the MG and MT samples was in accordance with the 16S rRNA gene amplicon sequencing results of the ARTIS population. Firmicutes and Bacteroides dominate the bacterial microbiota. The major bacterial community composition shifts observed in the 16S rRNA gene amplicon sequencing results were confirmed by the MG and MT reads. Chloroflexi contributed only 0.13% of the MG reads and 0.03% of the MT reads according to the Kraken2/Bracken taxonomy assignment, whereas the high relative abundance of Lactobacillus was also reflected in the MT and MG samples: 10.5% of the MG reads were assigned to Lactobacillus, which accounted for 14.6% of the MT reads.

Almost all eukaryotic transcripts were assigned to the genera Piromyces, Pecoramyces, and Neocallimastix (Supplementary Table S4), and all three belonged to Neocallimastigomycetes, a class of anaerobic gut fungi previously identified in the rumen and gut of other herbivores (95% of eukaryotic RNA reads, 19.0% of the total metatranscriptome, and 0.07% of the total metagenome). The remaining 5% of eukaryotic transcripts (1.1% of the total metatranscriptome) are found in soil or plants, such as the pathogenic genus Melamspora, which infects Salix species (willow tree) that are consumed by ZHGs. Thus, these are not considered true members of the gut microbiota, but most likely originated from the passage of food through the gastrointestinal tract. Sequence reads matching the Archaea were present at rather low abundance in the metagenome (0.4%) but showed much higher abundance in the metatranscriptome (19.3%). They originate from methylotrophic methanogens (genera “Candidatus Methanoplasma,” 6.8, “Candidatus Methanomethylophilus,” 3.8%, Methanomassiliicoccus, 0.4%), and the more widespread hydrogenotrophic methanogens (genus Methanobrevibacter, 8.3%), as indicated in Supplementary Table S5. Metatranscriptome analysis revealed a substantial contribution of anaerobic fungi and archaea to the activity of the gut microbiome, accounting for 38.3% of the assignable metatranscriptome.

The gorilla microbiome harbours an extensive repertoire of carbohydrate-active enzymes

Shifting our perspective from a taxonomic to a functional one, we characterised the carbohydrate-active enzyme repertoire of the ZHG gut microbiome. We utilised the classification of the Carbohydrate-Active enzyme database (CAZyDB), which provides a sequence-based family classification linking the sequence to the specificity and 3D structure of the enzymes that assemble, modify, and break down oligo- and polysaccharides (Lombard et al., 2014). We identified genes and transcripts belonging to 268 CAZyme families. CAZyDB classifies enzymes into six superfamilies (Supplementary Table S6). Most CAZymes identified in the gorilla gut metagenome (MG) were glycoside hydrolases (GHs, 53.7%). The percentual contribution of each GH family to the total of GH families in the CAZy metagenome and metatranscriptome is shown in Figure 5. The percentages mentioned in the text below represent the percentual contribution of each family to the total CAZyome. Cellulases and hemicellulases were not found among the abundant GH families in the metagenome, accounting for more than 2% of the metagenome. GH13 (4.6%) is involved in starch breakdown. GH2 (3.2%), GH3 (3.2%), GH20 (1.9%), GH1 (1.8%), GH29 (1.6%), GH43 (1.5%), GH92 (1.5%), and GH97 (1.3%) contain oligosaccharide-degrading enzymes. GH78 (1.7%) and GH23 (1.7%) are debranching enzymes. However, in the metatranscriptome, cellulase and hemicellulases were observed within the abundant GH families (Figure 5). We refer to “abundant” as more than 2% of all GH families in the metatranscriptome. Families GH5 (3.5% MT, 0.98% MG), GH9 (2.1% MT, 0.3% MG), GH48 (2.1% MG, 0.01% MT), and GH6 (1.7% MT, not detected in MG) contain cellulases. GH11 (2.6% MG, 0.01% MT) and GH10 (1.4% MT, 0.3% MG) contain endo-hemicellulases. Abundant GH families in the metatranscriptome that are involved in the degradation of oligosaccharides are GH1 (4.6% MT, 1.8% MG) and GH18 (1.9% MT, 0.5% MG). GH57 (1.1% MT, 0.5% MG) and GH77 (3.1% MT, 0.9% MG) contain debranching enzymes.

Figure 5. Glycoside hydrolyase families in the gut metagenome and metatranscriptome of zoo-housed gorillas. Percentages represent the percentual contribution of each glycoside hydrolase (GH) family to the total number of GH families found by (A) shotgun metagenomics (DNA) and (B) metatranscriptomics or RNA-seq (mRNA) on faecal samples of a ZHG.

Mucin degradation by gut microbes requires the expression of GH33 (sialidases) to access mucin glycans (Glover et al., 2022). A total of 0.84% of the metagenomic CAZyome is accounted for by GH33, but no transcripts were observed. Thus, we observed no direct producers of mucin-degrading activity in the ZHG gut microbiome at the time of sampling. Other GH families known to be involved in mucin degradation are GH16, GH29, GH95, GH20, GH2, GH35, GH42, GH98, GH101, GH129, GH89, GH85, and GH84 (Glover et al., 2022); however, these are involved in the general degradation of oligosaccharides, and their presence cannot be considered direct evidence for mucin-degrading activity. Glycosyltransferases (GT) contribute 32.0% of the MG CAZyome and 20.8% of the MT CAZyome (but are not considered further here, since they are involved in polysaccharide synthesis instead of their breakdown). Carbohydrate esterase (CE) families aid in the deacetylation of plant polysaccharides (9.9% MT, 11.5% MG; Supplementary Figure S11). Carbohydrate-binding modules (CBM) target their substrates (mostly cellulose and hemicellulose) and promote prolonged interaction, thereby potentiating the enzymatic activities of other CAZymes. CBMs contributed 3.7% to MG CAZyome and 10.3% to MT CAZyome (Supplementary Figure S12).

Notably, the MT CAZyome did not predict the activity of CAZymes in the MG samples. No (hemi)cellulases were found among the abundant GH families in the MG samples, but the metatranscriptome indicated that they are among the most active families of the CAZyome. The most represented GH family in the metagenome, GH13 (starch degrading), contributed only 0.96% of the total MT CAZyome. Based on metagenome analysis alone, we would have underestimated the (hemi)cellulase activity of the zoo-housed microbiome, just as we would have underestimated the contribution of anaerobic fungi and archaea. The presence and activity of debranching enzymes and oligosaccharide-active enzymes allow for further consumption of degraded plant material by hemicellulases, or breakdown of oligosaccharides directly available in the diet. Interestingly, for families such as GH6 or GH48, none or extremely little MG reads were found, but significant levels of MT reads were detected. We further investigated two of the cellulase families for which this is the case, GH6 and GH48, with ShortBRED, a system for profiling protein families of interest at very high specificity in meta-omic sequencing data (Kaminski et al., 2015).

The cellulase families GH48 and GH6 originate from anaerobic gut fungi

The GH48 family contains cellulases, and its transcripts were detected in very low quantities in the metagenome. GH6 was found in MT samples but not in MG samples. From the 1,405 known GH48 DNA sequences, 110 corresponding protein sequences were deposited in UniProt, of which eight were detected using ShortBRED in our MT data (Supplementary Table S7) and none in the MG data. Within the GH48 family, we detected two transcripts encoding proteins of presumed fungal origin. The first of these transcripts resembles H2BPU6, a putative cellulase of Neocallimastix patriciarum that was previously detected in the rumen and functionally described by transcriptomic analysis (Li et al., 2014). The latter maps to Q8J1E3, a Cellulase Cel48A protein known from Piromyces sp. strain E2 isolated from the gut of an Indian elephant. Cellulase 48 is the major cellulosome component of Piromyces (Steenbakkers et al., 2002). Both species are members of the class Neocallimastigomycetes, which were previously shown to account for 19.0% of taxonomically assignable transcripts. The other GH48 transcripts found in ShortBRED were ascribed to an uncultured actinobacterium from the soil and uncultured bacteria from the rumen. Cellobiohydrolase from GH6 was first identified in Piromyces rhizinflatus.

Hydrogenosomal activity and symbiosis with methanogens

Neocallimastigomycetes possess hydrogenosomes instead of mitochondria (Czajkowski et al., 2020). The malic enzyme (ME) is considered a marker of hydrogenosomal activity (Mentel et al., 2008). Of the 16 malic enzyme proteins attributed to Neocallimastigomycetes in UniProt, transcripts from the ME of Neocallimastix californiae and Piromyces sp. (strain E2) were detected (A0A1Y2CXE5 in all samples and A0A1Y3NKG2 in three out of four samples). Other proteins involved in hydrogenosomal metabolism are succinate – CoA ligase subunit alpha (Q7Z941, originating from Neocallimastix patriciarum), with the beta subunit being detected in two out of four samples from another genome (P53587, Neocallimastix frontalis). Also expected to be in the hydrogenosome are homoacotinase, Aconitate hydratase, NADH dehydrogenase [ubiquinone] flavoprotein 1, Arg5,6 arginine biosynthetic enzyme, Alanine-tRNA ligase, Dynamin-type G domain-containing protein, and MSF1-domain-containing protein. An overview of the results is provided in Supplementary Table S8.

From the scientific literature it is known that the hydrogen produced in the hydrogenosome of anaerobic gut fungi is converted to methane by methanogenic archaea (Li et al., 2021). Methyl-coenzyme M reductase (MCR) catalyses the final step in methanogenesis and is considered a marker gene for methanogenesis (Chen et al., 2020). The alpha, beta, gamma, C, and D subunits of MCR were annotated as COG IDs (COG4058, COG4054, COG4057, COG4056, and COG4055, respectively). CPM values of the collective of these COG IDs were 14 to 145 times higher in the metatranscriptome than in the corresponding metagenome samples, indicating an overrepresentation of methanogenic activity in MT samples, in line with the Kraken2/Bracken taxonomy assignment. For methanogenesis, the abundance of genes (metabolic capacity) in the metagenome does not reflect its contribution to transcriptional activity. Our data indicated that anaerobic fungi and their symbiotic methanogens contribute substantially to enzymatic activity in the zoo-housed gut microbiome of gorillas (Supplementary Tables S4 and S5).

Discussion

In this study, we determined the gut microbiota composition of a population of zoo-housed western lowland gorillas using 16S rRNA gene V3–V4 amplicon sequencing. Beta diversity analysis showed that the ZHG microbiota had a significantly distinct composition compared to the WG microbiota. Their gut microbiota composition showed increased alpha diversity in terms of species richness compared to previously obtained data from wild western lowland gorillas, as indicated by the significantly lower Chao1 estimator, ACE, and Fisher diversity index values. The increased alpha diversity of the ZHG microbiome may be explained by the fact that all the WG samples were collected during periods of low dietary diversity. As observed in previous studies on primate zoo-housed microbiota, we found that the characteristic members of the gut microbiota of native gorillas vanished. Nine bacterial genera were significantly reduced in the ZHGs. The association of some of these genera with the consumption of high-fibre foods (Flexilinea and Olsenella) consumed during periods of low dietary diversity may explain this (Davenport et al., 2014; Hicks et al., 2018).

A striking observation was the relatively high abundance of Lactobacillus species, known to colonise the mammalian gut, in the microbiota of ZHGs, particularly in those from ARTIS. An anatomical explanation may be given here: gorillas share a gut anatomy that is similar to that of other monogastric herbivores. A typical type of gut anatomy is a simple stomach, in which the proximal part consists of non-glandular squamous cell epithelium. In horses, it has been shown that Lactobacillus species can colonise this area of the stomach and form biofilm-like structures (Yuki et al., 2000). Although this phenomenon has not been studied or observed previously in gorillas, based on the similarity in anatomy of the gorilla and stomach of monogastric herbivores, we hypothesise that the stomach of gorillas may be susceptible to colonisation by Lactobacillus species as well. Colonisation of the non-glandular part of the stomach is believed to be beneficial to the host because it prevents colonisation by pathogenic microorganisms (Yuki et al., 2000). This hypothesis of Lactobacillus colonisation of the gorilla stomach would clarify the presence of these Lactobacillus species in our compositional data, but does not explain why the relative abundances are so high compared to the compositional data of WGs. Again, diet most likely explains this observation. While the WG diet almost exclusively consists of high-fibre fruits and plant material, the ZHG diet inevitably contains high levels of oligosaccharides (Gänzle et al., 2012). They are mainly present in the upper part of the mammalian intestinal tract (Gänzle et al., 2012). The combined effect of easy adhesion to the upper intestinal epithelium and the high availability of oligosaccharides in this part of the gastrointestinal tract may explain why Lactobacillus species can flourish in the intestinal tract of zoo-house gorillas. However, based on the relative abundance alone, we cannot draw confident conclusions regarding the absolute abundance and distribution of bacteria in the gastrointestinal tract.

Annotation of metagenomic and metatranscriptomic CAZyome led to the identification of 268 distinct Carbohydrate-Active enzyme (CAZyme) families. None of the GH families abundant in the metagenome contained cellulases or hemicellulases. Instead, starch or other oligosaccharide-degrading and debranching enzymes are the most represented GH families in the metagenome. However, the metatranscriptome reveals the activity of multiple cellulases and hemicellulases, next to these oligosaccharide-degrading enzymes and debranching enzymes. From a methodological point of view, it is important to note that we only resorted to metagenomics for functional profiling and concluded that the ZHG microbiota hardly displayed potential for cellulose and hemicellulose degradation. However, metatranscriptomics revealed the activity of cellulases and hemicellulases, and that of Neocallimastigomycetes and methanogens in the gut microbiota of gorillas, which we would not have identified or largely underestimated by metagenomic analysis alone. A first possible explanation for this difference is that, due to the short half-life of mRNA, transcripts active in the upper intestinal tract may be less represented in the metatranscriptome than in the hindgut. Oligosaccharides are mainly digested in the upper intestinal tract, which may explain the overrepresentation of oligosaccharide-degrading enzymes in the metagenome compared with the metatranscriptome. A second explanation could be that some of the (hemi)cellulases are of eukaryotic origin. Interestingly, the activity of at least two families of cellulases (GH6 and GH48) can be attributed to members of the Neocallimastigomycetes, which are underrepresented (0.07%) in the metagenome but contribute 19.0% to the total metatranscriptome. Splicing effects may explain the lack of detection of eukaryotic genes in the metagenome, since methods for metagenomic analysis are tailored to the bacterial microbiome and thus do not account for introns and exons in assembly and contig annotation or lack reference genomes from eukaryotes in their databases (Lind and Pollard, 2021). We presume that this is the reason why no shotgun metagenomics studies have found or mentioned the presence of Neocallimastigomycetes in the gorilla microbiome (Hicks et al., 2018; Campbell et al., 2020).

In fact, only one previous study utilising ITS1 sequencing has discussed the presence of Neocallimastigomycetes in the gorilla gut (Schulz et al., 2018). The authors of this study recognised that anaerobic fungi are well-known members of the ruminant intestinal tract and are essential for the breakdown of plant material (Gruninger et al., 2014; Czajkowski et al., 2020). ITS1 sequencing results suggested that these fungi were indeed members of the gorilla gut microbiota (Schulz et al., 2018). Later studies utilised 16S rRNA gene amplicon sequencing (which does not detect eukaryotes) or shotgun metagenomics, and did not identify these anaerobic gut fungi. Similarly, in our shotgun metagenomics data, all analyses except for the taxonomic assignment with Kraken2/Bracken did not detect these anaerobic fungi at all. However, our RNA-seq results show that they are highly active in the gorilla gut and that they contribute to cellulase activity. Fungal isolates from the rumen microbiome degraded approximately 70% of dried leaf blades and stems in in vitro digestion studies (Akin et al., 1990). Thus, Neocallimastigomycetes may be major players in the degradation of plant cell wall material by the gorilla microbiome and are preserved in the zoo-housed microbiome.

The isolation of Neocallimastigomycetes and novel cellulases from the gorilla gut microbiome could be of great biotechnological interest for biomass degradation. Their polysaccharide-degrading efficiency is the highest in symbiosis with methanogenic archaea (Bauchop and Mountfort, 1981; Li et al., 2021), which is also overrepresented in the metatranscriptome compared to the metagenome. In human microbiome data, metagenome analysis has been found to underestimate methanogenesis when compared to metatranscriptome data (Franzosa et al., 2014, 2018). It should be noted that not all methanogenic activity found in our data originated from methanogens that are symbionts to anaerobic fungi, as these methanogens can also be symbionts to hydrogen-producing bacteria in the gut microbiota (Ma et al., 2020). Based on our analyses, we could not quantify bacterial and fungal contributions to hydrogen production. However, it has been found that anaerobic fungi and methanogens together degrade polysaccharides at a higher rate than a consortium of bacteria and methanogens (Ma et al., 2020), highlighting the biotechnological potential of Neocallimastigomycetes.

The methanogenic activity identified in our study in the gut microbiome of ZHGs originated from hydrotrophic and methylotrophic methanogens. We identified three genera of Thermoplasmata as the origin of methylotrophic methanogenesis and the genus Methanobrevibacter as the origin of hydrogenotrophic methanogenesis. Because methanogenesis from carbon dioxide consumes four molecules of hydrogen per molecule of methane and methanogenesis from methanol requires only one molecule, methyl-reducing methanogens should have an energetic advantage over hydrogenotrophic methanogens at low hydrogen partial pressures (Vanwonterghem et al., 2016). It remains to be investigated whether the relatively high percentage of methylotrophic methanogens is a specific characteristic of ZHGs and whether this is linked to a lower hydrogen pressure.

Our study shows how metatranscriptomics can uncover the contribution of eukaryotes and archaea to the herbivorous gut microbiome. It should be noted that we carried out metagenome and metatranscriptome analyses of four faecal samples from only one ZHG. Ideally, future research would include metatranscriptomic data of WG microbiomes, although we recognise that this is challenging because faecal samples must be collected as soon as possible after defecation due to the short lifetime of RNA.

In addition, shotgun metagenomic analysis workflows should be improved to detect eukaryotic and prokaryotic genes. In particular, the underrepresentation of gorilla-specific microbes in databases, including the HUMAnN3 database used in this study, poses limitations to the number of unassignable reads. One solution is to map all unassignable reads using BLAST analysis. However, this is computationally very intense if performed for millions of sequences. The prediction of eukaryotic and prokaryotic genes requires a different approach (i.e., to detect splicing sites in eukaryotic genes), which starts by correctly classifying metagenomic contigs. While some promising tools are being developed to address this, for instance, Whokaryote (Pronk and Medema, 2022), these need to be integrated into larger analysis frameworks.

In conclusion, the zoo-housed microbiome is altered in terms of diversity, composition, and expected function. However, hemi(cellulases), other carbohydrate-active enzymes targeting polysaccharides, and enzymes originating from Neocallimastigomycetes remain active in the zoo-housed microbiome. The abundance of fibre-degrading Chloroflexi was strongly reduced in the microbiota of ZHGs; however, the herbivore-characteristic activity of the microbiome was conserved. These results are promising considering conservation efforts and require further investigation to “rewild” the zoo-housed microbiome. Our study highlights the contribution of eukaryotes and archaea to the gorilla gut microbiome, thereby showing the added value of RNA sequencing and the importance of further improvements in metagenomics analysis pipelines for the detection of eukaryotic genes in shotgun metagenomics data.