Five phylogenetically different bacterial strains isolated from a soil-derived lignocellulolytic microbial consortium were selected. This selection was based on parameters established by Puentes-Téllez and Salles (2018), after screening of 60 synthetic five-species consortia varying in their functional metabolic diversity and drawn from a pool of 18 bacterial species. The selected strains (mixture 48) when grown together on PB (after 96 hrs of cultivation), could degrade up to 50% of the total lignocellulose fraction, with specific degradation percentages of 43.01, 39.65, 69.18% of lignin, cellulose and hemi-cellulose, respectively (Puentes-Téllez and Salles, 2018). Based on previous 16S rRNA gene sequencing, these strains were identified, at the genus level, as Stenotrophomonas, Paenibacillus, Microbacterium, Chryseobacterium, and Brevundimonas.

In order to analyze the temporal growth dynamics of the synthetic bacterial consortium, we cultivated the bacterial mixture (five strains) on four different substrates as sole carbon source varying in their complexity (glucose -GLC, carboxymethylcellulose -CMC, xylan -XYL and sugarcane bagasse -SCB). Reagents like GLC, CMC and XYL were obtained from Sigma–Aldrich. The SCB was previously milled into small pieces (≤1 mm), washed twice with 70% of ethanol and distilled water, and then dried at 55°C for 48 h. Firstly, each bacterial strain was cultivated in 5 ml of LB broth for 48 h. The cultures were centrifuged (6,500 g × 10 min) and the cells were washed with 1 ml of 0.9% NaCl. Absorbance (OD_{600 nm}) for all the strains was adjusted with 0.9% NaCl at 1.0 in order to ensure similar starting densities. Subsequently, 25 μl of each strain (relation 1:1) were inoculated into 25 ml of sterile mineral salt medium containing 0.5% of each substrate, vitamins and trace element solutions (Jiménez et al., 2014b). Triplicate flasks per time point were incubated at 28°C in aerobic conditions (180 rpm) and samples were taken at 12, 48, 96, and 192 h for subsequent analysis. Microbial growth was evaluated by absorbance (OD_{600 nm}) for each time point. Absorbance in the SCB-derived cultures was measured directly at the upper layer, after substrate sedimentation. Two controls, i.e., one without substrate and the other without microbial source were set up (Supplementary Figure S1).

Sugar and protein concentrations, as well as the enzymatic activities of the consortium cultivated on XYL and SCB were measured in the extracellular fraction. To obtain these fractions, 10 ml of each culture at each time point were centrifuged (6,500 g × 8 min). The supernatant was clarified by passing through a 0.22 μm filter. Quantity of sugars and proteins was obtained by 3,5-dinitrosalicilic acid (DNS) method (Miller, 1959) and the Quick Start^TM Bradford Protein Assay (Biorad, Hercules, CA, United States), respectively. Enzymatic activities related with PB degradation were evaluated using pNP-labeled substrates, i.e., p-nitrophenyl α-D-glucopyranoside (pNPαGlu), p-nitrophenyl β-D-xylopyranoside (pNPβXyl), p-nitrophenyl β-D-galactopyranoside (pNPβGal) and p-nitrophenyl α-D-mannopyranoside (pNPαMan). The reaction mixtures consisted of 200 μl of 7.5 mM of each p-nitrophenol-glycoside (diluted in 50 mM of Tris-HCL pH 7.5) and 150 μl of each supernatant. The mixtures were incubated at 37°C for 4 h, after which the reactions were stopped on ice. Two negative controls were used for all assays: (i) reaction mixture without pNP-substrate; (ii) reaction mixture using sterile water. Enzymatic activities were determined from the measured absorbance units using a standard calibration curve. Activity was detected by the presence of yellow color and the amount of para-nitrophenol liberated was measured by absorbance at 410 nm (Maruthamuthu et al., 2017). All the analyses were performed using three biological replicates and standard deviations (σ) were obtained.

