Bacterial composition along the digestive tract of the Horned Screamer (Anhima cornuta), a tropical herbivorous bird

Sample collections

An adult male Horned Screamer, 2.5 Kg in weight, 79.0 cm total length and 210.0 cm wingspan, was hunted in San Pablo de Urama wetlands, Carabobo state, Venezuela (approximate coordinates: 10°31′14.7″N, 68°23′7.4″W) under the hunting permit for scientific purposes provided by the Ministerio del Poder Popular para el Ecosocialismo y Aguas (MINEC No 1727) and with approval of the IVIC Animal Ethics Committee (COBIANIM Dir-0884/1517). In our study area in Venezuela, we only found two individuals; and the second flew once we caught it. We returned some days after and did not see any other individual. As explained in the introduction, the Horned Screamer populations in Venezuela are low and primarily found in couples, hence using only one individual.

In the field, the animal was ventrally incised from throat to vent, to expose the gastrointestinal tract (GIT). After being removed from the body cavity, the GIT was dissected. The different organs (proventriculus, gizzard, small intestine, large intestine and cecum) were kept frozen in liquid nitrogen until arrival to the IVIC laboratory where the tissues were stored at −80 °C until DNA extraction.

In the IVIC lab, the GIT was defrosted and samples (200 mg) from tissues and content were sterile placed in the PowerBead tubes to perform the DNA extraction using the PowerSoil DNA Isolation kit (MO BIO, Inc., Carlsbad, CA, USA). The digestive tract was dissected, and samples from tissues and GIT content were taken from different organs and sections, including proventriculus, gizzard, small intestine, large intestine and cecum. The proventriculus and the small intestine were considered for this study as upper GIT, and the large intestine and cecum as lower GIT for analyses purposes (Fig. 1).

Sequence data processing and statistical analyses

