Comparative metagenomic discovery of the dynamic cellulose-degrading process from a synergistic cellulolytic microbiota

To reveal the dynamic process of cellulose biodegradation and explore more potential cellulases, a microbiota (FPDM) with cellulose-degrading ability was cultivated, and different stages of filter paper degradation were compared. Ion chromatography and comparative metagenomic sequencing revealed that the diversity of FPDM enhanced as the hydrolysate length diversity increased. Sporocytophaga and Cohnella dynamically dominated the synergistic degradation of cellulose in early-intermediate and intermediate-final periods, respectively. Moreover, 363 declining shifting and 231 progressive shifting unannotated genes were speculated to participate in the catabolism of cellulose to cellodextrin/cello-oligosaccharide and to cellobiose, respectively. Based on the dynamic changes in hydrolysates, community structure and gene abundance, a dynamic cellulose-degrading pathway of FPDM was predicted. Our work should provide a new perspective for subsequent identification of key cellulolytic strains and enzymes and clarification of the mechanism of cellulose biodegradation.


Introduction
Currently, the development of cellulose-derived renewable biofuels is of great significance for recycling agricultural waste and relieving the shortage of fossil fuels (Rodionova et al. 2017). However, the biodegradation of cellulose is difficult due to its recalcitrant nature and the relatively low activity of currently available hydrolytic systems, which is considered to be a major obstacle to industrial-scale production of fuel from cellulose (Jayasekara and Ratnayake 2019). Therefore, it is urgent to further clarify the molecular mechanism of cellulose biodegradation and explore efficient cellulose hydrolysis systems.
Recently studies revealed that cellulose decomposition by pure cultures of the cellulolytic strain was Ming Yang and Jingjing Zhao have equally contributed to this work. less efficient than in mixed cultures (Dumova and Kruglov 2009;Jiménez et al. 2014;Wongwilaiwalin et al. 2010). In natural environments, cellulolytic microbes work together with other microorganisms that utilize soluble sugars (cellobiose and glucose) to catabolize insoluble cellulosic matter (Deng and Wang 2016;Dumova and Kruglov 2009). The cellulolytic microbes could secrete several cellulases, including endoglucanase, cellobiohydrolase and beta-glucosidase (Juturu and Wu 2014), which catalyze the hydrolysis of the beta-1,4-glycosidic bonds in a cellulose molecule (Herrera et al. 2019;Lee et al. 2016). Catabolite repression of cellulases can be relieved due to the subsequent consumption of cellulose degradation products (Dumova and Kruglov 2009). Therefore, the de novo synthesis of cellulases and increased biomass of cellulolytic microbes drive continuous decomposition of cellulose. However, due to the diversity of natural cellulose substrates and habitats, these cellulolytic communities are biologically diverse, which hampers the elucidation of their structure and function (Wilson 2011).
As a result of the rapid development of next generation sequencing technology in recent years, metagenomic sequencing has become a powerful strategy for exploring the structure of cellulolytic microbiota and mining functional cellulose degradation related strains and enzymes (Knight et al. 2018). For instance, metagenome analysis of a rice strawadapted microbial consortium enriched from compost revealed a distinctive set of microorganisms and carbohydrate-active enzyme (CAZyme) genes distributed within these bacteria phyla, suggesting a synergistic system efficient in catabolism of cellulosic components in the compost habitat . With the help of culture-independent ''omics'' techniques, the metabolic ''blue-print'' of a proficient cellulose-degrading microbiota enriched from soil sample was predicted, and the putative operation of cellulosome and free-enzyme within the community was demonstrated (Zhou et al. 2014). Despite the exciting findings from metagenomic researches of cellulolytic microbiota, there is still much we do not understand about the cellulose-degrading mechanism and the function of numerous putative proteins, which can be partly attribute to the lack of knowledge about the dynamic changes of the microbial community at different cellulose-degrading stages (Mande et al. 2012).
To further clarify the cellulose degrading mechanism and explore more potential cellulases, we cultured a stable soil-derived microbiota (FPDM) capable of efficiently degrading filter paper in the present study. The dynamic changes in filter paper hydrolysates and the microbiota composition at different degradation stages were respectively analyzed by ion chromatography (IC) and metagenomic shotgun sequencing, revealing the synergistic action of the core functional strains. Notably, based on shorttime series expression miner (STEM) algorithm (Ernst and Bar-Joseph 2006), the function of some key functionally unannotated genes (FUGs) were predicted and a more distinct metabolic pathway involved in cellulose degradation were mapped.

