Fungiculture in Termites Is Associated with a Mycolytic Gut Bacterial Community

Understanding functional capacities of gut microbiomes is important to improve our understanding of symbiotic associations. Here, we use peptide-based functional annotation to show that the gut microbiomes of fungus-farming termites code for a wealth of enzymes that likely target the fungal diet the termites eat. Comparisons to other termites showed that fungus-growing termite guts have relatively more fungal cell wall-degrading enzyme genes, whereas wood-feeding termite gut communities have relatively more plant cell wall-degrading enzyme genes. Across termites with different diets, the dominant biomass-degrading enzymes are predominantly coded for by the most abundant bacterial taxa, suggesting tight links between diet and gut community compositions.

T ermites are widespread in tropical, subtropical, and warm temperate regions (1) and form a diverse group of more than 3,000 described species in 281 genera and seven families (2)(3)(4)(5). They have major impacts on their environments (1), and this success has been attributed to their capacity to use nutritionally imbalanced, recalcitrant food sources, allowing for colonization of otherwise inaccessible niches (6). Different termites forage on distinct substrates, including soil, wood, dung, and fungus (7,8), decomposed through intricate interactions with complex gut microbial communities (6,9). In most termites, the main role of gut microbiota is believed to be the digestion of lignocellulose (10,11), but gut microbes also play key roles in nitrogen fixation (12)(13)(14), microbial defense (see, for example, reference 15), and immune regulation (16,17), which have major importance for the evolutionary history of the symbioses.
Approximately 30 million years ago, the basal higher termite subfamily Macroter-mitinae engaged in a mutualistic association with Termitomyces fungi (18,19) and have a distinct composition of the gut microbiota (20)(21)(22). Termitomyces decomposes plant material within external fungus gardens (combs) (23,24), but the gut still remains central in the association because plant substrate is macerated and mixed with asexual Termitomyces spores in a first gut passage prior to comb deposition (25). After Termitomyces breaks down the plant substrate, the termites ingest mature parts of the comb in a second gut passage (25), where gut microbes may contribute enzymes for final digestion of any remaining plant components (24). This division of labor is consistent with gut bacteria being of importance mainly when the comb material passes through the termite gut in a second passage (cf. references 23, 24, and 26), but recent work has suggested that partial lignin breakdown may also be accomplished during this first gut passage in Odontotermes formosanus (27). A set of microbes distinct from the gut microbiota of other termites persists in the fungus-growing termite guts, but limited work has examined functional implications of these differences (8,24). It has been hypothesized that it was associated with the more protein-rich fungal diet (20) and/or to break down chitin and other fungal cell wall components (8,24,28). Termitomyces domestication exposed fungus-growing termite gut communities to large quantities of fungal cell wall glucans (composed of D-glucose monomers), chitin (glucosamine polymer), and glycoproteins (see, for example, reference 29). Their breakdown requires a combination of carbohydrate-active enzymes (CAZymes; www.cazy.org) (30,31) and fungus-growing termite gut bacteria indeed encode glycoside hydrolase (GH) families of enzymes that may cleave chitin (GH18, GH19, and GH20), ␤-glucan (GH55, GH81, and GH128), and ␣-mannan (GH38, GH76, GH92, GH99, and GH125) (8,24).
In nature, bacteria are the major chitin degraders and its hydrolysis has been correlated with bacterial abundances in, e.g., soil communities (32). In fungus-growing termites, Bacteroidetes and Firmicutes bacteria appear to be the main producers of CAZymes putatively producing mycolytic enzymes, i.e., enzymes that lyse the fungal cell wall (8,24). These studies remained preliminary, however, because they were based on either an unassembled low-coverage metagenome (8) or had limited functional predictions (24). Here, we sequenced the gut metagenome of the fungus-growing termite Odontotermes sp. and performed in silico analyses to elucidate its fungal and plant cell wall-degrading capacities at deeper functional levels (i.e., to EC numbers when possible), assigned putative enzyme functions to gut community members using peptide-based functional annotation where prediction of function was confirmed by more than one method. To investigate the link between termite diet and gut community composition, we compared our findings to metagenomes from the fungusgrowing termite Macrotermes natalensis (24) and seven non-fungus-growing termite species feeding on plant material at different degrees of decomposition: the dung feeder Amitermes wheeleri (33), the two wood feeders Nasutitermes corniger and Microcerotermes parvus, a litter feeder Cornitermes sp., the two humus feeders Termes hospes and Neocapritermes taracua, and the soil feeder Cubitermes ugandensis (34). We reveal that the difference in gut community composition is associated with the presence of a mycolytic microbiota, providing insights into digestion and the role of gut communities in the fungus-growing termite symbiosis.