Total microbial RNA extraction from the consortium was carried out along the batch of culture on XYL and SCB. For each time point (12, 48, 96, and 192 h) three biological replicates were taken. Briefly, 10 ml of RNAprotect Bacteria Reagent (Qiagen, Hilden, Germany) were added to 5 ml of microbial culture. After 5 min (at room temperature) the mixture was centrifuged at 5,000 g for 10 min. The cell pellets were lysed using proteinase K (50 mg/ml) and lysozyme (15 mg/ml), and the RNA extraction was performed using the RNeasy Mini Kit (Qiagen, Hilden, Germany) following the manufacturer’s instructions. The quality and quantity of the extracted total RNA was checked by non-denaturing agarose gel electrophoresis and absorbance ratios 260/280 nm. Depletion of ribosomal RNA (using Ribo-Zero^TM Kit), mRNA libraries preparation and sequencing (Illumina NextSeq 2 × 150 bp; 1 lane) were performed at LGC Genomics (Berlin, Germany).

A total of 116,205,704 quality filtered paired-end reads (>20 bp) from 21 samples were used as input for de novo transcript assembly using Trinity (Haas et al., 2013). Transcript abundance estimation and normalization was performed using scripts included in the Trinity package. The script abundance_estimates_to_matrix.pl was used to create a matrix with the estimated RNA-Seq fragment raw counts that were used for differential expression analysis with edgeR (Robinson et al., 2010). To enable direct in-between sample comparison of transcript abundances, raw abundances were converted to TPM (transcripts per million transcripts) and TMM (trimmed mean of M-values) normalized in RSEM (RNA-Seq by Expectation Maximization) (Dillies et al., 2012). The script analyze_diff_expr.pl was used to extract those transcripts that were at least fourfold differentially expressed at a FDR (false discovery rate) corrected P-value cut off of 0.001 in any of the pairwise sample comparisons, followed by hierarchical clustering based on the Pearson correlation matrix or pairwise comparisons of samples using edgeR (Robinson et al., 2010). Moreover, total transcript contigs were converted to potential coding regions using the TransDecoder function developed within Trinity. For functional classification, the proteins were annotated using the GhostKOALA (Kanehisa et al., 2016) and dbCAN web-platforms (Yin et al., 2012). In order to unveil the relative expression dynamics of lignocellulolytic proteins, and in accordance with previous omics studies (Simmons et al., 2014; Zhou et al., 2014; Berlemont and Martiny, 2015; Jiménez et al., 2015a; 2016; Abot et al., 2016; Antunes et al., 2016; Zhu et al., 2016), we selected thirty carbohydrate-active enzyme (CAZy) families (Lombard et al., 2014) that are highly relevant in the deconstruction of PB. All transcripts affiliated within those CAZy families were manually annotated by BLASTp against the NCBI-nr database. The taxonomic affiliation of each transcript was performed by manually taking the most frequent BLASTp hit between the top-50. Finally, the TMM-normalized TPM values were added to each transcript classification and summed up to obtain the values per each CAZy family (Pfreundt et al., 2016). Relative expression of the CAZy families that segregated significantly between the time points were identified using random forest analysis with 1000 trees followed by the Boruta algorithm for feature selection (average z-scores of 1000 runs > 3) (Cardenas et al., 2015). All metatranscriptomic data are publically accessible on the MG-RAST (Meyer et al., 2008) server under the IDs mgm4734165.3 to mgm4734173.3; and mgm4735448.3 to mgm4735459.3.

The growth of the consortium was evaluated on four carbon sources with different structural complexity. The results showed that the consortium is able to grow using GLC, XYL, and SCB, but not CMC as a sole source of carbon (Figure 1A). In the GLC culture, a maximal growth was observed at 48 h, gradually decreasing toward the end of the batch (192 h). Regarding the XYL culture, we observed the exponential growth phase after 12 h, followed by the start of the stationary phase at 96 h. In the case of SCB, the exponential phase started at an early growth stage (between 0 to 12 h), and apparently was maintained until 96 h of the culture. Subsequent analyses were performed on the cultures with bacterial growth in substrates with higher complexity (XYL and SCB cultures).