Materials and methods
Enrichment of the filter paper-adapted microbiota from soil samples Lignocellulose-rich soil samples were collected from woodland soil in the campus of Dalian Polytechnic University, China (38°58 0 59.2176'' N, 121°31 0 55.1964'' E). These soil samples were respectively suspended in sterilized water and dispersed onto a double-layered filter paper placed on a solid medium, then pressed with a flat-ended glass rod to ensure sufficient contact between the soil and filter paper and incubated at 30°C. After 15-day fermentation, a piece of filter paper containing the cellulosedegrading microbiota was cut and transferred onto the same filter paper medium and incubated at 30°C for 15 days. Then the second batch microbiota culture was transferred again onto the same new medium. This subculture process was repeated for 10 times, and a stable microbiota named FPDM with an efficient cellulolytic ability was established. Further fermentation was conducted in double-layered filter paper solid medium and liquid filter paper medium, respectively, and performed in biological triplicate. The doublelayered filter paper solid medium applied a circular Whatman No. 1 filter paper (Whatman, Tokyo) as the only carbon source. The filter paper was presoaked with sterilized water and then placed on the agar medium with the inorganic components. The inorganic components of the solid media contained per liter: 0.5 g each of KCl, NaNO 3 , and MgSO 4 7H 2 O; 0.01 g of FeSO 4 7H 2 O; and 1 g of K 2 HPO 4 , pH 7.0. To prepare the liquid medium, 5 g per liter filter paper (cut into pieces 5 mm in length and width) was added into the same inorganic solution as the solid medium. The supernatant was collected from 100 mL microbiota FPDM culture after being cultivated in liquid medium at 30°C with shaking (180 rev/min) for 3, 8, and 15 days for the subsequent IC analysis. The FPDM culture after being cultivated in solid medium for 3, 8, and 15 days was also sampled for the subsequent metagenomic analysis.

Ion chromatography
Considering the cellulose hydrolysates were difficult to be detected from the FPDM solid culture, hydrolysates from FPDM liquid culture were determined by high-performance anion-exchange chromatography using an ICS-5000 IC system (Dionex, USA) equipped with a pulsed amperometric detector and a CarboPac PA100 column. The samples were analyzed by gradient elution. A 1.0 mL/min of mobile phase flow rate with a 25 lL of injection volume was adjusted at 30°C of column temperature.