RESULTS
Taxonomic composition of fungus-growing termite gut microbiotas. We assigned bacterial taxonomies to metagenome contigs by searching for the closest matches of protein-coding genes on each contig against the NR database in NCBI and compared the relative abundance of the contigs in each group to assess the composition of termite gut microbiotas. M. natalensis and Odontotermes sp. were distinct from other higher termites primarily being relatively richer in Bacteroidetes (Fig. 1A), corroborating previous work (20,35), and termites in the same feeding group tend to be similar in gut microbiota composition ( Fig. 1B; see Table S1 in the supplemental material). Taxonomic compositions at the phylum level were consistent with previous 16S rRNA amplicon surveys of fungus-growing termites (8,20,21). Firmicutes and Bacteroidetes dominated in Odontotermes sp. (24 and 25%, respectively) and M. natalensis (24 and 30%, respectively) ( Fig. 1A; see also Table S1 in the supplemental material), comparable to an average 35 and 32% abundance in 16S-rRNA studies of Macrotermes subhyalinus and Odontotermes sp. in the Ivory Coast (21) and Odontotermes yunnanensis from Southwest China (8). Spirochaetes were relatively abundant in wood-feeding termites (46 to 50%) and Cornitermes sp. (32%), and they were low in relative abundance in fungus-growing termites (3 to 6%). At the genus and family levels, we identified the presence of genera that are part of the core gut microbial community in fungus-growing termites (e.g., Alistipes, Treponema, Dysgonomonas, Desulfovibrio, Ruminococcaceae, and Lachnospiraceae) (21). Alistipes and Bacteroides were most abundant in fungus-growing termites, representing 5 to 8% and 2 to 4% total abundances, respectively, sharply contrasting with other higher termites (on average, 0.1 and 0.5%, respectively) (Table S1). Treponema (Spirochaetes) was low in relative abundance in fungus growers and termites feeding on decaying plant material (3% and 6%, respectively), except for the litter feeder Cornitermes sp. (32%), while it was the most abundant taxon in wood feeders (42 to 45%) (Table S1). Comparisons to composition estimates from 16S rRNA studies (21,22,34) revealed that some genera were underrepresented in the metagenomes. Alistipes, for example, were found in much higher relative abundances (12% in average) from classifications using 16S rRNA than protein-coding genes (5 to 8%) (Table S3). Some taxonomic groups classified in 16S rRNA surveys, such as the TG3 phylum (20,34) were not detected in the metagenomes (Table S1).
Fungus and plant cell wall-targeting enzymes. To gain insights into the functional capacity for carbohydrates degradation, we first identified carbohydrate-active enzyme (CAZyme) families (30,31) and classified the genes by their substrate target and thus putative enzyme function by assigning EC numbers using peptide-based functional annotation (36)(37)(38) (Table 1; see also Table S2 in the supplemental material). We focus our presentation and comparisons to enzymes, for which the prediction of function was confirmed by more than one method (for details, see Materials and Methods). Principle component analysis (PCA) of glycoside hydrolase (GH) family compositions (Fig. 2B) support that gut microbial enzyme capacities are similar for termites with similar diets. Fungus-growing termite guts were thus systematically different in GH family gene composition to the guts of termites feeding on plant biomass ( Fig. 2B; see also Table S2 in the supplemental material), reflecting the differences in cell wall composition of plant and fungus material. Genes encoding GH enzymes that putatively target fungal cell wall components (detailed based on EC numbers in the paragraph below) were more abundant in fungus-growing termites, with the most marked differences being for GH125 (34-fold) and GH92 (14-fold) that likely target ␣-mannose (Table S2). Genes encoding GH families GH17, GH128, GH55, and GH87 that contain (among others) ␤-1,3-glucan-targeting enzymes and GH18 and GH19 that include several chitinases showed 2-to 12-fold higher relative abundance in fungus-growing termites (Table S2). In contrast and as expected, many GH family encoding genes primarily targeting plant cell wall components were low in relative abundance in fungus-growing termites but higher in wood-feeding termites (Table S2). Examples include GH94, GH10, and GH5 containing, e.g., cellulases, and GH74, GH120, GH39, GH26, GH10, and GH11 that contain hemi-cellulases, the genes of which were 3to 5-fold more abundant in wood feeders than in other termites (Table S2).
The bacteria encoding the most abundant cell wall-targeting enzymes. To gain further insight into the functional contribution of gut microbiota members to carbohydrate degradation, we grouped enzymes targeting polysaccharides from fungus and plant cell wall components by their taxonomy (Fig. 3). Clostridiales and Bacteroidales contributed most of the cell wall-degrading enzyme genes in fungus-growing termites (83% in M. natalensis and 68% in Odontotermes sp.). The two orders may, however, differ in what enzymes they contribute. The majority of the endo-1,2-␣-mannanases, endo-1,6-␣-mannanases, and exo-1,6-␣-mannosidases (60%) that target ␣-mannan of the fungus cell wall were coded for by members of Bacteroidales, whereas the exo-␣mannosidases were mainly contributed by members of Clostridiales (Fig. 3). Members from both orders contribute genes for endo-1,3-␤-glucanase and exo-1,3-␤-glucanase that target ␤-glucan in the fungus cell wall, but only Bacteroidales contribute endo-1,6-␤-glucosidase.
In wood feeders, genes encoding enzymes for cellulose and hemicellulose cleavage were coded for by all abundant members of the gut microbiotas (Fig. 3). Spirochaetales contributed 32% of the cell wall-degrading enzyme genes (Fig. 3), but notably, a large proportion of these genes (36%) could not be assigned to bacterial orders. Genes putatively encoding chitinases and ␤-N-acetylhexosaminidases for chitin degradation were contributed by both Clostridiales and Bacteroidales and, to a lesser extent, Spirochaetales in wood feeders. In termite species that were neither fungus nor wood feeding, most of the fungus cell wall-degrading enzyme genes were also contributed by Clostridiales, Bacteroidales, and Spirochaetales (17 to 92%), but the relative abundances of these enzyme genes were far lower than those observed in fungus and wood feeders (Table 1 and Fig. 3).

