Anaerobic fungi in the tortoise alimentary tract illuminate early stages of host-fungal symbiosis and Neocallimastigomycota evolution

Abstract

Anaerobic gut fungi (AGF, Neocallimastigomycota) reside in the alimentary tract of herbivores. While their presence in mammals is well documented, evidence for their occurrence in non-mammalian hosts is currently sparse. Culture-independent surveys of AGF in tortoises identified a unique community, with three novel deep-branching genera representing >90% of sequences in most samples. Representatives of all genera were successfully isolated under strict anaerobic conditions. Transcriptomics-enabled phylogenomic and molecular dating analyses indicated an ancient, deep-branching position in the AGF tree for these genera, with an evolutionary divergence time estimate of 104-112 million years ago (Mya). Such estimates push the establishment of animal-Neocallimastigomycota symbiosis from the late to the early Cretaceous. Further, tortoise-associated isolates (T-AGF) exhibited limited capacity for plant polysaccharides metabolism and lacked genes encoding several carbohydrate-active enzyme (CAZyme) families. Finally, we demonstrate that the observed curtailed degradation capacities and reduced CAZyme repertoire is driven by the paucity of horizontal gene transfer (HGT) in T-AGF genomes, compared to their mammalian counterparts. This reduced capacity was reflected in an altered cellulosomal production capacity in T-AGF. Our findings provide insights into the phylogenetic diversity, ecological distribution, evolutionary history, evolution of fungal-host nutritional symbiosis, and dynamics of genes acquisition in Neocallimastigomycota.

Fig. 1. AGF diversity and community structure in Tortoises.

A Community composition in the 11 tortoise fecal samples studied. The tortoise phylogenetic tree was downloaded from timetree.org. The pie chart to the right shows the total percentage abundance of the three tortoise-affiliated genera (NY54, NY36, and NY56) (green) versus other AGF genera (peach). AGF community composition for each tortoise sample is shown to the right as colored bars corresponding to the legend key. B Percentage occurrence (left) and percentage abundance (right) of the three tortoise-affiliated genera in previously studied cattle, white-tail deer, goats, horses, sheep, and other mammals, as well as in the 11 tortoise samples studied. The number of individuals belonging to each animal species is shown on the X-axis. Color code follows the key in (A). C Maximum likelihood phylogenetic tree constructed from the alignment of the D1/D2 region of the LSU rRNA genes and highlighting the position of the three tortoise-affiliated genera in relation to all previously reported cultured and uncultured AGF genera. Genera are color-coded by family or putative family, and the three tortoise-affiliated genera are shown in green boldface. D Distribution of sequence divergence within each genus. Data are derived from 129284, 7525594, and 586689 distances for NY36, NY54, and NY56 respectively. E AGF load (determined using qPCR and expressed as copy number/g fecal sample) in the 11 tortoise samples studied here in comparison to ten individual cattle, goats, sheep, and horses selected. Significance is shown above the boxplots and corresponds to the two-sided Student t-test p value. Source data for Fig. 1B–D are provided as a Source Data file. Boxplots in (B, D, E) extend from the first to the third quartile and the median is shown as a thick line in the middle. The whiskers extending on both ends represent variability outside the quartiles and are calculated as follows: Minimum whisker = minimum quartile − 1.5 × inter-quartile range; Maximum whisker = maximum quartile + 1.5 × inter-quartile range. All points outside the box and whiskers are outliers.

Fig. 2. Patterns of AGF alpha and beta diversity in the 11 tortoise samples studied in comparison to a subset of mammalian hosts previously studied (Dataset S2A).

