Transcriptome Analysis of Polyhydroxybutyrate Cycle Mutants Reveals Discrete Loci Connecting Nitrogen Utilization and Carbon Storage in Sinorhizobium meliloti

The ability of bacteria to store carbon and energy as intracellular polymers uncouples cell growth and replication from nutrient uptake and provides flexibility in the use of resources as they are available to the cell. The impact of carbon storage on cellular metabolism would be reflected in global transcription patterns. By investigating the transcriptomic effects of genetically disrupting genes involved in the PHB carbon storage cycle, we revealed a relationship between intracellular carbon storage and nitrogen metabolism. This work demonstrates the utility of combining transcriptome sequencing with metabolic pathway mutations for identifying underlying gene regulatory mechanisms.

While the regulatory mechanisms that govern PHB accumulation remain unclear, evidence suggests that the PHB cycle is linked to the regulation of other carbon storage mechanisms (exopolysaccharide [EPS] and glycogen production), as well motility and symbiosis. For example, a phbC mutation has been shown to result in decreased EPS levels (succinoglycan [EPSI] and galactoglucan [EPSII]) and increased glycogen levels, as well as a decreased ability to compete for nodule infection (9,17). Mutations in chvI lead to hypermotility and altered PHB accumulation (18).
Glycogen is produced in S. meliloti as a reserve carbon store, similar to PHB (15). There is not a clear mechanism for the allocation of carbon between glycogen synthesis and PHB synthesis in S. meliloti. In cyanobacteria, it has been proposed that carbon is diverted to glycogen (over PHB) during nitrogen stress (19). Unlike PHB, glycogen forms granules that are water soluble, and they are generally much smaller than PHB granules (20). Glycogen is synthesized through the production of ADP-glucose by pyrophosphorylase, GlgP, followed by the sequential addition of these molecules to a growing ␣-1,4-glucan chain by glycogen synthase (GlgA1) and glycogen branching enzyme (GlgB1) (21). Glycogen can be debranched to its individual units by the glycogendebranching enzyme GlgX1 to free up carbon during times of metabolic activity. Disruption of the glycogen synthase gene glgA1 results in decreased PHB levels and increased EPS levels compared to the wild type (15). Strains unable to synthesize glycogen are able to form nodules and fix nitrogen, but they exhibit a noticeable delay in nodulation and a reduction in competitive ability (15).
Previous transcriptomics studies of Ralstonia eutropha (Cupriavidus necator) examined expression differences between the wild type (strain H16) and PHB synthesis mutants, as well as between growth phases. The PHB synthesis mutants exhibited increased expression of various PHB cycle and related genes (phasins), appeared to be defective in fatty acid metabolism, and showed downregulation of two cbb operons that are involved in CO 2 assimilation (22). In another study, during PHB accumulation, repression of genes involved in central metabolism was observed in addition to active transcription of cbb and ␤-oxidation genes (23). This work also showed consistently high transcription of the PHB cycle-associated genes phaR, phaC1, and phaB1 and some phasin genes (23).
Previous studies have helped to establish knowledge of transcriptomic regulation in S. meliloti. Early work involved both macroarray and microarray investigations of symbiotic development, growth under micro-oxic conditions, acidic adaptation response, cell cycle and environmental signaling, Clr regulation, and the nitrogen stress response (24)(25)(26)(27)(28)(29)(30)(31)(32). More recently, transcriptome sequencing (RNA-seq) analysis has contributed to the general mapping of transcription start sites (TSS) and regulatory motifs. To our knowledge, the only specific application of RNA-seq in S. meliloti is to Hfq regulation (33,34). Here, to investigate the influence of intracellular carbon storage mechanisms on global transcriptional regulation, we performed RNA-seq transcriptomic analyses of individual S. meliloti mutants in which each of the key reactions of the PHB cycle are disrupted. Given the role of glycogen as an intracellular carbon store, we also assessed a glycogen synthase gene (glgA1) mutant. We reasoned that by systematically examining the transcriptional effects of multiple PHB cycle disruptions, it might be possible to differentiate core regulatory determinants of the PHB pathway from pleotropic noise that might arise from single gene knockouts. In addition, by using a glycogen synthase mutant, we can differentiate the effects specific to PHB storage mutants from effects caused by mutation in any carbon storage system. By comparing the transcriptomic profiles of the different mutants, we identified core loci of differential transcriptional activity that reside on the pSymA megaplasmid and contain numerous genes involved in nitrogen utilization, nitrogen fixation, and denitrification. Analysis of these genes and their regulatory regions suggests a relationship between intracellular carbon accumulation and the regulation of nitrogen metabolism genes. We are thus continuing the use of RNA-seq technology to delve more deeply into specific regulatory pathways in S. meliloti.

