Skip to main content

Transcriptomic analysis of Myxococcus xanthus csgA, fruA, and mrpC mutants reveals extensive and diverse roles of key regulators in the multicellular developmental process

Abstract

Background

The bacterium Myxococcus xanthus provides an important multicellular model for understanding stress responses. The regulatory proteins CsgA, FruA, and MrpC are essential to survive prolonged starvation by forming fruiting bodies, which are mounds containing hardy round spores formed from vegetative rods, but the genome-wide pathways affected by these proteins remain poorly understood. Only a fruA mutant transcriptome and MrpC ChIP-seq have been reported. We describe RNA-seq transcriptome analysis of csgA, fruA, and mrpC mutants relative to a wild-type laboratory strain midway during the starvation-induced developmental process, when mounds, but not spores, have formed.

Results

We show that CsgA, FruA, and MrpC broadly impact developmental gene expression, with over 60% of the genes differentially expressed in one or more mutants. Building upon previous investigations, we found that strongly regulated genes in the mrpC mutant correlate with MrpC DNA-binding sites located ~ 80 bp upstream of transcriptional start sites. We also confirmed that FruA directly or indirectly regulates many genes negatively, as well as many others positively. CsgA regulates indirectly and its strongest effects are positive. MrpC strongly stimulates fruA transcription and FruA accumulation, impacting many genes, but our results reveal that MrpC is also a strong negative or positive regulator of hundreds of genes independently of FruA. Indeed, we observed nearly every possible pattern of coregulation, unique regulation, and counterregulation by comparing the wild-type and mutant transcriptomes, indicating diverse roles of CsgA, FruA, and MrpC in the developmental gene regulatory network. The genes most strongly regulated were coregulated in two or three of the mutants. Each set of genes exhibiting differential expression in one or more mutants was analyzed for enrichment of gene ontology (GO) terms or KEGG pathways, and predicted protein-protein interactions. These analyses highlighted enrichment of pathways involved in cellular signaling, protein synthesis, energetics, and envelope function. In particular, we describe how CsgA, FruA, and MrpC control production of ribosomes, lipid signals, and peptidoglycan intermediates during development.

Conclusions

By comparing wild-type and mutant transcriptomes midway in development, this study documents individual and coordinate regulation of crucial pathways by CsgA, FruA, and MrpC, providing a valuable resource for future investigations.

Peer Review reports

Background

Myxobacteria are a diverse group of gram-negative bacteria found ubiquitously in soil. They exhibit social characteristics and prey en masse on other soil microbes and organic detritus, digesting them by secreting a cocktail of lytic enzymes [1]. They produce an array of bioactive secondary metabolites, including antimicrobials such as myxovirescin [2] and epothilones [3], the latter being a class of microtubule-targeting agents used in the treatment of breast cancer. Analysis of myxobacterial genomes reveals an abundance of genes predicted to code for proteins involved in signaling and transcriptional responses [4, 5]. In terms of gene regulation in this fascinating group of bacteria, much of our understanding is through studies of the developmental process of the model myxobacterium Myxococcus xanthus [6].

When nutrient sources are scarce, starving M. xanthus rod-shaped cells undergo development by producing signals, regulating over a thousand genes, coordinating their movements to build multicellular mounds, and differentiating into round spores to form fruiting bodies [1, 7, 8]. The spores are physically resilient and resistant to harsh environmental conditions, yet can germinate into new, motile rods when nutrients return. Interestingly, the majority of the starving rods undergo lysis, presumably to supply nutrients to the remaining cells [9, 10]. Other cells stay outside of fruiting bodies in a quiescent state seemingly similar to antibiotic persister cells [11, 12]. These “peripheral rods” may be part of a bet-hedging strategy, allowing growth if nutrients become available shortly, whereas spores can provide longer-term survival. For over half a century, scientists have sought to understand the signaling and gene regulatory network underpinning the control of cellular physiology, movement, and fate during M. xanthus fruiting body formation.

Early investigations identified mutants defective in producing extracellular signals important for development [13]. Ensuing work revealed three types of signal molecules and their roles. The A-signal is a mixture of peptides and amino acids produced by extracellular protease activity of M. xanthus rods shortly after the onset of starvation to measure whether the cell density is sufficient to proceed with the developmental process [14,15,16]. The E-signal is a combination of a branched-chain fatty acid and an ether lipid necessary for mound formation and sporulation, respectively [17,18,19,20]. The C-signal is a fragment of the CsgA protein [21, 22] and/or diacylglycerols produced by cardiolipin phospholipase activity of intact CsgA [23]. The C-signal acts later than the A-signal [24,25,26]. C-signaling is needed for both proper mound formation and sporulation [27,28,29]. Indeed, C-signaling appears to ensure that mounds form first, then spores, perhaps by gauging end-to-end contacts between rods as they align during mound building [30,31,32,33]. The alignment of rods enhances C-signaling [31] and this positive feedback loop is proposed to trigger sporulation [27]. Here, we examine the impact of a loss-of-function mutation in csgA on the M. xanthus transcriptome at 24 h poststarvation (PS), when mounds normally would have formed, but prior to sporulation [34].

Intracellular signals are also crucial for M. xanthus development. The concentrations of the second messengers guanosine penta- and tetra-phosphate [(p)ppGpp] [35, 36] and cyclic diguanylate (c-di-GMP) [37] increase in starving cells. Amino acid limitation leads to uncharged tRNAs bound to ribosomes, stimulating RelA to produce (p)ppGpp during the stringent response [38]. Ectopic production of (p)ppGpp is sufficient to turn on some genes required for fruiting body formation [39]. (p)ppGpp produced during the stringent response is necessary for the production of the extracellular A- and C-signals [40,41,42]. DmxB is a development-specific diguanylate cyclase that increases the c-di-GMP concentration, inducing genes for exopolysaccharide synthesis [37].

A signal-responsive gene regulatory network governs the expression of over a thousand genes during the developmental process of M. xanthus [6, 43, 44]. The gene regulatory network has interconnected cascades of signal-responsive transcription factors that impact the expression of target genes, whose products affect cellular metabolism, motility, and differentiation. Early in development, a cascade of enhancer-binding proteins (EBPs) expected to be phosphorylated by protein kinases in response to unknown starvation-dependent signals, regulates many target genes [45], including genes involved in (p)ppGpp accumulation [46, 47] and A- and C-signal production [48,49,50,51,52,53,54] (Fig. 1). The EBP cascade connects via MrpB ~ P to a second transcription factor cascade involving MrpC followed by FruA [50, 55, 56], which regulates additional target genes [6]. For example, MrpC was recently shown to directly repress dmxB and directly activate pmxA (a phosphodiesterase that degrades c-di-GMP) to decrease the c-di-GMP level during development [57].

Fig. 1
figure 1

Signal-responsive gene regulatory network leading to mound formation during M. xanthus development. A simplified model shows starvation initiating the enhancer-binding protein (EBP) cascade of transcription factors (blue boxes, and blue arrows show their effects). Each EBP is likely phosphorylated (~ P) by a protein kinase in response to a signal molecule. The EBP cascade controls production of the intracellular signal (p)ppGpp and the extracellular A- and C-signals, which positively feed back on Nla28 ~ P and ActB ~ P, respectively (black arrows). Nla28 ~ P activates transcription of the gene for MrpB ~ P, which in turn activates transcription of the gene for MrpC (orange boxes show Mrp proteins and orange arrows show their effects). MrpC positively feeds back on C-signal production and activates transcription of the gene for FruA. C-signaling and MXAN4899 ~ P appear to activate FruA posttranslationally, designated FruA* (green boxes show FruA/FruA* and green arrows show their effects). FruA* is depicted to positively regulate the fmg, dev, and fadIJ operons in combination with MrpC. FruA/FruA* and MrpC lead to mound formation by 24 h poststarvation. Alignment of cells during mound building feeds back positively on C-signaling. Addition of nutrients before or during mound formation causes rapid proteolysis of MrpC and halts development. See the text and [6] for references

Additional direct targets of the cascade involving MrpC followed by FruA include the fmg, dev, and fadIJ operons (Fig. 1). MrpC is a Crp/Fnr family DNA-binding protein whose expression is induced by the MrpAB two-component signal transduction system [55, 58], which in turn is induced by Nla28 ~ P [50], a part of the EBP cascade. MrpB ~ P is an EBP, but its cognate protein kinase is unknown. MrpA appears to act as a phosphatase of MrpB ~ P [58]. MrpC is a master transcriptional regulator of mound formation and sporulation [58,59,60], and serves as a checkpoint for developmental progression – addition of nutrients before or during the appearance of mounds causes rapid proteolysis of MrpC and blocks commitment to sporulation [34, 61] (Fig. 1). MrpC directly activates transcription of fruA [56]. FruA is similar to response regulators of two-component signal transduction systems [62], but FruA appears to be atypical since a cognate protein kinase has not been found, its receiver domain lacks residues normally important for phosphorylation, and addition of small-molecule phosphodonors does not change FruA DNA-binding activity [63]. C-signaling appears to activate FruA posttranslationally [64, 65] (Fig. 1, activated FruA is designated FruA*), but the mechanism of activation is unknown. FruA binds cooperatively with MrpC to the promoter regions of four fmg genes or operons, the dev operon, and the fadIJ operon [63, 66,67,68,69,70] (Fig. 1). Mutations that interfere with binding of FruA or MrpC in vitro decrease promoter activity in vivo, indicating that FruA or FruA*, together with MrpC, directly activate transcription of fmg and dev genes [71,72,73,74,75,76]. This regulatory paradigm was reinforced by ChIP-seq analysis of MrpC-binding sites at 18 h PS and the finding that 13 out of 15 DNA fragments with an MrpC-binding site were cooperatively bound by MrpC and FruA in vitro [77].

Recently, new paradigms involving regulation by MrpC and FruA emerged from a study of the exoA-I, exoL-P, and nfsA-H operons involved in spore coat biogenesis [70]. Transcript levels from these operons were greater in a fruA mutant than in a wild-type (WT) laboratory strain early in development, suggesting negative regulation by FruA. Later in development, FruA or FruA* appeared to positively regulate exoA-I and nfsA-H, but not exoL-P. A mutation in mrpC affected the transcript levels from each operon differently over developmental time. Furthermore, each promoter region was unique in terms of MrpC and FruA binding in vitro, and in terms of binding to a fusion protein containing the Nla6 DNA-binding domain, which presumably reflects Nla6 ~ P of the EBP cascade.

Beyond the complexity of the gene regulatory network revealed by gene-specific studies, DNA microarray experiments and more recently RNA-seq based transcriptomics have begun to show the remarkable rewiring at the heart of M. xanthus fruiting body formation. The first RNA-seq analysis followed development on nutrient-limited agar over a 96-h timecourse and identified 1415 genes expressed sequentially in 10 groups [43]. The authors highlighted changes in expression of genes for motility, glycogen and lipid body utilization as energy sources, gluconeogenesis and turnover of polysaccharides and proteins to provide precursors for macromolecular synthesis, secondary metabolite production, signal transduction, transcription, and translation.

A second RNA-seq study examined development on plastic coverslips in flow cells over a 72-h timecourse and used rRNA depletion [44]. These authors identified 1522 genes expressed in eight groups, including 512 genes found in the first study plus 1010 genes not observed previously. Together, the two studies identified 2425 developmental genes. Both studies identified signal transduction, and energy production and conversion, as major functional categories of regulated genes. Genes for cell wall/membrane/envelope biogenesis also emerged as a major category under developmental control in the second study. There were also new insights into regulation of the initiation of development and central metabolism.

A third RNA-seq study used submerged culture conditions in which cells form a biofilm on a plastic surface, as in the second study, but with incubation in static petri dishes rather than flow cells [78]. The timecourse was shorter, 24 h, and compared a WT laboratory strain with a fruA mutant. The analysis showed that > 4000 genes were up- or down-regulated in the WT strain and ~ 1350 genes differed in expression comparing the fruA and the WT strain at 24 h. In general, the results agreed with the two prior studies, showing upregulation of genes for signal transduction and secondary metabolite production during development of the WT strain, and downregulation of genes for translation and many metabolic processes. Interestingly, FruA positively regulated some genes for signal transduction and polyketide synthases involved in secondary metabolite production, but FruA negatively regulated other genes for signal transduction and some genes for fatty acid synthesis.

There remains much to discover about the breadth of developmental regulation by FruA, as well as by its transcriptional activator MrpC and its posttranslational activator CsgA. To begin filling this knowledge gap, we performed RNA-seq analysis of csgA, fruA, and mrpC mutants relative to a WT laboratory strain at 24 h PS, when mounds but not spores have formed under our submerged culture conditions [34]. We found that CsgA, FruA, and MrpC regulate thousands of genes, often in unison (coregulation) or individually (unique regulation), and much less often in opposition (counterregulation). Pathways involved in cellular signaling, protein synthesis, energetics, and envelope function are enriched among the genes regulated by CsgA, FruA, and MrpC. We focus on their regulation of three pathways – control of ribosome biogenesis by the stringent response to starvation, branched-chain amino acid (BCAA) metabolism important for developmental lipid signals, and peptidoglycan (PG) biosynthesis owing to the cell wall remodeling necessary during spore formation and germination.

Methods

Bacterial strains, growth, and development