In the XYL and SCB cultures, the concentration of the total proteins secreted by the consortium was less than 0.1 mg/ml. Regarding release of reducing sugars, maximum values (0.30–0.34 mg/ml) for the XYL culture were found between 48 and 192 h, during the exponential and stationary phases. In the SCB, we observed an increase of sugars during the entire growth period, with a maximum value at the end of the culture (0.28 ± 0.01 mg/ml; mean ± SD) (Figure 1B). Along the incubation time, we evaluated four types of enzymatic activities (target enzymes that can act on terminal ends of substrates) related with hemicellulose (β-xylosidase, β-galactosidases, α-mannosidase) and pectin (α-glucosidase) degradation processes. Thus, the consortium cultivated on XYL displayed maximal β-xylosidase and α-glucosidase activity values at 48 and 192 h, which was consistent with the sugar accumulation data (Figure 1C). Concerning the consortium cultivated on SCB, we observed a high β-xylosidase, α-glucosidase and β-galactosidases activity at 48 h versus a low α-mannosidase activity. The activities were relatively stable until the last measurement, except for β-xylosidase where a significant peak (p < 0,001, t-test of pairwise comparisons) was observed at 48 h (Figure 1D).

Total RNA extraction from the consortium cultivated on XYL and SCB was performed on triplicate systems along four time points of incubation. However, it was not possible to isolate sufficient RNA at the initial stage (12 h) of the XYL culture due to the low microbial activity and the presence of undigested xylan. Thus, a total of 21 samples were subjected to mRNA sequencing (9 from XYL and 12 from SCB cultures). Approximately 5.7 and 7.1 Gb of total mRNA sequencing data were obtained for the XYL and SCB samples, respectively. The approximately 116 million reads were assembled into 59,734 contigs. Given that N50 values discard read coverage information, we report the E90N50 value (1.1 kb) as a more useful indicator of transcriptome assembly quality than the N50 (665 bp) (Haas et al., 2013). Also note that the E90 number of transcripts (23,364) for which the E90N50 value is computed is roughly a third of the total number of transcripts assembled (68,738) and for which the N50 statistic was calculated. Based on the RNA-Seq fragment raw counts, we observed that 4,213 transcripts were differentially expressed in all pairwise combinations. As expected, the clustering of these transcripts showed that the XYL and SCB cultures had different expression profiles. For each, the profiles clustered by triplicates divided by treatment. Only the SCB time points 96 and 192 h formed an intermingled cluster based on the differential expression profile (Supplementary Figure S2), suggesting that after 96 h of growth on SCB the consortial expression patterns were ‘stable.’

Thirty CAZy families were selected in order to evaluate the expression of transcripts related with PB-degrading enzymes in the consortium. Here, we focus our analyses on the SCB cultures. A heat map using TMM-normalized TPM values showed a trend in the expression profiles of enzymes from the consortium along the incubation time (Figure 2 and Supplementary Table S1). In general terms, most of the selected CAZy families were highly expressed at 12 and 48 h compared with 96–192 h. For instance, xyloglucanases/endo-glucanases (GH16), endo-xylanases (families GH8 and GH10) and α-glucuronidases (GH67) were significantly expressed at 12 hrs. Moreover, lytic polysaccharide monooxygenases -LPMOs (AA10), exo-glucanases (GH6 and GH48), α-arabinosidase/β-xylosidase (GH43), next to other endo-glucanases (GH5 and GH9) and endo-xylanases (GH11) were highly expressed at 48 h compared with all other time points. In addition, enzyme families involved in pectin degradation (e.g., GH28 and PL10) and β-glucosidases (specifically GH3) were mostly expressed at 48 hrs. Regarding lignin degradation, the family AA2 (catalases-peroxidases) showed higher values of expression at 12 (117 ± 42) and 192 hrs (118 ± 26). In terms of relative expression values (i.e., TMM-normalized TPM) between the CAZy families in each time point, we observed that transcripts of enzymes of families GH16 (1115 ± 447) and GH10 (1005 ± 207) were highly expressed at 12 h (Figure 3A). Moreover, those from families GH11 (439 ± 42), GH10 (339 ± 2.5), GH43 (414 ± 44), GH3 (236 ± 31) and GH9 (163 ± 20) were the most expressed at 48 hrs. At the final stages of the cultures (96–192 h), the relative expression values of enzymes belonging to families GH43, GH13, GH11, GH3 and CE3 were higher in comparison with those of the other families.