DNA was extracted from these samples using the PowerSoil DNA Isolation kit (MO BIO, Inc., Carlsbad, CA, USA). DNA from the GIT content samples were normalized to 4 nM during 16S library prep and the hypervariable region V4 of the 16S ribosomal RNA gene (~291 bp) was amplified using the universal bacterial primers: 515F (5′GTGCCAGCMGCCGCGGTAA3′) and 806R (5′GGACTACHVGGGTWTCTAAT3′) in the Earth Microbiome Project (http://www.earthmicrobiome.org/emp-standard-protocols/16s/) (Caporaso et al., 2012) using previously reported conditions in projects from our group (Rodriguez-Barreras, Tosado-Rodriguez & Godoy-Vitorino, 2021; Ruiz Barrionuevo et al., 2022; Abarca et al., 2018). We used the Illumina MiSeq Reagent kit 2 × 250 bp to sequence the 16S amplicons. The 16S-rRNA reads were deposited in QIITA (Gonzalez et al., 2018) Bio project ID 12669, and the raw sequences are available in the European Nucleotide Archive ENA ID ERP135861.

Raw read pre-processing of the demultiplexed files in QIITA were processed with a Phred offset of 33, and default parameters. Forward reads were trimmed to 250 bp and a closed reference approach was selected for OTU picking using the SILVA reference database for taxonomy assignment with a minimum similarity threshold of 97% (Pruesse et al., 2007). The species table file in biom format was downloaded for downstream analyses using a locally run version of QIIME (Caporaso et al., 2010). Singletons and reads matching chloroplasts or mitochondria were removed from downstream analyses. Data were collected as previously described (Rodriguez-Barreras, Tosado-Rodriguez & Godoy-Vitorino, 2021; Ruiz Barrionuevo et al., 2022). Data analysis was done using a rarefaction level of 1,000 sequences per sample considering the gut sections: proventriculus n = 3, small intestine n = 1, large intestine n = 4 and cecum n = 2. Due to sequence quality and rarefaction, gizzard samples were removed and only 10 samples were analyzed. Data analyses was done by grouping samples into two main groups, upper GIT (proventriculus and small intestine n = 4) and lower GIT (large intestine and cecum n = 6). Community level analyses (beta diversity) were done by computing the pairwise Bray-Curtis distances between samples and plotted as 2D NMDS plot in microbiomeAnalyst (Dhariwal et al., 2017). Statistical significance between sample groups was assessed using the PERMANOVA test (Anderson, 2001).

Alpha diversity measures were estimated including richness (Observed, Chao 1), ACE Abundance-based coverage estimator (non-parametric) and Fisher diversity were plotted as boxplots and nonparametric statistical t-tests with Monte Carlo permutations were calculated. Barplots and heatmap revealing phyla and genus taxa were computed using the same parameters with MicrobiomeAnalyst (https://www.microbiomeanalyst.ca/). Linear Discriminant Analysis Effect Size (LEfSe) was done to identify significant features (p-value cutoff: 0.05; Log LDA Score: 2.0) discriminating lower and upper gastrointestinal tract samples (Segata et al., 2011). A supplementary file details the analyses (Supplemental Document).

Taxon-function analyses

Taxon-function attributions in the bird’s digestive and intestinal systems related to cellulose and hemicellulose fermentations were analyzed by the BURRITO framework (McNally et al., 2018). 16S rRNA OTU data were separated into two microbiome datasets: Upper GIT and Lower GIT. Attribution calculations were based on the original PICRUSt and a genomic content table characterized using the Greengenes OTU IDs. The association between taxonomy and function in both microbiome datasets (lower and upper GIT) were evaluated using their Greengenes OTU IDs and a custom KEGG BRITE hierarchical tree, which was composed using KO attributed to cellulose and hemicellulose metabolism, including starch and sucrose metabolism, galactose metabolism, glycolysis, pentose phosphate pathway, pyruvate, butanoate and propanoate metabolism. The summary levels of functional attributions for each taxon to each metabolic pathway involved in the fermentation of plant carbohydrates and average shares of each function attributed to each taxon were presented as a heatmap constructed using GraphPad Prism version 9 (GraphPad Software, San Diego, CA, USA). Metabolic network of the fermentation of plant carbohydrates were composed using KEGG BRITE hierarchical tree (Kanehisa & Goto, 2000) and identified KO functional data attributed to the fermentation of plant carbohydrates. The network was constructed using Cytoscape (Shannon et al., 2003).

Bacterial composition between upper and lower GIT

A total of 123,022 good quality 16S rRNA gene reads were used for downstream analyses from the 10 gut samples of the Horned Screamer. Data analyses were done using a rarefaction level of 1,000 reads.

Beta diversity analyses did not show significant differences in community structure or composition (Permanova p-value = 0.289), however the NMDS plot shows a distinction of upper GIT samples to the right of the plot x-axis, while most lower GIT samples converge to the left of the NMDS plot (Fig. 2A). Alpha diversity analyses among sites, showed nearly significant differences among GIT sites (Fig. 2B). The lower GIT samples had nearly significant higher alpha diversity (KW p-value: 0.064).

In terms of composition, we found seven different phyla in the Screamer’s GIT. The most abundant phyla were Firmicutes and Bacteroidetes. Other phyla like Proteobacteria, Cyanobacteria, Euryarchaeota, Verrucomicrobia and Actinobacteria were found in lower abundance (Fig. 3A). At a genus level, a wider range of bacteria was found with uncultured species predominant. Some of the genera include Prevotella, Anaerostipes, or Bacteroides (Figs. 3B and 4).

LEfSe analyses determined that at the phyla level, Proteobacteria most likely to explain differences between upper and lower GIT (p-value = 0.019). Bacteroidetes and Actinobacteria were more dominant in the upper GIT although there were no significant differences (Fig. 5A). At a genus level, Helicobacter, Vibrio, Enterobacter and Acinetobacter (from the Proteobacteria phylum), had higher significant abundance in the upper GIT (Fig. 5B). Blautia, Oribacterium, Roseburia, Ruminococcus, Intestinimonas, Marvinbryantia (Firmicutes), Desulfovibrio (Proteobacteria) and Parabacteroides (Bacteroidetes) are more dominant in the lower GIT (Fig. 5B).

Taxon-function attributions to metabolic pathways involved in the fermentation of plant carbohydrates

The microbiome functional capacity in the fermentation of plant carbohydrates was evaluated in the upper and lower GIT via the BURRITO framework (McNally et al., 2018). We identified 22 functions associated with starch and sucrose metabolism, galactose metabolism, glycolysis, pentose phosphate pathway, pyruvate, butanoate, and propanoate metabolism (Fig. 6A), and 28 species contributors. In both regions, Firmicutes was found as the main drivers of the activity of metabolic enzymes involved in the fermentation of plant carbohydrates. Bacteroidetes and Proteobacteria mainly contributed to many metabolic reactions in the upper GIT (Figs. 6B and 6C).

In the lower GIT, only two genera of Bacteroidetes, such as Bacteroides, contributed to galactose metabolism and Rikenellaceae—to galactose metabolism and butanoate metabolism. On the other hand, among Proteobacteria identified in the lower GIT, Desulfovibrionaceae contributed to butanoate metabolism and propanoate metabolism correspondingly (Figs. 6B and 6C). Actinobacteria, including Corynebacterium, Micrococcus, and Propionibacterium, contributed to pyruvate metabolism, galactose metabolism, butanoate metabolism, and propanoate metabolism only in the upper GIT (Figs. 6B and 6C). In contrast, Methanomassiliicoccaceae_vadinCA11 (Euryarchaeota) contributed with many metabolic reactions only in the lower GIT (Figs. 6B and 6C). Cerasicoccaceae (Verrucomicrobia) and Cyanobacteria_YS2 also contributed to many metabolic reactions in both anatomic regions (Figs. 6B and 6C).

Metabolic functions associated with gut microbiome on the fermentation of plant carbohydrates

It is recognized that the central metabolic hub in the process of fermentation of carbohydrates is glycolysis. Cellulose (CL), a monomeric polymer of glucose, was converted into cellodextrin and cellobiose by endoglucanase (K01179) and cellulose 1,4-beta-cellobiosidase (K01225), followed by breaking cellobiose into glucose by beta-glucosidase (K01188) in all regions (Fig. 6C, i). Identification of possible taxonomic drivers of hemicellulose (hCL) fermentation, a more complex heterogeneous polysaccharide, did not produce any conclusive results with some exceptions, such as the microbial involvement in the conversion of galactose and xylose in both regions. In particular, the digestion of galactose to glucose was associated with alpha-glucosidase activities (K01187), while the digestion of xylose to fructose-6-phosphate, erythrose-6-phosphate, and xylulose-5-phosphate through the pentose phosphate pathway was concomitant with the activity of transketolase (K00615) (Fig. 6C, i).

We also analyzed the functional capacity of microbiota in yielding end-products acetate, propanoate, butanoate, and acetoacetate (Fig. 6C, ii). Thus, propanoate produced from propanoyl-CoA by acetyl-CoA synthetase (K01895) was predicted for both regions and associated with Firmicutes, Bacteroidetes, Proteobacteria, Cyanobacteria, Euryarchaeota, and Verrucomicrobia. In the upper GIT, Prevotella and Ruminococcus, and in the lower GIT Ruminococcus were the most robust contributors for this metabolic reaction (Fig. 6C, ii). The contribution of Actinobacteria for propanoate production was predicted only in the upper GIT (Fig. 6C, ii). Contribution to butanoate and acetoacetate production in the upper GIT was associated with Firmicutes, Bacteroidetes, Proteobacteria, and Verrucomicrobia, where the strongest contributors were Ruminococcus, Peptostreptococcaceae, and Bacteroidales_S24-7 (Fig. 6C, ii). In the lower GIT, butanoate production was concomitant with Firmicutes, Proteobacteria, Cyanobacteria and Euryarchaeota, and Ruminococcus was the most robust contributor to this metabolic reaction (Fig. 6C, ii). Next, we identified that the production of acetate by acetyl-CoA hydrolase (K01067) in both queried regions was associated with Firmicutes, Cyanobacteria, and Verrucomicrobia. In addition, Actinobacteria and Proteobacteria also contributed to this metabolic reaction in the upper GIT. Examining the shares in acetate production revealed that Coprococcus was the strongest contributor in both regions, while Oscillospira in the upper GIT and Phascolarctobacterium, Cerasicoccaceae Methanomassiliicoccaceae_vadinCA11 in the lower GIT (Fig. 6C, ii). Conversion of fructose-6-phosphate to glyceraldehyde-3-phosphate and 3-phosphoglycerate were associated with the activity of transaldolase (K00616) and glyceraldehyde 3-phosphate dehydrogenase (K00134) and fructose-bisphosphate aldolase (K01623) (Figs. 6C, i and iii). Several glycolytic enzymes such as L-lactate dehydrogenase (K00016), 2,3-bisphosphoglycerate-dependent phosphoglycerate mutase (K01834), enolase (K01689), glucose-6-phosphate isomerase (K01810), pyruvate kinase (K00873), and phosphoglycerate kinase (K00927) were predicted to be associated with either Firmicutes or Bacteroidetes (Fig. 6C, iii).