The M. xanthus strains used in this study were the laboratory WT strain DK1622 [79], the csgA mutant DK5208 (csgA::Tn5-132 Ω205) [80], the fruA mutant DK5285 (fruA::Tn5 lac Ω4491) [25], and the mrpC mutant SW2808 (ΔmrpC) [58]. Three biological replicates were performed for each strain. M. xanthus were grown in CTTYE liquid medium (1% Casitone, 0.2% yeast extract, 10 mM Tris-HCl [pH 8.0], 1 mM KH2PO4-K2HPO4, 8 mM MgSO4 [final pH 7.6]) with shaking at 32 °C. CTT agar (1.5%) lacking yeast extract was used for growth on solid medium, supplemented with 40 µg/ml of kanamycin sulfate or 12.5 µg/ml of oxytetracycline as needed. Starvation-induced development was performed under submerged culture conditions [81] as described previously [34]. Briefly, M. xanthus from log-phase CTTYE cultures were collected by centrifugation and resuspended in MC7 (10 mM morpholinepropanesulfonic acid [MOPS, pH 7.0], 1 mM CaCl2) at a concentration of 1,000 Klett units. An aliquot (1.5 mL) of cell suspension plus 10.5 mL of MC7 was added to an 8.5-cm-diameter plastic petri dish and incubated at 32 °C. Cells adhere to the bottom of the dish, forming a biofilm. After 24 h, the MC7 overlay was removed by aspiration and replaced with a solution to inhibit RNase activity (0.5 mL 5% phenol [pH 7] in ethanol plus 4.5 mL of MC7 buffer). The biofilm was scraped from the bottom of the dish using a sterile cell scraper, aspirated into a 15-mL centrifuge tube, flash frozen in liquid nitrogen, and stored at -80 °C.

RNA extraction and library preparation

Frozen samples were thawed on ice and centrifuged at 8000 × g for 10 min. The supernatant was discarded and RNA was isolated from the cell pellet using the hot-phenol extraction method and digested with DNase I (Roche) as described previously [60]. RNA quality was determined using an Agilent Bioanalyzer system and RNA samples with an RNA Integrity Number (RIN) > 8 were subject to rRNA depletion as described previously [44]. Briefly, 2–3 µg of RNA was mixed in hybridization buffer (0.1 M Tris-HCl [pH 7.0], 0.2 M NaCl) with 76 DNA oligonucleotides (0.17 µM/oligo) that anneal to M. xanthus 16 S and 23 S rRNAs, the mixture was incubated in a thermocycler at 95 °C for 2 min, then ramped down to 45 °C at a rate of 0.1 °C/s, followed by additional incubation at 45 °C for 5 min. The annealed samples were treated with Hybridase™ Thermostable RNase H (Epicentre) in order to hydrolyze RNA-DNA hybrids and digested with DNase I (Roche) to remove remaining oligos. Agencourt RNAClean XP beads (Beckman Coulter) were used for cleanup of the rRNA-depleted RNA samples. qPCR verification was performed to confirm successful (i.e., ~ 90%) rRNA depletion. The KAPA RNA HyperPrep kit (KAPA Biosystems) was used as described in the product manual to generate adapter-ligated libraries for multiplex sequencing. The Qubit 1X dsDNA HS (High-Sensitivity) Assay Kit (ThermoFisher Scientific) was used to measure library concentration. Libraries (20 ng each) were pooled and sequenced in one lane of an Illumina HISEQ 4000 Rapid Run flow cell in the 50 bp, single-end read format. The library from one biological replicate of the fruA mutant yielded only ~ 0.05–0.1% as much sequence as each of the other libraries, so those sequences were not included in our analyses.

RNA-seq transcriptome and differential gene expression analyses