RESULTS AND DISCUSSION
Rm1021 shows large-scale significant differences in the transcription of carbon, nitrogen, and symbiosis-associated genes under nitrogen-limited conditions. (i) Establishment of conditions for differential PHB accumulation. It has been established that S. meliloti accumulates substantial amounts of PHB during growth on YM medium, which is carbon rich due to high levels of mannitol and contains limiting nitrogen (16). We first determined that increasing the amount of yeast extract 20-fold resulted in the absence of detectable PHB accumulation, likely because of relief of the unbalanced carbon/nitrogen ratio. We therefore adopted the original and yeast extractenhanced versions of YM medium as PHB-accumulating and -nonaccumulating medium conditions and termed these conditions nitrogen limited and balanced, respectively, in line with the ratio of carbon and nitrogen levels. Nitrogen-limited conditions consistently resulted in the wild-type strain accumulating substantial PHB at stationary phase, whereas balanced conditions always resulted in undetectable PHB levels. Glycogen is also known to be produced in YMB medium, which is the base medium that was modified for the nitrogen-limited conditions in this work (15).
Gene expression profiles were determined at mid to late log phase under both conditions (see Fig. S1A in the supplemental material). Comparison of the gene expression profiles of parental strain Rm1021 under these two conditions revealed within-strain (wild type to wild type) differential expression of 2,335 genes in total, with 1,096 that were upregulated under nitrogen-limited conditions, showing increased transcript abundance, and 1,239 that were downregulated. It should be cautioned, however, that this effect is not likely to be solely related to carbon storage, as the differences in the yeast extract concentration likely influenced numerous metabolic pathways.
An overview of the differential expression in Rm1021 under balanced and nitrogenlimited conditions across the coordinates of the genome (Fig. 2) demonstrates an even distribution of effect across the three replicons, the megaplasmid pSymA, the chromid pSymB, and the chromosome. Key genes of interest are indicated and discussed in detail below.
(ii) Differentially expressed carbon storage genes. Of note, the acetoacetyl-CoA synthetase-encoding gene acsA2 (accession no. SMc00774; see below), was the only PHB cycle gene differentially expressed under the two conditions. Hence, the PHB cycle genes are not transcriptionally regulated in response to environmental changes that affect PHB accumulation, thus implying that regulation of PHB cycle activity and PHB accumulation occurs at the translational or posttranslational level. We found that two glycogen synthesis gene transcripts, glgA1 and glgX1, differed between the conditions, with increased transcript abundance under nitrogen-limited conditions, indicating that they may be transcriptionally regulated in response to the extracellular carbon and/or nitrogen supply.
The downregulation of the acsA2-encoded acetoacetyl-CoA synthetase raises the possibility that this enzyme represents a key control point in the PHB cycle. However, since acsA2 mutants are not known to hyperaccumulate PHB, it may alternatively reflect the lower levels of available acetoacetate substrate, considering that flux through the degradation pathway would be reduced under these conditions. Why this would be different for acsA2 than other PHB cycle genes is not known and bears further investigation.
The two phasin genes phaP1 (accession no. SMc00777) and phaP2 (accession no. SMc02111) were upregulated under nitrogen-limited conditions, while there was no significant difference in aniA (accession no. SMc03880), phaR in other species, under the two conditions. This was not unexpected, as PhaP1 and PhaP2 are well established as key proteins in the formation of PHB granules (17). In R. eutropha, phasin genes were also found to be upregulated in response to mutation of the PHB synthesis genes (22). Previous studies have implicated aniA in the regulation of PHB and EPS synthesis under microaerobic conditions (35). The lack of difference in aniA transcription suggests that AniA's activity is not transcriptionally regulated, as AniA has been shown to repress the transcription of phasin genes in several species, although this has not been directly shown in S. meliloti, and here the phasin genes are upregulated (36)(37)(38). Insertional inactivation of aniA in S. meliloti Rm41 results in decreased PHB accumulation, demonstrating a role in the regulation of PHB production (35). Recently, PhaR in Bradyrhizobium diazoefficiens was shown to be a global regulator of excess carbon allocation and symbiosis through control of the fixK2 gene (37). The increase in phasin gene expression under nitrogen-limited conditions observed here could be mediated by AniA or another transcriptional regulator.
(iii) Differentially expressed nitrogen-associated genes. The most strongly upregulated genes under nitrogen-limited conditions were those involved in amino acid and nitrogen uptake and utilization. This most certainly reflects the very different availability of nitrogenous compounds under the two conditions because of the yeast extract concentration differences. Included among the genes with the strongest transcript level upregulation were several that encode amino acid uptake and binding proteins, as well as nitrate uptake genes. The amino acid uptake genes include the aap operon (SMc02118 to SMc02121) and the liv (bra) operon (SMc01946 to SMc01951), both involved in branched-chain amino acid uptake, as well as the hut (SMb21163 to SMb21166) and his (SMc00669 to SMc00674) operons, involved in histidine and proline uptake (39,40). These transcriptional changes are all consistent with transcriptional regulation (induction or removal of repression) of amino acid uptake genes in response to the scarcity of nitrogenous compounds. Interestingly, this finding has similarities to another study that detected significant differences in the transcription of amino acid uptake genes in a cyaJ overexpression mutant of S. meliloti (32) (further discussed below). Under nitrogen-limited conditions, exogenous levels of amino acids are expected to be much lower than under balanced conditions, which would influence the expression of genes involved in the uptake and utilization of amino acids.
Studies have demonstrated that when S. meliloti is grown under nitrogen-limited conditions, characterized by a low glutamine/␣-ketoglutarate ratio, the cell activates the nitrogen stress response, mediated by GlnD through GlnK and GlnB (41). Activated GlnK and GlnB contribute to the activation of NtrC to increase the catabolism of nitrogen-containing compounds, freeing up ammonium. Both glnK and glnB (SMc03806 and SMc00947) were upregulated under nitrogen-limited conditions, while ntrC and glnD were not. glnK was the fourth most upregulated transcript under nitrogen-limited conditions, suggesting that it could be a key component of the observed cellular response. This finding suggests that the nitrogen stress response is transcriptionally regulated through the upregulation of glnK and glnB, whose protein products then trigger catabolism (through the activation of NtrC) of intracellular nitrogenous com-pounds to free nitrogen stores in the cell. Previous work with Azospirillum brasilense has implicated NtrB and NtrC in the regulation of PHB synthesis in response to the ammonia concentration in the medium (42). Consistent with the activation of NtrC, we observed upregulation of the following known transcriptional targets of NtrC: glnII, which encodes glutamine synthetase (GS; SMb20745), glnE (SMc02368), which encodes an adenylyltransferase that inactivates GS, and amtB (SMc03807), which encodes an ammonium transporter.
The glx operon (SMc02610 to SMc02612), which codes for glutamate synthase, was also upregulated, while glxA (SMc02609), the regulator of the operon, was downregulated. This suggests that glxA acts as an inhibitor of glx transcription. Immediately upstream of the glx operon is the sox operon (soxGA1BD), which was starkly upregulated under nitrogen-limited conditions. The sox operon encodes sarcosine oxidase, which degrades sarcosine, an intermediate in glycine synthesis and degradation (43).
The GSIII-encoding gene glnT (SMc02613) was also upregulated, although it had previously been reported to be expressed only when GSI and GSII are absent (44). The expression of GSI-encoding glnA did not change under the two conditions, and out of three putative GSs, one showed no difference (SMc00762), one was upregulated (SMc01594), and the other was downregulated (SMc02352). This suggests areas in which experimental investigation would be fruitful, as these results are not consistent with the existing body of literature or in some cases contradict previous findings, as discussed above.
The dctA gene (SMb20611), located on the chromid, showed a strong decrease in transcription under nitrogen-limited conditions. DctA is a dicarboxylic acid transporter used by the cell during symbiosis to transport in dicarboxylic acids, such as malate, provided by the host plant. This downregulation may be a response to high intracellular carbon levels, as further import of organic acids into the cell would increase the intracellular carbon content. It is possible that this is not specific to dctA but instead is due to a global repression effect on multiple carbon uptake pathways, even those that would normally be inactive.
EPSs are critically important to the interaction of S. meliloti with its plant host, leading to the development and infection of nitrogen-fixing root nodules. The synthesis of EPSs has long been associated with carbon-sufficient conditions under which growth is hindered by low levels of other essential nutrients such as nitrogen (47). These conditions also promote PHB accumulation. In S. meliloti Rm1021, this is manifested in a mucoid colony morphology on agar plates (48). Consistent with this, we observed that several genes involved in the succinoglycan and galactoglucan EPS synthesis pathways were significantly upregulated under nitrogen-limited conditions. This includes many of the succinoglycan synthesis genes in the 35-kb region bounded by SMb20932 and exoP (SMb20961). Interestingly, the regulatory genes exsF (SMb20934), exsI (SMb20935), and exsB (SMb20940) were downregulated. Similarly, many of the genes in the 27-kb galactoglucan synthesis region between wgeH (SMb21307) and wgaJ (SMb21327) were significantly upregulated. This suggests that regulation of EPS accumulation is transcriptional and is in response to growth conditions. Although Rm1021 contains a null expR mutation, we observed upregulation of sinI (SMc00168) and sinR (SMc00170), which encode components of a quorum-sensing system that is known to positively regulate EPS synthesis (49,50).
The operon containing msbA2 (SMb21191) was observed to be significantly upregulated. MsbA2 has been suggested to be involved in EPS production and is essential for the establishment of a symbiosis between S. meliloti and alfalfa (51). Mutations in msbA2 disrupt the infection leading to root nodule symbiosis, and the operon has also been shown to be regulated by the exoR-exoS-chvI regulatory system, which is known to also be required for succinoglycan and galactoglucan synthesis (28,52). Interestingly, the expression of exoS (SMc04446) was significantly downregulated in Rm1021 under nitrogen-limited conditions.
(v) Other metabolism. Surprisingly, genes for phosphate and phosphonate uptake and utilization (SMb21174, SMb21175, SMb21176, and SMb2177) and purine metabolism (SMb21286, SMb21287, and SMb21288), were strongly downregulated. Rm1021 has a mutation in pstC that causes the cell to constantly respond as if it were in a state of phosphate stress via constitutive expression of the pho operon (53). Downregulation of the pho operon and purine metabolism genes is counterintuitive, as the extracellular phosphate and purine concentrations would be low in nitrogen-limited medium because of the decreased amount of yeast extract present. Also of note is the downregulation of many genes involved in the energy currency of the cell, such as those encoding components of the electron transport chain (SMb21488 and SMb21489), as well as those involved in core functions of transcription, translation, and homologous recombination (e.g., SMc00565, SMc01290, SMc01291, SMc02101, SMc01318, SMc01319, SMc01378, SMc01379, SMc03779, SMc01243, SMc00335, and SMc02692). This may reflect a decrease in central metabolism exhibited by the cell in response to nitrogen limitation.
The rhizobactin 1021 synthesis operon on pSymA (SMa2400 to SMa2410) was downregulated under nitrogen-limited conditions. Rhizobactin is a siderophore that captures iron for transport into the cell (54). This operon is known to be repressed under high-iron conditions, as might be expected under balanced conditions (55). The downstream genes rhrA and rhtA (SMa2412 and SMa2414) were similarly downregulated. RhtA is a known outer membrane receptor protein for rhizobactin 1021, while RhrA is a positive transcriptional regulator of rhrA and the rhizobactin operon (55). Both conditions have equal molarity of FeCl 3 · 6H 2 O, but there is 20 times as much yeast extract, which is known to contain trace amounts of iron, under balanced conditions as under nitrogen-limited conditions. It is therefore surprising to observe decreased transcription of the rhizobactin synthesis operon under nitrogen-limited conditions. Previous studies provide an explanation for this phenomenon, as it has been suggested that the differential regulation of the gene for rhizobactin and other iron uptake genes is attributable to a general stress response (56,57).
(vi) Coregulation of carbon and nitrogen pathways. Although none of the PHB cycle genes were upregulated under nitrogen-limited conditions, the xdhA2 and xdhB2 genes showed increased expression. Notably, these genes are involved in hypoxanthine degradation and are immediately adjacent to the bdhA gene, which encodes the key PHB degradation enzyme D-3-hydroxybutyrate dehydrogenase. In rhizosphere colonization and infection, it has been hypothesized that a regulatory association exists between PHB degradation and the degradation of nitrogen-containing metabolites (hypoxanthine) (58). Consistent with this idea, expression of the predicted LysR-type regulator (SMb20847) downstream of xdhA2-xdhB2 was also upregulated, raising the possibility that it is somehow involved in the regulation of this operon. It was previously reported that SMb20847 in a nifA mutant background was upregulated compared to that in the wild-type background during symbiosis, which establishes this potential regulator's involvement in symbiosis (26). Increased expression of these genes under nitrogen-limited conditions would contribute to increased breakdown of hypoxanthine in the cell to free up nitrogen for use in critical cell activities.
Another previous study of interest examined the phosphate starvation response by comparing the transcriptomic profiles of Rm1021, Rm2011, and a phoB mutant under limited-phosphate conditions (57). That study identified 234 genes that were differentially expressed in the phoB mutant and wild-type strains (57). In comparison to our data set, 155 of these genes overlap those that were differentially expressed in the parental strain under nitrogen-limited conditions (Fig. S1B). The common genes include iron uptake genes (rhizobactin genes and others), EPS synthesis genes, xanthine dehydrogenase genes (xdhA2 and xdhB2), phosphate and other ABC transporter genes, ribosomal protein genes, flagellar and chemotaxis protein genes, arginine transport genes, and 40 hypothetical genes.
Comparing the three transcriptomes collectively, only 17 genes were commonly differentially expressed (Fig. S1B). These genes include four EPS synthesis genes (exoF1, exoY, exoN, and exoW), a cytochrome gene (cycF), and two flagellar protein genes (flgB and flgG). It is difficult to interpret the relevance of this set of overlapping genes, which may speak more to the limited depth of microarray data than a lack of biological importance.
These results demonstrate that there are numerous differentially regulated genes under balanced versus nitrogen-limited conditions. Common themes include expression changes in genes related to nitrogen uptake and denitrification, nitrogen fixation, and nitrogen stress response and PHB-associated genes. This is in contrast with the observed lack of transcriptional change in genes directly associated with the two main carbon storage systems, PHB and glycogen. These findings demonstrate widespread transcriptional changes in nitrogen metabolism genes in the cell in response to a change in the environmental nitrogen profile.
Genetic disruption of the PHB cycle leads to strong transcriptional downregulation of loci on pSymA. (i) PHB cycle mutants. By applying transcriptomics to preconstructed phbA, phbB, bdhA, acsA2, and glgA1 mutants, as well as newly generated phbC, phbAB, and phaZ mutants (see Materials and Methods), we analyzed global patterns of gene expression at multiple steps in the PHB pathway. We also included a glycogen synthase gene (glgA1) mutant as another carbon storage system mutant for comparison. For each mutant, we first confirmed the absence or reduction of the transcript corresponding to the mutated gene, according to the deletion or insertion nature of the mutation. For all of the genes, we then compared the transcript abundance in each mutant with that in the parental strain, which we refer to as betweenstrain differential expression. We are using this term to differentiate transcript abundance comparisons of the same strain under different conditions (within-strain differential expression) and comparisons between different strains under the same conditions, i.e., parental to mutant (between-strain differential expression). Earlier we discussed within-strain differential expression, the log 2 -fold change between balanced and nitrogen-limited conditions in the wild type. In this section, use of the term differential expression can be assumed to refer to the difference between strains.
When evaluating between-strain differential expression under each of the two conditions for each mutant, some general trends were recognized. The number of transcripts ranged from 120 (55 up/65 down) in the phbC mutant to Ͼ1,500 (696 up/812 down) in the glgA1 mutant, both under balanced conditions (Fig. 3A). Interestingly, the acsA2 and glgA1 mutants both appeared to show larger numbers of differ- entially expressed genes under balanced conditions and a larger difference in the number of differentially expressed genes between the two conditions than the other mutants, which tended to be similarly affected by the two conditions. We directly examined the relationship between the differentially expressed transcripts in each mutant under the two conditions (Fig. 3B). For each mutant, the overlap of expression under the two conditions was relatively small. This underlines the relevance of examining transcriptional responses under the two different environmental conditions, only one of which results in PHB accumulation.
The genome-wide distribution of log 2 -fold changes in genes for each mutant and background was determined. There was consistency in the frequency distribution for each mutant between the conditions (Fig. S1C). However, when focusing on the genes showing the largest effect size (largest positive and negative log 2 -fold changes), substantial differences were observed in the mutants between conditions (Fig. 3C). Under nitrogen-limited conditions, a larger number of genes show extreme log 2 -fold values (in terms of both up-and downregulation). The most notable effect was the increase in genes showing negative log 2 -fold values. This trend excludes the acsA2, glgA1, and bdhA mutants, which still showed a high frequency of strong effect but less of a notable difference between the two conditions.
(ii) Expression patterns in each of the mutant strains. The patterns of expression in the eight mutants were compared and visualized by using a clustered heat map ( Fig. 4A; enlarged in Fig. S2). The clustering reflects the clear distinction between nitrogen-limited and balanced conditions. Under both conditions, the acsA2 and bdhA degradation pathway mutants were tightly clustered together. Under balanced conditions, the phbA, phbB, and phbAB synthesis mutants were more closely associated with each other than with phbC. Surprisingly, the phbC and phaZ mutants had similar transcriptional patterns. The clustering of the phbC and phaZ mutants under balanced conditions possibly reflects their unique location in the PHB cycle, at the juncture between PHB synthesis and degradation.
In contrast, under nitrogen-limited conditions, the phbA, phbB, and phbC mutants clustered together while the phaZ and phbAB mutants formed a loose cluster with the glgA1 mutant. The discordant responses of the phbA, phbB, and phbAB mutants under nitrogen-limited conditions were unexpected, considering that phbA and phbB are organized in a single operon and the mutants were expected to behave similarly. The phbA and phbAB mutants would be expected to accumulate acetyl-CoA similarly, as the phbA mutation is nonpolar to phbB, but the resulting transcriptional differences do not correlate with this. A potential hypothesis could be that acetoacetyl-CoA is produced by other enzymes in the cell, acting as the substrate for the PhbB enzyme in the phbA mutant. This is unlikely, as the phbA mutant has been shown to be unable to produce detectable accumulation of PHB despite its deletion mutation being nonpolar on phbB (62).
This initial analysis of the RNA-seq data suggests important differences in the patterns of gene expression in the mutants under the two conditions.
(iii) Commonly differentially expressed genes are overrepresented in distinct regions of pSymA. To better visualize patterns of expression, the log 2 -fold changes in genes in each mutant background were mapped to their chromosomal locations and compared ( Fig. S3 and S4). There was a clear distinction in the transcriptional response pattern of the mutants under nitrogen-limited conditions and that seen under balanced conditions. Strong common transcriptional changes were most evident under nitrogen-limited conditions in loci on pSymA (Fig. S3). Given the remarkable effect observed under nitrogen-limited conditions, we focused our attention on these common genes. The differentially expressed genes on pSymA are clustered in three distinct regions, the largest being a 120-kb segment bounded by SMa1077 and SMa1297, with 75% of the genes differentially expressed in at least one of the mutants. The two smaller regions, 37 and 7 kb, respectively, are separated from each other by 40 kb and bounded by SMa0631 and SMa0697 (50% of the genes differentially expressed in at least one mutant) and by SMa0760 and SMa0766 (100% of the genes differentially expressed in at least one mutant). Intriguingly, the transcriptional response in this region was exclusive to the PHB cycle mutants and did not include the glgA1 mutant, which implies that the response is specific to mutation in the PHB cycle.
To connect the discovery of these loci to our initial analysis in the Rm1021 background, we took a closer look at the differential expression in Rm1021 under balanced and nitrogen-limited conditions. We focused on the 120-kb region, as it has the largest number of genes, at 124. Of the genes in the 120-kb region, the expression of 52 differed under the two conditions. Most of these genes were upregulated in the wild type under nitrogen-limited conditions, and many of them were strongly differentially expressed in the mutants compared to the wild type (Fig. S5). The genes upregulated in the wild type under nitrogen-limited conditions tend to be downregulated in the mutants compared to the parental strain under the same conditions. We hypothesize that this may be due to an inability or impairment of the mutants to respond to nitrogen limitation. This would be an interesting line of inquiry for future work.
The 120-kb region contains a smaller number of genes that were downregulated than of genes that were upregulated in Rm1021 under nitrogen-limited conditions. The nine downregulated genes were SMa1099, SMa1182, SMa1183, SMa1184, SMa1213, SMa1214, SMa1269, SMa1272, and SMa1273. These genes encode a hypothetical protein with sequence similarity to the cytochrome c oxidase genes, the nos operon (although the nosR regulator was upregulated), fixQ1/fixP1, and the nor operon. These downregulated genes may reflect the metabolism of the cell under nitrogen-limited conditions, where denitrification would be unlikely to be supported. In contrast to upregulated genes, the strongly downregulated genes tend to show similar responses in the parental and mutant backgrounds (Fig. S5).
The 120-kb region is rich in genes that encode proteins involved in nitrogen utilization, symbiotic nitrogen fixation, and denitrification. This includes the FixLJ two-component regulatory system that regulates N 2 fixation in symbiosis in response to oxygen levels, through its adjacent Fix cluster, including the FixK1 transcriptional regulator. Also included are the NirK nitrite reductase, Nap periplasmic nitrate reductase, Nos nitrous oxide reductase, and Nor nitric oxide reductase operons, all of which have been implicated in denitrification (63). No genes in this region have previously been associated with PHB metabolism or with a regulator known to act on targets other than nitrogen utilization genes.
The 37-kb region from SMa0631 to SMa0697 is located downstream of the fixI2S2 operon and contains mostly conserved/hypothetical genes, including SMa0662, which is a putative Fnr/Crp regulator. This region also includes the arc operon, which encodes arginine transport and metabolism proteins, consistent with the effects of nitrogen levels on PHB regulation. Also in this region is SMa0669, which is predicted to encode an HlyD family protein that contains a biotin binding domain. Biotin has previously been associated with PHB accumulation (64).
The smallest pSymA region exhibiting differential regulation was a single sevengene operon, SMa0760 to SMa0766, which spans the region from fixT2 to fixP2. This operon is in a region that includes nod genes, which are involved in the plant-microbe signaling that precedes symbiosis. The fix genes are involved in nitrogen fixation and specifically encode parts of the cytochrome c oxidase protein. fixNOQP1 mutants exhibit reduced nitrogen fixation, but redundancy of the fixNOQP operon prevents the abrogation of nitrogen fixation (65). The involvement of fixT2 here is also of interest. FixT is an antikinase protein that can affect the transcription of its regulator, FixK, by affecting the phosphorylation of FixL (66).
All three regions contain numerous genes related to nitrogen metabolism. A consistent pattern among all of these regions in the mutants was lower expression in the mutant than in the wild type under nitrogen-limited conditions, while many of these genes were upregulated compared to the level in the parental strain under balanced conditions. On the basis of the strong transcriptional response in nitrogen-associated genes, we hypothesize that obstruction of the PHB cycle disrupts the ability of the cell to respond normally to nitrogen limitation, indicating that carbon storage is involved in the regulatory response to nitrogen limitation.
There are previously discovered and described regulators in these regions. It has been established that the FixL/FixJ two-component regulatory system regulates fixK1 transcriptionally in S. meliloti and is central to a signal cascade leading to nitrogen fixation (67). Further work has indicated that Hfq is involved in the regulation of fixLJ transcription (68). FixK1 is known to downregulate its own transcription through fixT1 activation (69). The transcription of fixK1 itself was downregulated in the mutants under nitrogen-limited conditions, suggesting that there is inhibition of fixK1, which then prevents it from initiating the regulatory cascade in the region. The second copy, fixK2, is found in the 7-kb region. The fixK2 gene in Bradyrhizobium japonicum has been studied more extensively, but little is known about the role of fixK2 in S. meliloti. NnrR, which has been shown to act as a transcriptional activator of other genes in the region involved in the cellular response to NO, is also present in the region (30).
The localized effect of gene expression changes in each of these regions is common among the PHB mutant strains that we tested but, importantly, is not shared by the glycogen synthase mutant, indicating that the effect appears to be specific to disruption of the PHB carbon storage pathway rather than disruption of carbon storage in general. Some of the mutants do not produce PHB, while others accumulate significant amounts of PHB. The shared effect illuminates the difference between the function of intracellular PHB and the presence of PHB in the cell. If the cell were responding to the presence of the PHB polymer only, the effect would differ substantially between synthesis and degradation mutants. Instead, what was observed here is that disruption on either side of the PHB cycle resulted in a similar effect on transcription in the pSymA regions indicated. This suggests that the effect is related to a fully functional PHB cycle rather than the physical presence of a PHB granule.
A discovered motif is related to the transcriptional response in pSymA. (i) Identification of putative motifs associated with differential expression. Hierarchicalcluster analysis of the expression in the mutants produced two major clusters of genes, with one cluster further subdivided into five more clusters (Fig. 4A). The genes in each cluster represent those that show transcriptional patterns that are common to the mutants, possibly because of shared regulation. To identify potential regulatory sequences associated with the observed differential expression, de novo motif discovery was performed. Specifically, regions 200 bp upstream of each gene within each cluster were analyzed with MEME. All of the motifs discovered by MEME with E values below 1 ϫ 10 Ϫ15 are reported (Fig. S6).
Several potential regulatory motifs were identified upstream of the genes in each expression cluster (Fig. S6). Cluster A showed the strongest transcriptional response in all of the mutants under nitrogen-limited conditions (Fig. 4B). The most significant (E value of 3.3 ϫ 10 Ϫ80 ) and biologically intriguing motif in this cluster is shown in Fig. 4C. This motif was compared to a library of known motifs by using Tomtom, which identified a close match to the previously identified Fnr/Crp family FixK transcription factor binding motif (Fig. 4C). We therefore will refer to the cluster A motif as the discovered motif. The first transcription factor and motif of this family was discovered in Escherichia coli and termed Fnr (fumarate and nitrate reductase). The Fnr regulatory system plays a role in anaerobic regulatory responses (67,70,71). More recent work has found Fnr family genes in Rhizobiales and termed this subgroup of proteins FixK. FixK proteins are a subset of the Fnr/Crp regulatory family, and their identified motifs, while varied, generally show the same pattern as Fnr proteins. Interestingly, multiples of these homologs have been previously identified in S. meliloti, as FixK1, FixK2, and NnrR, and all are involved in the regulation of nitrogen fixation during symbiosis and the response to NO (67,72). The FixK consensus sequence has been identified in previous work as TTGAT-N 4 -ATCAA (69). NnrR also has a consensus binding sequence, identified as TTG-N8-CAAA, that has been previously hypothesized to share homology with FixK and may, in fact, overlap the FixK consensus sequence (30). Both FixK1 and NnrR are homologs of the Fnr transcriptional activator in E. coli (72). This suggests that a system closely resembling the Fnr regulatory system in E. coli may be responsible for the observed transcriptional response in the PHB cycle mutants under nitrogen-limited conditions.
The density of the discovered motif across the entire genome was mapped after the discovered motif was identified across the genome by using the coordinates of each motif location and compared against the expression profile of the PHB synthesis mutants and the glgA1 mutant under nitrogen-limited conditions (Fig. 5). This showed a strong correlation between the density of occurrences of the discovered motif in the genome and the transcriptional effect seen in the PHB cycle mutants under nitrogenlimited conditions. That is, the region with the highest motif density overlaps the identified 120-kb region of pSymA associated with the strong decrease in transcript abundance. This finding supports the idea that the hypothesized mechanism of Fnrmediated regulation is responsible for the observed transcriptional response in the pSymA regions. This correlation of density and transcriptional effect does not apply to the glgA1 mutant, suggesting that this region is regulated in response to PHB and not glycogen accumulation. In addition, it has been previously demonstrated that bacterial transcription factors tend to act locally, which suggests that the relevant Fnr/Crp regulator may be located in the region itself (73,74).
(ii) The discovered motif colocates with differential transcription in the PHB cycle mutants under nitrogen-limited conditions. The overrepresentation of the discovered motif in regions with a strong transcriptional effect is compelling evidence of a putative regulatory system that is influenced by the PHB cycle. To further investigate this, we identified all of the instances of the predicted motif across the 120-kb region and examined the log 2 -fold changes in the expression of the associated downstream genes and operons (Fig. 6). The transcriptional effects in the three regions were strongest and most consistent between the PHB synthesis (phbA, phbB, phbC, and phbAB) mutants under nitrogen-limited conditions. It can be qualitatively observed that the genes showing the strongest effect size occurred directly downstream of a predicted motif. We then manually curated genes into four groups (A, genes with motif upstream; B, genes without a motif upstream; C, operons with a motif upstream; D, operons without a motif upstream) to perform a significance analysis of the effect of a motif on the expression of the downstream gene. The average log 2 -fold change in genes with a motif upstream was found to be significantly lower than that in all of the other groups (Fig. 7A). In addition, we observed a striking trend in operons with a motif upstream where the first gene was strongly downregulated but each gene further down the operon had a decreased effect size. To study this effect further, we then classified genes into five groups (A, first-order genes with a motif upstream; B, secondorder genes with a motif upstream of the operon; C, third-order genes with a motif upstream of the operon; E, operons with no motif upstream; Z, first-order genes without a motif upstream) for further significance analysis. The average log 2 -fold change in first-order genes with a motif upstream was found to be lower than that of all of the other groups (Fig. 7B). A further reduction of the effect size was seen among groups A, B, and C. This trend was observed only in operons in the three regions that had a motif upstream and not in operons with no motif (Fig. 7C). This interesting effect of gene expression intensity across operons is a known phenomenon and has been discussed to some degree in the literature (75). To determine the common location of the motif, the distance from the open reading frame (ORF) to the closest upstream motif was mapped (Fig. 7D). The motif was found to occur most often within the first 50 bp upstream of an ORF, consistent with a transcriptional regulatory function. When the locations of the motifs are compared to recently published TSS data on S. meliloti, there is close alignment of many of the motifs with mRNA TSS (mTSS) positions (Fig. S1D) (33). These analyses strongly support the hypothesis that the discovered motif represents the sites for transcriptional regulator binding that influence the expression of genes affected in the PHB cycle mutants. for transcripts in PHB synthesis mutants shows signal decay in operons with a predicted motif upstream (each dot represents a gene, and colored connected dots denote genes in the same operon). (D) First-order genes with a predicted motif upstream tend to be Ͻ200 bp away from the predicted motif.

