The microbial load and alpha diversity of the gut microbiome increase after fungus acquisition
Microbial load, quantified with ddPCR of 16S rRNA gene copies, was significantly affected by timepoint (ANOVA: F2,80 = 4.632, p = 0.013, Cohen’s F = 0.34) which furthermore was dependent on sample type, evident by a significant interaction term (F3,80 = 17.79, p < 0.001, Cohen’s F = 0.82). Pairwise examination concluded that while both Post and Field samples had significantly higher loads than Pre, Post and Field where indistinguishable in workers (TukeyHSD - Pre vs Post: diff = -1.680, p = 0.003; Pre vs Field: diff = -2.231, p < 0.001; Post vs Field: diff = -0.551, p = 0.938), yet only Pre and Field were distinguishable for soldiers (Pre vs Post: diff = -2.325, p = 0.230; Pre vs Field: diff = -3.101, p = 0.0288; Post vs Field: diff = − 0.777, p = 0.995) (Fig. 2A). However, while colony maturity in general increased microbial load for workers and soldiers, it decreased the microbial load of fungus comb (Post vs Field: diff = 5.371, p < 0.001) (Fig. 2A). Timepoint also had a significant effect on Chao1 richness (ANOVA: F2,84 = 185.7, p < 0.001, Cohen’s F = 2.10) that was dependent on caste (F2,84 = 19.75, p < 0.001, Cohen’s F = 0.84). Richness increased across all timepoints in workers (Pre vs Post: diff = -232.840, p < 0.001, Pre vs Field: diff = -331.452, p < 0.001, Post vs Field: diff = -98.612, p < 0.001) but only between Pre and Post/Field in soldiers (Pre vs Post: diff = -289.877, p < 0.001, Pre vs Field: diff = -225.300 p < 0.001, Post vs Field: diff = 64.577, p = 0.743) (Fig. 2B). As for microbial load, Chao1 richness decreased from Post to Field fungus combs (Post vs Field: diff = 150.723, p < 0.001) (Fig. 2B).
Multivariate analysis of the metabarcoding data using PERMANOVA allowed us to infer that 24.2% of the variation in community compositions could be explained by timepoint (F2,84 = 16.84, R2 = 0.2427, p < 0.001) (Fig. 3A). Pairwise comparisons further demonstrated that all timepoints differed from each other but that Post and Field were more similar (Pre vs Post: F2,84 = 15.49, R2 = 0.1999, p < 0.001: Pre vs Field: F2,84 = 17.98, R2 = 0.2645, p < 0.001; Post vs Field: F2,84 = 10.35, R2 = 0.1355, p < 0.001).
Metabolic potential changes associated with the acquisition of the fungal cultivar
To infer the metabolic potential of the fungus-farming termite gut microbiome with regard to carbohydrate metabolism, we investigated the catalogue of enzymes involved in the breakdown of complex carbohydrates from the CAZy database using dbCAN4 (v.4.0.0). We identified a total of 292 unique CAZyme families across timepoints and symbiotic partners. The number of CAZyme families in the gut microbiome increased by approximately 1/3 with the acquisition the fungal cultivar, from 204 (Pre) to 270 (Post) or 277 (Field). This clarified that fungus acquisition did alter the potential function of the microbiome, although most were maintained with colony founding (Fig. 3B). The complementarity observed between termite, fungal, and bacterial CAZyme contributions is consistent with previous findings22, where there is also some overlap in carbohydrate metabolism between symbiotic partners. The fungus and termite shared respectively 35.3% and 18.5% of CAZyme families identified across time points.
Given the clear expansion in a portion of the carbohydrate active functions of the gut microbiome, we further examined the relative abundances of CAZyme families. Heatmaps of the logged relative abundance profiles of CAZyme families revealed six distinct groupings by timepoint (Fig. 3C). Families in group one were present at all timepoints in relatively high abundance. Group two was variably present in Pre but always present in Post and Field colonies, while group three was absent from Pre but almost always present in Post and Field. Group four was also absent from Pre but only variably present in Post and Field, suggesting functional maturation of metabolic potential in the microbiome post-fungus acquisition. Group five families were, as group one, consistently present across timepoints but in relatively lower abundances, and finally group six was always present in Pre but only variably present in Post and Field, potentially indicating loss of function during microbiome maturation. Therefore, a core set of functions appear to always be maintained through colony development, making up the majority of carbohydrate metabolism in the gut, while horizontal transmission of the fungus drives a general increase in abundance of other putative functions.
The clear groupings of families, based on relative abundances, were as a whole consistent with functional shifts indicated by substrates predictions from dbCAN4 (Fig. 3D). Despite evident shifts, when exploring the range of substrates that the predicted enzymes putatively target, we saw a lack of enzymes targeting sialic acid, fructose, and gellan in the Pre timepoint, whilst Post and Field variably lacked enzymes targeting xanthan and chitooligosaccharides (Fig. 3D). Notably, enzymes targeting lignin reduced in relative abundance after fungal acquisition (TukeyHSD of ANOVA - Pre vs Post: diff = 9.221x10− 05, p < 0.001; Pre vs Field: diff = 1.085x10− 04, p < 0.001; Post vs Field: diff = 1.630x10− 05, p < 0.639). The absence of predicted substrate from many of the CAZymes meant that we could only establish a partial picture of how acquisition of the fungal cultivar affects substrate targets. However, it was evident that enzymes targeting the majority of plant (e.g., cellulose, pectin, hemicellulose, and lignin) and fungal (e.g., α-glucan, β-glucan, chitin, and mannan) cell wall components were abundant at all timepoints suggesting ample ability to aid in their breakdown from colony inception.
We then compared Shannon diversity and Observed richness of CAZyme potential across timepoints revealing that time had a significant effect on both CAZyme diversity and richness with Pre being substantially lower for both (ANOVA – Shannon: F2 − 16 = 76.51, Cohen’s F = 3.09, p < 0.001; Observed: F2 − 16 = 1281, Cohen’s F = 12.66, p < 0.001) (Fig S1A). Matriline also had a significant effect on richness but not Shannon diversity; however, to a lesser extent (Shannon: F2 − 16 = 1.660, p = 0.221; Observed: F2 − 16 = 79.34, Cohen’s F = 3.15, p < 0.001). Pairwise comparisons of Shannon diversity and Observed richness with TukeyHSD showed Pre to be lower than both Post and Field, whilst the latter two did not exhibit significantly different Shannon diversity despite significantly different Observed richness (Fig S1, see Table S5 for details). This suggests that acquisition of the fungus not only shifts the microbiome to a greater diversity of bacteria, as was also evident from the ASV analysis, but also in encoded CAZymes.
Analysing the Bray Curtis distances between CAZyme relative abundance profiles through multivariate compositional analysis revealed that timepoint had a significant effect and explained a substantial component of the observed variation (PERMANOVA: F2,16 = 13.20, R2 = 0.4677, p < 0.001) whilst matriline was not (F2,33 = 0.9839, p = 0.4039). The pairwise comparisons further revealed that Pre stood out as particularly distinct, as indicated by larger R2 values (Pre vs Post: F1,14 = 22.62, R2 = 0.5456, p = 0.001: Pre vs Field: F1,10 = 8.739, R2 = 0.3201, p = 0.001; Post vs Field: F1,8 = 7.592, R2 = 0.2340, p = 0.035). Matriline remained insignificant in all comparisons (Table S6). These clear groupings by timepoint were also evident from PCoA clustering that further illustrated that the reduced explanatory power of time in comparisons of Pre and Field was likely two samples that substantially deviated from the rest of the Field timepoint samples (Fig. 4A). These two samples were consistently more abundant in multiple CAZyme families than other Field samples (Fig. 3C). The composition analysis established that the majority of the maturation in CAZyme potential takes place between Pre and Post stages in colony life, but likely continues until colonies reach maturity.
Enzymatic activity of termite guts shifts with fungus acquisition
Relative abundance of CAZymes can inform the potential metabolic profiles of an environment but is not always reflective of actual enzymatic activity. Therefore, we profiled the enzymatic activity of termite guts for nine key enzymes, which metabolise plant and fungal compounds, through AZCL enzymatic activity assays. Multivariate compositional analysis of Bray Curtis distances between enzyme activity profiles demonstrated that timepoint explained almost all of the observed variation (PERMANOVA: F2,33 = 383.7, R2 = 0.9588, p < 0.001). Similar to CAZyme family profiles, pairwise comparisons showed that timepoints were significantly different from each other, and that Pre again stood out as most different (Pre vs Post: F1,22 = 648.2, R2 = 0.9672, p < 0.001: Pre vs Field: F1,22 = 523.6, R2 = 0.9597, p < 0.001; Post vs Field: F1,22 = 47.90, R2 = 0.6853, p < 0.001). The extreme explanatory power of timepoint was also evident from PCoA clustering into three distinct groups with little dispersion (Fig. 4A). As for the predicted CAZyme relative abundances, this indicates that the majority of enzyme activity difference is between Pre and Post stages of colony life but that activity continues to change until colony maturity.
The extreme variance between Pre and both Post and Field appeared to be primarily driven by enzyme activities being absent or very low prior to fungus acquisition (Fig. 4B). Since genes encoding these CAZymes were consistently present (Fig S3), this is conceivably due to metabolic dormancy of the gut microbiota at this point in time. For enzymes that were active in Pre, Post and Field colonies, timepoint had a major effect (ANOVAs: Alpha-amylase: F3,33 = 52.77, p < 0.001, Cohen’s F = 1.79; Endo-cellulase: F2,33 = 11.83, p < 0.001, Cohen’s F = 0.85). Tukey HSD pairwise comparisons showed that Alpha-amylase enzymatic activity was higher in Pre than Post (p < 0.001) and Field (p < 0.001), while Post was higher than Field (p = 0.0030). Pre and Post were comparable in Endo-cellulase activity (p = 0.514) and both were higher than in Field colonies (Pre: p < 0.001 and Post: p = 0.0033) samples (Fig. 4B). This could potentially be explained by an increased need for bacteria-derived cellulases prior to the establishment of the fungal combs where plant catabolism is mostly taken over by Termitomyces. However, given that the AZCL assays were conducted on whole termite guts, including any fungus comb they had consumed, some enzyme activity could be attributed to the consumed fungus comb in addition to that from the gut bacteria.
The termite gut microbiome provides a diverse range of nitrogen cycling potential
Similarly, to explore the nitrogen cycling potential of the fungus farming termite gut microbiomes, we investigated genes found in the bacterial gut microbiome, termite host, and fungus. These were curated with NCycDB with DIAMOND as part of the NCyc tool50. We identified 55 unique genes from all seven families categorised in NCycDB. Complementarity across members of the symbiosis was high but the gut microbiome contributed a substantial number of otherwise unavailable nitrogen cycling genes (Fig. 5A). Termitomyces and the termite host shared 43.6% and 21.8% of genes across timepoints, respectively, and these were mainly involved in organic degradation, organic synthesis, and assimilatory nitrate reduction (Fig. 5A).
When comparing the logged relative abundances of NCyc profiles of the gut microbiome across timepoints, heatmaps and hierarchical clustering revealed four distinct groups (Fig. 5B). Group one was present in high abundance across all timepoints, group two was mostly present in Post and variably in Field but absent in Pre, while group three was always present in Pre but more variably in abundance in Post and Field. Within group three, there was a subgroup that for the most part was absent in Field. Finally, group four, similarly to group one, was consistently present albeit in lower relative abundance. Grouping genes by broad metabolic pathways showed that most pathways were well represented across timepoints with a substantial proportion of the nitrogen cycling gene abundance being involved in organic degradation and synthesis, and denitrification, with minor variation in the relative abundances across genes involved in these pathways. Notably, genes involved in nitrogen fixation were present in high abundance at all timepoints. The most variable pathway across timepoints was Anaerobic Ammonium Oxidation (anammox), which was essentially absent before fungus acquisition. Nitrification was essentially absent across all timepoints. The patterns identified here are similar to those observed for the CAZymes: the majority of nitrogen cycling-related potential was present across timepoints, with only a few losses or gains post-fungus acquisition. This potential covered most of the nitrogen cycling pathway and notably also nitrogen fixation, for which microbial genes were present at all timepoint and are lacking in the termites and in the fungus.
We then explored nitrogen cycling gene diversity with Shannon diversity and Observed richness. Similar to CAZyme potential, timepoint significantly affected both nitrogen cycling gene diversity and richness (ANOVA: Shannon: F2,16 = 13.20, Cohen’s F = 1.28, p < 0.001; Observed: F2,16 = 182.08, Cohen’s F = 4.77, p < 0.001) (Fig S1B). Shannon diversity was consistent across matrilines, but richness was not (Shannon: F2,16 = 0.055, p = 0.6469; Observed: F2,16 = 17.07, Cohen’s F = 2.11, p < 0.001). Pairwise comparisons of Shannon diversity with TukeyHSD demonstrated Field to be lower than both Pre and Post, whilst Post is the highest in Observed richness followed by Pre and then Field (Fig S1, see Table S5 for details).
Multivariate PERMANOVA analysis further demonstrated timepoint explained more than 40% of the variation in nitrogen cycling gene profile compositions. Timepoint, similar to that of CAZyme potential (PERMANOVA: F2,16 = 10.65, R2 = 0.4320, p < 0.001), with compositions being consistent across matrilines (F2,16 = 0.6165, p = 0.6982). Pairwise comparisons demonstrated that timepoints were significantly different from each other (Pre vs Post: F1,14 = 17.58, R2 = 0.5039, p = 0.001: Pre vs Field: F1,10 = 7.991, R2 = 0.3157, p = 0.004; Post vs Field: F1,8 = 5.825, R2 = 0.2059, p = 0.027). These groupings were also reflected in PCoA clustering (Fig. 5C), again with the two distinctly different samples (9b and 15b) having higher abundances of several genes (Fig. 5B). Maternal lineage remained insignificant across comparisons (p > 0.05) (Table S6).