Ruminiclostridium thermocellum M3 strain was isolated and enriched from horse manure by Sheng et al. (2016) and deposited in the Microbiology Laboratory of the School of Environment and Chemical Engineering, Heilongjiang University of Science and Technology. The seed was stored in a constant temperature incubator at −20°C and cultivated in an anaerobic bottle (filled with nitrogen) containing modified ATCC 1191 (MA) medium. The main components of culture medium were K₂HPO_4, 1.5 g/L; MgSO₄·7H₂O, 0.2 g/L; (NH₄)₂SO_4, 1.0 g KCl, 0.2 g/L; L-cysteine, 0.5 g/L; KH₂PO_4, 3.0 g/L; CaCl₂•2H₂O, 0.025 g/L; NaCl, 1.0 g/L; Yeast, 1.5 g/L; and Avicel, 5.0 g/L. The temperature of the culture was maintained at 60°C at 120 rpm.

The R. thermocellum samples for whole genome sequencing analysis were cultured in MA medium for 24 h at 60°C, centrifuged for 5 min at 4°C, 12,000 × g, then the cell pellet was washed twice with normal saline (NS), and crushed with the FastPrep-24 instrument in lysing matrix B tubes (MP Biomedical) for 40 s to release the genomic DNA from cells. The extraction of genomic DNA was performed with the Tiangen bacterial DNA mini kit (Tiangen Biotech Co. Ltd., Beijing, China) according to the manufacturer’s protocol. The harvested DNA was detected using agarose gel electrophoresis and then quantified by Qubit 4.0 (ThermoFisher, Q33226). The genomic sequencing was conducted by the Pacbio sequencing platform with de novo assembly (SMRT portal; Berlin et al., 2015).

The genome annotation and gene function were predicted in the Gene Ontology (GO) database, the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, the Clusters of Orthologous Groups (COG) database, and the Non-Redundant Protein databases (NR). The whole genome BLAST search (E-value below 1^e−5, minimal alignment length percentage above 40%) was performed with the above four databases. The prediction of carbohydrate-active enzymes was conducted with the Carbohydrate-Active enZYmes Database (Lombard et al., 2013).

