Transcriptome dynamics of the Myxococcus xanthus multicellular developmental program

*For correspondence: [email protected] Competing interests: The authors declare that no competing interests exist. Funding: See page 22 Received: 20 July 2019 Accepted: 04 October 2019 Published: 14 October 2019 Reviewing editor: Taˆm Mignot, CNRS-Aix Marseille University, France Copyright Mun ˜oz-Dorado et al. This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited. Transcriptome dynamics of the Myxococcus xanthus multicellular developmental program Jose ´Mun ˜oz-Dorado 1 *, Aurelio Moraleda-Mun ˜oz 1 , Francisco Javier Marcos-Torres 1 , Francisco Javier Contreras-Moreno 1 , Ana Belen Martin-Cuadrado 2 , Jared M Schrader 3 , Penelope I Higgs 3 , Juana Pe ´rez 1 1 Departamento de Microbiologı ´a, Facultad de Ciencias, Universidad de Granada, Granada, Spain; 2 Departamento de Fisiologı ´a, Gene´tica y Microbiologı ´a, Universidad de Alicante, Alicante, Spain; 3 Department of Biological Sciences, Wayne State University, Detroit, United States Abstract 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. DOI: https://doi.org/10.7554/eLife.50374.001 Introduction Myxococcus xanthus is a soil-dwelling d-proteobacterium that exhibits a complex multicellular life cycle with two phases: growth and starvation-induced development (Mun ˜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 (Pe ´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øgaardAndersen, 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 »10 5 –10 6 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 Mun ˜oz-Dorado et al. eLife 2019;8:e50374. DOI: https://doi.org/10.7554/eLife.50374 1 of 27 TOOLS AND RESOURCES 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; Mu ¨ller et al., 2010;Mu ¨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. The developmental program is directed by sophisticated, but not completely defined, genetic regulatory networks, which are coupled to a series of intraand extra-cellular cues (Kroos, 2017; Mun ˜oz-Dorado et al., 2016). The first cue is starvation, which triggers accumulation of cyclic-diGMP and, via the stringent response, guanosine pentaand 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. Results and discussion 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 (R 2 Figure 1. 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. DOI: https://doi.org/10.7554/eLife.50374.002 Mun ˜oz-Dorado et al. eLife 2019;8:e50374. DOI: https://doi.org/10.7554/eLife.50374 2 of 27 Tools and resources Microbiology and Infectious Disease correlation >0.98), with the exception of 24 hr samples (R 2 correlation = 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). As a first data validation step, we determined the expression profiles of two different developmental genes using b-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 b-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 b-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). 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 (R 2 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) upor 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 Table 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. Sample name #Gb #mapped reads #rRNA-reads #clean reads (Non-rRNA) %rRNA rate Coverage (x) R 2 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_1 3.09 33003343 32498394 504949 98.47 5.52 0.80 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_48_2 5.14 34693098 34020090 673008 98.06 7.36 WT_72_1 6.43 63625946 62530758 1095188 98.28 11.98 0.99 WT_72_2 5.20 51866662 50779793 1086869 97.90 11.89 WT_96_1 6.33 63255801 62039033 1216768 98.08 13.31 0.99 WT_96_2 6.11 61187496 60080963 1106533 98.19 12.11 DOI: https://doi.org/10.7554/eLife.50374.003 The following source data is available for 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. DOI: https://doi.org/10.7554/eLife.50374.004 Source data 2. 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. DOI: https://doi.org/10.7554/eLife.50374.005 Mun ˜oz-Dorado et al. eLife 2019;8:e50374. DOI: https://doi.org/10.7554/eLife.50374 3 of 27 Tools and resources Microbiology and Infectious Disease 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 log 2 normalized RPKM value relative to the mean of the entire RNA-Seq trajectory for a given gene. This relative expression profile is log 2 (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. 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. Figure 2. Validation of the RNA-Seq transcription patterns for genes spiA (MXAN_RS20760) (A) and fmgE (MXAN_RS16790) (B). 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 b-galactosidase specific activity determination. DOI: https://doi.org/10.7554/eLife.50374.006 The following source data is available for figure 2: Source data 1. Comparison of expression profiles of several developmental genes described in the literature and included in the DGs with the RPKM profiles from this study. DOI: https://doi.org/10.7554/eLife.50374.007 Mun ˜oz-Dorado et al. eLife 2019;8:e50374. DOI: https://doi.org/10.7554/eLife.50374 4 of 27 Tools and resources Microbiology and Infectious Disease Figure 3. 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 mm, and 5 mm, respectively. (B) Relative expression levels of the genes in panel A during the indicated hours after chemical-induction of sporulation (Mu ¨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. DOI: https://doi.org/10.7554/eLife.50374.008 The following source data and figure supplements are available for figure 3: Source data 1. RPKM values for developmentally controlled genes and distribution of the genes in the ten developmental groups. DOI: https://doi.org/10.7554/eLife.50374.013 Source data 2. Genes with no reads at time 0 hr. DOI: https://doi.org/10.7554/eLife.50374.014 Source data 3. Previously reported developmental genes and identification in the 10 developmental groups shown in Figure 3A. DOI: https://doi.org/10.7554/eLife.50374.015 Source data 4. Genes involved in development that are not included in the DGs shown in Figure 3A and reasons for their exclusion. DOI: https://doi.org/10.7554/eLife.50374.016 Source data 5. Genes included for comparison of developmental and sporulation transcriptomes (Figure 3B). DOI: https://doi.org/10.7554/eLife.50374.017 Source data 6. Tally of enrichment of sporulation transcriptome genes (Mu ¨ller et al., 2010) found in each of the DGs (Figure 3B). DOI: https://doi.org/10.7554/eLife.50374.018 Source data 7. Tally of genes observed in both the developmental (this study) and sporulation [Mu ¨ller et al., 2010 BMC Genomics 11:264] transcriptomes as supporting data for Figure 3C. DOI: https://doi.org/10.7554/eLife.50374.019 Figure 3 continued on next page Mun ˜oz-Dorado et al. eLife 2019;8:e50374. DOI: https://doi.org/10.7554/eLife.50374 5 of 27 Tools and resources Microbiology and Infectious Disease 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 (Mu ¨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) (Mu ¨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 Figure 3 continued Figure supplement 1. Relative expression profiles of significantly regulated genes at the indicated hours after induction of starvation. DOI: https://doi.org/10.7554/eLife.50374.009 Figure supplement 2. Tally of reliable genes. DOI: https://doi.org/10.7554/eLife.50374.010 Figure supplement 3. Number of genes composing 50% of the transcriptome throughout the developmental program. DOI: https://doi.org/10.7554/eLife.50374.011 Figure supplement 4. A: Relationship between M. xanthus developmental phases, developmental groups, and a known functional interaction network. DOI: https://doi.org/10.7554/eLife.50374.012 Mun ˜oz-Dorado et al. eLife 2019;8:e50374. DOI: https://doi.org/10.7554/eLife.50374 6 of 27 Tools and resources Microbiology and Infectious Disease 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 upor 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) (Mu ¨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. Aand S-motility genes exhibit different developmental expression profiles The two phases of the M. xanthus lifecycle (predatory growth and development) depend on Aand S-motility engines and their associated regulatory proteins (Pe ´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 Aand 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). 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 Mun ˜oz-Dorado et al. eLife 2019;8:e50374. DOI: https://doi.org/10.7554/eLife.50374 7 of 27 Tools and resources Microbiology and Infectious Disease Figure 4. 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 Figure 4 continued on next page Mun ˜oz-Dorado et al. eLife 2019;8:e50374. DOI: https://doi.org/10.7554/eLife.50374 8 of 27 Tools and resources Microbiology and Infectious Disease 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). 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 straightand 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 (b-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. 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). 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 Figure 4 continued 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. DOI: https://doi.org/10.7554/eLife.50374.020 Mun ˜oz-Dorado et al. eLife 2019;8:e50374. DOI: https://doi.org/10.7554/eLife.50374 9 of 27 Tools and resources Microbiology and Infectious Disease 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). 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 standalone 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 (Mun ˜oz-Dorado et al., 1991; Pe ´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 s 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 (Yagu ¨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. Figure 9 continued 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. DOI: https://doi.org/10.7554/eLife.50374.025 Mun ˜oz-Dorado et al. eLife 2019;8:e50374. DOI: https://doi.org/10.7554/eLife.50374 16 of 27 Tools and resources Microbiology and Infectious Disease Figure 10. 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. Figure 10 continued on next page Mun ˜oz-Dorado et al. eLife 2019;8:e50374. DOI: https://doi.org/10.7554/eLife.50374 17 of 27 Tools and resources Microbiology and Infectious Disease 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. Materials and methods Key resources table Reagent type (species) or resource Designation Source or reference Identifiers Additional information 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 (Km r )W4521 Continued on next page Figure 10 continued 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. DOI: https://doi.org/10.7554/eLife.50374.026 Mun ˜oz-Dorado et al. eLife 2019;8:e50374. DOI: https://doi.org/10.7554/eLife.50374 18 of 27 Tools and resources Microbiology and Infectious Disease Continued Reagent type (species) or resource Designation Source or reference Identifiers Additional information Strain, strain background (Myxococcus xanthus) DK4294 Kroos et al., 1986 MXAN_RS16790; MXAN_ 3464 fmgE::Tn5-lacZ; Tn5 lac (Km r )W4406 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 10 8 cells/ml (optical density at 600 nm [OD 600 ] of 1), and then harvested and resuspended in TM buffer (10 mM Tris-HCl [pH 7.6]; 1 mM MgSO 4 ) to a calculated density of 4.5 10 9 cells/ml (OD 600 of 15). For each time replicate, 200 ml 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 (OD 600 of 1) were harvested by centrifugation as above, resuspended in RNA Protect Bacteria Reagent, and processed in the same manner. RNA extraction 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 ml 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 mg of RNA and incubated at room temperature for 10 min. The reaction was stopped by adding the Mun ˜oz-Dorado et al. eLife 2019;8:e50374. DOI: https://doi.org/10.7554/eLife.50374 19 of 27 Tools and resources Microbiology and Infectious Disease 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 ml of RNasefree 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 mg of RNA in a final reaction volume of 20 ml. 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 ml. 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 mg of DNA, with a minimum concentration of 200 ng/ml. Figure 11. 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, CheYlikeresponse 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. DOI: https://doi.org/10.7554/eLife.50374.027 The following source data is available for figure 11: Source data 1. Regulatory elements included in the developmental groups. DOI: https://doi.org/10.7554/eLife.50374.028 Mun ˜oz-Dorado et al. eLife 2019;8:e50374. DOI: https://doi.org/10.7554/eLife.50374 20 of 27 Tools and resources Microbiology and Infectious Disease 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 R 2 correlation coefficient of the two time course replicates. Log 2 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 b-galactosidase activity For quantitative determination of b-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-Mun ˜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. Microscopy To observe swarm and fruiting bodies, 10 ml 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. Supporting data 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_1biological_replicate_1) and SAMN10135984 (WT_72_2-biological_replicate_2); 96 hr: Mun ˜oz-Dorado et al. eLife 2019;8:e50374. DOI: https://doi.org/10.7554/eLife.50374 21 of 27 Tools and resources Microbiology and Infectious Disease SAMN10135985 (WT_96_1-biological_replicate_1) and SAMN10135986 (WT_96_2biological_replicate_2). Acknowledgements 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. Additional information Funding Funder Grant reference number Author Spanish Government CSD2009-00006 Jose Munoz-Dorado Spanish Government BFU2016-75425-P Aurelio Moraleda-Mun ˜oz National Institutes of Health R35GM124733 Jared M Schrader National Science Foundation IOS 1651921 Penelope I Higgs Ministerio de Educacio ´n, Cultura y Deporte Salvador de Madariaga Program Jose Munoz-Dorado Juana Pe ´rez The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication. Author contributions Jose´ Mun ˜oz-Dorado, Juana Pe´ rez, 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; Aurelio MoraledaMun ˜oz, Francisco Javier Marcos-Torres, Francisco Javier Contreras-Moreno, Acquisition of data; Revising the article critically for important intellectual content; Final approval of the version to be published; Ana Belen Martin-Cuadrado, Jared M Schrader, 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; 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 Author ORCIDs Jose´ Mun ˜oz-Dorado https://orcid.org/0000-0001-7765-5687 Jared M Schrader https://orcid.org/0000-0002-5728-5882 Decision letter and Author response Decision letter https://doi.org/10.7554/eLife.50374.033 Author response https://doi.org/10.7554/eLife.50374.034 Additional files Supplementary files .Transparent reporting form DOI: https://doi.org/10.7554/eLife.50374.029 Mun ˜oz-Dorado et al. eLife 2019;8:e50374. DOI: https://doi.org/10.7554/eLife.50374 22 of 27 Tools and resources Microbiology and Infectious Disease Data availability 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). The following dataset was generated: Author(s) Year Dataset title Dataset URL Database and Identifier Martin-Cuadrado, Ana-Belen 2018 Myxococcus xanthus DK1622 transcriptome (TaxID: 34) https://www.ncbi.nlm. nih.gov/bioproject/ PRJNA493545 NCBI BioProject, PRJNA493545 References Apelian D, Inouye S. 1993. A new putative sigma factor of Myxococcus xanthus.Journal of Bacteriology 175: 3335–3342. DOI: https://doi.org/10.1128/jb.175.11.3335-3342.1993,PMID: 8501037 Bhat S, Boynton TO, Pham D, Shimkets LJ. 2014. Fatty acids from membrane lipids become incorporated into lipid bodies during Myxococcus xanthus differentiation. PLOS ONE 9:e99622. DOI: https://doi.org/10.1371/ journal.pone.0099622,PMID: 24906161 Bode HB, Ring MW, Schwa ¨r G, Altmeyer MO, Kegler C, Jose IR, Singer M, Mu ¨ller R. 2009. Identification of additional players in the alternative biosynthesis pathway to isovaleryl-CoA in the myxobacterium Myxococcus xanthus.ChemBioChem 10:128–140. DOI: https://doi.org/10.1002/cbic.200800219,PMID: 18846531 Bretl DJ, Kirby JR. 2016. Molecular mechanisms of signaling in Myxococcus xanthus Development. Journal of Molecular Biology 428:3805–3830. DOI: https://doi.org/10.1016/j.jmb.2016.07.008,PMID: 27430596 Bretscher AP, Kaiser D. 1978. Nutrition of Myxococcus xanthus, a fruiting myxobacterium. Journal of Bacteriology 133:763–768. PMID: 415048 Campoy S, Fontes M, Padmanabhan S, Corte´ s P, Llagostera M, Barbe´ J. 2003. LexA-independent DNA damagemediated induction of gene expression in Myxococcus xanthus.Molecular Microbiology 49:769–781. DOI: https://doi.org/10.1046/j.1365-2958.2003.03592.x,PMID: 12864858 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. DOI: https://doi.org/10.1126/science.aad2001,PMID: 26965631 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.Journal of Bacteriology 189:169–178. DOI: https://doi. org/10.1128/JB.00792-06,PMID: 17056749 Curtis PD, Shimkets LJ. 2008. Metabolic pathways relevant to predation, signaling, and development. In: Whitworth D. E (Ed). Myxobacteria: Multicellularity and Differentiation. ASM Press. de Hoon MJ, Imoto S, Nolan J, Miyano S. 2004. Open source clustering software. Bioinformatics 20:1453–1454. DOI: https://doi.org/10.1093/bioinformatics/bth078,PMID: 14871861 Diodati ME, Gill RE, Plamann L, Singer M. 2008. Initiation and early developmental events. In: Whitworth D. E (Ed). Myxobacteria: Multicellularity and Differentiation. ASM Press. Dworkin M, Gibson SM. 1964. A system for studying microbial morphogenesis: rapid formation of microcysts in Myxococcus xanthus.Science 146:243–244. DOI: https://doi.org/10.1126/science.146.3641.243,PMID: 141 85314 Fawcett P, Eichenberger P, Losick R, Youngman P. 2000. The transcriptional profile of early to middle sporulation in Bacillus subtilis.PNAS 97:8063–8068. DOI: https://doi.org/10.1073/pnas.140209597,PMID: 1086 9437 Foster HA, Parish JH. 1973. Ribosomes, ribosomal subunits and ribosomal proteins from Myxococcus xanthus. Journal of General Microbiology 75:391–400. DOI: 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. PNAS 115:2484–2489. DOI: https://doi.org/ 10.1073/pnas.1716441115 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. DOI: https://doi.org/10.1186/1471-2180-13-279,PMID: 24304716 Mun ˜oz-Dorado et al. eLife 2019;8:e50374. DOI: https://doi.org/10.7554/eLife.50374 23 of 27 Tools and resources Microbiology and Infectious Disease Giglio KM, Caberoy N, Suen G, Kaiser D, Garza AG. 2011. A cascade of coregulating enhancer binding proteins initiates and propagates a multicellular developmental program. PNAS 108:E431–E439. DOI: https://doi.org/ 10.1073/pnas.1105876108,PMID: 21670274 Gohara DW, Yap MF. 2018. Survival of the drowsiest: the hibernating 100S ribosome in bacterial stress management. Current Genetics 64:753–760. DOI: https://doi.org/10.1007/s00294-017-0796-2,PMID: 2 9243175 Goldman BS, Nierman WC, Kaiser D, 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, et al. 2006. Evolution of sensory complexity recorded in a myxobacterial genome. PNAS 103: 15200–15205. DOI: https://doi.org/10.1073/pnas.0607335103,PMID: 17015832 Hagen DC, Bretscher AP, Kaiser D. 1978. Synergism between morphogenetic mutants of Myxococcus xanthus. Developmental Biology 64:284–296. DOI: https://doi.org/10.1016/0012-1606(78)90079-9,PMID: 98366 Harms A, Maisonneuve E, Gerdes K. 2016. Mechanisms of bacterial persistence during stress and antibiotic exposure. Science 354:aaf4268. DOI: https://doi.org/10.1126/science.aaf4268,PMID: 27980159 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. DOI: https://doi.org/10.1016/j.jinorgbio.2011.11.022,PMID: 22196016 Higgs PI, Hartzell P, Holkenbrink C, Hoiczyk E. 2014. Myxococcus xanthus vegetative and developmental cell heterogeneity. In: Yang Z, Higgs P. I (Eds). Myxobacteria: Genomics, Cellular and Molecular Biology. Caister Academic Press. Hodgkin J, Kaiser D. 1977. Cell-to-cell stimulation of movement in nonmotile mutants of Myxococcus.PNAS 74: 2938–2942. DOI: https://doi.org/10.1073/pnas.74.7.2938,PMID: 16592422 Hoiczyk E, Ring MW, McHugh CA, Schwa ¨r G, Bode E, Krug D, Altmeyer MO, Lu JZ, Bode HB. 2009. Lipid body formation plays a central role in cell fate determination during developmental differentiation of Myxococcus xanthus.Molecular Microbiology 74:497–517. DOI: https://doi.org/10.1111/j.1365-2958.2009.06879.x,PMID: 1 9788540 Holkenbrink C, Hoiczyk E, Kahnt J, Higgs PI. 2014. Synthesis and assembly of a novel glycan layer in Myxococcus xanthus spores. The Journal of Biological Chemistry 289:32364–32378. DOI: https://doi.org/10.1074/jbc.M114. 595504,PMID: 25271164 Inouye S. 1990. Cloning and DNA sequence of the gene coding for the major sigma factor from Myxococcus xanthus.Journal of Bacteriology 172:80–85. DOI: https://doi.org/10.1128/jb.172.1.80-85.1990,PMID: 2104614 Inouye S, Nariya H, Mun ˜oz-Dorado J. 2008. Protein Ser/Thr kinases and phosphatases in Myxococcus xanthus. In: Whitworth D. E (Ed). Myxobacteria: Muticellularity and Differentiation. ASM Press. Islam ST, Mignot T. 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. DOI: https://doi. org/10.1016/j.semcdb.2015.10.033,PMID: 26520023 Jakobsen JS, Jelsbak L, Jelsbak L, Welch RD, Cummings C, Goldman B, Stark E, Slater S, Kaiser D. 2004. Sigma54 enhancer binding proteins and Myxococcus xanthus fruiting body development. Journal of Bacteriology 186:4361–4368. DOI: https://doi.org/10.1128/JB.186.13.4361-4368.2004,PMID: 15205438 Kaiser D. 1979. Social gliding is correlated with the presence of pili in Myxococcus xanthus.PNAS 76:5952– 5956. DOI: https://doi.org/10.1073/pnas.76.11.5952,PMID: 42906 Keseler IM, Kaiser D. 1997. sigma54, a vital protein for Myxococcus xanthus.PNAS 94:1979–1984. DOI: https:// doi.org/10.1073/pnas.94.5.1979,PMID: 9050890 Konovalova A, Wegener-Feldbru ¨gge S, Søgaard-Andersen L. 2012. Two intercellular signals required for fruiting body formation in Myxococcus xanthus act sequentially but non-hierarchically. Molecular Microbiology 86:65– 81. DOI: https://doi.org/10.1111/j.1365-2958.2012.08173.x,PMID: 22834948 Korp J, Vela Gurovic MS, Nett M. 2016. Antibiotics from predatory bacteria. Beilstein Journal of Organic Chemistry 12:594–607. DOI: https://doi.org/10.3762/bjoc.12.58,PMID: 27340451 Kottel RH, Bacon K, Clutter D, White D. 1975. Coats from Myxococcus xanthus: characterization and synthesis during myxospore differentiation. Journal of Bacteriology 124:550–557. PMID: 809426 Kroos L, Kuspa A, Kaiser D. 1986. A global analysis of developmentally regulated genes in Myxococcus xanthus. Developmental Biology 117:252–266. DOI: https://doi.org/10.1016/0012-1606(86)90368-4,PMID: 3017794 Kroos L. 2017. Highly Signal-Responsive gene regulatory network governing Myxococcus Development. Trends in Genetics 33:3–15. DOI: https://doi.org/10.1016/j.tig.2016.10.006,PMID: 27916428 Laub MT, McAdams HH, Feldblyum T, Fraser CM, Shapiro L. 2000. Global analysis of the genetic network controlling a bacterial cell cycle. Science 290:2144–2148. DOI: https://doi.org/10.1126/science.290.5499.2144, PMID: 11118148 Leonardy S, Freymark G, Hebener S, Ellehauge E, Søgaard-Andersen L. 2007. Coupling of protein localization and cell movements by a dynamically localized response regulator in Myxococcus xanthus.The EMBO Journal 26:4433–4444. DOI: https://doi.org/10.1038/sj.emboj.7601877,PMID: 17932488 Li Y, Sun H, Ma X, Lu A, Lux R, Zusman D, Shi W. 2003. Extracellular polysaccharides mediate Pilus retraction during social motility of Myxococcus xanthus.PNAS 100:5443–5448. DOI: https://doi.org/10.1073/pnas. 0836639100,PMID: 12704238 Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R, 1000 Genome Project Data Processing Subgroup. 2009. The sequence alignment/Map format and SAMtools. Bioinformatics 25:2078–2079. DOI: https://doi.org/10.1093/bioinformatics/btp352,PMID: 19505943 Mun ˜oz-Dorado et al. eLife 2019;8:e50374. DOI: https://doi.org/10.7554/eLife.50374 24 of 27 Tools and resources Microbiology and Infectious Disease Li H, Durbin R. 2009. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25: 1754–1760. DOI: https://doi.org/10.1093/bioinformatics/btp324,PMID: 19451168 Lobedanz S, Søgaard-Andersen L. 2003. Identification of the C-signal, a contact-dependent morphogen coordinating multiple developmental responses in Myxococcus xanthus.Genes & Development 17:2151–2161. DOI: https://doi.org/10.1101/gad.274203,PMID: 12923062 Lux R, Shi W. 2005. A novel bacterial signalling system with a combination of a ser/Thr kinase cascade and a his/ Asp two-component system. Molecular Microbiology 58:345–348. DOI: https://doi.org/10.1111/j.1365-2958. 2005.04856.x,PMID: 16194223 Mauriello EM, Mignot T, Yang Z, Zusman DR. 2010. Gliding motility revisited: how do the myxobacteria move without flagella? Microbiology and Molecular Biology Reviews 74:229–249. DOI: https://doi.org/10.1128/ MMBR.00043-09,PMID: 20508248 McBride MJ, Zusman DR. 1989. Trehalose accumulation in vegetative cells and spores of Myxococcus xanthus. Journal of Bacteriology 171:6383–6386. DOI: https://doi.org/10.1128/jb.171.11.6383-6386.1989,PMID: 250 9436 McLaughlin PT, Bhardwaj V, Feeley BE, Higgs PI. 2018. MrpC, a CRP/Fnr homolog, functions as a negative autoregulator during the Myxococcus xanthus multicellular developmental program . Molecular Microbiology 109:245–261. DOI: https://doi.org/10.1111/mmi.13982 Mercier R, Mignot T. 2016. Regulations governing the multicellular lifestyle of Myxococcus xanthus.Current Opinion in Microbiology 34:104–110. DOI: https://doi.org/10.1016/j.mib.2016.08.009,PMID: 27648756 Moraleda-Mun ˜oz A, Carrero-Le´ rida J, Pe´ rez J, Mun ˜oz-Dorado J. 2003. Role of two novel two-component regulatory systems in development and phosphatase expression in Myxococcus xanthus.Journal of Bacteriology 185:1376–1383. DOI: https://doi.org/10.1128/JB.185.4.1376-1383.2003,PMID: 12562808 Morett E, Segovia L. 1993. The sigma 54 bacterial enhancer-binding protein family: mechanism of action and phylogenetic relationship of their functional domains. Journal of Bacteriology 175:6067–6074. DOI: https://doi. org/10.1128/jb.175.19.6067-6074.1993,PMID: 8407777 Mu ¨ller FD, Treuner-Lange A, Heider J, Huntley SM, Higgs PI. 2010. Global transcriptome analysis of spore formation in Myxococcus xanthus reveals a locus necessary for cell differentiation. BMC Genomics 11:264. DOI: https://doi.org/10.1186/1471-2164-11-264,PMID: 20420673 Mu ¨ller FD, Schink CW, Hoiczyk E, Cserti E, Higgs PI. 2012. Spore formation in Myxococcus xanthus is tied to cytoskeleton functions and polysaccharide spore coat deposition. Molecular Microbiology 83:486–505. DOI: https://doi.org/10.1111/j.1365-2958.2011.07944.x,PMID: 22188356 Mu ¨ller S, Strack SN, Ryan SE, Shawgo M, Walling A, Harris S, Chambers C, Boddicker J, Kirby JR. 2016. Identification of functions affecting Predator-Prey interactions between Myxococcus xanthus and Bacillus subtilis.Journal of Bacteriology 198:3335–3344. DOI: https://doi.org/10.1128/JB.00575-16,PMID: 27698086 Mun ˜oz-Dorado J, Inouye S, 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. DOI: https://doi.org/10. 1016/0092-8674(91)90372-6,PMID: 1835671 Mun ˜oz-Dorado J, Higgs PI, Elı ´as-Arnanz M. 2014. Abundance and complexity of signaling mechanisms in myxobacteria. In: Yang Z, Higgs P. I (Eds). Myxobacteria: Genomics Cellular and Molecular Biology. Caister Academic Press. Mun ˜oz-Dorado J, Marcos-Torres FJ, Garcı ´a-Bravo E, Moraleda-Mun ˜oz A, Pe´ rez J. 2016. Myxobacteria: moving, killing, feeding, and surviving together. Frontiers in Microbiology 7:781. DOI: https://doi.org/10.3389/fmicb. 2016.00781,PMID: 27303375 Nan B, McBride MJ, Chen J, Zusman DR, Oster G. 2014. Bacteria that glide with helical tracks. Current Biology 24:R169–R173. DOI: https://doi.org/10.1016/j.cub.2013.12.034,PMID: 24556443 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. Molecular Microbiology 52:273–283. DOI: https://doi.org/10.1111/j.1365-2958.2003.03972.x, PMID: 15049826 Nariya H, Inouye S. 2003. An effective sporulation of Myxococcus xanthus requires glycogen consumption via Pkn4-activated 6-phosphofructokinase. Molecular Microbiology 49:517–528. DOI: https://doi.org/10.1046/j. 1365-2958.2003.03572.x,PMID: 12828646 Nariya H, Inouye M. 2008. MazF, an mRNA interferase, mediates programmed cell death during multicellular Myxococcus development. Cell 132:55–66. DOI: https://doi.org/10.1016/j.cell.2007.11.044,PMID: 18191220 O’Connor KA, Zusman DR. 1991a. Analysis of Myxococcus xanthus cell types by two-dimensional polyacrylamide gel electrophoresis. Journal of Bacteriology 173:3334–3341. DOI: https://doi.org/10.1128/jb.173.11.33343341.1991,PMID: 1904431 O’Connor KA, Zusman DR. 1991b. Behavior of peripheral rods and their role in the life cycle of Myxococcus xanthus.Journal of Bacteriology 173:3342–3355. DOI: https://doi.org/10.1128/jb.173.11.3342-3355.1991, PMID: 1904432 O’Connor KA, Zusman DR. 1991c. Development in Myxococcus xanthus involves differentiation into two cell types, peripheral rods and spores. Journal of Bacteriology 173:3318–3333. DOI: https://doi.org/10.1128/jb. 173.11.3318-3333.1991,PMID: 1904430 Pe ´rez J, Castan ˜eda-Garcı ´a A, Jenke-Kodama H, Mu ¨ller R, Mun ˜oz-Dorado J. 2008. Eukaryotic-like protein kinases in the prokaryotes and the myxobacterial kinome. PNAS 105:15950–15955. DOI: https://doi.org/10.1073/pnas. 0806851105,PMID: 18836084 Mun ˜oz-Dorado et al. eLife 2019;8:e50374. DOI: https://doi.org/10.7554/eLife.50374 25 of 27 Tools and resources Microbiology and Infectious Disease