Skip to main content
Download PDF

New insights into the efficient secretion of foreign protein in Bacillus subtilis via Ribo-seq and RNA-seq integrative analyses

BMC Microbiology volume 24, Article number: 537 (2024) Cite this article

Abstract

Background

As an important prokaryotic model organism, Bacillus subtilis has been widely used in the industrial production of a variety of target products. The efficient secretion of target products has always been the main purpose of industrial microbial technology. The modification of gene regulatory networks is an important technical means to construct a factory of microbial cells that efficiently secretes target products. However, the regulatory network of the efficient expression of foreign genes in B. subtilis has not been studied at the translation level.

Results

In this study, Ribo-seq and RNA-seq technology were used to study the changes in differentially expressed genes during the efficient secretion of the protease PB92 by B. subtilis WB600, and the results revealed the gene regulatory network related to efficient secretion of foreign protein. The results revealed that the correlation between the differentially expressed genes of B. subtilis at the transcription and translation levels was only 0.5354. Forty-one common (transcription and translation) and 436 unique (translation) key differential gene sets that may be related to the efficient secretion of foreign proteins were revealed. KEGG enrichment analysis of these key gene sets revealed that they were involved mainly in the cell motility and central metabolic regulatory network of B. subtilis.

Conclusion

Our study provides important guidance for the construction of cell factories and metabolic networks for the efficient secretion of target products by B. subtilis.

Peer Review reports

Introduction

As a model organism of Bacillus subtilis, the expression system has been widely used in the production of industrial enzymes, the synthesis of metabolites, and basic research [1,2,3]. Compared with other expression systems, it has the advantages of a fast growth rate, strong protein secretion ability, clear genetic background, and food-grade safe microorganisms, which are ideal host bacteria for the expression of foreign genes [1,2,3].

As a microbial cell factory, B. subtilis can be considered a sustainable biochemical reactor [4]. At present, an effective strategy to improve the yield and production capacity of target compounds is to use metabolic engineering to rationally modify and regulate cell factories [5]. Two methods are commonly used: (i) according to the theory of metabolic regulation, fermentation process conditions (pH, temperature, aeration capacity, medium composition, etc.) are optimized to change the metabolic balance in bacteria and maximize the accumulation of target metabolites [6]; and (ii) strain-modification strategies based on gene knockout and overexpression can change the characteristics of microbial gene regulatory networks and fundamentally change the metabolism of microorganisms [7]. In response to changes in fermentation process, the optimal fermentation process conditions for the target metabolite can be determined via univariate optimization [8]. However, for the modification of strains, it is necessary to understand in depth and mine the gene regulatory network of cell factories to further modify the relevant metabolic pathways [8, 9].

Gene regulatory networks dynamically modulate the levels of mRNAs and proteins in response to cellular requirements, exerting influence across various stages of gene expression, including transcription, mRNA processing, transport, translation, and protein degradation [10]. To understand the intricate molecular mechanisms underlying biological processes and phenotypes, researchers have increasingly relied on next-generation sequencing techniques such as RNA-seq over the past two decades [11]. However, while RNA-seq can provide valuable insights into gene transcription, it may not fully capture the complete landscape of gene expression due to the influence of factors such as mRNA stability, degradation, and translation efficiency [11, 12]. In the gene regulatory network, translation regulation is also one of the key components of gene expression regulation, accounting for more than all other regulatory processes combined (> 50%) [13]. To overcome these limitations, ribosome profiling, a cutting-edge method that involves deep sequencing of mRNA regions shielded by ribosomes, has emerged as a powerful tool for investigating translation regulation on a genome-wide scale [14, 15]. By mapping the positions of ribosomes along entire ribosome profiles, Ribo-seq enables researchers to examine translation events more directly than with RNA-seq, providing a more precise and comprehensive understanding of protein expression [16].

Therefore, the purpose of this study was to explore the gene regulatory network in which the efficient expression of foreign gene in B. subtilis is related at both the transcriptional and translational levels. Prior to this, the protease PB92 gene has been shown to be overexpressed in the B. subtilis expression system, exhibiting good extracellular protease activity and an optimal pH of extreme alkalinity, which has aroused great research interest in academia and industry and has been applied in modern industry [17,18,19,20]. In this study, the extracellular enzyme activity of the protease PB92 was taken as the target product, and integrative analysis of the transcriptome (RNA-seq) and translatome (Ribo-seq) revealed the gene regulatory network of B. subtilis when it efficiently expressed foreign protein. We hope our findings will provide a theoretical basis and reference for the efficient expression of foreign genes in B. subtilis and the modification of gene regulatory networks in the future.

Materials and methods

Strain and fermentation

The gene encoding the protease PB92 derived from Bacillus alcalophilus (GenBank accession no. WP_094423791.1) spans a length of 1143 nucleotides and was initially cloned and inserted into the pMD19-T vector, which is currently housed in our laboratory [20]. To construct the recombinant plasmid pBE-09-PB92, we employed the In-Fusion HD Cloning Kit (Takara Inc., Dalian, China), which incorporates the robust pShuttle-09 promoter and BslB signal peptides [21]. Subsequently, electroporation was employed to introduce the recombinant plasmid into B. subtilis WB600 [22]. We screened for a single colony on a plate supplemented with 50 mg/mL kanamycin and 1% skim milk powder and observed the formation of a transparent circle as the criterion. The selected colony was used to inoculate 5 mL of LB broth, which served as the seed culture and was then incubated at 37 °C. For subsequent fermentation, a 500-mL flask containing fresh superrich medium (10 g/L lactose, 60 g/L peptone, 10 g/L NH4Cl, 8 g/L skim milk powder, 5 g/L NaCl, 5 g/L MgSO4·7H2O, and 5 g/L K2HPO4) was inoculated with 5% of the seed culture and further incubated at 37 °C. At 16 predetermined time points ranging from 0 to 96 h after inoculation, samples of the fermentation broth were collected. These samples were subjected to centrifugation at 12,000 ×g for 20 min to obtain the enzyme mixture, and the protease activity was evaluated using azocasein as a substrate at each time point, following the method described by Miao et al. [20, 21].

Sample preparation

Three sets of parallel experiments were performed according to the same method as above. For fermentation for 12 h (P1), 48 h (P2), and 60 h (P3), (3 groups × 3) fermentation broth samples were collected sequentially, and some samples were used for RNA sequencing and protease activity determination. To stop translation, chloramphenicol was added to the remaining sample volume at a final concentration of 100 µg/mL, and the mixture was incubated at 37 °C for 2 min and then centrifuged at 4500 rpm for 10 min at 4 °C, after which the supernatant was discarded [23, 24]. The bacteria were washed once with 10 mL of precooled wash solution (10 mM MgCl2, 100 mM NH4Cl, 20 mM Tris pH 8.0, and 1 mM chloramphenicol) and centrifuged at 3000 rpm at 4 °C for 5 min, after which the supernatant was discarded. The samples were quickly frozen in liquid nitrogen at -80 °C. Nine samples (12 h, 48 h, and 60 h) were RNA-sequenced, and 6 samples were Ribo-sequenced from the two groups (12 h and 60 h) with the greatest difference in extracellular protease activity.

