A genomic perspective on the potential of termite-associated Cellulosimicrobium cellulans MP1 as producer of plant biomass-acting enzymes and exopolysaccharides

Isolation of symbiotic cellulolytic bacteria

Drywood termites (Cryptotermes domesticus) were collected from rotten tree trunks and bagasse in Nghe An Province, Vietnam. Termites were surface-sterilized using 70% ethanol to remove contamination and then washed with sterile distilled water. The head of each termite was removed using forceps; the entire guts were removed, crushed with glass rods, and subsequently inoculated into 1 ml broth mineral medium M1 (NaNO₃ 2.5 g; NaCl 0.1 g; KH₂PO₄ 2 g; MgSO₄ 0.2 g; CaCl₂.6H₂O 0.1 g, pH 7.0 in a liter) containing 1% carboxymethylcellulose (CMC) or filter paper as a sole carbon source (Gupta, Samant & Sahu, 2012). These cultures were then incubated for 14 days in an incubating sharker at 30 °C. To isolate the cellulolytic bacteria, the growing cultures were spread on the M1 plate medium (K₂HPO₄ 1 g; NaNO₃ 2.5 g; KCl 2 g; peptone 2 g; MgSO₄ 0.5 g; CMC 10 g; agar 15 g; pH 7.0 in a liter). All bacterial isolates were subsequently purified by re-streaking on the M2 agar plate. Confirmation of the cellulolytic ability of pure isolates was performed on the solid medium by covering the Petri dishes with Congo-red dye (Teather & Wood, 1982). The colonies showing yellow-colored halo zones by Congo-red staining were considered as positive cellulolytic bacteria and the clear zones were measured.

Physiological, biochemical, and 16S rRNA sequencing analysis

The shape and size of strain MP1 were determined by scanning electron microscope (SEM) JSM-5410 (JEOL, Tokyo, Japan). Gram staining was performed using the conventional methodology and confirmed using the KOH test (Powers, 1995). The effects of different conditions (ranges of pH, NaCl concentration and temperature, carbon and nitrogen sources) on the growth were investigated as previously described (Kamlage, 1996).

