Differential gene expression analysis and functional enrichment of DEGs

To understand the molecular changes in N. hansenii SI1 associated with growth on glycerol as the main carbon and energy source, we performed comparative transcriptomic analysis. Gene expression analysis identified 424 significantly differentially expressed genes (13% of total), with 219 downregulated and 205 upregulated (Fig. 1a). Initial screening revealed few annotated genes in this list; subsequent RAST annotation-based analysis provided deeper insights into functional assignments.

Fig. 1

Analysis of gene expression and functional characterisation of DEGs between the glycerol and the glucose cultures of N. hansenii SI1. a Volcano plot representing gene expression changes. Shown are − log10 of adjusted p-values versus log2 fold changes (FC). The repressed and the induced DEGs are colored in blue and orange, respectively. Gene symbols are plotted only for the DEGs. The dashed horizontal line indicates the adjusted p-value threshold of 0.01. The dashed vertical lines indicate the |log2FC|= 1 thresholds. b, c Results of functional enrichment analysis of DEGs based on RAST annotation. Displayed are enriched RAST subcategories. The enriched functional bins of the downregulated and the upregulated DEGs are displayed as blue and orange bars, respectively. Shown are the results of the one-sided Fisher’s exact test (one-tailed). For clarity, only the bins that scored the adjusted p-value below 0.5 are presented. The red dashed vertical line indicates the adjusted p-value threshold of 0.01

The results of the gene set enrichment analysis yielded insights into functions overrepresented among the DEGs. Only three RAST functional categories demonstrated significant enrichment across both the downregulated and upregulated genes. These categories, ordered by ascending adjusted p-value, were “respiration” and “phages, prophages, transposable elements, plasmids” for the downregulated DEGs (Online Resource 1: Fig. S1a) and “sulphur metabolism” in the case of the upregulated DEGs (Online Resource 1: Fig. S1b). More detailed functions were exposed at the subcategory level (Fig. 1b, c). The significantly overrepresented RAST subcategories among the repressed DEGs included “electron donating reactions” and “phages, prophages” (Fig. 1b). For the induced DEGs, these were “riboflavin, FMN, FAD”, “inorganic sulphur assimilation”, “protein degradation”, “capsular and extracellular polysaccharides”, and “lysine, threonine, methionine, and cysteine” (Fig. 1c). Finally, at the subsystem level, several functions were enriched, and among them, significantly enriched subsystems included “NADH ubiquinone oxidoreductase”, “respiratory complex I”, and “cobalt–zinc–cadmium resistance” for the downregulated DEGs. For the upregulated DEGs, these were “inorganic sulphur assimilation”, “cysteine biosynthesis”, “riboflavin, FMN and FAD metabolism”, “riboflavin, FMN and FAD metabolism in plants”, “riboflavin to FAD”, “proteasome bacterial”, “dTDP-rhamnose synthesis”, and “rhamnose containing glycans”. These results guided further, detailed analysis of the enriched functional subcategories.

Gene expression changes in the central carbohydrate metabolism pathways

During growth on glycerol, oxidation of glucose in the periplasm, along with its uptake and phosphorylation in the cytosol, exhibited mild repression (Fig. 2). The Embden–Meyerhof–Parnas pathway was downregulated. Conversely, both the Entner–Doudoroff (ED) pathway and the pentose phosphate pathway (PPP) were induced, with significant changes observed in the gene encoding 6-phosphogluconate dehydrogenase (gndA, log2FC = 1.0, adj. p-value = 5.22 * 10−18; Fig. 2 and Online Resource 1: Fig. S2).

Fig. 2