RNA sequencing and data analysis

Total RNA was isolated via the TRIzol-based method (Life Technologies, CA, USA) [25]. Following RNA quality assessment, rRNA depletion was carried out via the Illumina MRZB12424 Ribo-Zero rRNA Removal Kit (Illumina, San Diego, CA, USA). First-strand cDNA synthesis was subsequently performed via ProtoScript II Reverse Transcriptase (New England BioLabs, Ipswich, MA, USA) at 25 °C for 10 min, 42 °C for 15 min, and 70 °C for 15 min. Second-strand cDNA synthesis was then conducted via NEBNext Second Strand Synthesis Reaction Buffer and dNTP mix (New England BioLabs, Ipswich, MA, USA) at 16 °C for 1 h. The resulting cDNA was purified via Agencourt AMPure XP beads (Beckman Coulter, Brea, CA) and subjected to end repair via NEBNext End Repair Reaction Buffer and Enzyme Mix (New England BioLabs, Ipswich, MA, USA) at 20 °C for 30 min followed by 65 °C for 30 min. Sequencing adapters were ligated via NEBNext Adaptor for Illumina (New England BioLabs, Ipswich, MA, USA) at 20 °C for 15 min. Subsequently, the second-strand cDNA was degraded via the USER enzyme mix (New England BioLabs, Ipswich, MA, USA) at 37 °C for 15 min, and the resulting product was purified via Agencourt AMPure XP beads (Beckman Coulter, Brea, CA). In the final step of our experimental workflow, the index-coded samples underwent clustering via the cBot cluster generation system. This process involved the use of NEBNext Q5 Hot Start HiFi PCR Master Mix from New England Biolabs (Ipswich, MA, USA). Once the clusters were successfully generated, we sequenced them on the Illumina NovaSeq 6000 platform, employing the paired-end sequencing strategy with 150 base reads.

