Genome-Wide Analysis of ResD, NsrR,
ABSTRACT Upon oxygen limitation , the Bacillus subtilis ResE sensor kinase and its cognate ResD response regulator play primary roles in the transcriptional activation of genes functioning in anaerobic respiration . 
The nitric oxide ( NO ) - sensitive NsrR repressor controls transcription to support nitrate respiration . 
In addition , the ferric uptake repressor ( Fur ) can modulate transcription under anaerobic conditions . 
However , whether these controls are direct or indirect has been investigated only in a gene-specific manner . 
To gain a genomic view of anaerobic gene regulation , we determined the genome-wide in vivo DNA binding of ResD , NsrR , and Fur transcription factors ( TFs ) using in situ DNase I footprinting combined with chromatin affinity precipitation sequencing ( ChAP-seq ; genome footprinting by high-throughput sequencing [ GeF-seq ] ) . 
A significant number of sites were targets of ResD and NsrR , and a majority of them were also bound by Fur . 
The binding of multiple TFs to overlapping targets affected each individual TF 's binding , which led to combinatorial transcriptional control . 
ResD bound to both the promoters and the coding regions of genes under its positive control . 
Other genes showing enrichment of ResD at only the promoter regions are targets of direct ResD-dependent repression or antirepression . 
The results support previous findings of ResD as an RNA polymerase ( RNAP ) - binding protein and indicated that ResD can associate with the transcription elongation complex . 
The data set allowed us to reexamine consensus sequence motifs of Fur , ResD , and NsrR and uncovered evidence that multiple TGW ( where W is A or T ) sequences surrounded by an A - and T-rich sequence are often found at sites where all three TFs competitively bind . 
IMPORTANCE Bacteria encounter oxygen fluctuation in their natural environment as well as in host organisms . 
Hence , understanding how bacteria respond to oxygen limitation will impact environmental and human health . 
ResD , NsrR , and Fur control transcription under anaerobic conditions . 
This work using in situ DNase I footprinting uncovered the genome-wide binding profile of the three transcription factors ( TFs ) . 
Binding of the TFs is often competitive or cooperative depending on the promoters and the presence of other TFs , indicating that transcriptional regulation by multiple TFs is much more complex than we originally thought . 
The results from this study provide a more complete picture of anaerobic gene regulation governed by ResD , NsrR , and Fur and contribute to our further understanding of anaerobic physiology . 
Bacteria live in environments where many chemical and physical parameters constantly change . 
In addition , their life is affected by the presence and activity of other living organisms . 
These effects are sometimes detrimental to bacteria but are often manageable by orchestrating global gene expression in response to each environmental change . 
One such change they encounter in diverse environments and in living hosts is a fluctuation of oxygen concentration . 
Bacillus subtilis adapts to oxygen limitation by undergoing nitrate respiration or anaerobic fermentation that generates ATP by substrate-level phosphorylation . 
Gene regulation required for the adaptation from aerobic to anaerobic conditions is mainly controlled at the transcriptional level ( 1 ) . 
The signal transduction system composed of the ResD response regulator and the ResE histidine sensor kinase plays essential roles in nitrate respiration ( 2 ) and aerobic respiration ( 3 ) . 
Due to the dual roles of ResDE in aerobic and anaerobic respiration , B. subtilis is endowed with an additional layer of regulation for controlling genes that function in nitrate respiration when cells encounter oxygen-limited conditions in the presence of nitrate . 
This regulatory mechanism is executed by two [ 4Fe-4S ] - containing transcription factors ( TFs ) , namely , oxygen-sensitive Fnr and nitric-oxide ( NO ) - sensitive NsrR . 
Fnr plays a major role in autoregulation through the narK-fnr operon promoter ( 2 ) , as well as in the activation of the narGHJI ( respiratory nitrate reductase ) operon ( 4 , 5 ) . 
The NsrR repressor functions in upregulating ResD-dependent transcription of the nasDEF nitrite reductase operon in response to NO ( 6 ) . 
NO is generated from nitrite as a by-product during nitrate respiration ( 6 ) ; hence , B. subtilis likely senses nitrate availability through the presence of NO . 
In addition , nitrite reductase converts nitrite to ammonium , thus mitigating NO accumulation that is harmful to cells . 
NO detoxification could also be carried out under microaerobic conditions by flavohemoglobin ( 7 ) , which is encoded in B. subtilis by another ResD/NsrR-controlled gene , hmp ( 8 ) . 
We have shown in vivo that nasD and hmp transcription are highly repressed by NsrR during anaerobic fermentative growth ( in the absence of NO ) , but that the expression of these genes is induced under conditions that support nitrate respiration ( in the presence of NO ) ( 6 ) . 
Studies using electrophoretic mobility shift assays ( EMSAs ) and in vitro transcription confirmed that NsrR is an NO-sensitive repressor ( 9 , 10 ) . 
Electron paramagnetic resonance and resonance Raman spectroscopies demonstrated that NO interaction with iron in the [ 4Fe-4S ] cluster causes dinitrosylation of iron , leading to the dissociation of NsrR from the nasD promoter DNA ( 11 ) . 
The direct modification of the [ 4Fe-4S ] cluster by NO has also been reported with an NsrR ortholog ( 12 ) . 
To gain insight into the genome-wide function of NsrR in transcription , a microarray analysis was conducted ( 9 ) . 
The results showed that hmp and nasD are the genes most highly regulated by NsrR . 
In addition , we identified other genes moderately controlled by NsrR , some of which belong to the Fur regulon involved in iron homeostasis ( 13 ) . 
Although in vivo transcription assays confirmed the negative effect of NO on NsrR repression of the newly identified genes , an EMSA showed that NsrR only weakly binds to this class of promoter DNA in an NO-insensitive manner ( 9 ) . 
Hence , NsrR might indirectly control the transcription of these genes . 
Alternatively , efficient NsrR binding could require ternary DNA structure and/or other TFs to assist NsrR binding to DNA , both of which were lacking in the in vitro studies ( 9 ) . 
The following in vivo protein-DNA binding approach using chromatin affinity precipitation ( ChAP ) - chip revealed that NsrR interacts with some but not all of the genes identified by the preceding microarray results ( 14 ) . 
Some of the DNA sites directly targeted by NsrR were also bound by ResD and/or Fur ( 14 ) . 
However , the resolution of ChAP-chip is not high enough to precisely localize the exact binding sites of each TF . 
In this study , we revisited the genome-wide interaction of the three TFs -- NsrR , ResD , and Fur -- in cells cultured under anaerobic fermentative conditions to resolve each TF-binding sequence by adapting the high-resolution mapping method called genome footprinting by high-throughput sequencing ( GeF-seq ) . 
GeF-seq involves in situ DNase I digestion of genomic DNA followed by ChAP-seq , which was successfully used by members of our group to delimit AbrB-binding sites ( 15 ) . 
A comparison of binding sites targeted by each TF in the wild type and in mutant strains lacking the other two TFs has revealed that the binding of ResD , NsrR , and Fur could be cooperative or competitive depending on the targeted promoters and the other TFs binding to nearby or to overlapping sites . 
By using a single-nucleotide-resolution analysis of GeF-seq data , we were able to determine that all three TFs competitively bind within a 40-bp region of promoter DNAs that contain sequences similar to the consensus sequence of each TF . 
We also showed that DNA-bound ResD distributions differ depending on ResD function in the transcriptional control of each ResD-targeted gene . 
RESULTS
Rationale . 
Our previous ChAP-quantitative PCR ( qPCR ) showed that ResD , NsrR , and/or Fur often affects the binding of the other TF ( s ) ( 14 ) . 
To gain further understanding of how binding by the multiple TFs to promoter DNA affects individual TF-DNA interactions , we carried out GeF-seq both in the wild type and in mutants that do not produce the other two TFs . 
A comparison of TF binding between the wild type and double mutant strains is an advantageous approach . 
First , the binding of multiple TFs to nearby but not overlapping sequences might lead to cross-linking between His12-tagged and untagged TFs , which could contribute to protection of the DNA region against DNase I cleavage . 
Thus , affinity purification of the DNA during ChAP generates a longer DNA that includes more than a single TF-binding site . 
In fact , a search for NsrR-and ResD-binding motifs by MEME ( multiple expectation maximization for motif elicitation ) analysis of GeF-seq data from wild-type extracts often identified the Fur-binding motif ( data not shown ) . 
We could partially circumvent the problem by examining each TF-binding target in mutant strains . 
Second , each TF affects , positively or negatively , DNA binding by other TFs . 
Thus , a comparison of binding distributions of each TF in the absence or presence of other TFs likely provides us with more valuable information , such as the existence of cooperative and competitive binding , which is useful to understand the biological importance of combinatorial control exerted by multiple TFs . 
Therefore , we carried out GeF-seq using both the wild-type and double mutant strains . 
More specifically , NsrR , ResD , or Fur binding was determined in the resD fur , nsrR fur , or resD nsrR mutant cells , respectively . 
All GeF-seq experiments were carried out in cells cultured anaerobically in 2 yeast extract-tryptone medium ( YT ) supplemented with 0.5 % glucose and 0.5 % pyruvate to support anaerobic fermentation where NsrR is active as a repressor ( 6 ) . 
ResD and NsrR binding is localized at the hmp and nasD promoters . 
Previous ChAP-qPCR revealed that binding of ResD to nasD increases in the nsrR mutant background ( 14 ) , suggesting that NsrR efficiently competes with ResD . 
However , we were unable to distinguish whether this competition is due to either an overlapping or nearby sequence targeted by the two TFs . 
The GeF-seq results of ResD and NsrR binding to hmp ( Fig. 1A ) and nasD ( Fig. 1B ) at a single nucleotide resolution not only solved the question but also provided a new finding that the DNA-binding pattern of ResD correlates with its function in transcriptional control . 
Protein-binding regions shown in the left panels are expanded and shown in the right panels in Fig. 1 . 
NsrR interaction with the hmp promoter in the resD fur mutant was at an intensity comparative to that in the wild-type strain , but the binding peak in the mutant became more prominent within the region between 9 to 24 ( relative to the transcription start site ) ( Fig. 1A ) , which corresponds to the NsrR consensus sequence ( NsrR-1 ) predicted by bioinformatics analysis ( 16 ) . 
Given that Fur does not bind to hmp in either the wild type or the resD nsrR mutant , NsrR , not Fur , is likely the TF that outcompetes ResD for binding to the hmp promoter under anaerobic fermentation conditions . 
In addition , the GeF-seq data from the mutant revealed that ResD interacts with two distinct regions upstream of the core promoter and the transcription start site , whereas very low association of ResD with the upstream region was detected in the wild-type cells . 
The upstream ResD-binding site ( 93 to 35 ) corresponds well with the sites previously determined by in vitro DNase I footprinting ( 76 to 40 ) ( 17 ) and hydroxyl radical footprinting ( 81 to 47 ) ( 18 ) . 
The downstream ResD-binding site , which covers the 10 sequence and the transcription start site , overlaps with the NsrR-binding site . 
The downstream ResD site has not been detected in the in vitro binding studies ( 17 , 18 ) . 
We interpret this result as meaning that the upstream ResD-binding site is required for the activation of hmp transcription as demonstrated in the earlier study ( 17 ) , and the downstream site was not detected in the in vitro footprinting experiments because ResD binding to this site requires another protein , very likely RNA polymerase ( RNAP ) , as described later . 
Both NsrR - and ResD-binding profiles to nasD showed similarities to as well as differences from those to the hmp promoter region . 
Our previous work showed that the nasD promoter carries a sequence similar to the predicted NsrR consensus sequence in hmp ( 6 ) . 
This nasD sequence functions as the operator for NsrR based on the effect of base substitutions on in vitro NsrR-DNA binding and in vivo transcription ( 9 ) . 
Binding of NsrR to nasD , as detected by GeF-seq , covered the previously identified NsrR-binding site ( Fig. 1B ) . 
However , the binding intensity was lower in the resD fur mutant unlike that for hmp . 
The requirement of ResD for efficient binding of NsrR to nasD was also demonstrated by ChAP-qPCR ( 14 ) . 
As observed with hmp , ResD was more enriched at the nasD promoter in the absence of NsrR . 
At the hmp promoter , the upstream site showed greater interaction with ResD than the downstream site , whereas for nasD , the promoter region ( corresponding to the downstream site of hmp ) had a higher intensity of ResD binding , and upstream of 35 ( corresponding to the upstream site of hmp ) in the mutant was recognized , at best , only as a shoulder of the peak . 
In summary , the results clearly demonstrated a high-resolution view of NsrR interference in ResD interaction with the hmp and nasD promoter DNA . 
Previous studies suggested that the transcription of hmp and nasD is controlled by ResD and NsrR through similar mechanisms , but this study revealed that there are differences between each promoter with respect to ResD and NsrR . 
Two different ResD-binding profiles reflect the regulatory roles of ResD . 
ResD has been known to act as a direct transcriptional activator but has not been identified as a direct transcriptional repressor . 
GeF-seq revealed two different ResD-binding profiles ( Fig. 2 ) , which suggested the possibility that ResD likely plays diverse roles in transcriptional control . 
Table 1 lists representative promoters where ResD strongly binds , as well as how ResD regulates the transcription of each gene based on previously published and unpublished results together with those obtained from this study . 
Seventeen genes ( 18 sites ) showed ResD binding at both their promoters and throughout the entire coding regions ( Table 1 , group A ) . 
The binding of ResD was clearly sharper and stronger to the promoters than to the coding regions ( Fig. 2A to D ) . 
As evident with the nasDEF and the cydABCD operons , the trail of binding ends at the - independent terminator of the operon , suggesting that a portion of the ResD population that binds to the promoter remains associated with the RNAP elongation complex . 
ResD-binding regions in the second class of genes are limited to the promoter DNA ( Table 1 , group B , and Fig. 2E to G ) . 
We examined whether ResD plays distinct roles in transcription between genes that show different binding profiles . 
As shown in Table 1 , most of genes in group A , to which ResD binds both the promoter and coding regions , function in aerobic or anaerobic respiration and are controlled positively by ResD . 
These genes include hmp ( 8 ) , nasD ( 19 ) , narG ( 4 ) , ctaA ( heme A synthase gene ) ( 20 ) , ctaB ( major heme O synthase gene ) ( 21 ) , qox ( cytochrome aa3 quinol oxidase gene ) ( 22 ) , and cydA ( cytochrome bd ubiquinol oxidase gene ) ( 23 , 24 ) . 
Therefore , we speculated that ResD activates the transcription of all genes in this class . 
The speculation was tested using transcriptional lacZ fusions to as-yet-uncharacterized promoters that belong to this class ( Fig. 3 ; see also Fig . 
S1 in the supplemental material ) . 
Transcriptional control of these genes by ResD has not been reported , and the biological functions of most of the genes remain to be uncovered . 
As expected , ResD indeed activated yxiE and yozB transcription ( Fig. 3A and B ) , as well as ctaO ( minor heme O synthase gene ) ( 25 ) , ydbL , and yfmQ ( Fig . 
S1 ) . 
A previous study using affinity purification of RNAP followed by mass spectrometry identified ResD as an RNAP-associated protein ( 26 ) . 
The high affinity of ResD to RNAP might explain why ResD travels with elongating RNAP through coding regions . 
Consistent with this finding , a previous EMSA showed that ResD binding to the nasD promoter is enhanced in the presence of RNAP ( 10 ) . 
Altogether with these results , we concluded that genes exhibiting ResD associations with both promoter and coding regions are directly activated by ResD . 
Next , we chose five genes -- glpF , yjlC , ywcE , ymfC , and ytcP -- representing the second class of ResD-bound loci to examine whether and how ResD controls the transcription of these genes ( Fig. 3C to F ) . 
glpF and the downstream glpK constitute an operon encoding glycerol uptake facilitator and glycerol kinase , respectively . 
The third gene in the operon , glpD ( glycerol-3-phosphate dehydrogenase gene ) , is transcribed from a gene-specific promoter as well as the glpF promoter ( 27 ) . 
A previous transcriptomic analysis showed that glpFK transcription is repressed under anaerobic conditions compared with that under aerobic conditions and that aerobic glpD expression is upregulated in the resDE mutant ( 1 ) . 
These results suggested that ResD negatively controls the transcription of the glpF operon , which was confirmed by data showing that glpF transcription is derepressed in the resD mutant ( Fig. 3C ) . 
Taken together with the GeF-seq results , we concluded that ResD is a direct repressor of the glpF operon . 
Similarly , the transcription of the yjlC-ndh operon was downregulated under anaerobic conditions and increased in the resDE mutant under aerobic conditions ( 1 ) . 
The regulatory loop of this operon was previously reported ( 28 ) . 
A higher NAD concentration leads to repression of the yjlC-ndh operon by the Rex repressor that responds to a decrease in the NADH/NAD ratio ( 29 ) , thereby leading to a reduction in NADH dehydrogenase encoded by ndh and a consequent elevation in the NADH/NAD ratio ( 28 ) . 
The operon has a relatively long untranslated leader sequence , and the transcription initiation site resides 257-bp upstream of the initiation codon . 
The Rex-binding site was identified within the leader region by deletion analysis and EMSA ( 28 ) . 
Our GeF-seq results showed that ResD binds to a sequence that includes the 10 site of the yjlC promoter ( 27 ) ( see Data Set S1 ) . 
ResD was able to repress the expression of lacZ fused to the yjlC promoter ( 138 to 29 ) that lacks the Rex-binding site ( 145 to 166 ) , thus functioning in a Rex-independent manner ( Fig. 3D ) . 
ywcE , whose transcription is under negative control by AbrB , encodes a holin-like protein , and the ywcE mutant forms spores with a reduced outer coat ( 30 ) . 
ResD repressed ywcE ( Fig. 3E ) at least under anaerobic conditions . 
ResD also repressed the transcription of ymfC ( GntR family TF ) ( Fig. 3F ) , but the effect of the resD mutation on ymfC transcription was not as strong as those of other genes listed in Fig. 3 . 
ytcP expression was not detected under the culture conditions tested regardless of the presence or absence of ResD . 
We previously showed that ResD functions as an antirepressor of Fur for ykuN ( flavodoxin gene ) that also belongs to this class ( 14 ) . 
In summary , ResD plays the role of a repressor ( or antirepressor ) in transcriptional control of the second class of genes , but not as a direct activator , by binding to promoter regions . 
To the best of our knowledge , this is the first report to show that ResD directly represses transcription . 
Some ResD-binding promoters interact with NsrR and/or Fur . 
Previous studies including ChAP-chip analyses showed that ResD and NsrR , by directly binding to hmp and nasD , positively and negatively regulate the transcription of these genes . 
Furthermore , ResD and NsrR modulate Fur-dependent repression of ykuN , which involves the direct binding of three TFs . 
The questions that remained to be answered were whether there are more promoter regions targeted by a combination of the three TFs and , if there are , whether and how the coordinated or simultaneous interactions affected the transcription of these genes . 
Table 1 classifies ResD-binding genes by the ResD-binding profiles described above and by their interactions with NsrR and/or Fur . 
NsrR does not bind to the group A genes listed in Table 1 either in the wild type or in resD fur mutant strains , except for hmp and nasD as described above . 
This unique features of hmp and nasD are attributed to the physiological importance of NO-sensitive NsrR repression , the mechanism by which the RNAP-ResD complex efficiently competes for promoter binding with NsrR only when NO is present . 
Table 1 , group B , shows that ResD-binding patterns to promoter DNA somewhat correspond to whether NsrR also binds to the DNA . 
Genes in group Ba to which ResD binds in both the wild-type and mutant strains generally do not interact with NsrR or Fur , except yvaF , which exhibits higher binding by NsrR in the absence of ResD and Fur . 
ykuN and fbpC , which belong to group Bb , are enriched for both ResD and NsrR in the wild-type strain but not in the mutants . 
By contrast , Fur binding is less significantly affected by the nsrR and resD mutations . 
Previous studies showed that fbpAB and fbpC are Fur-regulated RNA chaperones for the fsrA small RNA ( sRNA ) that functions in the iron-sparing response ( 31 ) . 
In comparison to fbpAB where NsrR and ResD bind weakly at most ( data not shown ) , all three TFs clearly bind to the fbpC promoter ( Table 1 ) . 
FsrA-dependent control of resA expression was uncovered in a previous study ( 31 ) . 
Given that resD is under the control of the resA operon promoter that is indirectly activated by ResD itself ( 32 ) , it is attractive to envision the involvement of FbpC in the ResDE autoregulatory pathway ; however , the role of the RNA chaperones , particularly FbpC , has not been fully established . 
Genes in the last group , Bc , could interact with NsrR and/or Fur ( except yppF ) , but the binding intensity of each TF is dramatically decreased in the wild-type cells , indicating the competitive nature of binding . 
An untranslated region ( UTR ) named S842 was identified upstream of yppF ( 27 ) , and we found that the yppF ResD-binding site is within the 146-bp region coding for S842 . 
S842 is induced by anaerobiosis ( 27 ) , thus ResD binding of the S842 site might be classified as belonging to group A instead of Bc , although it is difficult to distinguish the ResD-binding profile because of the small size of S842 . 
Transcriptional control of group Bc genes such as ymfC by each TF is often subtle and complex ( Fig. 3 and data not shown ) . 
This is likely caused by the overlap of each TF-binding site . 
For example , NsrR and/or Fur likely represses transcription by occupying the site available in the resD mutant . 
In summary , GeF-seq and transcription assays enabled us to identify genes that are direct targets of ResD , which could be further classified into different groups using in vivo binding profiles , ResD 's role in transcription , and whether or how NsrR/Fur affects mutual DNA binding . 
The three TFs show combinatorial binding to DNA . 
Genome-wide association studies of the TFs showed that ResD , NsrR , and Fur bind to a significant number of overlapping sites ( Data Set S1 ) . 
These sites not only are in promoter regions but also localize outside promoter DNA . 
Each TF binds to a larger number of the overlapping sites when the other TFs are missing compared with that in the wild-type strain , suggesting that these TFs participate in competitive interactions more than cooperative interactions at overlapping binding sites ( Table 2 ) . 
Fur shares target sites with ResD and/or NsrR but almost exclusively binds to those associated with both ResD and NsrR . 
Such a trend was observed in both the wild type and the mutant backgrounds . 
These TFs use overlapping binding sequences of target genes ( Data Set S1 ) . 
The combinatorial binding by the TFs is classified into three groups based on how the binding of the TFs is mutually affected ( Table 2 ) . 
The first group of genes includes ykuN and fbpC , which were classified into group Bb in Table 1 . 
The TFs cooperatively bind to either the promoter DNA ( ykuN , fbpC , and exlX [ extracellular endoglucanase gene ] ) or the coding region of ppsB ( plipastatin synthase gene ) . 
The binding of ResD and NsrR to these genes requires Fur , which is in good agreement with ChAP-qPCR results of ykuN ( 14 ) , whereas Fur is still able to bind to these genes in the absence of ResD and NsrR , albeit with less intensity ( fbpC and ppsB ) or with altered binding distribution ( ykuN , which will be discussed later and in Fig. 5 ) . 
The second group of genes includes those classified into group Bc in Table 1 . 
Each TF only binds to these genes , or the binding intensities increase when the other two TFs are missing . 
This competitive binding of the three TFs was also observed at two other sites , the promoter region of yngD ( nrnB ) that encodes nano-RNase B and the ydjA ( BsuM DNA restriction system [ 33 ] ) coding region . 
The three TFs also competitively bind to yet-unassigned genes located between yoaM-yozS and ynfE-xynC , as described later . 
As Fur binding of pps is enhanced in the mutant , although the intensity value remains lower than the defined threshold for peak detection , pps likely belongs to this group . 
This group of genes generally shows higher intensities of NsrR binding compared with those of ResD and Fur . 
The third group shows binding competition that is different from the second group . 
NsrR binding in the wild-type background is weak or not detected but enhanced in the resD fur double mutant . 
Fur binding also increases in the resD nsrR mutant despite the observation that intensity values do not reach the threshold in the case of yvaF , spoIIIC , and spoIVCB . 
As ResD binding does not significantly change between the wild type and the mutant , ResD likely has a higher binding affinity than NsrR and Fur . 
NsrR and Fur likely compete with each other for binding in the absence of ResD . 
In summary , GeF-seq in the wild type and in double mutants has made it possible for us to classify all three TF binding characteristics into three groups based on the effect of mutual DNA interactions . 
Search for consensus sequences for Fur , ResD , and NsrR binding . 
As considerable numbers of sites were bound by ResD/NsrR/Fur , and the TFs often bind to overlapping sites with different affinities depending on the site , it is not easy to identify each TF-binding motif , particularly for ResD and NsrR . 
A search for the consensus sequences was conducted using data collected from GeF-seq in the mutant backgrounds . 
Results using the MEME motif search ( 34 ) and BiPad ( a two-block de novo DNA motif-finding tool ) ( 35 ) are shown in Data Sets S2 to S4 . 
Fur . 
Previous studies by Helmann 's group determined the consensus sequence of Fur-binding sites ( 36 , 37 ) . 
The binding sequence was identified as a 7-1-7 inverted repeat ( TGATAATNATTATCA , where N is any nucleotide ) where one Fur dimer interacts . 
The sequence was further confirmed in a study of condition-dependent transcriptome analysis in which Fur operators were characterized using a MEME motif search ( 27 ) . 
Here , we carried out a MEME search using all Fur-binding regions identified by GeF-seq and found a common TGANAA motif ( Data Set S2 ) . 
As MEME did not detect the 7-1-7 motif in all Fur-binding sites , we repeated the search using subsets of genes divided into four groups by binding strength ( Data Set S2 ) . 
The highest binding sites ( intensity rank 1 to 25 ) were found to harbor a 7-1-7 motif similar to the previously identified sequence , except that the fourth T in the upstream 7 was only moderately conserved ( 17 of 25 ) and the fourth A of the downstream 7 was the least conserved ( 8 of 25 ) in the newly identified motif . 
Thus , the common motif in genes of the highest intensity rank contains a modified 7-1-7 ( TGANAATNATTNTCA ) . 
When Fur-binding sites that belong to the next three ranks ( 26 to 50 , 51 to 75 , and 76 to 98 ) were used , MEME identified shorter consensus motifs that overlap the upstream half-site in the 7-1-7 motif ( TGANAA ) ( Data Set S2 ) . 
The result might suggest that either Fur binds as a monomer or the downstream half motif was not detected , because MEME motif search does not identify bipartite motifs separated by gaps of various length . 
Therefore , we used BiPad , which is capable of predicting various pairs of bipartite motifs with different gap lengths . 
Changes of spacing between the halfsites from 0 to 2 when incorporated into the search successfully detected target sequences in all potential Fur-binding sites ( Fig. 4 ; see also Data Set S2 ) . 
Fifty-five Fur-binding sites contain the modified 7-1-7 motif identified by MEME . 
An additional 27 sites and 16 sites were detected to have a common 7-2-7 and 7-0-7 sequence motif ( Fig. 4 ) , respectively , although the downstream half of the 7-0-7 configuration is not well conserved . 
The high-intensity binding sites mostly carry the original 7-1-7 motif as expected , and often carry two overlapping 7-1-7 motifs to generate a 21-bp binding site , a sequence similar to the 19-bp sequence proposed earlier ( 36 ) . 
We also noticed tandem TGANAA-like sequences in both orientations near the 7-1-7 motif , suggesting cooperative binding to multiple sites to compensate for weaker affinities . 
This result might indicate that a Fur dimer interacts with each half-site of the Fur operator in a flexible manner . 
ResD . 
The ResD consensus sequence was determined using all the ResD-binding sites identified in the nsrR fur mutant by GeF-seq . 
A MEME motif search identified a common 10-1-10 motif . 
TGANW6 or TNTGANW4 ( where W is A or T ) was observed particularly in high-intensity ResD-binding sites as tandem repeats with 1-bp spacing ( Fig. 4 ; see also Data Set S3 ) . 
A directly repeated TGANAANW3 motif was also visually detected . 
The alignment of these logos is very complex , and the overall A T richness of the sequence where ResD binds makes it difficult to conclude which motif fits best with an individual binding site . 
However , it is likely that a 10-bp motif with 1-bp spacing functions primarily as the ResD target sequence . 
The data support the hydroxyl radical footprinting results suggesting that ResD likely binds tandemly at the same surface of the DNA helix in the hmp and nasD promoters ( 18 ) . 
The consensus sequence obtained here is somewhat similar to the ResD consensus sequences identified earlier based on in vitro studies ( 38 ) . 
The current work supported the previously detected in vitro ResD-binding architecture but also demonstrated the importance of TGA for ResD binding as well as for Fur-operator interaction . 
The GA sequence in TGANAA of the Fur operator site is important for the recognition by Fur and functions as a discriminatory core sequence distinguishing the Fur operator from the sites recognized by other Fur homologs ( Per and Zur ) ( 37 ) . 
The G base of TGA in the ResD-binding site is less conserved compared with that in the Fur site , which might act to favor Fur over ResD . 
Unexpectedly , another C-rich motif was identified by MEME analysis ( Data Set S3 ) , and interestingly , the reverse complement sequence to this motif was identified as one of the common motifs for NsrR targets ( Data Set S4 ) . 
At the moment , it is unknown whether this sequence is involved in ResD - and NsrR-binding recognition or serves as a binding site for other DNA-binding proteins , including TFs . 
NsrR . 
A MEME motif search determined that 9 of 26 NsrR-binding sites carry a similar motif ; however , this motif is not present in more than half of the NsrR-binding sites ( Data Set S4 ) , including hmp . 
Accordingly , the automatic search did not detect the motif in the nasD NsrR-binding site that is similar to the hmp site . 
As described above , all previous in vivo and in vitro work and a computational study indicated that a partial dyad symmetry sequence of 8-1-8 ( ATRTATYTTAAATAtat , where R is G or A and Y is C or T [ the last three bases are not conserved between NsrR-binding sites in nasD and hmp ] ) is the cis site required for NO-sensitive repression by NsrR ( 6 , 9 , 16 ) . 
Base substitutions of every residue in the sequence revealed that T and A at the fourth and fifth positions in the upstream half-site are important , and A at the fifth position in the downstream half-site is the most critical residue for in vivo and in vitro NsrR activity . 
In addition , the deletion of the center T , but not its substitution to G , resulted in the loss of NsrR activity ( 9 ) . 
The downstream part of the NsrR motif identified by MEME ( Fig. 4 ; see also Data Set S4 ) , namely , ATATTATGTAAAC , does not show dyad symmetry ; however , it could be read as ATATYATGTAAAC , which somewhat resembles RTATYTTA AATAt ( italicized letters show important residues for NsrR activity ) in the NsrR-binding site of nasD , and importantly , the critical distance between TA and A ( TAN8A ) is conserved . 
The results suggest that the consensus sequence emerging from the analysis presented above might function as an NsrR operator site , at least for certain genes ( Data Set S4 ) . 
In conclusion , comparison of the proposed consensus sequences for Fur , ResD , and NsrR brought to light that the A T-rich sequences have a moderate similarity to each other , which could explain the overlapping DNA contacts of these TFs . 
Interestingly , the proposed Fur and NsrR consensus sequences contain TGA and TGT , respectively , whereas the ResD consensus sequence has TGA and TGW . 
Any variation departing from each consensus sequence likely determines which TF more strongly occupies these sites to control transcription . 
Consensus sequences are identified within a narrow range of DNA where the three TFs competitively bind within promoters . 
ResD and NsrR bind to narrow regions within the ykuN and fbpC promoters only in the presence of the other two TFs , while others such as the ymfC promoter are sites where ResD/NsrR/Fur ( or NsrR/ResD ) competitively bind ( Table 2 ; see also Data Set S1 ) . 
As a proof of concept , we determined whether the proposed consensus sequences ( Fig. 4 ) could be identified in the overlapping binding sites of these promoters . 
Fur enrichment to the ykuN promoter is 30 - to 50-fold stronger than for NsrR and ResD ( Fig. 5A ) , respectively , which mirrors the presence of two strong Fur consensus sequences within the overlapping binding site . 
We reported previously that efficient binding of NsrR and ResD to ykuN depend on each other , and neither protein is enriched at ykuN in the fur mutant ( 14 ) . 
The GeF-seq results of NsrR and ResD binding are consistent with the previous findings . 
On the other hand , the resD nsrR mutation led to a shift of the Fur-binding site to the core promoter region . 
Because of the presence of the strong Fur consensus sequences , it was difficult to pinpoint NsrR and ResD consensus sequences , although some similarity to their consensus sequences is detected within the binding site . 
Fur binding is relatively weak compared with ResD and NsrR binding at the ymfC , yoeC , ydhB , yobR , and yclA promoters and is only visibly detected at pps . 
NsrR was enriched at these promoters to a greater extent than ResD or Fur . 
When the binding region of ymfC was expanded ( Fig. 5B ) , we clearly observed that the three TFs bound within a 40-bp segment . 
Furthermore , the consensus sequences for Fur , ResD , and NsrR proposed in Fig. 4 were detected within a 25-bp sequence that carries multiple TGAs and TGTs . 
The results demonstrated that the high-resolution mapping is a powerful tool for identifying each consensus sequence for multiple TFs that bind to overlapping targets , if the DNA affinity of each is comparable . 
ResD and NsrR bind to sequences outside promoter regions . 
ResD and NsrR mostly bind to promoter sites and ResD often binds to entire coding regions of positively controlled genes , but they occasionally show a narrow peak of binding in coding regions . 
Such examples are recombination sites of SP ( yodU-yotN and ykoA-ypgP ) and sigK intervening ( skin ) element ( 39 ) ( spoIVCB-spoIVCA and yqaB-spoIIIC ) , to which both ResD and NsrR bind ( Table 2 ) . 
The B. subtilis genome carries 10 prophagelike sequences identified as A T-rich regions ( 40 ) . 
Interestingly , the SP prophage and skin are inserted into protein-coding regions ( 41 , 42 ) . 
Upon DNA damage ( caused by mitomycin C or UV ) , SP is excised from the genome and begins lytic development . 
The skin element interrupts the sigK gene encoding one of the sporulation mother cellspecific sigma factors . 
At late sporulation , two parts of genes interrupted by skin , spoIVCB and spoIIIC , are connected by site-specific recombination in the mother cells to generate the intact sigK gene . 
Recent work showed that the site-specific recombination of the SP region also occurs during sporulation , which generates the functional spore polysaccharide gene spsM from truncated yodU and ypqP ( 43 ) . 
The sequences that ResD and NsrR bind cover the sites of recombination catalyzed by recombinases SprA ( for SP ) ( 43 ) and SpoIVCA ( for sigK ) ( 39 ) ( Data Set S1 ) . 
ResD and NsrR bind to the upstream site of ypqP that encodes the C-terminal part of SpsM ( see Fig . 
S2 ) . 
The NsrR-binding sequence partially overlaps with attR , and the ResD-binding sequence covers half of the inverted repeat of the SP recombination site . 
As expected from the relatively low intensity of both NsrR and ResD binding to yodU encoding the N-terminal part of SpsM , both binding sites show less similarity to each TF 's consensus sequence . 
NsrR again binds to the site that overlaps with attL , and ResD binds to the other half of the inverted repeat , except on the opposite strand . 
Excision of the phage will leave an intact NsrR-binding site , while the ResD-binding region is excised with the phage . 
It remains to be determined whether the binding of ResD and/or NsrR to these sites plays any role in site-specific recombination during sigK generation and/or the life cycle of SP phage . 
The TFs , particularly Fur , bind within a narrow segment of coding-region DNA or at the 3 = end of certain genes . 
In addition to prophage-like insertion sites , GeF-seq provided evidence that the TFs bind to other locations outside promoter regions . 
For example , Fur binds to 32 and 48 promoter DNAs in the wild-type and resD nsrR mutant strains , respectively , but Fur also binds to 65 and 50 coding regions in the wild type and the resD nsrR mutant , respectively . 
These sites include resE , hemX , gltT , and kinB . 
ResD and NsrR only bind to coding regions of ppsB ( plipastatin synthase gene ) in the wild type , whereas Fur binds in both the wild type and the mutant ( Table 1 ) . 
NsrR interacts with the ydjA-coding region in the resD fur mutant , but Fur only weakly binds to the site , while ResD shows no interaction . 
Although yodU , spoIVCB , and spoIIIC belong to this category , these sites function as targets of site-specific recombination as described above . 
Why Fur binds to many coding regions and how this binding affects gene expression are unknown , although a road-block mechanism has been reported in the binding of the CodY TF to a coding region ( 44 ) . 
Among the binding sites at the 3 = ends of genes , the regions between ynfE and xynC , as well as between yoaM and yozS , contain sites where the three TFs bind in the mutant strains ( Table 2 ) . 
In these cases , NsrR appears to bind more strongly than ResD or Fur , a pattern similar to that for ymfC described above ( Fig. 5B ) . 
Consistent with this finding , we were able to detect the consensus sequences of all three TFs within a 30-bp region where all three bind ( Fig. 5C and D ) . 
We noticed that the binding sites likely reside within upstream regions of unassigned genes in the strain 168 genome . 
Small hypothetical proteins were annotated to both regions in the genomes of other B. subtilis strains according to NCBI . 
Although both genes carry putative 10 sequences , a Shine-Dalgarno sequence is missing ( the putative open reading frame located downstream of yoaM ) or weak ( the one located downstream of ynfE ) , and whether the regions code for protein or RNA needs to be investigated . 
DISCUSSION
Genome-wide transcriptome analysis , such as RNA-seq , is a powerful method to identify genes controlled by each TF . 
However , the method is not always applicable to identifying genes that are directly regulated by TF-DNA contact . 
Particularly , it is difficult to determine direct or indirect control when a TF exhibits less sequencespecificity toward its target and/or is influenced by multiple TFs binding within the target 's vicinity . 
Here , we investigated genome-wide binding of Fur , ResD , and NsrR . 
NsrR was originally considered a repressor that binds to DNA in a sequence-specific manner ( 9 , 10 , 16 ) , but later studies raised the possibility that it interacts with promoters lacking the putative consensus sequence identified in nasD and hmp ( 9 , 14 ) . 
We used GeF-seq to detect in vivo binding of the three TFs . 
Compared with our previous ResD and NsrR ChAP-chip data ( 14 ) , GeF-seq enabled us to identify binding sites with a much higher resolution . 
The results confirmed the previous finding that the method resolves protein-DNA binding at a resolution similar to that obtained by in vitro DNase I footprinting ( 15 ) . 
For example , Fur binding identified by GeF-seq achieved a comparable resolution to that obtained by in vitro DNase I footprinting ( 37 ) ( Data Set S5 ) . 
The three TFs bind cooperatively within a 75-bp DNA sequence ( Fig. 5A ) and competitively when targeting a segment smaller than 40 bp ( Fig. 5B to D ) . 
We examined whether these DNA sites contain each consensus sequence . 
One such site , the ykuN promoter , is extremely complex . 
As marked with arrows for Fur , direct and inverted repeats of four TGAW3 motifs were detected as two clusters downstream of 35 and the transcription start site ( Fig. 5A ) . 
Identical sequences ( TGAAAATCATTATCA ) of the 7-1-7 Fur consensus motif reside within both TGAW3 clusters . 
In addition , as TGANAANW3 might be a ResD-binding sequence , we were unable to draw a reasonable conclusion to describe the ResD/NsrR-DNA binding location . 
On the other hand , we successfully identified all consensus sequences in the TF-binding sites of ymfC , pps , yobR , and two unassigned genes downstream of ynfE and yoaM ( Fig. 5B and data not shown ) . 
The highest NsrR-binding intensity among the three TFs is a common characteristic to these binding sites . 
Although this characteristic applies to genes , including ydhB and yclA , the NsrR consensus sequence was not detected in these genes . 
Given that more than half of the NsrR-binding sites do not carry the putative consensus sequence , the entire picture of NsrR-operator binding is not fully uncovered . 
This study showed that a part of the ResD population likely associates with RNAP during the elongation of certain genes . 
Although ResD interaction with RNAP was previously reported ( 26 ) , our work revealed that a majority ( if not all ) of the genes to which ResD binds across coding regions are activated by ResD and play roles in aerobic/anaerobic respiration . 
It is not clear whether one such gene , ndk , has any role in respiration , although a possible involvement of Ndk ( nucleotide diphosphate kinase ) in energy production was previously indicated by a study in which Ndk was shown to function in ATP generation from phosphoenolpyruvate ( PEP ) in an oxygen-independent manner via a phosphotransfer network in Pseudomonas fluorescens ( 45 ) . 
Ndk that interacts with the inner mitochondrial membrane was shown to couple nucleotide transfer with respiration ( 46 ) . 
The functions of yxiE and yozB remain to be investigated . 
yxiE is induced by phosphate starvation independently from the PhoPR two-component regulatory system and B , which is a rare exception among phosphate starvation-induced genes ( 47 ) . 
Transcription of ctaO , ydbL , and yfmQ is higher under anaerobic than under aerobic conditions ( 27 ) , and YdbL and YfmQ , like CtaO , are predicted to localize to the cell membrane . 
yizD transcription is activated by oxygen limitation ( 27 ) , and the divergently transcribed small gene was identified as a homolog of RsaE , an sRNA widely conserved in Gram-positive bacteria ( 48 ) . 
The B. subtilis RsaE ( now renamed RoxS ) is activated by ResDE in response to NO ( 49 ) . 
Genes involved in redox homeostasis are upregulated in the roxS mutant , which include cydA , ykuN , and resA that were a focus of our study . 
Therefore , ResD negatively regulates these genes through RoxS , while it upregulates cydA , ykuN , and resA as a direct activator , an antirepressor , and indirectly via an as-yet-unidentified regulatory pathway , respectively . 
The small basic protein YizD ( 55 amino acids [ aa ] , pI 9.7 ) might function as a RoxS chaperone . 
All our data in earlier studies led us to believe that ResD - and NsrR-dependent transcriptional control of hmp and nasD operate similarly , but the current study revealed that the assumption is not fully supported and that the proposed mechanism will need to be revised based on the results from future experimentation . 
One unexpected observation is that ResD-RNAP binding around the nasD transcription start site was much stronger than binding to the upstream sites previously identified by in vitro footprinting experiments ( Fig. 1B ) . 
On the other hand , ResD interacts with the hmp upstream regulatory site as well as the transcription start site ( Fig. 1A ) . 
The difference might be due to the core promoter ( a suboptimal 10 for hmp and a strong TG 10 extended promoter motif [ TG 10 ] element for nasD ) . 
ResD binding at the upstream site might be released upon recruitment of RNAP to the nasD promoter , while the 
RNAP-ResD complex might be stalled at the promoter due to the stable RNAP/TG 10 promoter DNA complex before promoter escape . 
The upstream ResD-binding sites in nasD overlap a CodY-binding sequence ( 50 ) , and a TnrA-binding sequence resides between the CodY and NsrR sites ( 51 ) ( Fig. 1B ) . 
TnrA activates the nasD operon only under poor nitrogen conditions and not under the anaerobic fermentation conditions used in this study ( 19 ) . 
Because the CodY-binding site overlaps with the ResD-binding site , we examined whether CodY has any role in anaerobic nasD expression . 
For example , CodY might play a positive role by binding to the site , either as an activator , coactivator , or an antirepressor . 
nasD expression was upregulated by the nsrR mutation , but the expression in the codY mutant was as low as that in the wild type . 
nasD-lacZ expression in the codY nsrR mutant is similar to that in the nsrR mutant ( data not shown ) , indicating that CodY does not play a major role in nasD transcription under anaerobic conditions . 
However , we observed a reproduc-ible minor growth defect in the absence of codY . 
MATERIALS AND METHODS
Bacterial strains and culture conditions . 
Strains and plasmids used in this study are listed in Table S1 in the supplemental material . 
All B. subtilis strains used in this study are derivatives of strain 168 unless otherwise stated . 
Strains OC0010 , ORB8238 , and ORB8440 expressing C-terminally His12-tagged nsrR , resD , and fur at their native loci were constructed previously ( 14 ) . 
These strains produce the 12 His-tagged proteins that are functional in vivo based on the effects of each TF in transcriptional control of target genes as previously described ( 14 ) . 
Strains expressing His12-tagged nsrR in the resD fur double mutant background ( ORB9091 ) , His12-tagged resD in the nsrR fur mutant background ( ORB9092 ) , and His12-tagged fur in the nsrR resD mutant background ( ORB9093 ) were also constructed . 
To this end , chromosomal DNA of HB2501 ( fur : : neo ) was used to transform OC0010 to generate ORB8277 ( nsrR-His12-tet fur : : neo ) , which was then used for transforming with MH5260 ( resD : : cat ) chromosomal DNA , resulting in ORB9091 ( nsrR-His12-tet fur : : neo resD : : cat ) . 
Similarly , HB2501 and ORB6179 ( nsrR : : cat ) chromosomal DNA was used to construct ORB9092 ( His12-resD fur : : neo nsrR : : cat ) , and ORB6179 and LAB2511 ( resD : : spc ) DNA was used to construct ORB9093 ( fur-His12 nsrR : : cat resD : : spc ) . 
The successful construction of His-tagged genes was confirmed by PCR , followed by sequencing . 
To determine the regulatory roles of the TFs that interact with specific promoter regions , transcriptional lacZ fusions were used . 
The lacZ fusions used are those integrated at the native loci or ectopic loci ( thrC or amyE ) ( see Table S1 ) . 
- Galactosidase activities were measured in the wild type and in mutant backgrounds that lack resD , fur , and/or nsrR as indicated . 
B. subtilis strains were grown anaerobically in 2 YT supplemented with 0.5 % glucose and 0.5 % pyruvate ( fermentative conditions ) where NsrR is active ( 6 ) . 
Genome footprinting by high-throughput sequencing . 
GeF-seq was carried out as described in the previous report on genome-wide AbrB-binding in B. subtilis ( 15 ) with some modifications . 
In short , OC0010 ( nsrR-His12 ) , ORB8238 ( resD-His12 ) , or ORB8440 ( fur-His12 ) was grown anaerobically at 37 °C in 2 YT with 0.5 % glucose and 0.5 % pyruvate until the end of exponential growth ( T0 ) and cells were treated with formaldehyde . 
For the mutant backgrounds , ORB9091 ( nsrR-His12 resD : : cat fur : : neo ) , ORB9092 ( resD-His12 nsrR : : cat fur : : neo ) , and ORB9093 ( fur-His12 resD : : spc nsrR : : cat ) were used and cultured similarly . 
Cells were harvested and treated with 5 mg/ml lysozyme in isotonic SMM buffer ( 0.02 M maleic acid , 0.5 M sucrose , and 0.02 M MgCl2 , pH 6.5 adjusted with NaOH ) in the presence of EDTA-free protease inhibitor cocktail and 1 mM phenylmethylsulfonyl fluoride ( PMSF ) before DNase I and RNase A treatment . 
DNase I concentrations were adjusted to the level that generates fragments less than 100 bp in size based on agarose gel electrophoresis . 
The reactions were terminated by the addition of UT buffer ( 0.1 M HEPES [ pH 7.5 ] , 0.01 M imidazole , 8 M urea , 0.5 M NaCl , 1 % Triton X-100 , 10 mM - mercaptoethanol , 1 mM PMSF ) . 
The samples were then sonicated on ice using an Astrason Ultrasonic Processor XL ( Misonix ) for 10 min ( 4 s `` on '' and 10 s `` off '' at output level 5 ) . 
After removing cell debris , protein-DNA complexes were affinity-purified using Dynabeads His-Tag isolation and pulldown ( Thermo Fisher Scientific ) and were reverse cross-linked by incubating at 65 °C overnight . 
The DNA library for sequencing by the Illumina Genome Analyzer IIx ( GAIIx ) or Illumina HiSeq2000 was generated using the NEBNext DNA Sample Prep reagent kit ( New England BioLabs ) according to the manufacturer 's instructions with some modifications . 
To reduce the loss of small-sized DNA fragments , all DNA purification steps before adapter ligation were carried out using a QIAquick nucleotide removal kit ( Qiagen ) . 
In addition , after adapter ligation to the DNA fragments , the size selection step that was done previously was eliminated to recover all sizes of DNA fragments bound by target protein . 
The DNA fragments were then purified using a QIAquick PCR purification kit ( Qiagen ) and were amplified using PCR to obtain at least 1 fmol of DNA library . 
The amount of DNA was determined using an Agilent 2100 Bioanalyzer with a high-sensitivity DNA kit ( Agilent ) . 
The sequence of the library was then determined by 69-bp or 100-bp paired-end sequencing using Illumina GAIIx or Illumina HiSeq2000 , respectively . 
GeF-seq data analysis . 
The mapping of read sequences to a reference genome was carried out using the mpsmap software ( 52 ) . 
The detection of depth peaks and visualization of the results were carried out with methods previously described ( http://metalmine.mydns.jp/maps/gefseq ) with some modifications ( 15 ) . 
In the present analysis , the positions of both ends of DNA fragments had been determined based on paired-end read information , using pair_map software we developed ( http : / / metalmine.mydns.jp / maps/pair _ map ) . 
Regions between the starting positions of paired reads within the range of 25 to 100 bp were collected as peak candidates . 
To detect the majority of potential peaks , the regions where signal intensities exceeded threshold value ( 50 for NsrR binding in the wild-type strain and 100 for others ) at more than half of the nucleotides between them were extracted . 
The detected peaks include background peaks derived from preferential DNA digestion of A T-rich regions by DNase I and thus are consistent with the local G C content at every 50 bp as shown in Fig. 1 . 
The maximum read count was used to represent the intensity of each peak . 
The distribution of intensities of primary peaks shows a normal distribution caused by the background peaks derived from the G C content , and the equation was evaluated using Gaussian curve fitting of the IGOR Pro software ( Hulinks Inc. ) . 
To extract the peaks caused by protein binding from the background , P values of each peak were calculated and all peaks with P values lower than the threshold ( indicated in the supplemental material ) were assigned as protein-binding peaks . 
All positive peaks were then reconfirmed by visual inspection using IMC software ( In Silico Biology , Inc. ) . 
Ambiguous peaks were removed manually in this process . 
Finally , to normalize the read number variation among experiments , read counts at each position of the genome were multiplied by 500 million and then divided by the total read count . 
Motif analysis . 
NsrR - , ResD - and Fur-binding DNA motifs were analyzed by the BiPad web server ( http://bipad.cmh.edu ) for modeling bipartite sequence elements or the MEME program ( http://meme - suite.org/tools/meme ) . 
To search bipartite motifs for Fur , BiPad was used with the following settings : left half-site , gap range lengths , right half-site , and the iteration cycles were set to 7 , 0 to 2 , 7 , and 100 , respectively ( 35 ) . 
Using the results , sequence logos of the bipartite motifs with 0 , 1 , and 2 gaps were generated by WebLogo ( 53 ) . 
Accession number ( s ) . 
Sequencing data in this study were submitted to the DDBJ Sequence Read Archive ( DRA ) under accession number DRA005609 . 
SUPPLEMENTAL MATERIAL
Supplemental material for this article may be found at https://doi.org/10.1128/JB .00086 -17 . 
ACKNOWLEDGMENTS
We thank Adriano O. Henriques , Tsutomu Sato , Boris Belitsky , and Linc Sonenshein for providing B. subtilis strains . 
We thank Kazuo Kobayashi for sharing unpublished transcriptomic results used in Table 1 . 
We also thank Peter Zuber and Naotake Ogasawara for critically reading the manuscript . 
This work was supported by grants from the National Science Foundation ( grant number MCB1157424 to M.M.N. ) , the Advanced Low Carbon Technology Research and Development Program ( ALCA ) of the Japan Science and Technology Agency ( JST ) ( grant number JP15K07359 to S.I. ) , and KAKENHI of the Japan Society for the Promotion of Science ( JSPS ) ( grant number JP26430199 to K.N. ) . 
T.Q. was supported by a Partners in Science Program Grant from the M. J. Murdock Charitable Trust .