The relative expression of the rpoA gene (housekeeping single-copy gene throughout bacteria) was used as a marker to evaluate the metabolic activity of each member the consortium along the incubation time on SCB, thus providing a depiction of the relative contribution of each strain in the CAZy expression profile described above. The taxonomic affiliation and the TMM-normalized TPM values of the rpoA transcripts showed that Brevundimonas and Paenibacillus were the most active members at 12 and 48 h, respectively. In contrast, Chryseobacterium and Stenotrophomonas were highly active at 96–192 h compared with the other three consortium members. Microbacterium showed low activity along the entire incubation time (Figure 3B). In addition, we characterized all transcripts encoding enzymes of the thirty CAZy families to species (Figure 4). This analysis indicated that Brevundimonas and Paenibacillus were major contributors to the expression of enzymes involved in cellulose degradation. Specifically, Brevundimonas was relevant in the expression of β-glucosidases/galactosidases (GH1 and GH3) at 12 and 48 h (transcript IDs 14 and 15). The expression of enzymes of CAZy families GH10, GH11, GH30, GH43, GH51, and GH67, Paenibacillus and Chryseobacterium were the most important drivers of the degradation of xylan/arabinoxylan along incubation time (transcript IDs 17–26). In addition, endo-xylanases (GH8), endo-xyloglucanases (GH74), β-galactosidases (GH2), xyloglucanases/mannosidases (GH44) and β-endomannases (GH26) were mainly expressed by Paenibacillus at 48 h. Pectin degradation was mainly associated with the expression of enzymes of families GH28 and PL10 from Brevundimonas at 48 h (transcript IDs 47 and 50). The expression of catalases-peroxidases (AA2) was carried out by Paenibacillus (at 12 h), Chryseobacterium and Brevundimonas (at 192 h). Taxonomic affiliation, TMM-normalized TPM values (average and standard deviations) and the functional annotation of the 50 most important transcripts within the 30 CAZy families can be found in the Supplementary Table S2.

Out of the thirty CAZy families evaluated, a total of 23 were highly expressed at this stage (row Z-score > 0 in at least one biological replicate). By using a random forest analysis and the Boruta statistical test, eleven CAZy families were significantly more expressed at 12 h (Figure 2), suggesting that this is the most functionally diverse phase compared with the other time points. In this stage, transcripts encoding enzymes of CAZy families GH8, GH10 and GH16 were significantly expressed compared with other time points. The family GH16 contains enzymes that can be assigned to five subgroups according to their substrate specificities, including xyloglucan transglucosylases/hydrolases (XTHs), (1,3)-β-galactanases, (1,4)-β-galactanases/κ-carrageenases, “non-specific” (1,3/1,3; 1,4)-β-D-glucan endohydrolases and (1,3; 1,4)-β-D-glucan endohydrolases (Strohmeier et al., 2004). Specifically, in our study, the high expression of this family at 12 h was due to one transcript that was affiliated to a 1,3-1,4-endoglucanase from Paenibacillus (78% identity) (Supplementary Table S2). The relevance of GH16 proteins has been acknowledged in previous studies. For instance, proteins from this family have been found in high relative abundances in the metagenomes of three aerobic PB-degrading microbial consortia (Jiménez et al., 2016). In addition, Zhou et al. (2014) reported that one of the microbial members (Clostridium -FC2) of the cellulolytic consortium F1RT has the capacity to produce cellulosomal enzymes belonging to family GH16. In Trichoderma reesei and T. harzianum these types of enzymes have been found highly expressed (or induced) in presence of cellulosic substrates (e.g., cellulose, avicel and sophorose) (Häkkinen et al., 2012; Ferreira Filho et al., 2017). On the other hand, proteins from GH16 have a putative structural and evolutionary connection with endoxylanases from family GH11 (Strohmeier et al., 2004). Based on the lack of experimental support on the cellulolytic activity of GH16 transcripts, we can speculate that, in our system, these enzymes may have an important role in transformation of hemicellulose (e.g., xyloglucan) present in the SCB. Based on the ability of our consortium to grow on GLC, XYL and SCB, but not on CMC, we suggest that the expression of enzymes involved in cellulose degradation (e.g., endo-glucanases) could be regulated by the presence natural cellulose of other type of PB polysaccharides present in the SCB. CMC could be used as an inductor of endoglucanases (Xiao et al., 2005). However, this polymer has a different structure compared to the natural cellulose found in PB.