Transcriptomic changes in the central carbohydrate metabolism pathway between N. hansenii SI1 cultures grown in the glycerol and the glucose medium. a The pathway of glycerol metabolism and glycolysis/gluconeogenesis. Genes are colored according to log2 fold change in expression between the glycerol and the glucose cultures. IM, inner membrane; OM, outer membrane. b Expression of the glycerol operon (glp). c Glycerol dehydrogenase subunits of two loci. Transcripts mean FPKM values are shown in either gray or green for cells grown in the glucose or the glycerol medium, respectively. Bars represent the means from 3 replicated cultures. Thin black bars denote standard error. Stars in all panels denote statistically significant changes (called by DESeq2; adjusted p-value < 0.01 and |log2FC|  ≥ 1)

The pathway for glycerol uptake and catabolism was upregulated (Fig. 2a), with strong induction observed in genes within the glp cluster (Fig. 2b). The genes on the pathway from DHAP to glucose-6-phosphate, such as tpiA, fba, tal, and glpX, were upregulated, where, for the last one, this change was significant (log2FC = 1.0, adj. p-value = 1.23 * 10−8; Fig. 2a).

The genes encoding copies of large and small subunits of glycerol dehydrogenase are located in two loci in the genome of N. hansenii SI1 (sldAB1, sldAB2). The sldAB2, of which expression is the highest of the two gene sets, was mildly downregulated in the glycerol medium (Fig. 2a and 2c). The expression of sldAB1, which is lower than that of sldAB2, is significantly induced in the glycerol medium (log2FC = 1.9 for both subunits, adj. p-value = 1.23 * 10−8 (sldA1) and 4.62 * 10−8 (sldB1)).

The sldAB1 genes are located in close proximity to a gene cluster which encloses predicted tagatose kinase, tagatose 6-phosphate 4-epimerase, sorbitol dehydrogenase, and putative metabolite transport protein (Online Resource 1: Fig. S3a). These genes are all significantly upregulated (Online Resource 1: Fig. S3b).

The metabolism of DHA produced by glycerol dehydrogenase is not completely clear. Based on the functional predictions, the genome of N. hansenii SI1 does not encode dihydroxyacetone kinase, as is the case with the strains of the Komagataeibacter genus and Novacetimonas maltaceti LMG 1529 (based on UniProt database search, data not shown). One possibility is that phosphorylation of DHA is conducted by glycerol kinase (glpK), since this activity was found for its ortholog in K. xylinus (Acetobacter xylinum in the original work (Weinhouse and Benziman 1976)). This gene appears to be the sole candidate for this function in the N. hansenii SI1 genome.

Gene expression changes in the TCA cycle, acetoin metabolism, and respiration

The tricarboxylic acid (TCA) cycle was overall downregulated (Fig. 3a), which agreed with the results of RAST enrichment analysis (Fig. 1b). Many genes of this pathway were significantly repressed. These were such genes as maeA1 (log2FC =  − 1.7, adj. p-value = 1.44 * 10−43), maeA3 (log2FC =  − 1.9, adj. p-value = 9.15 * 10−23), acnA (log2FC =  − 1.0, adj. p-value = 3 * 10−86), icd (log2FC =  − 1.1, adj. p-value = 6.48 * 10−37), hicd (log2FC =  − 1.0, adj. p-value = 2.4 * 10−28), aarC (log2FC =  − 1.22, adj. p-value = 5.72 * 10−36).

Fig. 3

Changes in the expression of genes involved in the TCA cycle and respiration in N. hansenii SI1. a Transcriptomic changes in the TCA cycle. Genes are colored according to log2 fold change in expression between glycerol and glucose cultures. Stars denote statistically significant changes (called by DESeq2; adjusted p-value < 0.01 and |log2FC| ≥ 1). IM, inner membrane; OM, outer membrane. b Mean expression levels of the subunits of cytochrome ba(3) ubiquinol oxidase (UOX). c Mean expression levels of the subunits of NADH-quinone oxidoreductase operon (nuo, complex I). Transcripts’ mean FPKM values are shown either gray or green for cells grown in either the glucose or glycerol medium, respectively. Bars represent the means from 3 replicated cultures. Thin black bars denote standard error. Stars in all panels denote statistically significant changes (called by DESeq2; adjusted p-value < 0.01 and |log2FC|  ≥ 1)

