Time-Course Transcriptome Analysis of Bacillus subtilis DB104 during Growth

Bacillus subtilis DB104, an extracellular protease-deficient derivative of B. subtilis 168, is widely used for recombinant protein expression. An understanding of the changes in gene expression during growth is essential for the commercial use of bacterial strains. Transcriptome and proteome analyses are ideal methods to study the genomic response of microorganisms. In this study, transcriptome analysis was performed to monitor changes in the gene expression level of B. subtilis DB104 while growing on a complete medium. Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis, K-mean cluster analysis, gene ontology (GO) enrichment analysis, and the function of sigma factors were used to divide 2122 differentially expressed genes (DEGs) into 10 clusters and identified gene functions according to expression patterns. The results of KEGG pathway analysis indicated that ABC transporter is down-regulated during exponential growth and metabolic changes occur at the transition point where sporulation starts. At this point, several stress response genes were also turned on. The genes involved in the lipid catabolic process were up-regulated briefly at 15 h as an outcome of the programmed cell death that postpones sporulation. The results suggest that changes in the gene expression of B. subtilis DB104 were dependent on the initiation of sporulation. However, the expression timing of the spore coat gene was only affected by the relevant sigma factor. This study can help to understand gene expression and regulatory mechanisms in B. subtilis species by providing an overall view of transcriptional changes during the growth of B. subtilis DB104.