Genome component prediction of R. thermocellum M3 strains included the coding gene, repetitive sequences, signal peptide, genomic islands, prophage, lipoprotein, prophage, and clustered regularly interspaced short palindromic repeat sequences (CRISPR). The GeneMarkS program was conducted to retrieve the coding genes (http://topaz.gatech.edu/; Yang et al., 2020). The tRNA was predicted by the Aragorn program (Laslett and Canback, 2004), the rRNA was predicted by the RNAmmer program (Lagesen et al., 2007), and the miscRNA was predicted by Infernal (v1.1.2; Mostajo Berrospi et al., 2019). Repeat Modeler was used to predict the repeat sequence de novo of the assembly results, and the RepeatMasker program was used for identifying repetitive elements in nucleotide sequences.

Genomic islands were predicted by the Island Path-DIOMB program (Hsiao et al., 2003). The CRISPR was predicted by the CRISPR recognition tool (CRT; Grissa et al., 2007). Prophages were predicted using PhiSpy (Hähnke et al., 2009). RepeatMasker was used to identify the location and frequency of repeats on the genome (Saha et al., 2008). NCBI Blast+ was used to compare the protein sequences with CDD, KOG, COG, NR, NT, PFAM, Swissprot, TrEMBL, and other databases to obtain the functional annotation information.

The core genes, specific genes, the gene family phylogenetic tree, single nucleotide polymorphism (SNP), and genome visualization were analyzed to reveal the result of comparative genomics (Zhong et al., 2018). The genome sequences of R. thermocellum DSM 2360, R. thermocellum ATCC 27405, R. thermocellum DSM 1313, and R. thermocellum AD2 were obtained from the NCBI database. Genomic alignments among R. thermocellum M3 and other genomes of R. thermocellum were performed using the MUMmer and LASTZ tools (Kurtz et al., 2004). Core genes and specific genes were analyzed using the CDHIT rapid clustering of similar proteins software with a threshold of 50% pairwise identity and a 0.7 length difference cutoff for amino acids (Li et al., 2001, 2002; Li and Godzik, 2006). The relationships between five R. thermocellum strains were analyzed, and the results were represented in a Venn diagram. NCBI Blast+ was used to compare the predicted 16S rRNA sequence with the NCBI 16S database to obtain its homologous strain information and a phylogenetic tree was constructed using mega software (Hall, 2013).

The whole genome of R. thermocellum M3 was sequenced and analyzed with regard to the predictions of coding genes. The identified total size of the genome R. thermocellum M3 was 3,602,270 bp with 39% GC content using SPAdes, and the number of coding genes was 3,195 with an average gene length of 973.56 bp (Supplementary Figure S1). The major characteristics of the genome of the R. thermocellum M3 and four R. thermocellum are summarized to acquire more perceptions concerning the genetic information (Table 1). The final assemblies indicated that the genomes of five R. thermocellum were similar in size and G + C contents, genome annotation yielded 3,062, 3,196, 2,949, 2,959, and 3,077 genes for strains DSM 2360, ATCC 27405, DSM 1313, AD2, and M3, respectively.

The result of genome annotation indicated that a series of genes encoding virulence factors (VFDB), antibiotic resistance (CARD), and pathogen-host interactions (PHI-base) were present in R. thermocellum (Supplementary Table S1). Meanwhile, R. thermocellum strains had several genes related to carbohydrate degradation, biosynthesis, and modification enzymes. The whole-genome visualization map clearly identified the composition, location, and function of the R. thermocellum M3 genome, which is demonstrated in Figure 1, with 100% coverage of sequencing. Detailed information includes analysis of CDS, non-coding RNA, COG, functional classification of genes, and the size of genome and GC%. Additionally, among the 3,077 predicted genes of R. thermocellum M3, only 2,071 CDSs were assigned any COG, which accounted for 67.3% of the predicted genes, and an additional 171 CDSs were assigned to group S with an unknown function. There are 137 genes involved in amino acid production and transformation (group E) in the aligned COG database genes, accounting for 6.62% of the COG annotations; the number of carbohydrate transport and metabolism-related genes (group G) was 130, accounting for 6.28% of the COG-annotated genes; the quantity of energy production and conversion-related genes (group C) was 107, taking up 5.17% of the COG-annotated genes; and only 16 genes fall into group Q (secondary metabolites biosynthesis, transport, and catabolism; see Supplementary Figure S2).

Gene Ontology and KEGG databases were used to obtain the elucidation of the “character” of the coding gene in the bacteria from a macroscopic perspective. In the genome of R. thermocellum M3, 2021 genes were annotated in the GO database (Table 2). In the biological process, the metabolic and cellular processes account for the highest proportion; in molecular function, it is mainly related to catalytic activity and binding; in the cellular component, it is closely compared to cells and cell parts, which indicated that strain M3 has more proteins involved in metabolism, cell composition, and enzymatic catalysis. These results are in accordance with the biochemical characteristics of M3, which possesses an outstanding catalytic capacity for cellulose substrate. In the KEGG database, 1,246 genes were annotated, which can be divided into five branches according to the metabolic pathways involved in the genes: cellular processes, environmental information processing, genetic information processing, metabolism, and organismal systems. As shown in Supplementary Figure S3, aligning these annotated genes into metabolic pathways yielded a total of 147 metabolic pathways. It can be clearly seen that the predominant pathway was carbohydrate metabolism and overview with 329 and 229 unigenes, followed by amino acid metabolism and energy metabolism with 197 and 153 unigenes.

Meanwhile, the predicted numbers of SignaIP-TM and SignaIP-noTM were 24 and 106, respectively, among the total Signa proteins. A total of four CRISPR arrays were found by CRT (Supplementary Table S2). The number of repeat counts in four arrays was 51, 90, 134, and 144. Moreover, several cas or cas-like genes were found in their neighborhoods. This suggests that R. thermocellum M3 has a defense against phage contamination, as CRISPR is very important in prokaryotes and is involved in resisting foreign phages and plasmids and recognizing and silencing invading functional elements. The results of gene-island (GI) prediction obtained eight GI with an average G + C content of about 36.4%, which is slightly lower than the G + C content of the M3 genome. It is worth noting that the G + C content of GI5 (33.7%) is significantly different than that of other GI, which indicates GI5 may be an exogenous sequence by horizontal transfer. The main components of the exogenous sequence are the transposase protein and the hypothetical protein (Supplementary Table S3).

To reveal the genetic and evolutionary relationships between R. thermocellum M3 and other typical R. thermocellum, the analysis in view of the core-pan gene of the whole genome sequence was conducted. The Gene family boxplot (Supplementary Figure S4) revealed that the pan-genome trend of R. thermocellum is open with the increasing number of R. thermocellum strains sequenced. The open pan-genome indicated that there is a significant capacity for discovering novel genes with the evolution and development of strains. Among the five R. thermocellum strains, the number of genes in the pan-genome is 3,042 and includes a core gene set (2,544 genes; Figure 2). Dispensable genes and unique genes existed across all genomes of five R. thermocellum strains. Among them, strain ATCC 27405 had the largest number of specific genes (267 genes), followed by strain M3 (192 genes). The high resemblance of five R. thermocellum strains was reflected by the large proportion of core genes.

To identify the functional classes of the R. thermocellum pan-genome, the COG database was used to classify the functional genes. In the pan-genome of R. thermocellum, the most core, dispensable, and specific gene clusters fell in the metabolism category (Figure 3). Meanwhile, the results indicated that both gene clusters for general function prediction and unknown functions were also more abundant. Compared with the dispensable and specific gene clusters, the majority of genes in the core gene clusters were involved in translation, ribosomal structure, biogenesis (J), cell wall/membrane/envelope biogenesis (M), carbohydrate transport and metabolism (G), amino acid transport, and metabolism (E). By comparison, the majority of genes in the no-core gene clusters was concerned with housekeeping functions, for example, replication, recombination and repair (L), and defense mechanisms (V). Single nucleotide polymorphism (SNP) represents the variation situation of bacteria in the evolutionary process. We constructed the phylogenetic tree based on SNP, which revealed the similarity in bacterial strain variation in adaptation to the natural environment. Compared with R. thermocellum ATCC 27405, M3 showed stronger evolutionary relationships with the other three R. thermocellum strains (Figure 4). Orthologous protein linear analysis based on MCScanX software was performed to further understand the differences in protein homology and amino acid arrangement between M3 and the other four R. thermocellum strains (Figure 5). Between M3 and other three R. thermocellum strains (DSM 1313, DSM 2360 and AD2), it was found that the proteins are not only homologous but also have a good linear relationship in sequence, while R. thermocellum ATCC 27405 was related to M3 but showed a large number of inversions. At the same time, we performed a genome-linear analysis of M3 and R. thermocellum ATCC 27405 (Supplementary Figure S5). The results were similar to the results of the orthologous protein analysis, and it was more clearly evident that nearly half of the genes were inversions.

The results of the CAZyme system analysis indicated that the multi-modular enzyme system of R. thermocellum M3 consisted of 75 dockerin and eight cohesin, which construct scaffolding structural proteins offering a large number of binding sites for cellulase. The quantity of enzyme protein and composition across dissimilar CAZy families in R. thermocellum M3 were analyzed and compared to those in the other four R. thermocellum to evaluate the inclination for lignocellulose saccharification. In the genome of R. thermocellum, most genes fell into glycoside hydrolases (GHs), carbohydrate binding molecules (CBMs), and glycosyltransferases (GTs), whereas a few genes were annotated in carbohydrate esterases (CEs), polysaccharide lyases (PLs), and auxiliary activities (AAs). To be specific, 271 cellulase proteins were detected in R. thermocellum M3 using dbCAN2 (DIAMOND algorithm; Table 3). Most proteins were detected to be GHs (108 candidates), with GH124 (n = 71), GH9 (n = 17), and GH5 (n = 13) being the most abundant families. In addition, there were 75 carbohydrate-binding module (CBM) proteins, 55 glycosyl transferases (GTs) proteins, 22 carbohydrate esterases (CEs) proteins, seven polysaccharide lyases (PLs) proteins, and four auxiliary activities (AAs) proteins. It is worth noting that there were nine CBM genes and 24 GH genes that had not been identified in R. thermocellum before (Supplementary Table S4).

Four encoded proteins associated with cellobiose degradation were found among the cellulase-encoded proteins of strain M3 by Conserved Domain Database (CDD) analysis (Supplementary Figure S6), including BglA (Protein: PROKKA_02182), BglB (Protein: PROKKA_00862), and BglX (Protein: PROKKA_01050 and PROKKA_02059). BglA encoded glycoside hydrolase family 1 (GH1) protein (β-glucosidase). BglB encoded glycoside hydrolase family 1 (GH1) protein (6-phospho-β-glucosidase). BglX (PROKKA_01050) encoded glycoside hydrolase family 3 (GH3) protein (β-glucosidase). BglX (PROKKA_02059) encoded glycoside hydrolase family 3 (GH3) protein (β-glucosidase). The result of protein sequence alignment indicated that the ratio of similarity of pairwise comparison of proteins encoded by three genes (PROKKA_02182, PROKKA_01050, and PROKKA_02059), which correlated with β-glucosidase, were lower than 28.39%. To further demonstrate the diversity among the three Bgls, we constructed phylogenetic trees of Bgls from different strains by MEGA software. Not surprisingly, the phylogenetic tree based on β-glucosidase showed that each of the three Bgls was located in one of three different clusters (Figure 6). Furthermore, a gene encoding β-galactosidase (PROKKA_03012) was present in the M3 genome, which was not found in the other wild R. thermocellum before.

In addition, 11 coding proteins about ABC sugar transport protein and two coding proteins (PROKKA_01086 and PROKKA_02112) about cellobiose phosphorylase coming from glycoside hydrolase family 94 (GH94) were detected in the genome of R. thermocellum M3 (Supplementary Table S5). Two group coding proteins (MglA and UpgA) relating polysaccharide transportation (cellobiose, fructose, and arabinose) ware found in 11 ABC sugar transport proteins by CDD. MglA (PROKKA_00080 and PROKKA_01980) is the ATPase component of ABC sugar transporter, which mainly promotes sugar transport; the primary function of UpgA (PROKKA_01330 and PROKKA_02933) is the transmembrane of sugar substrates, which is a permease component of the ABC transporter. Meanwhile, two unique ABC sugar transport proteins (PROKKA_02931 and PROKKA_02932), which had not been reported in other R. thermocellum, were found in strain M3.

Ruminiclostridium thermocellum, as a typical genus of thermophagic cellulosic-degrading bacteria, was a high-utility candidate in lignocellulosic biomass refinement by means of a powerful cellulosome system (Akinosho et al., 2014; Sheng et al., 2016). However, the activity of β-glucosidase (BglA) was low in wild-type R. thermocellum strains (Shinoda et al., 2019), and due to the adsorption of cellulosome to cellulose, only a small fraction of BglA was available on cellulosome. In nature, the efficient decomposition of cellulosic biomass strongly depends on the fast adaptation of R. thermocellum to enable the regulation of its cellulosomal enzymes for a specific substrate composition. Compared to other R. thermocellum, such as genus R. thermocellum ATCC 27045 and AD 2, genus R. thermocellum M3 not only has high cellulosic saccharification ability but also possesses the properties of cellubiose-resistance, which contributes to M3 being a potential industrial option for cellulosic biofuels refine.

In general, the wild type R. thermocellum containing β-glucosidase A (BglA) and β-glucosidase B (BglB), gene BglX (β-glucosidase) was first found in R. thermocellum in this study. In general, β-glucosidase A was sensitive to temperature and thus inactivated under the optimum temperature of R. thermocellum (Yoav et al., 2019). Different from β-glucosidase A, BglB encodes a novel thermostable β-glucosidase B, which is more thermally stable than β-glucosidase A, whereas the biosynthesis of β-glucosidase B is repressed by cellobiose. Moreover, the β-glucosidase B of R. thermocellum was inhibited by low glucose concentration (Waeonukul et al., 2012). Different from BglA and BglB, the gene BglX (Protein:PROKKA_01050) may encode β-xylosidase or β-glucosidase, since R. thermocellum cannot utilize xylose as a carbon source as the previous report (Sheng et al., 2016) suggests that the BglX (Protein: PROKKA_01050) gene mainly plays a role in the hydrolysis of cellobiose in strain M3 and that the BglX (PROKKA_01050) gene encodes an enzyme with β-glucosidase activity predominantly. It is reported that the β-glucosidase encoded by BglX was a periplasmic cellulase that hydrolysis cellobiose into glucose (Souto et al., 2021). In addition, a signal peptide that anchors β-glucosidase BglX to the periplasmatic space was found in the conserved domain of BglX, which might help to reduce the cellobiose concentration in a restricted area of the cell surface (Souto et al., 2021). Kim also found that the β-glucosidase of Aspergillus aculeatus can bind to the yeast surface without any modification (Kim et al., 2013), which suggests that the β-glucosidase of strain M3 could be directly connected to the cell membrane to hydrolyze cellobiose without the involvement of the CBM module.

It is believed that the genome of R. thermocellum does not encode any β-galactosidase (GH42; Ravachol et al., 2016). Surprisingly, GH42 β-galactosidase (PROKKA_03012) was found in the genome of R. thermocellum M3. Some studies found that in addition to β-glucosidase, β-1,4 glucosidic can also be cleaved by β-galactosidase (Morita et al., 2008; Yang et al., 2018). Therefore, GH42 might be another key factor in cellobiose hydrolysis by R. thermocellum M3. Similar to GH42, GH116 was not found in R. thermocellum in the previous study; the GH116 family was reported as a thermophilic β-d-glucosidase, which was found in animals, plants, archaea, and bacteria (Sansenya et al., 2015; Rohman et al., 2019). The GH116 protein from thermophilic bacterial was reported for high hydrolytic activity toward β-1,3- and β-1,4-linked gluco-oligosaccharides and 4-nitrophenyl β-D-glucopyranoside (4NPGlc) artificial substrate (Sansenya et al., 2015); therefore, GH116 protein might be another key factor involved in the hydrolysis of cellobiose by R. thermocellum M3.

The combination of ABC transport protein and cellobiose phosphorylase is a common strategy concerning cellobiose transportation in R. thermocellum (Parisutham et al., 2017). Five putative cellodextrin-specific ABC transporters, labeled as CbpA-D and Lbp, had been identified in R. thermocellum DSM 1313 (Nataf et al., 2009); Yan et al. found that only CbpB plays a key role in cellobiose transport using the functional verification of CbpA-D by genetic inactivation (Yan et al., 2022). We identified four ABC sugar transporters that specifically transport polysaccharides by gene annotation and conserved domain database analysis, but the high affinity for binding to cellobiose needs further proof in future works.

It is believed that the cellulase system of R. thermocellum-cellulosomes contains cellulose-binding modules (CBMs), which leads to catalytic activity varying greatly in different regions. As a result, there are higher local cellobiose concentrations at particular sites. However, the β-glucosidase of wild type R. thermocellum was not bound to the CBM which led to the ineffective hydrolysis of local cellobiose (Yoav et al., 2019). Therefore, the construction of an artificial chimeric cohesion-containing scaffold in which binding the β-glucosidase to the cellulosome and mimicking the enzymatic synergism of native cellulosome systems was feasible to enhance the hydrolysis efficiency of cellobiose (Fierobe et al., 2002, 2005; Morais et al., 2010). Different from other R. thermocellum strains, nine CBM family genes were unique to strain M3, including CBM23, CBM36, CBM37, CBM40, CBM47, CBM53, CBM61, CBM70, and CBM75. It is worth noting that the CBM37s were initially discovered in R. albus on the basis of adhesion-defective strains that lacked specific surface proteins. The CBM37s were collectively shown to exhibit a broad specificity pattern, which indicated a mechanism for binding the parent enzymes to cellulosic substrates (Himmel et al., 2010). It is reported that CBM37 was responsible for anchoring substrate with enzymes to the cell surface. For example, CBM37 was reported as the mode that binds cellobiohydrolase and endoglucanase to the bacterial cell surface by the C terminus of glycoside hydrolases (Ezer et al., 2008). In addition, binding the substrate with enzymes (Himmel et al., 2010) might enhance the linkage between the cellulosome and enzymes related to the cellobiose degradation of R. thermocellum M3.

Meanwhile, CBM13 and CBM15 were also found in R. thermocellum M3. CBM13s acquire a larger variety of carbohydrate binding specificities including endo-β-glucanase (EC 3.2.1.6; Ferrer et al., 1996), α-galactosidase (EC 3.2.1.22; Holan et al., 1993), and some other glycoside hydrolases (White et al., 1995). Different from CBM13, CBM35 was reported to have conserved ligand specificity, which is often appended to plant cell wall-degrading enzymes (Montanier et al., 2009) and xylan-degrading enzymes (Chen et al., 1995), but it is often found in β-galactosidase (EC 3.2.1.23; Ali-Ahmad et al., 2017). Therefore, it is suggested that the β-galactosidase of strain M3 is mainly related to CBM13 and CBM35 modules during the degradation of cellulosic feedstocks.

AA6 and AA8 families were first identified in R. thermocellum in this study. In general, members of the AA8 family (cellobiose dehydrogenase, CDH) can be isolated or appended to a CBM. Proteins contain iron reductase domains and may generate reactive oxygen species that could contribute to the non-enzymatic degradation of cellulose chains by the generation of highly reactive hydroxyl radicals (OH•) via Fenton’s reaction (Levasseur et al., 2013). It is reported that the poor cellobiose availability of the substrate was a limiting factor to CDH activity; on the surface of R. thermocellum, cellobiose was accumulated in a restricted area, which was a sufficient substrate for the CDH. Good contact with CDH might play an important role in reducing the concentration of cellobiose in the environment (Kittl et al., 2012) and potentially releasing the strain M3 from the feedback inhibition of cellobiose.