A Box and whisker plots showing the distribution of 4 alpha diversity measures (observed number of genera (Sobs), Shannon, Simpson, and Inverse Simpson) for the different animal species. Results of two-sided Wilcoxon signed rank test for pairwise comparison of tortoise (pink) alpha diversity indices to mammals (cyan; cattle (n = 25), deer (n = 24), goats (n = 25), horses (n = 25), and sheep (n = 25)) are shown above the boxplots. * Denotes p values of 0.01 < p < 0.05 and ** denotes 0.001 < p < 0. Boxplots in Fig. 1B, D, E extend from the first to the third quartile and the median is shown as a thick line in the middle. The whiskers extending on both ends represent variability outside the quartiles and are calculated as follows: Minimum whisker = minimum quartile − 1.5 × inter-quartile range; Maximum whisker = maximum quartile + 1.5 × inter-quartile range. All points outside the box and whiskers are outliers. B, C Principal coordinate analysis (PCoA) plot based on the phylogenetic similarity-based index weighted Unifrac. The percentage variance explained by the first two axes is displayed on the axes, and ellipses encompassing 95% of variance are displayed. Samples and ellipses are color-coded by host class (B), and host species (C). Some of the circles representing tortoise samples might not be apparent due to overlap with other data points. D Double principal coordinate analysis (DPCoA) biplot based on the phylogenetic similarity-based index weighted Unifrac. The percentage variance explained by the first two axes is displayed on the axes, and ellipses encompassing 95% of variance are displayed. Samples and ellipses are color-coded by host species. AGF genera are shown as black empty circles, and the three tortoise-affiliated genera are labeled. Source data are provided as a Source Data file.

Fig. 3. Bayesian phylogenomic maximum clade credibility (MCC) tree of Neocallimastigomycota with estimated divergence time.

The isolate names are color-coded by host class as shown in the legend. Strains belonging to the two T-AGF genera are shown in boldface, and the taxa label is shown to the right. All clades above the rank of the genus are fully supported by Bayesian posterior probabilities. The 95% highest probability density (HPD) ranges (blue bars) are denoted on the nodes, and the average divergence times are shown. Geological timescale is shown below.

Fig. 4. CAZyome composition difference between tortoise-sourced and mammalian-sourced (n = 54) strains.

A Box and whisker plots for the distribution of the total number of GHs (top), CEs (bottom left), and PLs (bottom right) identified in the transcriptomes (mammalian sourced, cyan; tortoise sourced, pink). Only CAZy families with >100 total hits in the entire dataset are shown, and CAZy families that were significantly more abundant in mammalian versus tortoise transcriptomes are shown in red text. Boxplots in Fig. 1B, D, E extend from the first to the third quartile and the median is shown as a thick line in the middle. The whiskers extending on both ends represent variability outside the quartiles and are calculated as follows: Minimum whisker = minimum quartile − 1.5 × inter-quartile range; Maximum whisker = maximum quartile + 1.5 x inter-quartile range. All points outside the box and whiskers are outliers. Two-sided Wilcoxon test adjusted p values for the significance of difference in CAZyome composition for the families in red text are shown to the right, along with the values for GH5, GH13, GH16, GH43, and PL1 sub-families. B Principal coordinate analysis (PCoA) biplot based on the GH families composition in the studied transcriptomes. The % variance explained by the first two axes is displayed on the axes, and strains are color-coded by AGF genus, as shown in the figure legend to the right, while GH families are shown as smaller cyan spheres with black borders. Source data are provided as a Source Data file.

Fig. 5. Comparative cellulosomal analysis between representatives of the two tortoise-affiliated genera (genus Astrotestudinimyces, strain B1.1; and genus Testudinimyces, strain T130A) and one mammalian affiliated strain (Orpinomyces joyonii strain AB3).

A Maximum likelihood mid-point rooted phylogenetic tree showing the relationship between scaffoldin ScaA protein homologs identified in Orpinomyces joyonii strain AB3 (12 copies denoted AB3_1 through AB3_12 and shown in purple text), Astrotestudinimyces strain B1.1 (2 copies denoted B1.1_and B1.1_2 and shown in brown boldface text), and Testudinimyces strain T130A (2 copies denoted T130A_and T130A_2 and shown in orange boldface text) in comparison to a reference set of 319 Neocallimastigomycota ScaA homologs retrieved from Uniprot. All reference ScaA homologs are shown with their Uniprot ID followed by the AGF strain name color-coded by genus, as shown in the legend. B Comparison of the percentage distribution of functions (as predicted by NCBI Conserved Domain database) encoded by cellulosomal peptides (all predicted peptides harboring a non-catalytic dockerin domain in the two tortoise affiliated genera (genus Astrotestudinimyces, strain B1.1; and genus Testudinimyces, strain T130A) and the mammalian affiliated strain (Orpinomyces joyonii strain AB3) and destined to the extracellular milieu (as predicted by DeepLoc)). The total number of peptides is shown above each column. C CAZyome composition of the predicted cellulosome in the three strains compared. Source data for Fig. 5B, C are provided as a Source Data file.