The SPARTA (Simple Program for Automated reference-based bacterial RNA-seq Transcriptome Analysis) workflow [82] was used to analyze raw data from the HISEQ read output. SPARTA is turnkey software for analysis of RNA-seq data sets that conducts read trimming and adapter removal with Trimmomatic, performs quality analysis of the data sets with FastQC, maps the reads to the reference with Bowtie, and counts transcript or gene feature abundance with HTSeq. In brief, raw data files were subjected to trimming of low-quality bases and removal of adapter sequences using Trimmomatic 0.30 [83] with a 4 bp sliding window, cutting when the read quality dropped below 15 or read length was less than 36 bp. Trimmed reads were aligned to the M. xanthus DK1622 genome (GenBank: CP000113.1) using Bowtie [84]. Aligned reads were then counted per gene feature in the M. xanthus DK1622 genome using the HTSeq software suite [85]. Data were normalized by estimating effective library sizes using robust regression within the DESeq package [86]. Statistical analysis and differential gene expression was performed in RStudio [87] by fitting a negative binomial model to each set of conditions and testing for differences utilizing the DESeq package. Differentially expressed genes were identified as genes with differential gene expression ≥ 2-fold, and a False Discovery Rate (FDR)-adjusted p-value ≤ 0.05. Venn diagrams were generated using the eulerr web tool (https://eulerr.co/).

Bioinformatic analyses

Heatmaps were prepared with the pheatmap package in R [88]. Horizontal bar and lollipop graphs were prepared in R with the ggplot2 package [89]. The horizontal stacked bar graph was prepared in Microsoft Excel. Comparative analysis of MrpC ChIP-seq [77] and Cappable-Seq [57] findings were performed with custom R scripts, with combined histogram and density plots prepared using the ggplot2 and ggExtra libraries in R. Briefly, for each gene with an identified ChIP-seq MrpC peak, the peak position was compared to Cappable-Seq transcriptional start sites (TSSs), keeping all TSSs within 200 bp of the ChIP-Seq peak. Each gene was assigned to a TSS based on the highest count abundance. Scatterplots were prepared using the matplotlib library in Python [90]. STRING cluster analysis and visualization used Cytoscape 3.10.1 [91] with the stringApp package. GO term overrepresentation analysis was performed with the ShinyGO 0.77 webserver [92] and the Database for Annotation, Visualization and Integrated Discovery (DAVID) [93, 94]. We used all 4610 genes that were differentially expressed in one or more mutants as the background for both methods. The output lists for ShinyGO and DAVID were compared and GO terms from DAVID which were > 95% in agreement with those from ShinyGo were used. Pathway enrichment analysis was performed using the clusterProfiler [95] and pathview [96] packages in R. Correlation analyses were prepared using the ggpubr library in R [97].

Results

RNA profiles show tight intra-strain clustering and distinct inter-strain differences

To identify overall differences in transcript levels between the strains tested and check the quality of biological replicate transcript counts, the pcaExplorer library in R was used to perform supervised principal component analysis [98]. Nearly all the variance across the strains could be explained by the first three principal components, with the vast majority by PC1 (81%) followed distantly by PC2 and PC3 (9.5% and 8%, respectively) (Additional file 1). Dimensional regression and cross-sample comparison of the replicates showed low differences within strain replicates and higher differences between strains. When plotted against the first and second principal components, replicates clustered within strains (with the fruA mutant having the weakest intrastrain clustering) and distinctly between strains. The main variance was that of the wild-type (WT) strain and the csgA mutant on one hand versus the mrpC and fruA mutants on the other. We conclude that the transcript datasets are in high agreement within strain replicates and show distinct transcriptional profiles between the strains.

Expression of known and predicted regulatory targets in the csgA, fruA, and mrpC mutants is mostly as expected

The formation of M. xanthus fruiting bodies depends on the products of genes regulated by CsgA, FruA, and MrpC, with a loss of their function leading to deficient fruiting body formation under starvation conditions [13, 58, 62, 65]. To determine whether our RNA-seq data agree with previous work, we cataloged developmentally regulated genes whose expression was known to be altered in the mutants (Table 1). For example, expression from the promoter of the dev operon was known to be reduced in all three mutants (designated c- f- m-) relative to the WT strain (see Table 1 for the first gene of the operon, devI, and for references), in agreement with our data (see Additional file 2 for other genes in the dev operon). Expression of the fmgA, fmgD, and fmgE genes, and the fmgBC and fadIJ operons, was likewise reported previously to exhibit the c- f- m- pattern, and our data largely agree (Table 1 and Additional file 2).

Table 1 Comparison of RNA-seq results with regulation reported previously

Our RNA-seq data also agreed in most other cases with the patterns of gene regulation in all three mutants reported previously. For example, exoA and exoL were downregulated in the csgA mutant and upregulated in the fruA and mrpC mutants (c- f + m+) based on RT-qPCR analysis [70]. Our RNA-seq data agreed with the exoA and exoL results (Table 1) and showed similar results for the other genes in the exoA-I and exoL-P operons (Additional file 2). The nfsA gene exhibited a different pattern of regulation (c- f + m-) than exoA and exoL using RT-qPCR [70], and our RNA-seq data agreed (Table 1). Interestingly, our data revealed differential regulation of genes in the nfsA-H operon (Additional file 2). While nfsA-C transcript levels showed the c- f + m- pattern, those of nfsD-E were similar in the mrpC mutant and the WT strain (c- f+), and those of nfsF-H were upregulated in the mrpC mutant (c- f + m+). Suboperonic promoters upstream of nfsD and nfsF [57] may account for differential regulation of these genes in the mrpC mutant (see Discussion). Expression of mrpA and sdeK was reported to be < 2-fold regulated in csgA and fruA mutants, and reduced in an mrpC mutant (m-), and our RNA-seq data agreed (Table 1) and showed similar results for mrpB of the mrpAB operon (Additional file 2). The expression of 11 other genes has been measured in two of the three mutants. Our RNA-seq data agreed with prior work except in two instances (Additional file 2).

We conclude that expression of known regulatory targets was altered in the mutants in our RNA-seq dataset as expected from previous work, with only a few exceptions. The high degree of agreement indicates our data can be mined for new insights with confidence. One new insight is that genes in the nfsA-H operon are differentially regulated in the mrpC mutant.

Differential gene expression in the mrpC mutant correlates with position of MrpC DNA-binding sites upstream of transcriptional start sites

To explore whether our RNA-seq data could serve to augment previous findings, we hypothesized that genes differentially expressed in the mrpC mutant (relative to the WT strain) would correlate with the position of MrpC DNA-binding sites identified by ChIP-seq [77] relative to developmental TSSs identified by Cappable-seq [57]. In general, transcriptional activators typically bind slightly upstream of or partially overlapping with the promoter region bound by RNA polymerase (-50 to + 1) and repressors bind within the promoter region [104,105,106], but we nevertheless included MrpC-binding sites between − 200 to + 200 bp relative to a TSS in our analysis. As multiple potential TSSs were assigned to many of the genes [57], we simplified our analysis by taking the TSS with the highest number of reads for a given gene, and determining the nearest MrpC-binding site peak coordinate [77] relative to that TSS. We observed an abundance of genes with an MrpC-binding site upstream of a TSS, with an above-average concentration spanning from − 115 to + 5 bp and a maximum at -80 bp (Fig. 2a). We then segregated the strongly expressed genes (log10 TSS max count ≥ 1.45, which is ~ 28 reads) that were downregulated in the mrpC mutant from those upregulated in the mutant (Additional file 3), and plotted their differential expression against the position of the MrpC-binding site peak relative to the TSS (Fig. 2b and c). We found 42 genes strongly downregulated in the mrpC mutant (i.e., log2 fold difference < -3) that had an MrpC-binding site peak located between − 120 and − 40 (Fig. 2b and Additional file 3), consistent with strong transcriptional activation by MrpC bound upstream of RNA polymerase. In comparison, a plot of genes upregulated in the mrpC mutant shows a broader distribution of MrpC-binding site locations (Fig. 2c and Additional file 3). Interestingly, nine of the strongly upregulated genes had an MrpC-binding site located between − 105 and − 60, perhaps due to MrpC acting as a repressor by binding to a site that prevents a transcriptional activator from binding or functioning in a subsequent step of activation (see Discussion). These results illustrate the value of combining our RNA-seq data with previous data to discern general features of transcriptional activation and repression by MrpC.

Fig. 2
figure 2

Comparative analysis of M. xanthus MrpC ChIP-seq data collected at 18 h poststarvation [77], 18-h Cappable-seq data [57], and our 24-h RNA-seq data. (A) Combined histogram and density plot for ChIP-seq peak distances from the most abundant transcriptional start site (TSS) of the nearest M. xanthus gene. Distance segments corresponding to density levels > 15% of the maximum density are colored in red in the histogram. Strongly expressed genes (log10 TSS max count ≥ 1.45; Additional file 3) [57] downregulated (B) or upregulated (C) in the mrpC mutant relative to the wild-type strain based on our RNA-seq data, plotted against the ChIP-seq peak [77] distances from TSS shown in panel A. Points with log2 fold differences < -3 or > 3 (horizontal dashed lines in panels B and C, respectively) and within the density peak range (vertical dashed lines) are colored red. The point corresponding to dmxB is indicated in panel C

Large numbers of genes are up- or down-regulated in the mutants

To determine the extent of positive versus negative effects on developmental gene expression in the csgA, fruA, and mrpC mutants, we compared the log2 fold difference between each mutant and the WT strain, sorting by increasing difference in expression (Fig. 3a, left). The numbers of differentially expressed genes were similar for the csgA and mrpC mutants, with 1339 and 1410 genes upregulated, and 1455 and 1450 downregulated, respectively (Additional file 4). They were lower for the fruA mutant, with 1140 genes upregulated and 1010 downregulated, consistent with FruA/FruA* acting downstream of MrpC and CsgA in a regulatory pathway (Fig. 1), but MrpC and CsgA acting in one or more other pathways as well.

Fig. 3
figure 3

Patterns of differential gene expression in mutants relative to wild type. (A) Stacked bar graph of the numbers of genes binned by log2 fold difference. Numbers of genes and color coding are in the table. (B) Venn diagrams of genes downregulated or upregulated in one or more mutants. The number of genes in each segment is given with the percentage relative to the number of genes in all segments. (C) Bar graphs of the numbers of genes with the indicated pattern of regulation in two (left) or all three (right) mutants. Abbreviations are c, f, and m for csgA, fruA, and mrpC mutants, and + or - for up- or down-regulated in mutants relative to wild type

Most of the upregulated genes were increased in expression near the lower end of our cutoff, between 2- to 4-fold (i.e., log2 fold difference 1 to 2) (Fig. 3a, right). This was particularly noteworthy for the csgA mutant, with 89.8% of the upregulated genes being in this grouping. Only 8 genes, or 0.6% of the upregulated genes in the csgA dataset, increased in expression > 8-fold. This is much less than the 25.8% and 28.9% of upregulated genes with > 8-fold increase for the fruA and mrpC mutants, respectively. Furthermore, 20.3% of the downregulated genes in the csgA dataset decreased in expression more than 8-fold, highlighting a much greater role of CsgA as a strong positive regulator of gene expression than as a strong negative regulator. Conversely, FruA has a greater role as a strong negative regulator since nearly twice as many genes increased expression > 8-fold in the fruA mutant compared with the number that decreased expression more than 8-fold. It is important to note that, as expected [56], fruA was downregulated strongly in the mrpC mutant (64-fold, Additional file 2), and FruA is undetectable by immunoblot in the mrpC mutant at 18–30 h PS [65], which presumably accounts for some of the effects on gene expression observed in the mrpC mutant. Yet, the numbers of genes strongly upregulated (408 genes) or downregulated (322 genes) in the mrpC mutant (using > 8-fold difference as a cutoff) are larger than the numbers in the fruA mutant (294 upregulated, 166 downregulated), indicating that MrpC is a strong negative and positive regulator of many genes independently of its effect on the FruA level.

We conclude that CsgA, FruA, and MrpC are negative and positive regulators of > 1000 genes during development (Fig. 3a). CsgA-dependent regulation is indirect and its strongest effects are positive, with the bulk of its negative effects being comparatively weaker. Regulation by transcription factors FruA and MrpC may be direct or indirect. FruA strongly affects hundreds of genes, especially as a negative regulator. MrpC not only affects genes by strongly stimulating fruA transcription, but it also strongly regulates hundreds of other genes negatively or positively apart from its effect on the FruA level.

CsgA, FruA, and MrpC are coregulators and counterregulators of many developmental genes

To determine the numbers of coregulated, counterregulated, and uniquely regulated genes among the mutants, we compiled lists of the differentially expressed genes observed in only one mutant, two mutants, or all three mutants, categorized according to their regulatory pattern (Additional files 4 and 5). We found that 382 genes were downregulated and 157 were upregulated in all three mutants (Fig. 3b), indicating a larger role of CsgA, FruA, and MrpC as positive coregulators of gene expression than as negative coregulators. Positive coregulation was known for a few genes from previous work (Table 1) and our RNA-seq data showed that downstream genes in the same operons were regulated similarly (Additional file 2), as described above (i.e., genes exhibiting the c- f- m- pattern of regulation in the mutants). In contrast, negative coregulation by CsgA, FruA, and MrpC (i.e., a c+ f + m+ pattern) is a novel finding of this study.

Among the genes coregulated in two of the three mutants, the greatest number were in the fruA and mrpC mutants, with 421 upregulated and 94 downregulated (Fig. 3b). The strong dependence of FruA on MrpC [56, 65] (Additional file 2) and cooperative binding of FruA and MrpC to DNA [63, 66,67,68,69, 77] likely contribute to this high frequency of coregulation. Surprisingly large numbers of genes were coregulated by CsgA and either FruA or MrpC, underscoring regulation not anticipated from the pathway depicted in Fig. 1. For the csgA and fruA mutants, 275 genes were downregulated and 114 were upregulated (Fig. 3b). Since FruA depends strongly on MrpC [56, 65] (Additional file 2), identifying genes impacted in the csgA and fruA mutants, but not the mrpC mutant, was particularly surprising (see Discussion). For the csgA and mrpC mutants, 216 genes were downregulated and 161 were upregulated (Fig. 3b). Perhaps regulation of these genes reflects partial dependence of CsgA on MrpC (Fig. 1). The CsgA protein level [58] and the csgA (MXAN_1294) transcript level (Additional file 5) are reduced in the mrpC mutant. How an elevated level of CsgA in the WT strain would positively regulate gene expression apart from FruA (i.e., < 2-fold regulation in the fruA mutant) is unknown.

While the general observation was that CsgA, FruA, and MrpC tend to co-positively or co-negatively regulate target genes as described above, another surprising observation was that sizable numbers of genes were counterregulated across two (Fig. 3c, left) or three (Fig. 3c, right) of the mutants. Genes counterregulated in csgA, fruA, and/or mrpC mutants have been described in several previous studies (Table 1 and Additional file 2), but the extent of counterregulation genome-wide was unknown. Strikingly, 249 (40%) of the genes regulated by CsgA and MrpC (henceforth called the CsgA/MrpC dataset) were counterregulated (Fig. 3c, left), including two previously reported c-di-GMP-associated genes, MXAN_4232 and MXAN_2902 (Additional file 2). In contrast, only 49 (11%) and 101 (16%) of the genes in the CsgA/FruA and FruA/MrpC datasets, respectively, were counterregulated (Fig. 3c, left). Hence, CsgA and MrpC exert opposing regulatory effects on many more genes than does FruA in combination with CsgA or MrpC. Interestingly, FruA regulation rarely opposed that of both CsgA and MrpC in the CsgA/FruA/MrpC dataset, where 218 (28.8%) of the genes were counterregulated (Fig. 3c, right), but only 33 (4.4%) were upregulated in the fruA mutant and downregulated in the csgA and mrpC mutants (i.e., c- f + m-, like nfsA-C and pkn1 in Table 1 and Additional file 2), and no genes exhibited the c+ f- m+ pattern (Fig. 3c, right). The tendency of FruA to coregulate rather than counterregulate with CsgA and MrpC is consistent with FruA/FruA* acting downstream of MrpC and CsgA in the regulatory pathway depicted in Fig. 1.

We were also surprised to find large numbers of uniquely-regulated genes in each mutant, including both up- and down-regulated genes (Fig. 3b). The numbers of uniquely-regulated genes in the csgA and mrpC mutants (973 and 861, respectively) were much larger than in the fruA mutant (339), consistent with FruA/FruA* acting downstream of CsgA and MrpC to regulate many genes as depicted in Fig. 1. However, finding any uniquely-regulated genes in the fruA mutant is surprising given the strong dependence of FruA on MrpC [56, 65] (Additional file 2) (see Discussion).

We conclude that CsgA, FruA, and MrpC are primarily positive coregulators, consistent with the pathway depicted in Fig. 1. The most common pair of coregulators is FruA and MrpC, presumably due to dependence of FruA on MrpC and/or the two transcription factors binding cooperatively to DNA. Typically, the pair act negatively when CsgA does not appear to be involved (i.e., < 2-fold expression difference in the csgA mutant). Conversely, pairwise coregulation by CsgA and FruA was primarily positive, perhaps related to CsgA-mediated C-signaling leading to FruA* formation and transcriptional activation of late developmental genes.

Genes most impacted are coregulated in two or three of the mutants

The previous section focused on the numbers of genes exhibiting particular patterns of regulation in the mutants, without regard to the strength of regulation. In this section, we focus on the strength of regulation or more precisely the impact (strength and statistical power) of the regulation. To determine the genes most impacted in the csgA, fruA, and/or mrpC mutants based on their strength of differential expression (log2 fold difference) and statistical power (-log10 FDR-adjusted p-value), we prepared a volcano plot for each mutant and colored the dots (representing the genes) according to whether the gene was differentially expressed in one (green), two (red), or all three (blue) mutants (Fig. 4a). We observed two general trends. First, the bulk of the most impacted downregulated genes were differentially expressed in all three mutants, whereas the upregulated ones were differentially expressed in two or three of the mutants. Second, for both up- and down-regulated genes, as the strength of differential expression increased and the -log10 FDR-adjusted p-value increased, the genes tended to be differentially expressed in two or three of the mutants. Both trends indicate that the genes most impacted are regulated in two or three of the mutants.

Fig. 4
figure 4

Scatter plots of differential gene expression in mutants relative to wild type. (A) Volcano plots for the mutants. The top ten strongest differentially expressed and statistically significant genes are illustrated with larger points. Colors indicate differential expression in one (green), two (red), or all three (blue) mutants. (B) Differential expression of genes in the following pairs of mutants: fruA and mrpC (red), csgA and fruA (green), or csgA and mrpC (blue). (C) Differential expression of genes in all three mutants. Axes show strength of differential expression in fruA and mrpC mutants, and color shows strength of differential expression in the csgA mutant

We also prepared lists of the ten most impacted up- or down-regulated genes for each mutant (Additional file 6). Among the downregulated lists, the genes were distinct for the csgA mutant, with the exception of katE, which was also in the mrpC mutant list. In contrast, six of ten genes were the same for the fruA and mrpC mutants. One of six genes (MXAN_7102) was not downregulated in the csgA mutant (i.e., exhibited the f- m- pattern), but all the other most impacted downregulated genes exhibited the c- f- m- pattern, highlighting the strong positive coregulation by CsgA, FruA and MrpC (Additional file 6, bold). Among the upregulated lists, only MXAN_2809 was on more than one list. This gene exhibited the most prevalent upregulated pattern (f + m+) for the fruA (9/10) and mrpC (5/10) mutant lists, followed by c+ f + m+ (1/10 and 4/10 for fruA and mrpC, respectively). We conclude that FruA and MrpC are strong negative coregulators, often independently of CsgA, which we examine further below. For the csgA mutant, the upregulated patterns varied. Finally, we note that many genes (25/60) on the lists are likely in operons with other genes on the lists.

We examined the differential expression strength for genes affected in two of the three mutants and found that it varied depending on the pairing. For example, of the 421 genes upregulated in both the fruA and mrpC mutants, the strength of upregulation varied over a broad range (Fig. 4b, top panel, upper right quadrant) and was highly correlated via Spearman analysis (rs = 0.82, Additional file 7, top right graph). In agreement with our conclusion above that FruA and MrpC are strong negative coregulators, often independently of CsgA, we note that 91 genes were strongly upregulated (log2 fold difference > 5) in both the fruA and mrpC mutants (Fig. 4b, top panel, upper right quadrant). Comparatively, the number of genes downregulated in the fruA and mrpC mutants was much smaller at 94 (Fig. 4b, top panel, lower left quadrant). The strength of downregulation varied over a broad range (as for upregulation), but only 11 genes were strongly downregulated (log2 fold difference < -5) in both mutants. The strength of downregulation was highly correlated in the two mutants (rs = 0.67, Additional file 7, top left graph), albeit slightly lower than for upregulation. The highly correlated strength of regulation may reflect the strong dependence of FruA on MrpC [56, 65] (Additional file 2) and/or cooperative binding of FruA and MrpC to DNA at regulatory sites [63, 66,67,68,69, 77]. Notably, we observed no strongly counterregulated genes (log2 fold difference > 5 in one mutant and < -5 in the other) (Fig. 4b, top panel, upper left and lower right quadrants).

For genes affected in the csgA and fruA mutants, downregulation tended to be stronger than upregulation (Fig. 4b, middle panel, lower left and upper right quadrants, respectively). Both included many genes more strongly affected in the csgA mutant than in the fruA mutant (i.e., many dots closer to the x-axis than to the diagonal). The correlations were moderate to nonexistent (rs = 0.49 to 0.11 for down- and up-regulated genes, respectively).

For genes affected in the csgA and mrpC mutants, there was also skewing toward stronger downregulation than upregulation (Fig. 4b, bottom panel, lower left and upper right quadrants, respectively). Both showed many genes more strongly affected in one of the mutants (i.e., many dots off the diagonal), so the correlations were low to nonexistent (rs = 0.33 to -0.07 for down- and up-regulated genes, respectively).

We also examined the differential expression strength for genes affected in all three mutants and observed some similar trends as mentioned above. Like the 421 genes upregulated in the fruA and mrpC mutants, for the 157 genes upregulated in all three mutants, the strength of upregulation in the fruA and mrpC mutants varied over a broad range (Fig. 4c, upper right quadrant) and was highly correlated (rs = 0.76, Additional file 7, bottom right graph). As the strength of upregulation increased in the fruA and mrpC mutants, that in the csgA mutant tended to decrease (Fig. 4c, upper right quadrant). For the 382 genes downregulated in all three mutants, the strength of downregulation in the fruA and mrpC mutants also varied over a broad range (Fig. 4c, lower left quadrant). However, many genes were more strongly affected in one of the mutants (i.e., many dots off the diagonal), so the correlation was low (rs = 0.38, Additional file 7, bottom left graph). In this case, as the strength of downregulation increased in the fruA and mrpC mutants, that in the csgA mutant also tended to increase (Fig. 4c, lower left quadrant).

Altogether, our analysis of gene regulatory strength indicates that FruA and MrpC are strong negative coregulators of many genes independently of CsgA. In conjunction with CsgA, FruA and MrpC negatively coregulate less genes, but still in a highly correlated fashion. More genes are positively coregulated by CsgA, FruA, and MrpC, but interestingly the correlation between strength of regulation by FruA and MrpC is low (see Discussion).

CsgA, FruA, and MrpC are regulators of genes involved in cellular signaling, protein synthesis, energetics, and envelope function

To determine if specific functional families of genes are preferentially regulated by CsgA, FruA, and MrpC, we performed Gene Ontology (GO) term overrepresentation analysis on the various lists of up- and down-regulated genes (Additional file 4) using two separate platforms. ShinyGO queries multiple annotation databases, such as Ensembl, to calculate enriched GO terms for a given list of genes [92]. An FDR-adjusted p-value is determined for each fold-enrichment score, with a cutoff of less than 0.05 used for our analysis. We queried the gene lists related to Fig. 3a (i.e., all genes up- or down-regulated in each mutant, designated as csgA all, fruA all, and mrpC all in Additional file 4) and the shorter gene lists related to the groupings depicted in Fig. 3b and c (designated csgA only, fruA only, mrpC only, csgA/fruA, csgA/mrpC, fruA/mrpC, and csgA/fruA/mrpC, also in Additional file 4). As an output comparison, we repeated this process using DAVID [93, 94]. Like ShinyGO, DAVID utilizes several annotation databases, including Entrez Gene and Uniprot, to ascertain statistically enriched GO terms. The breadth of databases utilized by DAVID is wider, but the outputs of both tools should be largely in agreement given the similar nature of the approaches. The output GO term lists were compared to one another, and GO terms that were statistically enriched and in agreement between the two analytical methods are reported in Fig. 5 and Additional file 8.

Fig. 5
figure 5

Statistically enriched Gene Ontology (GO) groupings identified by both ShinyGO and DAVID analyses. (A) GO groupings enriched in each individual mutant regardless of gene presence in other datasets (designated ALL in Additional file 8). (B) GO groupings enriched considering only genes uniquely regulated in a single mutant (designated ONLY in Additional file 8). (C) GO groupings enriched considering genes with particular patterns of co- or counter-regulation. Gene regulatory patterns are abbreviated as explained in the Fig. 3c legend. In all three panels, the RichFactor is defined as the number of differentially expressed genes over the total number of genes in the GO term, the numbers of differentially expressed genes are to the right of the dots, and a Benjamini score of < 0.05 was used as a statistical cutoff

When looking at all the genes up- or down-regulated in each individual mutant (i.e., the csgA all, fruA all, and mrpC all lists in Additional file 4), several observations stood out. First, each of the mutants had downregulated genes enriched for ribosomal proteins and rRNA (Fig. 5a, left), with the csgA mutant having the most at 45 genes, nearly the entire GO term list (RichFactor nearly 1). Hence, CsgA, FruA, and MrpC are positive regulators of many genes necessary for protein synthesis during development. Second, each mutant also had downregulated genes enriched for one other GO term. The mrpC mutant was enriched for “defense response to virus” genes, including two type I CRISPR-Cas systems: the dev operon spanning MXAN_7266−7259 and MXAN_7014–7020. The latter was weakly impacted only by mrpC, with average log2 fold difference near − 1 (Additional file 5), while the former was affected in each mutant, though more strongly so in the mrpC and fruA mutants, with average log2 fold differences of about − 5 (Additional file 2). It should be noted that the M. xanthus GO term list for “defense response to virus” genes is limited and does not include every gene within these CRISPR-Cas clusters. This accounts for the RichFactor = 1 with only 11 genes (Fig. 5a, left) (absent from the GO term list are MXAN_7266, 7014, 7017, and 7018 in the two type I systems, and the entire type III system MXAN_7276–7283). The fruA and csgA mutants were enriched for genes related to energy production and conversion, with the former controlling the electron transport chain via NADH dehydrogenase subunit expression and the latter affecting energy storage via fatty acid biosynthesis. Third, only genes upregulated in the csgA mutant yielded GO terms that met our criteria, and all three GO terms were related to signaling pathways: phosphorelay signal transduction systems, sensor kinases, and protein Ser/Thr kinases (Fig. 5a, right). Between ~ 40–50% of all genes in these pathways were present within the dataset (i.e., RichFactor ~ 0.4–0.5), highlighting CsgA’s role as a negative regulator of signaling pathways at 24 h PS.

Our analyses of genes uniquely regulated in a single mutant (i.e., the csgA only, fruA only, and mrpC only lists in Additional file 4) yielded GO terms distinct from those described above for the “all” lists. Interestingly, genes uniquely downregulated in the csgA mutant were enriched in four GO terms related to cellular shape change (cell division, peptidoglycan biosynthesis, cell wall organization, and regulation of cell shape), while genes uniquely downregulated in the mrpC mutant included most genes involved in histidine biosynthesis (Fig. 5b, left). These results suggest that CsgA and MrpC play unique roles in positive regulation of genes involved in cellular shape change and histidine biosynthesis, respectively, during development. Genes uniquely upregulated in the csgA mutant were enriched in the broad cellular compartment category of “integral component of membrane” (Fig. 5b, right).

Our analyses of co- and counterregulated gene lists yielded mostly the same GO terms as those described above for the “all” lists, albeit with slightly to considerably lower gene numbers and RichFactors. For example, analysis of coregulated genes with the c- f- m- pattern yielded the “ribosome” GO term with 30 genes and RichFactor = 0.63 (Fig. 5c and Additional file 8), slightly lower than the 31 downregulated genes and RichFactor = 0.65 observed for the mrpC mutant and considerably lower than for the csgA and fruA mutants (Fig. 5a, left), whose “all” lists have an additional 12 “ribosome” genes with the coregulated c- f- pattern (Fig. 5c). Likewise, the 7 “defense response to virus” genes of the dev operon (devI, aka MXAN_7266, is absent from the GO term list, as noted above) are coregulated (c- f- m-) (Fig. 5c and Additional files 2 and 8), but the mrpC all list has 4 additional genes with this GO term (Fig. 5a, left), and the fruA all list has 1 gene with the NADH dehydrogenase GO term in addition to the 8 with the coregulated c- f- pattern (Fig. 5a and c). Only counterregulated genes with the c+ m- pattern yielded GO terms that met our criteria. Among the upregulated genes on the csgA “all” list that yielded GO terms related to signaling pathways (Fig. 5a, right), considerably smaller numbers were also downregulated in the mrpC mutant and yielded two of the same GO terms (Fig. 5c). These 12 signal transduction system genes and 9 sensor kinase genes are negatively regulated by CsgA, but positively regulated by MrpC.

Overall, our GO term overrepresentation analysis supports broad alteration of cellular signaling, protein synthesis, energetics, and envelope function by CsgA, FruA, and/or MrpC at the developmental timepoint prior to sporulation.

CsgA, FruA, and MrpC counterregulate genes corresponding to particular predicted STRING protein-protein interaction networks

As an alternative approach to glean functional information about genes counterregulated in the csgA, fruA, and/or mrpC mutants, we utilized the Cytoscape software platform [91, 107] to map out and visualize putative STRING interaction networks. Several of the predicted STRING interaction networks are shown in Fig. 6. These and other networks are viewable in Cytoscape (Additional file 9). Several networks corresponding to counterregulated genes grouped broadly in sugar metabolism: polysaccharide biosynthetic process/glycosyl transferase proteins to the c + f- m- or c- f + m + regulatory pattern (Fig. 6a), maltose/maltodextrin import proteins to f + m- (Fig. 6b), and galactose metabolism proteins to c- f+ (Fig. 6d). Additionally, we identified a network involved in the non-ribosomal synthesis of siderophores corresponding to the f + m- regulatory pattern (Fig. 6b) and networks of Ser/Thr kinases, histidine kinases, and response regulators corresponding primarily to c + m- (Fig. 6c). Much additional work will be required to determine the significance of the predicted STRING interaction networks whose corresponding genes are counterregulated by CsgA, FruA, and MrpC, as well as the many other networks viewable in Cytoscape (Additional file 9).

Fig. 6
figure 6

Select STRING interaction networks for counterregulated genes. (A) Genes differentially expressed in all three mutants. Regulatory patterns: c + f- m- (gray), c- f + m + (yellow), or c- f- m + (light blue). (B-D) Genes differentially expressed in two mutants. Regulatory patterns: (B) f + m-, (C) c + m- (green) or c- m + (magenta), or (D) c- f +. Numbers indicate MXAN_#### gene designations (aka old locus ID). See the Fig. 3C legend for regulatory pattern abbreviations

Gene set enrichment analysis identifies several enriched pathways controlled in combination by CsgA, FruA, and MrpC

To determine if certain cellular pathways were affected in the csgA, fruA, and/or mrpC mutants, we performed KEGG pathway gene set enrichment analysis using the gseKEGG function in the clusterProfiler R-package on our segmented or combined gene lists. The segmented lists of genes were those differentially expressed relative to the WT strain in only one mutant or coregulated in two mutants or all three. The combined list included all genes from the segmented lists. We identified 18 unique pathways across five gene lists with a q-value ≤ 0.05 (Additional file 10). The KEGG pathways mirror the findings of our GO term analysis, with a strong emphasis on protein synthesis (ribosome and aminoacyl-tRNA biosynthesis), energetics (fatty acid/branched-chain amino acid metabolism, TCA cycle, and oxidative phosphorylation), and cellular shape (peptidoglycan biosynthesis). We investigated three KEGG pathways further: the ribosomal pathway (mxa03010) given its strong downregulation across all mutants, BCAA metabolism (mxa00280) due to its upstream importance for production of key developmental lipids, and PG biosynthesis (mxa00550) owing to the cell wall remodeling necessary for cellular shape change during development.

We identified several points upstream of rRNA and ribosomal protein production that are impacted by loss of CsgA, FruA, and/or MrpC function. Figure 7a summarizes previous knowledge of the M. xanthus ribosomal pathway and indicates whether expression of genes corresponding to the proteins shown is up- or down-regulated in mutants, based on our RNA-seq data (Fig. 7b and Additional file 5). We observed downregulation of socE (MXAN_0731) and upregulation of relA (MXAN_5201), nla18 (MXAN_3692), and nla4 (MXAN_2516) in the indicated mutants (Fig. 7a). All these effects potentially increase (p)ppGpp levels in the mutants, which could explain the observed downregulation of most rRNAs and genes for ribosomal proteins (Fig. 7). Since downregulation of the ribosomal pathway requires not only elevated (p)ppGpp levels, but also elevated levels of uncharged tRNA, we examined expression of genes involved in tRNA charging. Of the 25 M. xanthus aminoacyl-tRNA synthetase genes, 17, 10, and 5 were downregulated in the csgA, fruA, and mrpC mutants, respectively (Fig. 7b). Thus, we propose that in the WT strain, CsgA, FruA, and MrpC positively regulate aminoacyl-tRNA synthetase genes and negatively regulate relA, nla18, and nla4, respectively, while MrpC positively regulates socE, and all these effects decrease (p)ppGpp levels, diminishing the stringent response to starvation and preserving ribosomes for translation of developmental mRNAs into proteins.

Fig. 7
figure 7

Ribosomal pathway genes differentially expressed in mutants relative to wild type. (A) Control of ribosome biogenesis by the stringent response to starvation and summary of mutant effects. Starvation inhibits aminoacyl-tRNA synthetases due to amino acid limitation, leading to uncharged tRNA in ribosomes and causing RelA to synthesize (p)ppGpp. SocE inhibits the RelA-dependent stringent response [108], whereas Nla18 [46] and Nla4 [47] stimulate (p)ppGpp accumulation. (p)ppGpp and one or more DksA orthologs [109] may inhibit RNA polymerase transcription of genes for rRNAs, diminishing ribosome biogenesis. Gene regulatory patterns are indicated in parenthesis using abbreviations explained in the Fig. 3c legend. (B) Heatmaps of differential expression of genes for stringent response proteins, aminoacyl-tRNA synthetases, and ribosomal proteins and RNA in csgA, fruA, and mrpC mutants relative to wild type based on our RNA-seq data (Additional file 5)

We found that genes coding for many enzymes involved in the metabolism of BCAAs are downregulated in the csgA, fruA, and/or mrpC mutants. Figure 8a shows the affected part of the BCAA metabolic pathway and indicates whether expression of genes corresponding to the enzymes shown is up- or down-regulated in mutants, based on our RNA-seq data (Fig. 8b and Additional file 5). In relation to this pathway (Fig. 8a), a critical early discovery was that loss of the branched-chain keto acid dehydrogenase (BCKAD) subunits (bkd/esg, MXAN_4264 and MXAN_4265) blocks development, and the defect can be rescued by providing isovaleric acid, which the mutant presumably converts to isovaleryl-CoA, a primer for biosynthesis of long branched-chain fatty acids such as the abundant iso15:0 species [17]. In the WT strain, BCKAD converts the leucine derivative 2-ketoisocaproate (not shown) to isovaleryl-CoA, which fatty acid synthase elongates to iso15:0 (Fig. 8a, vertical pathway at left). The iso15:0 fatty acid is a precursor of TG-1 triglyceride, and either can rescue development of a bkd/esg mutant, suggesting they are developmental signals [18,19,20]. A second pathway for the formation of isovaleryl-CoA was identified that relies on MvaS (MXAN_4267) conversion of acetyl-CoA and acetoacetyl-CoA to (S)-3-hydroxy-3-methylglutaryl-CoA (Fig. 8a, dashed arrows), followed by subsequent activity of LiuC (MXAN_3757), BCKAD, and MXAN_4266 [110, 111]. Notably, isoleucine can contribute to isovaleryl-CoA production via the MvaS-dependent pathway since isoleucine can be converted to (S)-2-methyl-butanoyl-CoA and lead to acetyl-CoA (Fig. 8a, center vertical pathways can form a loop). Strikingly, the genes corresponding to nearly all the enzymes involved in converting branched-chain keto acids (BCKAs) to iso15:0 are downregulated in the csgA, fruA, and/or mrpC mutants, with a few exceptions (see Discussion), most notably the initial conversion of BCAAs to BCKAs being upregulated in the csgA mutant (Fig. 8). Overall, our results imply that in the WT strain, CsgA, FruA, and MrpC positively regulate genes whose products generate isovaleryl-CoA, iso15:0, and TG-1, consistent with observations related to the abundance of lipid bodies and lipid signals in the WT strain and csgA, fruA, and mrpC mutants during development (see Discussion). Separate from synthesis of isovaleryl-CoA, iso15:0, and TG-1, genes for energy metabolism are positively regulated by CsgA, FruA, and MrpC, presumably increasing isobutyryl-CoA from valine catabolism, eventually being converted to succinyl-CoA via multiple enzymatic steps (Fig. 8a, vertical pathway at right).

Fig. 8
figure 8

Branched-chain amino acid (BCAA) metabolic genes differentially expressed in mutants relative to wild type. (A) BCAA degradation leading to biosynthesis of lipid signals involved in fruiting body formation and summary of mutant effects. BCAAs are converted to branched-chain keto acids (BCKAs), which in turn are converted to the indicated acyl-CoA derivatives. Isovaleryl-CoA can feed into fatty acid biosynthesis, including that of the abundant iso15:0 species, which is a precursor of TG-1 triglyceride (vertical pathway at left). The iso15:0 fatty acid and TG-1 are lipid signals involved in fruiting body formation. (S)-2-methyl-butanoyl-CoA (derived from Ile) can lead to acetyl-CoA and acetoacetyl-CoA, which via MvaS activity yields (S)-3-hydroxy-3-methylglutaryl-CoA that can lead to isovaleryl-CoA (i.e., the two center vertical pathways form a loop). Isobutyryl-CoA (derived from Val) can lead to succinyl-CoA (vertical pathway at right) that feeds into metabolism and energy production and conversion pathways. Asterisks indicate enzymatic steps for which corresponding M. xanthus gene homologs were not apparent in the KEGG database. Gene regulatory patterns are indicated in parenthesis using abbreviations explained in the Fig. 3c legend. (B) Heatmaps of differential expression of genes for BCAA degradation and fatty acid biosynthesis in csgA, fruA, and mrpC mutants relative to wild type. Since many enzymes participate in the synthesis of fatty acids, the corresponding M. xanthus genes are not shown separately in the vertical pathway at the left in panel A. Rather, they are summarized as “Fatty acid biosynthesis” and their typical regulatory pattern (c- f-) is indicated. The individual genes are listed in functional groups at the right in panel B and the heatmap shows their differential expression based on our RNA-seq data (Additional file 5)

Finally, we draw attention to the PG biosynthetic pathway, in which MrpC negatively regulates genes involved in the initial and penultimate steps, but CsgA positively regulates genes involved in intermediate steps. Figure 9a shows the affected portion of the PG pathway and indicates whether expression of genes corresponding to the enzymes shown is up- or down-regulated in mutants, based on our RNA-seq data (Fig. 9b and Additional file 5). Genes coding for enzymes involved in the production of the four principal inputs to PG biosynthesis, UDP-GlcNAc-enolpyruvate (MXAN_4909, MurA), (D-Ala)2 (MXAN_5741 and MXAN_5951, Ddl enzymes), m-diaminopimelic acid (m-DAP) (MXAN_5054, DapF; MXAN_0971 and MXAN_4877), and undecaprenyl phosphate (Und-P) (MXAN_0609), were upregulated in the mrpC mutant, as were genes for class A penicillin-binding proteins (MXAN_5181 and MXAN_5911, PbpAs), which are glycosyltransferase-transpeptidases that synthesize glycan chains using the lipid II precursor and then cross-link the peptide side chains at penultimate steps of the pathway. Hence, MrpC negatively regulates genes important for PG synthesis, an appropriate response during starvation-induced developmental growth arrest. Conversely, genes corresponding to enzymes involved in intermediate steps of the pathway were downregulated in the csgA mutant, including genes for MurB-G and MraY (MXAN_5602–5604 and MXAN_5606–5609) that convert UDP-GlcNAc-enolpyruvate to the lipid II precursor. Previous investigations provide insight into the negative regulation by MrpC that we observe, and suggest explanations for the novel positive regulation of PG biosynthesis by CsgA (see Discussion).

Fig. 9
figure 9

Peptidoglycan (PG) biosynthesis genes differentially expressed in mutants relative to wild type. (A) Overview of PG biosynthesis and summary of mutant effects. Amino acids, sugar nucleotides, and phospholipids are converted into UDP-GlcNAc-enolpyruvate, (D-Ala)2, m-DAP, and Und-P, which are the four principle inputs for synthesis of the lipid II precursor, which is used by glycosyltransferase-transpeptidases to synthesize the glycan chains and cross-link the peptide side chains of PG. Gene regulatory patterns are indicated in parenthesis using abbreviations explained in the Fig. 3c legend. (B) Heatmaps of differential expression of genes for enzymes involved in the initial (synthesis of inputs), intermediate (lipid II synthesis), and completion (glycan chain synthesis and peptide side chain cross-linking) stages of PG biosynthesis in csgA, fruA, and mrpC mutants relative to wild type based on our RNA-seq data (Additional file 5)

Discussion

Our results show that loss-of-function mutations in csgA, fruA, or mrpC broadly impact the M. xanthus transcriptome at 24 h PS, a developmental timepoint when the WT strain would have formed mounds, but prior to the transition from rods into spores. We discovered that thousands of genes are up- or down-regulated at least 2-fold in each mutant relative to the WT strain. In agreement with a recent developmental timecourse comparison of fruA mutant and WT transcriptomes [78], our results indicate that FruA is a strong negative regulator of gene expression. Conversely, we found that CsgA is primarily a strong positive regulator. We also found that MrpC is both a strong negative and positive regulator independently of its effect on the FruA level. While only a few different regulatory patterns across the three mutants had been reported previously (Table 1 and Additional file 2), we observed nearly every possible pattern of coregulation, unique regulation, and counterregulation (Fig. 3), revealing much more complicated roles of CsgA, FruA, and MrpC in the developmental gene regulatory network than has been appreciated. In particular, our results show gene set enrichment for pathways involved in cellular signaling, protein synthesis, energetics, and envelope function in the csgA, fruA, and/or mrpC mutants. We drilled deeper into the effects on three pathways – control of ribosome biogenesis by the stringent response to starvation, BCAA metabolism important for developmental lipid signals, and PG biosynthesis since the cell wall undergoes remarkable remodeling during the rod to spore conversion. We discuss our observations related to these pathways, as well as our other results, below.

New insights from comparing genome-wide datasets

Our results add to a growing number of genome-wide studies that provide opportunities for new insights to emerge. For example, our RNA-seq data uncovered three distinct regulatory patterns of genes in the nfsA-H operon (Additional file 2). The patterns differ with respect to the effect in the mrpC mutant (i.e., m- for nfsA-C, < 2-fold regulation for nfsD-E, and m + for nfsF-G). Only the nfsA transcript level had been measured previously in the mrpC mutant [70]. We note that Cappable-seq identified developmental TSSs upstream of nfsA, nfsD, and nfsF [57]. ChIP-seq identified an MrpC-binding site centered at -29 relative to the TSS upstream of nfsA [77]. Downregulation of nfsA-C in the mrpC mutant implies positive regulation by MrpC, which was proposed to involve transcriptional activation by MrpC binding to the site centered at -29, cooperatively with phosphorylated Nla6 binding to a site centered at -42.5, based on in vitro DNA-binding results [70]. We propose that suboperonic promoters upstream of nfsD and nfsF account for < 2-fold regulation of nfsD-E and upregulation of nfsF-G in the mrpC mutant, respectively, due to other mechanisms controlling transcriptional initiation. Premature termination and/or differential stability of transcripts from the three promoters may also contribute to the distinct effects of the mrpC mutation on genes in the nfs operon.

Comparison of the genome-wide Cappable-seq [57] and MrpC ChIP-seq [77] datasets with our RNA-seq differential gene expression for the mrpC mutant relative to the WT strain showed MrpC binding between − 120 and − 40 relative to TSSs correlates with strong transcriptional activation (Fig. 2b and Additional file 3). This position of binding is common for CRP/Fnr superfamily members, such as MrpC, and the mechanism of transcriptional activation usually involves direct interaction with subunit α or another subunit of RNA polymerase to facilitate its recruitment to the promoter [106]. Alternatively or in addition, activation can involve other steps in the initiation of transcription [104, 105]. Establishing whether MrpC follows these paradigms or utilizes novel regulatory mechanisms, and determining what factors influence MrpC’s regulatory strength at particular promoters, warrant further investigation. Interestingly, our dataset comparison also suggested that MrpC binds upstream of promoters to strongly repress transcription (Fig. 2c and Additional file 3). Analysis of mrpC regulation provides an example of how this might work. MrpC negatively autoregulates by outcompeting a transcriptional activator, MrpB (presumably phosphorylated), for binding to overlapping sites upstream of the mrpC promoter [55]. Likewise, MrpC binds upstream of the dmxB promoter and strongly represses transcription (Fig. 2c and Additional file 3), perhaps by outcompeting the binding of a transcriptional activator that remains to be identified [57].

Comparison of our RNA-seq data with other transcriptomic studies related to M. xanthus development also yields some new insights. The most similar study was the recent developmental timecourse comparison of fruA mutant and WT transcriptomes [78]. RNA-seq analysis showed that > 4000 genes were up- or down-regulated in a WT strain during the first 24 h PS, and expression of ~ 1350 genes differed in a fruA mutant relative to a WT strain at 24 h. Using similar (though not identical) methods of submerged culture development, rRNA depletion, and statistical analysis, we found that 2150 genes differed in expression between our fruA mutant and WT strain at 24 h. Since there may be a temporal difference in development between the studies, and McLoon et al. [78] found that ~ 2950 genes differed in expression at 12 h, we performed pairwise comparisons of our 24-h data with their 12- and 24-h data. While all three datasets had high broad agreement in terms of regulation by FruA for both clusters II and IV, which were down- and up-regulated in their fruA mutant, respectively (Additional file 11), the Spearman correlation for strength of differential expression was better for the two 24-h datasets (Additional file 12). Of the genes predicted to be FruA-dependent by McLoon et al. [78] and observed in our study, only 43 are solely controlled by FruA at 24 h, while 276 are controlled by FruA and either CsgA or MrpC, and 268 require all three for normal expression (Additional file 11). The differences and similarities between the work of McLoon et al. [78] and two other developmental timecourse transcriptomic studies [43, 44] have been discussed previously [78]. Common themes included regulation of signaling, translation, and metabolism. Our results show that much of this regulation involves combinations of CsgA, FruA, and/or MrpC.

Diverse regulatory roles of CsgA, FruA, and MrpC

Our results reveal much more complicated roles of CsgA, FruA, and MrpC in the developmental gene regulatory network than was evident from prior work. Early studies of fmg and dev genes suggested that CsgA, FruA, and MrpC are positive coregulators, and our RNA-seq data was in nearly complete agreement (Table 1 and Additional file 2). The early studies support the model that cooperative binding of FruA and MrpC to DNA immediately upstream of the fmg and dev promoters stimulates transcription [63, 66,67,68,69]. CsgA was proposed to activate FruA by an unknown mechanism, forming FruA*, to enable or enhance transcription [63,64,65] (Fig. 1). Our RNA-seq data identified 382 genes that are downregulated in the csgA, fruA, and mrpC mutants (Fig. 3b, left). Some of these genes showing the c- f- m- regulatory pattern likely follow the paradigm of fmg and dev gene regulation involving cooperative binding of FruA/FruA* and MrpC, given that cooperative binding appears to occur commonly [77]. However, the strong dependence of FruA on MrpC [56, 65] (Additional file 2) also predicts coregulation in the fruA and mrpC mutants. Importantly, the strength of regulation in the two mutants is expected to be similar, whether caused by cooperative binding or the dependence of FruA on MrpC, unless in the latter case a small amount FruA/FruA* in the mrpC mutant exerts an effect on gene expression.

Our results suggest that a small amount FruA/FruA* in the mrpC mutant may indeed exert an effect on gene expression. Although the fmg and dev genes exhibited similar strength of downregulation in the fruA and mrpC mutants (Table 1 and Additional file 2), this was not the case for many of the other genes with the c- f- m- regulatory pattern (Additional file 7, bottom left graph). Many of these genes (i.e., those above the diagonal) were more strongly downregulated in the fruA mutant, perhaps due to the complete absence of FruA/FruA*, whereas in the mrpC mutant a small amount of FruA/FruA* may positively regulate gene expression in the absence of MrpC. This hypothesis can be tested by measuring gene expression in an mrpC fruA double mutant.

Other genes with the c- f- m- regulatory pattern (i.e., those below the diagonal in the bottom left graph of Additional file 7) are more strongly downregulated in the mrpC mutant than in the fruA mutant (e.g., fadIJ). For these genes, the complete absence of MrpC and FruA/FruA* (or a small amount of the latter) in the mrpC mutant affects gene expression differently than the absence of just FruA (not MrpC) in the fruA mutant. Here too, measuring gene expression in an mrpC fruA mutant may help elucidate the regulatory mechanism. In the case of the fadIJ operon, the fadI transcript level was very low in csgA, fruA, and mrpC single mutants at 24 h PS, and FruA and MrpC bound to the promoter region both independently and with weak cooperativity [70], so it is likely that independent binding of MrpC in the fruA mutant weakly stimulates transcription. Understanding the regulatory mechanisms of other genes with the c- f- m- pattern will likewise require analysis of FruA and MrpC binding to the promoter regions. We note that FruA ChIP-seq could be performed and analyzed in combination with the existing MrpC ChIP-seq [77], Cappable-seq data [57], and our RNA-seq data to predict regulatory mechanisms genome-wide not only for genes with the c- f- m- pattern, but for genes with other regulatory patterns as well (see below).

A novel finding of our study is negative coregulation by CsgA, FruA, and MrpC (i.e., the c + f + m + pattern). Interestingly, the 157 genes with this pattern (Fig. 3b, right) exhibit much more similar strength of regulation in the fruA and mrpC mutants than the 382 genes with the c- f- m- pattern (Additional file 7, bottom graphs). We infer that FruA/FruA* and MrpC act independently as negative regulators less frequently than as positive regulators. This may reflect strongly cooperative DNA binding of the two proteins to repress transcription. In agreement, we note that the 421 genes with the f + m + pattern (Fig. 3b, right) also exhibit very similar strength of regulation in the fruA and mrpC mutants (Additional file 7, top right graph). Since these genes were < 2-fold upregulated in the csgA mutant, we infer that unactivated FruA (not activated FruA*, whose formation requires CsgA) binds cooperatively with MrpC to repress these genes.

A surprising finding of our study is hundreds of genes differentially expressed only in one of the three mutants (Fig. 3b). Previous work had uncovered only the mrpAB operon and the sdeK gene as uniquely downregulated in the mrpC mutant (Table 1 and Additional file 2). Our RNA-seq data adds 375 genes positively regulated and 483 genes negatively regulated by MrpC independently of CsgA and FruA (Fig. 3b). These are likely a combination of direct and indirect effects of MrpC acting independently as an activator and a repressor of transcription. ChIP-seq identified 1608 putative MrpC binding sites at 18 h PS, but in 13 of 15 cases tested there was evidence for cooperative binding with FruA [77]. To estimate how frequently MrpC acts independently as a direct regulator, we determined the number of genes uniquely regulated in the mrpC mutant and with an MrpC-binding site [77] between − 200 to + 200 bp relative to a TSS [57]. We found only 75 down- and 62 up-regulated genes (Additional file 13), suggesting that MrpC acts independently as a direct regulator infrequently. Among the downregulated genes, we found mrpA but not mrpB (Additional files 2 and 13), which is cotranscribed with mrpA [58], illustrating that our gene number estimates are low since downstream genes in operons are not included. Even so, most of the independent effects of MrpC appear to be indirect.

Presumably, all the effects of CsgA are indirect, since CsgA is not known to bind DNA. Finding 347 down- and 626 up-regulated genes unique to the csgA mutant (Fig. 3b) shows the profound effects of CsgA on the developmental gene regulatory network apart from FruA and MrpC. How CsgA exerts these effects emerges as an important question from our study.

It was most surprising to find 158 down- and 181 up-regulated genes unique to the fruA mutant (Fig. 3b), given the strong dependence of FruA on MrpC [56, 65] (Additional file 2). However, as mentioned in the preceding paragraphs, a small amount of FruA/FruA* in the mrpC mutant could explain its differences from the fruA mutant, which completely lacks FruA/FruA*. In the csgA mutant, FruA accumulates normally early in development and then persists rather than decreasing as in the WT strain [64], but presumably FruA is not converted to FruA*. Previously, expression of only one gene had been shown to depend on FruA, but not on CsgA. The dofA gene was downregulated in a fruA mutant [78], but unchanged in a csgA mutant [112] (Additional file 2) and FruA appears to directly activate dofA transcription by binding to two sites between − 82 and − 42 relative to the TSS [113]. Our RNA-seq data suggests that MrpC is not involved in dofA regulation (Additional file 2), so we infer that FruA (not FruA*) acts independently to activate dofA transcription. To estimate the frequency of this regulatory mechanism among the 338 other genes uniquely regulated in the fruA mutant, FruA ChIP-seq could be performed and analyzed in combination with the Cappable-seq data reported recently [57].

Our RNA-seq data revealed that 617 genes are counterregulated in the csgA, fruA, and mrpC mutants at 24 h PS. Our data agree with recent reports of exoA, exoL, nfsA, pmxA, pkn1, MXAN_6957, and MXAN_7024 being counterregulated, though our data did not agree with the previously reported counterregulation of spiA and MXAN_2902 (Additional file 2). Importantly, the developmental timepoint sampled can impact the comparison between transcript levels in the WT strain and mutants. For example, exoA exhibits the c- f + m + regulatory pattern at 24 h, but shows the c- pattern (i.e., loses counterregulation) at 30 h [70]. The developmental conditions and the particular strains may also influence the timing of developmental regulatory patterns. For example, our 24-h data better matched 12-h data from other studies in several instances (Additional file 2). In any case, counterregulation suggests regulatory complexity, and consistent with this notion, the regulation of exoA, exoL, and nfsA was shown to involve not only CsgA, FruA, and MrpC, but transcription factor Nla6-P and CRISPR-Cas component DevI as well [70]. Further investigation of counterregulated genes promises to be intriguing.

In summary, our work highlights diverse regulatory roles of CsgA, FruA, and MrpC late during M. xanthus development in submerged culture. All three are extensive negative and positive regulators, and we observed nearly every possible pattern of coregulation, unique regulation, and counterregulation.

CsgA, FruA, and MrpC broadly regulate genes involved in cellular signaling, protein synthesis, energetics, and envelope function during development

Prior timecourse transcriptomic studies showed that cellular signaling, translation, and metabolism are broadly altered during M. xanthus development [43, 44, 78], and our results show that much of this regulation involves CsgA, FruA, and MrpC. Our GO term overrepresentation analysis showed that CsgA is a negative regulator of many signaling pathways and MrpC positively regulates a subset of those pathways (Fig. 5). CsgA, FruA, and MrpC are positive regulators of many genes necessary for ribosome biogenesis and translation, based on both our GO term and KEGG pathway gene set enrichment analyses (Fig. 5a and Additional files 8 and 10). Interestingly, this positive regulation opposes the initial downregulation of genes involved in translation during development [43, 78], as discussed further below. Our GO term analysis identified CsgA and FruA as positive regulators of genes involved in energy production and conversion (Fig. 5a), and our KEGG pathway analysis revealed a strong emphasis on regulation of cellular energetics and metabolism (Additional file 10). Developmental timecourse studies showed that after an initial phase of downregulation of genes involved in energy production and many metabolic pathways, transcriptional rewiring of metabolic pathways occurs, presumably to utilize changing carbohydrate, lipid, and amino acid resources as development proceeds [43, 44, 78]. Our results show that CsgA, FruA, and MrpC participate broadly in the rewiring of metabolic and energy production and conversion pathways at 24 h PS, when rods have formed mounds. Some rods have shortened by 24 h, although round spores have not yet formed [19, 114]. Cell shortening requires envelope remodeling. Both our GO term and KEGG pathway analyses indicated that CsgA in particular regulates genes implicated in cell envelope structure and function (Fig. 5b and Additional files 8 and 10). Consistently, transcripts from genes involved in cell wall/membrane/envelope biogenesis increased by 24 h [44]. Some of these genes (e.g., nfsA-C, nfsF-H, and pkn1) were upregulated in a fruA mutant relative to a WT strain [78]. Our RNA-seq data agree with the prior studies and show that CsgA and MrpC also regulate envelope function during development (Additional file 2).

Control of ribosome biogenesis, lipid signals, and peptidoglycan biosynthesis by CsgA, FruA, and MrpC during development

Our RNA-seq data suggests that CsgA, FruA, and MrpC regulate transcript levels at 24 h PS to decrease (p)ppGpp levels (Fig. 7), thus diminishing the stringent response to starvation and promoting ribosome biogenesis for developmental gene expression. A simple prediction of this model is that (p)ppGpp levels would be elevated in a csgA mutant relative to the WT strain. However, (p)ppGpp levels were 2- to 4-fold lower in a csgA mutant than in a WT strain at 24 h PS, and enzymatic activity of CsgA appeared to induce the stringent response [108], which raises (p)ppGpp levels. These findings suggest that the role of CsgA in the regulatory network is more complex than can be predicted from the effects of a csgA mutation on transcript levels. For example, posttranscriptional regulatory mechanisms may be important. In any case, positive and negative regulation of the stringent response and (p)ppGpp levels must be balanced during development to ensure that some cells within mounds form spores.

Our results suggest that CsgA, FruA, and MrpC positively regulate genes whose products generate isovaleryl-CoA, iso15:0, and TG-1, which have been implicated in the formation of lipid bodies and lipid signals. Coincident with cell shortening early in development, membrane phospholipids appear to be converted to lipid bodies composed primarily of triacylglycerols (many derived from iso15:0) and their ether derivatives (with TG-1 being the most abundant) [115, 116]. The lipid bodies appear to be carbon storage organelles for later generation of energy to complete sporulation. Mutations in csgA and fruA severely impair cell shortening and lipid body formation, whereas a deletion of mrpC causes milder defects [115]. These defects likely result from failure to upregulate genes for BCAA metabolism (Fig. 8). In the WT strain, TG-1 and related iso15:0-derived ether lipids in lipid bodies may be released from lysing cells to function as developmental signals [18,19,20]. The levels of these lipids are lower in csgA and fruA mutants than in the WT strain at 24 h [116, 117] and 72 h PS [18], but the levels in an mrpC have not been reported. Our data suggest that ether lipids may be lower in the mrpC mutant than in the WT strain, since MrpC, like CsgA and FruA, positively regulates genes involved in many steps of BCAA metabolism (Fig. 8). However, elbD (MXAN_1528) is upregulated in the mrpC mutant (Fig. 8). ElbD is primarily responsible for ether lipid biosynthesis, including TG-1, at the end of the pathway [118]. Increased elbD expression in the mrpC mutant might restore ether lipid levels, at least partially. Another intriguing exception to the downregulation of genes involved in BCAA metabolism in the csgA, fruA, and mrpC mutants is the upregulation of MXAN_2987 in the csgA mutant (Fig. 8). MXAN_2987 codes for BCAA transaminase, which catalyzes oxidative deamination of BCAAs at the beginning of the pathway. We speculate that CsgA negatively regulates the first step of the pathway to conserve BCAAs for protein synthesis and/or for allosteric activation of downstream enzymes in the pathway (e.g., as observed for leucine activation of glutamine dehydrogenase) [119].

We found that MrpC negatively regulates genes involved in the initial and penultimate steps of PG biosynthesis, but CsgA positively regulates genes involved in intermediate steps (Fig. 9). The negative regulation by MrpC is an appropriate response to starvation-induced growth arrest, since elongation of the PG sacculus ceases. Presumably, the levels of PG components increase in the starving WT strain. This may further increase csgA expression beyond that caused by starvation, since addition of the PG components N-acetylglucosamine (GlcNAc), diaminopimelate (DAP), and D-alanine (D-Ala) to cells induces csgA expression even in the presence of nutrients [29]. The starvation-dependent increase in csgA expression may be a direct effect of MrpC, since an MrpC-binding site centered at -407 relative to the TSS was detected by ChIP-seq [77]. Both the CsgA protein [58] and the csgA (MXAN_1294) transcript levels are reduced in the mrpC mutant (Additional file 5). Hence, in the starving WT strain, PG components and MrpC may elevate the CsgA level, which would positively regulate genes involved in intermediate steps of PG synthesis based on our RNA-seq data (Fig. 9). This positive regulation by CsgA, together with the negative regulation by MrpC of genes involved in the initial and penultimate steps of PG biosynthesis, may lead to the accumulation of the lipid II precursor as well as the smaller PG components. Later in development, as rods become spores, PG degradation presumably releases PG intermediates. We speculate that PG degradation and earlier regulation by MrpC and CsgA of genes involved in PG biosynthesis create a store of PG intermediates in spores that are used to rebuild PG during germination. In support of this hypothesis, chemically-induced myxospores treated with nutrients elongate into rods in the presence of fosfomycin [120], an antibiotic that inhibits MurA at the beginning of the PG biosynthetic pathway (Fig. 9a) and thus prevents synthesis of PG intermediates. Similar experiments with starvation-induced fruiting body myxospores have not been reported. Our hypothesis predicts that starvation-induced myxospores treated with nutrients would also elongate into rods in the presence of fosfomycin.

Conclusions

This work provides a wealth of information about the impact of mutations in csgA, fruA, and mrpC on the M. xanthus developmental transcriptome. The developmental process is a delicate balance, as the lytic fate of the majority of rods, the torpor of hardy spores in fruiting bodies, and the alternate fate of peripheral rods that remain outside fruiting bodies, all require great inputs of energy. An intricate network of signals and both positive and negative regulators of gene expression are required to ensure proper temporal and spatial determination of cell fate. Splendid work has been undertaken to study the regulation of individual genes as well as temporal progression of the developmental transcriptome. Our work expands on previous findings and will serve as a springboard for additional targeted investigations into the considerable gaps in knowledge that remain.

Data availability

The sequence data has been deposited in the NCBI SRA database and can be accessed at https://www.ncbi.nlm.nih.gov/sra/PRJNA1177521. The BioProject accession number is PRJNA1177521.

References

  1. Yang Z, Higgs P, editors. Myxobacteria: genomics, cellular and molecular biology. Norfolk, UK: Caister Academic; 2014.

    Google Scholar 

  2. Gerth K, Irschik H, Reichenbach H, Trowitzsch W. The Myxovirescins, a family of antibiotics from Myxococcus virescens (Myxobacterales). J Antibiot (Tokyo). 1982;35(11):1454–9.

    Article  CAS  PubMed  Google Scholar 

  3. Altmann K-H, Kinghorn AD, Höfle G, Müller R, Prantz K. The epothilones: an outstanding family of anti-tumor agents: from soil to the clinic. Springer Science & Business Media; 2009.

  4. Goldman BS, Nierman WC, Kaiser D, Slater SC, Durkin AS, Eisen J, et al. Evolution of sensory complexity recorded in a myxobacterial genome. Proc Natl Acad Sci USA. 2006;103(41):15200–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Whitworth DE, Sydney N, Radford EJ. Myxobacterial genomics and Post-Genomics: A review of genome biology, genome sequences and related ‘omics studies. Microorganisms. 2021;9(10):2143.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Kroos L. Highly signal-responsive gene regulatory network governing Myxococcus development. Trends Genet. 2017;33:3–15.

    Article  CAS  PubMed  Google Scholar 

  7. Bretl DJ, Kirby JR. Molecular mechanisms of signaling in Myxococcus xanthus development. J Mol Biol. 2016;428:3805–30.

    Article  CAS  PubMed  Google Scholar 

  8. Munoz-Dorado J, Marcos-Torres FJ, Garcia-Bravo E, Moraleda-Munoz A, Perez J. Myxobacteria: moving, killing, feeding, and surviving together. Front Microbiol. 2016;7:781.

    Article  PubMed  PubMed Central  Google Scholar 

  9. Lee B, Holkenbrink C, Treuner-Lange A, Higgs PI. Myxococcus xanthus developmental cell fate production: heterogeneous accumulation of developmental regulatory proteins and reexamination of the role of MazF in developmental Lysis. J Bacteriol. 2012;194(12):3058–68.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. O’Connor KA, Zusman DR. Reexamination of the role of autolysis in the development of Myxococcus xanthus. J Bacteriol. 1988;170(9):4103–12.

    Article  PubMed  PubMed Central  Google Scholar 

  11. O’Connor KA, Zusman DR. Behavior of peripheral rods and their role in the life cycle of Myxococcus xanthus. J Bacteriol. 1991;173(11):3342–55.

    Article  PubMed  PubMed Central  Google Scholar 

  12. O’Connor KA, Zusman DR. Development in Myxococcus xanthus involves differentiation into two cell types, peripheral rods and spores. J Bacteriol. 1991;173(11):3318–33.

    Article  PubMed  PubMed Central  Google Scholar 

  13. Hagen DC, Bretscher AP, Kaiser D. Synergism between morphogenetic mutants of Myxococcus xanthus. Dev Biol. 1978;64:284–96.

    Article  CAS  PubMed  Google Scholar 

  14. Kuspa A, Plamann L, Kaiser D. Identification of heat-stable A-factor from Myxococcus xanthus. J Bacteriol. 1992;174:3319–26.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Kuspa A, Plamann L, Kaiser D. A-signalling and the cell density requirement for Myxococcus xanthus development. J Bacteriol. 1992;174:7360–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Plamann L, Kuspa A, Kaiser D. Proteins that rescue A-signal-defective mutants of Myxococcus xanthus. J Bacteriol. 1992;174:3311–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Toal D, Clifton S, Roe B, Downard J. The Esg locus of Myxococcus xanthus encodes the E1a and E1b subunits of a branced-chain Keto acid dehydrogenase. Mol Microbiol. 1995;16:177–89.

    Article  CAS  PubMed  Google Scholar 

  18. Ring MW, Schwar G, Thiel V, Dickschat JS, Kroppenstedt RM, Schulz S, et al. Novel iso-branched ether lipids as specific markers of developmental sporulation in the myxobacterium Myxococcus xanthus. J Biol Chem. 2006;281(48):36691–700.

    Article  CAS  PubMed  Google Scholar 

  19. Bhat S, Ahrendt T, Dauth C, Bode HB, Shimkets LJ. Two lipid signals guide fruiting body development of Myxococcus xanthus. MBio. 2014;5(1):e00939–13.

    Article  PubMed  PubMed Central  Google Scholar 

  20. Ahrendt T, Dauth C, Bode HB. An iso-15:0 O-alkylglycerol moiety is the key structure of the E-signal in Myxococcus xanthus. Microbiology. 2016;162(1):138–44.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Kim SK, Kaiser D. C-factor: a cell-cell signaling protein required for fruiting body morphogenesis of M. xanthus. Cell. 1990;61:19–26.

    Article  CAS  PubMed  Google Scholar 

  22. Lobedanz S, Sogaard-Andersen L. Identification of the C-signal, a contact-dependent morphogen coordinating multiple developmental responses in Myxococcus xanthus. Genes Dev. 2003;17(17):2151–61.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Boynton TO, Shimkets LJ, CsgA. Drosophila Sniffer and human HSD17B10 are cardiolipin phospholipases. Genes Dev. 2015;29:1903-14.

  24. Konovalova A, Wegener-Feldbrugge S, Sogaard-Andersen L. Two intercellular signals required for fruiting body formation in Myxococcus xanthus act sequentially but non-hierarchically. Mol Microbiol. 2012;86(1):65–81.

    Article  CAS  PubMed  Google Scholar 

  25. Kroos L, Kaiser D. Expression of many developmentally regulated genes in Myxococcus depends on a sequence of cell interactions. Genes Dev. 1987;1:840–54.

    Article  CAS  PubMed  Google Scholar 

  26. Kuspa A, Kroos L, Kaiser D. Intercellular signaling is required for developmental gene expression in Myxococcus xanthus. Dev Biol. 1986;117:267–76.

    Article  CAS  PubMed  Google Scholar 

  27. Kim SK, Kaiser D. C-factor has distinct aggregation and sporulation thresholds during Myxococcus development. J Bacteriol. 1991;173:1722–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Kruse T, Lobedanz S, Berthelsen NM, Sogaard-Andersen L. C-signal: a cell surface-associated morphogen that induces and co-ordinates multicellular fruiting body morphogenesis and sporulation in Myxococcus xanthus. Mol Microbiol. 2001;40(1):156–68.

    Article  CAS  PubMed  Google Scholar 

  29. Li S-F, Lee B, Shimkets LJ. CsgA expression entrains Myxococcus xanthus development. Genes Dev. 1992;6:401–10.

    Article  CAS  PubMed  Google Scholar 

  30. Kim SK, Kaiser D. Cell motility is required for the transmission of C-factor, an intercellular signal that coordinates fruiting body morphogenesis of Myxococcus xanthus. Genes Dev. 1990;4:896–905.

    Article  CAS  PubMed  Google Scholar 

  31. Kim SK, Kaiser D. Cell alignment required in differentiation of Myxococcus xanthus. Science. 1990;249:926–8.

    Article  CAS  PubMed  Google Scholar 

  32. Kroos L, Hartzell P, Stephens K, Kaiser D. A link between cell movement and gene expression argues that motility is required for cell-cell signaling during fruiting body development. Genes Dev. 1988;2:1677–85.

    Article  CAS  PubMed  Google Scholar 

  33. Gomez-Santos N, Glatter T, Koebnik R, Swiatek-Polatynska MA, Sogaard-Andersen L. A TonB-dependent transporter is required for secretion of protease PopC across the bacterial outer membrane. Nat Commun. 2019;10(1):1360.

    Article  PubMed  PubMed Central  Google Scholar 

  34. Rajagopalan R, Kroos L. Nutrient-regulated proteolysis of MrpC halts expression of genes important for commitment to sporulation during Myxococcus xanthus development. J Bacteriol. 2014;196:2736–47.

    Article  PubMed  PubMed Central  Google Scholar 

  35. Manoil C, Kaiser D. Accumulation of guanosine tetraphosphate and guanosine pentaphosphate in Myxococcus xanthus during starvation and myxospore formation. J Bacteriol. 1980;141:297–304.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Manoil C, Kaiser D. Guanosine pentaphosphate and guanosine tetraphosphate accumulation and induction of Myxococcus xanthus fruiting body development. J Bacteriol. 1980;141(1):305–15.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Skotnicka D, Smaldone GT, Petters T, Trampari E, Liang J, Kaever V, et al. A minimal threshold of c-di-GMP is essential for fruiting body formation and sporulation in Myxococcus xanthus. PLoS Genet. 2016;12(5):e1006080.

    Article  PubMed  PubMed Central  Google Scholar 

  38. Boutte CC, Crosson S. Bacterial lifestyle shapes stringent response activation. Trends Microbiol. 2013;21(4):174–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Singer M, Kaiser D. Ectopic production of guanosine penta-and tetraphosphate can initiate early developmental gene expression in Myxococcus xanthus. Genes Dev. 1995;9:1633–44.

    Article  CAS  PubMed  Google Scholar 

  40. Harris BZ, Kaiser D, Singer M. The guanosine nucleotide (p)ppGpp initiates development and A-factor production in Myxococcus xanthus. Genes Dev. 1998;12(7):1022–35.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Crawford EW Jr., Shimkets LJ. The Myxococcus xanthus socE and csgA genes are regulated by the stringent response. Mol Microbiol. 2000;37(4):788–99.

    Article  CAS  PubMed  Google Scholar 

  42. Konovalova A, Lobach S, Sogaard-Andersen L. A RelA-dependent two-tiered regulated proteolysis cascade controls synthesis of a contact-dependent intercellular signal in Myxococcus xanthus. Mol Microbiol. 2012;84(2):260–75.

    Article  CAS  PubMed  Google Scholar 

  43. Munoz-Dorado J, Moraleda-Munoz A, Marcos-Torres FJ, Contreras-Moreno FJ, Martin-Cuadrado AB, Schrader JM, et al. Transcriptome dynamics of the Myxococcus xanthus multicellular developmental program. eLife. 2019;8:e50374.

    Article  PubMed  PubMed Central  Google Scholar 

  44. Sharma G, Yao AI, Smaldone GT, Liang J, Long M, Facciotti MT, et al. Global gene expression analysis of the Myxococcus xanthus developmental time course. Genomics. 2021;113(1 Pt 1):120–34.

    Article  CAS  PubMed  Google Scholar 

  45. Giglio KM, Caberoy N, Suen G, Kaiser D, Garza AG. A cascade of coregulating enhancer binding proteins initiates and propagates a multicellular developmental program. Proc Natl Acad Sci USA. 2011;108:E431–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Diodati ME, Ossa F, Caberoy NB, Jose IR, Hiraiwa W, Igo MM, et al. Nla18, a key regulatory protein required for normal growth and development of Myxococcus xanthus. J Bacteriol. 2006;188:1733–43.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Ossa F, Diodati ME, Caberoy NB, Giglio KM, Edmonds M, Singer M, et al. The Myxococcus xanthus Nla4 protein is important for expression of stringent response-associated genes, ppGpp accumulation, and fruiting body development. J Bacteriol. 2007;189(23):8474–83.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Giglio KM, Zhu C, Klunder C, Kummer S, Garza AG. The enhancer binding protein Nla6 regulates developmental genes that are important for Myxococcus xanthus sporulation. J Bacteriol. 2015;197(7):1276–87.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Sarwar Z, Garza AG. Two-Component signal transduction systems that regulate the temporal and spatial expression of Myxococcus xanthus sporulation genes. J Bacteriol. 2016;198(3):377–85.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Ma M, Garza AG, Lemon DJ, Caro EA, Ritchie L, Ryan C, et al. Identifying the gene regulatory network of the starvation-induced transcriptional activator Nla28. J Bacteriol. 2022;204(12):e0026522.

    Article  PubMed  Google Scholar 

  51. Gronewold TM, Kaiser D. The act Operon controls the level and time of C-signal production for Myxococcus xanthus development. Mol Microbiol. 2001;40(3):744–56.

    Article  CAS  PubMed  Google Scholar 

  52. Gronewold TM, Kaiser D. act operon control of developmental gene expression in Myxococcus xanthus. J Bacteriol. 2002;184(4):1172-9.

  53. Jelsbak L, Givskov M, Kaiser D. Enhancer-binding proteins with an FHA domain and the sigma54 regulon in Myxococcus xanthus fruiting body development. Proc Natl Acad Sci USA. 2005;102:3010–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Volz C, Kegler C, Muller R. Enhancer binding proteins act as hetero-oligomers and link secondary metabolite production to myxococcal development, motility, and predation. Chem Biol. 2012;19(11):1447–59.

    Article  CAS  PubMed  Google Scholar 

  55. McLaughlin PT, Bhardwaj V, Feeley BE, Higgs PI. MrpC, a CRP/Fnr homolog, functions as a negative autoregulator during the Myxococcus xanthus multicellular developmental program. Mol Microbiol. 2018;109:245–61.

    Article  CAS  PubMed  Google Scholar 

  56. Ueki T, Inouye S. Identification of an activator protein required for the induction of fruA, a gene essential for fruiting body development in Myxococcus xanthus. Proc Natl Acad Sci USA. 2003;100(15):8782–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Kuzmich S, Blumenkamp P, Meier D, Szadkowski D, Goesmann A, Becker A, et al. CRP-Like transcriptional regulator MrpC curbs c-di-GMP and 3’,3’-cGAMP nucleotide levels during development in Myxococcus xanthus. mBio. 2022;13:e00044–22.

    Article  CAS  PubMed Central  Google Scholar 

  58. Sun H, Shi W. Genetic studies of mrp, a locus essential for cellular aggregation and sporulation of Myxococcus xanthus. J Bacteriol. 2001;183(16):4786–95.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Sun H, Shi W. Analyses of Mrp genes during Myxococcus xanthus development. J Bacteriol. 2001;183(23):6733–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Higgs PI, Jagadeesan S, Mann P, Zusman DR. EspA, an orphan hybrid histidine protein kinase, regulates the timing of expression of key developmental proteins of Myxococcus xanthus. J Bacteriol. 2008;190:4416–26.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Hoang Y, Kroos L. Ultrasensitive response of developing Myxococcus xanthus to the addition of nutrient medium correlates with the level of MrpC. J Bacteriol. 2018;200:e00456–18.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Ogawa M, Fujitani S, Mao X, Inouye S, Komano T. FruA, a putative transcription factor essential for the development of Myxococcus xanthus. Mol Microbiol. 1996;22(4):757–67.

    Article  CAS  PubMed  Google Scholar 

  63. Mittal S, Kroos L. A combination of unusual transcription factors binds cooperatively to control Myxococcus xanthus developmental gene expression. Proc Natl Acad Sci USA. 2009;106:1965–70.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Ellehauge E, Norregaard-Madsen M, Sogaard-Andersen L. The FruA signal transduction protein provides a checkpoint for the Temporal co-ordination of intercellular signals in Myxococcus xanthus development. Mol Microbiol. 1998;30:807–17.

    Article  CAS  PubMed  Google Scholar 

  65. Saha S, Patra P, Igoshin O, Kroos L. Systematic analysis of the Myxococcus xanthus developmental gene regulatory network supports posttranslational regulation of FruA by C-signaling. Mol Microbiol. 2019;111:1732–52.

    Article  CAS  PubMed  Google Scholar 

  66. Campbell A, Viswanathan P, Barrett T, Son B, Saha S, Kroos L. Combinatorial regulation of the Dev Operon by MrpC2 and FruA during Myxococcus xanthus Development. J Bacteriol. 2015;197:240–51.

    Article  PubMed  Google Scholar 

  67. Lee J, Son B, Viswanathan P, Luethy P, Kroos L. Combinatorial regulation of FmgD by MrpC2 and FruA during Myxococcus xanthus development. J Bacteriol. 2011;193:1681–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Mittal S, Kroos L. Combinatorial regulation by a novel arrangement of FruA and MrpC2 transcription factors during Myxococcus xanthus development. J Bacteriol. 2009;191:2753–63.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  69. Son B, Liu Y, Kroos L. Combinatorial regulation by MrpC2 and FruA involves three sites in the FmgE promoter region during Myxococcus xanthus development. J Bacteriol. 2011;193:2756–66.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Saha S, Kroos L. Regulation of late-acting operons by three transcription factors and a CRISPR-Cas component during Myxococcus xanthus development. Mol Microbiol. 2024;121:1002–20.

    Article  CAS  PubMed  Google Scholar 

  71. Viswanathan K, Viswanathan P, Kroos L. Mutational analysis of the Myxococcus xanthus omega4406 promoter region reveals an upstream negative regulatory element that mediates C-signal dependence. J Bacteriol. 2006;188:515–24.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  72. Viswanathan P, Kroos L. cis elements necessary for developmental expression of a Myxococcus xanthus gene that depends on C signaling. J Bacteriol. 2003;185(4):1405-14.

  73. Viswanathan P, Ueki T, Inouye S, Kroos L. Combinatorial regulation of genes essential for Myxococcus xanthus development involves a response regulator and a LysR-type regulator. Proc Natl Acad Sci USA. 2007;104:7969–74.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  74. Yoder D, Kroos L. Mutational analysis of the Myxococcus xanthus omega4400 promoter region provides insight into developmental gene regulation by C signaling. J Bacteriol. 2004;186:661–71.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Yoder D, Kroos L. Mutational analysis of the Myxococcus xanthus omega4499 promoter region reveals shared and unique properties in comparison with other C-signal-dependent promoters. J Bacteriol. 2004;186:3766–76.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  76. Yoder-Himes D, Kroos L. Regulation of the Myxococcus xanthus C-signal-dependent omega4400 promoter by the essential developmental protein FruA. J Bacteriol. 2006;188:5167–76.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  77. Robinson M, Son B, Kroos L. Transcription factor MrpC binds to promoter regions of many developmentally-regulated genes in Myxococcus xanthus. BMC Genomics. 2014;15:1123.

    Article  PubMed  PubMed Central  Google Scholar 

  78. McLoon AL, Boeck ME, Bruckskotten M, Keyel AC, Sogaard-Andersen L. Transcriptomic analysis of the Myxococcus xanthus FruA Regulon, and comparative developmental transcriptomic analysis of two fruiting body forming species, Myxococcus xanthus and Myxococcus stipitatus. BMC Genomics. 2021;22(1):784.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  79. Kaiser D. Social gliding is correlated with the presence of pili in Myxococcus xanthus. Proc Natl Acad Sci USA. 1979;76:5952–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  80. Shimkets LJ, Asher SJ. Use of recombination techniques to examine the structure of the Csg locus of Myxococcus xanthus. Mol Gen Genet. 1988;211:63–71.

    Article  CAS  PubMed  Google Scholar 

  81. Kuner JM, Kaiser D. Fruiting body morphogenesis in submerged cultures of Myxococcus xanthus. J Bacteriol. 1982;151:458–61.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  82. Johnson BK, Scholz MB, Teal TK, Abramovitch RB. SPARTA: simple program for automated reference-based bacterial RNA-seq transcriptome analysis. BMC Bioinformatics. 2016;17:1–4.

    Article  Google Scholar 

  83. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for illumina sequence data. Bioinformatics. 2014;30(15):2114–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  84. Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10(3):R25.

    Article  PubMed  PubMed Central  Google Scholar 

  85. Anders S, Pyl PT, Huber W. HTSeq–a Python framework to work with high-throughput sequencing data. Bioinformatics. 2015;31(2):166–9.

    Article  CAS  PubMed  Google Scholar 

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

    Article  PubMed  PubMed Central  Google Scholar 

  87. Team R, RStudio, Boston MA. 2020;http://www.rstudio.com/.

  88. Kolde R, pheatmap. Pretty Heatmaps. R package version 1012. 2018;https://github.com/raivokolde/pheatmap

  89. Wickham H. Ggplot: elegant graphics for data analysis. New York: Springer-; 2016.

    Book  Google Scholar 

  90. Hunter J, Matplotlib. A 2D graphics environment. Comput Sci Eng. 2007;9(3):90–5.

    Article  Google Scholar 

  91. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  92. Ge SX, Jung D, Yao R. ShinyGO: a graphical gene-set enrichment tool for animals and plants. Bioinformatics. 2020;36(8):2628–9.

    Article  CAS  PubMed  Google Scholar 

  93. Sherman BT, Hao M, Qiu J, Jiao X, Baseler MW, Lane HC, et al. DAVID: a web server for functional enrichment analysis and functional annotation of gene lists (2021 update). Nucleic Acids Res. 2022;50(W1):W216–21.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  94. Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44–57.

    Article  CAS  PubMed  Google Scholar 

  95. Yu G, Wang LG, Han Y, He QY. ClusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  96. Luo W, Brouwer C. Pathview: an R/Bioconductor package for pathway-based data integration and visualization. Bioinformatics. 2013;29(14):1830–1.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  97. Kassambara A. ggpubr: ‘ggplot2’ Based Publication Ready Plots. R package version 0.6.0. https://rpkgs.datanovia.com/ggpubr/. 2023.

  98. Marini F, Binder H. PcaExplorer: an R/Bioconductor package for interacting with RNA-seq principal components. BMC Bioinformatics. 2019;20(1):331.

    Article  PubMed  PubMed Central  Google Scholar 

  99. Bullock HA, Shen H, Boynton TO, Shimkets LJ. Fatty acid oxidation is required for Myxococcus xanthus development. J Bacteriol. 2018;200(10):e00572–17.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  100. Licking E, Gorski L, Kaiser D. A common step for changing cell shape in fruiting body and starvation-independent sporulation of Myxococcus xanthus. J Bacteriol. 2000;182(12):3553–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  101. Muller FD, Schink CW, Hoiczyk E, Cserti E, Higgs PI. Spore formation in Myxococcus xanthus is tied to cytoskeleton functions and polysaccharide spore coat deposition. Mol Microbiol. 2012;83(3):486–505.

    Article  PubMed  Google Scholar 

  102. Ueki T, Inouye S. Identification of a gene involved in polysaccharide export as a transcription target of FruA, an essential factor for Myxococcus xanthus development. J Biol Chem. 2005;280(37):32279–84.

    Article  CAS  PubMed  Google Scholar 

  103. Muller FD, Treuner-Lange A, Heider J, Huntley SM, Higgs PI. Global transcriptome analysis of spore formation in Myxococcus xanthus reveals a locus necessary for cell differentiation. BMC Genomics. 2010;11:264.

    Article  PubMed  PubMed Central  Google Scholar 

  104. Browning DF, Busby SJ. Local and global regulation of transcription initiation in bacteria. Nat Rev Microbiol. 2016;14(10):638–50.

    Article  CAS  PubMed  Google Scholar 

  105. Browning DF, Butala M, Busby SJW. Bacterial transcription factors: regulation by pick N mix. J Mol Biol. 2019;431(20):4067–77.

    Article  CAS  PubMed  Google Scholar 

  106. Lee DJ, Minchin SD, Busby SJ. Activating transcription in bacteria. Annu Rev Microbiol. 2012;66:125–52.

    Article  CAS  PubMed  Google Scholar 

  107. Otasek D, Morris JH, Boucas J, Pico AR, Demchak B. Cytoscape automation: empowering workflow-based network analysis. Genome Biol. 2019;20(1):185.

    Article  PubMed  PubMed Central  Google Scholar 

  108. Crawford EW, Shimkets LJ. The stringent response in Myxococcus xanthus is regulated by SocE and the CsgA C-signaling protein. Genes Dev. 2000;14(4):483–92.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  109. Garcia-Moreno D, Polanco MC, Navarro-Aviles G, Murillo FJ, Padmanabhan S, Elias-Arnanz M. A vitamin B12-based system for conditional expression reveals DksA to be an essential gene in Myxococcus xanthus. J Bacteriol. 2009;191(9):3108–19.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  110. Bode HB, Ring MW, Schwar G, Kroppenstedt RM, Kaiser D, Muller R. 3-Hydroxy-3-methylglutaryl-coenzyme A (CoA) synthase is involved in biosynthesis of isovaleryl-CoA in the myxobacterium Myxococcus xanthus during fruiting body formation. J Bacteriol. 2006;188(18):6524–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  111. Bode HB, Ring MW, Schwar G, Altmeyer MO, Kegler C, Jose IR, et al. Identification of additional players in the alternative biosynthesis pathway to isovaleryl-CoA in the myxobacterium Myxococcus xanthus. ChemBioChem. 2009;10(1):128–40.

    Article  CAS  PubMed  Google Scholar 

  112. Horiuchi T, Akiyama T, Inouye S, Komano T. Analysis of DofA, a fruA-dependent developmental gene, and its homologue, DofB, in Myxococcus xanthus. J Bacteriol. 2002;184(24):6803–10.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  113. Ueki T, Inouye S. Activation of a development-specific gene, dofA, by FruA, an essential transcription factor for development of Myxococcus xanthus. J Bacteriol. 2005;187(24):8504–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  114. Hoang Y, Franklin JL, Dufour YS, Kroos L. Cell density, alignment, and orientation correlate with C-signal-dependent gene expression during Myxococcus xanthus development. Proc Natl Acad Sci USA. 2021;118(45):e2111706118.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  115. Bhat S, Boynton TO, Pham D, Shimkets LJ. Fatty acids from membrane lipids become incorporated into lipid bodies during Myxococcus xanthus differentiation. PLoS ONE. 2014;9(6):e99622.

    Article  PubMed  PubMed Central  Google Scholar 

  116. Hoiczyk E, Ring MW, McHugh CA, Schwar G, Bode E, Krug D, et al. Lipid body formation plays a central role in cell fate determination during developmental differentiation of Myxococcus xanthus. Mol Microbiol. 2009;74(2):497–517.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  117. Ahrendt T, Wolff H, Bode HB. Neutral and phospholipids of the Myxococcus xanthus lipodome during fruiting body formation and germination. Appl Environ Microbiol. 2015;81(19):6538–47.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  118. Lorenzen W, Ahrendt T, Bozhuyuk KA, Bode HB. A multifunctional enzyme is involved in bacterial ether lipid biosynthesis. Nat Chem Biol. 2014;10:425–7.

    Article  CAS  PubMed  Google Scholar 

  119. Tomita T, Kuzuyama T, Nishiyama M. Structural basis for leucine-induced allosteric activation of glutamate dehydrogenase. J Biol Chem. 2011;286(43):37406–13.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  120. Zhang H, Mulholland GA, Seef S, Zhu S, Liu J, Mignot T, et al. Establishing rod shape from spherical, peptidoglycan-deficient bacterial spores. Proc Natl Acad Sci USA. 2020;117(25):14444–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We thank Dale Kaiser and Wenyuan Shi for providing M. xanthus strains. We thank Suchitha Ragunathan for contributions to library preparation and Sarah Wilson for generating differential expression heat maps in the early stages of data analysis.

Funding

This work was supported by the National Science Foundation (awards MCB-1411272 and IOS-1951025 to LK) and by salary support for LK from Michigan State University AgBioResearch.

Author information

Authors and Affiliations

Authors

Contributions

RR and LK designed the RNA-seq experiments. RR performed the RNA-seq experiments, analyzed the transcriptomes and differential expression, and contributed to writing the manuscript. MF performed the additional bioinformatic analyses, prepared figures, and was a major contributor in writing the manuscript. LK contributed to the design of bioinformatic analyses and figures, and was a major contributor in writing the manuscript.

Corresponding authors

Correspondence to Ramya Rajagopalan or Lee Kroos.

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.

Additional file 1

: Principal component analysis for transcript reads of M. xanthus wild-type (WT) strain DK1622 and the csgA, fruA, and mrpC mutants. A) Scree plot of the proportion of variance per principal component. B) Heatmap of Euclidean similarity scores within strains and between strains. A lower score denotes higher similarity. C) PCA plot for each replicate of the M. xanthus strains.