Introduction
Bacillus subtilis is a Gram-positive endospore-forming bacterium and is not considered pathogenic or toxigenic to humans and animals [1,2]. Because of its well-known genetic background and easy genetic manipulation, it is an ideal strain as a host for metabolic engineering or industrial enzyme. B. subtilis DB104, a protease-deficient derivative of B. subtilis 168, is widely used for recombinant protein expression [3].
Understanding bacterial gene expression is essential for the development of a strain or vector for industrial use. To date, microarray technology has been used to monitor bacterial gene expression. Patten et al. [4] described alternative sigma factor RpoS-dependent gene expression in Escherichia coli K-12 using a microarray, while Kitichalermkiat et al. [5] explored the effect of epigallocatechin gallate on Staphylococcus aureus gene expression using this technology. During growth in specific conditions, such as with glucose [6] and batch culture [7], transcript profiling was also conducted on B. subtilis using microarrays. Because the technique relies on hybridization, several disadvantages exist, such as background hybridization that could limit the accuracy and difficulty of designing a large number of probes. RNA-sequencing (RNA-seq) is another technique that provides information about gene expression [8]. RNA-seq based on next-generation sequencing is a widely used transcriptome analysis technology. However, it has limitations of information loss due to the The transcriptome sequences were mapped to the B. subtilis subsp. subtilis str. 168 genome (GenBank accession number: AL009126.3) using HTseq version 0.10.0 and Bowtie 1.1.2 software. All raw transcriptome data are publicly available in the NCBI Sequence Read Archive (https://ncbi.nlm.nih.gov/sra, accessed on 13 February 2023) database under the accession number PRJNA934277.

Analysis of DEGs
The expression levels of each gene were calculated using the fragments per kilobase of transcript per million fragments mapped (FPKM) value, based on the mapped genes and using Bowtie software. Next, the differentially expressed genes (DEGs) were extracted using edgeR software with a p-value < 0.05 and |fold change| ≥ 2. DEGs were obtained from five compared groups ( Afterward, all DEGs were grouped into 10 K-mean clusters based on the Z-score.

Functional Analysis of DEGs
All DEGs were enriched through the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (https://www.genome.jp/kegg, accessed on 20 December 2022). Meanwhile, the DEGs from each cluster that exhibited distinct expression patterns were analyzed with the Biological Networks Gene Ontology (BiNGO) tool to identify their function. The DEGs were annotated using then Gene Ontology (GO) database according to p-value < 0.03 and Benjamini-Hochberg-corrected p-value < 0.03.

Reverse Transcription-Quantitative Polymerase Chain Reaction (RT-qPCR) Confirmation
The primers used for RT-qPCR are listed in Table S2 (Supplementary File S2). A total of 8 genes (cotD, spoIIID, yqfX, spmH, gsiB, frlO, cheY, and gbsB) that were found to be most up-and down-regulated DEGs according to fold change were selected and validated by RT-qPCR. One µg of RNA was used to synthesize cDNA using random hexamer (Roche, Base, Switzerland) and M-MLV reverse transcriptase (Promega, Madison, SI, USA), following the manufacturer's protocols. RT-qPCR was conducted at 95 • C for 10 s, followed by 40 cycles of 95 • C for 5 s and 60 • C for 30 s using Taq Pro Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China) on a CFX Connect Real-Time System (Bio-Rad, Hercules, CA, USA). The RT-qPCR analysis was conducted in triplicate independently. The target genes were normalized against rpsJ gene as an internal control. The relative expression levels of the target genes were calculated with the 2 −∆∆Ct method [23].

Cell Growth and Sporulation
Growth was monitored while cultured in an LB medium to investigate time-related transcriptome changes in B. subtilis DB104. First, the growth curve was monitored by optical density measurements at OD 600 ( Figure 1). After 12 h, cell density declined for a while and then from 15 h increased again. Next, sporulation was calculated at various times by plate counting (Figure 1) and observed using an optical microscope (Supplementary File S1: Figure S1). Sporulation was approximately 1% until 12 h, and then rapidly progressed to 18% at 15 h. After 15 h, the sporulation rate constantly increased to 90% at 24 h. Therefore, B. subtilis DB104 growth could be divided into two parts: exponential growth and sporulation with a transition point at 12 h.

RNA Sequencing and Analysis of Differentially Expressed Genes (DEGs)
Cell samples from six-time points were obtained based on the growth curve: the middle of exponential growth (8 h), the end of exponential growth (10 h), the beginning of sporulation (12 h), the former part of sporulation (15 h), the latter part of sporulation (18 h), and the end of sporulation (24 h). At each time point, the total RNAs were extracted from three independent cultures and used for RNA-seq. The transcriptome sequences were mapped against the B. subtilis subsp. subtilis str. 168 reference genome (GenBank accession number: AL009126.3). A summary of the resulting sequence data is shown in Table S1 (Supplementary File S2). The filtered clean reads ranged from 21,238,576 to 35,756,134. The Q20 values were greater than 98.19%, Q30 values were greater than 94.16%, and the GC contents were greater than 44.58% in each sample.
To ensure the reliability of the transcriptome analysis, multidimensional scaling and the Pearson correlation coefficient were used to confirm the differences among all samples (Figure 2a,b). According to the results, the 10 h and 12 h groups had similar RNA expression patterns. The 15 h, 18 h, and 24 h groups were similar, whereas the 8 h group was relatively separated from the other groups. This suggests that the transcription profile of B. subtilis DB104 was altered, depending on the growth stage.
With a p-value < 0.05 and |fold change| ≥ 2 as the criteria, a total of 2122 DEGs were identified in all the comparison groups. The numbers of DEGs identified between the various time points during the growth of B. subtilis DB104 were analyzed using a Venn diagram ( Figure 2c). The number of genes that were up-or down-regulated in each compared group is shown in Figure 2d. The results show that more transcriptome changes occurred during exponential growth to the former part of sporulation than the latter part and end of sporulation. In general, the number of up-regulated DEGs was higher than that of down-regulated DEGs.

RNA Sequencing and Analysis of Differentially Expressed Genes (DEGs)
Cell samples from six-time points were obtained based on the growth curve: the middle of exponential growth (8 h), the end of exponential growth (10 h), the beginning of sporulation (12 h), the former part of sporulation (15 h), the latter part of sporulation (18 h), and the end of sporulation (24 h). At each time point, the total RNAs were extracted from three independent cultures and used for RNA-seq. The transcriptome sequences were mapped against the B. subtilis subsp. subtilis str. 168 reference genome (GenBank accession number: AL009126.3). A summary of the resulting sequence data is shown in Table S1 (Supplementary File S2). The filtered clean reads ranged from 21,238,576 to 35,756,134. The Q20 values were greater than 98.19%, Q30 values were greater than 94.16%, and the GC contents were greater than 44.58% in each sample.
To ensure the reliability of the transcriptome analysis, multidimensional scaling and the Pearson correlation coefficient were used to confirm the differences among all samples (Figure 2a,b). According to the results, the 10 h and 12 h groups had similar RNA expression patterns. The 15 h, 18 h, and 24 h groups were similar, whereas the 8 h group was relatively separated from the other groups. This suggests that the transcription profile of B. subtilis DB104 was altered, depending on the growth stage.

DEGs with the Highest Expression Levels and Fold Changes
To identify the most highly expressed genes, expression levels at all time points were calculated based on transcripts per millions mapped reads (TPM). The 10 most highly expressed DEGs are listed in Table 1. Genes encoding various spore coat proteins, such as cotX, cotY, and cotZ, showed high TPM values. Other extreme TPM values were observed for genes related to the spore of B. subtilis, such as gerE, safA, cgeA, and cgeB. Most of these genes were strongly expressed after 12 h. Among the genes listed in Table 1, only sdpC and skfA were highly expressed during exponential growth, and their expressions were down-regulated after 12 h. The DEGs with the highest fold-change values were also confirmed, and the top 10 of them are listed in Table 2. The most up-regulated genes were mostly related to sporulation and germination, such as the spore coat proteins cotD and cotG and the germination proteins yqfX and yqzD. In contrast, the most down-regulated genes were associated with a higher variety, such as response to stress (gsiB and gabT), chemotaxis (cheY), and metabolism (gabT, gbsB, gbsA, and xpt). Furthermore, the downregulated genes showed a much smaller fold-change value than the up-regulated genes. The highest expression changes did not appear after 18 h. Altogether, these results show that sporulation-related genes are highly expressed and regulated, and sharp changes in gene expression mainly occur up to the former part of sporulation (15 h) during B. subtilis DB104 growth. With a p-value < 0.05 and |fold change| ≥ 2 as the criteria, a total of 2122 DEGs were identified in all the comparison groups. The numbers of DEGs identified between the various time points during the growth of B. subtilis DB104 were analyzed using a Venn diagram (Figure 2c). The number of genes that were up-or down-regulated in each compared group is shown in Figure 2d. The results show that more transcriptome changes occurred during exponential growth to the former part of sporulation than the latter part and end of sporulation. In general, the number of up-regulated DEGs was higher than that of down-regulated DEGs.

DEGs with the Highest Expression Levels and Fold Changes
To identify the most highly expressed genes, expression levels at all time points were calculated based on transcripts per millions mapped reads (TPM). The 10 most highly expressed DEGs are listed in Table 1. Genes encoding various spore coat proteins, such as cotX, cotY, and cotZ, showed high TPM values. Other extreme TPM values were observed for genes related to the spore of B. subtilis, such as gerE, safA, cgeA, and cgeB. Most of these genes were strongly expressed after 12 h. Among the genes listed in Table 1, only sdpC and skfA were highly expressed during exponential growth, and their expressions were down-regulated after 12 h. The DEGs with the highest fold-change values were also confirmed, and the top 10 of them are listed in Table 2. The most up-regulated genes were mostly related to sporulation and germination, such as the spore coat proteins cotD and cotG and the germination proteins yqfX and yqzD. In contrast, the most down-regulated genes were associated with a higher variety, such as response to stress (gsiB and gabT), chemotaxis (cheY), and metabolism (gabT, gbsB, gbsA, and xpt). Furthermore, the downregulated genes showed a much smaller fold-change value than the up-regulated genes. The highest expression changes did not appear after 18 h. Altogether, these results show that sporulation-related genes are highly expressed and regulated, and sharp changes in gene expression mainly occur up to the former part of sporulation (15 h) during B. subtilis DB104 growth.

KEGG Pathway Enrichment Analysis
KEGG pathway analysis was performed to elucidate changes in pathways. DEGs in each compared group were divided into two groups consisting of up-or down-regulated cases and were analyzed separately. The total number of DEGs involved in each pathway was calculated. Table 3 shows the top 10 KEGG pathways with the number of DEGs involved. The DEGs were related to metabolic pathways (bsu01100), microbial metabolism Microorganisms 2023, 11, 1928 7 of 22 in diverse environments (bsu01120), biosynthesis of amino acids (bsu01230), and butanoate metabolism (bsu00650) in both cases, which suggests that metabolic pathway changes are significant throughout growth. In addition, it was noticed that the genes involved in ABC transporters (bsu02010) were down-regulated up to the beginning of sporulation (12 h) and up-regulated from the latter part of sporulation (18 h).

Cluster Analysis
Based on the Z-score, a K-mean cluster analysis was performed to further analyze the DEG expression patterns further. Ten distinct clusters were obtained, and their expression levels across growth were presented by a heatmap (Figure 3). Each cluster exhibits various expression patterns over time; in cluster 1, gene expression was decreased from 8 to 12 h, while it was increased thereafter until being equal to 8 h (n = 111 genes); in cluster 2, expression sharply decreased at 10 h and remained at that level (n = 177 genes); in cluster 3, gene expression increased slightly at 10 h and decreased with a nadir at 15 h, while it was slightly recovered thereafter (n = 133 genes); in cluster 4, the expression of the included genes was down-regulated to 12 h and then remained at a similar level (n = 265 genes); in cluster 5, the expression increased in a time-dependent manner with the zenith at 15 h (n = 336 genes); in cluster 6, 7, and 9, gene expression was up-regulated from 8 to 12 h, and then decreased at 15 h (n = 221, 330, and 116, respectively); in cluster 8, expression levels increased in a time-dependent manner with maximum expression at 18 h (n = 116 genes); and in cluster 10, genes were expressed in a zig-zag manner with the deepest parts at 12 and 18 h and the highest point at 15 h (n = 86 genes). The genes belonging to each cluster are listed in Supplementary File S3.

Functional Enrichment Analysis of Each Cluster
Gene ontology (GO) enrichment analysis was conducted to confirm the biological function of genes belonging to each cluster that showed different expression patterns over time. The genes in each cluster were classified into biological process (BP), molecular function (MF), and cellular component (CC). The GO terms that were significantly enriched and their relationships are shown in Tables S3-S11 (Supplementary File S2) and Figures S2-S10 (Supplementary File S1).

Functional Enrichment Analysis of Each Cluster
Gene ontology (GO) enrichment analysis was conducted to confirm the biological function of genes belonging to each cluster that showed different expression patterns over time. The genes in each cluster were classified into biological process (BP), molecular function (MF), and cellular component (CC). The GO terms that were significantly enriched and their relationships are shown in Tables S3-S11 (Supplementary File S2) and Figures S2-S10 (Supplementary File S1). The 177 genes in cluster 2 showed the high expression levels in the middle of exponential growth (8 h) but decreased expression from 10 h at the end of exponential growth ( Figure 3). DEGs belonging to this cluster were classified into 130 GO terms. The enriched GO terms are listed in Table S4 (Supplementary File S2), and their relationship is shown in Figure S2 (Supplementary File S1). Among them, the terms related to BP were the most common with 96, and mostly consisted of various metabolic processes (GO:0008152), such as ribonucleotide metabolic process (GO:0009259), amine catabolic process (GO:0009310), and organic acid metabolic process (GO:0006082). Terms, such as signaling (GO:0023052), secretion (GO:0046903), locomotion (GO:0040011), and localization (GO:0051179), were also found to be significantly enriched terms. Sixteen GO terms were related to CC and associated with the flagellum (GO:0019861). Four out of ten clusters showed a pattern where expression decreased from 15 h or 18 h (middle of sporulation). Among them, clusters 6 and 7 showed similar expression patterns, with expression at a low level in the exponential phase and then increased until sporulation began; however, differences in expression level were observed after 15 h. The 221 DEGs belonging to cluster 6 were classified into only 21 BP terms. The 330 DEGs belonging to cluster 7 were classified into 30 BP, 8 MF, and 2 CC terms. The enriched terms in each cluster are listed in Tables S8 and S9 (Supplementary File S2), and their relationships are shown in Figures S7 and S8 (Supplementary File S1), respectively. The BP-enriched DEGs of clusters 6 and 7 were mainly related to developmental process (GO:0032502), such as sporulation (GO:0043934) and cell differentiation (GO:0030154), and some other metabolic processes. The expression pattern of cluster 9 differed in that it had a high expression level from exponential growth, although it decreased from 15 h. As a result of annotating 116 DEGs belonging to cluster 9 into the GO database, there was no significant term that satisfied the p-value < 0.03, and the Benjamini-Hochberg-corrected p-value < 0.03. Unlike other clusters, DEGs of cluster 5 increased expression up to 15 h and then decreased from 18 h. The 336 DEGs belonging to cluster 5 were classified into nine BP terms and one CC term. The enriched terms are listed in Table S7 (Supplementary File S2), and their relationship is shown in Figure S6 (Supplementary File S1). BP terms were focused on developmental processes (GO:0032502), such as sporulation (GO:0043934) and cell differentiation (GO:0030154).

Genes That Are Up-Regulated in the Middle of Sporulation: Clusters 1 and 8
Genes belonging to cluster 1 were down-regulated at 12 h, but their expression increased from 15 h, the former part of the sporulation (Figure 3). The 111 DEGs of cluster 1 were classified into 53 GO terms. The enriched terms are listed in Table S3 (Supplementary File S2), and their relationship is shown in Figure S2 (Supplementary File S1). In cluster 1, most DEGs were enriched in BP terms, including metabolic processes such as tryptophan metabolic process (GO:0006568) and lipoteichoic acid biosynthetic process (GO:0070395) in addition to cell wall organization or biogenesis (GO:0071554), localization (GO:0051179), and small molecule metabolic process (GO:0044281). Of the seven MF terms, terms related to transporter activity (GO:0005215) were the most common, appearing in four. Terms enriched in this cluster were classified with relatively high p-values. In contrast, terms classified as enriched in cluster 8 had low p-values and were concentrated in developmental processes (GO:0032502). DEGs belonging to cluster 8 steadily increased expression from the beginning to the latter part of sporulation. The 347 genes belonging to this cluster were enriched in 15 BP terms. The enriched terms are listed in Table S10 (Supplementary File S2), and their relationship is expressed in Figure S9 (Supplementary File S1). In addition to the subterm of the developmental process (GO:0032502), these terms were also found to be significant in terms related to the organic cation transport (GO:0015695).

Genes That Are Transiently Turned on at the Former Part of Sporulation: Cluster 10
The genes belonging to cluster 10 were specifically turned on at 15 h (the former part of sporulation) and then down-regulated. The enriched terms are listed in Table S11 (Supplementary File S2), and their relationship is shown in Figure S10 (Supplementary File S1). The 86 DEGs in this cluster were enriched mainly for BP terms related to metabolic process, such as lipid catabolic process (GO:0016042) and organic acid catabolic process (GO:0016054), or monosaccharide transport (GO:0015749).

Sigma Factors
Various alternative sigma factors exist in B. subtilis and are transcribed by recognizing the promoter of the specialized gene set. By analyzing the expression patterns of genes regulated by each alternative sigma factor, their expression patterns were confirmed with clusters revealed through K-mean cluster analysis. The genes regulated by each sigma factor were classified using the DBTBS database [24]. σ A is a housekeeping sigma factor that regulates gene expression during the exponential phase and early sporulation in B. subtilis. The gene that encodes σ A , sigA, was expressed highly during whole growth, especially at 12 h, where sporulation started (Table 4). Genes controlled by σ A were found to be distributed across all clusters. In B. subtilis, σ B is a sigma factor that regulates various general stress genes. Clusters 6, 7, and 9, which have the highest expression at the transition point of 12 h, commonly contain many genes controlled by this sigma factor. However, the gene that encodes σ B , sigB, was expressed constantly during growth (Table 4). σ D is a sigma factor that regulates genes related to flagellar gene expression, chemotaxis, autolysis, and motility in B. subtilis. Cluster 2, whose expression was suppressed at the end of exponential growth, contained many genes controlled by σ D , including sigD. sigD (the gene for σ D ) is a component of the fla-che operon, which consists of more than 30 genes and can selfactivate. Many genes comprising this operon were included in cluster 2, including the genes encoding the rod protein flgC and the hook proteins flgE, cheA, cheB, and cheW, which comprise the flagellar of B. subtilis. In addition, flhO and flhP, which do not belong to this operon but encode the proteins that make up the flagellar rod, were also included in this cluster. Cluster 4, which exhibits a pattern where expression decreases as sporulation starts, also contained many genes controlled by σ D , similar to cluster 2. Among the fla-che operon genes, genes in cluster 4 included the genes encoding the hook proteins flgL and flgK and the rod protein flgB. In addition, the genes encoding flagellar assembly protein genes fliS and fliT, the flagellin protein hag, and the flagellar stators motA and motB are not part of the fla-che operon and are controlled by σ D , belonging to this cluster. Among the genes controlled by σ D , genes related to autolysis, such as lytB, lytC, and lytF, were included in cluster 4. Genes regulated by the sigma factor associated with sporulation were included in the cluster with the highest expression at 15 h or 18 h. The genes controlled by σ E and σ F , which are early sporulation-related sigma factors, were included in cluster 7. The DEGs controlled by σ E that regulate early sporulation genes in mother cells were also included in cluster 6. In contrast, genes controlled by σ K , a sigma factor related to late sporulation gene expression in the forespore, and σ G , sigma factor related to late sporulation gene expression in the mother cell, were predominantly included in clusters 8 and 5, respectively. The genes that encode early sporulation-specific sigma factors, sigE and sigF, were up-regulated at 10 h (Table 4). However, the genes that encode late sporulation-specific sigma factor, sigG, were up-regulated from 10 h. The other late sporulation-specific sigma factor σ K , which consists of two truncated genes, sigKc and sigKn [25], are highly expressed from 12 h. General stress protein genes, such as yceD, yceE, yceF, and yceG, are controlled by σ W , an extra-cytoplasmic function sigma factor and were included in clusters 4 or 9. The gene for σ W , sigW, was also included in cluster 9, and its expression was down-regulated from 15 h (Table 4). Among the regulons of σ L , known to be involved in degradative enzyme expression, only the acetoin dehydrogenase operon genes (acoABCL) were classified as cluster 8, and the rest did not belong to any cluster. The gene that encodes σ L , sigL, was up-regulated until 10 h and then down-regulated. Several other factors, such as σ I , σ M , σ X , σ Y , and YlaA, along with the sigma factors described above, regulate gene expression in B. subtilis. The dhb operon genes (dhbACEBF and mbtH) are controlled by σ I and were included in cluster 9. Some flavodoxin genes (fldN, ykuO, and fldP), which are controlled by the same sigma factor, were included in cluster 10. σ X regulon has dlt operon genes, such as dltABCDE, belonging to cluster 1 and the efe operon genes (efeUBO) and yce operon genes (yceDEFGH) belonging to cluster 4. In addition, some of the genes controlled by σ X belong to cluster 7 or cluster 6. YlaC, known to regulate genes associated with oxidative stress, regulates the yla operon to which it belongs, and genes belonging to this operon (ylaACD) were included in cluster 7 (the ylaB gene in yla operon was not mapped in this study).

Flagella
Among the movement of bacteria, swimming and swarming are the most representative. B. subtilis, a flagellated bacterium, can swim in liquids and swarm on solid surfaces through flagella. Therefore, the analysis of the flagellar gene could help to understand the mobility of this strain. The expression values of some flagellar structure genes, such as fliO, fliP, flgB, fliD, and flgD, are shown in Figure 4a-c and Table S12 (Supplementary File S2), which were classified as cluster 2 or 4 through K-mean clustering analysis. These genes remained down-regulated from the end of exponential growth (10 h) or the transition point (12 h). On the other hand, yvzB showed a different expression pattern from other flagellar genes, maintaining its expression level until 12 h, when sporulation began, and then down-regulated from the former part of sporulation (15 h

Flagella
Among the movement of bacteria, swimming and swarming are the most representative. B. subtilis, a flagellated bacterium, can swim in liquids and swarm on solid surfaces through flagella. Therefore, the analysis of the flagellar gene could help to understand the mobility of this strain. The expression values of some flagellar structure genes, such as fliO, fliP, flgB, fliD, and flgD, are shown in Figure 4a-c and Table S12 (Supplementary File S2), which were classified as cluster 2 or 4 through K-mean clustering analysis. These genes remained down-regulated from the end of exponential growth (10 h) or the transition point (12 h). On the other hand, yvzB showed a different expression pattern from other flagellar genes, maintaining its expression level until 12 h, when sporulation began, and then down-regulated from the former part of sporulation (15 h).

Sporulation and Spore Coat
The expression patterns of the gene involved in each sporulation stage are shown in Figure 5. In stage 0, cannibalism delays entering the sporulation process (in the next section). Most genes related to sporulation stage 0 and I were up-regulated at the exponential growth (8 h and 10 h). When septum is formed, many genes expressed in stage II were mostly up-regulated at 10 h (the end of exponential growth). Most of the genes expressed when stage III is completed also showed similar expression patterns to those of stage II but remained at a highly expressed level until 12 h (the beginning of sporulation). Most genes expressed during stages VI and V, in which the cortex and spore coat are formed and assembled, were up-regulated after 12 h, and they remained consistently high until 24 h when sporulation was completed.

Sporulation and Spore Coat
The expression patterns of the gene involved in each sporulation stage are shown in Figure 5. In stage 0, cannibalism delays entering the sporulation process (in the next section). Most genes related to sporulation stage 0 and I were up-regulated at the exponential growth (8 h and 10 h). When septum is formed, many genes expressed in stage II were mostly up-regulated at 10 h (the end of exponential growth). Most of the genes expressed when stage III is completed also showed similar expression patterns to those of stage II but remained at a highly expressed level until 12 h (the beginning of sporulation). Most genes expressed during stages VI and V, in which the cortex and spore coat are formed and assembled, were up-regulated after 12 h, and they remained consistently high until 24 h when sporulation was completed. After spore engulfment was completed, σ E first induced the expression of the spore coat proteins in the mother cell. The expression value of representative spore coat proteins of B. subtilis is shown in Figure 6 and Table S12 (Supplementary File S2). The σ E -induced spore coat protein genes (cotE, cotJA, cotJB, and cotJC) and spore morphogenic genes (spoIVA and spoVID) became up-regulated from 10 h at the end of exponential growth ( Figure 6; the cotE gene was not mapped in this study). Afterwards, most of the spore coat protein genes, such as cotA, cotD, cotF, cotH, cotM, cotT, cotV, cotW, cotY, cotZ, cotB, cotG, cotS, and cotX, and the spore coat gene transcription regulator gene gerE were up-regulated by σ K from 15 h, the former part of sporulation ( Figure 6). After spore engulfment was completed, σ E first induced the expression of the spore coat proteins in the mother cell. The expression value of representative spore coat proteins of B. subtilis is shown in Figure 6 and Table S12 (Supplementary File S2). The σ E -induced spore coat protein genes (cotE, cotJA, cotJB, and cotJC) and spore morphogenic genes (spoIVA and spoVID) became up-regulated from 10 h at the end of exponential growth ( Figure 6; the cotE gene was not mapped in this study). Afterwards, most of the spore coat protein genes, such as cotA, cotD, cotF, cotH, cotM, cotT, cotV, cotW, cotY, cotZ, cotB, cotG, cotS, and cotX, and the spore coat gene transcription regulator gene gerE were up-regulated by σ K from 15 h, the former part of sporulation ( Figure 6).
Microorganisms 2023, 11,1928 14 of 23 Figure 6. Expression patterns of the genes involved in spore coat expression. Gene expression is visualized as the Z-score. The red and blue colors represent the up-regulated and down-regulated gene, respectively. E represents the σ E -dependent genes. K represents the σ K -dependent genes.

Cannibalism
At the beginning of sporulation in B. subtilis (stage 0), cannibalistic behavior occurs to delay entering this energy-consuming process. The skf (sporulation killing factor) and sdp (sporulation delay protein) operons become involved in cannibalism. Among the skf operon genes, skfA, which encodes the sporulation killing factor itself, was strongly expressed in 10 h at the end of exponential growth and 12 h at the beginning of sporulation then decreased (Figure 4d). However, except for skfA, the remaining skf operon genes showed a Figure 6. Expression patterns of the genes involved in spore coat expression. Gene expression is visualized as the Z-score. The red and blue colors represent the up-regulated and down-regulated gene, respectively. E represents the σ E -dependent genes. K represents the σ K -dependent genes.

Cannibalism
At the beginning of sporulation in B. subtilis (stage 0), cannibalistic behavior occurs to delay entering this energy-consuming process. The skf (sporulation killing factor) and sdp (sporulation delay protein) operons become involved in cannibalism. Among the skf operon genes, skfA, which encodes the sporulation killing factor itself, was strongly expressed in 10 h at the end of exponential growth and 12 h at the beginning of sporulation then decreased (Figure 4d). However, except for skfA, the remaining skf operon genes showed a different pattern in which expression decreased from 12 h. In addition, the expression of those genes was significantly lower than that of skfA (Supplementary File S2: Table S13). Among the sdpABC operon, only the sdpC gene belonged to DEGs. Its expression was highest at 10 h and decreased sharply from 12 h when sporulation began (Figure 4e). Conversely, the expression of the sdpRI immunity operon continued to be up-regulated from 12 h to 24 h. The sdpA and sdpB genes of the sdpABC operon and the sdpI and sdpR genes of the sdpRI operon showed very-low expression values compared to sdpC (Supplementary File S2: Table S13).

RNA-Seq Data Validation by RT-qPCR
To validate the accuracy of transcriptome data, RT-qPCR analysis was conducted. Table 2 shows the top 10 up-and down-regulated DEGs according to fold change. Therefore, eight genes included in Table 2 (cotD, spoIIID, yqfX, spmH, gsiB, frlO, and cheY) were selected for validation by RT-qPCR. RNA samples from six-time points were used as templates. The results are shown in Figure 7, indicating that the data from RNA-seq and RT-qPCR analysis were consistent.

Discussion
This study investigated gene expression changes in B. subtilis DB104 during growth in a LB complete medium through RNA-seq. The simultaneous application of several

Discussion
This study investigated gene expression changes in B. subtilis DB104 during growth in a LB complete medium through RNA-seq. The simultaneous application of several analysis methods (DEG analysis, KEGG analysis, K-mean cluster analysis, GO analysis, and sigma factor function analysis) allowed the expression patterns of genes and their function to be revealed. We described gene expression changes using RNA-seq analysis, providing genome-wider information compared to using microarray analysis [6,7].
When bacterial growth becomes limited due to nutrient depletion, exponential growth ceases and the bacteria enter the stationary phase. The period between these two phases is called a transition point [26,27]. This transition point has been also observed in B. subtilis, and various metabolic pathway changes occur here [6,7]. In the present study, 12 h was found to be the transition point at which exponential growth ceases and the stationary phase starts (Figure 1). In addition, most of the highest fold changes in DEGs were found at this time point (Table 2). In addition, KEEG analysis revealed that alterations occurred in metabolic pathways around the transition point. The metabolic pathway (bsu01100) was up-regulated until 12 h, and some of the DEGs involved in this pathway started to become down-regulated from 15 h (Table 3). Previous studies showed that the metabolic rate decreased when the stationary phase started in the growth of B. subtilis [28,29]. The other enriched metabolic pathways, such as microbial metabolism in diverse environments (bsu01110) and the biosynthesis of secondary metabolites (bsu01110), were up-regulated to 12 h.
ABC transporters are responsible for transporting a wide range of substrates, such as amino acids, sugars, and organic acids, across the cell membrane through ATP hydrolysis and can perform various functions. Among the pathways enriched by KEGG analysis, only the DEGs involved in ABC transporters (bsu02010) were down-regulated before the transition point and subsequently up-regulated. A possible reason for this may be that a specific ABC transporter is required to initiate sporulation. The ABC transporter gene ftsEX, whose transcription levels were down-regulated at 10 h in this study (Supplementary  File S2: Table S13), plays a role in spore formation initiation [30]. In addition, the other ABC transporter genes oppA, oppB, oppC, oppD, and oppF, which are also involved in the initiation of sporulation [31,32], were down-regulated at 12 h (Supplementary File S2: Table S13). The fact that this system is related to stress responses, such as osmotic stress, could be another reason. When bacteria reach the stationary phase, they undergo various metabolism changes and become resistant to several stresses [7,33]. Corresponding to this, the genes for the osmoprotectant ABC transporters opuBA, opuBB, and opuBC were up-regulated at 18 h but were down-regulated at 12 h.
Diverse transcription regulators exist in B. subtilis. Among them, various alternative sigma factors, which allow RNA polymerase to recognize the promoter sequence rapidly and specifically, are essential for sporulation, stress response, chemotaxis, and motility [34]. Therefore, the sigma factors are a significant tool for understanding the gene expression of B. subtilis. To identify the different expression patterns of DEGs and their functions, K-mean clustering was coupled with GO and sigma factor analysis. In this study, the DEGs were divided into 10 clusters whose expression profiles could be redivided into five patterns in accordance with growth stage. Because gene expression was examined according to six-time points, the 10 clusters were divided again, as described above (Section 3.6).
The first group of clusters, which contains only cluster 2, has a pattern specifically expressed only in the middle of exponential growth. The DEGs belonging to cluster 2 were enriched in most terms. The relationship between them was also complicated, as shown in Figure S3 (Supplementary File S1), indicating that genes for maintaining various biological functions are highly expressed in exponential growth and subsequently turned off. Genes transcribed by σ A -or σ D -dependent RNA polymerase were the most prevalent among the 177 DEGs in cluster 2. B. subtilis σ A , similar to its E. coli counterpart σ 70 , is the primary sigma factor involved in the transcription of most housekeeping genes [35]. This sigma factor also participates in the transcription of the early sporulation genes [36,37] and stress-inducible genes [38,39]. It has been found that the expression of sigA (gene for σ A ) remains constantly high during the whole growth of B. subtilis DB104 (Table 4). In addition, the primary sigma factor, σ A in B. subtilis, usually transcribes genes with other sigma factors [40,41]. This could be why σ A -dependent genes were evenly distributed across all clusters in the transcriptome data. Many genes included in clusters 1 to 10 were transcribed by σ A -dependent promoter. However, this result does not correspond with that of a previous study, which found that the expression of σ A -dependent genes is turned off in the stationary phase [6]. In contrast, Blom et al. [7] showed results consistent with this research and observed that the SigA regulon is overrepresented in the exponential growth phase but could be expressed after the transition point. Furthermore, results from this study suggest that σ A -dependent genes have various expression patterns because σ A works both alone and jointly with other sigma factors. B. subtilis σ D controls genes related to flagellar [42], autolysis [43], chemotaxis, and motility [44]. σ D -dependent genes were included predominantly in cluster 2 in addition to cluster 4 in the second cluster group (genes that are down-regulated from the beginning of sporulation). With corresponding data, the GO term "flagellum (GO:0019861)" was enriched in clusters 2 and 4. As previously stated, yvzB, the flagella structure gene, began to decline after the transition point, but remained highly expressed until sporulation began (Figure 4a-c). This may be because this gene is a putative filament gene, and it is unknown whether it contains the proteins that make up the filament structure. Mukherjee and Kearns [45] suggested that flagella synthesized in a mother cell could function in a granddaughter cell. The results from this study showed that flagellar structural genes are down-regulated during the stationary phase in B. subtilis DB104, but further research is needed to determine whether this weakens the swarming/swimming mobility. Meanwhile, among the genes controlled by σ D , autolysis-related genes were included in cluster 4 rather than cluster 2, indicating that autolysin expression continued later than that of flagella proteins. The sigD (gene for σ D ) was also included in cluster 2, and the expression level was down-regulated from 10 h (Table 4). This result indicates that σ D produced in the middle of exponential growth remains and functions afterward. The DEGs belonging to clusters 3 and 4 showed a high expression level in exponential growth, similar to genes belonging to cluster 2. However, the difference is that they maintained a high expression afterward and then showed a low expression from the beginning of sporulation. They were also classified as diverse and complex terms (Tables S5-S7 in Supplementary File S2, Figures S4 and S5 in Supplementary File S1). Therefore, these results show that genes for maintaining various biological functions are highly expressed during exponential growth, and some are upregulated even after exponential growth is over.
An extracellular RNA polymerase sigma factor of B. subtilis, σ W , is involved in the response to cell envelope stress, such as to antimicrobial compounds [46,47]. σ B is the general stress sigma factor of B. subtilis and is also involved in the transcription of general stress proteins induced by heat, salt, and oxygen stress [48]. σ W -dependent genes were distributed in clusters 4 and 9, and σ B -dependent genes were distributed in clusters 6, 7, and 9. Except for cluster 4, these clusters belong to the third group of clusters, which contains clusters 5, 6, 7, and 9 (genes that are turned on at the beginning of sporulation and then turned off in the middle). The genes for these stress-related sigma factor, sigB and sigW, were also up-regulated until 12 h and then down-regulated from 15 h ( Table 4). The time point 12 h, where sporulation starts, is also a transition point in this study (Figure 1), and it was found that, if the growth of B. subtilis enters the transition point, stress responses may occur, even if the stresses are not induced from outside factors [7,49]. This study supports this view, showing that the σ W -or σ B -dependent genes that respond to stress are primarily turned on at the beginning of sporulation. Interestingly, cluster 4, whose expression was decreased from the transition point, also contains σ W -dependent genes: iron transporter genes from the efeUOB operon [50], general stress genes from the yceCDEFGH operon [51], and phage shock protein gene and unknown function genes from the pspA-ydjGHI operon [52]. It can be assumed that these genes have different expression patterns among the stress response genes because of the growth conditions used in the study. The results show that a spontaneous stress response may occur at the transition point, but even a general stress gene may not be induced at this time.
For the spore-forming bacteria B. subtilis, nutrition limitation or high population density can trigger the complex process of sporulation [53]. During the growth of B. subtilis in this study, sporulation started from 12 h and progressed rapidly up to 18 h and then proceeded to 90% in 24 h ( Figure 1). As expected, the GO term "sporulation (GO:0043934)" was enriched in clusters that belong to the third group, clusters 5, 6, and 7. Furthermore, in the GO analysis of these clusters, sporulation was significantly enriched with a relatively low p-value. Consistent with this, σ E , σ F , σ G , and σ K , which are related to gene expression during sporulation, are also included in clusters 6 and 7. The sporulation of this strain could be described as eight stages (stages 0-VII) based on morphological features [54,55]. In this study, we compared these sporulation stages with transcriptomic expression changes and analyzed them with the growth stage ( Figure 5). Stages 0 and I are vegetative growth and axial filamentation, respectively. However, they are no longer distinguishable because of any mutant arrests at this stage [55], so we considered them together. In these stages, the cells start the sporulation with cannibalistic behavior (see the discussion below). The master regulator Spo0A is essential in transitioning from vegetative growth to the sporulation stage. The activation via the phosphorylation process from KinA/KinB to Spo0F, Spo0B, and Spo0A of the spo0A gene transcribed by the σ H -dependent promoter, the first activated sigma factor, signals the start of sporulation. In addition, phosphatase genes, such as the Rap family (rapA, B, E, and H) [56] and Spo0E [57], the cell division protein gene racA that helps chromosomes properly segregate [58], and the genes like sirA [59] and sda [60] that encode proteins for regulating the appropriate chromosome number, are expressed at these stages. All their expression levels were high at 8 h and/or 10 h, indicating that stages 0-I correspond to those before spores are formed. The activated Spo0A, Spo0A~P, regulates the spoIIA operon, which starts stage II. Genes, such as phosphatase gene spoIIE that dephosphorylates SpoIIAA [61], cell-division protein ftsZ that is required for septum formation [62], mother cell-specific sigma factor gene sigE, and protease gene spoIIGA that activates Pro-SigE [63], are expressed in stage II, including genes of spoIIA operon (anti-anti-sigma factor gene spoIIAA, anti-sigma factor gene spoIIAB, and forespore-specific sigma factor gene sigF). All of these genes expressed highly at 12 h, showing that stage II, at which asymmetric division occurs, matched the end of exponential growth. At stage III, the cell-wall hydrolase genes spoIID and spoIIP that are required for engulfment [64,65] are expressed. The membrane fission protein gene fisB, the spoIIQ gene, and the spoIIIA operon gene (spoIIIAA-spoIIIAH) that are required for the activation of σ G [66] are also transcripted at this stage, where engulfment is processed. Their up-regulated expressed level from 10 h to 12 h indicates that stage III spans the end of exponential growth and the beginning of sporulation. Cortex and spore coat are produced and assembled at stages VI and V. Genes, such as spore coat morphogenetic protein gene spoIVA, DNA recombinase gene spoIVCA that is required for σ K gene, and mother cell-specific sigma factor gene sigKc/sigKn, are related to cortex and spore coat. Most of these genes were highly expressed after 15 h (the formal part of sporulation), but some were up-regulated from 10 h (end of exponential growth). This result shows that the genes involved in stages VI and V could be expressed earlier. So, we decided to discuss spore coat genes in depth.
σ E and σ F are the sigma factors related to early sporulation genes in the mother cell and forespore, respectively. At the same time, σ K and σ G are the sigma factors related to late sporulation genes in the mother cell and forespore, respectively. Additionally, many σ K -dependent genes were also included in cluster 8 of the fourth cluster group (genes that are up-regulated in the middle of sporulation). Notably, the expression of spore coat genes was clearly distinguished according to which sigma-factor-dependent RNA polymerase transcribes ( Figure 6). The B. subtilis spore coat consists of two layers, an inner coat and outer coat. According to the sporulation process, the spore coat proteins are sequentially expressed in the mother cell [67]; however, it is not completely clear how many proteins make up the spore coat and how they are assembled. However, the results show that, no matter where the proteins exist, their transcription is affected only by σ E or σ K , and there is no difference in the expression time point between them.
Cluster 8, where the GO term "acetoin metabolic process (GO:0045149)" was enriched (Supplementary File S2: Table S10), included some of the genes that are controlled by σ L , the acetoin dehydrogenase operon genes acoABCL. σ L is known to be involved in degradative enzyme gene expression [68]. The acetoin is one of the extracellular products of B. subtilis [69]. Since it is reused as a carbon source after the cells enter to stationary phase and sporulation [70,71], the operon genes are up-regulated after 15 h. On the other hand, other genes belonging to σ L regulon were not classified into any cluster, and even sigL (gene for σ L ) was down-regulated from 12 h (Table 4). Similar to σ L , the regulons of other alternative sigma factors in B. subtilis, such as σ I , σ M , σ X , σ Y , and YlaA, did not show a consistent expression pattern (Supplementary File S1: Figure S11). Moreover, not only were the TPM values of sigI, sigM, sigX, sigY, and ylaA (gene for σ I , σ M , σ X , σ Y , and YlaA, respectively) much lower than those of other sigma factors (Table 4), there were also fewer genes belonging to their regulons. These results suggest that these sigma factors are minor compared to the other sigma factors described above (σ A , σ B , σ D , σ E , σ F , σ G , σ K , and σ W ), and the role of other transcription regulators is essential for the regulation of the genes controlled by them.
Of particular interest, GO terms related to the lipid catabolic process were enriched only in cluster 10, which has a unique expression pattern. From these results, the genes related to the lipid catabolism of B. subtilis are expressed specifically throughout the entire sporulation when grown in a rich medium. The programmed cell death (PCD) system could be the reason for this specific expression. In B. subtilis, as in other eukaryotes or prokaryotes, a PCD system exists to delay sporulation, an energy-consuming and irreversible process [72]. The step-by-step expression of the gene involved in this system was well identified in our data. As shown in Figure 4c,d, the expression of skfA and sdpC, which play an essential role in the PCD system, peaked at 12 h and 10 h, respectively. In the Spo0A-active cell, genes encoded in the skf operon are expressed, and killing factors are produced and exported [73,74]. In previous transcriptome analysis conducted in LB with glucose [6], where sporulation was suppressed, there was no skfA expression changes during B. subtilis growth, indicating the PCD system occurs with spore formation. On the other hand, in both the Spo0A-active and -inactive cells, the sdp operon is turned on. SdpC (product of sdpC), a constituent protein of the sdp operon, increases the expression of the fadNAE operon, which encodes lipid catabolism enzymes through SdpR, which delays sporulation [75]. The sdpR and fadNAE operon gene expressions peaked at 12 h and 15 h, respectively (Supplementary File S2: Table S13). Therefore, it can be assumed that regulating the lipid catabolic process through the cascade of this sequential expression prevents sporulation triggers caused by a lack of energy.

Conclusions
This transcriptome study provided a novel insight into the dynamic transcriptional response of each gene in B. subtilis DB104 during growth. The highest expressed genes and fold changes were identified around the transition point. Furthermore, the gene expression patterns were identified and DEG functions of were established through Kmean cluster analysis, GO enrichment analysis, the roles of alternative sigma factors, and their simultaneous application. The spore formation process is the most regulated in B. subtilis DB104, and most transcriptomic and metabolic changes occur with the entry of sporulation, named "the transition point". Because sporulation is an essential program in this organism, many processes during growth, such as metabolism changes, stress-induced responses, and cannibalism, overlap with sporulation. The gene expression involved in each sporulation stage was compared and matched with a growth curve. Additionally, this research showed that the timing of expression of the spore coat gene is determined only by the sigma factor, regardless of its location. Our results also suggest that, among the sigma factors of B. subtilis, σ A , σ B , σ D , σ E , σ F , σ G , σ K , and σ W play significant roles. Knowledge of transcriptional changes during the growth of B. subtilis DB104 will help to elucidate the expression and regulatory mechanisms of genes under various growth conditions and for the industrial application of this strain.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/microorganisms11081928/s1. Supplementary File S1. Figure S1: Observation of sporulation by optical microscopy at each time point, Figure S2: Gene ontology (GO) term networks of differentially expressed genes involved in cluster 1, Figure S3: Gene ontology (GO) term networks of differentially expressed genes involved in cluster 2, Figure S4: Gene ontology (GO) term networks of differentially expressed genes involved in cluster 3, Figure S5: Gene ontology (GO) term networks of differentially expressed genes involved in cluster 4, Figure S6: Gene ontology (GO) term networks of differentially expressed genes involved in cluster 5, Figure S7: Gene ontology (GO) term networks of differentially expressed genes involved in cluster 6, Figure S8: Gene ontology (GO) term networks of differentially expressed genes involved in cluster 7, Figure S9: Gene ontology (GO) term networks of differentially expressed genes involved in cluster 8, Figure S10: Gene ontology (GO) term networks of differentially expressed genes involved in cluster 10, Figure S11: Expression patterns of the genes included in each regulon. Supplementary File S2; Table S1: Summary of RNA-seq and the reads mapped to the genome of Bacillus subtilis subsp. subtilis str. 168, Table S2: Primers used for RT-qPCR analysis. Table S3: Gene ontology (GO) enrichment analysis of cluster 1, Table S4: Gene ontology (GO) enrichment analysis of cluster 2, Table S5: Gene ontology (GO) enrichment analysis of cluster 3, Table S6: Gene ontology (GO) enrichment analysis of cluster 4, Table S7: Gene ontology (GO) enrichment analysis of cluster 5, Table S8: Gene ontology (GO) enrichment analysis of cluster 6, Table S9: Gene ontology (GO) enrichment analysis of cluster 7, Table S10: Gene ontology (GO) enrichment analysis of cluster 8, Table S11: Gene ontology (GO) enrichment analysis of cluster 10, Table S12: Transcripts per millions mapped reads (TPM) values of DEGs involved in each structure, Table S13

Data Availability Statement:
The datasets supporting the conclusions of this article are included within the article as well as their additional files. The raw sequence data have been submitted to the Sequence Read Archive (SRA, https://www.ncbi.nlm.nih.gov/sra, accessed on 13 February 2023) under Bioproject accession number PRJNA934277 and reference Biosample accession number SAMN33269663-80.

Conflicts of Interest:
The authors declare no conflict to interest.