text:
Transcriptome dynamics of the Myxococcus xanthus multicellular developmental program
José Muñoz-Dorado ,
Aurelio Moraleda-Muñoz,
Francisco Javier Marcos-Torres,
Francisco Javier Contreras-Moreno,
Ana Belen Martin-Cuadrado,
Jared M Schrader,
Penelope I Higgs,
Juana Pérez,
Universidad de Granada, Spain;
Universidad de Alicante, Spain;
Wayne State University, United States;
Tools and Resources Oct 14, 2019
The bacterium Myxococcus xanthus exhibits a complex multicellular life cycle. In the presence of nutrients, cells prey cooperatively. Upon starvation, they enter a developmental cycle wherein cells aggregate to produce macroscopic fruiting bodies filled with resistant myxospores. We used RNA-Seq technology to examine the transcriptome of the 96 hr developmental program. These data revealed that 1415 genes were sequentially expressed in 10 discrete modules, with expression peaking during aggregation, in the transition from aggregation to sporulation, or during sporulation. Analysis of genes expressed at each specific time point provided insights as to how starving cells obtain energy and precursors necessary for assembly of fruiting bodies and into developmental production of secondary metabolites. This study offers the first global view of developmental transcriptional profiles and provides important tools and resources for future studies.
Myxococcus xanthus is a soil-dwelling δ-proteobacterium that exhibits a complex multicellular life cycle with two phases: growth and starvation-induced development (Muñoz-Dorado et al., 2016). When nutrients are available, cells divide to produce a community known as swarm. Swarms are predatory (although not obligate) and can digest prokaryotic and eukaryotic microorganisms (Pérez et al., 2016). Upon starvation, cells in the swarm enter a developmental program, during which they migrate into aggregation centers and climb on top of each other to build macroscopic structures termed fruiting bodies. To form fruiting bodies, starving cells glide on solid surfaces by using two motility systems, known as A- (adventurous) and S- (social) motility, which allow individual cell movement or group movement that requires cell-cell contact, respectively (Mauriello et al., 2010; Nan et al., 2014; Islam and Mignot, 2015; Chang et al., 2016; Schumacher and Søgaard-Andersen, 2017). After completion of aggregation (24 hr post-starvation), cells differentiate into environmentally resistant myxospores, which are embedded in a complex extracellular matrix (Figure 1). Each fruiting body contains ≈105–106 myxospores. Interestingly, only ≈10% of the starving population become myxospores (O'Connor and Zusman, 1991a) as most cells (around 60%) undergo programmed cell death, most likely to provide the rest of the population enough nutrients to successfully build fruiting bodies (Wireman and Dworkin, 1977; Nariya and Inouye, 2008). The remaining cells differentiate into a persister-like state, termed peripheral rods (PR), which surround the fruiting bodies (O'Connor and Zusman, 1991a; O'Connor and Zusman, 1991b; O'Connor and Zusman, 1991c). While PRs are morphologically similar to vegetative cells, myxospores are coccoid and are surrounded by a thick coat mainly consisting of polysaccharides (Kottel et al., 1975; Müller et al., 2010; Müller et al., 2012; Holkenbrink et al., 2014). Myxospores can germinate when nutrients are available, and collective germination of myxospores from a fruiting body generates a small swarm that facilitates cooperative feeding.
Schematic of the M. xanthus developmental program.
The time line indicates aggregation and sporulation phases. M. xanthus cells (yellow rods) aggregate into mounds (arrows indicate gliding to aggregation centers) and then differentiate into resistant spores (gray circles) to produce mature fruiting bodies. Peripheral rods (gray rods) remain outside of the fruiting bodies as a distinct differentiated state. Cells undergoing lysis are depicted with dashed lines.
The developmental program is directed by sophisticated, but not completely defined, genetic regulatory networks, which are coupled to a series of intra- and extra-cellular cues (Kroos, 2017; Muñoz-Dorado et al., 2016). The first cue is starvation, which triggers accumulation of cyclic-di-GMP and, via the stringent response, guanosine penta- and tetraphosphate [(p)ppGpp] inside the cells. These global signals somehow activate four master cascade modules (Nla24, Mrp, FruA, and bacterial enhancer-binding proteins [bEBPs]), which interconnect to control the correct timing of gene expression (Kroos, 2017). Proper progression of development requires intercellular communication, wherein cells produce and transmit five sequential extracellular signals, named A, B, C, D, and E (Bretl and Kirby, 2016).
Although much knowledge has been generated in the last 40 years about the M. xanthus developmental cycle, especially with respect to signaling and gene regulatory networks, we are far from having an overall picture of all the events that occur during aggregation and sporulation. Several partial transcriptome analyses from developmental samples based on microarrays have been published, which were restricted to a few genes related to bEBPs (Jakobsen et al., 2004; Diodati et al., 2008), two-component systems (TCSs) (Shi et al., 2008), A-signaling genes (Konovalova et al., 2012), or lipid metabolism (Bhat et al., 2014). Here, we used RNA-Seq technology to measure global changes in transcript abundance at seven time points during M. xanthus development, which represents a substantial step forward compared to previous analyses. We found that at least 19.6% of M. xanthus genes (1415/7229) had statistically significant changes in transcript abundance during development. These data and analyses provide, for the first time, a comprehensive view of the transcriptional regulatory patterns that drive the multicellular developmental program of this myxobacterium, offering an essential scaffold for future investigations.
Transcriptome analysis of the developmental program by RNA-Seq
Global gene expression patterns were examined by RNA-Seq analysis of the wild-type M. xanthus strain DK1622 developed on nutrient limited CF agar plates. RNA was harvested from two independent biological replicates at 0, 6, 12, 24, 48, 72, and 96 hr of development, reverse transcribed to cDNA, and sequenced by Illumina methodology. On average, 54.72 million read pairs and a coverage of 591X was obtained. After removing the ribosomal sequences, the genome coverage varied from 5.52 to 14.18X (median of 10.49X), enough to provide an adequate coverage of the mRNA fraction. The two sample-replicates showed a high degree of concordance in gene expression (R2 correlation >0.98), with the exception of 24 hr samples (R2correlation = 0.80), which may be related to lack of synchrony between cells in the transition from aggregation to sporulation. The median of both values was utilized for further analysis (Table 1 and Table 1—source data 1).
Statistical analysis of the M. xanthus DK1622 transcriptome raw data.
Data for each of the replicas at 0, 6, 12, 24, 48, 72 and 96 hr of development are shown.
#Gb
#mapped reads
#rRNA-reads
#clean reads (Non-rRNA)
%rRNA rate
Coverage (x)
R2 correlation
WT_0_1 5.70 61906718 60713758 1192960 98.07 13.05 0.99
WT_0_2 5.62 61117544 59826984 1290560 97.89 14.12
WT_6_1 5.02 54468281 53756387 711894 98.69 7.79 1.00
WT_6_2 4.80 52436258 51759879 676379 98.71 7.40
WT_12_1 5.04 54054132 53271826 782306 98.55 8.56 0.99
WT_12_2 5.37 57646891 56798096 848795 98.53 9.29
WT_24_2 5.38 53962216 53018166 944050 98.25 10.33
WT_48_1 6.32 62796702 61500435 1296267 97.94 14.18 0.99
WT_72_2 5.20 51866662 50779793 1086869 97.90 11.89
Table 1—source data 1
Number of reads for each ORF of Myxococcus xanthus at 0, 6, 12, 24, 48, 72 and 96 hr of development.
Download elife-50374-table1-data1-v1.xlsx
RPKM values of the developmental time course and correlation scores.
RPKM values reported here were calculated from the total number of non-tRNA/rRNA containing reads. The old (MXAN_) locus tags, new gene identifiers (MXAN_RS), gene name and predicted functions or pathways in which they have been previously implicated are included. The number of missing data points and fold change were used as criteria for the developmental gene analysis. DG or reason that genes were not included in the DGs is indicated.
As a first data validation step, we determined the expression profiles of two different developmental genes using β-galactosidase transcriptional reporters [spiA::Tn5-lacZ (strain DK4322) and fmgE::Tn5-lacZ (strain DK4294)] (Kroos et al., 1986) from cells developed under the same conditions used in this study. Comparison of these β-galactosidase activities to the RNA-Seq data indicated the patterns were similar and in agreement with those obtained with other strategies (Figure 2). Moreover, the expression profiles of many genes that have been previously characterized from β-galactosidase reporter activity, qRT-PCR, or microarray analyses were compared with those obtained with our data. This analysis revealed a general agreement (Figure 2—source data 1).
Validation of the RNA-Seq transcription patterns for genes spiA (MXAN_RS20760) (A) and fmgE (MXAN_RS16790) (B).
β-galactosidase specific activity (SA) of the strains harboring lacZ fusions to the respective genes (blue lines) compared to RNA-Seq RPKM values (red lines) at each developmental time point (h). Error bars indicate standard deviations for β-galactosidase specific activity determination.
Comparison of expression profiles of several developmental genes described in the literature and included in the DGs with the RPKM profiles from this study.
References are in Figure 3—source data 1.
Gene expression profiles organize into 10 developmental groups (DGs)
To identify developmentally regulated transcripts with similar expression patterns, genes containing measured RPKM (reads per kilobase pair of transcript per million mapped reads) values for all time points were further analyzed. First, all genes with <50 reads and/or high replicate variability between the two replicate datasets (R2 correlation <0.7) were removed (Table 1—source data 2). 1557/7229 (21.5 %) genes passed quality criteria filters. Some of these genes (142; 9.1%) were not significantly (>2 fold) up- or down regulated during the developmental program, suggesting these genes could be considered constitutively expressed. 1415 significantly regulated genes (90.9% of the genes passing quality control) were then analyzed for clusters of similar expression patterns. Briefly, initial evaluation of several clustering methods (see Materials and methods for details) revealed kmeans clustering with 6–12 clusters produced the best clustering of genes with similar expression profiles. Refinement of kmeans clusters by visual inspection indicated 10 DGs best explained the number of unique gene profiles (Figure 3A, Figure 3—figure supplement 1, and Figure 3—source data 1). Clusters were organized with peak expression profiles corresponding to progression through the developmental program (Figure 3A). The relative expression profiles presented here in the heat maps are a log2 normalized RPKM value relative to the mean of the entire RNA-Seq trajectory for a given gene. This relative expression profile is log2(RPKMi time-point x/RPKMi average), where i is a given gene, x is a given time-point, and RPKM average is the average RPKM value of all time-points. Although the sensitivity of the RNA-Seq technology allows detection of genes with low expression levels, some genes were excluded in this analysis after removing those with less than 50 reads. This is especially important for genes that are expressed during growth at levels not detectable with this technology. Only 13 genes were found in this situation (Figure 3—source data 2), which have not been included in the analyses shown below.
The relative expression profiles of M. xanthus genes observed during the developmental program compared to those previously observed during chemical-induction of sporulation.
(A) Relative expression profiles of significantly regulated genes at the indicated hours after induction of starvation. Genes were clustered into 10 developmental groups based on the time of peak expression and then organized according to the temporal progression of development. Developmental group number and the phase of the developmental program (with photographs of aggregates under the dissecting microscope and cells under the scanning electron microscope) are indicated to the left of the heat map. In the photographs, red, blue and yellow bars represent 2 mm, 100 µm, and 5 µm, respectively. (B) Relative expression levels of the genes in panel A during the indicated hours after chemical-induction of sporulation (Müller et al., 2010). The position of individual genes in panel B is matched to panel A. Relative expression levels for panels A and B are indicated by color code according to the legend. DGs significantly represented in up-, down-, or not-regulated sporulation gene sets are indicated by the probability values in yellow, blue or gray, respectively (C) Comparison of up-, down-, or not-regulated starvation-induced development and chemical-induced sporulation gene tallies. Tally of glycerol-induced sporulation genes up- (top), down- (middle) or not-regulated (bottom) that are significantly enriched in the indicated DGs.
RPKM values for developmentally controlled genes and distribution of the genes in the ten developmental groups.
Genes with >2 fold change in mRNA levels and >0.7 R2 correlation between the two replicate developmental time courses are listed with their corresponding RPKM mRNA level measurements, following the same order shown in Figure 3A. The old (MXAN_) locus tags, new gene identifiers (MXAN_RS), gene name and predicted functions or pathways in which they have been previously implicated are included.
Genes with no reads at time 0 hr.
Since no data existed at 0 hr for these genes (and therefore not included in Figure 3A), we compared correlation scores of the developmental time points and found that 7 out of 13 genes had high correlations and passed the developmental regulation criteria (highlighted in blue). Four of these genes are most likely forming operons with other genes included in DG4 and DG7 based on the genomic environment and relative expression profiles (written in red). The profiles have been represented assuming that relative time at 0 hr = 0 RPKM. Genes excluded because of low correlation values are highlighted in yellow.
Previously reported developmental genes and identification in the 10 developmental groups shown in Figure 3A.
Genes have been highlighted by color (see color code) to depict known or implied roles in developmental processes and/or membership of protein families involved in development. References are the same as in Figure 3—source data 1.
Genes involved in development that are not included in the DGs shown in Figure 3A and reasons for their exclusion.
Genes included for comparison of developmental and sporulation transcriptomes (Figure 3B).
Developmental transcriptome genes (this study) include those assigned to DGs and those which passed requirements for all data points with >0.7 R2 correlation (considered constitutive). Sporulation transcriptome data were from Müller et al. (2010), in which vegetative cells were induced to artificially sporulate by addition of glycerol to 0.5 M. The sporulation gene data set included genes significantly down-, up-, or not-regulated. Up-regulated genes were classified as up1 (peak expression within 2 hr of sporulation induction) or up2 (peak expression at 2–4 hr of sporulation induction). Only genes present in both studies are included in this table.
Tally of enrichment of sporulation transcriptome genes (Müller et al., 2010) found in each of the DGs (Figure 3B).
Only genes identified as having reliable expression patterns in both developmental and sporulation transcriptome studies (Figure 3—source data 5) were included. Spore transcriptome genes enriched or depleted in any DG were tested for significance (p<0.01) using a chi square test [as per Müller et al. (2010) generating probability (p) values that this enrichment was due to chance. Sporulation transcriptome up-, down- or not-regulated gene categories significantly over-enriched in any DGs are highlighted in yellow, genes significantly depleted from any DG category are highlighted in blue. Text in red: sporulation up-regulated genes are enriched in DGs 7 and 8.
Tally of genes observed in both the developmental (this study) and sporulation [Müller et al., 2010 BMC Genomics 11:264] transcriptomes as supporting data for Figure 3C.
Total number indicates genes with expression profiles considered reliable in both data sets. Developmental DG1 was considered down-regulated genes, while up-regulated genes were considered all genes in DGs 2–10. Developmental constitutive genes were those whose expression level was not more than 2-fold different from vegetative cells at any time point.
DGs 1 and 2 contained genes that were immediately down-regulated relative to growth conditions, DGs 3–5 correlated with the aggregation phase, DGs 6–7 correlated with the transition from aggregation to sporulation, and DGs 8–10 correlated with the sporulation phase (Figure 3A). The number of genes attributed to each cluster ranged from 5.5% (DG3) to 14.1% (DG8) (Figure 3—figure supplement 2). Similar proportions of genes exhibited peak expressions in growth (21.8%), aggregation (23.8%), and transition from aggregation to sporulation (18.8%) phases. The final phase accounted for peak expressions of 35.6% of genes, consistent with the significant morphological and physiological rewiring that must occur during spore differentiation and preparation for extended quiescence.
When the number of genes by RPKM values that account for 50% of all mRNA expression at each time point was compared, it was observed that this number was low (~300) during growth, and at 48 hr it steadily increases to 762 by the end of sporulation (Figure 3—figure supplement 3). This suggests that the transcriptome is becoming distributed across a broader number of genes during development.
As a first step in analysis of the DGs, genes that have been previously described to affect M. xanthus development were identified. Of the ~280 characterized genes, 167 were included in the DGs. Most showed peak expression patterns at a time point that matched with the developmental phase where they have been reported to function, some of which have a well-defined role on development (Figure 3—source data 3 and Figure 3—figure supplement 4). A notable exception to this are genes that are thought to function in very early stages during development, but are included in DG10. Examples of these genes are sasS, which is involved in A signaling (Yang and Kaplan, 1997); csgA, which encodes the C signal (Shimkets and Rafiee, 1990); and romR, which is involved in polarity control of motility (Leonardy et al., 2007). However, DG10 genes also showed high expression at 6 hr, consistent with early activities and raising the possibility of additional functions at the latest stage of development. Finally, not all known developmental genes appeared in the 10 DGs defined here (Figure 3—source data 4). Of the 108 missing genes, 28 genes encoded TCSs and 20 serine/threonine protein kinases (STPKs).
Genes identified through chemical induction of sporulation mainly map to DG7 and DG8
A core sporulation transcriptome was previously defined using an artificial method for inducing spore differentiation (Müller et al., 2010). Using this method, myxospores can be induced in cells growing in rich broth culture by addition of chemicals such as 0.5 M glycerol (Dworkin and Gibson, 1964). Chemical-induced sporulation bypasses the requirement for starvation, motility (aggregation), and alternate cell fates (Higgs et al., 2014). Comparison of these two transcriptome sets determined that 1388 genes passed quality criteria in both transcriptome studies (Figure 3B, Figure 3—source datas 5 and 6, Figure 3C, and Figure 3—source data 7). 92% (218/237) of genes significantly up-regulated (>2 fold) in the sporulation transcriptome were also significantly up-regulated in the developmental transcriptome. These genes were significantly over-represented in DGs 7 and 8, with low probability (p) that this association is due to random chance (p=1×10−12 and p=2×10−11, respectively). DGs 7 and 8 peak expression at the transition to sporulation and sporulation phases, respectively. Co-regulated genes in these groups include those predicted to be involved in generation of sugar precursors for spore coat (see below). Likewise, genetic loci involved in spore coat synthesis (exo) and surface polysaccharide arrangement (nfs) (Müller et al., 2012; Holkenbrink et al., 2014) fell in DG8.
Interestingly, DGs 9 and 10, with peak expression profiles corresponding to the latest developmental time points, were not well represented in the up-regulated sporulation transcriptome. DG9 genes were significantly over-represented in the down-regulated sporulation transcriptome (p=0007). The DG9 cluster represents genes that were expressed during growth in rich media, down-regulated in response to starvation, and later up-regulated during the final phases of development (Figure 3A). Genes in this cluster appear to be involved in transport, respiration, and transcriptional regulation, and may be required during rapid growth and perhaps in the final maturation phases of sporulation. It is likely that the chemical-induced sporulation transcriptome may not have included very late induced sporulation genes as RNA was harvested until four-hours after induction. Chemical-induced spores are heat and sonication resistant at this stage, but final maturation may continue past this point. Finally, DG10 contained genes that were over-represented (p=7×10−10) in the not-regulated sporulation transcriptome gene set. We speculate this pool of genes may be present in PRs, a cell fate enriched late in the developmental program and not present during chemical induction of sporulation. Genes in this DG encode several STPKs, the bEBP Nla26, and several proteins involved in secondary metabolite (SM) biosynthesis (see below).
Good correlation was observed in genes that were expressed during growth and down-regulated in both transcriptome sets, with genes from DGs 1, 2, 5, and nine being over-represented in the down-regulated sporulation transcriptome set (Figure 3C). Of the large pool of genes that were not significantly regulated in the sporulation transcriptome (754), a relatively small number (12%) were also not significantly up- or down-regulated in the developmental transcriptome. This pool of genes likely represents constitutively expressed genes that serve as good normalization markers (Figure 3—source data 5) and includes housekeeping genes such as the transcription termination factor Rho (MXAN_RS11995), the cell cytoskeletal protein MreB (MXAN_RS32880) (Müller et al., 2012; Treuner-Lange et al., 2015; Fu et al., 2018), the gliding motility and sporulation protein AglU (MXAN_RS14565) (White and Hartzell, 2000), and the CheA homolog DifE (MXAN_RS32400), which is required for exopolysaccharide production and social motility (Yang et al., 2000).
Analysis of developmental regulated genes
All genes included in the 10 DGs were individually analyzed to find out which processes were affected during development, which may explain the different events that occur during aggregation and sporulation. Here, we have focused on six different processes.
A- and S-motility genes exhibit different developmental expression profiles
The two phases of the M. xanthus lifecycle (predatory growth and development) depend on A- and S-motility engines and their associated regulatory proteins (Pérez et al., 2014; Islam and Mignot, 2015; Mercier and Mignot, 2016; Schumacher and Søgaard-Andersen, 2017). Our data have revealed that many of the motility genes were developmentally regulated and that the expression profiles of the distinct A- and S-machinery genes clearly differed during development (Figure 4). A-motility genes were up-regulated during early development, while S-motility genes first decreased at 6 hr, and then returned to growth levels during aggregation (except for pilA), suggesting that A motility is preferentially used by cells during aggregation. During sporulation, expression of genes encoding both motility engines decreased (Figure 4), although some of them peaked again at this stage. The peak of some A-motility genes during sporulation is consistent with the repurposing of certain A-motility proteins to function in spore coat assembly (Wartel et al., 2013). At the end of development, the expression levels of the motility genes remain high. However, it should be reminded that although mature fruiting bodies are static structures, the PRs surrounding them are motile (O'Connor and Zusman, 1991b).
Developmental expression levels of M. xanthus motility proteins.
Schematic representation of the focal adhesion motor complexes necessary for A motility (top), the type IV pili motor complexes necessary for S motility (bottom), and the proteins involved in controlling polarity of both engines (polarity control module; center). The developmental expression levels (RPKM) of significantly regulated motility genes at the indicated times (h) of development are depicted. Gene expression profiles are colored to match the proteins shown in the schematic. Proteins depicted in gray represent genes that were not included in the developmental groups.
Gene expression patterns are consistent with use of glycogen and lipid bodies as energy sources
During growth, M. xanthus does not appear to consume sugars as carbon or energy sources (Watson and Dworkin, 1968; Bretscher and Kaiser, 1978). Instead, pyruvate, amino acids, and lipids are efficiently utilized which directly enter the tricarboxylic acid (TCA) cycle (Bretscher and Kaiser, 1978). It has been long debated as to whether M. xanthus utilizes a fully functioning glycolytic pathway (Curtis and Shimkets, 2008). It was speculated that the pathway may be utilized primarily in the gluconeogenic direction to produce sugar precursors necessary for spore coat production (Youderian et al., 1999; Chavira et al., 2007; Getsin et al., 2013). It is unknown how these pathways contribute to energy production during development when starving cells must synthesize energy currencies (i.e. ATP) over a period of at least three days.
Analysis of the DGs revealed that most genes involved in energy generation (pyruvate dehydrogenase complex, TCA cycle and oxidative phosphorylation) were found in DG2 (Figure 5A). Therefore, although they are down-regulated at 6 hr, they are later up-regulated at lower levels than during growth. In contrast, many genes encoding enzymes of the glycolytic/gluconeogenic pathway were up-regulated during development, with most reaching maximum expression levels at the completion of aggregation (24–48 hr) (Figure 5B). This up-regulated group included genes encoding homologs of glucokinase (glkC) and phosphofructokinase (pfkA), which are specific for the glycolytic pathway. These observations are consistent with a transcriptional rewiring of metabolic pathways during the developmental program, perhaps to take advantage of changing carbon/energy sources. For instance, developmental up-regulation of the glycolytic pathway genes may allow developing cells to obtain energy from sugars released from cells undergoing developmental lysis or from glycogen. Glycogen accumulates during late stationary phase/early development, and then disappears prior to sporulation (Nariya and Inouye, 2003). Enzymes predicted to be involved in synthesis of glycogen, such as GlgC, were found in DG4, while enzymes involved in utilization of glycogen, such as trehalose synthase, GlgP, and MalQ appeared to be constitutively expressed. The observation that these competing pathways show overlapping expression profiles suggests that regulation of glycogen production/consumption is likely regulated post-transcriptionally, as has been demonstrated by phosphorylation of PfkA by Pkn4 (Nariya and Inouye, 2003).
Relative developmental expression profiles of genes involved in energy generation.
(A) Genes encoding protein homologs for the pyruvate dehydrogenase complex, TCA cycle, and oxidative phosphorylation proteins. (B) Genes necessary for glycolysis/gluconeogenesis. Developmental time points in hours are indicated above each panel. Relative expression levels for panels A and B are indicated by color code according to the legend. For simplicity, the MXAN_RS designation was omitted from the locus tag of each gene.
Lipid bodies also accumulate in cells prior to sporulation. They mainly consist of triacylglycerides derived from membrane phospholipids as cells shorten in length, which are later used as an energy source (Hoiczyk et al., 2009; Bhat et al., 2014). However, the profiles of genes involved in both straight- and branched-chain primer synthesis and elongation of fatty acids observed in this study (Figure 6A and B) suggest that some level of lipid synthesis occurs during development. Moreover, genes involved in the alternative pathway to produce isovaleryl-CoA (Bode et al., 2009) were induced (Figure 6A and B). These lipids might be either incorporated into lipid bodies, be responsible for changes in lipid composition of the membranes of myxospores and/or PRs, and/or be used as precursors for SMs production. On the other hand, genes involved in lipid degradation reached maximum expression at 24–48 hr (β-oxidation) or even at later times (by other pathways) (Figure 6C). These fatty acid degradation enzyme profiles suggest that lipids are preferably consumed during sporulation.
Genes involved in synthesis and degradation of lipids.
(A) Simple representation of the M. xanthus branched fatty acid metabolic pathways depicting leucine degradation and alternative mevalonate-dependent routes. (B) Relative developmental expression profiles of the genes involved in straight-chain and branched-chain fatty acid biosynthesis as designated to the right. (B1) Straight-chain fatty acid primer synthesis; (B2) Branched-chain fatty acid primer synthesis of isovaleryl-CoA via leucine degradation (bkd genes); (B3) Branched-chain fatty acid primer synthesis of isovaleryl-CoA via the alternative pathway (mevalonate); B4: Fatty acid elongation. (C) Lipid degradation via β oxidation (C1) and other pathways (C2). Relative expression levels for panels B and C are indicated by color code according to the legend and developmental time points in hours are indicated above each panel. The MXAN_RS designation was omitted from the locus tag of each gene.
These data offer an overall picture of the central metabolism during development, which reinforces the notion that macromolecules recycled from growth phase or released from lysing cells can be directly used to yield energy, but are also used to synthesize glycogen and lipids that are stored for later consumption.
Amino acid and sugar precursors required for developmental macromolecule synthesis may be released by protein and polysaccharide turnover and gluconeogenesis
In addition to energy, starving cells need a source of sugar precursors to synthesize developmentally specific polysaccharides required for motility (Li et al., 2003), fruiting body encasement (Lux and Shi, 2005), spore coat synthesis (Kottel et al., 1975; Holkenbrink et al., 2014), and spore resistance (McBride and Zusman, 1989). It has been suggested that these sugars are derived from gluconeogenesis (Youderian et al., 1999). Consistently, our data have revealed that genes encoding enzymes specific for gluconeogenesis, such as phosphoenolpyruvate carboxykinase was in DG2, and GlpX (fructose-1,6-bisphosphatase) were present during growth and throughout development (Table 1—source data 2). Thus, these observations, as well as those presented above (Figure 5B), indicate that gluconeogenesis likely contributes to sugar precursor production at various stages during development. Moreover, the observation that four glycosyl hydrolases were specifically up-regulated during development (Figure 7A) suggests the cells may recycle vegetative polysaccharides or scavenge polysaccharides released from cells induced to lyse. These released free monomers could be synthesized into specific developmental polysaccharides by the series of glycosyl transferases that are also developmentally up-regulated (Figure 7B).
Developmental expression profiles of genes involved in production of polysaccharides and proteins.
Relative expression profiles of genes predicted to be necessary for polysaccharide hydrolysis (A), polysaccharide synthesis (B), and encoding proteases and peptidases (C). Developmental time points in hours are indicated above each panel and relative expression levels are indicated by color code according to the legend at the bottom. The MXAN_RS designation was omitted from the locus tag of each gene.
The developmental program is mainly triggered by amino acid starvation, yet developmentally specific proteins need to be newly synthesized. A source of these amino acids could be released from lysed cells and turnover of proteins that are not required during development. Consistently, 60 proteases appeared in the 10 DGs. 16 of these were specific to growth and down-regulated after 6 hr of starvation, while the rest were sequentially up-regulated at different time points (Figure 7C). All of these proteases likely release amino acids, but some of them may additionally function in regulatory processes, as has been reported for the protease PopC (DG7), which is thought to function in C-signal generation by cleavage of CsgA p25 to generate the p17 fragment (Lobedanz and Søgaard-Andersen, 2003; Rolbetzki et al., 2008).
Genes involved in secondary metabolism are developmentally up-regulated
M. xanthus produces multiple SMs, some of which (i.e. myxovirescine and myxoprincomide) facilitate predation (Xiao et al., 2011; Müller et al., 2016). However, we have found that although genes responsible for their biosynthesis are expressed during growth, their expression increases during development (Figure 8A–D). The M. xanthus genome codes for 18 nonribosomal peptide synthetases (NRPS), 22 polyketide synthases (PKS), and six mixed PKS/NRPS genes located in regions predicted to be involved in SMs synthesis (Korp et al., 2016). Of these, 10 were developmentally up-regulated, and 9/10 were found in DGs 8, 9 and 10 (Figure 8C–E). Moreover, 129 genes assigned to the DGs were located in those SM regions (Korp et al., 2016), and 51 of these were also identified in DGs 8, 9, and 10 (Figure 3—source data 1).
Developmental expression levels of genes involved in myxovirescine (antibiotic TA), myxoprincomide and other secondary metabolite biosynthesis.
(A) Schematic of the myxovirescine gene cluster. Names of genes depicted here are ta followed by the capital letter written with each arrow. (B) and (C) Developmental expression levels (RPKM) of significantly regulated myxovirescine genes plotted against the indicated developmental time points in hours. Gene expression profiles are colored to match the genes depicted in panel A. Genes depicted in gray were not included in the developmental groups. (D) Expression profiles of genes involved in myxoprincomide biosynthesis. (E) Profiles of others NRPS, PKS and PKS/NRPS not included in panels C and D. The MXAN_RS designation was omitted from the locus tag of each gene.
Together, these observations suggest that SMs may play previously unrecognized roles during development. For instance, they may be used to protect developing cells from other microbes in the soil, kill competitors to yield nutrients, or used as signaling molecules. In agreement with the data presented above concerning chemical-induced sporulation, an intriguing possibility is that PRs may produce SMs to defend spores inside the fruiting bodies or to release nutrients from prey to promote germination.
Translation may be developmentally rewired
It was previously reported that the protein composition of ribosome complexes purified from growing cells versus myxospores was different (Foster and Parish, 1973), but this intriguing observation has not been further pursued. Data from this study suggest a role for translation regulation in control of M. xanthus development. 76 genes involved in translation fit the criteria for inclusion in the DGs, which represents 5.4% of the genes included in them. Most of the genes involved in translation (i.e. encoding ribosomal proteins, initiation, elongation and termination factors, ribosome maturation and modification proteins, and several aminoacyl-tRNA ligases) are down-regulated by 6 hr, likely as the result of the stringent response (Starosta et al., 2014). Interestingly, most of them were subsequently up-regulated to a level similar to, or even higher than that observed during growth (Figure 9A and B). Some genes, including several aminoacyl-tRNA ligases, were mainly expressed during sporulation (Figure 9A and B).
Developmental expression profiles of genes involved in protein production.
(A) Relative expression levels of genes involved in translation or ribosome assembly. (B) Relative expression levels of genes encoding ribosomal proteins. Relative expression levels for panels A and B are indicated by color code according to the legend at the bottom, and developmental time points in hours are indicated above each panel. The MXAN_RS designation was omitted from the locus tag of each gene. (C) Developmental expression levels (RPKM) of the paralogous genes encoding protein S4 plotted against developmental time points in hours.
The observation that many ribosomal proteins are differentially regulated during development (Figure 9B) suggests that their relative ratio varies during development, yielding ribosomes with altered protein composition. It is remarkable the differences in expression profiles exhibited by duplicated genes of ribosomal proteins. M. xanthus encodes two paralogs for proteins S1, S4, S14, L28, and L33 (Table 1—source data 2). Only one of the two paralogous genes for ribosomal proteins S1 and S14 was found in the DGs (Figure 3—source data 1). Most notably, the two S4 paralogs were in the DGs, and exhibit similar RPKM values during growth. However, the RPKM for MXAN_RS16120 was 10-fold higher than that for MXAN_RS32955 at 48 hr of development (Figure 9C). Neither of the two paralogs of the 50S subunit were found in the DGs (Table 1—source data 2). These data confirm the previous results of Foster and Parish (1973), and provide a new overall perspective on the changing composition of translational machinery during development.
While it is possible that some of the transcriptional changes observed in the translational machinery could be related to ribosomal hibernation in myxospores and/or PRs (Yoshida and Wada, 2014; Harms et al., 2016; Gohara and Yap, 2018), we speculate that regulation of translational machinery may play an important role in directing the developmental program. It is known that some bacteria build alternative ribosomes to improve fitness on different growth conditions by altering the core ribosomal protein stoichiometry and differential expression of paralogous ribosomal proteins (Foster and Parish, 1973; Nanamiya et al., 2004; Hensley et al., 2012; Prisic et al., 2015; Slavov et al., 2015). In addition, some ribosomal proteins also play extraribosomal functions, such as control of transcription or mRNA decay (Warner and McIntosh, 2009).
A large interconnected regulatory network controls development
M. xanthus encodes a large repertoire of signaling/regulatory proteins presumed necessary to direct and coordinate its multicellular lifecycle in response to extra- and intra-cellular cues. Examples include one-component systems (OCS; transcriptional regulators that contain a sensing domain), TCS genes (sensor histidine kinase [HK] and response regulator [RR] proteins connected by phosphorelay), alternative sigma factors, and STPKs. Many of these proteins have been characterized in this bacterium, and in some cases, individual signal-transduction pathways and their exact role in controlling development have been defined (Inouye et al., 2008; Schramm et al., 2012; Muñoz-Dorado et al., 2014; Rajagopalan et al., 2014). However, the data presented here provide for the first time a view of the entire developmental cycle as an integrated system. We observed developmental up- or down-regulation of a larger number of regulatory genes than previously reported (Muñoz-Dorado et al., 2014; Rajagopalan et al., 2014). A significant number of these genes showed expression patterns that suggested they play important roles in modulation of aggregation and/or sporulation.
21 OCS genes were developmentally regulated, 9 of which were down-regulated. The remaining were differentially up-regulated during development (Figure 10A), including the repressor LexA (Campoy et al., 2003), and the regulators SasN (Xu et al., 1998) and MrpC (Sun and Shi, 2001; Ueki and Inouye, 2003; McLaughlin et al., 2018).
Developmental expression profiles of genes involved in transcriptional regulation and signal transduction.
Relative expression levels of genes encoding one-component regulators (A), two-component signal transduction proteins (B), sigma factors (C), and serine/threonine protein kinases (F) are depicted. Relative expression levels for panels A, B, C, and F are indicated by color code according to the legend at the bottom, and developmental time points in hours are indicated above each panel. The MXAN_RS designation was omitted from the locus tag of each gene. Expression levels (RPKM) of genes encoding the major sigma factors (D) and the subunits of the RNA polymerase (E) plotted against developmental time points in hours.
Of the 272 encoded TCS genes, we found 47 included in the DGs: 23 HKs, 10 hybrid HKs (HyHKs, containing HK and RR modules in the same polypeptide), and 14 RRs. 6 TCS genes were shut down during development (Figure 10B). Only 16 of 47 TCS genes have been previously characterized, thus pinpointing additional candidates for further characterization. Additionally, 5 of the 14 RRs belong to the group of bEBPs (plus MXAN_RS07605, which exhibits an architecture similar to bEBPs, but contains a GAF instead of a receiver domain). bEBPs function to activate expression from sigma 54 dependent promoters (Morett and Segovia, 1993). Besides the bEBPs, FruA is the only other RR found in the DGs that contains a DNA-binding domain. Out of the remaining 8 RRs, five were stand-alone receiver domains (CheY-like), two contain a putative diguanylate cyclase output domain, one of which has been characterized (Skotnicka et al., 2016), and one is RomR, which modulates motility (Leonardy et al., 2007). As many of the HKs and RRs that were developmentally regulated are orphans, the results presented here may help to identify cognate pairs.
Regarding other transcription factors, nine sigma factor genes were found to be developmentally regulated in a sequential fashion (Figure 10C), including those encoding the major sigma factor RpoD (Inouye, 1990) and RpoN (sigma 54) (Keseler and Kaiser, 1997). The expression profiles of these two sigma factor genes clearly differ during development. While rpoD was down-regulated at 6 hr, up-regulated during aggregation, and then down-regulated during sporulation, rpoN is up-regulated throughout development (Figure 10D). The expression profiles of rpoN and bEBP genes are consistent with previous results demonstrating that a bEBP cascade is triggered upon starvation (Giglio et al., 2011). Additionally, sigC, which encodes a group II sigma factor (Apelian and Inouye, 1993), was identified in DG7. The remaining six are predicted to encode extracytoplasmic functions (ECF) sigma factors, including rpoE1, involved in motility (Ward et al., 1998). It is noteworthy that the genes encoding the four subunits of the core RNA polymerase are developmentally regulated with profiles similar to that observed for the ribosomal proteins (Figure 10E).
M. xanthus encodes 99 STPKs, 11 of which have been reported to be pseudokinases, because they lack at least one of the three required catalytic residues (Muñoz-Dorado et al., 1991; Pérez et al., 2008). 22 STPK genes were included in the DGs (Figure 10F). Interestingly, seven encode pseudokinases while 15 encode predicted active kinases (Figure 3—source datas 1 and 3).
Together, these observations suggest that a high number of regulators directs the developmental program of M. xanthus, with some acting simultaneously and others sequentially to perfectly modulate the different events that occur through development. In other model bacteria such as Caulobacter crescentus, more than 19% of the genes have been found to be developmentally regulated (Laub et al., 2000), a similar percentage to that found in M. xanthus, 57% of which are under direct control of the five master regulators identified in this bacterium (Zhou et al., 2015). In Bacillus subtilis, out of 4100 genes, 520 were dependent on Spo0A, but not on σF, while 66 were dependent upon both regulatory proteins (Fawcett et al., 2000). And in the case of Streptomyces coelicolor, 1901 genes (24% of the ORFs) exhibit differences when the substrate mycelium differentiates to a multinucleated mycelium, with a large number of transcriptional regulators involved (Yagüe et al., 2013). Although the number of regulators included in the DGs is larger in M. xanthus than in other model bacteria, it remains to be elucidated how many of them directly modulate development.
Toward elucidation of a complete developmental gene-regulatory circuit
M. xanthus is an extraordinary bacterium with a large coding potential and a complex lifecycle. The developmental cycle consists of two consecutive events: aggregation and sporulation, during which cells segregate into three different cell fates. Moreover, this process lasts over three days and is triggered by starvation. During this extended period, starving cells must glide to build fruiting bodies and synthesize numerous macromolecules exclusive to fruiting bodies and myxospores. Consistent with this complexity, the transcriptomic analyses presented here revealed that 1415 genes are developmentally regulated with a high degree of confidence, exhibiting expression times that peak at either growth, aggregation and/or sporulation. Analysis of individual genes and the processes in which they participate has shed some light about how cells regulate the expression of motility genes to allow the cells to reach the aggregation centers, or how they rewire metabolism to both obtain energy and monomers to build new macromolecules or to synthesize SMs. Moreover, these data have also revealed that the translational and transcriptional machinery is deeply altered to modulate the different events of development, offering new insights that require to be experimentally pursued to determine the function of all these regulators.
Although the role of translation in the regulation of the developmental cycle has not yet been addressed, a large number of transcriptional regulators that function during development have been identified (Figure 3—figure supplement 4, and Figure 11—source data 1). The data presented here corroborate the results obtained by other myxobacteriologists. For instance, our analyses are in good agreement with an established model in which starvation activates four genetic regulatory networks: bEBP, Mrp, FruA, and Nla24 modules, with the latter being activated by cyclic-di-GMP (Kroos, 2017). However, the high number of transcriptional regulators and signaling proteins found in this study to be developmentally regulated at the mRNA level is much higher than expected. As shown in Figure 11, many uncharacterized regulators peak at either aggregation or sporulation (Figure 11—source data 1). Those that are expressed at the same time point may be interconnected to properly modulate specific events and guarantee the proper sequential expression of genes. The analysis of the expression levels of the different regulators has revealed that fruA and mrpC exhibit maximum RPKM values of 7046 and 4415, respectively, while none of the others reach 1000 (except for two sigma factors) (Figure 11—source data 1). This may explain why FruA and MrpC are considered master regulators of development, while many of the rest have not been identified as playing crucial roles in the lifecycle of M. xanthus. Undoubtedly, the developmental mRNA expression profiling presented here will act as a blueprint for the complete elucidation of the M. xanthus developmental regulatory program. Now that the changes in gene expression are measured, identifying the regulatory inputs of each promoter will be critical to understand the complete genetic circuitry controlling development.
New signaling proteins that are developmentally regulated in M. xanthus identified in this study.
Code used to distinguish among types of regulators is indicated in the upper part, where OCS indicates one-component systems; HK, histidine kinase: hyHK, hybrid histidine kinase; CheY, CheY-likeresponse regulator; RR, response regulator; ECF, ECF sigma factor; STPK, active Ser/Thr protein kinase; PseudoK, pseudokinase. Numbers inside each symbol indicate the number of each type of regulator that have not been previously identified as being developmentally regulated. Information about proteins depicted here is shown in Figure 11—source data 1.
Figure 11—source data 1
Regulatory elements included in the developmental groups.
Time point of maximum expression is highlighted in green when genes where previously identified to play a role in development; in yellow if these genes have been previously identified but their role in development have not been well established; and in blue if they have been first identified in this study. References for genes previously identified and shown with names here are found in Figure 3—source data 1.
Download elife-50374-fig11-data1-v1.xlsx
Reagent type (species)
or resource
Strain, strain background (Myxococcus xanthus) DK1622 Kaiser, 1979; Goldman et al., 2006 Genome: NC_008095.1 Wild-type strain used to obtain RNA
Strain, strain background (Myxococcus xanthus) DK4322 Kroos et al., 1986 MXAN_RS20760; MXAN_4276 spiA::Tn5-lacZ; Tn5 lac (Kmr) Ω4521
Strain, strain background (Myxococcus xanthus) DK4294 Kroos et al., 1986 MXAN_RS16790; MXAN_ 3464 fmgE::Tn5-lacZ; Tn5 lac (Kmr) Ω4406
Chemical compound, drug DNase I Sigma-Aldrich Cat No./ID: AMPD1
Chemical compound, drug SuperScript III Reverse Transcriptase Invitrogen Cat No./ID: 18080044
Chemical compound, drug Escherichia coli DNA polymerase New England Biolabs D1806
Chemical compound, drug E. coli RNAse H Invitrogen Cat No./ID: 18021071
Chemical compound, drug Proteinase K Ambion Cat No./ID: 25530–015
Chemical compound, drug E. coli DNA ligase New England Biolabs Cat No./ID: M0205L
Chemical compound, drug lysozyme Roche Cat No./ID: 10837059001
Commercial assay or kit RNeasy Midi Kit Qiagen Cat No./ID: 75142
Commercial assay or kit RNA Protect Bacteria Reagent Qiagen Cat No./ID: 76506
Software, algorithm BWA software Li and Durbin, 2009
Software, algorithm SAMtools Li et al., 2009
Software, algorithm Artemis v.16.0.0 Rutherford et al., 2000
Software, algorithm Cluster 3 Software Package de Hoon et al., 2004
Preparation of cells for RNA-Seq experiment
M. xanthus strain DK1622 (Kaiser, 1979; Goldman et al., 2006) was used in this study. Cells were grown in CTT liquid medium (Hodgkin and Kaiser, 1977) at 30°C with vigorous shaking (300 rpm) to 3.0 × 108 cells/ml (optical density at 600 nm [OD600] of 1), and then harvested and resuspended in TM buffer (10 mM Tris-HCl [pH 7.6]; 1 mM MgSO4) to a calculated density of 4.5 × 109 cells/ml (OD600 of 15). For each time replicate, 200 μl aliquots of concentrated cell suspension were spotted onto thirteen separate CF agar plates (Hagen et al., 1978). Two replicates of cells were harvested from plates at 6, 12, 24, 48, 72 and 96 hr (samples WT_6, WT_12, WT_24, WT_48, WT_72 and WT_96, respectively), and the obtained pellets were transferred immediately into 0.5 ml of RNA Protect Bacteria Reagent (Qiagen). Cells were then incubated at room temperature for 5 min, harvested by centrifugation at 5000 g for 10 min (4°C), and stored at −80°C after removal of the supernatant. For the t = 0 samples (sample WT_0), two replicates of 30 ml of the original liquid culture (OD600 of 1) were harvested by centrifugation as above, resuspended in RNA Protect Bacteria Reagent, and processed in the same manner.
To isolate RNA, frozen pellets were thawed and resuspended in 1 ml of 3 mg/ml lysozyme (Roche) and 0.4 mg/ml proteinase K (Ambion) prepared in TE buffer (10 mM Tris-HCl; 1 mM ethylenediaminetetraacetic acid [EDTA], pH 8.0) for cell lysis. Samples were incubated 10 min at room temperature. RNA extraction was carried out using the RNeasy Midi Kit (Qiagen) and each sample was eluted in 300 µl of RNase-free water. The concentration of RNA was measured using a NanoDrop ND-2000 spectrophotometer (NanoDrop Technologies, USA). To remove DNA, each RNA sample was supplemented with 1 unit of DNase I (from the DNA Amplification Grade Kit of Sigma) per µg of RNA and incubated at room temperature for 10 min. The reaction was stopped by adding the stop solution included in the kit and incubating 10 min at 70°C. The obtained RNA was precipitated with 1/10 vol of 3 M sodium acetate and 3 volumes of ethanol, and resuspended in 50 µl of RNase-free water. The quality of the total RNA was verified by agarose gel electrophoresis, and the concentration was determined using NanoDrop as indicated above.
Double stranded copy DNA synthesis
First strand DNA was synthesized using SuperScript III Reverse Transcriptase (Invitrogen) starting with 5 μg of RNA in a final reaction volume of 20 µl. In the next step, second-strand DNA was synthesized by adding 40 units of Escherichia coli DNA polymerase (New England Biolabs), 5 units of E. coli RNAse H (Invitrogen), 10 units of E. coli DNA ligase (New England Biolabs), 0.05 mM (final concentration) of dNTP mix, 10x second-strand buffer (New England Biolabs), and water to 150 µl. After 2 hr at 16°C, the reaction was stopped with 0.03 mM EDTA (final concentration). The obtained DNA was purified and concentrated using the DNA Clean and Concentrator Kit of Zymo Research according to manufacturer's instructions. The final product was eluted in DNA elution buffer from the kit to reach, at least, a yield of 2 µg of DNA, with a minimum concentration of 200 ng/µl.
Sequencing and transcriptomic data analysis
The cDNA from two biological replicates of each condition (see above) was used for sequencing using the Illumina HiSeq2000 (100 bp paired-end read) sequencing platform (GATC Biotech, Germany). Sequence reads were pre-processed to remove low-quality bases. Next, reads were mapped against M. xanthus DK1622 ribosomal operon sequences using BWA software with the default parameters (Li and Durbin, 2009). Remaining reads were subsequently mapped to the genome sequence with the default parameters and using the pair-end strategy. SAMtools (Li et al., 2009) was used to convert resulting data into BAM format. Artemis v.16.0.0 (Rutherford et al., 2000) was used for the visualization of the sequence reads against the M. xanthus genome. Once the transcripts were mapped to the genome, the average median value for each condition was used in further analyses (Table 1 and Table 1—source data 1 and 2).
Developmental gene analysis
Genes with fewer than 50 reads in a given time point were removed from analysis. RPKM values of the remaining genes were then compared across the developmental time-points. Developmental expression was characterized by a > 2 fold change in RPKM values across the time course and a > 0.7 R2 correlation coefficient of the two time course replicates. Log2 fold-change calculations were performed using the Cluster 3 Software Package (de Hoon et al., 2004). Genes passing these criteria are presented in Figure 3—source data 1. Randomization of the time point RPKM values within each replicate data set yielded a false discovery rate of 3.97% based on five randomized simulations that scrambled the order of the time points across all genes in the two datasets. Genes passing the developmental expression criteria were hierarchically clustered using pearson correlation, spearman rank, euclidian distance, and kmeans clustering. By visually inspecting the clusters we found that kmeans clustering gave the best clustering of genes with similar expression profiles.
Assay of β-galactosidase activity
For quantitative determination of β-galactosidase activity during development, strains containing lacZ fusions (Figure 2) were cultured and spotted onto CF plates as described above. Cell extracts were obtained at different times by sonication and assayed for activity as previously reported (Moraleda-Muñoz et al., 2003). The amount of protein in the supernatants was determined by using the Bio-Rad Protein Assay (Bio-Rad, Inc) with bovine serum albumin as a standard. Specific activity is expressed as nmol of o-nitrophenol (ONP) produced per min and mg of protein. The results are the average and associated standard deviation from three independent biological replicates.
To observe swarm and fruiting bodies, 10 µl of cell suspension prepared as mentioned above were spotted onto CTT (for growth) or CF (for development) agar plates and incubated at 30°C. Observation was on an Olympus SZX7 dissecting microscope. For scanning electron microscopy, samples obtained from CF and CTT agar plates were fixed with glutaraldehyde vapors for 24 hr at room temperature and then postfixed in aqueous 1% osmium tetroxide for 1 hr at 4°C, washed three times in buffer, and poststained for 2 hr in buffered 0.5% uranyl acetate. Dehydration was accomplished by a graded series of ethanol. Samples were then critical-point dried and sputter coated with gold. Photographs were taken in a Zeiss DSM950 scanning electron microscope.
The transcriptome sequencing data (raw-reads) was submitted to NCBI SRA under the Bioproject accession number: PRJNA493545. SRA accession numbers for each of the replicas are as follows: 0 hr: SAMN10135973 (WT_0_1-biological_replicate_1) and SAMN10135974 (WT_0_2-biological_replicate_2); 6 hr: SAMN10135975 (WT_6_1-biological_replicate_1) and SAMN10135976 (WT_6_2-biological_replicate_2); 12 hr: SAMN10135977 (WT_12_1-biological_replicate_1) and SAMN10135978 (WT_12_2-biological_replicate_2); 24 hr: SAMN10135979 (WT_24_1-biological_replicate_1) and SAMN10135980 (WT_24_2-biological_replicate_2); 48 hr: SAMN10135981 (WT_48_1-biological_replicate_1) and SAMN10135982 (WT_48_2-biological_replicate_2); 72 hr: SAMN10135983 (WT_72_1-biological_replicate_1) and SAMN10135984 (WT_72_2-biological_replicate_2); 96 hr: SAMN10135985 (WT_96_1-biological_replicate_1) and SAMN10135986 (WT_96_2-biological_replicate_2).
The transcriptome sequencing data (raw-reads) was submitted to NCBI SRA under the Bioproject accession number: PRJNA493545. SRA accession numbers for each of the replicas are as follows: 0 h: SAMN10135973 (WT_0_1-biological_replicate_1) and SAMN10135974 (WT_0_2-biological_replicate_2); 6 h: SAMN10135975 (WT_6_1-biological_replicate_1) and SAMN10135976 (WT_6_2-biological_replicate_2); 12 h: SAMN10135977 (WT_12_1-biological_replicate_1) and SAMN10135978 (WT_12_2-biological_replicate_2); 24 h: SAMN10135979 (WT_24_1-biological_replicate_1) and SAMN10135980 (WT_24_2-biological_replicate_2); 48 h: SAMN10135981 (WT_48_1-biological_replicate_1) and SAMN10135982 (WT_48_2-biological_replicate_2); 72 h: SAMN10135983 (WT_72_1-biological_replicate_1) and SAMN10135984 (WT_72_2-biological_replicate_2); 96 h: SAMN10135985 (WT_96_1-biological_replicate_1) and SAMN10135986 (WT_96_2-biological_replicate_2).
Martin-Cuadrado
Ana-Belen
(2018) NCBI BioProject
ID PRJNA493545. Myxococcus xanthus DK1622 transcriptome (TaxID: 34).
https://www.ncbi.nlm.nih.gov/bioproject/PRJNA493545
Apelian D
Inouye S
(1993) A new putative sigma factor of Myxococcus xanthus
Bhat S
Boynton TO
Pham D
Shimkets LJ
(2014) Fatty acids from membrane lipids become incorporated into lipid bodies during Myxococcus xanthus differentiation
Bode HB
Ring MW
Schwär G
Altmeyer MO
Kegler C
Jose IR
Singer M
Müller R
(2009) Identification of additional players in the alternative biosynthesis pathway to isovaleryl-CoA in the myxobacterium Myxococcus xanthus
ChemBioChem 10:128–140.
https://doi.org/10.1002/cbic.200800219
Bretl DJ
Kirby JR
(2016) Molecular mechanisms of signaling in Myxococcus xanthus Development
Journal of Molecular Biology 428:3805–3830.
Bretscher AP
Kaiser D
Nutrition of Myxococcus xanthus, a fruiting myxobacterium
Journal of Bacteriology 133:763–768.
Campoy S
Fontes M
Padmanabhan S
Cortés P
Llagostera M
Barbé J
(2003) LexA-independent DNA damage-mediated induction of gene expression in Myxococcus xanthus
Molecular Microbiology 49:769–781.
Chang YW
Rettberg LA
Treuner-Lange A
Iwasa J
Søgaard-Andersen L
Jensen GJ
(2016) Architecture of the type IVa pilus machine
Science 351:aad2001.
https://doi.org/10.1126/science.aad2001
Chavira M
Cao N
Le K
Riar T
Moradshahi N
McBride M
Lux R
Shi W
(2007) Beta-D-Allose inhibits fruiting body formation and sporulation in Myxococcus xanthus
https://doi.org/10.1128/JB.00792-06
Curtis PD
Myxobacteria: Multicellularity and Differentiation
Metabolic pathways relevant to predation, signaling, and development, Myxobacteria: Multicellularity and Differentiation, ASM Press.
de Hoon MJ
Imoto S
Nolan J
Miyano S
(2004) Open source clustering software
https://doi.org/10.1093/bioinformatics/bth078
Diodati ME
Gill RE
Plamann L
Initiation and early developmental events, Myxobacteria: Multicellularity and Differentiation, ASM Press.
Dworkin M
Gibson SM
(1964) A system for studying microbial morphogenesis: rapid formation of microcysts in Myxococcus xanthus
https://doi.org/10.1126/science.146.3641.243
Fawcett P
Eichenberger P
Losick R
Youngman P
(2000) The transcriptional profile of early to middle sporulation in Bacillus subtilis
Foster HA
Parish JH
(1973) Ribosomes, ribosomal subunits and ribosomal proteins from Myxococcus xanthus
Journal of General Microbiology 75:391–400.
https://doi.org/10.1099/00221287-75-2-391
Fu G
Bandaria JN
Le Gall AV
Fan X
Yildiz A
Mignot T
Zusman DR
Nan B
(2018) MotAB-like machinery drives the movement of MreB filaments during bacterial gliding motility
Getsin I
Nalbandian GH
Yee DC
Vastermark A
Paparoditis PC
Reddy VS
Saier MH
(2013) Comparative genomics of transport proteins in developmental Bacteria: Myxococcus xanthus and Streptomyces coelicolor
BMC Microbiology 13:279.
Giglio KM
Caberoy N
Suen G
Garza AG
(2011) A cascade of coregulating enhancer binding proteins initiates and propagates a multicellular developmental program
PNAS 108:E431–E439.
Gohara DW
Yap MF
(2018) Survival of the drowsiest: the hibernating 100S ribosome in bacterial stress management
Current Genetics 64:753–760.
Goldman BS
Nierman WC
Slater SC
Durkin AS
Eisen JA
Eisen J
Ronning CM
Barbazuk WB
Blanchard M
Field C
Halling C
Hinkle G
Iartchuk O
Kim HS
Mackenzie C
Madupu R
Miller N
Shvartsbeyn A
Sullivan SA
Vaudin M
Wiegand R
Kaplan HB
(2006) Evolution of sensory complexity recorded in a myxobacterial genome
Hagen DC
(1978) Synergism between morphogenetic mutants of Myxococcus xanthus
Developmental Biology 64:284–296.
Harms A
Maisonneuve E
Gerdes K
(2016) Mechanisms of bacterial persistence during stress and antibiotic exposure
Science 354:aaf4268.
https://doi.org/10.1126/science.aaf4268
Hensley MP
Gunasekera TS
Easton JA
Sigdel TK
Sugarbaker SA
Klingbeil L
Breece RM
Tierney DL
Crowder MW
(2012) Characterization of zn(II)-responsive ribosomal proteins YkgM and L31 in E. coli
Journal of Inorganic Biochemistry 111:164–172.
Higgs PI
Hartzell P
Holkenbrink C
Hoiczyk E
Myxobacteria: Genomics, Cellular and Molecular Biology
Myxococcus xanthus vegetative and developmental cell heterogeneity, Myxobacteria: Genomics, Cellular and Molecular Biology, Caister Academic Press.
Hodgkin J
(1977) Cell-to-cell stimulation of movement in nonmotile mutants of Myxococcus
McHugh CA
Bode E
Krug D
Lu JZ
(2009) Lipid body formation plays a central role in cell fate determination during developmental differentiation of Myxococcus xanthus
Kahnt J
(2014) Synthesis and assembly of a novel glycan layer in Myxococcus xanthus spores
The Journal of Biological Chemistry 289:32364–32378.
https://doi.org/10.1074/jbc.M114.595504
(1990) Cloning and DNA sequence of the gene coding for the major sigma factor from Myxococcus xanthus
Journal of Bacteriology 172:80–85.
https://doi.org/10.1128/jb.172.1.80-85.1990
Nariya H
Muñoz-Dorado J
Myxobacteria: Muticellularity and Differentiation
Protein Ser/Thr kinases and phosphatases in Myxococcus xanthus, Myxobacteria: Muticellularity and Differentiation, ASM Press.
Islam ST
(2015) The mysterious nature of bacterial surface (gliding) motility: a focal adhesion-based mechanism in Myxococcus xanthus
Seminars in Cell & Developmental Biology 46:143–154.
https://doi.org/10.1016/j.semcdb.2015.10.033
Jakobsen JS
Jelsbak L
Welch RD
Cummings C
Goldman B
Stark E
Slater S
(2004) Sigma54 enhancer binding proteins and Myxococcus xanthus fruiting body development
(1979) Social gliding is correlated with the presence of pili in Myxococcus xanthus
Keseler IM
(1997) sigma54, a vital protein for Myxococcus xanthus
Konovalova A
Wegener-Feldbrügge S
(2012) Two intercellular signals required for fruiting body formation in Myxococcus xanthus act sequentially but non-hierarchically
Molecular Microbiology 86:65–81.
Korp J
Vela Gurovic MS
Nett M
(2016) Antibiotics from predatory bacteria
Beilstein Journal of Organic Chemistry 12:594–607.
https://doi.org/10.3762/bjoc.12.58
Kottel RH
Bacon K
Clutter D
White D
Coats from Myxococcus xanthus: characterization and synthesis during myxospore differentiation
Kroos L
Kuspa A
(1986) A global analysis of developmentally regulated genes in Myxococcus xanthus
(2017) Highly Signal-Responsive gene regulatory network governing Myxococcus Development
Trends in Genetics 33:3–15.
https://doi.org/10.1016/j.tig.2016.10.006
Laub MT
McAdams HH
Feldblyum T
Fraser CM
Shapiro L
(2000) Global analysis of the genetic network controlling a bacterial cell cycle
Leonardy S
Freymark G
Hebener S
Ellehauge E
(2007) Coupling of protein localization and cell movements by a dynamically localized response regulator in Myxococcus xanthus
https://doi.org/10.1038/sj.emboj.7601877
Sun H
Ma X
Lu A
Zusman D
(2003) Extracellular polysaccharides mediate Pilus retraction during social motility of Myxococcus xanthus
Handsaker B
Wysoker A
Fennell T
Ruan J
Homer N
Marth G
Abecasis G
1000 Genome Project Data Processing Subgroup
(2009) The sequence alignment/Map format and SAMtools
Lobedanz S
(2003) Identification of the C-signal, a contact-dependent morphogen coordinating multiple developmental responses in Myxococcus xanthus
https://doi.org/10.1101/gad.274203
(2005) A novel bacterial signalling system with a combination of a ser/Thr kinase cascade and a his/Asp two-component system
Mauriello EM
(2010) Gliding motility revisited: how do the myxobacteria move without flagella?
McBride MJ
(1989) Trehalose accumulation in vegetative cells and spores of Myxococcus xanthus
McLaughlin PT
Bhardwaj V
Feeley BE
(2018) MrpC, a CRP/Fnr homolog, functions as a negative autoregulator during the Myxococcus xanthus multicellular developmental program
Molecular Microbiology 109:245–261.
https://doi.org/10.1111/mmi.13982
Mercier R
(2016) Regulations governing the multicellular lifestyle of Myxococcus xanthus
Current Opinion in Microbiology 34:104–110.
https://doi.org/10.1016/j.mib.2016.08.009
Moraleda-Muñoz A
Carrero-Lérida J
Pérez J
(2003) Role of two novel two-component regulatory systems in development and phosphatase expression in Myxococcus xanthus
https://doi.org/10.1128/JB.185.4.1376-1383.2003
Morett E
Segovia L
(1993) The sigma 54 bacterial enhancer-binding protein family: mechanism of action and phylogenetic relationship of their functional domains
Müller FD
Heider J
Huntley SM
(2010) Global transcriptome analysis of spore formation in Myxococcus xanthus reveals a locus necessary for cell differentiation
BMC Genomics 11:264.
Schink CW
Cserti E
(2012) Spore formation in Myxococcus xanthus is tied to cytoskeleton functions and polysaccharide spore coat deposition
Strack SN
Ryan SE
Shawgo M
Walling A
Harris S
Chambers C
Boddicker J
(2016) Identification of functions affecting Predator-Prey interactions between Myxococcus xanthus and Bacillus subtilis
Inouye M
(1991) A gene encoding a protein serine/threonine kinase is required for normal development of M. xanthus, a gram-negative bacterium
Cell 67:995–1006.
Elías-Arnanz M
Myxobacteria: Genomics Cellular and Molecular Biology
Abundance and complexity of signaling mechanisms in myxobacteria, Myxobacteria: Genomics Cellular and Molecular Biology, Caister Academic Press.
Marcos-Torres FJ
García-Bravo E
(2016) Myxobacteria: moving, killing, feeding, and surviving together
Frontiers in Microbiology 7:781.
https://doi.org/10.3389/fmicb.2016.00781
Oster G
(2014) Bacteria that glide with helical tracks
Current Biology 24:R169–R173.
Nanamiya H
Akanuma G
Natori Y
Murayama R
Kosono S
Kudo T
Kobayashi K
Ogasawara N
Park SM
Ochi K
Kawamura F
(2004) Zinc is a key factor in controlling alternation of two types of L31 protein in the Bacillus subtilis ribosome
(2003) An effective sporulation of Myxococcus xanthus requires glycogen consumption via Pkn4-activated 6-phosphofructokinase
(2008) MazF, an mRNA interferase, mediates programmed cell death during multicellular Myxococcus development
O'Connor KA
(1991a) Analysis of Myxococcus xanthus cell types by two-dimensional polyacrylamide gel electrophoresis
(1991b) Behavior of peripheral rods and their role in the life cycle of Myxococcus xanthus
(1991c) Development in Myxococcus xanthus involves differentiation into two cell types, peripheral rods and spores
Castañeda-García A
Jenke-Kodama H
(2008) Eukaryotic-like protein kinases in the prokaryotes and the myxobacterial kinome
Jiménez-Zurdo JI
Martínez-Abarca F
Millán V
(2014) Rhizobial galactoglucan determines the predatory pattern of Myxococcus xanthus and protects Sinorhizobium meliloti from predation
Environmental Microbiology 16:2341–2350.
(2016) Bacterial predation: 75 years and counting!
Environmental Microbiology 18:766–779.
Prisic S
Hwang H
Dow A
Barnaby O
Pan TS
Lonzanida JA
Chazin WJ
Steen H
Husson RN
(2015) Zinc regulates a switch between primary and alternative S18 ribosomal proteins in Mycobacterium tuberculosis
Rajagopalan R
Sarwar Z
Developmental gene regulation, Myxobacteria: Genomics Cellular and Molecular Biology, Caister Academic Press.
Rolbetzki A
Ammon M
Jakovljevic V
(2008) Regulated secretion of a protease activates intercellular signaling during fruiting body formation in M. xanthus
Developmental Cell 15:627–634.
Rutherford K
Parkhill J
Crook J
Horsnell T
Rice P
Rajandream MA
Barrell B
(2000) Artemis: sequence visualization and annotation
Bioinformatics 16:944–945.
https://doi.org/10.1093/bioinformatics/16.10.944
Schramm A
Lee B
(2012) Intra- and Interprotein Phosphorylation between Two-hybrid Histidine Kinases Controls Myxococcus xanthus Developmental Progression
Schumacher D
(2017) Regulation of cell polarity in motility and cell division in Myxococcus xanthus
Annual Review of Microbiology 71:61–78.
https://doi.org/10.1146/annurev-micro-102215-095415
Shi X
Huntley S
Hamann N
Hedderich R
(2008) Bioinformatics and experimental analysis of proteins of two-component systems in Myxococcus xanthus
Rafiee H
(1990) CsgA, an extracellular protein essential for Myxococcus xanthus development
Skotnicka D
Smaldone GT
Petters T
Trampari E
Liang J
Kaever V
Malone JG
(2016) A minimal threshold of c-di-GMP is essential for fruiting body formation and sporulation in Myxococcus xanthus
PLOS Genetics 12:e1006080.
https://doi.org/10.1371/journal.pgen.1006080
Slavov N
Semrau S
Airoldi E
Budnik B
van Oudenaarden A
(2015) Differential stoichiometry among core ribosomal proteins
Cell Reports 13:865–873.
Starosta AL
Lassak J
Jung K
Wilson DN
(2014) The bacterial translation stress response
FEMS Microbiology Reviews 38:1172–1201.
(2001) Analyses of mrp genes during Myxococcus xanthus development
Macia E
Guzzo M
Hot E
Faure LM
Jakobczak B
Espinosa L
Alcor D
Ducret A
Keilberg D
Castaing JP
Lacas Gervais S
Franco M
(2015) The small G-protein MglA connects to the MreB actin cytoskeleton at bacterial focal adhesions
The Journal of Cell Biology 210:243–256.
Ueki T
(2003) Identification of an activator protein required for the induction of fruA, a gene essential for fruiting body development in Myxococcus xanthus
Ward MJ
Lew H
Regulation of motility behavior in Myxococcus xanthus may require an extracytoplasmic-function sigma factor
Warner JR
McIntosh KB
(2009) How common are extraribosomal functions of ribosomal proteins?
Molecular Cell 34:3–11.
Wartel M
Thutupalli S
Bergam P
Brun YV
Shaevitz J
(2013) A versatile class of cell surface directional motors gives rise to gliding motility and sporulation in Myxococcus xanthus
PLOS Biology 11:e1001728.
Watson BF
Comparative intermediary metabolism of vegetative cells and microcysts of Myxococcus xanthus
Journal of Bacteriology 96:1465–1473.
Hartzell PL
(2000) AglU, a protein required for gliding motility and spore maturation of Myxococcus xanthus, is related to WD-repeat proteins
Wireman JW
Developmentally induced autolysis during fruiting body formation by Myxococcus xanthus
Xiao Y
Wei X
Ebright R
Wall D
(2011) Antibiotic production by Myxobacteria plays a role in predation
Xu D
Yang C
Myxococcus xanthus sasN encodes a regulator that prevents developmental gene expression during growth
Yagüe P
Rodríguez-García A
López-García MT
Martín JF
Rioseras B
Sánchez J
Manteca A
(2013) Transcriptomic analysis of Streptomyces coelicolor differentiation in solid sporulating cultures: first compartmentalized and second multinucleated mycelia have different and distinctive transcriptomes
Tong L
(2000) Myxococcus xanthus dif genes are required for biogenesis of cell surface fibrils essential for social gliding motility
(1997) Myxococcus xanthus sasS encodes a sensor histidine kinase required for early developmental gene expression
Yoshida H
Wada A
(2014) The 100S ribosome: ribosomal hibernation induced by stress
Wiley Interdisciplinary Reviews: RNA 5:723–732.
https://doi.org/10.1002/wrna.1242
Youderian P
Lawes MC
Creighton C
Cook JC
Mutations that confer resistance to 2-deoxyglucose reduce the specific activity of hexokinase from Myxococcus xanthus
Zhou B
Schrader JM
Kalogeraki VS
Abeliuk E
Dinh CB
Pham JQ
Cui ZZ
Dill DL
(2015) The global regulatory architecture of transcription during the Caulobacter cell cycle
Tâm Mignot
Reviewing Editor; CNRS-Aix Marseille University, France
Gisela Storz
Senior Editor; National Institute of Child Health and Human Development, United States
Lotte Sogaard-Andersen
Reviewer; Max Planck Institute for Terrestrial Microbiology, Germany
[Editors' note: a previous version of this study was rejected after peer review, but the authors submitted for reconsideration and the revised version was accepted for publication. The first decision letter after peer review is shown below.]
Thank you for submitting your work entitled "Global transcriptome analysis of the Myxococcus xanthus multicellular developmental program" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and a Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Lotte Sogaard-Andersen (Reviewer #1); Patrick Eichenberger (Reviewer #2); Roy Welch (Reviewer #3).
Our decision has been reached after consultation between the reviewers. Based on these discussions and the individual reviews below, we regret to inform you that your work cannot be considered for publication in eLife. All reviewers agree that the dataset is high quality and will constitute a valuable resource for the field. However, in the present manuscript, the current data are largely used to validate previous results in developmental gene regulation and as such, do not reveal "an integrated regulatory network". New analyses could take the manuscript in this direction and ways to perform them are suggested in each of the individual reviews. Given that the modification will likely take more than the two-months period that eLife grants for revisions, we can only offer to reconsider a deeply modified version that would take these comments into account. Please consider that if you decide to re-submit, the manuscript will have to go through another round of reviews before it can be considered for publication.
This manuscript describes a global analysis of gene expression changes during development in Myxococcus xanthus using RNAseq. During development, M. xanthus cells segregate into three cell types. The data presented here contains the combined changes for all three cell types and no attempts were made to distinguish between cell types. Previously, two transcriptome analyses from developmental samples based on microarrays were published (Diodati et al., Nla18, a key regulatory protein required for normal growth and development of Myxococcus xanthus. J Bacteriol. 188, 1733-1743 (2006) and Shi et al., 2008. Given the technological developments in the past decade to map transcriptomes, the data presented here represent a clear qualitative improvement and represent a rich and important novel resource for the community. Moreover, the data are carefully analyzed and provide a number a testable hypotheses for future studies.
1) Introduction third paragraph: Maybe briefly mention that Diodati et al. and Shi et al. previously published transcriptome analyzes of developing cells and how your data represent are different.
2) Subsection "Transcriptome analysis of the developmental program by RNA-Seq" paragraph two: I did not find microarray analyzes in the Bath et al., 2014 reference. Please clarify.
3) In the same paragraph: It is appreciated that the authors validate the quality of the RNAseq data using the two lacZ fusions. Later in the text, when expression profiles are discussed for genes that are important for development, it would be helpful to include a reference to those papers and mention if the listed genes were shown to be up- or down-regulated. Also, many more genes than the ones listed are important for development. To further validate the RNAseq data, I was wondering if it would be worthwhile to prepare a list with genes that are known to be important for development, their expression pattern from qRT-PCR or gene fusions and the expression profile reported here.
4) Subsection "Gene expression profiles organize into 10 developmental groups": The criteria for including genes in the analysis should be more clearly described: It is my understanding that in order for a gene to be included, it needs to have >50 reads at all seven time points. If this is correct, it seems that this would for instance exclude genes that are not expressed at T=0 or strongly downregulated genes. Please clarify.
5) Throughout the data presentations: It is not clear (to this reviewer) what is meant by "relative expression profiles". If the developmental samples are compared to the T=0 samples, then the T=0 induction ratio for the T=0 samples should be 0.00 but it is not. Please clarify precisely what is shown in the expression profiles.
6) Figure 3, 7, 8 and 9: Standard deviations cannot be derived from n<3; please remove SD from the figures.
7) Figure 8AB: The logic underlying the order in which genes are listed in these two figures is not clear. Maybe it would be a good idea to list them according to function, e.g. in 8B L subunits would be listed together and S-subunits would be listed together.
8) Subsection "A large interconnected regulatory network controls development" paragraph three: DmxB is the response regulator-diguanylate cyclase responsible for the increase in c-di-GMP during development and the dmxB gene was previously shown to be upregulated during development (Skotnicka et al., 2016). Is the dmxB gene among the upregulated response regulators?
This is a comprehensive study of the 96-hour program of fruiting body formation in the bacterium Myxococcus xanthus. Gene expression during the development program was investigated by collecting RNA at 7 successive time points and performing RNA-sequencing on two biological replicates. Based on clustering analyses, genes whose expression varied during development were assigned to 10 developmental groups. The composition of each group was discussed in view of 40 years of research on the model organism. This rich dataset will be invaluable to research groups working on Myxococcus and δ-proteobacteria. It might also hold some value to researchers working on other bacteria with well-characterized developmental programs (B. subtilis, Streptomyces, C. crescentus), but comparison to these systems has not been explored in the current version of this paper. In addition, the following points require clarification.
Subsection "Transcriptome analysis of the developmental program by RNA-Seq": "two independent biological replicates" and "(R2 correlation>0.98), with the exception of 24-h samples (R2 correlation=0.8)". Why only two and not three? Even though the correlation was high for most samples, the fact that there was lower correlation at one time point might have justified the addition of a third replicate. Also, "high replicate variability between the two replicate datasets (R2 correlation <0.7) were removed". A third replicate would help for that category as well.
"After removing the ribosomal sequences (about 98% of the reads)" This seems like a waste of resources. I wonder why the authors did not use a method to pull down rRNA before sequencing? If the authors had done so, they would have saved money for including a third biological replicate.
Second paragraph of subsection "Transcriptome analysis of the developmental program by RNA-Seq": "two independent biological replicates". The authors should insist on how the use of RNA-seq extends previous knowledge acquired from microarray experiments. The impact of the present paper is lessened if the data reported are just a confirmation of previous work.
The authors perform a global analysis of the Myxococcus xanthus developmental transcriptome, taking timepoints throughout the biofilm's starvation stress response past the self-organization of fruiting bodies and the differentiation of cells into spores (96 hours). They then compare these data to prior studies, both single gene and high throughput, to determine if the changes in RNA levels for different genes at different time points either match previous data or agree with contemporary functional interaction models. Their comparison serves both as data validation and as the primary source of their conclusions. Because the authors choose to largely limit their manuscript in this way, they miss several opportunities to perform a more varied and compelling set of analyses, and ultimately fails to deliver on their impact statement; the work, as presented, does not "reveal a genetic regulatory network"; with significant modification, it might.
Many interesting questions that could be addressed by these data are left unanswered. There is enough new knowledge to merit publication, but it will require additional analysis and substantial reorganization. Because the changes are significant, there is more than one way to make them, but at least three core problems need to be addressed in one way or another.
Problem 1: Analysis of the primary expression data set is not explained in sufficient detail for the reader to understand and reproduce it. Statistics were sometimes provided in the text with incomplete descriptions of analytical methods. For example, what specific methods were used to make the log2FC comparisons? How were criteria established for determining which genes fell into each DG? Were alternative methods tried before settling on these 10 DG clusters, and could alternative methods produce a different number of clusters and/or parsing of genes within them?
Proposed method to address Problem 1: Additional analyses should be performed on the primary data set to address the following questions: Would application of statistical methods, such as multiple range tests, confirm the stability of the 10 DGs, and are the differences in alternative methods statistically relevant? How much of the transcriptome is involved at each time point (i.e. what percentage of the genome is regulated, and how does that percentage change across the time points)? Are there genes that seem to be unique to each time point (i.e. are there any genes that are significantly up- or down-regulated only at one time point and may therefore be particularly important at that specific time point)? Are there other important expression patterns that can be revealed through different statistical methods, such as a functional PCA (i.e. what is the dimensionality of the time-course data)?
Problem 2: The main focus of the middle part of the manuscript (Figures 3-9) seems to be the confirmation of an existing 'consensus' interaction network, and the authors accomplish this by providing a supporting narrative. Although this narrative is interesting, the choice of genes to focus on and the alternative explanations for inconsistent data make this part of the manuscript subjective. Developmental data already exists for the majority of genes in the 'consensus' interaction network (Figure 1B); sometimes these data are from individual gene studies and sometimes they are from high-throughput studies (also, the authors should carefully qualify their claims of novelty and double check some references). This prior work can be used to generate hypotheses that can be tested using the their data set.
Proposed method to address Problem 2: The hypothesis could be something like "given the consensus functional interaction network and data from prior studies, we expect the expression patterns of this set of genes to follow the interaction network in the following way…". The hypothesis must include genes that do not meet the authors' expression cutoffs. Also, crucial conditional statements and plausible alternative explanations like the one stated in the final sentence of subsection "A‐and S‐motility genes exhibit different developmental expression profiles" can't be left to the end of a section, because they effectively negate everything that comes before it. These statements must be addressed throughout the presentation of relevant results and within the context of prior work, rather than as an afterthought.
Problem 3: Comparison of starvation-induced to glycerol-induced sporulation expression is a very interesting idea; an analysis might provide meaningful insight regarding the similarities and differences between these seemingly related events. The authors only begin to perform this analysis.
Proposed method to address Problem 3: The extent and nature of differences in gene expression between starvation-induced and glycerol-induced sporulation must be characterized and quantified. A superficial scan of Figure 2B reveals whole regions of different expression patterns. For example, what are all of the additional repressed genes in DG 9 during glycerol-induced sporulation? Could some of these genes provide insight into the differences between the two kinds of sporulation? Could some of these data be used to address the 'spore versus peripheral rod' alternative hypothesis proposed in subsection "A‐and S‐motility genes exhibit different developmental expression profiles"? Could some of these data be used to support the 'consensus' functional interaction network?
There are more minor issues involving the consistent use of abbreviations and nomenclature, claims of novelty, appropriate references, and the accuracy of concluding statements, but these can wait for the next round of revisions.
Major Points: Text
The text should be divided into three main sections: a comprehensive analysis of the gene expression data generated for this study (see point 1 above), a detailed comparison of these data to a 'consensus' functional interaction network (see point 2 above), and a detailed comparison of these data to the glycerol-induced sporulation time course (see point 3 above). All validation work (i.e. B-gal) should be put in supplementary materials, including figures that parse the expression data to support the narrative.
Major Points: Figures
The figures should be reorganized to better help the reader navigate the manuscript; at present, their structure and sequence do not match the text in a linear sequence. In the following discussion, the current figures are referred to as 'old Figure X' and proposed figures are referred to as 'new Figure X'.
New Figure 1A should be real images of M. xanthus development rather than a cartoon. It should include images that represent each of the time points taken so that the reader can see what development looks like under a microscope.
Old Figure 1B should be moved later in the manuscript because that is when it is discussed in the text.
New Figure 2 should have the cartoon of M. xanthus development running vertically on the left side of the heatmap accurately spaced to represent events, along with representative times. A timeline is also represented horizontally, but the authors use developmental stages and times interchangeably throughout the text, and so it would be helpful to have one diagram that runs along one axis showing clearly which times refer to which DGs and which developmental stages. The times should also match new Figure 1, so that the reader can easily go back and see what each stage actually looks like. Old Figure 2B should be moved later in the manuscript because the comparison to sporulation data is not yet discussed in the text.
New Figure 3 should include at least some set of analyses (see above).
Old Figures 3 through 9 should be moved to supplementary figures.
New Figure 4 should have a diagram showing the interaction of genes like one in Figure 1B. The Results section must be clear regarding which genes have transcription profiles that support the authors' hypotheses and which ones don't, and the figure should also include any genes that didn't make the expression cutoff to be included in the analysis. The interaction network diagram in old Figure 1B is designed to horizontally match to the cartoon showing development in old Figure 1A – this is a good idea for new Figure 4, but the matching should be more obvious, and should include developmental stages and DG groupings as in new Figure 2.
New Figure 5 should be the comparison between starvation induced development/sporulation vs glycerol-induced sporulation transcription profiles. Emphasis should be focused on genes whose expression profiles are similar between the two data sets, and genes whose expression profiles are different. Perhaps some of the genes from new Figure 4 could be identified as involved in development, sporulation, or both. There may also need to be a new Figure 6 providing additional statistical analysis, similar to new Figure 3.
[Editors' note: what now follows is the decision letter after the authors submitted for further consideration.]
Thank you for submitting your article "Transcriptome dynamics of the Myxococcus xanthus multicellular developmental program" for consideration by eLife. Your article has been reviewed by two peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Gisela Storz as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Lotte Sogaard-Andersen (Reviewer #2).
The authors have made considerable changes in the latest revision. Most importantly, the statistical analyses involved in the parsing and evaluation of expression data are now described clearly and early in the manuscript. Of course the authors' data are not perfect but, for these kinds of large data sets, variations between experimental replicates must be expected. For example, the relatively weak correlation between replicates of the 24 hour time point does not diminish the overall impact of the manuscript, even though the authors can only speculate about an explanation. The inclusion in supplementary materials of alternate heatmaps for different numbers of Developmental Groups (Figure 2A—figure supplement 1) is also very helpful to a reader who may well use these data for the authors' primary stated purpose – as "important tools and resources for future studies." The comparison between starvation and glycerol induced spores is particularly interesting.
[Editors' note: the author responses to the first round of peer review follow.]
Our decision has been reached after consultation between the reviewers. Based on these discussions and the individual reviews below, we regret to inform you that your work cannot be considered for publication in eLife. All reviewers agree that the dataset is high quality and will constitute a valuable resource for the field. However, in the present manuscript, the current data is largely used to validate previous results in developmental gene regulation and as such, it does not reveal "an integrated regulatory network". New analyses could take the manuscript in this direction and ways to perform them are suggested in each of the individual reviews. Given that the modification will likely take more than the two-months period that eLife grants for revisions, we can only offer to reconsider a deeply modified version that would take these comments into account. Please consider that if you decide to re-submit, the manuscript will have to go through another round of reviews before it can be considered for publication.
First, I would like to thank the three reviewers for their valuable comments. Undoubtedly, their suggestions have improved the manuscript. This kind of manuscript, with such a huge amount of data, are difficult to handle and focus, and I really think they have provided many inputs for improvement. Accordingly, we have done a great effort to address all their suggestions and concerns. As it can be observed in the new version, we have reorganized and rewritten most of the manuscript to make clear what it was known and what we consider are our most relevant contributions. Moreover, we have done new analyses, and accordingly, a considerable number of new figures and source data are included in the revised version. Furthermore, I would like to mention that, from our point of view, the most important part of the manuscript is to release all the raw data to the public. We consider that these data are very useful, and that they will be a great tool for many researchers. Therefore, I hope we have satisfactorily answered the reviewers' concerns. Alternatively, I hope we have satisfactorily justified why we have not attended a few demands.
These authors did not use microarrays, but they analyze data from other researchers to determine which genes involved in lipid metabolism are regulated during development. I think we have included this reference now in a more proper manner.
These authors did not use microarrays, please see response to previous point.
For simplicity, we do not mention in the text all genes that have been previously described to be important for development. All those genes, and their references, are included in the revised version in Figure 3A—source data 1 and 3. Moreover, we have added a Supplementary file 1, where we have compared the profiles obtained by microarrays, qRT-PCR and/or lacZ fusions of some studied genes (some of them shown in Figure 3A—figure supplement 4) with the profiles obtained in our study, to further validate our data.
I hope it is now clearer that we used three criteria to consider a gene as developmentally regulated: >50 reads, R2 correlation >0.7, and fold change >2 (subsection "Gene expression profiles organize into 10 developmental groups (DGs" and "Sequencing and transcriptomic data analysis"). We are aware that using filters will exclude genes that may be important for development. On the other hand, the use of these filters give a high degree of confidence for a gene to be considered as being developmentally regulated. And regarding the concern about exclusion of genes with no reads during growth, we agree with the reviewer that this deserve a special analysis. Therefore, we have analyzed genes in this situation. We have now mentioned these genes in the manuscript and a new table has also been included in the new version (Figure 3A—source data 2).
We have explained what "relative expression profiles" shown in the heat maps means (subsection "Developmental gene analysis").
We have eliminated the error bars when only two samples were used (see all figures).
We always ordered genes according to the developmental time-point when they reach maximum expression levels. Nevertheless, we have redone the panel for ribosomal proteins for two reasons: 1) differences between all genes are not very high, and 2) the reviewer is right to mention that this specific panel will be more understandable if we separate proteins of the two ribosomal subunits and we order them by their numbers (see new Figure 9).
Sorry for this error. We have included this gene in the text and in Figure 10.
We have included a brief comparison with the transcriptome of the three bacteria the reviewer mentions (subsection "A large interconnected regulatory network controls development").
We do not totally agree with the reviewer about the requirement of a third replicate, and we have several reasons that justify our disagreement. One reason is that not all the M. xanthus cells in the population are doing the same thing at the same time point (development is not synchronic, and there are particularly differences between the center and the edge of the spot), and for this reason we think that a third replicate would only introduce variability without adding much information. With a third replicate, many genes would be excluded because they will not accomplish the criteria we have used to consider a gene as being developmentally regulated, and they could be quite interesting for further studies, because they fit these restrictive criteria in two replicates. Actually, a few hours delay in one of the samples would make this study nearly useless. Actually, we think the R2 correlation is lower at 24 h than at other time points because some cells already started the sporulation process while other did not complete aggregation, yet (in the other time points they are either aggregating or sporulating). We have added a sentence in the manuscript to mention why variability might be higher at this specific time point. In fact, I must say that we decided to sequence two replicates instead of three for these reasons. We think this manuscript has generated new data (as those related to motility, transcription and translation) that will undoubtedly contribute to this field, because it is an excellent scaffold for future studies. Moreover, although a third replicate would offer a better statistical support, two biological replicates have been widely accepted by many authors and it is operationally fine. A correlation value of 0.8 at 24 h is still quite acceptable and does not justify the necessity of a third replicate. Finally, I would like to mention that I contacted with the journal about the need of a third replicate and they agreed that this is not necessary.
There is a plethora of kits based on different methodologies to remove the rRNA from the RNA-Seq library constructions. However, it is still out on the best method for ribosomal removal and several publications have shown extensive side-by-side tests of one method versus another and have declared different winners. An additional challenge is presented by the high GC presented by M. xanthus genome (68.89%), in which rRNA removal is less efficient. Another factor to consider is that, at the time of sequencing, since probe based rRNA removal methods are sequence specific, the availability of specific probes for Myxococcus were not available. Summarizing, at the end of the rRNA removal process, it is frequently found high amounts of rRNA (50% or more) with the putative mRNA lost in the laboratory manipulations. As this work present a time-series with two replicates, rationale indicates that the less manipulation of the samples, the better. Therefore, we chose the option to keep the rRNA. On the analysis side, ribosomal DNA content was filtered out so that differences in rRNA reads across samples do not affect alignment rates and skew subsequent normalization of the data.
We have rewritten many parts of the manuscript to clarify how our data contribute to the previous knowledge of M. xanthus development. Moreover, we have written a short comparison between our transcriptome work (global) and (partial) published microarrays transcriptomes (Introduction section).
We have rewritten many parts of the manuscript to emphasize on all the new knowledge our data offer. We have also softened our claim of finding a genetic regulatory network.
First, I would like to mention that classification of the genes in developmental group is just one of the many results we present in this manuscript. We consider that this classification was especially important to facilitate the analyses, to visualize how the transcriptome is changing upon starvation, and to validate our data with data already published, including the transcriptome of chemically induced myxospores. Therefore, definition of the 10 DGs is one of our data, but not the most relevant one. Once the groups are defined, the current paper focuses on many single gene expression profiles of developmental genes, including those that were expected to be developmentally regulated according to work already published in the last 40 years, and many of the new ones identified here, which have not been demonstrated to be important for development. Our work provides a lot of information that undoubtedly will have to be experimentally pursued to know the exact role of each gene in development. It seems that we failed in the narrative of this last part, and we have written the new version of the manuscript to show the analyses in a clearer manner to make readers understand all the novelties our data and analyses provide. Nevertheless, and going back to the definition of the 10 clusters, we have explained in more details how the 10 DGs were obtained, and why we chose 10 instead of other number (subsections "Gene expression profiles organize into 10 developmental groups (DGs)"; "Sequencing and transcriptomic data analysis"; and "Developmental gene analysis"). We have also included a new figure (Figure 3A—figure supplement 1) to show that the election of 10 groups instead of other number was the most appropriate one.
Proposed method to address Problem 1: Additional analyses should be performed on the primary data set to address the following questions: Would application of statistical methods, such as multiple range tests, confirm the stability of the 10 DGs, and are the differences in alternative methods statistically relevant?
As mentioned above, we have explained in more details how the 10 DGs were obtained, and why we chose 10 instead of other number. We have also included a new figure to show that the election of 10 groups instead of other number was the most appropriate one.
How much of the transcriptome is involved at each time point (i.e. what percentage of the genome is regulated, and how does that percentage change across the time points)?
We have included two new figure supplements to respond this question (Figure 3A—figure supplement 2 and 3).
Are there genes that seem to be unique to each time point (i.e. are there any genes that are significantly up- or down-regulated only at one time point and may therefore be particularly important at that specific time point)?
I think the response to this question may already be in the manuscript. For instance, all genes included in DG1 are only required during growth. Additionally, all genes in DG3 peak at 6 h. In other groups and time points, we can find some genes that peak at a specific time point. Nevertheless, we have mentioned in the new version of the manuscript some genes that are involved in regulation of gene expression that peak at specific time points or specific process (aggregation, sporulation or transition from one to the other) (Figure 11).
Are there other important expression patterns that can be revealed through different statistical methods, such as a functional PCA (i.e. what is the dimensionality of the time-course data)?
Although it could be interesting, according to all the responses given above we think this is out of the scope of this manuscript. I agree with the reviewer that many things can be done with the data we have obtained, and in fact, as we mention in the manuscript, our data represent a tool and resource to work in many directions in the future. I think this suggestion (and some others mentioned) is out of this scope of the manuscript we try to publish now. Even more, I encourage Roy, as well as other myxobacteriologists, to use our data once they are public to prepare other(s) manuscript(s) that will undoubtedly help to understand the M. xanthus developmental cycle.
We only partially agree with this this comment. In every section we have dealt with in the manuscript, we have tried to combine what it was already known with the complexity that our data provide. It is clear that we failed to clearly differentiate these two types of data, known and novel, and our novelty seems to be hidden in the manner we presented the data. For this reason, we have rewritten all this sections to make clear what it was already known, what we have found, and how everything matches to explain M. xanthus developmental cycle. I would also like to mention that the previous figure addressed a streamlined functional network (and new Figure 3A—figure supplement 4). This functional network diagram illustrates those genes for which known molecular mechanisms exist.
I hope the new presentation of the data satisfy this reviewer.
We think this could also be interesting, but analysis here is simply a hypothesis generator, and without experimental validation, there is no direct point to it. It could also be an interesting point to pursue in the future, but not in this manuscript. We have nevertheless addressed the differences between the two transcriptome sets more precisely and pointed out which gene sets have the potential to reveal differences in spore maturation and alternate cell types (we have rewritten the comparison in the revised version).
We disagree with this comment. As mentioned above, these kinds of data can be handled and many different ways, and while we are sure it would generate an interesting paper doing what the reviewer proposes, we believe that the data should be presented in the manner we have opted to do it.
Reorganization of the text and Figures was undertaken to address this issue as is indicated in the text.
We have decided to maintain the cartoon of Figure 1A (Figure A in the revised version). Maybe myxobacteriologists will understand microscopic figures, but a general reader can obtain a clearer information in a cartoon where many details can be drawn that do not appear in an image (for instance, cells that lyse cannot appear in a microscopic image). Nevertheless, and as microscopic images may also be helpful, we have included them in Figure 3A. As addressed above, Figure 1B (functional network) seems to have been problematic for not including ALL genes known to be involved in development. We have deleted previous Figure 1B and recreated a new figure, which is now placed later in text as Figure 3A—figure supplement 4.
This point includes several questions. Regarding the presentation of the situation of the developmental program at each time point, we think that Figure 3A (when it is clear what happens during development depicted in Figure 1)is the moment to include pictures taken under dissecting and electron microscopes. Therefore, we have included pictures about the aggregation stage of the cells and their morphology in the new Figure 3A. Regarding the proposal to move Figure 3B to later, we feel the figure would lose a lot of the information it contains in both panels. Each gene is represented at the same position in panels A and B. Therefore, it is possible to visualize what happens with each individual gene (and group) in both situations, development by starvation and spore formation by glycerol. These two panels have to be shown side by side to provide all the information we want to show. Therefore, we have maintained panels A and B in Figure 3. Moreover, as the section about comparison between both types of spores has been moved earlier in the text, no other figure is mentioned before Figure 3B.
New Figure 3 should include at least some set of analyses mentioned (see above).
We have included a new Figure (Figure 3A—figure supplement 1) to explain why we defined 10 DGs, as mentioned above.
We respectfully disagree again in this point. We consider all the information we show in these figures to be relevant and therefore insist they cannot be supplemental. We have a huge amount of information and it can be handled in many different manners. We consider the one we have chosen is the best one. Therefore, we have maintained Figures 4-10, with some modifications in some of them.
We have redone Figure 1, eliminating panel B from this figure. We agree with the reviewer that this position could be premature. We have moved this panel to supplementary data as Figure 3A—figure supplement 4, to show this information closer to a new and more complex regulatory network unveiled by our data.
Our feeling is that the reviewer overlooked the information we included as supplemental with the pie charts for this part (Figure 3B—source data 1, 2 and 3), since much of what he is asking for is already there, including the stat analysis. Nevertheless, we have rearranged the section and included a new panel as Figure 3C. We recommend to examining new Figure 3B—source data 1 and 2 and Figure 3C—source data 1.
José Muñoz-Dorado
Departamento de Microbiología, Facultad de Ciencias, Universidad de Granada, Granada, Spain
Substantial contributions to conception and design, acquisition of data, and analysis and interpretation of data; Drafting the article and revising it critically for important intellectual content; Final approval of the version to be published
jdorado@ugr.es
Aurelio Moraleda-Muñoz
Acquisition of data; Revising the article critically for important intellectual content; Final approval of the version to be published
Francisco Javier Marcos-Torres
Francisco Javier Contreras-Moreno
Ana Belen Martin-Cuadrado
Departamento de Fisiología, Genética y Microbiología, Universidad de Alicante, Alicante, Spain
Acquisition of data and analysis and interpretation of data; Revising the article critically for important intellectual content; Final approval of the version to be published
Jared M Schrader
Department of Biological Sciences, Wayne State University, Detroit, United States
Penelope I Higgs
Analysis and interpretation of data; Drafting the article and revising it critically for important intellectual content; Final approval of the version to be published
Juana Pérez
Spanish Government (CSD2009-00006)
Jose Munoz-Dorado
Spanish Government (BFU2016-75425-P)
National Institutes of Health (R35GM124733)
National Science Foundation (IOS 1651921)
Ministerio de Educación, Cultura y Deporte (Salvador de Madariaga Program)
This work has been supported by the Spanish Government (grants CSD2009-00006 to JMD and BFU2016-75425-P to AMM [70% funded by FEDER]), and by NIGMS of the National Institutes of Health under award number R35GM124733 to JMS. JPT and JMD were granted with fellowships of the Salvador de Madariaga Program to stay at Wayne State University for four months. PIH was funded by a grant from NSF IOS 1651921. We also want to thank Prof. Lee Kroos for providing strains.
Gisela Storz, National Institute of Child Health and Human Development, United States
Tâm Mignot, CNRS-Aix Marseille University, France
Lotte Sogaard-Andersen, Max Planck Institute for Terrestrial Microbiology, Germany
Received: July 20, 2019
Accepted: October 4, 2019
Version of Record published: October 14, 2019 (version 1)
© 2019, Muñoz-Dorado et al.
Myxococcus xanthus
transcriptome