In our study, two transcripts (IDs 17–18) contributed to a large extent to the high expression level of enzymes from family GH10 at 12 h (Figure 4). These were affiliated to endo-1,4-β-xylanases from Paenibacillus (>96% identity). Proteins from the family GH10 and GH16 were also strongly expressed and secreted in a corn stover-adapted microbial consortium at 24 and 72 h, respectively (Zhu et al., 2016). The high expression level of endo-xylanases in the consortium at 12 h suggests that xylan/arabinoxylan is attacked at early growth stages. One of the most important features of this stage was the high metabolic activity of Brevundimonas (Figure 3B). This organism may participate in the expression of endo-glucanases and endo-xyloglucanases from families GH5 and GH74, respectively. It is key to note that proteins from family GH5 can breakdown different plant polymers [e.g., cellulose, xyloglucan, and galacto(gluco)mannan]. Proteins derived from organisms affiliated with Brevundimonas species have been found in high relative abundance in two consortia cultivated with once-used PB (Jiménez et al., 2016). These results suggest that Brevundimonas may play an important role in xylan, cellulose and lignin degradation. Moreover, the expression of α-glucuronidases – which cleave 4-O-methyl-α-D-glucuronic acid (MeGlcA) from xylose residues when they are attached to the terminal, non-reducing end of xylo-oligosaccharides (Dodd and Cann, 2009) – by Chryseobacterium and Paenibacillus was significantly higher at 12 h, indicating that the consortium can carry out this process at the early stages of cultivation. Interestingly, the removal of MeGlcA could expose the xylan to attack by endo-xylanases and β-xylosidases (Rogowski et al., 2014). Two transcripts (ID 41 and 44) affiliated to AA2 and AA6 (and derived from Paenibacillus and Brevundimonas, respectively) were comparatively highly expressed at 12 h (Figure 4). We speculate that the expression of these transcripts is probably related to cell protection from oxidative damage and the transfer of electrons.

At this stage, the low cell densities and the presence of accessible PB polysaccharides could enhance the fast growth of Paenibacillus and Brevundimonas. Here, competition and niche partitioning events are likely to occur (Jiménez et al., 2017). Thus, based on the high expression of Brevundimonas-derived transcripts affiliated to α-fucosidases (GH29) and β-galactosidases (GH35) at this stage (Figure 4), we suggest that this species may consume preferably fucose and galactose, whereas Paenibacillus could have more access to arabinose and mannose. Competition for the most abundant sugars (glucose and xylose) could be a valid explanation for the lower activity of Brevundimonas in the next stages. At an early phase, Microbacterium could consume easily available oligo- or monosaccharides present in the SCB. However, we suggest that the low metabolic activity of Microbacterium could be directly related with an intrinsic characteristic of the organism (e.g., oligotrophy), niche competition and possibly a low metabolic versatility. Probably, in the consortium, Microbacterium could confer functional stability, redundancy and robustness, thus increasing the metabolic diversity and providing a potential buffer against external disturbance (Awasthi et al., 2014). In addition, Microbacterium species could secrete molecules that synergistically stimulate the degradation capacity of other PB-degrading bacterial species.

This is the key stage for the degradation of the main PB polysaccharides. Indeed, we observed raised extracellular enzymatic activities at 48 h, with a significant peak of β-xylosidase (Figure 1D). Here, efficient hydrolysis of cellulose can be carried out by a synergistic action of non-processive endo-glucanases (GH5 and GH9), which produce new ends at random within the polysaccharide chain, and processive exo-glucanases (GH6 and GH48), which remain attached to the substrate and split off cellobiose from such free ends. Then, cellobiose can be converted to glucose by the action of β-glucosidases expressed by Brevundimonas. An interesting synergistic interaction is carried out at this stage, where consumption of cellobiose by Brevundimonas potentially prevents metabolite repression of endo/exo-glucanases produced by Paenibacillus (Lynd et al., 2002; Berlemont and Martiny, 2013). Moreover, LPMOs from Paenibacillus may play a relevant role in the cellulose degradation at this stage. These LPMOs are metalloenzymes (copper-dependent) that catalyze the oxidative cleavage of (1,4)-linked glycosidic bonds of crystalline PB polysaccharide surfaces (Horn et al., 2012; Agger et al., 2014).