DNA preparation and metagenomic sequencing
Totally 1 g of cellulolytic microbiota FPDM samples adhere to the filter paper were collected for nucleic acid extraction. Metagenomic DNA was isolated using a PowerSoil Kit (MoBio, USA) following the manufacturer's recommendations. Nanodrop-1000 (Thermo Scientific, USA) and gel electrophoresis were used to determine the DNA quality. High-quality DNA was then stored at -20°C before library construction.
DNA libraries were constructed using a NEB-NextÒ Ultra TM DNA Library Prep Kit for Illumina (New England Biolabs, USA) under the manufacturer's protocols. First, 1 lg of each metagenomic DNA sample was fragmented to a size of * 350 base pairs (bp) using an ultrasonicator (Covaris, USA). Next, the fragmented DNA was end-polished, A-tailed, and ligated with Illumina adaptors for further polymerase chain reaction (PCR) amplification. The PCR products were then purified using the AMPure XP system (Beckman Coulter, USA); libraries were profiled for size distribution using an Agilent 2100 Bioanalyzer (Agilent Biotechnologies, USA) and quantified by quantitative PCR . Finally, the high-quality libraries were sequenced on an Illumina NovaSeq 6000 platform (Novogene, China) using the paired-end 150 bp sequencing strategy.
Metagenomic data preprocessing, assembly, and annotation Raw reads that contained a [ 15 bp of adapter overlap, [ 40 low-quality bases (Q-value \ 38) or [ 10% ''N'' bases were removed using SOAPnuke (Chen et al. 2018). The remaining clean reads from each sample were assembled independently into scaftigs without ''N'' bases, using SOAPdenovo2 with a k-mer length of 55 bp (Luo et al. 2012;Qin et al. 2014). The clean reads from each sample were then mapped against the corresponding scaftigs using SOAPaligner (Li et al. 2009). Subsequently, all unmapped reads were mixed and assembled using SOAPdenovo2 with the same parameters. Open reading frames (ORFs) were predicted from scaftigs longer than 500 bp using MetaGeneMark (Sunagawa et al. 2015;Zhu et al. 2010). To obtain a nonredundant gene catalog, ORFs of more than 100 bp were clustered at 95% sequence identity using CD-HIT (Li and Godzik 2006;Sunagawa et al. 2015). SOAPaligner was used to align the clean reads with the nonredundant gene catalog, and only genes supported by at least three mapped reads were designated as unigenes for further analysis (Qin et al. 2010). The unigene abundance was calculated for each sample as the number of reads mapped to the unigene and normalized by unigene length (Villar et al. 2015). The unigenes were mapped against the National Center for Biotechnology Information (NCBI) nonredundant protein (NR) database (ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/) and the available fungi genomes from EnsemblFungi database (https://fungi.ensembl.org) for taxonomic annotation using blastp in DIAMOND (Buchfink et al. 2015), with a maximum e-value of 1e-5. MEGAN (Huson et al. 2011) was applied to determine the taxonomic level of each unigene using the lowest common ancestor-based algorithm. The abundance of each taxonomic group was estimated by summing the abundance of unigenes assigned to the same feature . The relative abundance of sequence from eukaryotes, virus and unclassified organisms were also calculated, and showed that the sum of these relative abundances is only 0.052%. Therefore, downstream analyses focused only on prokaryotic taxa, genes and genomes. The Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Kanehisa et al. 2008) was used to evaluate the function and metabolic pathways of the unigenes via blastp in DIAMOND. To mine the carbohydrate-active genes, the unigenes were annotated and selected using a hidden Markov model-based database for automated CAZyme annotation (dbCAN) (Lombard et al. 2014;Yin et al. 2012). To analyze the shift pattern of relative abundance, the STEM algorithm was employed to cluster the relative abundance variant profiles of the genus-level taxa, functional genes and the CAZyme families respectively.

Metagenomic contigs binning
In order to recover high quality metagenome-assembled genomes (MAGs), all clean reads were mixed and assembled by SOAPdenovo2, and the resulting contigs were then binned using three different binning software tools with Binning module in MetaWRAP pipeline (v1.2.1) (Uritskiy et al. 2018), including MetaBAT2 (Kang et al. 2019), MaxBin 2.0 (Wu et al. 2016), and CONCOCT (Alneberg et al. 2014). The resulting bins were collectively processed to produce consolidated MAGs with both Bin_refinement (criterion: completeness [ 70%; contamination \ 10%) and Reassemble_bins modules in MetaWRAP. The bins resulting from the Bin_refinement module were qualified with Salmon using the Quant_bins module of MetaWRAP. Bins obtained from each relative step above were subjected to quality control by estimating completeness and contamination using CheckM (Parks et al. 2015). Gene prediction and functional annotation were carried out using Annotate_bins in MetaWRAP based on Prokka gene database (Seemann 2014). CAZyme annotation for each bin was performed using the same method to that for unigenes as described above. These MAGs were then uploaded to MiGA (Rodriguez-R et al. 2018) for a complementary analysis to determine the most likely taxonomic classification and novelty rank of the bins at genuslevel resolution.

Statistical analysis
The alpha diversity was characterized by calculating the Shannon-Wiener and Simpson indices using Mothur (Schloss et al. 2009). The R package ''agricolae'' (http://tarwi.lamolina.edu.pe/*fmendiburu) was used to determine significant differences (P \ 0.05) in alpha diversity indices and their change rates between treatments. The beta-diversity for Bray-Curtis dissimilarities was calculated using the ''beta_diversity.py'' program in QIIME (Caporaso et al. 2010), and visualized by nonmetric multidimensional scaling (NMDS), which was performed using the ''WGCNA'', ''stat'', and ''ggplot2'' packages in the R software. Bar plots and heat maps were created using the R packages ''ggplot2'' and ''pheatmap'', respectively. Pairwise comparisons of the taxonomic and functional profiles between treatments were carried out using Metastats (White et al. 2009).

Prediction of gene function and construction of cellulolytic metabolic pathway
Based on the abundance variation patterns, FUGs were clustered together with annotated cellulolytic genes (ACGs) and were thus assigned a certain function similar to those of ACGs from the same pattern. Furthermore, a dynamic pathway of cellulose degradation was predicted for FPDM based on the results of IC, annotation and prediction of gene function, and taxonomical annotation. The detailed steps were as follows: (1) FUGs and ACGs that displayed no significant difference (P [ 0.05) in relative abundance between treatments, had a low average relative abundance (\ 0.01%) in all treatments, or were not assigned at the genus level were discarded. Moreover, FUGs derived from cellulolytic bacteria were screened out for further analysis.
(2) The remaining FUGs and ACGs were mixed together and clustered based on the abundance variation patterns using STEM algorithm. Abundance variation patterns with fewer than 10 ACGs were removed. FUGs clustered with a group of ACGs (C 60%) that had only one certain function, were presumed to have the same function with them.
(3) To establish the cellulolytic pathway for FPDM, the degradation products determined by IC were first mapped to the cellulolytic pathway according to the order of metabolism. Subsequently, in combination with the pathways and reactions described in the KEGG database, ACGs and putative cellulolytic genes (PCGs) were mapped to the cellulolytic pathway of FPDM based on the annotated and putative functions, respectively.

Isolation, identification and functional analysis of dominant strains
The microbiota FPDM was cultured in liquid filter paper medium at 30°C, and a tenfold serial dilution of the FPDM culture with sterile saline was spread on the solid medium plates consisting of 2% glucose, 0.2% yeast extract, 0.05% MgSO 4 , and 1.5% agar at pH 7.0, and incubated at 30°C for 8 days. Colonies were randomly chosen from the plates and purified by repeating the above-mentioned spread plate method for 5 times. The isolated strains were identified using full-length 16S rRNA sequence analysis as described previously (Hayat et al. 2012).
To confirm synergistic cellulose-degrading capability of two dominant strains, Sporocytophaga sp. CX and Cohnella sp. CX were simultaneously inoculated into liquid filter paper medium with 2% (v/v) inoculum volume for a single strain, and co-fermented at 30°C for 3 days and 8 days, respectively. As a control, a single Sporocytophaga sp. CX strain was cultivated in liquid filter paper medium with 2% (v/v) initial inoculum volume at 30°C for 3 days and 8 days, respectively. The filter paper culture was observed and photographed everyday.
Sporocytophaga sp. CX and Cohnella sp. CX were respectively cultured in liquid medium using filter paper and cellodextrin as the sole carbon source for 48 h. To measure the cellulose-degrading activity, enzyme assay of the culture supernatant was the performed. The endo-1,4-beta-glucanase activity was measured by 3,5-dinitrosalicylic acid (DNS) method as described previously (Miller 1959), and carboxymethyl cellulose (CMC) was used as the substrate (Younesi et al. 2016). The cellobiosidase activity was assayed using p-nitrophenyl-beta-D-cellobioside (pNPC) (Song et al. 2016), microcrystalline cellulose (Avicel PH101) (Kołaczkowski et al. 2020) and phosphoric acid swollen cellulose (PASC) (Dos Santos et al. 2020) as substrate, respectively. The betaglucosidase activity was tested with p-nitrophenylbeta-D-glucopyranoside (pNPG) as substrate as reported before (Song et al. 2016). The cellobiose dehydrogenase activity was assayed according to Sulej et al. (2013). All the assays were performed in triplicate and the mean values ± standard deviation were presented.

Results and discussion
Sequencing data, genome assembly and annotation of microbiota FPDM Overall, * 198.11 million high quality clean pair-end reads (mean * 22.01 million per sample) were generated from the nine microbiota FPDM samples, * 187.79 million (94.79%) of these reads were used to assemble for each corresponding sample and obtained a total of 82,261 scaftigs ([ 500 bp) with a total length of 232.89 megabases (Mb) and an N50 length of 37.35 kilobases (Kb). Subsequently, the unused reads were mixed and assembled, which resulted in 3,473 scaftigs. Considering each assembly, an average of 29,436 ORFs per sample were predicted, and a total of 44,829 unigenes were obtained by clustering all these ORFs, of which 32,614 (72.81%) could be assigned a putative function based on NCBI NR database. In total, 172 CAZyme families and carbohydrate-binding modules (CBMs) were identified, and 1,956 (4.36%) unigenes were assigned to them (Table S1). Specially, almost all of these CAZymes and/or CBMs related genes were present in the members of Firmicutes, Bacteroidetes, Proteobacteria and Actinobacteria (Table S2).
To perform a comprehensive characterization of microbiota FPDM, the MAGs were retrieved from the metagenomic sequences using MetaWRAP pipeline, yielding six near-complete draft bacterial genomes with sizes from 3.41 to 8.06 Mb, GC contents from 36.44 to 69.45%, and estimated genome completeness over 96%. A total of 32,642 ORFs were predicted from these MAGs, of which 20,844 were functionally annotated based on Prokka gene database. The taxonomic assignments via determining average amino acid identity (AAI) based on the type material genomes and NCBI Genome database (Prokaryotes) of MiGA showed that the bins labeled from Bin.1 to Bin.6 were assigned to Achromobacter (Proteobacteria), Herbaspirillum (Proteobacteria), Sporocytophaga (Bacteroidetes), Cohnella (Firmicutes), Microbacterium (Actinobacteria) and Paenibacillus (Firmicutes), respectively (Table S3). CAZyme analysis showed that a lot of CAZyme genes were detected from these bins across the process of cellulose degradation, especially in Bin.3 and Bin.4 (Table S4).
The diversity of microbiota FPDM increased with the increment of filter paper hydrolysate complexity The filter paper degradation process by microbiota FPDM could be divided into three stages: the early stage (3-day cultivation, FPDM.3), the intermediate stage (8-day cultivation, FPDM.8) and the final stage (15-day cultivation, FPDM.15). During this process, the morphology of filter paper changed from the center light yellow strain lawn to all over spread sticky or even translucent strain lawn (Fig. 1). This observation was similar to that reported by Dumova and Kruglov (2009), who observed a slight degradation area in filter paper after the microbiota sample was cultured for 7 days. To further analyze the filter paper hydrolysates at different degradation stages, IC analysis was performed. As shown in Fig. 2, a large number of short-chain saccharide molecules with a degree of polymerization (DP) greater than 15 were detected in FPDM.3 culture. After 8 days of continuous degradation, the short-chain molecules were broken down into cello-oligosaccharides and cellobiose with DP values of 2-6. At FPDM.15 stage, the cello-oligosaccharides and cellobiose contents decreased, and the end product, glucose, was observed. All the results suggest that the product complexity, i.e., the length diversity of polymer chain in hydrolysates, increased during the early-intermediate period of cellulose degradation.
The disturbances in the diversity and structure of microbiota FPDM at different cellulose degradation stages were assessed using Shannon-Wiener and Simpson indices. As shown in Table 1, the species diversity of the microbiota FPDM significantly increased in the intermediate and final stages compared with the early stage (P \ 0.05). In addition, the increasing rate of species diversity in the earlyintermediate period was significantly faster than that in the intermediate-final period (P \ 0.05). The 2-dimensional NMDS plot showed that the microbiota FPDM.3 was well separated from FPDM.8 and FPDM.15, which were relatively close to each other (Fig. 3), suggesting that the structure of FPDM was Fig. 1 The degradation process of filter paper by microbiota FPDM. The microbiota was cultivated in filter paper medium at 30°C for 3, 8 and 15 days Fig. 2 Ion chromatography analysis of the filter paper hydrolysates in the culture supernatant. The supernatant of the 0-day culture was used as a control, and the microbiota FPDM was cultivated in filter paper medium at 30°C for 3, 8, and 15 days. Abbreviations: G represents glucose, C2-C6 refers to cello-oligosaccharides with DP from 2 to 6 disturbed in the early-intermediate period but relatively stable in the intermediate-final period, which coincided well with the alpha diversity results. Considering the changes of filter paper hydrolysates during the degradation process, it could be deduced that the diversity of microbiota FPDM increased with the increment of the length diversity of polymer chain in hydrolysates. Therefore, we hypothesized that the complexity of filter paper hydrolysates affected the way strains interacted with greater cooperation according to the study of Deng and Wang (2016).
Different shift patterns of relative abundance reveal the synergistic action of strains in microbiota FPDM To reveal the ''jungle rule'' in microbiota FPDM due to cooperation and competition for carbon sources in the process of cellulose degradation, the relative abundance of strains in microbiota FPDM and their shifts were analyzed. As shown in Fig. 4a, Bacteroidetes, Firmicutes and Proteobacteria were the main phyla; they accounted for over 96% of the abundance in FPDM samples from different cellulose degradation stages, which was consistent with previous studies (Berlemont et al. 2014;Rosnow et al. 2016). At the genus level, the genera with a total relative abundance greater than 0.5% included Sporocytophaga (40.82%), Cohnella (34.85%), Herbaspirillum (5.87%), Achromobacter (4.99%), Paenibacillus (2.00%), and Microbacterium (0.54%). As illustrated in Fig. 4b, there were 6 relative abundance-shifting patterns among microbial taxa in FPDM. These patterns could be roughly divided into two types: progressive type (pattern 8, 11, 12 and 15) and declining type (pattern 0 and 4). The declining shift microorganisms were deduced to be involved in the early stage of cellulose degradation (cellulose to cellodextrin). Notably, almost all declining taxa belonged to Bacteroidetes, with Sporocytophaga as the representative genus, which has been proved to glide on solid surfaces of cellulose substrates and initiate cellulose degradation (Berg et al. 1972;Liu et al. 2014;Taillefer et al. 2018). On the other hand, progressive shift taxa might play important roles in catabolizing metabolites such as cello-oligosaccharides, cellobiose and glucose that are produced in the intermediate-final degradation stage. Most of the progressive shift taxa were derived from Firmicutes and Proteobacteria. Cohnella, the most dominant genera from Firmicutes, has been reported to contain a set of genes encoding cellulases and possess the ability to utilize cellulose, hemicellulose, and cellobiose (Kämpfer et al. 2006;Eida et al. 2012;Khianngam et al. 2012;Hameed et al. 2013). Therefore, it could be deduced that these carbon source utilizing features of progressive shift and declining shift taxa in FPDM were consistent with IC results, which revealed that cellulose were processively hydrolyzed into short-chain saccharides (DP [ 15), cello-oligosaccharides and cellobiose, and glucose during different degradation stages. Furthermore, all these progressive and declining microorganisms have been identified in cellulolytic microbiota from other samples, such as industrial samples, the gut of native insects (Manfredi et al. 2015), peat-straw (Dumova and Kruglov 2009) and cow manure (Qiao et al. 2019), however, these works focused on mining cellulolytic strains and cellulases, giving few information about the dynamic cooperation among stains during cellulose hydrolysis process.
At present, most researches explored cellulolytic enzymes mainly in annotated gene database produced by omic analysis (Wilhelm et al. 2018;Moraes et al. 2018;Cui et al. 2019;Tomazetto et al. 2020). In this work, 4,031 genes with high relative abundance (C 0.01%) and clear species information were still functionally unannotated. To further explore the FUGs potentially related to cellulose catabolism, cluster analysis was performed using the STEM algorithm based on different abundance variation patterns of genes in different cellulose degradation stages. As shown in Table 2, four typical patterns (0, 4, 11 and 14), in which more than 10 cellulolytic genes were annotated, were selected for clustering. Annotated genes clustered to the declining pattern 0 and 4 usually encode endo-type cellulases (GH5,(8)(9)(10)44,48 and 74) and other non-hydrolytic CAZyme modules (CBM4 and 6) derived from cellulolytic bacteria, such as Sporocytophaga. Therefore, it could be speculated that 363 FUGs from pattern 0 and 4 might be related to early-intermediate stage of cellulose b Fig. 4 Relative abundances and their change patterns of the dominant microbial taxa. a Bar plot of the relative abundances of dominant microbial taxa. Each bar represents the average relative abundance of a treatment. Taxa with relative abundances greater than 0.2% are shown. b Heat map of genus relative abundance of microbiota FPDM at different cellulosedegrading stages. The tree on the left is clustered based on the change pattern of the relative abundance of genus, and the color intensity (log scale) of each panel is proportional to a z-value calculated by visualizing the relative abundance of a genus in the corresponding sample, with reference to the color bar (top right). Each bracketed number following the genus name refers to the change pattern profiles of the relative abundance (middle right   5 Heat map of the relative abundance of cellulosedegrading related CAZyme and CBM families from microbiota FPDM. The tree on the left is clustered based on the change patterns of the relative abundances of CAZyme families. Each panel shows a z-value calculated by normalizing the relative abundance and translating it into color intensity (log scale), with reference to the color bar (top right). Each bracketed number followed CAZyme family refers to the change pattern profiles of the relative abundance (middle right). ''*'' refers to the relative abundance of CAZyme family more than 0.1%, and ''**'' refers to the relative abundance more than 1%. The different shapes after bracketed numbers refer to the cellulose-degrading related enzymes and the EC numbers (bottom right). Specifically, FPDM.3, FPDM.8 and FPDM.15 refer to microbiota FPDM cultivated in fiber medium for 3, 8, and 15 days respectively, and each time points sampled in biological triplicate with the postfix name of ''1'', ''2'', and ''3'' The cellulose-binding function has been demonstrated in one case on amorphous cellulose and beta-1,4-xylan Some of these modules also bind beta-1,3-glucan, beta-1,3-1,4-glucan, and beta-1,4-glucan The cellulose-binding function has been demonstrated in one case on amorphous cellulose and beta-1,4-xylan Some of these modules also bind beta-1,3-glucan, beta-1,3-1,4-glucan, and beta-1,4-glucan   Eida et al. (2012) Dynamic shifts in taxa and CAZymes make the cellulose catabolism pathway more distinct Currently, extensive efforts have been made for clarifying the cellulose catabolism pathway in cellulolytic microbial consortia (Raman et al. 2011;Fosses et al. 2017). Zhou et al. (2014) predicted how a structured and cooperative microbial community functions in biomass conversion process by analyses of reconstructed bacterial draft genomes from a soilderived cellulose-degrading microbial community. Based on microbial draft genomes, Cui et al. (2019) predicted a mutualistic interaction among the dominant microbes regarding the cellulolytic process in the niche. However, all these works built the cellulolytic pathway at a certain culture stage of the microbiota. Few works focused on the dynamically synergistic action at both the gene and genera levels throughout the cellulose degradation process. Using temporal metagenomic analysis, the microbial composition that thrives in hypersaline environments in the tropics has been investigated (Couto-Rodriguez and Montalvo-Rodriguez 2019), which sparked our interest in constructing a more clear cellulose catabolism pathway based on the dynamic analysis of the filter paper hydrolysates, the community structure and composition, and the clustered functional genes of microbiota FPDM. As shown in Fig. 6, two dominant genera, Sporocytophaga and Cohnella dynamically participated in the synergistic degradation process and had different functional tendencies in different stages. In the early-intermediate stage, cellulose was first hydrolyzed into shortchain cellulose, cellodextrin and cellobiose by endo-1,4-beta-glucanase (EC3.2.1.4, GH5,[8][9][10]44,48,and 74), cellulose 1,4-beta-cellobiosidase (nonreducing end) (EC3.2.1.91, GH5 and 9), cellulose 1,4-betacellobiosidase (reducing end) (EC3.2.1.176, GH48) and PCG1 secreted from progressive shift microorganisms (e.g. Sporocytophaga). In the intermediatefinal stages, beta-glucosidase (EC3.2.1.21, GH1,3,5,9,30,39,and 116), exo-beta-1,4-glucosidase Modules of ~120 residues. The cellulose-binding function has been demonstrated in one case on amorphous cellulose and beta-1,4-xylan. Some of these modules also bind beta-1,3-glucan, beta-1,3-1,4-glucan, and beta-1,4-glucan Cohnella bc Eida et al. (2012) FPDM.3.3_16369 NA CBM6 Modules of ~120 residues. The cellulose-binding function has been demonstrated in one case on amorphous cellulose and beta-1,4-xylan. Some of these modules also bind beta-1,3-glucan, beta-1,3-1,4-glucan, and beta-1,4-glucan Was the variation pattern of gene abundance clustered by STEM algorithm, the significant patterns of 0, 4, 11 and 14 with more than 10 annotated cellulolytic genes were displayed, and the top 5 of relative abundance of putative cellulolytic genes and candidate cellulolytic genes were listed, respectively b Represents the species that carries the cellulose-degrading related genes detected in microbiota FPDM, including cellulases and exoglucanases c Represents the species that carries the cellobiose-degrading related genes detected in microbiota FPDM, including beta-glucosidases and phosphocellobiase Putative cellulolytic genes were highlighted in gray (EC3.2.1.74, GH1, 3, 5, and 9) and PCGs secreted from both progressive and declining shift microorganisms continued to hydrolyze short-chain cellulose and cellodextrin (DP [ 15) into shorter cellooligosaccharides and cellobiose. Then, cello-oligosaccharides and cellobiose were completely hydrolyzed to glucose and glucose 6-phosphate by beta-glucosidase (EC3.2.1.21, GH1,3,5,9,30,39,and 116), exobeta-1,4-glucosidase (EC3.2.1.74, GH1, 3, 5, and 9), phosphocellobiase (EC3.2.1.86, GH1 and 4) and PCG2 produced from declining shift microorganisms (e.g. Cohnella).
To further identify the biological function of two dominant genera, Sporocytophaga sp. CX and Cohnella sp. CX were isolated from microbiota FPDM and co-cultured using filter paper as the substrate. The activities of endo-1,4-beta-glucanases (EC3.2.1.4, GH5, 8-10, 44, 48 and 74) (250 ± 32 U/L) could be detected in supernatant of Sporocytophaga sp. CX, and Cohnella sp. CX showed the activities of betaglucosidases (EC3.2.1.21, GH1,3,5,9,30,39,and 116,and EC3.2.1.74,GH1,3,5, and 9) (14.6 ± 0.3 U/L) and cellobiose dehydrogenases (EC1.1.99.18, AA3) (39.0 ± 3.3 U/L) (Fig. 6). Activities of other cellulases were not detected probably due to the low secretion level. As shown in Fig. 7, when Sporocytophaga sp. CX was cultivated for 8 days, filter paper was partly degraded. However, when this strain was co-fermented with Cohnella sp. CX, a completely hydrolysis of filter paper could be observed. All these results confirmed the synergistic work of the two strains deduced by our metagenomic analysis. Fig. 6 Cellulolytic pathway of microbiota FPDM. Different colored dots represent different metabolites. Arrows indicate both the pathways and directions of metabolism. The numbers in colored boxes refer to the EC number of the enzymes. Numbers in brackets represent the number of genes identified. Data highlighted in gray refer to the enzyme activity, and values represent means ± standard error. PCG1 and PCG2 represent putative cellulolytic genes with different functions. Orange and light green wedges represent the pathways of the declining shift taxa (e.g. Sporocytophaga) and the progressive shift taxa (e.g. Cohnella) respectively. The width of the wedges roughly represents the species' contribution rate to the metabolic pathway parallel to the corresponding position

Conclusions
The cellulolytic function of microbiota FPDM was clarified by the dynamic analysis of hydrolysates and metagenomes at different cellulose-degrading stages. The biodiversity of FPDM gradually enhanced as the length diversity of polymer chain in cellulose hydrolysates increased. FPDM-mediated cellulose degradation was accomplished by the collaboration and competition of two types of taxa, represented by Sporocytophaga and Cohnella respectively. Based on the deduced cellulolytic function of 432 annotated and 594 unannotated genes, a dynamic cellulose-degrading pathway of FPDM dominated by two core genera was predicted. Our work systematically depicted the synergistic and dynamic action of different taxa and their enzymes in microbiota during cellulose-degrading process. Significantly, based on the temporal analysis of cellulose-degrading microbiota samples, functional strains and proteins belong to the same shift pattern could be clustered, which should be valuable for mining key cellulolytic strains and enzymes and rebuilding more efficient cellulose-degrading system (Fig. 8). Fig. 7 The degradation process of filter paper by dominant strains isolated from microbiota FPDM. The strains were cultivated in liquid culture medium containing round filter paper at 30°C. Sporocytophaga sp. CX was fermented for Acknowledgments FY, XL, and XC were supported by grants from the National Natural Sciences Foundation of China (31671796, 31771907, and 31801469), and the Liaoning BaiQianWan Talents Program is also gratefully acknowledged.
Author contributions MY: designed the strategy for metagenomic analysis, developed the scripts and wrote the manuscript. JZ: contributed with reagents, materials, and analyses, performed the experiments, and wrote the manuscript. YY, XC: contributed with reagents, materials, and performed the experiments. FY, XL: drafted, wrote and revised the manuscript. All authors agreed with the submitted version of the paper. All authors read and approved the final manuscript.
Funding This work was supported by the National Natural Sciences Foundation of China (31671796, 31771907, and 31801469), and the Liaoning BaiQianWan Talents Program is also gratefully acknowledged.
Data availability The metagenomic reads generated in this study have been submitted to NCBI Sequence Read Archive (SRA) under the accession number SRP255666. A previous version of this paper was published as preprint in Research Square (https://doi.org/10.21203/rs.3.rs-22654/v1) (Yang et al. 2020).