DISCUSSION
Diet is a major driver of taxonomic and functional composition of gut microbial communities (39)(40)(41). Our characterization of cell wall degrading enzyme from higher termites confirms that the functional profiles of gut microbial communities are tightly linked with termite diet, which is consistent with previous work that has focused on community compositions (21,22,35), and preliminary functional associations with microbiota structure in the Macrotermitinae (8,24). Fungus-growing termite diets consist of fungal hyphae and plant material that is partly degraded by the symbiotic fungus (22). The gut microbiota thus not only functionally complements final plant biomass decomposition by providing oligosaccharide-targeting enzymes (8,24) but also provides key enzymes for the digestion of fungal biomass. The mycolytic potential of the gut microbiota of the South African Odontotermes sp. and M. natalensis (this study) was comparable to O. yunnanensis from Southwest China (8), suggesting conserved functions across space and time and highlighting the robust link between taxonomy and function of the intimate interactions between gut community members and termite hosts.
Fungus-growing termite gut mycolytic enzymes were primarily coded for by Clostridiales and Bacteroidales, which also dominate the core microbiota of the termite subfamily (20,21,35,43,44). The ancestor of the Macrotermitinae likely had a bacterial gut microbiota similar to those of lower termites (but without protists) (6,9). Fungiculture exposed the gut microbiota to larger amounts of fungal biomass than the ancestral lignocellulolytic diet, likely resulting in gut microbiotas that became more similar to those observed in extant cockroaches (20,21). This is likely a product of several factors. First, bacterial strains present in the Macrotermitinae ancestor that could utilize a fungal diet were likely selected for and consequently increased in relative abundance (e.g., Alistipes, Dysgonomonas, and members of the Ruminococcaceae) (20,21). Many mycolytic microbes were conceivably already present in termite guts prior to the origin of fungiculture, since the ancestral termites fed on partially degraded plant substrates containing fungal biomass. Second, bacteria that do not contribute to the breakdown of fungal material in the Macrotermitinae ancestor were likely outcompeted/selected against in the fungal-material-rich environment. Third, novel lineages adopted from other termites or the environment were likely coopted when fungiculture evolved. This is consistent with recent work demonstrating rampant horizontal transmission of gut bacterial lineages associated with termites (45). Collectively, this led to a gut microbiota with relatively more mycolytic enzymes (i.e., ␣-mannanases, ␤-1,3/1,6-glucanases, and chitinases) and relatively fewer lignocellulolytic enzymes (i.e., cellulase and hemicellulose).
In sharp contrast to the fungus growers, lignocellulose-degrading enzymes dominated the gut microbiota of wood-feeding termites (N. corniger and M. parvus), as expected from the requirement for the breakdown of recalcitrant plant components, which in fungus farmers is handled by Termitomyces (24; although lignin cleavage may be initiated during the first gut passage [27]). Interestingly, however, the relatively high abundance of chitinolytic enzymes of Spirochaetales origin in wood-feeding termites suggests that the decaying wood these species feed on harbors fungal biomass that gut bacteria, or the termite host, likely utilize. However, fungal cell wall-targeting enzymes such as ␤-1,3-glucanase may also serve to protect against fungal infections (cf. reference 46). Experimental work targeting the expression of these enzymes in the presence of fungal pathogens, ideally combined with a varying fungal biomass diet content, have the potential to shed light on their relative defensive and dietary roles.
An important caveat of our study is that we were limited by comparisons of bacterial abundances and encoded enzymes between whole-gut DNA extractions and metagenomes in the two fungus-growing termites with P3 compartment of wood-and litter-feeding termites and the lumen metagenome of the dung feeder. This may have biased our comparative analyses somewhat, but we believe that at least the P3 comparison to whole guts is reliable for the following reasons. First, the P3 compartment is expected to contain the vast majority of microbial cells of termite guts (e.g., Ͼ97% in the wood feeder N. corniger [47]). Consequently, even if relative abundances may differ across gut sections (cf. reference 34), the overall community structure and functions should be primarily driven by P3-residing bacteria. Nevertheless, the contribution of enzyme genes from some member such as Clostridiales, which is consistently more abundant in the P1 compartment (48), may have been underestimated when only the P3 region was sequenced. Second, comparisons of community structure between whole gut 16S rRNA from fungus-growing termite (21) with P3 community analyses (20) found nearly identical community structure (see Fig. 4 in reference 26). More importantly, although community composition of the gut fluid of N. corniger (33) largely resemble hindgut and P3 compartment analyses (20,34,47,48), we advocate that the comparisons to the lumen metagenome of A. wheeleri is taken with a grain of salt for two main reasons. First, lumen and gut wall-associated microbiomes may differ substantially (49,50). Second, although we normalize our comparisons by metagenome size, the low coverage of A. wheeleri most likely undersamples the true composition of the metagenome, thus probably precluding the identification of low-abundance bacteria and enzyme genes (Table 1). Third, although the impact is expected to be minor (cf. reference 42), discrepancies could be impacted by sequencing protocols (454 versus Illumina) due to potential sequencing biases associated with GC content, sequencing errors, and differences in read length.
The shotgun metagenomics data showed distinct compositional profiles of cell wall-degrading enzyme genes in termite guts consistent with expected gut community functions. However, without experimental data of the bacteria densities and complete bacteria genomes, our comparisons cannot reflect enzyme quantities and in situ activities (cf. reference 51). A large fraction of the contigs also do not have a taxonomic assignment (Table S1), and although identification of all abundant genera was consistent with previous 16S rRNA studies (21,22,34), some bacteria were absent or underrepresented in the metagenomes compared to 16S rRNA surveys and vice versa (Table S3). The lack of identification of the presence of members of the TG3 phylum in the metagenome is likely an artifact, as the phylum lies within the Fibrobacteres in the NCBI sequence database we used for taxonomic identification. The genera Chitinispirillum and Ruminiclostridium we identified in the metagenomes were absent in 16S rRNA analyses (21,22,34) because there were no rRNA gene references for these two genera in the database at the time they were analyzed and published (52). Also, the lack of appropriate reference genomes for termite-associated strains limits the classification of the metagenome to genus level, which is likely the cause for the underrepresentation of the genus Alistipes. Some Alistipes sequences may not have been identified because they do not match to the reference genome which are mostly isolated from the human microbiome (53)(54)(55)(56). Thus, limited by reference-dependent methods that require closely related sequences in the public databases, the contribution of cell wall-degrading enzymes of some taxa is likely to be underestimated. Deeper metagenome sequencing that can enable binning of sequences to improve taxonomic classification and functional predictions, coupled with functional studies of the gutcompartment-specific expression of bacterial enzymes, are thus warranted.
Similarities in host diet have been shown to drive convergence in the functional potential of gut microbes in other organisms (57,58). Selection for particular physiological traits may, however, not necessarily be directly linked to specific phylogenetic groups of microbes. Exploring communities associated with diverse fungus-growing hosts (and their associated fungi) would allow us to explicitly test for convergent evolution of mycolytic microbial communities, even if these were likely comprised of different microbial consortia. A number of other insects utilize fungus material as a nutrient source, including fungus-growing ants (59,60), some Drosophila species (61), the Malaysian mushroom-harvesting ant Euprenolepis (62), Sirex wood wasps (63,64), and ambrosia beetles (65)(66)(67). Different bacteria are likely to be involved, but predictions would be that predominantly fungal diets should select for microbial communities with comparable mycolytic capacities. Recent work on mycophagous Drosophila supports that this may be so. Bost et al. (68) demonstrated that gut bacteria in mycophagous Drosophila are implicated in fungal cell wall metabolism, with cysteine and methionine metabolism enzymes originating from Bacteroidetes, Firmicutes, and Proteobacteria gut microbiota members, phyla that are also abundant and mycolytic in fungus-growing termites.
The prevalence of microbial communities with ample mycolytic capacities in the guts of fungus-growing termite species suggests that the adoption of a fungal diet has been associated with a functional and compositional change in gut microbial communities at the onset of fungiculture in termites. The high relative abundance of fungal cell wall-degrading enzymes indicates adaptations to the decomposition of fungal diet, which is consistent with this capacity being absent or less in termites with predominantly plant-based diets. An exception is wood-feeding termites, for which wooddegrading fungi may also comprise an appreciable component of the termite diet. To improve our understanding of the link between the digestive function and the gut microbes in fungus-growing termites, further work will be needed to elucidate whether the functional capacities of the gut microbiota reflect the amount of fungal biomass in the diet and differences in properties of the fungal species fed on. Estimates of bacterial densities and enzyme activities are also needed to fully understand the capacity of biomass degradation in termite gut. Metagenomic analysis, including more fungusgrowing termite species, improved coverage, and longer reads, paired with referencefree taxonomy binning methods and functional annotation, will also largely enhance our understanding of the property of gut microbes in a functional angle.