Other genes associated with the TCA cycle, such as those encoding L-lactate dehydrogenase and alcohol dehydrogenase subunits, were significantly repressed (lutA (log2FC =  − 1.4, adj. p-value = 1.44 * 10−43), lutB (log2FC =  − 1.2, adj. p-value = 8.55 * 10−41), adhA1 (log2FC =  − 1.4, adj. p-value = 3.02 * 10−80), adhB1 (log2FC =  − 1.4, adj. p-value = 2.21 * 10−79), adhS (log2FC =  − 1.0, adj. p-value = 4.44 * 10−13); Fig. 3a). In contrast, the genes encoding pyruvate phosphate dikinase, the subunits of pyruvate dehydrogenase, and aldehyde dehydrogenase were induced, where for ppdk this change was significant ((log2FC = 1.1, adj. p-value = 1.17*10−27); Fig. 3a).

Additionally, we focused on another pathway associated with the TCA cycle and connected to acetoin metabolism. The two genes of this operon, budB and alsD, are significantly repressed in the glycerol medium (log2FC =  − 3.1, adj. p-value = 9.82*10−231 and log2FC =  − 2.6, adj. p-value = 1.53*10−91, respectively). The products of this operon, acetolactate synthase (budB) and alpha-acetolactate decarboxylase (alsD), are involved in acetoin biosynthesis. Our transcriptomic results show that the putative pathway of acetoin metabolism is generally downregulated when the cells are grown on the glycerol medium (Online Resource 1: Fig. S4).

Finally, we investigated the response of genes associated with the respiratory apparatus. Oxidation of substrates in the periplasm by membrane-bound dehydrogenases is coupled to the reduction of ubiquinone to ubiquinol. Our results indicate mild downregulation of the genes encoding the subunits of the ba3-type ubiquinol oxidase (Fig. 3b). Cytoplasmic enzymes assimilate these oxidation products and further metabolise substrates to generate NADH. The expression of the nuo operon of the NADH dehydrogenase (complex I) was strongly repressed, with most genes showing significant expression changes (Fig. 3c).

The results obtained at this stage have shown that, when the N. hansenii SI1 cells are grown on the glycerol medium, TCA cycle progression, organic acid production, and cytoplasmic (NADH) respiration are repressed.

Transcriptional activation of riboflavin biosynthesis pathway

In agreement with the results of RAST enrichment analysis, we found that the riboflavin biosynthesis pathway was upregulated, with the majority of genes in this pathway showing significant expression changes (Fig. 4a). The highest expression changes were displayed by the rib genes (log2FC > 1 and adj. p-value < 0.01 for most of these genes), of which ribD, ribC, ribBA, and ribE likely form a transcriptional unit (Fig. 4b). A predicted FMN riboswitch is located upstream of the rib cluster, which may regulate FMN/FAD biosynthesis (Fig. 4b).

Fig. 4

The predicted riboflavin biosynthesis pathway in N. hansenii SI1. a Transcriptomic changes in the riboflavin biosynthesis pathway (based on Solopova et al. 2020). Genes are colored according to log2 fold change in expression between glycerol and glucose cultures. Stars denote statistically significant changes (called by DESeq2; adjusted p-value < 0.01 and |log2FC| ≥ 1). b Gene loci encoding the predicted FMN riboswitch and the riboflavin synthesis enzymes in the genome of N. hansenii SI1. The gene identifiers and symbols correspond to the following predicted functions: ribD—bifunctional diaminohydroxyphosphoribosylaminopyrimidine deaminase (EC 3.5.4.26)/5-amino-6-(5-phosphoribosylamino)uracil reductase (EC 1.1.1.193); ribH—6,7-dimethyl-8-ribityllumazine synthase (EC 2.5.1.78); ribE—riboflavin synthase (EC 2.5.1.9); ribBA—bifunctional GTP cyclohydrolase II (EC 3.5.4.25); 3,4-dihydroxy-2-butanone 4-phosphate synthase (EC 4.1.99.12); ribCF—bifunctional riboflavin kinase (EC 2.7.1.26)/FAD synthase (EC 2.7.7.2); rutF—FMN reductase; SI1_01887—5,6-dimethylbenzimidazole synthase