The raw sequencing data were subjected to a series of filtering steps based on the following criteria: (1) removal of reads containing ≥ 10% unidentified nucleotides (N); (2) removal of reads with > 50% bases having phred quality scores ≤ 20; and (3) elimination of reads aligned to the barcode adapter via FASTP (version 0.18.0) [26]. The quality-trimmed reads were subsequently mapped to the reference genome of B. subtilis subsp. subtilis str. 168 (https://www.ncbi.nlm.nih.gov/genome/665?genome_assembly_id=300274) via Bowtie2 (version 2.2.8), with zero mismatches allowed. Reads mapped to ribosomal RNA were then excluded from further analysis [27, 28]. The filtered reads were subsequently aligned to the reference genome via Bowtie2 to identify known genes, and gene expression levels were quantified via the RSEM algorithm [29]. To assess reproducibility among samples, the correlation coefficient was calculated via the R package psych (version 1.8.4) [30], with values closer to one indicating greater reproducibility. To account for variations in gene length and sequencing depth, the gene expression levels were further normalized via the fragments per kilobase of transcript per million (FPKM) method [29]. This normalization approach reduces biases caused by differences in gene length and total sequencing reads. Differential expression analysis was carried out via the edgeR package (http://www.r-project.org/), considering genes with fold changes ≥ 2 and a false discovery rate-adjusted P value (FDR) < 0.05 as differentially expressed genes (DEGs) [31].

Ribo sequencing and data analysis

The sample was subjected to bacterial lysis by adding bacterial lysis buffer, which was then mixed and incubated on ice for 10 min. The resulting lysate was centrifuged at 17,000 ×g for 10 min at 4 °C, and the supernatant was collected. To 400 µL of lysate, 10 µL of RNase I (NEB, Ipswich, MA, USA) and 6 µL of DNase I (NEB, Ipswich, MA, USA) were added, followed by incubation for 45 min at room temperature with gentle mixing on a Nutator mixer. Nuclease digestion was halted by adding 10 µL of SUPERase·In RNase inhibitor (Ambion, Austin, TX, USA). Size exclusion columns (Illustra MicroSpin S-400 h Columns; GE Healthcare; catalogue no. 27-5140-01) were prepared by equilibrating them with 3 mL of polysome buffer through gravity flow, followed by centrifugation at 600 ×g for 4 min at room temperature. Subsequently, 100 µL of digested ribosome footprints (RFs) were loaded onto the column and centrifuged at 600 ×g for 2 min. Elution was performed by adding 10 µL of 10% (wt/vol) SDS, and RFs larger than 17 nt were isolated via the RNA Clean and Concentrator-25 kit (Zymo Research; R1017). To remove residual rRNA sequences, short (50–80 bases) antisense DNA probes complementary to rRNA sequences were added to the solution containing RFs. Subsequently, RNase H (NEB, Ipswich, MA, USA) and DNase I (NEB, Ipswich, MA, USA) were introduced to digest the rRNA and residual DNA probes. The purified RFs were further processed via magnetic beads (Vazyme, Nanjing, Jiangsu, China). With the obtained RFs, Ribo-seq libraries were constructed via the NEBNext® Multiple Small RNA Library Prep Set for Illumina® (catalogue no. E7300S, E7300L). This involved adding adapters to both ends of the RFs, followed by reverse transcription and PCR amplification. PCR products ranging from 140 to 160 bp in size were selectively enriched to generate a cDNA library, which was subsequently subjected to sequencing via the Illumina HiSeqTM X10 platform by Gene Denovo Biotechnology Co. (Guangzhou, China).

The FASTP tool [26] was employed to filter out low-quality reads. Specifically, raw reads were scrutinized, and any reads that contained more than 50% low-quality bases or more than 10% N bases were discarded. Additionally, adapter sequences were trimmed from the reads. For alignment to specific databases, such as the ribosomal RNA database, GenBank, and Rfam database [27], the Bowtie2 short-read alignment tool was utilized. Reads that mapped to rRNAs, transfer RNAs, small nuclear RNAs, small nucleolar RNAs, and miRNAs were subsequently eliminated from further analysis. Through read length distribution statistics, we retained only reads with a length of 20–40 bp, that is, RFs, and we used those RFs that met the expected length for subsequent analysis. RSEM software was used to calculate the number of reads at the Ribo-seq level in the open reading frame (ORF) region encoding the gene and convert it to the FPKM value to obtain the expression of the gene at the translation level [29]. RiboDiff was used to identify differential translational efficiency genes (FPKM in Ribo-seq/FPKM in RNA-seq) across sample groups [32], and genes with a fold change ≥ 2 and FDR < 0.05 in a comparison were considered significant differential translation efficiency genes. The enrichment analysis of the Pearson correlation coefficient and differentially expressed genes at the translation level was performed in the same way as for the transcriptome.

Combination analysis

The transcription and translation levels of the differentially expressed genes (fold changes ≥ 2, FDR < 0.05) were plotted with a scatter diagram, and the Pearson correlation coefficient was calculated. The VennDiagram package in R was used for Venn plotting. Target gene sets were mapped to the KEGG database (https://www.genome.jp/kegg/) [33], gene numbers were calculated for every pathway, and significantly enriched pathways of the target gene sets compared with the genome background were defined by hypergeometric testing. With FDR < 0.05 as a threshold, pathways meeting this condition were defined as significantly enriched pathways in DTGs.

Quantitative real-time PCR validation and protein activity assays

In the P1 and P3 comparison groups, three replicate samples were selected from each group for qRT‒PCR experiments and protein activity assays. Total RNA was isolated via the TriQuick Reagent Kit (Solarbio, Beijing, China), and the procedure used was described in the instruction manual of the kit. cDNA was synthesized from 1 µg of total RNA with Servicebio’s RT First Strand cDNA Synthesis Kit (Solarbio, Wuhan, China) for qPCR. qPCR was performed with a Q2000B (LongGene, China) using 2× SYBR Green qPCR Master mix (Solarbio, Wuhan, China). In this study, all qRT‒PCR experiments were performed in triplicate via a Q2000B Thermal Cycler (Roche, USA) for detection in strict accordance with the manufacturer’s instructions. The relative expression levels of genes were calculated via the 2-ΔΔCt method [34]. All the primers used are listed in Table S1.

A total of 20 mL of fermentation broth was centrifuged for 10 min (13,000 rpm/min, 4.0 °C), the supernatant was discarded, and the mixture was resuspended in 1 mL of 0.1 M citrate phosphate buffer. The supernatants were harvested after the cells were disrupted by an ultrasonic cell pulverizer (Scientz-IID, Beyotime, Shanghai, China) and centrifuged at 13,000 rpm/min for 20.0 min at 4.0 °C. The activities of proline dehydrogenase (putB) (KTB1431, Abbkine, Wuhan, China) and fructokinase (scrK) (G0871F, Grace Biotechnology, Suzhou, China) were assayed via a kit. The activities of N-acetylornithine aminotransferase (ArgD) [35] and oxamate amidohydrolase (ywrD) [36] were determined spectrophotometrically. Using OriginPro 2018 C statistical analysis data, differences between two independent samples were considered statistically significant when P < 0.05 according to a t test.

Results

Fermentation and sample preparation of protease PB92

A foreign expression plasmid of the target gene PB92 was constructed after verifying the accuracy of the sequence for the recombinant plasmid pBE-09-PB92 (promoter pShuttle-09 and the SP BslB), which was expressed in B. subtilis WB600. The transformed B. subtilis WB600 appeared as a clear transparent proteolytic circle after overnight incubation on a plate containing 1.0% skim milk powder (Fig. 1A). The results showed that the protease PB92 could achieve efficient expression in B. subtilis. To determine the time point of sampling, we further investigated the relationship between the extracellular enzyme activity and fermentation time of the protease PB92. By detecting the extracellular activity of the protease PB92 at 16 fermentation time points (0–96 h), the results revealed that the extracellular protease activity of the protease PB92 rapidly increased from 10 h to 60 h of fermentation, and after 66 h, the extracellular protease activity plateaued, with the highest extracellular enzyme activity reaching 4687.93 ± 191.03 U/mL in B. subtilis (Fig. 1B). The extracellular protease activity of the 3 groups of sequenced samples was 130.63 ± 4.26 U/mL, 2363.20 ± 137.32 U/mL, and 4528.53 ± 121.09 U/mL. Compared with that of the P1 group, the extracellular protease activity of the P2 group increased by 18.09-fold. Compared with that of the other groups, the P3 group reached the highest activity point and increased by 34.67-fold compared with that of the P1 group, and the P3 group increased by 1.91-fold compared with the P2 group (Fig. 1B). The growth trend of the strains showed that the P2 group was the highest, and the OD 600 nm of the P3 strain decreased (Fig. 1C).

Fig. 1
figure 1

Preparation of sequenced samples. (A) Transparent zone formed by the protease PB92. CK and PB92 represent B. subtilis WB600 strains containing empty vectors and PB92 target genes, respectively. (B) Fermentation results of the protease PB92 recombinant strain. To determine protease activity, culture samples were collected at 16 time points from 4 to -96 h. (C) Strain growth data of the sequence sample. Three parallel samples were obtained at 12 h (P1), 48 h (P2), and 60 h (P3) of fermentation

Full size image

Sequencing data analysis of RNA-sequencing and Ribo-sequencing

After RNA sequencing, filtering and alignment, 3982, 3988, and 3996 genes were detected in the P1, P2, and P3 groups, respectively, accounting for more than 93.98% of the total number of reference genome genes (Table S2). Correlation coefficient analysis was performed to assess sample reproducibility between groups. The results revealed that each set of samples has a high correlation coefficient (R ≥ 0.91) (Fig. 2A), there were no abnormal samples. Cluster analysis between the samples also revealed that each group of samples was clustered individually into the same branch (Figure S1). These results show that the samples within each group exhibit good reproducibility. The differentially expressed genes between the three groups were statistically analysed, with fold change ≥ 2 and FDR < 0.05 as the screening thresholds. The results revealed that the genes that differed between the three comparison groups were all arranged in a volcanic arrangement (Fig. S2). The comparison groups P1 vs. P2, P1 vs. P3, and P2 vs. P3 revealed 605 (403 upregulated, 202 downregulated), 1089 (697 upregulated, 392 downregulated), and 367 (228 upregulated, 139 downregulated) significantly different genes, respectively. P1 vs. P3 was also the group comparison with the greatest difference among the three comparison groups (Fig. 2B).

Fig. 2
figure 2

RNA-sequencing and Ribo-sequencing data analysis (n = 3 in each group). (A) Pearson correlation between two samples of RNA-sequencing data; (B) histogram of differentially expressed genes with fold changes ≥ 2 and FDR < 0.05 between two comparison groups; (C) Pearson correlation between two samples of Ribo-sequencing data; (D) volcano map of Ribo-sequencing data; (E) histogram of differentially expressed genes with fold changes ≥ 2 and FDR < 0.05 between two comparison groups. Abbreviations: P1 = 12 h; P2 = 48 h; P3 = 60 h; “r” stands for samples for Ribo-sequencing

Full size image

After Ribo-sequencing and filtering, the alignment rate with the reference genome reached 86.0%. A total of 4226 and 4200 genes were detected in rP1 and rP3, respectively, covering more than 99.00% of the genes and exhibiting higher sequencing quality than RNA sequencing (Table S3). The results of the Pearson correlation coefficient analysis revealed that the correlation between three replicates in the two groups of samples was extremely high, and the value reached more than 0.87 (Fig. 2C). The differentially expressed genes in the rP1 vs. rP3 groups showed a volcanic distribution (Fig. 2D). A total of 1272 differentially expressed genes (fold change ≥ 2, FDR < 0.05) were found, of which 701 were upregulated and 571 were downregulated (Fig. 2E). This result suggests that there are more differentially expressed genes at the level of translation.

Analysis of the associations of differentially expressed genes at the transcription and translation levels

To further explore the changes in gene expression at both the transcription and translation levels in B. subtilis, we compared the relative expression of 4237 genes annotated via Ribo-seq and RNA-seq in the same sample. The results revealed that the relative expression of the genes measured at the two levels was not consistent, and there was a significant difference in the expression of some genes at the two levels (Fig. 3A). Furthermore, we calculated the Pearson correlation coefficient of differentially expressed genes between translation expression and transcript abundance in the two comparison groups. The results revealed that the Pearson correlation coefficient of P1 vs. P3 at the transcription and translation levels was only 0.5354 (Fig. 3B).

Fig. 3
figure 3

Analysis of the associations of differentially expressed genes at the transcription and translation levels. (A) Comparison of the relative expression of genes. (B) Scatterplot of correlation analysis for different genes

Full size image

Integrative analysis to determine the target gene sets

To further explore the major differential gene sets associated with the efficient expression of foreign protein, we performed Venn analysis between the three comparison groups at the transcriptional level. There were 54 differentially expressed genes common to the three comparison groups (Fig. 4A). The characterization and analysis of 8 different trends of 54 highly correlated differentially expressed genes revealed significant upregulation (P = 0.00022) and significant downregulation (P = 4.9E-06) of 22 and 19 genes, respectively, with increasing fermentation time. The other types of trend plots were not significant (P > 0.05), so these genes were culled (Fig. 4B). Fortunately, the 41 differentially expressed genes (Table S5) associated with the efficient expression of foreign protein revealed by RNA-seq also showed significant differences at the level of translation via Venn analysis (Fig. 4C). On the other hand, 436 (187 upregulated, 149 downregulated) gene sets with significant differences only at the translation level were also of interest (Fig. 4D). Combined with the in vitro experimental indices between the two groups, the extracellular protease activity of the P3 group increased by 34.67 times compared with that of the P1 group, and these differential gene sets that were specific to the translation level may also be closely related to the efficient expression of foreign protein.

Fig. 4
figure 4

Determination of target gene sets. (A) Venn analysis of three comparison groups at the transcription level. (B) Trend analysis of the target gene set. (C) Venn analysis of all differentially expressed genes at the transcription level and the target gene set at the translation level. (D) Volcano map of uniquely different gene sets at the translation level. Abbreviations: P1 = 12 h; P2 = 48 h; P3 = 60 h

Full size image

Pathway enrichment analysis of the target gene set

To reveal the pathways of the key differential gene sets, the main differential gene sets were analysed by KO enrichment in units of the KEGG pathway. The results revealed the top 10 pathways with significant enrichment of 41 common differentially expressed genes at the transcription and translation levels, of which the first 9 pathways belong to metabolism: starch and sucrose metabolism; nitrogen metabolism; porphyrin metabolism; arginine and proline metabolism; 2-oxocarboxylic acid metabolism; oxidative phosphorylation; degradation of aromatic compounds; glycolysis/gluconeogenesis; biosynthesis of secondary metabolites; and the membrane phosphotransferase system (PST) pathway under the transport pathway (Fig. 5A). Sixteen significantly different genes were involved. The top 10 pathways were significantly enriched in 436 differentially expressed genes unique to the translation level, and the results revealed that 7 pathways were associated with metabolism: other glycan degradation, taurine and hypotaurine metabolism, propanoate metabolism, sulfur metabolism, arginine biosynthesis, lysine biosynthesis, and chlorocyclohexane and chlorobenzene degradation. There were two pathways involved in cell motility, flagellar assembly and bacterial chemotaxis, and a two-component system under the signal transduction pathway (Fig. 5B). Forty-seven significantly different genes were involved.

Fig. 5
figure 5

Top 10 significantly enriched pathways of the target differential gene sets. (A) Forty-one differentially expressed genes at both the transcription and translation levels; (B) 436 differentially expressed genes unique to the translation level

Full size image

Functional analysis of target genes

The pathways associated with the 16 major differentially expressed genes and the transformation process of the compounds were further revealed (Fig. 6). Two genes (putB and putC) with significantly upregulated expression were mainly involved in the regulation of the arginine and proline synthesis pathways (Fig. 6A). Only the membrane transport-related pathway (PTS) associated with the genes gmuA, gmuB, and gmuC was significantly downregulated, which was related mainly to the transmembrane transport of carbon sources (Fig. 6B). Carbohydrate metabolism involves two main pathways: glycolysis/gluconeogenesis and starch and sucrose metabolism. The former involves two differentially expressed genes, gmuD (downregulated) and acsA (upregulated) (Fig. 6C), and the latter involves five differentially expressed genes, gmuA, gmuB, gmuC, gmuD, and scrK. These genes were all downregulated, and these two pathways are related to glucose utilization (Fig. 6D). Energy metabolism involves two pathways, nitrogen metabolism and oxidative phosphorylation, which involve four downregulated differential genes, ctaA, ctaB2, narI, and narK, of which the ctaA and ctaB2 genes are involved in the oxidative phosphorylation pathway (Fig. 6E). These two genes also regulate porphyrin metabolism, which is related mainly to ATP synthesis (Fig. 6H). The narI and narK genes are involved in the nitrogen metabolism cycle and are ultimately closely related to the amino acid metabolism pathway (Fig. 6F). The biosynthesis of secondary metabolites pathway (bsdC, downregulated) and the degradation of aromatic compounds pathway (suhB, upregulated) also contain one gene each. The biosynthesis of secondary metabolites pathway was primarily involved in D-glucose-6P synthesis, which in turn participates in carbon metabolism (Fig. 6G). The pathway of degradation of aromatic compounds is associated with CoA, which further participates in the TCA cycle (Fig. 6I). Two significantly upregulated genes (argD and ilvD) associated with the 2-oxocarboxylic acid metabolism pathway were involved mainly in the synthesis of ornithine, valine, and isoleucine, which were in turn associated with the TCA cycle (Fig. 6J).

Fig. 6
figure 6

Functions of 16 differentially expressed genes and corresponding pathways. (A) Arginine and proline metabolism; (B) phosphotransferase system; (C) glycolysis/gluconeogenesis; (D) starch and sucrose metabolism; (E) oxidative phosphorylation; (F) nitrogen metabolism; (G) biosynthesis of secondary metabolites; (H) porphyrin metabolism; (I) degradation of aromatic compounds; (J) 2-oxocarboxylic acid metabolism. The expression of differentially expressed genes is represented by a heatmap; each column in the figure represents a sample, and different colours represent different amounts of expression of genes in the sample, with blue representing downregulation and red representing upregulation

Full size image

The pathways associated with the 47 major differentially expressed genes and the transformation process of the compounds were further revealed (Fig. 7). In the two-component system pathway, 16 differentially expressed genes were involved, of which 12 genes (comX, dnaA, bceB, dltC, dltD, cheB, tlpA, mcpA, cheA, cheW, eglS, and flgM) were significantly downregulated, while 4 genes (phoB, sacB, yodR, and yodS) were significantly upregulated (Fig. 7A). Among the cell motility pathways, flagellar assembly and bacterial chemotaxis were involved with 18 differentially expressed genes, and the flagellar assembly pathway (fliK, flgD, flgE, flhA, fliF, fliG, fliH, fliI, fliQ, fliJ, and flgM) and bacterial chemotaxis pathway (tlpA, mcpA, cheB, cheD, cheA, fliG, cheW, and cheC) showed a significant downwards trend (Fig. 7B and C). These 32 genes respond to changes in the external environment (orthophosphate, temperature, bacitracin, cAMP). In the amino acid metabolism pathway, arginine biosynthesis (yodQ, upregulated; ureB, rocG, and aspB, downregulated) and lysine biosynthesis (lysA, upregulated; murF and lysC, downregulated) were the main pathways, and arginine and lysine participated in the TCA cycle through the mediation of glutamate (Fig. 7D and E). Six differentially expressed genes (yngHB, yodR, and yodS, upregulated; pccB, ackA, and ytbE, downregulated) are involved in the propanoate metabolism pathway and participate in the TCA cycle through the production of succinyl-coA and acetyl-coA (Fig. 7F). Sulfur metabolism (tauC, yisZ, and yitB, upregulated; nrnA, downregulated) and taurine and hypotaurine metabolism (ywrD and rocG, upregulated; ackA, downregulated) participate in the TCA cycle through the production of acetate and isoleucine (Fig. 7G and H). In addition, sulfur metabolism can also be involved in glycolysis by synthesizing cysteine (Fig. 7H). Chlorocyclohexane and chlorobenzene degradation by the production of glycolate (Fig. 7J) and the degradation of other glycans (yesZ, upregulated; lytD, downregulated) by the degradation of carbon sources (Fig. 7I) are involved in glycolysis.

Fig. 7
figure 7

Functions of 47 differentially expressed genes and corresponding pathways. (A) Two-component systems; (B) flagellar assembly; (C) bacterial chemotaxis; (D) arginine biosynthesis; (E) lysine biosynthesis; (F) propionic acid metabolism; (G) taurine and hypotaurine metabolism; (H) sulfur metabolism; (I) degradation of chlorocyclohexane and chlorobenzene; (J) degradation of other glycans. The expression of differentially expressed genes is represented by a heatmap; each column in the figure represents a sample, and different colours represent different amounts of expression of genes in the sample, with blue representing downregulation and red representing upregulation

Full size image

Validation of partial genes and proteins

To validate the RNA-seq results, we randomly selected ten differentially expressed genes and examined their expression patterns at the two fermentation stages via q‒PCR. Among them, the differentially expressed genes argD, putB, gmuC, ctaA, and scrK had significant differences in transcription and translation levels (Fig. 8A and B), and the differentially expressed genes ywrD, fliI, dltC, cheC, and aspB had significant differences only in translation level, and there were no significant differences in transcription level (Fig. 8A and B). Five genes (argD, putB, gmuC, ctaA, and scrK) were differentially expressed (P < 0.05) at different stages according to q‒PCR. In addition, the expression of five genes tended to be similar (argD and putB, upregulated; gmuC, ctaA, and scrK, downregulated) to the results of RNA-seq and q-PCR (Fig. 8C). However, there was still no significant difference (P > 0.05) in the remaining five genes (ywrD, fliI, dltC, cheC, and aspB) between the results of RNA-seq and q-PCR (Fig. 8D). In addition, we further validated the protein activity of four differentially expressed genes (argD, putB, scrK, and ywrD). The results revealed significant differences in the protein activities of the four differentially expressed genes between the two groups (Fig. 8E) (P < 0.05). The activity of these four proteins showed a similar trend (argD, putB, and ywrD, upregulated; scrK, downregulated) to the Ribo-seq results (Fig. 8E and B). However, although the ywrD gene was verified by RNA-seq (Fig. 8A) and q-PCR (Fig. 8D), there was no significant difference (P > 0.05). However, there was significant variation in protein activity (Fig. 8E) (P = 0.0258). Therefore, the FPKM obtained from RNA-seq can be reliably used to determine gene expression, and Ribo-seq more accurately measures the protein activity of the gene of interest.

Fig. 8
figure 8

Validation of partial genes via q‒PCR (n = 9 in each group). (A) Heatmap analysis of 10 target genes at the transcriptional level in the sample. (B) Heatmap analysis of 10 target genes at the translational level in the sample. (C) The relative expression of five target genes with differences in transcription and translation levels was measured via q‒PCR. (D) The relative expression of five target genes with differences in only translation levels was measured via q‒PCR. (E) The protein activities of the four target genes were compared. Abbreviations: P1 = 12 h; P3 = 60 h; “r” stands for samples for Ribo-sequencing. Each column in the figure represents a sample, each row represents a gene, and different colours represent different amounts of gene expression in the sample. Using OriginPro 2018 C statistical analysis data, the difference between the two independent samples was statistically significant at P < 0.05

Full size image

Discussion

In this study, Ribo-seq and RNA-seq technology were used to study the changes in differentially expressed genes during efficient secretion of the protease PB92 via B. subtilis and revealed the differential gene regulatory network. However, several points are worth discussing.

First, RNA-seq technology can be used to study the expression, phenotype, and function of different genes at the transcription level and determine the number and level of genes expressed, and it can detect dynamic changes in transcription levels from different tissues or different states, which has a strong advantage in gene expression differences [37]. Previous studies revealed differentially expressed genes associated with osmotic pressure and glucose metabolism in B. subtilis [38, 39], and the metabolic pathway changes in B. subtilis during adenosine fermentation and isoprene production were analysed via RNA-seq technology [40, 41]. In this study, 54 differentially expressed genes related to the efficient expression of foreign protein in B. subtilis were also identified via RNA-seq via the target gene of the protease PB92 (Fig. 4A).

Second, genes can be transcribed into mRNAs, which does not necessarily mean that they can be translated; therefore, the correlation between the transcription level and metabolome is not high [13]. Maier et al. compared transcriptomes and proteomes across multiple species and reported that the correlation coefficient R2 between mRNA levels and the amount of protein they encode ranged from 0.01 to 0.5 [42]. Known as the “king of regulation”, translation regulation is the most important mode of regulation, accounting for more than half of the total regulation ratio [13]. Ribo-seq is closer to the true protein than RNA-seq is and may be more reflective of genetic differences in regulatory processes [13]. Our results revealed the correlation (0.5354) of the differentially expressed genes at the transcription and translation levels (Fig. 3). This result also proves that differential gene regulatory networks are also poorly correlated at the transcription and translation levels of B. subtilis. We theorize that the possible cause of this phenomenon is a translation pause involved in the translation of mRNAs in B. subtilis. We detected translation pauses in 1036 and 423 genes in two sets of samples (rP1 and rP3), respectively.

Third, in this study, the key gene set for efficient expression of foreign protein in B. subtilis was revealed via integrative analysis of Ribo-seq and RNA-seq data (Fig. 4). Two-component systems, flagellar assembly, and bacterial chemotaxis are involved in regulating a variety of physiological processes, including chemotaxis, substance metabolism, and transport [43,44,45]. Studies have shown that bacteria respond mainly to changes in the external environment through chemotaxis and two-component regulatory systems to respond to changes in their own motor state to obtain more favourable living conditions, which maybe include foraging, growing, and information exchanging, among others [43,44,45]. We speculate that this may be one of the key factors in the significant difference in foreign protein expression in B. subtilis, with bacteria manifesting greater chemotaxis when foreign protein is highly expressed. However, there have been no studies on the modification of the B. subtilis metabolic pathway to increase the expression of foreign proteins from a chemotaxis perspective, which may be related to the fact that the associated differentially expressed genes are only expressed at the translational level, while RNA-seq cannot be performed.

Finally, previous studies have shown that B. subtilis maximizes the utilization and allocation of carbon sources by coordinating the utilization and allocation of carbon sources in response to changes in the type and content of carbon sources in the environment through the central carbon metabolism pathway [46, 47]. This study also revealed that important central carbon metabolism pathways, such as the PTS system, glycolysis/gluconeogenesis, starch and sucrose metabolism, and other glycan degradation pathways, regulate the efficient expression of foreign protein of B. subtilis (Figs. 5 and 6). The function of the PTS system is to transport various carbon sources from the extracellular space to the intracellular space by phosphorylation of the phosphoric acid cascade, and glucose can be absorbed and utilized by cells mainly through the mediation of this system [48]. Studies have shown that by blocking the PTS pathway and strengthening the nonsugar phosphotransferase system pathway, the carbon metabolism overflow in B. subtilis can be weakened, and the yield of the product N-ethyl glucosamine can be increased from 6.5 g/L to 10.1 g/L [49]. Inhibition of the pfkA gene of the central carbon metabolism pathway can achieve co-utilization of glucose and xylose by B. subtilis and increase the yield of the product GlcNAc by 84.1% [50]. These results show that the rational modification of central carbon metabolism favours the expression of foreign proteins, which is consistent with our research results. The regulation of nitrogen metabolism in B. subtilis also involves a complex and fine regulatory network [46, 47]. This study revealed that nitrogen metabolism, 2-oxocarboxylic acid metabolism, arginine biosynthesis, lysine biosynthesis, and taurine and hypotaurine metabolism mediate the regulation of central nitrogen metabolism through the synthesis and degradation pathways of arginine, proline, lysine, and valine (Fig. 7). By knocking out the transcription factor CodY in B. subtilis, Dhali et al. removed the inhibitory effect of CodY on the synthesis of branched-chain amino acids (leucine), increasing the yield of surfactin by a factor of 10, reaching 2.3 g/L [51]. This result also reflects the importance of B. subtilis central nitrogen metabolism for high-density fermentation of foreign proteins. In addition, this study revealed that sulfate is converted to cysteine and acetate through sulfur metabolism and then participates in the central metabolic pathway. The e– produced by the sulfur metabolism pathway enters the electron transport chain and participates in energy metabolism [52]. Propionate is produced by propanoate metabolism after succinyl-CoA and acetyl-CoA and then participates in the TCA cycle. In addition, porphyrin metabolism is involved mainly in oxidative phosphorylation, which is the main way substances release energy in the body by synthesizing ATP through the respiratory chain [53]. In addition to their primary metabolism, microorganisms can secrete a variety of secondary metabolites, such as D-myo-inositol and other valuable growth factors, through the biosynthesis of secondary metabolites [54]. Interestingly, we also identified two pathways associated with the degradation of toxic and harmful substances, degradation of aromatic compounds, and chlorocyclohexane and chlorobenzene degradation [55], and we speculate that low doses of aromatic or halogenated compounds may also be involved as nutrients in the efficient expression of foreign proteins. In summary, we speculate that the main regulatory mechanism for the efficient expression of foreign genes in B. subtilis is to obtain more favourable nutrient resources and living conditions on the basis of cell chemotaxis, supply a biochemical reaction network dominated by the central metabolic pathway, maximize nutrient utilization, and improve the target product (Fig. 9).

Fig. 9
figure 9

Schematic diagram. A schematic showing the main regulatory mechanism for the efficient expression of foreign genes in B. subtilis: obtaining more favourable nutrient resources and living conditions on the basis of cell chemotaxis, supplying a biochemical reaction network dominated by the central metabolic pathway, maximizing nutrient utilization, and improving the target product

Full size image

Conclusions

Overall, new insights into the gene regulatory network of efficient secretion of foreign protein in B. subtilis were obtained via integrative analyses of multiomics joint analysis. To the best of our knowledge, this is the first study to reveal a gene regulatory network in B. subtilis via Ribo-seq and RNA-seq. The cell motility and central metabolic regulatory network constitute the key regulatory mode for the efficient expression of foreign proteins, and this new insight provides important reference value for the construction of cell factories and metabolic networks for the efficient secretion of target products by B. subtilis.

Data availability

The raw sequence data reported in this paper have been deposited in the Genome Sequence Archive [56] in National Genomics Data Center [57], China National Center for Bioinformation / Beijing Institute of Genomics, Chinese Academy of Sciences (GSA: CRA012098) that are publicly accessible at https://ngdc.cncb.ac.cn/gsa/s/5PKYxPYR. The other data that support the findings in this study were not deposited in an official repository, while these data are available from authors upon request.

References

  1. Gu Y, Xu XH, Wu YK, et al. Advances and prospects of Bacillus subtilis cellular factories: From rational design to industrial applications. Metab Eng. 2018;50:109–21.

    Article  PubMed  CAS  Google Scholar 

  2. Yzturk S, Yalık P, Yzdamar HT. Fed-batch biomolecule production by Bacillus subtilis: A state of the art review. Trends Biotechnol. 2016;34:329–45.

    Article  PubMed  Google Scholar 

  3. Zweers JC, Barak I, Becher D, et al. Towards the development of Bacillus subtilis as a cell factory for membrane proteins and protein complexes. Microb Cell Fact. 2008;7:10.

    Article  PubMed  PubMed Central  Google Scholar 

  4. Park SA, Bhatia SK, Park HA, et al. Bacillus subtilis as a robust host for biochemical production utilizing biomass. Crit Rev Biotechnol. 2021;41:31–4.

    Article  Google Scholar 

  5. Xiang MJ, Kang Q, Zhang DW. Advances on systems metabolic engineering of Bacillus subtilis as a chassis cell. Synth Syst Biotechnol. 2020;5:245–51.

    Article  PubMed  PubMed Central  Google Scholar 

  6. Sanchez S, Demain LA. Metabolic regulation of fermentation processes. Enzyme Microb Tech. 2002;31:895–906.

    Article  CAS  Google Scholar 

  7. Yup SL, Kim HU. Systems strategies for developing industrial microbial strains. Nat biotechnol. 2015;33:1061–72.

    Article  Google Scholar 

  8. Rychel K, Sastry AV, Palsson BO. Machine learning uncovers independently regulated modules in the Bacillus subtilis transcriptome. Nat Commun. 2020;11:6338.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  9. Faria JP, Overbeek R, Taylor RC, et al. Reconstruction of the regulatory network for Bacillus subtilis and reconciliation with gene expression data. Front Microbiol. 2016;7:25.

    Article  Google Scholar 

  10. Adeli K. Translational control mechanisms in metabolic regulation: Critical role of RNA binding proteins, microRNAs, and cytoplasmic RNA granules. Am J Physiol Endocrinol Metab. 2011;301:1051–64.

    Article  Google Scholar 

  11. Ansari SA, Dantoft W, Ruiz-Orera J, et al. Integrative analysis of macrophage Ribo-Seq and RNA-Seq data define glucocorticoid receptor regulated inflammatory response genes into distinct regulatory classes. Comput Struct Biotechnol J. 2022;20:5622–38.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  12. Zupanic A, Meplan C, Grellscheid SN, et al. Detecting translational regulation by change point analysis of ribosome profiling data sets. RNA. 2014;20:1507–18.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  13. Schwanhausser B, Busse D, Li N, Dittmar G, et al. Global quantification of mammalian gene expression control. Nature. 2011;473:337–42.

    Article  PubMed  Google Scholar 

  14. Schwanhäusser B, Busse D, Li N, et al. Genome-wide analysis in vivo of translation with nucleotide resolution using ribosome profiling. Science. 2009;324:218–23.

    Article  Google Scholar 

  15. Ingolia NT, Brar GA, Rouskin S, et al. The ribosome profiling strategy for monitoring translation in vivo by deep sequencing of ribosome-protected mRNA fragments. Nat Protoc. 2012;7:1534–50.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  16. Bagheri A, Astafev A, Al-Hashimy T, et al. Tracing translational footprint by Ribo-Seq: Principle, workflow, and applications to understand the mechanism of human diseases. Cells. 2022;11:2966.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  17. Siezen RJ, de Vos WM, Leunissen JA, et al. Homology modelling and protein engineering strategy of subtilases, the family of subtilisin-like serine proteinases. Protein Eng. 1991;4:719–37.

    Article  PubMed  CAS  Google Scholar 

  18. van der Laan JM, Teplyakov AV, Kelders H, et al. Crystal structure of the high-alkaline serine protease PB92 from Bacillus alcalophilus. Protein Eng. 1992;5:401–11.

    Google Scholar 

  19. Martin JR, Mulder FA, Karimi-Nejad Y, et al. The solution structure of serine protease PB92 from Bacillus alcalophilus presents a rigid fold with a flexible substrate-binding site. Structure. 1997;5:521–32.

    Article  PubMed  CAS  Google Scholar 

  20. Miao HB, Xiang X, Han NY, et al. Improving the thermostability of serine protease PB92 from Bacillus alcalophilus via site-directed mutagenesis based on semi-rational design. Foods. 2023;12:3081.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  21. Miao HB, Zhe YY, Xiang X, et al. Enhanced extracellular expression of a Ca2+- and Mg2+-dependent hyperthermostable protease EA1 in Bacillus subtilis via systematic screening of optimal signal peptides. JAFC. 2022;70:15830–9.

    CAS  Google Scholar 

  22. Miao HB, Jiang R, Han NY, et al. Enhanced extracellular expression of α-amylase DL3-4-1 in Bacillus subtilis via systematic screening of optimal signal peptides. Process Biochem. 2021;108:176–84.

    Article  CAS  Google Scholar 

  23. Arenz S, Wilson DN. Bacterial protein synthesis as a target for antibiotic inhibition. CSH Perspect Med. 2016;6(9):a025361.

    Google Scholar 

  24. Champney WS. Antibiotics targeting bacterial ribosomal subunit biogenesis. J Antimicrob Chemoth. 2020;75(4):787–806.

    Article  CAS  Google Scholar 

  25. Rio DC, Ares M Jr, Hannon GJ, Nilsen TW. Purification of RNA using TRIzol (TRI reagent). Cold Spring Harb Protoc. 2010;6:5439.

  26. Chen SF, Zhou YQ, Chen YR, et al. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34:884–90.

    Article  Google Scholar 

  27. Barbe V, Cruveiller S, Kunst F, et al. From a consortium sequence to a unified sequence: the Bacillus subtilis 168 reference genome a decade later. Microbiology. 2009;155:1758–75.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  28. Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9:357–9.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  29. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinform. 2011;12:323.

    Article  CAS  Google Scholar 

  30. Han XX. On statistical measures for data quality evaluation. J Geogr Inf Syst. 2020;12(3):178–87.

    Google Scholar 

  31. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Zhong Y, Karaletsos T, Drewe P, Sreedharan VT, Kuo D, Singh K, Wendel HG, Ratsch G. RiboDiff: detecting changes of mRNA translation efficiency from ribosome footprints. Bioinformatics. 2017;33(1):139–41.

    Article  PubMed  CAS  Google Scholar 

  33. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2002;28:27–30.

    Article  Google Scholar 

  34. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) method. Methods. 2001;25:402–8.

    Article  PubMed  CAS  Google Scholar 

  35. Morizono H, Cabrera-Luque J, Shi DS, et al. Acetylornithine transcarbamylase: a novel enzyme in arginine biosynthesis. J Bacteriol. 2006;188(8):2974–82.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  36. Hicks KA, Ealick SE. Biochemical and structural characterization of Klebsiella pneumoniae oxamate amidohydrolase in the uric acid degradation pathway. Acta Cryst. 2016;D72:808–16.

    Google Scholar 

  37. Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009;10(1):57–63.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  38. Hahne H, Mader U, Otto A, et al. A comprehensive proteomics and transcriptomics analysis of Bacillus subtilis salt stress adaptation. J Bacteriol. 2010;192:870–82.

    Article  PubMed  CAS  Google Scholar 

  39. Gutierrez-Rios RM, Freyre-Gonzalez JA, Resendis O, et al. Identification of network topological units coordinating the global expression response to glucose in Bacillus subtilis and its comparison to Escherichia coli. BMC Microbiol. 2009;9:176.

    Article  PubMed  PubMed Central  Google Scholar 

  40. Yu WB, Gao SH, Yin CY et al. Comparative transcriptome analysis of Bacillus subtilis responding to dissolved oxygen in adenosine fermentation. PLoS ONE. 2011;6:e20092.

  41. Shapiro JA. Revisiting the central dogma in the 21st century. Ann Ny Acad Sci. 2009;1178:6–28.

    Article  PubMed  CAS  Google Scholar 

  42. Maier T, Guell M, Serrano L. Correlation of mRNA and protein in complex biological samples. FEBS Lett. 2009;583:3966–73.

    Article  PubMed  CAS  Google Scholar 

  43. Ortet P, Whitworth DE, Santaella C, et al. P2CS: updates of the prokaryotic two-component systems database. Nucleic Acids Res. 2015;43:536–41.

    Article  Google Scholar 

  44. Keegstra JM, Carrara F, Stocker R. The ecological roles of bacterial chemotaxis. Nat Rev Microbiol. 2022;20:491–504.

    Article  PubMed  CAS  Google Scholar 

  45. Cremer J, Honda T, Tang Y, et al. Chemotaxis as a navigation strategy to boost range expansion. Nature. 2019;575:658–63.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  46. Fujita Y. Carbon catabolite control of the metabolic network in Bacillus subtilis. Biosci Biotech Bioch. 2009;73:245–59.

    Article  CAS  Google Scholar 

  47. Gunka K, Commichau FM. Control of glutamate homeostasis in Bacillus subtilis: acomplex interplay between ammonium assimilation, glutamate biosynthesis anddegradation. Mol Microbiol. 2012;85:213–24.

    Article  PubMed  CAS  Google Scholar 

  48. Reizer J, Bachem S, Reizer A et al. Novel phosphotransferase system genes revealed by genome analysis-the complete complement of PTS proteins encoded within the genome of Bacillus subtilis. Microbiology. 1999;145:3419–29.

  49. Gu Y, Deng JY, Liu YF, et al. Rewiring the glucose transportation and central metabolicpathways for overproduction of N-Acetylglucosamine in Bacillus subtilis. Biotechnol J. 2017;12:1700020.

    Article  Google Scholar 

  50. Wu YK, Chen TC, Liu YF, et al. CRISPRi allows optimal temporal control of N-acetylglucosamine bioproduction by a dynamic coordination of glucose and xylosemetabolism in Bacillus subtilis. Metab Eng. 2018;49:232–41.

    Article  PubMed  CAS  Google Scholar 

  51. Dhali D, Coutte F, Arias AA, et al. Genetic engineering of the branched fatty acidmetabolic pathway of Bacillus subtilis for the overproduction of surfactin C14 isoform. Biotechnol J. 2017;12:1600574.

    Article  Google Scholar 

  52. Mansilla MC, Albanesi D, Mendoza DD. Transcriptional control of the sulfur-regulated cysH operon, containing genes involved in L-cysteine biosynthesis in Bacillus subtilis. J Bacteriol. 2000;182:5885–92.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  53. Sousa J, Westhoff P, Methling K, et al. The absence of pyruvate kinase affects glucose-dependent carbon catabolite repression in Bacillus subtilis. Metabolites. 2019;9:216.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  54. Iqbal S, Begum F, Rabaan A, et al. Classification and multifaceted potential of secondary metabolites produced by Bacillus subtilis group: A comprehensive review. Molecules. 2023;28:927.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  55. Kumar AG, Manisha D, Rajan NN, et al. Biodegradation of phenanthrene by piezotolerant Bacillus subtilis EB1 and genomic insights for bioremediation. Mar Pollut Bull. 2023;194:115151.

    Article  Google Scholar 

  56. Chen TT, Chen X, Zhang SS, et al. The genome sequence archive family: Toward explosive data growth and diverse data types. GPB. 2021;19:578–83.

    PubMed  PubMed Central  Google Scholar 

  57. CNCB-NGDC members and partners. Database resources of the National Genomics Data Center. China National Center for bioinformation in 2022. Nucleic Acids Res. 2022;50:27–38.

    Article  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This work was financed by the colleges and universities in Yunnan Province serve key industrial science and technology projects (FWCY-ZNT2024009) and the National Natural Science Foundation of China (Grant No. 32260015).

Author information

Author notes

  1. Huabiao Miao and Xia Xiang contributed equally to this paper.

Authors and Affiliations

  1. School of Life Science, Yunnan Normal University, Kunming, 650500, China

    Huabiao Miao, Xia Xiang, Qian Wu & Zunxi Huang

  2. Engineering Research Center for Efficient Utilization of Characteristic Biological Resources in Yunnan, Ministry of Education, Kunming, 650500, China

    Huabiao Miao, Qian Wu & Zunxi Huang

  3. Key Laboratory of Yunnan for Biomass Energy and Biotechnology of Environment, Kunming, 650500, China

    Huabiao Miao, Qian Wu & Zunxi Huang

  4. Guangzhou Genedenovo Biotechnology Company Limited, Guangzhou, 510000, China

    Ling Cheng

Authors
  1. Huabiao Miao
    View author publications

    Search author on:PubMed Google Scholar

  2. Xia Xiang
    View author publications

    Search author on:PubMed Google Scholar

  3. Ling Cheng
    View author publications

    Search author on:PubMed Google Scholar

  4. Qian Wu
    View author publications

    Search author on:PubMed Google Scholar

  5. Zunxi Huang
    View author publications

    Search author on:PubMed Google Scholar

Contributions

HBM, XX designed the study, and performed the majority of experiments. HBM conducted all sampling, and extracted DNA. LC, QW analysed the data. HBM, XX, ZXH wrote the paper, and critically revised the manuscript. All authors discussed and approved the data and reviewed the manuscript.

Corresponding author

Correspondence to Zunxi Huang.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic Supplementary Material

Below is the link to the electronic supplementary material.

Supplementary Material 1

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Miao, H., Xiang, X., Cheng, L. et al. New insights into the efficient secretion of foreign protein in Bacillus subtilis via Ribo-seq and RNA-seq integrative analyses. BMC Microbiol 24, 537 (2024). https://doi.org/10.1186/s12866-024-03700-y

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12866-024-03700-y

Keywords

BMC Microbiology

ISSN: 1471-2180

Contact us