MATERIALS AND METHODS
Odontotermes sp. collection. Termites from an Odontotermes sp. colony (code Od127) were collected at the Experimental Farm of the University of Pretoria, Pretoria, South Africa (Ϫ25.742700, 28.256517). The species identity of this colony had been previously established as by mitochondrial gene COII barcoding (22). Fifty old major workers were sampled, and entire guts were dissected and pooled in a 1.5-ml Eppendorf tube and stored at -80°C until DNA extraction.
Gut microbiota DNA extraction. Guts were ground in liquid nitrogen, and DNA was extracted using DNeasy blood and tissue kits (Qiagen, Hilden, Germany) according to the manufacturer's description, except that a chloroform extraction step was included after incubation with protease K. After proteinase K digestion, 1 volume chloroform-isoamyl alcohol (24/1) was added. The tubes were incubated for 15 min on a slowly rotating wheel and then centrifuged at 3,000 ϫ g for 10 min. The supernatant was transferred to spin columns, and the remainder of the manufacturer's protocol was followed. The quality and purity of samples were determined using NanoDrop (Thermo Scientific, Wilmington, DE).
Metagenome sequencing and assembly. DNA was sheared to ϳ350-bp fragments, end repaired, A tailed, and ligated with Illumina paired-end adaptors (Illumina). The ligated fragments were selected from the desired size on agarose gels and amplified by ligation-mediated PCR, and libraries were sequenced with 150-bp read lengths on an Illumina HiSeq2500. The quality of raw sequencing reads was assessed before assembly. Reads containing the adaptor, Ͼ10% N, or Ͼ50% low-quality bases (Q-score Ն 5) were removed. To exclude sequences from Termitomyces and the termite hosts, qualitycontrolled reads were mapped to the M. natalensis and Termitomyces genomes (24) using the Burrows-Wheeler Aligner v0.7.15 (69) BWA-MEM algorithm; any aligned reads were filtered.
Clean reads were assembled by IDBA-UD v1.1.2 (70,71) with an iterative set up from k-mer size of 19 to 99 at step of 10 (-pre_correction -mink 19 -maxk 99 -step 10). Unassembled reads were picked out by mapping reads back to the initial assembly and assembled separately with the same setup. Redundancies of sequences from the same organism within the metagenome were removed by clustering all contigs at 95% identity with CD-hit v4.6.6 (72), and only the longest contig per cluster was kept. For comparison, high-quality reads of M. natalensis old major worker gut microbiota from Poulsen et al. (24) were also reassembled using the same procedure. Genes in both assemblies were predicted by Prodigal v2.6.3 (73) with metagenomics parameters (-c -p meta).
Relative abundances of contigs within metagenomes. Clean reads of fungus-growing termites gut microbiota were mapped to the contigs by Bowtie v2.2.9 (74) with subsequent masked duplication, taking the best match for each read. Contig coverage was first estimated by the contig length and the number of mapped reads per contigs. Mapped read numbers were scaled to an equivalent of 10 Gb of sequence per metagenome, after which the relative abundance of each contig was calculated as the coverage of the contig divided by the sum of coverages of all contigs (cf. reference 75). For non-fungusgrowing termites, coverage information of contigs was obtained from JGI (76,77), and the relative abundances were estimated as described above.
Metagenome taxonomic assignment. Taxonomic assignment of protein-coding genes was carried out using Diamond (18) alignment against the NR database in NCBI. Alignments with E values Ͼ1e-5 and sequences with identities Ͻ30% were removed. Taxonomic information of the top hit was assigned to each gene. Contigs referring to taxonomical levels were determined by a modified lowest common ancestor (LCA)-based algorithm implemented in MEGAN (78). Taxonomic classification supported by Ͻ10% of the genes on each contig were first filtered, and the LCAs for the taxonomic classification of the rest of the genes were assigned to the contig. The relative abundance of contigs belonging to the same taxonomic group was summed up to represent the taxonomic composition of that taxonomic group in the microbiota. Taxonomic composition at the phylum level for metagenomes was centered log ratio (CLR) transformed and then compared and visualized by principle component analysis (PCA) (79) using R v3.3.2 (80).
Carbohydrate-active enzyme analysis. We classified genes in CAZyme families by searching against CAZyme hidden Markov profiles in the dbCAN database (v7; CAZyDB accessed 31 July 2018 [81]) using HMMer v3.1b (E Յ 1e-5) (82). The longest matched profile was selected if CAZyme domain overlapped with peptide pattern recognition using HotPep (36)(37)(38). CAZyme family classifications that were supported by both approaches were kept. The composition of GH families in termite species was CLR transformed, and the differences between termite species were visualized using PCA (79) in R v3.3.2 (80). Enzyme functions of genes in each CAZyme family were determined by BLASTp searches (E Յ 1e-5) against the ExPASy enzyme records, CCD searches against the COG database (83), and GhostKOALA searches using the KEGG database online tool (84). Genes coding for enzymes related to fungal and plant cell wall degradation were selected, and their substrate targets and bacterial taxonomy were assigned and manually checked. To compare the fungal and plant cell wall-targeting enzymes across termite species, the relative abundance profiles were hierarchically clustered and are shown as a dendrogram in Fig. 2C. For each enzyme, the relative abundances were compared to the average across the nine termite species, and the fold change was log 2 transformed and is presented in a heatmap in Fig. 2C. Data availability. Clean reads and metagenome assembly have been submitted to the SRA and GenBank under BioProject accession numbers PRJNA476694 and PRJNA193472.