The ribCF gene is located in a different genome location than the rib cluster and forms a putative cluster with genes (mntC, mntD, and mntB) that are predicted to be involved in the methionine salvage pathway (Online Resource 1: Fig. S5a). These genes, similarly to ribCF, are significantly upregulated in the glycerol medium (Online Resource 1: Fig. S5b). Based on the genome annotation, other genes of this pathway are mildly upregulated (mtnK and mtnA) or downregulated (mtnP and mtnN).

Expression of sulphur assimilation and cysteine biosynthesis pathway

One of the significantly enriched RAST functional subcategories was “inorganic sulphur assimilation” (Fig. 1c). Upon examining the sulphur assimilation pathway, we observed its general activation in the culture grown on the glycerol medium (Fig. 5). Specifically, genes involved in sulphate transport and its reduction to sulphide were predominantly upregulated.

Fig. 5

Transcriptomic changes in the predicted pathway of sulphur assimilation and cysteine biosynthesis in N. hansenii SI1. Genes are colored according to log2 fold change in expression between glycerol and glucose cultures. Stars denote statistically significant changes (called by DESeq2; adjusted p-value < 0.01 and |log2FC| ≥ 1)

The genes associated with alkanesulfonate transport and metabolism were mildly induced, while the downstream portion of the cysteine biosynthesis pathway did not exhibit significant changes in expression (Fig. 5).

Protein degradation and stress response

The results of functional enrichment analysis revealed particular induction of the genes involved in protein degradation (Fig. 1c). Indeed, we observed that the genes encoding ATP-dependent proteases were mostly upregulated, and some of them (clpP1, clpP2, clpX1, lon) significantly (Fig. 6a). A comparable pattern was observed for the majority of genes encoding molecular chaperones (Fig. 6b). Here, such genes as groEL, groES, ibpA, and hrcA were significantly induced.

Fig. 6

Expression changes in genes associated with protein degradation and stress response. a ATP-dependent proteases. b Molecular chaperones. c DNA damage and the general stress-responsive genes. d Transcription regulators associated with oxidative stress response. e Genes coding for oxidative stress-responsive enzymes. f Glutaredoxins (grx) and thioredoxin A (trxA) genes. g Genes encoding methionine sulphoxide reductases (msr). h Suf cluster of genes involved in Fe–S cluster repair. Transcript mean FPKM values are shown in gray or green for cells grown in either the glucose or the glycerol medium, respectively. Bars represent the means from 3 replicated cultures. Thin black bars denote standard error. Stars in all panels denote statistically significant changes (called by DESeq2; adjusted p-value < 0.01 and |log2FC|  ≥ 1)

Further analysis revealed mild upregulation of genes linked to DNA damage and the general stress response (Fig. 6c).

The expression of homologs of transcription factors associated with the oxidative stress response was mostly induced in the culture grown on the glycerol medium (Fig. 6d). For iscR2 and oxyR1, this change was significant (log2FC = 1.25, adj. p-value = 1.29 * 10−13 and log2FC = 1.28, adj. p-value = 2.59 * 10−28, respectively). More severe upregulation was displayed by the genes encoding enzymes involved in the oxidative stress response (Fig. 6e). Here, such genes as ahpC, ahpF, hmp, and putA were significantly induced, where for hmp and putA the log fold change was greater than 2 (log2FC = 2.15, adj. p-value = 4.70 * 10−39 for hmp and log2FC = 2.19, adj. p-value = 1.64 * 10−43 for putA).

