ARTICLE https://doi.org/10.1038/s42003-020-1004-3 OPEN Integrative omics analysis of the termite gut system adaptation to Miscanthus diet identifies lignocellulose degradation enzymes 1234567890():,; Magdalena Calusinska 1 ✉, Martyna Marynowska1,2, Marie Bertucci1,3, Boris Untereiner1, Dominika Klimek1, Xavier Goux1, David Sillam-Dussès4, Piotr Gawron5, Rashi Halder 5, Paul Wilmes5, Pau Ferrer 1,6, Patrick Gerin 3, Yves Roisin2 & Philippe Delfosse1,7 Miscanthus sp. biomass could satisfy future biorefinery value chains. However, its use is largely untapped due to high recalcitrance. The termite and its gut microbiome are con- sidered the most efficient lignocellulose degrading system in nature. Here, we investigate at holobiont level the dynamic adaptation of Cortaritermes sp. to imposed Miscanthus diet, with a long-term objective of overcoming lignocellulose recalcitrance. We use an integrative omics approach combined with enzymatic characterisation of carbohydrate active enzymes from termite gut Fibrobacteres and Spirochaetae. Modified gene expression profiles of gut bacteria suggest a shift towards utilisation of cellulose and arabinoxylan, two main components of Miscanthus lignocellulose. Low identity of reconstructed microbial genomes to closely related species supports the hypothesis of a strong phylogenetic relationship between host and its gut microbiome. This study provides a framework for better understanding the complex lignocellulose degradation by the higher termite gut system and paves a road towards its future bioprospecting. 1 BioSystems and Bioprocessing Engineering, Luxembourg Institute of Science and Technology, Rue du Brill 41, L-4422 Belvaux, Luxembourg. 2 Evolutionary Biology and Ecology, Université Libre de Bruxelles, Avenue F.D. Roosevelt 50, 1050 Brussels, Belgium. 3 Laboratory of Bioengineering, Earth and Life Institute, Université Catholique de Louvain, Croix du Sud 2/L7.05.19, 1348 Louvain-la-Neuve, Belgium. 4 Laboratory of Experimental and Comparative Ethology EA444, University Paris 13-Sorbonne Paris Cité, 399 Avenue Jean-Baptiste Clément, 93430 Villetaneuse, France. 5 Luxembourg Centre for Systems Biomedicine, University of Luxembourg, 7 Avenue des Hauts-Fourneaux, L-4362 Esch-sur-Alzette, Luxembourg. 6 Chemical, Biological and Environmental Engineering, Universitat Autònoma de Barcelona, Edifici Q, 08193 Bellaterra (Cerdanyola del Vallès), Spain. 7 University of Luxembourg, Rectorate, 2 Avenue de l’Université, L-4365 Esch-sur-Alzette, Luxembourg. ✉email: magdalena.calusinska@list.lu COMMUNICATIONS BIOLOGY | (2020)3:275 | https://doi.org/10.1038/s42003-020-1004-3 | www.nature.com/commsbio 1 ARTICLE COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-020-1004-3 I n a world of finite biological resources, the agenda of the UN’s LM1 and LM3) of Cortaritermes intermedius were fed exclusively Sustainable Development Goals challenges scientific commu- with dried Miscanthus straw (Supplementary Figs. 1 and 2). This nity to develop transformative technologies that would enable Nasutitermitinae genus is known to feed on grass tussocks in its the replacement of petroleum-based raw materials and energy natural habitat10. Alteration of the termite gut microbiome (here with renewable bio-based feedstock. Plant biomass (lig- relative to bacterial communities in termite midgut and hindgut) nocellulose), being the most abundant and renewable natural was monitored at monthly basis during 9 months, by high- resource, could have many applications in different sectors1. throughput sequencing the V6–V8 regions of 16S rRNA gene Miscanthus sp. is a rhizomatous grass and owing to its adapt- (Fig. 1a). Quality-trimmed reads were assembled into 678 ability to various environmental conditions, it shows high operational taxonomic units (OTUs) assigned to 18 bacterial potential for sustainable production of lignocellulose over large phyla. Spirochaetae and Fibrobacteres were the most dominant, geographical range2. Considering its important agronomic as previously shown for plant fibre-feeding higher termites (e.g. advantages (e.g. high biomass yield per hectare, reduced soil ref. 11; Supplementary Data 1). By assessing bacterial community erosion and low fertiliser and pesticide requirements), it is sui- structures in control samples (colonies feeding on grass tussocks table for different biorefinery value chains, including bioethanol, in situ) and Miscanthus-fed microbiomes, we could observe biogas, food additives, ingredients for cosmetics, biopharmaceu- radical changes. Species richness and diversity were significantly ticals, bioplastics, biomaterials, organic fertilisers and animal higher (HOMOVA p < 0.001) before Miscanthus diet was initi- feed3. Yet, due to the high recalcitrance (resistance of the cell wall ated, possibly reflecting an adaptive selection for the most effi- components to enzymatic hydrolysis), its use is largely untapped4. cient microbial degraders facing lower complexity of carbon In living organisms, enzymatic hydrolysis of lignocellulose is sources in comparison to original diet (Fig. 1b). Further appli- mainly driven by carbohydrate active enzymes (CAZymes5). Gly- cation of linear discriminant analysis (LDA) effect size (LEfSe; coside hydrolases (GHs) are the primary enzymes that cleave gly- ref. 12) to two termite colonies demonstrated that nearly 140 cosidic linkages. Often, they are assisted by carbohydrate esterases bacterial OTUs were significantly enriched in control micro- (CEs), polysaccharide lyases (PLs) and other auxiliary enzymes biome, while roughly 13 were enriched in Miscanthus-fed (AAs). With its unique consortium of microorganisms, the termite microbiome (Supplementary Data 2). Out of the latter, six OTUs gut is considered as the most efficient lignocellulose degrading assigned to Fibrobacteres (mainly representing a genus exclu- system in nature6. Complete loss of gut cellulolytic flagellates in all sively containing termite Fibrobaceteres sequences) and Spir- evolutionary higher termites and acquisition of novel symbiotic ochaetae (associated with the termite Treponema cluster) were bacteria led to improved lignocellulolytic strategies. It allowed for particularly abundant, and on average they accounted for 55.39% diet diversification from mainly wood-restricted to e.g. dry grass ± 3.8 of the Miscanthus-fed microbiome, in comparison to the and other plant litter, herbivore dung and soil organic matter at 29.87% ± 1.8 average abundance in control microbiome. By different stages of humification7. Until now, most research has analysing bacteria naturally associated with Miscanthus diet, we focused on endogenous endoglucanases of termites and cellulases estimated the effect of immigration on the termite gut community from termite gut flagellates8. CAZymes from higher termite gut as negligible (Supplementary Fig. 3). All together, these results bacteria have only recently started receiving scientific attention9. demonstrated that diet change drives the development of Here, we investigated the higher termite gut system of Cortar- microbial consortium in a unique manner, yet food-associated itermes sp. (Nasutitermitinae subfamily) from French Guiana microbes cannot compete with highly specialised termite gut savannah. Using 16S rRNA gene amplicon profiling of termite gut microbiota for a niche. bacteria, we investigated the adaptation of two termite colonies to Miscanthus diet under laboratory conditions. Through the de novo metagenomic (MG) and metatranscriptomic (MT) recon- Comparison of de novo metatranscriptomics and metage- structions, we assessed the distribution of activities within gut nomics. The de novo MT was applied to nest LM1 and co- community, and we linked it to different bacteria, two main assembling reads from three samples (LM1_1 “control sample”, players being Spirochaetae and Fibrobacteres. Further analysis of LM1_2 and LM1_8, both representing “Miscanthus-fed micro- gene expression profiles proved microbial functional plasticity biomes”) yielded 603,579 open reading frames (ORFs), mainly (adaptation to changing environmental conditions through dif- representing partially reconstructed gene transcripts. The de novo ferential genes expression), and highlighted the abundance of gene MG reconstruction was also applied to colony LM1 at sampling transcripts involved in carbohydrate metabolism and transport. point LM1_8 and it yielded 211,724 ORFs annotated on 64,347 Analysis of the reconstructed community metagenome evidenced contigs for the total assembly size of 177.5 million base pairs CAZyme clusters targeting two main components of Miscanthus (Mbp). For both datasets, initial public database-dependent biomass, namely cellulose and (arabino)xylan. Most of these taxonomic classification of genes and gene transcripts pointed clusters were allocated to the reconstructed metagenome assem- to the abundance of Firmicutes (Supplementary Fig. 4), what bled genomes (MAGs) of Fibrobacteres and Spirochaetae origin. contrasted the results of community structure analysis. Sub- The de novo reconstruction of the host epithelial gut tran- sequent binning of MG contigs and phylum-level annotation of scriptome evidenced its contribution to lignocellulose degradation, the resulting bacterial bins allowed assigning correct taxonomic and its adaptation to Miscanthus diet. Based on the characterisa- origin, confirming the metagenomic abundance of Spirochaetae tion of purified bacterial CAZymes, we verified the in silico pre- and Fibrobacteres (Fig. 1c; Supplementary Fig. 5). Following dicted activities for many backbone-targeting (e.g. endocellulases mapping of RNA reads to the MG assembly (referred to as “RNA- and xylanases) and debranching enzymes (e.g. arabinofur- seq” analysis), we could confirm transcriptional dominance of anosidases). To finish, we discussed our findings in the context of these two bacterial phyla as well. Incomplete public databases and enzymes application in the developing biorefinery sector. extensive horizontal gene transfer were previously proposed as the origin of this misclassification9. Based on the classification of genes and transcripts to broad Results and discussion functional categories such as KEGG ontology profiles (KOs), Structural alteration of gut microbiome under Miscanthus diet. congruency between the de novo MG and MT reconstructions To examine enzymatic degradation of Miscanthus by the higher was high (Supplementary Fig. 6). However, out of the de novo termite gut system, two laboratory-maintained colonies (nests MT reconstructed gene transcripts of prokaryotic origin, only 2 COMMUNICATIONS BIOLOGY | (2020)3:275 | https://doi.org/10.1038/s42003-020-1004-3 | www.nature.com/commsbio COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-020-1004-3 ARTICLE a b Bacterial diversity Bacterial richness nest LM1 nest LM3 LM1_9 600 600 (sobs index) LM1_8 400 400 LM1_7 200 200 0 0 LM1_6 LM1_ 1 2 3 4 5 6 7 8 9 LM3_ 1 2 3 4 5 6 7 (invsimson index) LM1_5 30 30 nest LM1 LM1_4 20 20 Miscanthus feeding LM1_3 10 10 LM1_2 0 0 LM1_ 1 2 3 4 5 6 7 8 9 LM3_ 1 2 3 4 5 6 7 LM3_7 Time (months) Time (months) LM3_6 LM3_5 c nest LM1 d LM3_4 Relative MG abundance (%) 100 Relative MT abundance (%) nest LM3 LM3_3 80 LM3_2 60 LM3_1 control 40 LM1_1 20 0 20 40 60 80 100 Relative community abundance (%) 0 0.02 16S rRNA gene amplicon LM1_1 LM1_2 LM1_8 LM1_8 Spirochaetae Firmicutes Actinobacteria Fibrobacteres Proteobacteria Chlorobi Bacteroidetes Acidobacteria other Fig. 1 Structural composition of Cortaritermes sp. gut microbiome under Miscanthus sp. diet. a Clustering of samples based on the calculated Bray −Curtis index and phylum level taxonomic assignment of sequencing reads from the 16S rRNA gene amplicon study. b Bacterial richness and diversity indices before (highlighted in yellow on sub-figures a and b) and under Miscanthus. Boxes represent the interquartile range and error bars show the 95% confidence intervals (n = 3). c Relative metatranscriptomic (MT) and metagenomic (MG) reads abundance assigned at the phylum level. Taxonomic gene and gene transcript assignments were inferred from the metagenomic contigs binning and phylum-level bin classification. d Cortaritermes sp. colony (top) and termite workers under the protection of soldiers while feeding on Miscanthus fibres in laboratory conditions (bottom). 37.8% showed significant similarity to the de novo MG genes at herbivorous insects, this strategy was proposed as a mechanism the protein level (blastp e-value ≤ 10−5), sharing on average contributing to the success of termites14 and herbivorous ants16 76.04% of amino acid identity (Supplementary Fig. 4). Coherent in their marginal dietary niches. Importance of lignocellulose to our study, differences in functional gene profiles between MG degradation under Miscanthus diet was evidenced by increased and MT reconstructions have been previously underlined13. Even abundance of transcripts broadly assigned to cellulose and xylan in the context of the termite gut, some authors highlighted the processing KOs (Fig. 2; Supplementary Data 4). Multiple sugar value of the de novo MT assembly in retrieving highly expressed ABC transporters were up-regulated in the Spirochaetae meta- genes9,14. transcriptome, while they were nearly absent from the MG and MT reconstructions of Fibrobacteres origin. This observation Genomic potential and transcriptional adaptation of gut could suggest the governance of exogenous carbohydrates uptake microbes. Aggregation of 68.9 ± 1.8% of the de novo recon- and utilisation by Spirochaetae. structed gene transcripts into clusters of orthologous genes (COGs) pointed at functional microbiome stability at the differ- ent stages of feeding campaign (Supplementary Data 3). Con- Diversity and abundance of termite gut microbial CAZymes. sistently with previous reports11,14,15, cell motility and The de novo MT reconstruction yielded over 2000 of manually chemotaxis together with carbohydrate transport and metabolism curated transcripts assigned as CAZy-coding genes (cazymes). were the two most highly expressed gene categories. Recon- Out of these, 38.4% were further assigned to 55 GH families. The struction of (nearly) complete metabolic modules was quite de novo MG reconstruction resulted in close to 7000 different similar between Fibrobacteres and Spirochaetae. However, further cazymes, 43.6% of which were assigned to 86 GH families (Sup- comparative analysis using LEfSe (Fig. 2; Supplementary Data 4 plementary Figs. 6 and 7). Although there was a good correlation and 5) identified several biologically informative features differ- between the distributions of identified GHs to different GH entiating these two bacterial phyla. Both were capable of nitrogen families (Pearson r 0.83), roughly 150 cazymes identified in the fixation and glycogen synthesis, but the two pathways were de novo MT were also reconstructed from the assembled meta- enriched in Fibrobacteres. Expression of Amt ammonium trans- genome. Novelty of reconstructed cazymes was evidenced porters was highly up-regulated, and together with increased through sequence comparison to NCBI refseq database, and a abundance of gene transcripts involved in urea transport and metagenomic dataset previously generated for Nasutitermes sp.15 metabolism (restricted to Spirochaetae), it indicated nitrogen (Supplementary Fig. 4). In the latter case, average amino acids deficiency of a Miscanthus-fed termite colony. Both Spirochaetae identity for the 943 query hits equalled 65.4 ± 19.9% (blastp, and Fibrobacteres could also synthetise ten essential amino acids e-value ≤ 10−5), pointing to the diversity of microbial cazymes in that animals cannot synthetise de novo. Even though nitrogen guts of different termite species, even those phylogenetically provisioning by bacterial symbionts is not employed by all closely related. COMMUNICATIONS BIOLOGY | (2020)3:275 | https://doi.org/10.1038/s42003-020-1004-3 | www.nature.com/commsbio 3 ARTICLE COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-020-1004-3 a b Chemotaxis and motility Transporters Carbohydrate metabolism Nitrogen fixation and regulation Ribosome Genetic information processing other uncharacterised c carbohydrates uptake and metabolism - gut microbiome xylan cellulose, starch (hetero)xylans chitin β-mannan xylose xylose, uronates β-mannan ABC/PTS transp. ? ABC termite gut glucose glucose transp. lumen Fibrobacteres PTS Spirochaetae PTS glycolysis glycolysis transp. transp. glucose ribose synthesis glycogen pyruvate pyruvate reductive essential amino acids essential amino acids acetogenesis lactose, galactose chitin trehalose, α-mannan starch (amylose) cellulose, cellobiose CO CO2 acetate acetate (essential) (essential) acetate amino acids degradation amino acids glycogen glycogen glycogen host purine synthesis (storage) metabolism glutamine glutamine glycogen carbohydrates glutamate glutamate uric acid metabolism urea uptake and - termite biotin degradation NH3 synthesis biotin NH3 urea ABC transp. N2 fixation N2 fixation NH3 NH3 biotin N2 N2 Fig. 2 Functional characterisation of the termite gut system feeding on Miscanthus sp. a, b Tag clouds of enriched (LefSe LDA > 2, p < 0.05) KOs reconstructed from the de novo metatranscriptomics for the termite gut Fibrobacteres (a) and Spirochaetae (b) at LM1_8. Top 25 most abundant KOs are displayed. Size of the text reflects transcriptomics abundance of a specific KO. c Simplified metabolic reconstruction, with a focus on carbohydrate metabolism, for the termite gut lignocellulolytic system. Hypothetical pathways are indicated with dashed lines. Metabolic pathways enriched in Fibrobacteres and Spirochaetae (metatranscriptomes) are indicated with bold lines. Metabolites putatively shared between gut bacteria and the host are indicated with square boxes. Differential expression of specific cazymes at different stages of abundant but only partially reconstructed GH11 gene transcripts, the feeding experiment suggested quick acclimation to new what reduced initial over-dominance of this CAZy family by 3.3 (laboratory) conditions, also reflecting adaptation of gut microbes ± 0.9-fold (Supplementary Fig. 8). Similarly, only highly frag- to digest Miscanthus (Fig. 3). There were roughly 29.7% of mented GH11 genes were recovered from the reconstructed common cazymes transcripts between LM1_8 and control metagenome. We hypothesise that closely related Spirochaetae LM1_1 sample, while over 55% were shared between LM1_8 strains contain highly similar GH11 genes, possibly shared by and LM1_2 (both fed with Miscanthus). Along the experiment, horizontal gene transfer; likewise, it was shown for e.g. human gut GH5 (mainly subfamilies GH5_4 and GH5_2) was the most Bacteroidetes17. Their conserved nature may impede proper gene highly expressed family. Still, its cumulative expression nearly reconstruction from sequencing reads, similarly to other doubled under Miscanthus diet (Fig. 3b). Other abundant families structurally conserved genes such as 16S rRNA or transposase, included GH43, GH10 and GH11, all potentially involved in that typically do not reconstruct into larger genomic fragments18. (hetero)xylan degradation. The latter was previously shown as GH45 was only present in Fibrobacteres and it was the second largely expressed by the termite gut fibre-associated Spirochae- most expressed GH family under Miscanthus diet. According to tae9. Following manual curation, we removed three highly the CAZY database, nearly 95% of proteins assigned to this family 4 COMMUNICATIONS BIOLOGY | (2020)3:275 | https://doi.org/10.1038/s42003-020-1004-3 | www.nature.com/commsbio COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-020-1004-3 ARTICLE Fig. 3 Carbohydrate active enzymes (CAZymes) reconstructed from metatranscriptomic (MT) and metagenomic (MG) reads for the termite gut system. a Venn diagram showing the number of assigned glycoside hydrolase (GH) families for the termite gut microbiome and host gut epithelium. b, c Comparison of the gene expression profiles (de novo MT, log2 transformed) for gene transcripts assigned to the different GH families and at the different stages of the Miscanthus feeding experiment, for the gut microbiome (b) and the host gut epithelium (c). d Venn diagram showing the number of assigned GH families for Fibrobacteres and Spirochaetae, based on the de novo MG reconstruction. e Average CAZyme genes expression and cumulative gene expression at the different stages of the feeding experiment and analysed separately for Fibrobacteres and Spirochaetae. Lower panel is a zoom on the gene expression profiles with the outliers (highly expressed genes; in some cases representing only partially reconstructed genes) removed. Boxes represent the interquartile range and error bars show the 95% confidence intervals (n = number of transcripts annotated as glycoside hydrolases). f, g Number of genes (f) and cumulative abundance of the most abundant GH families (g RNA-seq log2 transformed) at the time point LM1_8 (the end of the Miscanthus feeding experiment), and visualised separately for Fibrobacteres and Spirochaetae. Transcripts abundance (g) is calculated based on the RNA mappings (RNA-seq) to the MG contigs. Shaded parts correspond to shared GH families. are of eukaryotic origin (Supplementary Fig. 10), and show peptide-based functional annotation21, we further assigned in endocellulase activity. Previously, GH45 CAZymes were char- silico specific functions (EC numbers) to 60.3 ± 2.2% of gene acterised from lower termite symbiotic protists19. Lytic poly- transcripts classified as GHs. In many cases, these predictions saccharide monooxygenases (LPMOs) typically assigned to AA9 were experimentally validated (Supplementary Data 6). We (fungal) and AA10 (predominantly bacterial enzymes) families20 confirmed β-glucosidase, endocellulase, endoxylanase and arabi- were neither present in the reconstructed metagenome nor nofuranosidase activities of several Spirochaetae CAZymes. We metatranscriptome. also characterised active endoxylanases and endomannanases from Fibrobacteres. Abundance of transcripts associated with endocellulase Expression and activities of GHs from Fibrobacteres and (EC:3.2.1.4) and endoxylanase (EC:3.2.1.8) increased under Spirochaetae. Based on the RNA-seq analysis, Fibrobacteres and Miscanthus diet (for both Fibrobacteres and Spirochaetae), while Spirochaetae expressed respectively 47.9 ± 14.8% and 45.6 ± those involved in chitin and starch (α-glucans) degradation 18.5% of their cazymes genomic content when the termite was fed decreased (Fig. 4; Supplementary Fig. 12). Endocellulase-assigned with Miscanthus. Total diversity of Spirochaetae cazymes was transcripts were nearly equally abundant between Fibrobacteres 2.3 ± 0.5-fold higher than for Fibrobacteres (Fig. 3d–g). Their and Spirochaetae, while abundance and diversity of endoxyla- cumulative transcriptional abundance was also higher for Spir- nases of Spirochaetae origin was much higher. Most of the ochaetae; however, calculated average gene expression was assigned endocellulases were classified as GH5_4 enzymes slightly higher for Fibrobacteres. This observation was consistent (Supplementary Fig. 8.). Phylogenetic reconstruction comprising across different GH families (Supplementary Fig. 11). As various the previously characterised CAZymes from this family revealed GH families are characterised with broad functionalities, using the presence of multiple protein clusters separately grouping COMMUNICATIONS BIOLOGY | (2020)3:275 | https://doi.org/10.1038/s42003-020-1004-3 | www.nature.com/commsbio 5 ARTICLE COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-020-1004-3 a Gut microbiome xyloglucan (XG) EC:3.1.1.72 AC AC degradation Termite host EC:3.2.1.21 EC:3.2.1.151 lignin 10-30% EC:3.2.1.177 EC:3.2.1.21 AA15 EC:2.4.1.49 lignin degradation LPMO + P EC:3.2.1.51 EC:3.2.1.4 PTS EC:3.2.1.22 cellobiose P cellulose EC:2.4.1.20 EC:3.2.1.86 EC:1.10.3.2 EC:3.2.1.21 40-60% X EC:3.2.1.91 + P arabinoxylan (AX) 20-40% + cellulose degradation degradation EC:3.2.1.37 hemicellulose AX EC:3.2.1.39 EC:3.2.1.8 XG βG β1-3 β1-4 EC:3.1.1.73 FA EC:3.2.1.55 EC:3.2.1.74 FA FA Glucose Galactose EC:3.2.1.4 EC:3.1.1.72 AC Xylose AC Fucose Arabinose EC:3.2.1.21 Lignin subunits FA Ferulic acid β-glucan (βG) AC Acetic acid degradation b c cellulose degradation other carbohydrate sources backbone-targeting enzymes debranching enzymesother carbohydrate sources EC:3.2.1.132 EC:1.10.30.2 EC:3.2.1.151 EC:3.2.1.177 AA15 LPMO EC:3.2.1.21 EC:3.2.1.28 EC:2.4.1.25 EC:2.4.1.25 EC:3.2.1.41 EC:3.2.1.14 EC:3.2.1.52 EC:3.2.1.14 EC:3.2.1.52 EC:2.4.1.49 EC:3.2.1.21 EC:2.4.1.20 EC:3.2.1.86 EC:3.2.1.37 EC:3.2.1.39 EC:3.2.1.74 EC:3.2.1.55 EC:3.2.1.22 EC:3.2.1.51 EC:3.1.1.72 EC:3.1.1.73 EC:3.2.1.4 EC:3.2.1.1 EC:3.2.1.1 EC:3.2.1.4 EC:3.2.1.8 trehalose chitin starch cellulose hemicellulose glycogen chitin starch Normalized cumulative d transcripts abundance Relative abundance (RNA-seq) at LM1_8 (log ) 12 Fibrobacteres 120 Spirochaetae Number of expressed genes (log2) 10 Fibrobacteres genes 100 Spirochaetae genes 8 80 6 60 4 40 2 20 0 0 3.2.1.4 2.4.1.49 3.2.1.21 2.4.1.20 3.2.1.86 3.2.1.8 3.2.1.37 3.2.1.151 3.2.1.39 3.2.1.74 3.2.1.55 3.2.1.22 3.2.1.177 3.1.1.72 3.1.1.73 3.2.1.41 3.2.1.14 3.2.1.52 3.2.1.132 3.2.1.1 2.4.1.25 cellulose hemicellulose debranching other carbon sources Fig. 4 Characterisation of the termite gut lignocellulose degradation strategies. a Simplified overview of enzymatic pathways involved in the degradation of main components of the Miscanthus biomass, based on enzymes (gene transcripts assigned an EC number) revealed in our study. Dashed lines indicate hypothetical pathways. Lignin subunits correspond to: p-hydroxyphenyl (H), guaiacyl (G), and syringyl (S). CAZymes gene expression profiles (de novo MT) at the different stages of the Miscanthus feeding experiment analysed for the termite gut epithelium (b) and the termite gut microbiome (c). Gene expression analyses were done separately for the termite gut epithelium and the gut microbiome; therefore, data on sub-figures b and c should only be compared within a single sub-figure. d Relative CAZymes gene transcripts abundance (RNA-seq) and gene numbers assigned to different enzymatic categories and analysed separately for Fibrobacteres and Spirochaetae for LM1_8 sample. Spirochaetae and Fibrobacteres GHs (Fig. 5a). Concurrent CAZyme representing Spirochaetae cluster IX. Purified protein inspection of reconstructed genomic fragments suggested the was shown to be an endocellulase acting on carboxymethylcellu- existence of different cazymes loci containing GH5_4 genes lose (CMC) and glucomannan (Fig. 5b). In addition, activity on (Supplementary Fig. 13). Interestingly, CAZymes previously xylan and arabinoxylan was confirmed. This gene was also one of characterised to possess single enzymatic activity (mostly the most highly expressed cazymes under Miscanthus diet, endocellulase and to a lower extent endoxylanase) grouped in hypothesising the importance and interest for bacteria to express upper part of the tree. Lower part of the tree mainly contained multi-functional enzymes. To our best knowledge, it represents multi-functional enzymes (single enzyme simultaneously the first GH5_4 CAZyme of Spirochaetae origin ever charac- acting on cellulose and xylan). Suggested enzymatic multi- terised, and first multi-functional enzyme of higher termite gut functionality was further confirmed for a selected GH5_4 prokaryotic origin. 6 COMMUNICATIONS BIOLOGY | (2020)3:275 | https://doi.org/10.1038/s42003-020-1004-3 | www.nature.com/commsbio COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-020-1004-3 ARTICLE a MA Spirochaetae MA cluster I MA Spirochaetae cluster II multi-functional enzymes Fibrobacteres cluster III mono-functional enzymes Spirochaetae Spirochaetae cluster IV cluster IX Fibrobacteres cluster V Spirochaetae cluster IV Spirochaetae cluster VII b c 100 Relative enzyme activity (%) Spirochaetae cluster VIII 50 0 20 30 40 50 60 Temperature (°C) Fig. 5 Characterisation of the GH5_4 family. a Unrooted neighbour-joining tree containing the de novo reconstructed genes from MG study (genes expressed under Miscanthus diet are highlighted in orange on the tree). Tree was cut in two parts along the dashed line. All GH5_4 characterised proteins were retrieved from the CAZY database and included on the tree. Clusters indicated with an arrow and designated as “MA” contain known multi-functional enzymes. The percentage of replicate trees in which the associated sequences clustered together in the bootstrap test (500 replicates) are shown next to the branches. Final alignment involved 157 amino acid sequences. Protein from the Spirochaetes cluster IX indicated with a grey arrow was heterologously produced and characterised. b Activity profiles for the heterologously produced and purified protein tested against CMC, glucomannan (galactomannan was negative, not displayed on the graph), xylan and arabinoxylan. Boxes represent the interquartile range and error bars show the 95% confidence intervals (n = 4). c Optimal temperature was assessed for glucomannan substrate. Error bars represent the standard deviation of a dataset (n = 3). MAGs reconstruction and carbohydrate utilisation gene clus- others alpha glucans, (arabino)xylans (AX), beta glucans, cel- ters. Reconstruction of draft microbial genomes enriched in lulose, chitin, galactans, mannans and xyloglucans were evi- Miscanthus-fed termite gut microbiome resulted in 20 MAGs denced in reconstructed Spirochaetae genomes. Largely with completeness above 50% and contamination below 10%, complete AX-targeting clusters were assigned to MAG_17 and including eight Fibrobacteres-assigned MAGs, six Spirochaetae MAG_1, both representing Treponema. All of them were (four of Treponema origin) and three novel MAGs representing transcriptionally up-regulated under Miscanthus diet (Supple- Proteobacteria phylum (Supplementary Data 7 and Supple- mentary Fig. 15), what is consistent with high AX content of mentary Fig. 14). Average nucleotide identity to MAGs from Miscanthus hemicellulose. Genes encoding for carbohydrate gut microbiomes of several higher termite species22 equalled transporters (Supplementary Data 8) and in some cases even roughly 77.1 ± 1.8%, confirming the novelty of our MAGs. the two-component system response regulators and chemotaxis Frequency of cazymes was higher in the Treponema genomes were adjacent to several cazymes clusters. Organisationally and they were also enriched in GHs. Over 36% of annotated similar to polysaccharide utilisation loci (PULs) systems cazymes were aggregated into 1096 gene clusters containing employed by Bacteroidetes24, cazymes clusters reconstructed in more than one CAZy encoding gene. Similar gene clusters were our study more resemble the concept of Gram-positive PULs recently discovered in gut microbiota of a wood-feeding termite recently proposed by Sheridan et al.25 in the context of human Globitermes brachycerastes23. Putative cellulose-utilisation gene gut Firmicutes. Taking into account high sequence similarity of clusters were the most highly expressed in Fibrobacteres MAGs, Spirochaeteae and Firmicutes cazymes, it is possible that whole suggesting its major contribution to cellulose degradation cazymes clusters were acquired from Firmicutes in the course of (Figs. 3 and 4). Endoxylanase-encoding genes were often co- evolution. Although not yet shown for Spirochaetae, Bacter- localised with endocellulases (Supplementary Data 7); however, oidetes PULs are often encoded within integrative and con- xylosidases were scarce, questioning the ability of Fibrobacteres jugative elements enabling their transfer among closely related to utilise complex xylans. Cazymes clusters targeting among species17. COMMUNICATIONS BIOLOGY | (2020)3:275 | https://doi.org/10.1038/s42003-020-1004-3 | www.nature.com/commsbio 7 ARTICLE COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-020-1004-3 Host functional gene expression profiles under Miscanthus and goat faeces34. In insects, GH45 cazymes were previously diet. Following taxonomic annotation, 16,416 gene transcripts of detected in the genomes of Phytophaga beetles35 and until now they eukaryotic origin were classified to Arthropoda. Based on the were not reported from sequenced termite genomes. presence of 271 conserved orthologous reference eukaryotic genes26, we estimated the completeness of our termite gut epi- thelial transcriptome (relative to midgut and hindgut) at 89.5%, Diet on Miscanthus: who does what? Miscanthus biomass is while the contamination with foreign mRNA was below 0.7%. mainly composed of cellulose (41.4 ± 2.9%), hemicelluloses Number of assembled ORFs was in line with the two published (25.8 ± 5.2%) including arabinoxylans, xyloglucans, β-glucans, termite genomes, Macrotermes natalensis (16,140 protein-coding and lignin (21.4 ± 3.6%) and other trace components4. Based on genes; ref. 27) and Zootermopsis nevadensis (15,459 protein cod- the annotation of GH profiles and provided there is virtually no ing-genes; ref. 28; Supplementary Data 9 and Supplementary foreign nucleic acid contamination in the reconstructed termite Fig. 16). However, it was much lower than for another sequenced transcriptome, termite on its own could digest amylose (starch), lower termite genome of Cryptotermes secundus (26,726 protein cellulose and/or cellobiose, lactose, galactose, chitin, mannan (α- coding-genes; accession number PRJNA432597). Similar to gut mannan presumably contained in the fungal cell wall) and microbiome, gene transcripts related to carbohydrate transport mannose, trehalose, other glycans (e.g. N-acylsphingosine) and and metabolism were abundant in the host transcriptome, bacterial cell wall components (Fig. 4a, b). Until now, hemi- showing importance of carbohydrates metabolism to the termite cellulose degradation seems to be conducted by gut microbes, and lifestyle (Supplementary Data 3). Over 10,000 transcripts were no putative hemicellulolytic genes were recovered from the assigned a KO number and a summary of complete and partially reconstructed termite transcriptome. For comparison, recent reconstructed metabolic pathways is provided in Supplementary discovery of multi-functional GH9 cellulases in Phasmatodea Data 4. Except for a previously partially reconstructed wood- suggested that some insects are capable to target xylan and feeding Nasutitermes takasagoensis transcriptome with 10,910 xyloglucan in addition to cellulose36. The same study evidenced detected transcripts, our reconstruction represents the most multi-functionality for two GH9 cellulases from Mastotermes complete transcriptome of a higher termite representing Nasuti- darwiniensis, suggesting that some endogenous GH9 members in termitinae subfamily29. termites might also hydrolyse hemicellulose. Slightly increased We further identified 170 CAZY gene transcripts assigned to transcriptional abundance of laccase-coding gene might indicate four main classes (GH, CE, AA and GT) and associated CBMs. termite ability to target lignin. Beyond relatively high lignin Glycoside hydrolases were encoded by 66 genes, and their content, recalcitrance of Miscanthus biomass is mainly enhanced diversity patterns were similar to those identified in other by other features, in particular acetylation and esterification with termites with sequenced genomes (Supplementary Data 9). They ferulic acid (FA). While acetyl esterase activity can be deduced were assigned to 19 different GH families, out of which five were from Fibrobacteres and Spirochaetae metatranscriptomes, the not represented in the gut microbiome (Fig. 3a–c). The highest ability to break ferulic linkages seems limited to Spirochaetae. transcriptional abundance was attributed to GH13 (typically Putative feruoyl esterase from CE1 family are contained within assigned as alpha glucanases; Supplementary Fig. 17), and it AX-targeting cazymes clusters in Spirochaetae MAGs, and were slightly decreased only towards the end of Miscanthus campaign. all up-regulated under Miscanthus diet (Supplementary Fig. 15). Based on rather constant expression profiles of chitodextrinase, Based on the diversity and expression patterns of GHs, different chitin utilisation by the host did not change significantly upon sugar transporters (mainly ABC and to a lesser extent PTS; Miscanthus feeding (Supplementary Fig. 9). This indicates Supplementary Data 8) and specific sugar isomerases, Spirochae- constant complementation of diet with nitrogen rich chitin, tae are able to utilise a wider range of Miscanthus-derived sugars originating from either necrophagy and/or cannibalism, or fungi- (including glucose, glucoronate, rhamnose, arabinose, mannose, colonised food stored in the nest, as previously proposed for other xylose, ribose and fucose) than Fibrobacteres (mainly glucose, higher termite species30. Transcriptional abundance of endocel- mannose and possibly ribose). Both bacterial phyla can target the lulases increased at later stages of Miscanthus feeding, suggesting backbone of cellulose, xylans and mannans (the latter is not a shift towards increased cellulose utilisation by host. In addition, abundant in Miscanthus biomass). Enrichment of Fibrobacteres a gene transcript sharing 54% of identity (at protein level) with a metatranscriptome in endoglucanases (both targeting 1–3 β and newly characterised cellulose- and chitin-targeting AA15 LPMO 1–4 β glycosidic bonds as present in β-glucans and cellulose, from Thermobia domestica31 was also identified. This insect has a respectively) shows its preference for carbohydrates with a remarkable ability to digest crystalline cellulose without microbial glucose-unit backbone. Fibrobacteres also express endoxylanases, assistance. Further blast analysis revealed a presence of homo- and we could confirm experimentally xylanase activity for one logous genes in the other termite genomes, including Z. GH11 enzyme (Supplementary Data 6). However, hardly nevadensis, M. natalensis and C. secundus, suggesting that next represented xylosidase-assigned gene transcripts and the absence to certain eukaryotes (e.g. crustaceans, molluscs, chelicerates, of any xylose transporters and other known genes involved in algae, and oomycetes) termites might be able to oxidatively cleave xylose utilisation would question the ability of Fibrobacteres to glycosidic bonds. Similar observation was also recently stated32. utilise xylans. Co-localisation of many endoxylanases together Interestingly, a few gene transcript of eukaryotic origin were not with potential endocellulases in the reconstructed Fibrobacteres assigned to Arthropoda, suggesting the presence of active small MAGs further suggest that termite gut Fibrobacteres mainly eukaryotes in the higher termite gut. One gene transcript was remove xylan polymers from Miscanthus fibres to better expose assigned to AA3 family and classified as putative cellobiose cellulose to the action of own endocellulases (Supplementary oxidoreductase (EC 1.1.99.18). This type of oxidase is involved in Fig. 13). By contrast, xylose isomerases and xylulose kinases were oxidative cellulose and lignin degradation in wood-decaying enriched in Spirochaetae metatranscriptome and both were highly fungi33. In addition, two other transcripts assigned to GH45 family expressed under Miscanthus diet. All reconstructed gene tran- (putative endoglucanases) were present in the reconstructed scripts were assigned to Spirochaetae, and together with enrich- transcriptome, and their expression slightly increased under ment of endoxylanase transcripts, it confirms the ability of these Miscanthus diet. Homology search (blastp) revealed their closest bacteria to degrade xylans, as recently proposed by Tokuda et al.9. similarity to Anaeromyces contortus, sp. nov. (Neocallimastigomy- Based on the in silico prediction of enzyme sub-cellular cota), an anaerobic gut fungal species recently isolated from cow localisations, most of the endoxylanases from Fibrobacteres and 8 COMMUNICATIONS BIOLOGY | (2020)3:275 | https://doi.org/10.1038/s42003-020-1004-3 | www.nature.com/commsbio COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-020-1004-3 ARTICLE Spirochaetae are exported outside the cell (Supplementary esterases, by removing cross-links between polysaccharides and Data 10), suggesting initial degradation of xylan backbone in lignin help separating lignin from the rest of biomass, offering an the extracellular space. Many Spirochaetae endocellulases are also alternative and/or complementation to currently applied indus- putative extracellular enzymes or anchored to the outer trial treatments45. In addition, ferulic acid and other hydro- membrane. In contrast, multiple Fibrobacteres endocellulases xycinnamic acids can have many applications in food and possibly lack signal peptide and are assumed to be localised in the cosmetic industries due to their antioxidant properties etc.46, thus cytoplasm. In general, as much as half of Fibrobacteres GHs are further extending the application range of Miscanthus biomass. predicted to be localised in cytoplasm. At the same time, three Approach-wise, experimental design undertaken in this study times more GHs are exported outside the cell by Spirochaetae. represents the enrichment strategy where a nature-derived This could indicate a rather selfish carbohydrates degradation microbial inoculum is grown in liquid batch cultures. Here, a strategy employed by termite Fibrobacteres, where cellulose fibres natural system of the termite gut was shown to progressively primarily detached by extracellular endocellulases are transported adapt yielding a consortium of microbes specialised in degrada- inside the cell for further breakdown. Selfish carbohydrate capture tion of Miscanthus biomass. Integrative omics combined with and degradation was previously proposed for Bacteroidetes in the protein characterisation provides a framework for better under- rumen37 and anaerobic digestion reactors38. By contrast, recent standing of complex lignocellulose degradation by the higher work done on a rumen isolate Fibrobacter succinogenes S85 termite gut system and paves a road towards its future indicates that enzymes involved in cellulose degradation are bioprospecting. localised on the cell surface39. By maximising intracellular cellulose breakdown termite gut Fibrobacteres would avoid being in competition with much more abundant Spirochaetae. Enrich- Methods Nest origin, laboratory maintenance and sampling. Initially, three colonies of ment of exbB (encoding for a biopolymer transport protein ExbB; grass-feeding higher termites from Nasutitermitinae (nests: LM1, LM2 and LM3) Fig. 2a) gene transcripts in Fibrobacteres metatranscriptome were identified in January 2017 in tropical savannah in French Guiana, in proxi- would suggest possible cellulose/cellodextrin uptake through a mity to Sinnamary town (radius of 5 km to GPS: N 05°24.195′ W 053°07.664′). mechanism similar to an experimentally demonstrated TonB- Termite nests were transported to the laboratory where colonies were maintained in separate glass containers at 26 °C, 12 h light and 12 h dark, and 90% humidity dependent transport of maltodextrins across outer membrane of conditions. Termite colonies were fed with dried Miscanthus grass winter harvest Caulobacter crescentus40, also discussed for F. succinogenes39. rich in recalcitrant lignocellulose47, for a period of up to 9 months. Colony LM1_2 died after few months and was excluded from extended analysis. Mature worker- caste individuals were sampled in regular monthly time intervals before (sample taken before Miscanthus diet corresponds to “wild-microbiome”) and following Conclusions. Retrieving lignocellulose-active enzymes from Miscanthus diet (samples correspond to “Miscanthus-adapted microbiome”). naturally evolved biomass-degrading systems, with the use of Termite specimens were cold immobilised, surface-cleaned with 80% ethanol and continuously improving high-throughput sequencing technolo- 1× phosphate-buffered saline and decapitated. Whole guts (here relative to the gies, presents a promising strategy to identify new enzymes with midgut and hindgut compartments) were dissected (n ≈ 30 per replicate, minimum potentially enhanced activities. Limited metatranscriptomic three replicates per sample) and preserved directly in liquid nitrogen. Additionally, for a sample selected for metagenomics analysis (LM1 time point 8 months; reports highlight high representation and overexpression of LM1_8), the hindgut luminal fluid was collected as previously described11. Samples CAZymes in termite digestomes (e.g. refs. 9,14). In higher termite were stored at −80 °C until further processing. Termite species were identified by guts, many lignocellulolytic steps are assisted by gut micro- morphology and by sequencing of the partial COII marker gene, as described organisms (e.g. cellulose deconstruction), while some are exclu- before48. sively attributed to hindgut bacteria (e.g. hemicellulose degradation). Cellulose degradation capacities of different lig- Extraction of nucleic acids. Total DNA and RNA were co-extracted from all nocellulose degrading environments, including the termite gut samples using the AllPrep PowerViral DNA/RNA Kit (Qiagen) following the system41, have extensively been studied in the past. Decomposi- manufacturer’s protocol. To assure the proper disruption of bacterial cell wall and termite gut epithelium cells, mechanical bead-beating step with 0.1 mm glass beads tion of hemicellulose and general mobilisation of different lig- at 20 Hz for 2 min was introduced to complement the chemical lysis. The eluents nocellulose components (breaking bonds between diverse plant were divided in two parts. First part was treated with 1 µL of 10 µg/mL RNase A polymers) have received comparably less scientific focus. (Sigma) for 30 min at room temperature. The second part was treated with TURBO Importantly, there is an increasing industrial interest in xylan- DNA-free kit (Invitrogen) according to the manufacturer’s protocol. The resulting processing enzymes42, regarding their application in biomass pure DNA and RNA fractions were quality assessed using agarose gel electro- phoresis and Bioanalyser RNA 6000 Pico Kit (Agilent). Nucleic acid concentration (wood) processing, pulp bio-bleaching, animal nutrition, food was quantified using Qubit dsDNA HS Assay and Qubit RNA HS Assay Kit additives, etc.43. (Invitrogen). DNA and RNA were stored at −20 °C and −80 °C, respectively. According to the recently established glycome profile of Miscanthus sp., next to glucose, xylose and arabinose are the 16S rRNA gene amplicon high-throughput sequencing and data analysis. The two main cell wall monosaccharides, both originating from bacterial 16S rRNA gene amplicon libraries were prepared using Illumina com- arabinoxylan fibres which are esther-cross-linked by ferulic patible approach as previously described49. Briefly, modified universal primers S- acid44. Consequently, cazymes specifically targeting (feruloylated) D-Bact-0909-a-S-18 (ACTCAAAKGAATWGACGG) and S-*-Univ-*-1392-a-A- arabinoxylan components were highly up-regulated under 15 (ACGGGCGGTGTGTRC, ref. 50), and Nextera XT Index Kit V2 (Illumina) Miscanthus diet, making them potentially suitable candidates were used along with Q5 Hot Start High-Fidelity 2× Master Mix (New England Biolabs) in a two-step polymerase-chain reaction (PCR). In the first step selective for industrial xylan-targeting applications. Many of these genes amplification of the 484 bp long fragments of bacterial 16S rRNA gene V6–V8 were combined in clusters with a set of complementary hydrolytic region was performed. llumina-compatible adapters and barcodes were attached in activities to degrade e.g. AXs. Nature-optimised synergy between the second step. Purified and equimolarly pooled libraries were sequenced along enzymes of the same CAZyme cluster could further provide the with PhiX control (Illumina) using MiSeq Reagent Kit V3-600 on in-house Illu- mina MiSeq Platform. Usearch v.7.0.1090_win64 software51 was used for quality basis to better define industrially relevant enzymatic cocktails. trimming, chimera check, singletons removal and assignment of the obtained Specific lignocellulose fractions could be selectively targeted to sequences to OTUs at 97% similarity level. Taxonomic affiliation of the resulting deliver desired products, with potential effects being fully OTUs was performed with SILVA database v.128 (ref. 52). Downstream analyses expectable and controllable38. It would allow for their fine- were performed with mothur53 and R environment54. Bacterial richness and diversity were calculated using respectively sobs and invsimpson indices. The tuned degradation for a variety of applications, e.g. oligoarabi- dissimilarity of bacterial community structures was calculated using the Bray noxylans for food industry (prebiotics), lignin fibres for −Curtis index. OTUs differentially abundant between the wild- and Miscanthus- biomaterials, glucomannans as food additives, etc. Feruloyl adapted microbiomes were assessed using the LEfSe approach12. COMMUNICATIONS BIOLOGY | (2020)3:275 | https://doi.org/10.1038/s42003-020-1004-3 | www.nature.com/commsbio 9 ARTICLE COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-020-1004-3 De novo metagenomics and data analysis. Sample LM1_8 was selected for detected using dbCAN (dbCAN-fam-HMMs.txt.v6; ref. 65) and CAZy database66. metagenomic sequencing in order to reconstruct genomes/larger genomic frag- Using the threshold of e-value of <1e−18 and coverage >0.35 recommended for ments of the dominant microbes in the Miscanthus-adapted microbiome. Pro- prokaryotic CAZymes resulted in the removal of a high number of putative karyotic DNA was enriched from the total hindgut DNA extract using NEBNext CAZymes; therefore, additional manual curation was performed to maximise the Microbiome DNA Enrichment Kit (New England BioLabs). Following sequencing, number of entries retained for further analysis. Additionally, gene transcripts over 170 Mbp raw reads were quality trimmed in CLC Genomics Workbench outliers (very partially reconstructed gene fragments with average MT expression v.9.5.2 (Qiagen), using a phred quality score of 20, minimum length of 50 and significantly exceeding the average expression of other genes assigned to the same allowing no ambiguous nucleotides, resulting in close to 154 million quality- group) were manually identified and removed as they were considered as chimeric trimmed paired reads. Quality-trimmed reads were assembled using the CLC’s de (additional Blast search was launched in each case). Homology to peptide pattern novo assembly algorithm in a mapping mode, using automatic bubble size and (Hotpep) was used to assign an EC class to the identified CAZymes (ref. 21). Sub- word size, minimum contig length of 1000, mismatch cost of 2, insertion cost of 3, cellular localisation of CAZymes was predicted using BUSCA web67. To link the deletion cost of 3, length fraction of 0.9, and similarity fraction of 0.95. The average degradation of the different lignocellulose fractions and subsequent sugar utilisa- contig abundance was calculated as DNA-RPKMs (reads per kilo base per million tion, we looked for the presence of suitable sugar transporters and also specific mapped reads). This type of normalisation allows for comparing contigs (genomic sugar isomerases and kinases. Eukaryotic CAZymes, including for other sequenced fragments) coverage (abundance) values, as it corrects differences in both sample termite genomes, were further searched with dbCAN2 and using dbCAN-fam- sequencing depth and contig length. ORFs were searched and annotated using the HMMs.txt.v8, with new AA families included. default pipeline integrated in the IMG/MER55. Information about KEGGs and Genes encoding for CAZymes of interest, selected based on their predicted COGs assignment was retrieved based on the IMG/MER annotations. Metabolic activities and their expression profiles, were further PCR amplified (Veriti™ 96 wells pathways/modules were reconstructed using the tool integrated in the online Thermal cycler, Applied Biosystems, Foster City, USA) and the resulting PCR version of the KEGG database56. Initially, contigs were binned using myCC57 what products were purified using a PCR purification kit (Qiagen, Hilden, Germany). If resulted in their separation into 35 phylum-level bins of relatively high con- any signal peptide was predicted (using LipoP version 1.0, http://www.cbs.dtu.dk/ tamination. This approach was undertaken to correctly assign phylum-level tax- services/LipoP/, ref. 68), it was removed before cloning to enhance cytoplasmic onomy to the resulting contigs. Subsequently, the bin refinement module integrated protein production. Purified PCR products were cloned into the pET52b(+) in MetaWRAP was used to fine-tune the resulting bins (MAGs) to the species/ plasmid and expressed in E. coli Rosetta (DE3) strain (Millipore Corporation, strain levels58. The completeness and contamination of the generated MAGs were Billerica, MA, USA), as previously described38. Cells were harvested by assessed with checkM59. Taxonomic affiliation was assessed with PhyloPhlan60. centrifugation (5000 × g, 4 °C, 15 min) and re-dissolved in a lysis buffer (50 mM Similarity to the previously reconstructed MAGs was verified with the FastANI61. NaH2PO4, 10 mM imidazole, 300 mM NaCl, pH8). Proteins were released by MAGs abundance within the reconstructed metagenome was calculated as average sonication, cell debris were removed by centrifugation (16,000 × g, 4 °C, 15 min) of contigs metagenomic abundance (DNA-RPKMs, see above) assigned to a spe- and subsequent filtration step (13-mm syringe filter, 0.2-µm-pore-size). Affinity tag cific MAG. Given the relatively high microbial diversity in the termite gut, only purification was achieved using a histidine tag located at the C terminus of a 54.6% of the resulting MG sequencing reads could map back to the reconstructed recombinant protein. NGC™ Medium-Pressure Liquid Chromatography system MG contigs, potentially mitigating the rate of functional gene discovery if solely (Bio-Rad) equipped with a HitrapTM column of 5 mL (Bio-Rad, Hercules, CA, relaying on the de novo MG reconstruction. Phylogenetic analyses were performed USA) was used to purify produced proteins. A constant flow rate of 5 mL/min was on MAFFT-aligned protein sequences62 using MEGA X63. applied. Initially, column was equilibrated with six column volumes (CV) of lysis buffer. Following equilibration, 290 mL of sample was injected and washed with three CVs of mixed buffer (97% of lysis buffer and 3% of elution buffer, the latter De novo metatranscriptomics, host transcriptomics and data analysis. For composed of 50 mM NaH2PO4, 250 mM imidazole, 300 mM NaCl, pH8). First step three selected samples (LM1-1, LM1-2 and LM1-8) the (meta)transcriptomic ana- of elution was achieved using ten CVs of a linear gradient of mixed buffer (from 3% lysis was performed using an optimised approach described earlier11. Ribo-Zero of elution buffer plus 97% of lysis buffer to 50% of elution buffer plus 50% of lysis Gold rRNA Removal Kit “Epidemiology” (Illumina) was used to enrich the sample buffer), followed by four CVs of 100% of elution buffer and finally one CV of 100% for prokaryotic and eukaryotic mRNA. Enriched mRNA was purified using Agen- lysis buffer to detach the remaining protein. Fractions of two milliliters were court RNAClean XP Kit and analysed with Bioanalyser RNA 6000 Pico Kit (Agilent). collected during the washing and elution steps, and were analysed on SDS-PAGE. In continuation, SMARTer Stranded RNA-Seq Kit (Clontech) was used according to The release of 4-nitrophenol (PNP assay) and/or reducing sugar (RS assay) was the manufacturer’s instructions to prepare sequencing libraries. Final libraries were used to determine the activity of recombinant proteins. Briefly, 50 µL of a purified quantified with High Sensitivity DNA Kit (Agilent) and KAPA SYBR FAST Uni- protein solution was incubated with 50 µL of substrate (respectively, 4-nitrophenol versal qPCR Kit. Libraries were pair-end sequenced at the Luxembourg Centre for derivatives were used for PNP assay and polysaccharides for RS assay) and 100 µL Systems Biomedicine (University of Luxembourg) using Illumina NextSeq 500/550 (PNP assay) or 25 µL (RS assay) of citrate phosphate buffer (pH7 0.1 M citric acid, High Output v2-300 Kit. Over 285 million raw reads were quality trimmed in CLC 0.2 M dibasic sodium phosphate). The targeted substrates included Genomics Workbench v.9.5.2, using a phred quality score of 20, minimum length of carboxymethylcellulose (CMC), arabinoxylan, galactomannan, glucomannan and 50 and allowing no ambiguous nucleotides, resulting in close to 214 million quality- xylan. Enzymatic reaction was carried out at 37 °C during 1 hour (PNP assay) or trimmed reads. Contaminating rRNA reads were removed using SortMeRNA 30 min (RS assay). The rate of release of 4-nitrophenol was instantly monitored at 2.0 software64. The resulting non-rRNA reads were used to perform de novo (meta) 405 nm using SPECORD 250 PLUS (Analytic Jena). The release of reducing sugars transcriptome co-assembly using the CLC assembly algorithm in mapping mode was determined following the Somogyi-Nelson method (refs. 69,70). All assays were with default parameters, except for minimum contig length of 200, length fraction of performed in triplicates. 0.90 and similarity fraction 0.95. As a result, nearly 2 million contigs were assembled. Obtained contigs were further submitted to IMG/MER for taxonomic and functional annotation55. Following the taxonomic assignment, 759,451 transcripts of putative prokaryotic origin were selected for further analysis. Initial IMG/MER taxonomic Statistics and reproducibility. Whenever relevant biological or technical repli- annotation resulted in over-representation of transcripts of putative Firmicutes cates were included in our study, this information is provided in specific sections of origin (Supplementary Fig. 4). As nearly no Firmicutes OTUs were detected using the “Methods” chapter. All statistical tests used are indicated and the reference is the 16S rRNA gene amplicon sequencing, transcripts were re-annotated based on the provided. Correlation was calculated using Pearson correlation coefficient. de novo assembled metagenome and contig binning, resulting in re-classification of virtually all Firmicutes-assigned contigs to Fibrobacteres and Spirochaetae. To complement the study and to characterise potential contribution of the termite host Reporting summary. Further information on research design is available in to Miscanthus digestion, transcripts of eukaryotic origin and taxonomically assigned the Nature Research Reporting Summary linked to this article. to Insecta (based on the IMG/MER annotation) were further evaluated for the completeness of the de novo reconstructed transcriptome with the BUSCO pipe- line26 and using the Eukaryota database (odb9). There were only two duplicated Data availability genes out of the total of 303 searched groups, suggesting that the level of possible The raw sequencing data used in this study were deposited in the SRA database under the contamination of non-Insecta origin was very low (below 0.7%). following project numbers: PRJNA587606 (16S rRNA amplicon sequencing), For both the de novo assembled metatranscriptome and termite transcriptome, PRJNA587423 (Metagenomics) and PRJNA587406 (Metatranscriptomics). OTU in order to determine the relative abundances of transcripts across studied samples, sequences were deposited in GenBank under project numbers PRJNA586754 and sequencing reads were mapped back to the annotated transcript sets, using the CLC PRJNA434195. COII nucleotide sequences are available in GenBank under accession “RNA-seq analysis” mode, with default parameters except for minimum similarity numbers MN803317−19. MG and MT assemblies and all other data underlying the of 0.95 over 0.9 of the read length, both strands specificity and 1 maximum number findings of this study are available from the corresponding author upon reasonable of hits per read. The mapping results were represented as TPMs (transcripts per request. million), what directly resulted in normalised reads counts. Received: 19 December 2019; Accepted: 13 May 2020; Identification of cazymes and enzymes characterisation. Genes (metagenomics) and gene transcripts (metatranscriptomics) encoding for microbial CAZymes were 10 COMMUNICATIONS BIOLOGY | (2020)3:275 | https://doi.org/10.1038/s42003-020-1004-3 | www.nature.com/commsbio COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-020-1004-3 ARTICLE References takasagoensis (Blattodea (Isoptera): Termitidae). Appl. Entomol. Zool. 51, 1. Bar-On, Y. M., Phillips, R. & Milo, R. The biomass distribution on Earth. Proc. 429–440 (2016). Natl. Acad. Sci. USA 115, 6506–6511 (2018). 30. Menezes, L. et al. Food storage by the savanna termite Cornitermes cumulans 2. Lewandowski, I., Clifton-Brown, J. C., Scurlock, J. M. O. & Huisman, W. (Syntermitinae): a strategy to improve hemicellulose digestibility? Microb. Miscanthus: European experience with a novel energy crop. Biomass-. Ecol. 76, 492–505 (2018). Bioenergy 19, 209–227 (2000). 31. Sabbadin, F. et al. An ancient family of lytic polysaccharide monooxygenases 3. Luo, H. et al. Total utilization of Miscanthus biomass, lignin and with roles in arthropod development and biomass digestion. Nat. Commun. 9, carbohydrates, using Earth abundant nickel catalyst. ACS Sustain. Chem. Eng. 756–756 (2018). 4, 2316–2322 (2016). 32. Tokuda, G. in Advances in Insect Physiology Vol. 57 (ed. Jurenka, R.) 97−136 4. van der Weijde, T., Dolstra, O., Visser, R. G. F. & Trindade, L. M. Stability of (Academic Press, 2019). cell wall composition and saccharification efficiency in Miscanthus across 33. Berlemont, R. Distribution and diversity of enzymes for polysaccharide diverse environments. Front. Plant Sci. https://doi.org/10.3389/fpls.2016.02004 degradation in fungi. Sci. Rep. 7, 222 (2017). (2017). 34. Hanafy, R. A., Johnson, B., Elshahed, M. S. & Youssef, N. H. Anaeromyces 5. Cantarel, B. L. et al. The Carbohydrate-Active enZymes database (CAZy): an contortus, sp. nov., a new anaerobic gut fungal species expert resource for Glycogenomics. Nucleic Acids Res. 37, D233–D238 (2008). (Neocallimastigomycota) isolated from the feces of cow and goat. Mycologia 6. Brune, A. Symbiotic digestion of lignocellulose in termite guts. Nat. Rev. 110, 502–512 (2018). Microbiol. 12, 168 (2014). 35. Busch, A., Danchin, E. G. J. & Pauchet, Y. Functional diversification of 7. Donovan, S. E., Eggleton, P. & Bignell, D. E. Gut content analysis and a new horizontally acquired glycoside hydrolase family 45 (GH45) proteins in feeding group classification of termites. Ecol. Entomol. 26, 356–366 (2001). Phytophaga beetles. BMC Evolut. Biol. 19, 100 (2019). 8. Tartar, A. et al. Parallel metatranscriptome analyses of host and symbiont gene 36. Shelomi, M., Wipfler, B., Zhou, X. & Pauchet, Y. Multifunctional cellulase expression in the gut of the termite Reticulitermes flavipes. Biotechnol. Biofuels enzymes are ancestral in Polyneoptera. Insect Mol. Biol. 29, 124–135 2, 25 (2009). (2020). 9. Tokuda, G. et al. Fiber-associated spirochetes are major agents of 37. Solden, L. M. et al. Interspecies cross-feeding orchestrates carbon degradation hemicellulose degradation in the hindgut of wood-feeding higher termites. in the rumen ecosystem. Nat. Microbiol. 3, 1274–1284 (2018). Proc. Natl. Acad. Sci. USA 115, E11996–E12004 (2018). 38. Bertucci, M. et al. Carbohydrate hydrolytic potential and redundancy of an 10. Cuezzo, C., Carrijo, T. F. & Cancello, E. M. Transfer of two species from anaerobic digestion microbiome exposed to acidosis, as uncovered by Nasutitermes Dudley to Cortaritermes Mathews (Isoptera: Termitidae: metagenomics. Appl. Environ. Microbiol. 85, e00895-00819 (2019). Nasutitermitinae). Austral. Entomol. 54, 172–179 (2015). 39. Raut, M. P., Couto, N., Karunakaran, E., Biggs, C. A. & Wright, P. C. 11. Marynowska, M. et al. Optimization of a metatranscriptomic approach to Deciphering the unique cellulose degradation mechanism of the ruminal study the lignocellulolytic potential of the higher termite gut microbiome. bacterium Fibrobacter succinogenes S85. Sci. Rep. 9, 16542 (2019). BMC Genomics 18, 681–681 (2017). 40. Neugebauer, H. et al. ExbBD-dependent transport of maltodextrins through 12. Segata, N. et al. Metagenomic biomarker discovery and explanation. Genome the novel MalA protein across the outer membrane of Caulobacter crescentus. Biol. 12, R60 (2011). J. Bacteriol. 187, 8300–8311 (2005). 13. Heintz-Buschart, A. & Wilmes, P. Human gut microbiome: function matters. 41. Brune, A. in Beneficial Microorganisms in Multicellular Life Forms (eds Trends Microbiol. 26, 563–574 (2018). Rosenberg, E. & Gophna, U.) 3−25 (Springer Berlin Heidelberg, 2011). 14. He, S. et al. Comparative metagenomic and metatranscriptomic analysis of 42. Sharma, M. Xylanases: an overview. Br. Biotechnol. J. 3, 1–28 (2013). hindgut paunch microbiota in wood- and dung-feeding higher termites. PLoS 43. Walia, A., Guleria, S., Mehta, P., Chauhan, A. & Parkash, J. Microbial ONE 8, e61126 (2013). xylanases and their industrial application in pulp and paper biobleaching: a 15. Rossmassler, K. et al. Metagenomic analysis of the microbiota in the highly review. 3 Biotech 7, 11 (2017). compartmented hindguts of six wood- or soil-feeding higher termites. 44. da Costa, R. M. F. et al. A cell wall reference profile for Miscanthus bioenergy Microbiome 3, 56–56 (2015). crops highlights compositional and structural variations associated with 16. Hu, Y. et al. Herbivorous turtle ants obtain essential nutrients from a development and organ origin. N. Phytol. 213, 1710–1725 (2017). conserved nitrogen-recycling gut microbiome. Nat. Commun. 9, 964 (2018). 45. Fazary, A. E. & Ju, Y.-H. Feruloyl esterases as biotechnological tools: current 17. Coyne, M. J., Zitomersky, N. L., McGuire, A. M., Earl, A. M. & Comstock, L. and future perspectives. Acta Biochim. Biophys. Sin. 39, 811–828 (2007). E. Evidence of extensive DNA transfer between Bacteroidales species within 46. Dilokpimol, A. et al. Diversity of fungal feruloyl esterases: updated the human gut. mBio 5, e01305–e01314 (2014). phylogenetic classification, properties, and industrial applications. Biotechnol. 18. Parks, D. H. et al. Recovery of nearly 8,000 metagenome-assembled genomes Biofuels 9, 231 (2016). substantially expands the tree of life. Nat. Microbiol. 2, 1533–1542 (2017). 47. Godin, B. et al. Chemical characteristics and biofuel potential of several vegetal 19. Franco Cairo, J. P. et al. Expanding the knowledge on lignocellulolytic and biomasses grown under a wide range of environmental conditions. Ind. Crops redox enzymes of worker and soldier castes from the lower termite Products 48, 1–12 (2013). Coptotermes gestroi. Front. Microbiol. 7, 1518 (2016). 48. Miura, T., Roisin, Y. & Matsumoto, T. Molecular phylogeny and biogeography 20. Cragg, S. M. et al. Lignocellulose degradation mechanisms across the Tree of of the nasute termite genus Nasutitermes (Isoptera: Termitidae) in the Pacific Life. Curr. Opin. Chem. Biol. 29, 108–119 (2015). Tropics. Mol. Phylogenet. Evolut. 17, 1–10 (2000). 21. Busk, P. K., Pilgaard, B., Lezyk, M. J., Meyer, A. S. & Lange, L. Homology to 49. Calusinska, M. et al. A year of monitoring 20 mesophilic full-scale bioreactors peptide pattern for annotation of carbohydrate-active enzymes and prediction reveals the existence of stable but different core microbiomes in bio-waste and of function. BMC Bioinform. 18, 214 (2017). wastewater anaerobic digestion systems. Biotechnol. Biofuels 11, 196 (2018). 22. Hervé, V. et al. Phylogenomic analysis of 589 metagenome-assembled 50. Klindworth, A. et al. Evaluation of general 16S ribosomal RNA gene PCR genomes encompassing all major prokaryotic lineages from the gut of higher primers for classical and next-generation sequencing-based diversity studies. termites. PeerJ Prepr. 7, e27929v27921 (2019). Nucleic Acids Res. 41, e1 (2013). 23. Liu, N. et al. Functional metagenomics reveals abundant polysaccharide- 51. Edgar, R. C. Search and clustering orders of magnitude faster than BLAST. degrading gene clusters and cellobiose utilization pathways within gut Bioinformatics 26, 2460–2461 (2010). microbiota of a wood-feeding higher termite. ISME J. 13, 104–117 (2019). 52. Yilmaz, P. et al. The SILVA and “all-species Living Tree Project (LTP)” 24. Zhang, M. et al. Xylan utilization in human gut commensal bacteria is taxonomic frameworks. Nucleic Acids Res. 42, D643–D648 (2014). orchestrated by unique modular organization of polysaccharide-degrading 53. Schloss, P. D. et al. Introducing mothur: open-source, platform-independent, enzymes. Proc. Natl. Acad. Sci. USA 111, E3708–E3717 (2014). community-supported software for describing and comparing microbial 25. Sheridan, P. O. et al. Polysaccharide utilization loci and nutritional communities. Appl. Environ. Microbiol. 75, 7537–7541 (2009). specialization in a dominant group of butyrate-producing human colonic 54. Team, R. C. R: A Language and Environment for Statistical Computing (R Firmicutes. Microb. Genomics 2, e000043 (2016). Foundation for Statistical Computing, 2019). 26. Simão, F. A., Waterhouse, R. M., Ioannidis, P., Kriventseva, E. V. & Zdobnov, 55. Chen, I.-M. A. et al. IMG/M v.5.0: an integrated data management and E. M. BUSCO: assessing genome assembly and annotation completeness with comparative analysis system for microbial genomes and microbiomes. Nucleic single-copy orthologs. Bioinformatics 31, 3210–3212 (2015). Acids Res. 47, D666–D677 (2018). 27. Poulsen, M. et al. Complementary symbiont contributions to plant 56. Kanehisa, M. & Goto, S. KEGG: Kyoto Encyclopedia of Genes and Genomes. decomposition in a fungus-farming termite. Proc. Natl. Acad. Sci. USA 111, Nucleic Acids Res. 28, 27–30 (2000). 14500–14505 (2014). 57. Lin, H.-H. & Liao, Y.-C. Accurate binning of metagenomic contigs via 28. Terrapon, N. et al. Molecular traces of alternative social organization in a automated clustering sequences using information of genomic signatures and termite genome. Nat. Commun. 5, 3636 (2014). marker genes. Sci. Rep. 6, 24175 (2016). 29. Kumara, R. P., Saitoh, S., Aoyama, H., Shinzato, N. & Tokuda, G. Metabolic 58. Uritskiy, G. V., DiRuggiero, J. & Taylor, J. MetaWRAP—a flexible pipeline for pathways in the mixed segment of the wood-feeding termite Nasutitermes genome-resolved metagenomic data analysis. Microbiome 6, 158 (2018). COMMUNICATIONS BIOLOGY | (2020)3:275 | https://doi.org/10.1038/s42003-020-1004-3 | www.nature.com/commsbio 11 ARTICLE COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-020-1004-3 59. Parks, D. H., Imelfort, M., Skennerton, C. T., Hugenholtz, P. & Tyson, G. W. Author contributions CheckM: assessing the quality of microbial genomes recovered from isolates, Experiments and molecular analyses were planned by M.C., M.M., M.B. and D.S.-D. and single cells, and metagenomes. Genome Res. 25, 1043–1055 (2015). carried out by M.C., M.M. and M.B. with technical support from B.U. and D.K. M.C. and 60. Segata, N., Börnigen, D., Morgan, X. C. & Huttenhower, C. PhyloPhlAn is a M.M. analysed the results. M.B. performed protein studies. Y.R., D.S.-D., X.G. helped new method for improved phylogenetic and taxonomic placement of with the collection of termites. R.H. and P.W. provided the support with Illumina microbes. Nat. Commun. 4, 2304–2304 (2013). NextSeq sequencing. P.Gawron and X.G. provided support with bioinformatics analysis. 61. Jain, C., Rodriguez-R, L. M., Phillippy, A. M., Konstantinidis, K. T. & Aluru, S. M.C., M.M. and M.B. wrote the manuscript. P.D., Y.R., D.S.-D., P.Gerin, P.F. participated High throughput ANI analysis of 90K prokaryotic genomes reveals clear in the planning and coordination of the study and in the manuscript correction. All species boundaries. Nat. Commun. 9, 5114 (2018). authors read and approved the final manuscript. 62. Katoh, K., Rozewicki, J. & Yamada, K. D. MAFFT online service: multiple sequence alignment, interactive sequence choice and visualization. Brief. Competing interests Bioinform. 20, 1160–1166 (2017). The authors declare no competing interests. 63. Kumar, S., Stecher, G., Li, M., Knyaz, C. & Tamura, K. MEGA X: molecular evolutionary genetics analysis across computing platforms. Mol. Biol. Evolut. 35, 1547–1549 (2018). Additional information Supplementary information is available for this paper at https://doi.org/10.1038/s42003- 64. Kopylova, E., Noé, L. & Touzet, H. SortMeRNA: fast and accurate filtering of 020-1004-3. ribosomal RNAs in metatranscriptomic data. Bioinformatics 28, 3211–3217 (2012). Correspondence and requests for materials should be addressed to M.C. 65. Yin, Y. et al. dbCAN: a web resource for automated carbohydrate-active enzyme annotation. Nucleic Acids Res. 40, W445–W451 (2012). Reprints and permission information is available at http://www.nature.com/reprints 66. Lombard, V., Golaconda Ramulu, H., Drula, E., Coutinho, P. M. & Henrissat, B. The carbohydrate-active enzymes database (CAZy) in 2013. Nucleic Acids Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in Res. 42, D490–D495 (2014). published maps and institutional affiliations. 67. Savojardo, C., Martelli, P. L., Fariselli, P., Profiti, G. & Casadio, R. BUSCA: an integrative web server to predict subcellular localization of proteins. Nucleic Acids Res. 46, W459–w466 (2018). 68. Juncker, A. S. et al. Prediction of lipoprotein signal peptides in Gram-negative bacteria. Protein Sci. 12, 1652–1662 (2003). 69. Somogyi, M. A new reagent for the determination of sugars. Int. J. Biol. Chem. Open Access This article is licensed under a Creative Commons 160, 61–68 (1945). Attribution 4.0 International License, which permits use, sharing, 70. Nelson, N. A photometric adaptation of the Somogyi method for the adaptation, distribution and reproduction in any medium or format, as long as you give determination of glucose. J. Biol. Chem. 152, 375–380 (1944). appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the Acknowledgements article’s Creative Commons license and your intended use is not permitted by statutory This research was funded through FNR 2014 CORE project OPTILYS (Exploring the regulation or exceeds the permitted use, you will need to obtain permission directly from higher termite lignocellulolytic system to optimise the conversion of biomass into energy the copyright holder. To view a copy of this license, visit http://creativecommons.org/ and useful platform molecules/C14/SR/8286517), and co-funded through the grant PDR licenses/by/4.0/. T.0065.15 from the Belgian F.R.S.-FNRS. We are grateful to Philippe Cerdan, Régis Vigouroux and the staff of the Laboratoire Environnement HYDRECO of Petit Saut (EDF-CNEH) for logistic support during the field work. © The Author(s) 2020 12 COMMUNICATIONS BIOLOGY | (2020)3:275 | https://doi.org/10.1038/s42003-020-1004-3 | www.nature.com/commsbio