(iii) Enrichment of Fnr homologs in the pSymA regions.
Recognizing that the Fnr-like motif is enriched in the pSymA regions of interest, we searched for potential Fnr-like regulators encoded in these regions. A Pfam search of the S. meliloti Rm1021 proteome for the main Fnr family domain (HTH_Crp2 [PF13545]) revealed 10 predictedgene-encoded proteins with this domain. The HTH_Crp2 domain is an annotation of a helix-turn-helix pattern found in the Crp protein of E. coli. Further searching by using the cNMP_binding (PF00027) domain identified another four predicted proteins with the cNMP_binding domain, including an HTH_Crp domain that had escaped annotation in the Pfam database. In total, there were 14 genes predicted to encode proteins in the Crp/Fnr family in S. meliloti Rm1021 (Table S2). Of these, four are in the 120-kb region of pSymA, which is a 15.9-fold overrepresentation over that expected on the basis of genome-wide occurrence. Additional Fnr-related family members reside in each of the two smaller regions, i.e., fixK2 in the 7-kb region and SMa0662 in the 37-kb region.
As discussed earlier, the 120-kb region includes the previously identified FixK1 and NnrR proteins but also contains SMa1141 and SMa1207, which are both putative Fnr/Crp regulatory proteins. Of these, FixK1 is a known regulator of nitrogen fixation in S. meliloti and NnrR is a regulator of the nap, nir, and nor operons (72,76). The fixK2 gene in the 7-kb region encodes another characterized Fnr regulator, but its specific role in the regulation of the fix operons in S. meliloti is not known. There is no information available regarding a regulatory role for SMa0662 in the 37-kb region.
While all six proteins were predicted to be structurally related to Fnr, they have one key difference in common. Unlike the originally identified Fnr protein in E. coli, they do not possess an N-terminal cysteine-rich region. This region allows Fnr to respond to oxygen (also NO) levels in the cell through an oxygen-labile Fe-S cluster in the N-terminal region. This cluster undergoes conversion upon oxygen binding, and the protein is unable to bind to DNA and act as a transcriptional activator (77). Since the six Fnr homologs in the pSymA regions are expected to be unable to detect cellular levels of oxygen, their function must be dependent on reaction to the presence of another compound or environmental condition or activation by another protein or be transcriptionally regulated. It has been previously established that fixK1 transcription is regulated by FixLJ (in the presence/absence of oxygen) and by FixT1 (69,78). NnrR induces a regulatory response in its associated genes in response to NO, but no mechanism has been elucidated (72).
Comparison and sequence alignment of the six Fnr/Crp proteins identified showed that five (SMa1207 excluded) contain a conserved, cyclic mononucleotide (cNMP) binding domain (Fig. S7). This is a trait shared with Crp, whose regulatory activity in E. coli is well known to be controlled by the binding of cyclic AMP (cAMP) to induce a catabolite repression response (79). The cNMP binding domain in the Fnr/Crp proteins suggests a regulatory role for cNMPs in the regulation of the pSymA regions through the involvement of cAMP or another cNMP, although this has not been studied in Fnr/Crp proteins in Rhizobiales. This would involve intracellular signaling to modulate the activity of these regulators as a reflection of the overall metabolic state of the cell. No experimental work has been published in support of this hypothesis. There is a known instance of catabolite repression in S. meliloti which is modulated through the incomplete phosphotransferase system to favor the use of succinate over many other carbon sources (80). Succinate-mediated catabolite repression is regulated separately from cAMP levels, but this does not discount the possibility of other regulatory systems with cAMP or other cyclic mononucleotides (81). Previous work with S. meliloti Rm2011 showed a decrease in the abundance of fixK1 when cAMP was overproduced by overexpressed cyaJ, indicating a plausible relationship between internal cAMP levels and fixK1 expression (32). No other genes in the three pSymA regions were affected in the cyaJ overexpression study, so the cAMP concentration in the wild type is not enough to trigger the transcriptional effect observed in these regions (SMa1082 excepted). Overall cAMP levels are reflective of the cell's metabolism and thus would be useful in modulating a regulatory system that responds to internal levels of carbon and/or nitrogen. A speculative model is that disruption of the PHB cycle interferes with cellular respiration and central metabolism, which may lead to alteration of the intracellular levels of cAMP, which could act as a ligand for an Fnr protein that regulates genes in the three pSymA regions identified.
Conclusions. We applied transcriptomic profiling to PHB cycle and glycogen pathway disruption mutants to identify and amplify the signal associated with the synthesis and degradation of PHB, as distinguished from individual gene mutation effects. Combined with computational genomic analysis, this revealed three distinct regions of pSymA that contain genes whose expression under nitrogen-limited conditions is impacted by PHB cycle disruption but not glycogen pathway disruption. These regions of pSymA contain genes involved in nitrogen uptake, nitrogen utilization, symbiotic nitrogen fixation, and denitrification, underlining the interrelation of PHB accumulation and nitrogen levels in the cell. We hypothesize that PHB cycle disruption results in an inability of the cell to respond appropriately to nitrogen limitation. Further examination of the genes whose expression was impacted in PHB cycle mutants led to the discovery of an associated motif that closely resembles the previously discovered FixK motif. The density of this motif was shown to correlate specifically with the genome regions with the greatest transcriptional change, and the presence of the motif upstream of genes in these regions correlates with the decreased gene expression level in the PHB synthesis mutants under nitrogen-limited conditions. These regions encode a high density of candidate Fnr/Crp family regulators that likely bind similar motifs. In the absence of any published findings that can explain the transcriptional effects seen, we speculate that intracellular cAMP levels in the PHB cycle mutants may be altered, leading to the modification of an Fnr/Crp family protein, which results in the downregulation of genes in three regions of pSymA through transcription factor binding of the motif discovered.
New questions that arise from this work include the nature of the mechanism of gene regulation in the regions identified; how gene regulation is influenced by PHB synthesis; whether intracellular cAMP levels are affected in the PHB cycle mutants; whether the local enrichment of proteins carrying the motif discovered is involved in fine-tuning of the regulatory response; and how this regulatory effect is connected to sensing of the carbon/nitrogen ratio. We hypothesize that one or more of the Fnr/Crp family proteins are required for some of the differential expression that is observed, and this can be tested through subsequent experiments.