Genomic DNA for 16S rRNA gene sequencing was prepared by phenol-chloroform extraction. The amplification of the 16S rRNA gene sequence of strain MP1 was performed by using the universal primer pair 27F (5′ -TAACACATGCAAGTCGAACG-3′) and 1429R (5′-GGTGTGACGGGCGGTGTGTA-3′) (Phi et al., 2010). A sequence similarity search was carried out using the BLAST program (http://blast.ncbi.nlm.nih.gov/Blast.cgi). The phylogenetic tree was computed by using the neighbor-joining method with 1,000 bootstrap in MEGA version 6.0 (Tamura et al., 2013). Numbers at nodes indicate percentages of 1000 bootstrap re-samplings and Bifidobacterium bifidum DSM 20456 (S83624) was used as the outgroup branch. The 16S rDNA gene sequence of strain MP1 was deposited onto the GenBank (NCBI) under accession number MW534740.

Genome sequencing, de novo assembly, and annotation

For library construction, DNA was extracted from a pure culture of a single colony of strain MP1 using G-spin™ Total DNA Extraction Mini Kit according to the manufacturer’s instructions. The quantity and quality of extracted DNA were evaluated by electrophoresis in 0.6% (w/v) agarose gel and NanoDrop spectrophotometer 2000 Thermo Scientific in order to construct a paired-end library. The constructed genome library was then sequenced using the Illumina platform (Illumina, California, USA). Quality control and read trimming were conducted using FastQC version 0.11.5 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc) and Trimmomatic version 0.36 (Bolger, Lohse & Usadel, 2014). The de novo genome assembly was made with SPAdes v.3.13 (Bankevich et al., 2012), which was then analyzed for its completeness using the Benchmarking Universal Single-Copy Orthologous (BUSCO) version 3 (https://gitlab.com/ezlab/busco).

The draft genome assembled into contigs was annotated using Prokaryotic Genomes Annotation Pipeline (PGAP; http://www.ncbi.nlm.nih.gov/genome/annotation_prok/) at NCBI. In addition, the CRISPRCasFinder was used to identify putative CRISPR loci and Cas cluster (Grissa, Vergnaud & Pourcel, 2007). Orthologous genes and Gene ontology (GO) were analyzed using clusters of orthologous genes (COGs) (Galperin et al., 2015) and InterProScan 5 (Jones et al., 2014), respectively. Virulence factors-encoding genes were identified using the Pathosystems Resource Integration Center (PATRIC) platform (Wattam et al., 2017). The graphical map of the circular genome was also generated using PATRIC. The draft genome sequence was deposited in the GenBank (NCBI) database under accession number: JAFGYF000000000.

Analysis of whole-genome similarity

To classify strain MP1 at the species level, whole-genome similarity including average nucleotide identity (ANI) calculation and digital DNA–DNA hybridization (dDDH) was performed. The ANI was calculated using the orthologous average nucleotide identity (OrthoANI) (Lee et al., 2016). The MP1 genome sequence data was uploaded to the Type (Strain) Genome Server (TYGS) for a whole-genome-based taxonomic analysis (https://tygs.dsmz.de). In silico dDDH and the branch lengths were and the Genome BLAST Distance Phylogeneny evaluated using Genome-to-Genome Distance Calculator (GGDC) (Meier-Kolthoff et al., 2013).

Comparative genomics and prediction of Carbohydrate-active enzyme

A genome-wide comparison of COGs between the assembled genome of Cellulosimicrobium cellulans MP1 and three other C. cellulans available in GenBank, including J36 (NZ_JAGJ01000000.1), LMG16121 (NZ_CAOI01000000.1), and ZKA 48 (NZ_QUMZ01000000.1) was implemented using the OrthoVenn web server with default parameters (E-value 1e−5 and inflation value 1.5) (Wang et al., 2015). The putative genes encoding CAZymes in the comparative genomes were predicted using the dbCAN2 meta server and classified by DIAMOND, HMMER, and Hotpep via CAZy and dbCAN databases, respectively. The top hits with e-value <1E−17, minimum homology rate >50%, and coverage >45% were considered to be homologs.

Determination of extracellular enzymatic activities

To verify the production of endoglucanase, C. cellulans MP1 was cultivated in the TN medium (rice bran 20 g; soya peptone 10 g, casein 10 g, NaCl 1 g, pH 7.0 in a liter) at 37 °C with vigorous sharking. At different time intervals, samples were taken, followed by centrifugation at 4 °C, 12,000 rpm for 10 min to remove bacterial cells and insoluble materials from the culture broth. About 0.5 ml of the crude enzyme solution was added into 0.5 ml of 0.05 M sodium phosphate, pH 7.0 buffer containing 1% of CMC. The mixture was then incubated at 30 °C for 30 min, which was stopped by adding 1 ml of 3,5′–dinitrosalicylic acid (DNS) reagent followed by boiling the reaction mixture at 100 °C for 5 min (Miller, 1959). As for xylanase production, strain MP1 was cultivated on M1 medium (K₂HPO₄ 1 g; KCl 2 g; NH₄Cl 2.5 g; yeast extract 2 g, MgSO₄ 0.5 g, sugar cane bagasse 5 g, pH 7.0 in a liter) at 37 °C for 4 days. At different time intervals, culture was centrifuged at 4 °C, 12,000 rpm for 10 min to obtain the crude enzyme. About 0.5 ml of the crude enzyme was mixed with 0.5 ml of the substrate solution and incubated at 50 °C for 5 min, followed by adding 1 ml DNS reagent. Xylanase activity was assayed using 1% birchwood xylan solution in acetate buffer, pH 5.0 as the substrate (Bailey, Biely & Poutanen, 1992). To quantify pectinase activity, the M2 medium (NH₄Cl 2.5 g; K₂HPO₄ 1 g; MgSO₄ 0.5 g; KCl 2 g; yeast extract 2 g, rice bran 10 g; lactose 10 g; pH 7.0 in a liter) was used. The pectinase assay was performed at 40 °C for 30 min as described previously (Mercimek Takcı& Turkmen, 2016). The quantitative enzyme assays of endoglucanase, xylanase, and pectinase were performed according to standard IUPAC procedures and expressed as international units (IU) (Ghose, 1987). One unit (IU) of enzymatic activity is defined as the amount of enzyme that releases 1 µmol reducing sugars per ml of culture supernatant per minute under assay conditions.

Levan determination

To test the ability to synthesize levan, the C. cellulans MP1 was cultured overnight in batch culture medium (yeast extract 7 g, (NH₄)₂SO₄ 1.5 g, K₂HPO₄ 2.5 g, pH 7.0 in a liter) at 37 °C under vigorous agitation. The overnight culture was transferred to a new batch culture medium supplemented with 100 and 200 g/l sucrose and adjusted to an optical density at 600 nm of 0.1. The culture was centrifuged at 10,000 rpm for 10 min at different time intervals. Cell-free supernatant was used to determine levan produced during fermentation (Gojgic-Cvijovic et al., 2019). Levan was harvested by adding three volumes of ice-cold ethanol. The mixture was kept at 4 °C for 12 h, centrifuged at 12,000 rpm at 4 °C for 20 min, and washed with 75% ethanol to remove the residual sugars. The obtained precipitate was hydrolyzed with 0.1 M HCl at 100 °C for 1 h. After neutralizing with 2 M NaOH, levan content was determined according to Somogyi and Nelson (Somogyi, 1945).

Cellulose-degrading potential and identification of the strain MP1

Of the 65 isolates able to grow on CMC plates as the sole carbon source, a total of 8 bacterial isolates produced variable zones of CMC clearance after Congo-red staining. Among those, MP1 was selected due to showing the maximum zone of clearance (9 ± 1, 1 mm). In support of this result, enzyme assays for cellulase activity on CMC and filter paper were found to be the highest for MP1 with 0.66 ± 0.15 IU/ml after 48 h and 0.33 ± 0.05 FPU/ml after 72 h of cultivation, respectively (Table S1). These results indicated that isolate MP1 is a potent cellulolytic bacterium for further study.

The phenotypic examinations indicated that MP1 grew well on LB agar after 2 days of incubation at 37 °C, producing colonies that were circular, smooth, convex, and pale yellow in colour. MP1 cells were Gram-positive, non-spore-forming rods, and nonmotile. The strain MP1 could utilize D-galactose, D-sucrose, D-raffinose, and amygdalin as sole carbon and energy (Table S2). It gave a positive test for catalase, starch hydrolysis, β-galactosidase, and nitrate reduction, whereas negative for oxidase, gelatinase, urease, indole, H₂S, acetoin. The 16S rDNA gene sequence of MP1 was aligned with the similar nucleotide sequences in the GenBank database in which a phylogenetic tree was then constructed. The neighbor-joining phylogenetic tree showed the close relationship between MP1 and related Cellulosimicrobium species and the highest homology to Cellulosimicrobium cellulans DSM 43879 (99.5%) and Cellulosimicrobium funkei A153 (99.6%) (Fig. 1). Further, the OrthoANI software was performed to determine the OrthoANI value between isolate MP1 and five closely related Cellulosimicrobium species. It revealed that MP1 shared high similarity to C. cellulans DSM 43879 (88.71%), C. funkei JCM 14302 (91.35%), and C. funkei JCM NRBC 104118 (91.35%), and low nucleotide similarity was observed with Promicromonospora sukumoe SAI-064 that were out of distinct sub-clade (Fig. 2A). This result confirmed that this strain was not considered to be a novel species.

To make identification more accurate at the species level, the whole-genome-based taxonomic analysis conducted by the Type Strain Genome Server (TYGS) platform suggested that Cellulosimicrobium sp. MP1 was closest to C. cellulans DSM 43879 with digital DNA-DNA hybridization (dDDH) values and differences in guanine-cytosine (GC) content of 57.5% (formula d₆) and 0.56%, respectively (Fig. 2B). As shown previously, the differences between C. funkei and C. cellulans in phenotypic characteristics are motility and ability to utilize raffinose (Hamada et al., 2016; Yoon et al., 2007). Based on the phenotypic characteristics and genome-wide comparison, this strain was identified as Cellulosimicrobium cellulans MP1. This bacterium was deposited at VAST-Culture Collection of Microorganisms (VCCM) with the accession number VCCM 14150.

Genome sequence and general features of the C. cellulans strain MP1

Briefly, the standard short insert library yielded 443,986,169 bases resulting in 4,487,842 mapped reads with about 92.2-fold sequencing depth. The draft genome of strain MP1 (4,580,223 bp with GC content of 73.9%) was assembled into a ring chromosome, comprising of 23 contigs and no plasmid was detected (Fig. 3). The genome was predicted to have 4,088 genes assigned for 3,964 protein-coding sequences (CDS), 61 rRNA sequences, 55 tRNA sequences, 3 ncRNA sequences (Table 1). Moreover, a total of three virulence factors present in C. cellulans MP1 included dihydroxy-acid dehydratase, CarD-like transcriptional factor, and two calmodulin-like proteins (Fig. 3).

A total of 79.5% (3,216 out of 4,046) of the protein-coding sequences were assigned to 21 out of 25 COG functional categories (Fig. 4A). Transcription (K: 365 protein-coding sequences), carbohydrate transport and metabolism (G: 307), amino acid transport and metabolism (E: 224), energy production and conversion (C: 199), and inorganic ion transport and metabolism (P: 180) were found to be the largest categories. By contrast, the least represent groups included cell cycle control, cell division, chromosome partitioning (D: 33), chromatin structure and dynamics (B: 2), and cell motility (N: 1).

GO analysis was used to provide a deeper understanding of the functional catalogs of strain MP1. A total of 3,165 genes were assigned to 43 subclasses, including 13 subclasses of the cellular component (CC) class, 10 subclasses of the molecular function (MF) class, and 20 subclasses of the biological process (BP) class (Fig. 4B). In detail, the CC class occupied the most genes (1638 genes; 45.3%), followed by the BP (1424 genes; 39.4%), and MF (553 genes; 15.3%) class. The most abundant pathways were cell (GO:0005623; 508 genes), cell part (GO:0044464; 508 genes), and membrane (GO:0016020; 340 genes), which could be considered as the main functional groups of genes belonging to the CC class. Within the BP class, the three most prevalent molecular functions were metabolic process (GO:0008152; 378 genes), cellular process (GO:0009987; 379 genes), and growth (GO:0040007; 281 genes) (Fig. 4B). Among all subclasses belonging to the MF class, catalytic activity (GO:0003824) and binding (GO:0005488) contributed the most genes, with 311 and 176, respectively.

Genome-wide comparison of COGs and Carbohydrate-active enzymes

Given that microbial evolution is due to vertical descent from a single ancestral gene leading to orthologous genes in different species, it is necessary to explore gene function, gene structure, and molecular evolution by using a genome-wide comparison of COGs in different strains. The COGs of C. cellulans MP1 were compared with four other strains, including J36, LMG16121, and ZKA48. It showed that C. cellulans MP1 comprised of 3405 COGs and 794 singletons. C. cellulans J36 and ZKA48 included 3604 and 3605 COGs, respectively, whereas the lowest COGs were observed in LMG16121 (Fig. 5A). Surprisingly, strain MP1 had the largest number of singletons (n = 794). Venn diagram denoted that a total of 2,539 COGs were commonly shared by all four strains of C. cellulans. Among the unique COGs observed in all strains, strain MP1 had the largest number of 20 (Fig. 5A). To support this result, the REALPHY phylogeny builder web tool was used to compare C. cellulans genomes. A maximum-likelihood phylogenetic tree showed that MP1 and LMG16121 formed a monophyletic clade, suggesting a strong evolutionary relationship (Fig. S1).

Using dbCAN 2 meta server, the putative genes encoded for CAZymes present in C. cellulans MP1 were screened to find out the genes responsible for cellulose and hemicellulose degradation. After removing the sequences that did not meet the filtering criteria, a total of 195 predicted CAZymes was identified corresponding to 4.8% of the total of 3,216 protein-coding sequences observed in this strain. Glycoside hydrolases (GHs) involved in the degradation of the most plant biomass such as cellulose and hemicellulose were predicted to be the most abundant subfamily with 99 enzymes. Next, 39 glycosyltransferases (GTs), 37 carbohydrate-binding modules (CBMs), 11 carbohydrate esterases (CEs), 6 enzymes with auxiliary activities (AAs), and 3 polysaccharide lyases (PLs) were detected (Fig. 5B).

The genome of C. cellulans MP1 was compared to closely related C. cellulans such as C. cellulans J36, C. cellulans LMG16121, C. cellulans ZKA48. Generally, C. cellulans MP1 had significantly more CAZyme domains found than other C. cellulans strains. Despite forming a monophyletic clade with strain LMG16121, MP1 was 56 CAZyme domains higher than LMG16121. Meanwhile, with the exception of LMG16121, the family numbers of GT, CBM, and AA family numbers were the same in all compared strains, suggesting the coexistence of these genes in breaking down cellulosic biomass (Fig. 5B).

Mining of plant biomass-acting enzymes in C. cellulans MP1 gene pool

CAZyme analysis suggested that 30 cellulose-related sequences were detected in the genome of C. cellulans MP1. The major families related to the degradation of cellulose are GH6, GH9, GH48, GH10, GH16, GH1, GH3, GH13, and GH64. Based on the annotation, 5 endoglucanases (two GH6 and three GH9), 3 exoglucanases (GH6, GH10, and GH48), and lichenase (GH16) were revealed in the genome of strain MP1 (Table 2). A total of five out of the eight annotated endoglucanases and exoglucanases (orf_454, orf_1616, orf_3244, orf_1607, orf_1611) contained a CBM2 domain appended to them, which have been known to bind to crystalline cellulose, insoluble chitin, and xylan (McLean et al., 2002). CBM2 was the major CBMs present in the genome of strain MP1. In the β-glucosidases family, three GH1 and six GH3 were considered as other important members for cellulose degradation and three GH13 had α-glucosidase activity. It is interesting to note that two pectate lyases such as PL4 (orf_102) and PL1 (orf_1326) belonging to the large class of PL were identified (Table 2). Given that pectin is involved in providing structural support for plant such as cell walls, wall integrity, and cell–cell cohesion, pectate lyases (EC 4.2.2.2) belonging to pectinase catalyze the eliminative cleavage of α-1,4-glycosidic bonds between C₄ and C₅ of pectin or pectic acid, producing unsaturated methyloligogalacturonates (Abbott, Gilbert & Boraston, 2010; Hugouvieux-Cotte-Pattat, Condemine & Shevchik, 2014). Pectate lyases are important symbiotic bacteria through facilitating their growth in presence of highly pectinolytic bacteria and under pectin-rich environments (Hugouvieux-Cotte-Pattat, Condemine & Shevchik, 2014).

For hemicellulose degrading genes mining, a total of 21 annotated proteins were deduced to involve in hemicellulose degradation. The major families responsible for the breakdown of hemicellulose were GH43, GH31, GH127, GH51, GH10, GH36, and GH4, and the total number of enzymes, including all families, was 21 (Table 2). Among them, the most abundant enzymes were attributed to GH43 (9 enzymes), followed by GH124 (3 enzymes) and GH31 (2 enzymes), indicating that strain MP1 might have great potentials for degradation of hemicellulosic backbones or debranding hemicellulose. Hemicellulose-related genes were predicted as mannan endo-1,4-β-mannosidase, arabinoxylan arabinofuranohydrolase, α-xylosidase, α-L-arabinofuranosidase, β-L-arabinofuranosidase, endo-1,4-β-xylanase, α-galactosidase, and bifunctional β-xylosidase/ α-arabinosidase (Table 2). Only two enzymes, mannan endo-1,4-β-mannosidase and α-galactosidase, have the ability to hydrolyse parts of mannan, the second major component in hemicellulose. Two endo-1,4-β-xylanase (orf_100 and orf_3772) and one bifunctional β-xylosidase/ α-arabinosidase (orf_3988) were coupled with either CBM2 or CBM9. Surprisingly, genes encoding arabinofuranosidase (orf_3999, orf_4000, orf_4001 and orf_4002), arabinoxylan arabinofuranohydrolase (orf_4003), and bifunctional β-xylosidase/ α-arabinosidase (orf_4004) are clustered in an operon. However, AA enzymes were not found to be related to cellulose and hemicellulose degradation.

To confirm the ability to produce endoglucanase, xylanase, and pectinase of strain MP1, the quantitative enzyme assay was monitored across different incubation periods (24 h, 48 h, 72 h, 96 h, 120 h) (Fig. 6). Using rice bran as an inducible substrate for endoglucanase production, the lowest endoglucanase activity of 0.15 ± 0.01 IU/ml was observed at 24 h, followed by an increase of activity within 48 h to 96 h. The highest endoglucanase production of about 3.17 ± 0.10 IU/ml was achieved after 72 h of incubation. The reduced activity of endoglucanase was notified as 1.40 ± 0.17 IU/ml at 120 h. Xylanase activity of strain MP1 was monitored on M1 medium supplemented with sugar cane bagasse at 37 ° C for 120 h. The result was a steep increase within 0 to 72 h of incubation times. The highest enzyme production was recorded during the stationary phase reaching maximum (1.84 ± 0.08 IU/ml) at 96 h, followed by a slight decrease after 120 h (Fig. 6C). Regarding to pectinase, the maximum pectinase yield was 0.16 ± 0.01 IU/ml after 48 h of incubation. Increasing incubation time was subjected to a significant decrease of pectinase activity. Despite the fact that a number of studies involved in the yields of endoglucanase, xylanase, and pectinase have been performed, it is hard to compare due to the effect of production conditions, substrate and assay conditions, and the way of defining the units.

Levan exopolysaccharide biosynthesis and degradation

Apart from cellulose-related genes, we were interested in the gene encoding for levansucrase (orf_22) since the ability of the genus Cellulosimicrobium to produce levan has not been studied yet. Levansucrase (EC 2.4.1.10) belonging to GH32 participates in synthesizing levan using high sucrose concentration as the main substrate (Raga-Carbajal et al., 2018). The C. cellulans MP1 ls gene consists of 1,851 nucleotides encoding a protein of 616 amino acids with the predicted molecular mass of 66.6 kDa. Alignment of the amino acid sequence of C. cellulans MP1 Ls shows the highest level of identity with Gluconacetobacter diazotrophicus (61%) and Microbacterium saccharophilum (58%) (Fig. S2). In addition, levansucrase is clustered alongside levanbiose-producing levanase (levB-orf_20) and levanases (sacC 1-orf_19 and sacC 2-orf_21) in an operon (Fig. 7A). This result suggested that levan produced by C. cellulans MP1 may be hydrolyzed into L-FOs such as levanbiose and levanases. It is noteworthy that the sacC1-levB-sacC2-ls operon was not present in the compared C. cellulans strains

The levan production in media with different initial sucrose concentrations was monitored as shown in Fig. 7B. Levan production depended on cultivation time and sucrose concentration. In the medium with 100 g/l sucrose, the maximum levan concentration was achieved at 20 h (9.9 ± 0.11 g/l), followed by a significant decrease after 24 h (5.5 ± 0.02 g/l). When sucrose concentration was increased to 200 g/l, the highest yield of levan was 14.8 ± 1.2 at 20 h. It was noteworthy that levan was accumulated in the stationary phase.