Additional file 2

: Expanded comparison of RNA-seq results with regulation reported previously. Includes the genes shown in Table 1 plus genes in the same operons and additional genes mentioned in the text

Additional file 3

: Strongly expressed genes that are differentially expressed in the mrpC mutant and have a ChIP-seq peak between -200 to +200 bp relative to a transcriptional start site (TSS). Genes shown as red dots in Figure 2b and 2c match the criteria in columns F and M of this table

Additional file 4

: Lists of genes with the indicated patterns of regulation in mutants based on RNA-seq results

Additional file 5

: List of all genes differentially expressed in one or more mutants. The table includes old and new gene locus identifiers, RNA-seq results, and annotation

Additional file 6

: Lists of top 10 most impacted up- or down-regulated genes in each mutant based on the absolute value of -log10 FDR-adjusted p-value*log2 fold difference

Additional file 7

: Spearman correlations of differentially downregulated (left) or upregulated (right) genes observed in only the fruA and mrpC mutants (top) or in the csgA, fruA, and mrpC mutants (bottom). A linear, smoothed trendline is overlaid for illustrative purposes

Additional file 8

: Expanded list of statistically enriched Gene Ontology (GO) groupings identified by both ShinyGO and DAVID analyses. Includes gene lists for GO groupings shown in Figure 5 plus closely related GO groupings that are not shown in Figure 5 because their gene lists overlap greatly with one that is shown, as indicated in column I

Additional file 9

: Cytoscape session file showing predicted STRING protein-protein interaction networks for proteins encoded by genes with each indicated pattern of regulation based on RNA-seq

Additional file 10

: KEGG pathways enriched among genes with the indicated patterns of regulation in mutants based on RNA-seq results

Additional file 11

: Comparison of the results for individual genes in RNA-seq based transcriptomic studies of fruA mutants

Additional file 12

: Spearman correlation analysis for genes observed to be differentially downregulated (left) or upregulated (right) in the fruA mutants of McLoon et al. (78) 12-h (top) and 24-h (bottom) datasets and this work. A linear, smoothed trendline is overlaid for illustrative purposes

Additional file 13

: Genes differentially expressed only in the mrpC mutant and with an MrpC ChIP-seq peak between -200 to +200 bp relative to a transcriptional start site (TSS).

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

A. Farrugia, M., Rajagopalan, R. & Kroos, L. Transcriptomic analysis of Myxococcus xanthus csgA, fruA, and mrpC mutants reveals extensive and diverse roles of key regulators in the multicellular developmental process. BMC Genomics 26, 355 (2025). https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12864-025-11417-z

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12864-025-11417-z

Keywords