MATERIALS AND METHODS
Bacterial strains, media, antibiotics, recombinant DNA, and bacterial genetics. The strains and plasmids used in this study are listed in Table 1. Bacteria were cultured under previously described conditions, in TY or LB with appropriate antibiotics (82). The antibiotics used with S. meliloti were neomycin (200 g/ml) and streptomycin (200 g/ml); those used with E. coli were kanamycin (25 g/ml) and ampicillin (100 g/ml). The deletion mutations of phbAB, phbC, and phaZ were carried out in a similar manner, by utilizing pK19mobsacB selection of homologous recombination on sucrose. For phbAB, phbC, and phaZ, 400-bp synthetic fragments that contained the 200 bp immediately up-and downstream of the corresponding ORF were used. These fragments were synthesized by BioBasic and delivered in purified pUC57. The inserts were subcloned from their vectors into pK19mobsacB as EcoRI/XbaI fragments, and the resulting constructs were transferred to Rm1021 by triparental conjugation. Neomycin selection resulted in integration of the plasmid into the chromosome via single crossover. Selection with 5% sucrose selected for resolution, resulting in revertants or double-crossover deletion mutants. The colonies were screened for neomycin sensitivity to confirm the loss of the plasmid. Confirmed clones were screened by colony PCR to differentiate between revertants and deletion mutants. The resulting strains were SmUW234 (ΔphbAB), SmUW235 (ΔphbC), and SmUW236 (ΔphaZ).
The RNA-seq analysis reported here confirmed that no transcript (deletion mutant) or no complete transcript (insertion mutant) corresponding to the deleted gene(s) of each of the mutants was detectable.
In preparation for RNA extraction, 5-ml cultures were grown to saturation in tryptone-yeast extract broth. The cells were washed twice and resuspended in an equivalent volume of 0.85% (wt/vol) NaCl. Aliquots (750 l) of washed cells were used to inoculate 250-ml baffled shake flask cultures containing 75 ml of broth. Nitrogen-limited medium consisted of 5 g of mannitol, 0.5 g of yeast extract, 1 g of K 2 HPO 4 , 0.2 g of MgSO 4 · 7H 2 O, 0.1 g of NaCl, 0.5 g of CaCl 2 · 2H 2 O, and 0.004 g of FeCl 3 · 6H 2 O/liter, adjusted to pH 7.2. Balanced medium was the same but had 20 times as much yeast extract (10 g).