Beyond antioxidant enzymes, genes encoding repair proteins for oxidised molecules were induced. One such group were the genes coding for three predicted glutaredoxins (grx) and thioredoxin A (trxA), where for the trxA gene the change in expression was significant (log2FC = 1.26, adj. p-value = 8.09*10−48; (Fig. 6f)).

Methionine sulphoxide reductases (msr), which reduce free and protein-based methionine sulphoxides to methionine, were also upregulated. Among the five predicted msr genes in the N. hansenii SI1 genome (msrA, msrB, msrC, msrP, and msrQ), only msrB exhibited significant upregulation (log2FC = 1.30, adj. p-value = 1.01 * 10−24; Fig. 6g).

Additionally, the suf cluster genes, involved in Fe–S cluster repair, displayed mild induction (Fig. 6h). Among its eight genes, only iscA, which is located upstream of the suf cluster, showed significant changes (log2FC = 1.18, adj. p-value = 2.79*10−45).

Metal homeostasis

Guided by the results of functional enrichment of RAST subcategories, we explored first the “cobalt-zinc-cadmium resistance” subsystems, both of which were significantly enriched among the downregulated DEGs (Fig. 1b and Online Resource 1: Fig. S1c).

Based on the N. hansenii SI1 genome annotation, two gene clusters (czcCBA1 and czcCBA2) encoding homologs of the cobalt-zinc-cadmium resistance system CzcCBA of Cupriavidus metallidurans CH34 (Janssen et al. 2010) were identified. These clusters exhibit high sequence similarity to the previously described cusCBASR cluster in K. xylinus E25 (Ryngajłło et al. 2019a). The czc2 cluster, similarly as in K. xylinus E25, is followed by two genes predicted to encode the CzcS/CzcR two-component system regulating the CzcCBA efflux system in Pseudomonas aeruginosa (Perron et al. 2004). Genes in both clusters are significantly downregulated during growth on the glycerol medium, except for czcS and czcR, which are only mildly repressed (Fig. 7a).

Fig. 7

Expression changes in genes encoding proteins involved in metal homeostasis. a Two gene clusters (czc1 and czc2) encoding homologs of the cobalt-zinc-cadmium resistance system. b Genes encoding homologs of SmtA and SmtB proteins involved in zinc and other essential ions homeostasis. c Copper resistance genes (copA1, copB, SI1_00647 [unknown function], and copA2) and a predicted merP gene encoding a mercury scavenger. d Genes encoding the most highly expressed TonB-dependent receptors. e Genes coding for homologs of the colicin I receptor, CirA. Transcript mean FPKM values are shown in gray or green for cells grown in either the glucose or glycerol medium, respectively. Bars represent the means from 3 replicated cultures. Thin black bars denote standard error. Stars in all panels denote statistically significant changes (called by DESeq2; adjusted p-value < 0.01 and |log2FC| ≥ 1)

Regarding zinc homeostasis, we found in the N. hansenii SI1 genome two genes in consecutive arrangement, which encode homologs of SmtA and SmtB of Synechococcus elongatus PCC 7942. SmtA functions as a metallothionein, while SmtB acts as a transcriptional repressor of smtA. These proteins are critical for maintaining zinc and other essential ions homeostasis (Huckle et al. 1993). In Synechococcus elongatus PCC 7942, SmtA sequesters and detoxifies four Zn2+ ions per molecule (Blindauer et al. 2002), while SmtB represses smtA expression under low Zn conditions (Huckle et al. 1993). In our analysis, we found that smtA and smtB were highly expressed in the glycerol medium and are both ranked as top differentially expressed genes (log2FC = 3.61, adj. p-value = 0 for smtA and log2FC = 3.41, adj. p-value = 0 for smtB, Fig. 7b).