Concerning hemicellulose degradation, we found that transcripts affiliated to Paenibacillus from the families GH11 (endo-xylanases) and GH43 (α-arabinosidases/β-xylosidases) were highly expressed when compared to the other time points. We found a shift in the expression of endo-xylanases – i.e., from GH10 at 12 h to GH11 at 48 h – which display an endo-action and thus hydrolyze β-(1,4)-linkages between adjacent xylose residues within the xylan/arabinoxylan chain. The enzymes of these families essentially differ in their structure and substrate specificity. For instance, GH10 xylanases have a (β/α)₈-fold and are much more versatile. In contrast, GH11 enzymes are the smallest xylanases with a β-jelly-roll architecture, are highly substrate-specific and do not tolerate high substitutions on the backbone (Biely et al., 1997; Paës et al., 2012). This suggests that a high level of enzymatic complementarity can occur in Paenibacillus within the consortium along the growth period. Regarding GH43, it can display dual action (Maruthamuthu et al., 2017), thus removing arabinose from the external-side chains attached to the xylan backbone or hydrolyze xylobiose into its monomeric units, releasing D-xylose (Rytioja et al., 2014). The family GH43 enzymes seem to play a key role in lignocellulose degradation, being highly enriched in two soil-derived microbial consortia cultivated on wheat straw (Jiménez et al., 2015a) and highly expressed and secreted between 72 to 168 h in PB-degrading microbial consortia cultivated under aerobic and anaerobic conditions (Abot et al., 2016; Zhu et al., 2016).

Although Paenibacillus was the most metabolically active organism and possibly the biggest contributor to the (hemi)cellulose deconstruction, Brevundimonas may play a key complementary role through the expression of pectinolytic enzymes (e.g., GH28 and PL10). Pectins comprise a diverse and complex group of polymers that function as structural components in PB (Mohnen, 2008). Thus, the action of pectinolityc enzymes could increase the sugar release from the (hemi)cellulosic fraction (Xiao and Anderson, 2013). Theoretically, different consortium members can contribute in distinctive ways and proportions to the SCB degradation process. However, we suggest that organisms that contribute more to the substrate transformation may profit from the most abundant sugars (glucose or xylose) (Jiménez et al., 2017). Based on this conceptualization, at this stage, Paenibacillus might indeed be the most abundant organism in terms of cell density. In other words, Paenibacillus invests a lot in the production of lignocellulolytic enzymes, and it might grow very well at the expense of the large proportion of sugars.

At this final stage, the expression of the majority of the CAZy families decreased. A similar trend occurs with the metabolic activity of Paenibacillus and Brevundimonas. These two organisms may be acting as biggest contributors to the degradation process at early stages (12 and 48 h) of the culture, At this stage, Paenibacillus probably entered stationary phase, which is reflected in the decrease of GHs gene expression. We suggest that at these time points the major fraction of the available (hemi)cellulose is already degraded. Here, not only carbon sources but also other essential nutrients, such as nitrogen, can become progressively limiting. Although production of sugars (by PB polysaccharide degradation) and their consumption occurs simultaneously in the consortium, we observed an increase in the concentration of reducing sugars along the time of incubation (Figure 1B). This suggests that, at the end of the culture, the environment becomes relatively more homogeneous. Thus, the increasing release of easily assimilable carbon sources is expected to enhance competition among the bacterial strains. Stenotrophomonas and Chryseobacterium showed higher normalized expression values of their rpoA gene at the end of the culture (192 h) (Figure 3B), suggesting that they become highly metabolically active when the concentration of monosaccharides is higher. We considered Stenotrophomonas to be a ‘sugar cheater’ in the consortium. That means that it profits from the sugars made available by the other community members, without effectively contributing to the PB degradation process (Jiménez et al., 2017). This is the case of Thermus thermophilus, which lacks the metabolic potential for PB polysaccharide hydrolysis, in a switchgrass-degrading consortium (D’haeseleer et al., 2013). Moreover, the exact functional roles of Microbacterium and Stenotrophomonas in the consortia remained unknown and further studies need to be performed in order to understand their impact in this system. However, we speculate that they could be involved in feedback regulation, cross-feeding, maintenance of functional stability and inhibition of metabolite repression. Finally, at this stage, lignin bioconversion was thought to take place by the expression of catalases-peroxidases (CAZy family AA2) from Chryseobacterium and Brevundimonas. Similar to patterns observed in composting (Antunes et al., 2016), where high expression of ligninolytic enzymes occurs at the end of the process.