RNA extraction.
Antibiotics were not used in cultures from which RNA was extracted, to avoid transcriptomic bias introduced by antibiotic presence. Cultures were shaken at 200 rpm at 30°C. Cells were harvested at mid-to late-log phase as indicated by optical density at 600 nm readings of 0.5 Ϯ 0.05 (Fig. S1). To harvest cells, 45 ml of culture broth was transferred to a sterile container and 5 ml of stop solution (5% phenol in absolute ethanol) was added. The cells were collected by centrifugation at 6,000 ϫ g and 4°C for 10 min, and the liquid was decanted. Pellets were frozen in liquid nitrogen and stored at Ϫ80°C for up to 2 months prior to RNA isolation.
RNA was extracted by the hot phenol method, much as previously described (83). Briefly, cells were lysed and total nucleic acids were initially quantified with a NanoDrop. Equivalent amounts of total nucleic acids were then compared with DNA standards of known quantity to determine the amount of contaminating genomic DNA. The samples were then treated with 1 U of RNase-free DNase I/g of DNA and 40 U of RiboLock/g of RNA. The RNA was cleaned up with the RNeasy minikit (Qiagen), and RNA quality was determined by electrophoresis on an agarose-formaldehyde denaturing gel stained with ethidium bromide. The RNA concentration was determined with a NanoDrop 1000 (Thermo Fisher), and subsamples from isolates with Ͼ1 g/l were run on agarose gels. If the 16S rRNA band was approximately half as bright as the 23S rRNA band, the samples were used for RNA-seq analysis.
RNA preparation for sequencing. rRNA depletion was performed with 1 g of total RNA and the Meta-Bacteria Ribo-Zero rRNA Removal kit (Mandel Scientific, Guelph, Ontario, Canada) at the Genome Quebec Innovation Center of McGill University. Final RNA was isolated with the RiboMinus Concentration Module and eluted with 17 l of Fragment, Prime, Finish buffer (Life Technologies, Inc.). The eluted RNA was then fragmented and prepared for cDNA synthesis. The TruSeq Stranded mRNA Sample Prep kit was used for cDNA synthesis, and then the quality of the library and fragment size were assessed with a LabChip, a LightCycler 480 II, and an Infinite M200 Fluorimeter (Roche, Mississauga, Ontario, Canada, and Tecan, Mannedorf, Switzerland).
Illumina sequencing and raw-read processing. The cDNA was sequenced (2 ϫ 100 bp) on an Illumina HiSeq 2000 instrument at the Genome Quebec Innovation Center of McGill University. There was an average of 21,592,337 reads with an average Phred quality score of 35 for each sample. The average number of base pairs per sample was 4,318,467,357. Trimmomatic was used to strip the barcode from the raw reads by using the default settings (84). TopHat was used to assemble reads to the reference genome; reads with more than two mismatches, indels, or gaps were not used (85). Cuffdiff was used to compare the normalized transcript count data (fragments per kilobase per million reads [FPKM] transcript  (86). Ten sequences aligned with each gene was required as the minimum threshold to allow testing; otherwise, no test was performed (significance displays as "no"). Only transcripts that were significantly different in the mutant and parental strains were considered in the following analysis. For the differential-expression data obtained, see http://www.bitbucket.org/ mdelow/phbrnaseq/. Motif discovery and comparison. The 200 bp upstream of every gene in the clusters delineated on the heat map were extracted, along with the 200 bp upstream of every gene in the S. meliloti Rm1021 genome, with RSAT (87). These files were used to perform de novo motif discovery with MEME (88). The top resulting motifs were mapped to the pSymA regions with FIMO (89). The motif outputs from MEME under a P value threshold of 1 ϫ 10 Ϫ15 were input into Tomtom for comparison to the prokaryotic motif databases (90). The Tomtom output results under a P value threshold of 0.001 were considered and are reported here.
Statistical analysis. Analysis of variance was performed by using the average log 2 -fold values for each of the manually curated groups. This was done in R with the stat.aov, TukeyHSD, and HSD.test functions.
Figure generation. Figures 2 to 7 were generated in part through data visualization in R. The scripts and data files used in R for all data visualization are available as discussed in the paragraph on data availability below.
Figures 2 and 5 used the starting coordinate of each ORF to anchor the log 2 -fold change (y-axis value) along the x axis based on genomic location.
For Fig. 4, S2, and S8, hierarchical clustering of log 2 expression values was performed with the hclust() function within R based on "complete linkage" clustering based on euclidean distances. To identify expression clusters, we used a semiautomated approach in which expression clusters were visually identified and fixed thresholds were selected to define cluster membership. At the highest level, the gene expression dendrogram split into two large clusters, A and B (see Fig. S8). Cluster A was identified as a biologically relevant cluster and was thoroughly investigated in this work. Expression cluster B was larger but had weaker expression patterns and subtle differences between subclusters, and we were unable to identify enriched motifs at the same level of significance as in cluster A. Therefore, it was further subdivided into subclusters with a lower threshold of 7, and these subclusters were analyzed further. Lastly, cluster 4 was also subdivided into two subclusters with noticeable differences in expression. Figure 5 visualized motif density along the genome. First, the FIMO tool was used to detect all instances of the discovered motif genome-wide. The coordinates of each motif location were then used to visualize motif density along the genome by using R's density() function with a scaling bandwidth parameter (b.w.) of 30,000. Figure 6 was generated in Geneious by using the Rm1021 annotation to visually map the locations of genes in the 120-kb region with the locations matching the "discovered" motif overlaid. Further editing was done in Inkscape to change the color of each gene arrow to correspond to the log 2 -fold change. A basic heat map of the average log 2 -fold change in the expression of genes in the 120-kb region in the PHB synthesis mutants under nitrogen-limited conditions was produced in R, and the colors of the tiles were mapped onto the gene arrows in the Geneious output. No changes in gene or motif locations were made.
The statistical analysis used for Fig. 7A and B is described above. Genes were manually curated into groups based on the direct presence of a motif upstream of the ORF (the motif was assumed to be functional in either direction because of its palindromic nature) without any other intervening genetic elements. Genes with multiple motifs upstream were not treated any differently. Motifs that occurred within genes were ignored. For Fig. 7A, operons were treated as a single group, with the presence of a motif being tested upstream of the first gene in the operon. For Fig. 7B, the first gene in an operon was separated into another group from the downstream genes in the operon, but the classification of these groups was still dependent on the presence of a motif upstream of the first gene in the operon. For Fig. 7D, we used a kernel density plot to visualize the distance between upstream motifs and the closely associated ORF. This was done by using the R plot(density) function with a line width of 2. Figure S1D and Fig. 6 were produced in Geneious. Figure S1D used the same reference as Fig. 6, with the TSS data layered on top (33). Figure S7 was produced by using the align utility on the protein sequences of the associated genes of interest.
Data availability. The raw RNA sequence data obtained in this study are available through the Sequence Read Archive (SRA) through BioProject accession number PRJNA369802 and associated SRA and BioSample accession numbers, titled "Analysis of the transcriptional response to mutations in PHB cycle genes in Sinorhizobium meliloti reveals disproportionately affected regions on pSymA containing nitrogen utilization genes." The BioSample accession numbers (36) with the analysis data and R script are available at https://bitbucket.org/mdelow/phbrnaseq.

ACKNOWLEDGMENTS
This work was financially supported by a Strategic Projects Grant and a Discovery Grant from the Natural Sciences and Engineering Research Council of Canada and by Genome Canada through the Applied Genomics Research in Bioproducts or Crops (ABC) program for the grant titled "Microbial Genomics for Biofuels and Co-Products from Biorefining Processes" awarded to David Levin and Richard Sparling. Maya D'Alessio was the recipient of an Ontario Graduate Scholarship.
We are grateful to Justin Zhang for assistance with initial RNA-seq data processing and Michael Mansfield for providing feedback on the manuscript.