We performed two in vitro batch fermentation experiments to understand how the IMMPs containing different amounts of α − (1 → 6) glycosidic linkages were broken down by human faecal microorganisms over time, and how the chemical structure of these compounds affected the functional dynamics of the microbial community during fermentation (Fig. 1). Experiment A included fermentation of IMMPs of varying percentage of α − (1 → 6) glycosidic linkages (27%, IMMP-27; 94%, IMMP-94) at three different time points. This was complemented by experiment B that was performed with IMMP with 96% α − (1 → 6) linkages (IMMP-96) and IMMP-27 after treatment with α-amylase and amyloglucosidase (IMMP-dig27). Furthermore, in experiment B an additional set of time points was evaluated to provide a more detailed understanding of microbial community dynamics. In both experiments a control blank without any IMMP substrate was included. We then performed metatranscriptome sequencing of all samples and assembled the resulting data into one reference metatranscriptome. Afterwards, machine learning techniques were applied to identify groups of similarly behaving bacteria and to identify consistent dynamic patterns in gene expression.Figure 1In total two experiments were performed. In Experiment A, 3 time points (0 h, 24 h, and 48 h) were sampled, with two replicates taken at each time point. In experiment B a 6 h and 12 h time point was added to allow closer monitoring of the degradation process. A blank without substrate was also included in experiment B, but with only one replicate.Quality control and statisticsThe metatranscriptomes were sequenced and subjected to a quality control before the data was further analysed (Fig. S1). As a result, 320 million reads (89% of the raw reads and 54% of all bases) passed the quality check and were used for contig assembly. In experiment A, the assembly yielded over 140,000 contigs, with more than 200,000 protein coding genes, and contained, on average, 81% of the input reads (range 71–85%) per sample. Read counts for experiment B were acquired by mapping to the same assembly obtained from experiment A (Table S1), and showed the same average mapping rate (81%, range 71–89%). After mapping, the biological replicates within each experiment showed a spearman correlation of on average 0.86 (range 0.78–0.93), indicating good reproducibility within sets of samples from the same treatment group (Table S2). The six samples from time point 0 showed a correlation of at least 0.79 (max 0.86), and all samples taken from batches incubated in the presence of prebiotics showed a correlation of at least 0.7 for time point 24 h, and all (except two) showed a correlation of at least 0.6 at time point 48 h, indicating similar development over all cultures.Of the ~ 200,000 protein coding genes, ~ 36,000 were full length genes. To ~ 144,000 protein sequences at least one Interpro domain (excluding “Coils” domain) could be identified. These included more than 24,000 genes with at least one (partial) Enzme Commission (EC) number.Community structure and expression patternsTaxonomic classification to at least the superkingdom of Bacteria was assigned to 190,000 of the 200,000 genes obtained from the RNA-assembly. Less than 3000 genes were assigned to eukaryotes and less than 2000 to Archaea. Of the bacterial groups, most genes were assigned to the orders Bacteroidales (> 67,000), Clostridiales (> 40,000), Lactobacillales (27,000) and Enterobacteriales (> 14,000). The genus with the highest number of assigned genes was the genus Bacteroides (> 54,000). Figure 2 shows the relative abundance of transcripts per assigned genus.Figure 2Average relative transcript expression of different genus level taxa in incubations sampled at time points 0 h, 6 h, 12 h, 24 h and 48 h. When the taxonomic assignment could not be made at genus level, the lowest classifiable taxonomy assignment was used for display. Low abundance genera are summarized as “Other taxonomic assignments” for display purposes.To identify bacterial gene expression patterns, we focused on RNA reads for which KEGG Orthology (KO) or EC identifiers could be assigned. The percentage of reads with defined KO or EC ranged from 42 to 83% for different samples. Most of the data with assigned KO or EC identifiers came from 22 bacterial groups, of which 12 could be assigned to a known genus, and only a small number of genes was assigned to minor groups (3%), unclassifiable sequences (3%), and sequences not classifiable beyond the superkingdom Bacteria (3.5%). In the activated inoculum at the start of the incubation (t0), unclassified Enterobacteriaceae were the most active group (Fig. 2), probably due to residual oxygen during the activation. However, once the incubation had started, the relative expression of Bacteroides increased in all treatment groups. In all samples combined across all treatments and time points, 39% of all expression data came from the genus Bacteroides and 27% from unclassified Enterobacteriaceae. Overall, relative abundances corresponded with those based on 16S rRNA gene analysis as described by Gu et al.11 (Fig. 3).Figure 3Correlation between the relative activity of the main bacterial groups based on metatranscriptome data, and their relative abundance based on 16S rRNA gene sequencing data11. In case when genus level assignment was ambiguous, unclassified fraction within the next higher taxonomic level was used.Global and IMMP specific co-occurrence of taxaIt is known that in microbial ecosystems bacterial taxa can occupy different niches and co-exist forming a complex network of co-dependencies. We wanted to assess whether, based on the metatranscriptome data, we could identify bacterial groups which co-occurred in our samples, with particular attention to the effect of specific IMMPs. We performed clustering analysis based on mRNA reads from all samples in our dataset to test for global co-occurrence patterns. We showed that clustering into nine groups was most stable, as it could be reproduced in multiple rounds of clustering. These resulting nine different groups showed different behaviours over all investigated conditions. An overview of organism assignment per cluster, with number of assigned genes and differentially expressed genes is provided in Table S3. Cluster 1 (Table S3) was present in all t0 samples but decreased or was absent at all other time points. This cluster consisted mostly of reads assigned to Ruminococcus and Lactococcus as well as reads that could be largely classified as contamination from the inoculum/preparation (e.g. Homo, Mus, Bos, unclassified Mammalia). The second and third cluster consisted mainly of genera including many genera known to contain many probiotic organism candidates, i.e. Bifidobacterium, Lactobacillus, and Enterococcus, and sequences, which could not be classified beyond a related higher taxon (e.g. unclassified Bifidobacteriaceae, unclassified Lactobacillaceae). These clusters also contained a related phage group (Myoviridae, mainly Lactobacillus phages), and an unrelated genus (Fusobacterium). The identified genera in cluster two and three showed increasing relative transcript abundance in all cultures supplied with IMMP substrates, whereas relative transcript abundance was decreased or undetected in the control cultures. The fourth cluster was predominated by E. coli and related higher order classifications (e.g. unclassified Enterobacteriaceae), together with other enterobacteria such as Enterobacter, Citrobacter and Klebsiella, and the unrelated Gram-positive genus Eubacterium. This cluster was mainly present in the samples without IMMP and declined in the samples with IMMP. The fifth cluster was predominated by Bacteroides and showed an increase with time in all incubations. This cluster also included Parabacteroides, Prevotella, Flavobacterium, and Desulfosporosinus. The sixth cluster consisted only of Clostridium/unclassified Clostridia, which showed some increase with time in all incubations. The seventh cluster contained Anaerostipes and related higher level taxon classifications (unclassified Clostridiales, unclassified Lachnospiraceae) and showed a similar pattern as cluster six. No clear pattern was seen for the eighth cluster consisting of Corynebacterium, Ethanoligenes, Odoribacter, and Sutterella. Finally, the ninth group consisted of different bacterial genera, some related to non-carbohydrate metabolizing bacteria (Acidaminococcus, Bilophila, Phascolarctobacterium), and some known gut symbionts like Veillonella and Megasphaera. This group was common in samples of incubations without any prebiotics at 48 h and was nearly absent in all the other samples.Detection of specific gene expression patternsBesides the co-occurrence of bacterial groups, the specific gene expression patterns within these groups were investigated based on the optimal gene clustering for all bacterial groups using DBSCAN. The clustering with the optimal tau was chosen for all bacterial groups, except for the genus Enterococcus, for which a suboptimal tau led to better cluster separation. As a result, the DBSCAN gene clustering analysis revealed the presence of three main expression patterns observed in nearly all bacterial clusters. These three patterns comprised in all cases at least 80% of all investigated genes, which were not considered noise. The first pattern was present in all incubations and was characterized by genes which were expressed only at t0, and not expressed at any of the later time points. The second pattern was found only in the control group and only at 48 h. The third and most common pattern found in all experimental groups included genes that were not expressed at t0 but showed upregulation at the later time points during incubation. This pattern was characteristic for genes assigned to the genera Enterococcus and Bacteroides, for which respectively 40% and 99% of differentially expressed genes were increasingly expressed over time in all treatment groups including the control group. Bifidobacterium/Lactobacillus and Clostridium also showed the same pattern, but only in the groups where IMMPs were present. Anaerobutyricum hallii (formerly Eubacterium hallii) showed the same gene expression pattern, but only in the group supplemented with IMMP-27.Considering the experimental set-up chosen for this study, the expression levels of genes assigned to a specific bacterial group indicate the contribution of this group to the utilization of the specified substrate, or its by-products. The high overall relative expression of genes assigned to bifidobacteria (and unclassified Bifidobacteriaceae), lactobacilli, enterococci, and unclassified Actinobacteria was positively associated with the presence of IMMPs (Fig. 2). Conversely, the expression of genes assigned to unclassified Proteobacteria, Prevotella, Sutterella, Acinetobacter, Eggerthella, Acidaminococcus, Streptococcus, Phascolarctobacterium, and Bilophila was negatively associated with the presence of IMMPs, as compared to the control group.General metabolic effects of IMMPWe wanted to further investigate the expression of genes and corresponding bacterial groups associated with the fermentation of different IMMPs. Our analysis of the metabolic clusters revealed that five bacterial groups found in the faecal inoculum, namely Bifidobacterium/Lactobacillus, Enterococcus, Bacteroides, Clostridium, and Anaerobutyricum hallii, showed a considerable upregulation of general metabolic pathways like glycolysis, nucleic acid or fatty acid biosynthesis, as compared to the gene expression at t0. When we compared metabolic patterns between different bacterial groups, the groups exhibited overall different metabolic patterns. Members of the genus Bacteroides showed at first a unique partial upregulation of vitamin B12 metabolism. An investigation of the cofactor requirements showed that vitamin B12 in Bacteroides is essential for methionine synthase and methylmalonyl-CoA mutase, the latter of which produces methylmalonyl-CoA from succinyl-CoA (Fig. 4A).Figure 4Overview of the metabolism of specific microbial groups observed in the samples taken during in vitro fermentation of different IMMPs by human faecal inoculum All samples show in general the same patterns for all organisms, besides for Anaerobutyricum hallii, which only showed expression in the samples with IMMP-dig27. The genus Enterococcus showed the same pattern as Bifidobacterium/Lactobacillus, but at lower relative transcript abundance. Grey indicates that certain genes were not differentially expressed within a pathway. 5-ALA 5-Aminolevulinate, AC Acetate, Ac-CoA Acetyl-CoA, BUT Butyrate, FORM Formate, FUM Fumarate, GLC Glucose, LAC Lactate, PROP Propionate, PROP-CoA Propanoyl-CoA, PYR Pyruvate, SUC Succinate, SUC-CoA Succinyl-CoA.Methylmalonyl-CoA mutase is involved in propionate biosynthesis, and our data showed that the whole pathway for propionate biosynthesis was, in fact, upregulated. The data further showed that many genes coding for proteins involved in iron scavenging were also upregulated (e.g. FecR). One of the genes coding for an enzyme with iron requirements was that encoding succinate dehydrogenase, which converts succinate into fumarate. This function, as well as all others in the tricarboxylic acid (TCA) cycle, showed upregulation in all samples tested. The genus Clostridium also showed an upregulation of genes involved in Vitamin B12 production (Fig. 4B), but the biosynthesis occurred via glutamate, whereas in the Bacteroides group it was produced via succinate (Fig. 4A). The genes in the pathway for propionate production were overall upregulated (production via acetyl-CoA, not succinyl-CoA), similar to the genes in lactate and butyrate production pathways. The only other enzyme requiring vitamin B12 in the microbiome was a multimer of propanediol dehydratase or glycerol dehydratase (ambiguous functional assignment), both involved in the breakdown of glycerol/glycerone phosphate to propanol/propionate/1,3-propanediol. However, a full upregulation of either pathway was not observed. The Bifidobacterium/Lactobacillus group and the Enterococcus group showed upregulation of genes related to production of lactate from pyruvate (Fig. 4C), and the Clostridium group also showed upregulation of genes encoding proteins involved in butyrate production, but it is unclear if butyrate would be directly produced from pyruvate, or derived from external acetate (Fig. 4B). Anaerobutyricum hallii, on the other hand, showed high gene expression related to converting lactate into butyrate, as also shown previously12 (Fig. 4D). In addition, our data indicated that formate might have been produced by proteins encoded by genes belonging to the Enterococcus and Bacteroides populations (Fig. 4A).Microbial groups directly involved in the degradation of the IMMPsIn order to gain insight into which bacterial groups were directly involved in degradation of different IMMPs, we used the KEGG reference pathway for starch and sucrose metabolism13. We surveyed our data for the expression of the genes encoding enzymes that are known to be involved in sucrose and starch metabolism. More specifically we focused on genes encoding enzymes from glycoside hydrolase family 13 (http://www.cazy.org/GH13_bacteria.html), as this family includes a number of bacterial proteins shown to be essential in degradation of similar compounds, such as isomaltooligosaccharides (IMOs)14. The majority of genes listed in the KEGG starch and sucrose metabolism pathway were detected in our transcriptome data, as well as some additional genes in glycoside hydrolase family 13 (EC 3.2.1.135, 3.2.1.68 and 3.2.1.11), which were not listed in the KEGG pathway, but which have been shown to be activated during the degradation of pullulan and dextran15,16,17,18. It is interesting to note that the relative contribution of these starch and sucrose metabolism genes to the total number of genes from each sample did not correlate with the presence or absence of IMMPs in the samples. The only exception was incubation with IMMP-dig27, in which starch and sucrose metabolism genes reached an expression of 10% at 12 h and about 12% at 48 h, whereas in other groups they ranged between 4 and 5% (Fig. S2). Despite of the similarities in the overall expression of the starch and sucrose metabolism genes in all samples, we could see differences in the relative abundance of genes coding for specific enzymes depending on the IMMP used, and the duration of the fermentation (Fig. S3).One of the aims of this study was to better understand the functional dynamics of the bacterial communities during IMMP degradation. Previously reported HPAEC and HPSEC analyses11 showed that the degradation of IMMP-94 and IMMP-96 occurred between 12 and 24 h of the incubation. At 24 h and 48 h we noted an increase in the expression of genes coding for enzymes that might be directly involved in the hydrolysis of α − (1 → 6) glycosidic linkages, namely EC 3.2.1.10—oligo-1,6-glucosidase, EC 3.2.1.11—dextranase, and EC 3.2.1.33—amylo-α-1,6-glucosidase (Fig. S4A,B). There was also an increase in the expression of genes coding for enzymes that can hydrolyse α − (1 → 4) glycosidic linkages, mainly EC 3.2.1.1—α-amylase, EC 3.2.1.20—α-glucosidase 4-α-glucanotransferase, and EC 2.4.1.25—4-α-glucanotransferase. Since IMMP-27 contains lower amounts of α − (1 → 6) linkages, its degradation also involves the activation of the same genes, however, the expression levels of the genes encoding enzymes which hydrolyse α − (1 → 6) linkages were much lower (Fig. S4A,B). Bacterial groups that contributed the most expression to the primary degradation of IMMP’s α − (1 → 6) linkages were Lactobacillus, Bifidobacterium and Bacteroides, all expressing the genes encoding EC 3.2.1.10 oligo-1,6-glucosidase and EC 3.2.1.11 dextranase. On the other hand, the metatranscriptomic data suggested that α − (1 → 4) linkages were hydrolysed mainly by proteins encoded by Bacteroides, unclassified Bacteroidales, unclassified Enterobacteriaceae, Lactobacillus and Bifidobacterium via EC 3.2.1.1 alpha-amylase and EC 2.4.1.1 glycogen/amylophosphorylase (Fig. S5). Based on the transcript data, genes belonging to the genera Bifidobacterium and Lactobacillus were mainly expressed in the degradation of IMMP-94 and IMMP-96 at 24 h (Fig. 5, and Fig. S5). These genera showed also expression during degradation of IMMP-27 and the IMMP-dig27, but their relative contributions were much lower (Fig. 5, and Fig. S5). The breakdown of IMMPs at 24 h and 48 h was otherwise dominated by expression from Bacteroides, with the exception of IMMP-dig27 at 48 h, which showed a high level of expression of genes assigned to unclassified Enterobacteriaceae. Figure 5 summarizes our model of IMMP degradation and shows the proposed specialized role of lactobacilli and bifidobacteria in hydrolysis of α − (1 → 6) linkages. It also reveals the important contribution of Bacteroides as both primary and secondary degraders of IMMPs and their by-products.Figure 5Overview over the main degradation pathways starting from dextran. Colours in the top panel indicate the main contributors to a reaction. The bottom panel shows the overall expression in reads per kilobase per million (RPKM) per organism at time points 24 and 48 h for all conditions. If an enzyme could not be identified by its associated EC number, then the KO or CAZy identifier used for identification are given in brackets. Only enzymes, which were significantly differentially expressed over at least one time series, are displayed. Please note the different scales between the first and the last enzymes in the bottom panels, due to differences in expression scale.Comparison to other metatranscriptomic projectsWe identified few limitations in our study, such as the fact that our starting culture was a mix of faecal samples from several donors for which we lack the initial gut microbiota profiles. Another limitation was the use of in-vitro batch fermentation system which does not allow to control for nutrient depletion and by-products accumulation. Therefore, to validate our findings we compared our results with other in-vivo and in-vitro gut metatranscriptome projects. We selected all available Illumina single-end metatranscriptomic bioprojects at the NCBI, as well as two relevant bioprojects with paired-end Illumina data. For each project we profiled taxonomic composition and compared the projects using PCA on the inverted Bray–Curtis distances. As it can be seen in Fig. 6, our data showed stronger clustering than a dysbiotic microbiome project (PRJNA396964, from graft-vs-host disease patients19). Our project data was also closely positioned with clusters of PRJNA416988 project20, another in-vitro study, where samples were derived from a single person and RNA was extracted after 7 days in a bioreactor and 24 h in a Hungate tube. The second most closely positioned project was PRJNA509512, where samples were derived from healthy volunteers21. Another in-vitro project, PRJNA59378722, with also a roughly comparable setup, was displayed further away, indicating a variation in the in-vitro setups comparable to in-vivo data.Figure 6PCA based on Bray–Curtis distance between samples from different bioprojects. The current study (light blue) groups with samples from two other studies (green, purple), which are derived from either an in-vitro model based on a single volunteer or direct sequencing of healthy volunteer samples, respectively. The samples displayed in orange are derived from a graft-vs-host disease study and include disturbed microbiomes.