Significant induction in expression was observed for the homologs of copper resistance proteins CopA and CopB from Pseudomonas syringae pv. tomato (Fig. 7c). In Pseudomonas species, the cop systems consist of four structural proteins (CopABCD), whose main role is the expulsion of copper from the cell (Hofmann et al. 2021). While N. hansenii SI1 lacks homologs of CopC/CopD, copA1, copB, and an uncharacterised gene (SI1_00647) form a cluster. Phobius predicts SI1_00647 encodes a periplasmic protein, potentially analogous to CopC, though functional validation is needed.

Apart from CopA and CopB proteins, the N. hansenii SI1 genome encodes a copper-exporting P-type ATPase homologous to Staphylococcus aureus NCTC 8325 CopA (copA2), which was significantly upregulated in the glycerol medium (Fig. 7c). Adjacent to copA2, we identified a gene encoding a MerP-like mercury scavenger from Shigella flexneri, showing mild induction (Fig. 7c).

For cobalt, we observed mild upregulation of genes encoding homologs of cobalt transporter subunits (Online Resource 1: Fig. S6a).

Previously reported TonB-dependent transporters (TBDTs (Ryngajłło et al. 2024)) were predominantly downregulated in the glycerol medium (Fig. 7d and Online Resource 1: Fig. S6b). Notable exceptions included the SI1_00389 gene (significant induction), annotated as colicin I receptor (cirA homolog of E. coli), and a second cirA copy (cirA2) with slight upregulation (Fig. 7e).

Phages and prophages

The results of RAST functional analysis highlighted the “phages, prophages” subcategory as significantly enriched among the downregulated DEGs (Fig. 1b). A total of six genes were assigned to this subcategory, encoding phage packaging machinery, neck, and capsid proteins. Despite being expressed at low levels, these genes exhibited significant downregulation in the glycerol medium cultures (Online Resource 1: Fig. S7).

We further compared these findings with the predictions from the PHASTEST program, which identified two prophage regions in the N. hansenii SI1 chromosome. The first region is intact and spans positions 1,089,369–1,117,668 bases, while the second region (designated as questionable) occupies positions 1,990,601–2,010,067 bases. Notably, the genes of the RAST subcategory overlap both prophage loci. By analysing transcription at single-nucleotide resolution, we observed that both regions exhibited predominantly repressed expression in the glycerol medium, with most genes showing statistically significant changes (Online Resource 1: Fig. S8 and Online Resource 2: Table S1).

Expression of cellulose and acetan synthesis pathways

The cellulose synthesis pathway was generally downregulated in the glycerol medium cultures, with the bcsC gene showing significant expression changes (log2FC =  − 1.1, adj. p-value = 1.76 * 10−22; Fig. 8a, b). Similarly, the genes from the other two cellulose synthase (CS) operons (bcsII and bcsIII) were downregulated (Online Resource 1: Fig. S9a, b).

Fig. 8

Gene expression changes in the cellulose and acetan-like II biosynthesis pathway in N. hansenii SI1. a The pathway. The genes are colored according to log2 fold change in expression between glycerol and glucose cultures. b Expression of the bcsI cellulose synthase cluster. c Expression of the acetan type II gene cluster. For b–c, transcripts’ mean FPKM values are shown in either gray or green for cells grown in the glucose or glycerol medium, respectively. Bars represent the means from three replicated cultures. Thin black bars denote the standard error. Stars in all panels denote statistically significant changes (called by DESeq2; adjusted p-value < 0.01 and |log2FC|  ≥ 1)

In contrast, strong and statistically significant induction was observed for the majority of genes involved in acetan-like polymer biosynthesis. Key examples include aceA, galE, and other genes of the acetan type II cluster (Fig. 8a, c), as well as chromosomal rml1 genes responsible for dTDP-rhamnose biosynthesis (Fig. 8a and Online Resource 1: Fig. S10c). The genes of the rml2 cluster, which are located on the p1 plasmid, were mildly downregulated (Fig. 8a and Online Resource 1: Fig. S10d).