Influence of Glucose Availability and CRP Acetylation on the Genome-Wide Transcriptional Response of Escherichia coli: Assessment by an Optimized Factorial Microarray Analysis

Background: While in eukaryotes acetylation/deacetylation regulation exerts multiple pleiotropic effects, in Escherichia coli it seems to be more limited and less known. Hence, we aimed to progress in the characterization of this regulation by dealing with three convergent aspects: the effector enzymes involved, the master regulator CRP, and the dependence on glucose availability. Methods: The transcriptional response of E. coli BW25113 was analyzed across 14 relevant scenarios. These conditions arise when the wild type and four isogenic mutants (defective in deacetylase CobB, defective in N(ε)-lysine acetyl transferase PatZ, Q- and R-type mutants of protein CRP) are studied under three levels of glucose availability (glucose-limited chemostat and glucose-excess or glucose-exhausted in batch culture). The Q-type emulates a permanent stage of CRPacetylated, whereas the R-type emulates a permanent stage of CRPdeacetylated. The data were analyzed by an optimized factorial microarray method (Q-GDEMAR). Results: (a) By analyzing one mutant against the other, we were able to unravel the true genes that participate in the interaction between ΔcobB/ΔpatZ mutations and glucose availability; (b) Increasing stages of glucose limitation appear to be associated with the up-regulation of specific sets of target genes rather than with the loss of genes present when glucose is in excess; (c) Both CRPdeacetylated and CRPacetylated produce extensive changes in specific subsets of genes, but their number and identity depend on the glucose availability; (d) In other sub-sets of genes, the transcriptional effect of CRP seems to be independent of its acetylation or deacetylation; (e) Some specific ontology functions can be associated with each of the different sets of genes detected herein. Conclusions: CRP cannot be thought of only as an effector of catabolite repression, because it acts along all the glucose conditions tested (excess, limited, and exhausted), exerting both positive and negative effects through different sets of genes. Acetylation of CRP does not seem to be a binary form of regulation, as there is not a univocal relationship between its activation/inhibitory effect and its acetylation/deacetylation stage. All the combinatorial possibilities are observed. During the exponential growth phase, CRP also exerts a very significant transcriptional effect, mainly on flagellar assembly and chemotaxis (FDR = 7.2 × 10−44).


INTRODUCTION
It is well known that microbial cells continuously sense their external and internal environment, with this information being filtered and integrated to trigger responses that allow them to best adapt to these stimuli (Snitkin and Segrè, 2008;Zhao et al., 2016). Beyond the "short-term adaptations" provided by the allosteric regulations working alongside the extensive set of protein post-translational modifications, changes in the transcriptional responses are the main mechanisms for achieving the "long-term adjustments" required (Shimizu, 2014). So, there are good reasons to analyse changes in the global transcriptome of microorganisms. To this end, one technique that is widely used is the high-density array of oligonucleotides (microarrays).
In fact, microarrays enable the analysis of the complete Escherichia coli transcriptome by monitoring the labeled hybridization of around ten thousand probes. Moreover, transcriptome analysis also provides information about the non-coding RNAs (anti-sense RNAs and microRNAs), which have increasing importance as additional layer of regulation at post-transcriptional level, not only in eukaryotes but also in bacteria (Delihas, 2015;Cech et al., 2016;Tronnet et al., 2016).
Here, we will focus on the analysis of changes in the transcriptome of some E. coli strains (wild-type and isogenic mutants) when these bacteria are challenged with different levels of a typical carbon source such as glucose. The analyses include three well-defined scenarios: the case of unrestricted glucose availability (exponential phase of batch culture), the case of glucose deprivation (stationary phase of batch culture), and the stage of transition between these extremes (glucose-limited chemostat culture). Of course, it must be noted that multiple experimental studies of this type have been done over the last 10 years (see Table 1). Also, some reviews are available (Fic et al., 2009;Hu et al., 2010;Bernal et al., 2014;Choudhary et al., 2014;Drazic et al., 2016;Wolfe, 2016). What then justifies the present study?
Inspection of Table 1 reveals some important aspects that justify the need to perform the present study. Hence, only a few studies simultaneously integrate the three biological scenarios employed herein. Moreover, our study aims not only to analyse the transcriptional effect of the glucose availability, but also to establish the extent to which these effects are specifically regulated through acetylation/deacetylation of the CRP protein.
In fact, this protein is the main protein responsible for catabolite repression (Deutscher, 2008). Together with RcsB and RpoB, CRP is one of the three main transcriptional regulators known in E. coli as influenced by acetylation (Zhang et al., 2009).
There are not many previous studies addressing catabolite repression from the point of view of its regulation by acetylation/deacetylation of CRP. Often, when studies referred to regulation by acetylation, they were not actually addressing the problem of glucose availability (Lima et al., 2011;Ma and Wood, 2011;Weinert et al., 2013a;Zhang et al., 2013;Baeza et al., 2014). Moreover, the mechanism of lysine-acetylation/deacetylation is also present in different proteins of E. coli as a type of posttranslational modification (Ma and Wood, 2011;Liu et al., 2014;Castaño-Cerezo et al., 2015;Fraiberg et al., 2015;Drazic et al., 2016).
Although most of the studies in Table 1 used microarrays, it is well known that this tool has major limitations (Guarnaccia et al., 2014;Yang et al., 2014;Chrominski and Tkacz, 2015). Many of these limitations derive from the current algorithms applied, such as the Empirical Bayes and Benjamini-Hochberg algorithms. To avoid these problems, we have developed an alternative algorithm called Q-GDEMAR (Quantile-Gaussian deconvolution of microarrays), which provides a better sensitivity of detection together with a low false discovery rate (FDR) .
Further, all the studies in Table 1 only considered pair-wise comparisons. The microarray data in our analyses, however, will be arranged as a factorial design, thus allowing us to account for the interaction between the variables. Furthermore, unlike the studies shown in Table 1 that "summarized" their microarrays (Irizarry et al., 2003), to optimize our study, the microarray data will be analyzed here at the probe-level. In any case, of special importance is the fact that all the technical variants that we will use for the microarray analyses in the present investigation have recently been discussed and tested successfully .
Note in Table 1 the occurrence of additional sources of variability, which affect the comparability of results among the studies. Thus, due to the different strains used by different authors, significantly different regulatory features can be expected, even though all the strains assayed belong to the E. coli taxon (Yoon et al., 2012;Castaño-Cerezo et al., 2015;Vital et al., 2015;Monk et al., 2016). Moreover, the presence of different trademarks in the microarray studies in Table 1 might imply a higher prevalence of false responses (cross-reactivity) in some cases due to the shorter length of probes in some kits (Chou et al., 2004). Importantly, the interpretations are conditioned not only by the type of computing algorithm used, but also by the application of different cut-off criteria (for example, see in Table 1 the use of t-student arbitrary values together with non-adjusted p values for multi-comparisons; p adjusted by Bonferroni; p adjusted by FDR). Taken together, these factors inexorably lead to the identification of different sets of differential genes and hence produce different biological inferences.
In brief, herein we advance in the characterization of regulatory acetylation mechanisms exerted by the DNA-binding dual regulator CRP (CRP: cAMP receptor protein) and their relation with glucose availability. To this end, we re-analyse two sets of previously published data, but unlike in the original reports, herein the transcriptional responses are crossed with the different conditions of glucose availability, and the microarray data are dealt with using our optimized post-processing approach by factorial Q-GDEMAR. In addition, the criteria applied to consider the control conditions are also reformulated.

Strains and Conditions Tested
Two sources of microarray data are re-analyzed herein. The first set concerns the global transcriptional responses to defective acetylation-deacetylation E. coli mutants. This data set has been gathered by Castaño-Cerezo et al. (2014), and was retrieved for our study from the GEO database (access code: GSE62094). The second set concerns the overall transcriptional responses of E. coli to mutations in the crp gene, which codifies the DNA-binding dual regulator CRP. This information has been gathered by Écija-Conesa and de Diego (2007) (GEO database, access code: 96955) 1 and published latter (Davis et al., 2018). These set of data are highly consistent since they came from the same laboratory, uses the same host E. coli strain, culture media, operational conditions, and analytics. Moreover, the analyses of both data sets are independent one of the other.
In the first case, the analysis encompasses the comparison of transcriptional responses to six experimental conditions. These scenarios result from testing three variants of the E. coli BW25113 strain (wild type, patZ, cobB) when cultured under two modes of bioreactor operation (batch and chemostat). The cobB mutant refers to a strain defective in the cobB gene (ycfY), which encodes the protein-lysine deacetylase and desuccinylase enzyme (Uniprot: P75960). This enzyme belongs to the family of sirtuins and is NAD + -dependent. By contrast, the mutant patZ refers to a strain defective in the patZ gene (yfiQ), which encodes the peptidyl-lysine acetyl-transferase enzyme (Uniprot: P76594). The samples coming from the exponential phase of the batch (i.e., culture without glucose limitation) correspond to a specific growth rate µ max = 0.62 h −1 . The samples coming from the chemostat (glucose-limited culture) correspond to a Dilution Rate = 0.23 h −1 and were taken at the steady state achieved after five residence times.
In the second case, the analysis encompasses the comparison of the transcriptional responses to eight experimental conditions. These scenarios result from testing four other E. coli BW25113 genetic variants [wild type; N (a cloning control); crp(Rmutant), crp (Q-mutant)] under two stages of batch culture (exponential and stationary phase). In fact, in one of these mutants, the strain E. coli BW25113 was engineered to produce a modified CRP protein in which the Lysine (K) 100, a critical residue for the CRP activity (Baeza et al., 2014), was substituted to codify an Arginine (R-mutant). In the second type of mutant, the K100 was substituted to codify a glutamine (Q-mutant). Importantly, the objective of the mutant crp(R-type) is to emulate a constitutive stage of deacetylation in Lysine 100 of CRP, whereas the objective of the mutant crp(Q-type) is to emulate a stage of constitutive acetylation of Lysine 100 in the same protein.
In this last series of experiments, an E. coli strain BW25113 defective in the crp gene was used as host of the cloning. The authors indicate that since they could not insert the modified genes at the native position, the variants of the crp gene (Qtype and R-type) were cloned on the paaH gene. This gene, which codes for 3-hydroxyadipyl-CoA dehydrogenase enzyme, is located within the paaABCDEFGHIJK operon. So that the cloned variants could be transcribed independently on the hosting operon promoter, the gene variants were cloned together with the promoter corresponding to the native crp gene. Moreover, to account for the effects caused by the change in the position, a third type of mutant that works as cloning control, called N, was considered. Thus, the "normal" crp gene and its associated promoter were cloned at the position of the paaH gene as well.
In both series of experiments, the bacterial cultures were done at 37 • C, using minimal defined medium (MM9) with glucose 1 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE96955 10 mM (pH = 7.4) and full aerobic conditions. However, while in the first series each condition was run in triplicate, in the second series each condition was run in duplicate. The samples corresponding to the exponential phase were taken at an optical density = 0.5, while the samples corresponding to the stationary phase samples were taken at an optical density = 1.8. For the microarray determinations, the authors indicate that RNA was extracted and purified from the culture samples, subjected to current quality control procedures, and finally run on Gene Chip E. coli Genome 2.0 (Affymetrix, USA) (for details see the original articles). Importantly, this microarray kit comprises 10,207 probes, from which only 8,662 have a gene assigned.

Factorial Microarray Analysis
The log 2 -normalized data of the microarray are subjected herein to a post-processing elaboration by the Q-GDEMAR method . This performs a computational deconvolution of the central region of the data distribution. The parametric characterization of this region in terms of a Gaussian distribution provides narrower limits to the genes whose expression fluctuates only stochastically. The comparison of these limits with the overall data distribution ultimately allows us to determine, with greater sensitivity and lower FDR, which probes are being differentially hybridized. The specific protocol for 2 × 2 factorial microarray analysis and the justification for why the probes are not subjected to median-polishing summarization have been described elsewhere .

Ontology Analysis
This was performed by the software DAVID 6.8 (Huang et al., 2009;Jiao et al., 2012). This is a free web service (https://david. ncifcrf.gov/). We have used the default settings provided by the program, while contrasts were performed against the list corresponding to the whole genome of E. coli as background.

Substrate Availability and Genome-Wide Transcription
In previous analyses of Castaño-Cerezo et al. (2014), the focus was placed on the study of two mutant strains ( cobB and patZ). There, the wild-type strain was considered only as the "control" against which the effects of the mutations were contrasted. However, besides assessing the performance of the mutants by using an optimized microarray post-processing, we also analyse herein the performance of the wild-type strain itself (see Table 2).
In Table 2, a ratio value significantly >1 means that genes corresponding to the condition considered at the numerator are up-regulated (or equivalently, that the genes corresponding to the condition considered at the denominator are down-regulated). Conversely, a ratio value significantly <1 means that genes corresponding to the condition considered at the denominator are up-regulated (or equivalently, that the genes corresponding to the condition considered at the numerator are downregulated).
With these criteria of interpretation in mind, several aspects of Table 2 are noticeable. First, we note the small number of differentially expressed genes observed in each condition in spite of the high number of probes present in the microarray (n = 8,662). Second, we note the low values of FDR  achieved throughout the analysis by the optimized microarray post-processing used herein. Note that the FDR values computed cover the entire set of associated genes .
Third, most of the genes operating during the exponential growth phase (where the wild-type strain achieved a µ max = 0.62 h −1 ) seem to continue being transcribed even during the stationary phase (where the growth of the strain has ceased). In fact, only six genes show a significant enhanced transcription during the exponential growth phase (second row, Experiment B), whilst 56 genes appear dysregulated during the stationary phase (first row, Experiment B). The small group of genes dysregulated during the exponential growth phase codify proteins involved in the metabolism of sulfur compounds such as methionine, cysteine, S-adenosyl methionine (cysD, cysP, metE genes), the novo biosynthesis of pyrimidines and glutamine metabolism (pyrD, carA genes), the biosynthesis of phenylalanine, tyrosine, and tryptophan (aroG gene).
Importantly, by the ontology analysis of the 56 differential genes detected for the stationary phase with the wild-type strain we identified as significant classes a set of four genes that codify periplasmic proteins with dodecin-like structure (ECs5165, ECs0884, ECs5164, ECs5175; FDR = 2.1 × 10 −2 ) and another set of five genes which codify enzymes involved in Figure 1).
The biological role of the dodecin-like proteins is not wellknown. However, based on their structural similarity with the dodecin molecule, these proteins could also have a protective role due to their capacity to sequester molecules containing oxidized FIGURE 1 | Arginine catabolic routes mapping the sub-set of five enzymes codified by the genes dysregulated in the microarray data corresponding to the stationary phase of batch culture with the wild type E. coli strain (first row, Experiment B, Table 2). The enzymes are indicated by their Enzyme Commission (EC) number, while the names of the genes that codify them are indicated in red. The diagram is a simplified adaptation based on KEGG Pathways (map 00330) after Ontology DAVID 6.8 software processing (FDR = 8 × 10 −4 ).
Frontiers in Microbiology | www.frontiersin.org flavin-moieties (Grininger et al., 2009). On the other hand, arginine, in addition to its structural role as amino acid, is also required as N donor for the biosynthesis of the messenger nitric oxide (NO). Moreover, through its transformation to ornithine, arginine is used for the biosynthesis of polyamines (putrescine and spermidine). Thus, is inferred that these two activities of arginine might be affected during the stationary phase.
Other groups of genes also appeared up-regulated during the stationary phase although they did not achieve the threshold of significance in the ontology analysis. In fact, we detected a group of five genes belonging to the hyaABCDE operon that codify the different subunits of the enzyme hydrogenase 1, two genes related with biofilm formation (bssR, and bsmA), several genes related with stress response (gadC, ytfQ, csiE, ymjC, yjbJ), and three genes related with the carbohydrates usage (araF, ytf R, sgcC). Noticeably, around 50% of the total differential genes detected at the stationary phase correspond to hypothetical proteins (see Table S2).
Consistent with our claim about the persistence of most of the genes operating under the exponential phase in further culture phases, from Experiment A in Table 2 it can be seen that only 104 differential genes appear to be up-regulated due to the glucose limitation in "chemostat 1" (Dilution rate = µ = 0.23 h −1 , a representative fixed point within the interval of the progressive deceleration phase in batch culture). As was previously commented, only 56 genes appear to be differentially up-regulated at the stationary phase, when the growth of the strain has ceased.
Despite that most of the genes operating at the exponential phase persist further, the sets of differential genes detected in the two analyzed transitions seem to be characteristic of each step. In fact, only 10 out of 160 total up-regulated genes (i.e., 6.2%) were common between the glucose-limited and the stationary phase. These common genes are the following: c2666, yjfN, c1304, prpB, c1608, ycfQ, bsmA/yjfO, c5447, mocA/ygfJ, ydcS. There, the gene ydcS codifies a protein with activity of polyhydroxybutyrate (PHB) synthase, the bsmA/ycfO gene is related to peroxide resistance, biofilm stress, and motility, while the mocA/ygfJ gene codifies the CTP:molybdopterin cytidylyltransferase enzyme.
Concerning the 104 genes detected in the wild type, the transcription of which increased as dependent on the glucose limitation (see Table S1), 48 genes codify proteins with no annotated functions (hypothetical proteins). Moreover, two genes correspond to non-coding RNA (ryeA and cyaR). The cyaR gene is a direct target of CRP and, by its negative effect on the OmpX translation, indirectly enhances the expression of fimbria proteins. The expression of cyaR also negatively influences the translation of the protein codified by the nadE gene, which is involved in the NAD salvage pathway. Importantly, on the remaining 56 differential genes detected, the ontology analysis indicates that the only significantly enriched term is the one corresponding to propanoate metabolism (FDR = 1 × 10 −5 ). The genes involved in this pathway are induced by CRP and belong to the operon prpBCDE (prpB [methyl-isocitrate lyase], prpC [2-methylcitrate synthase], prpD [2-methylcitrate dehydratase], and prpE [propanoate-CoA ligase]). This finding is in agreement with the trend observed in another study that used RNA sequencing (Vital et al., 2015). Propanoate metabolism is related to the degradation routes of the amino acids cysteine, methionine, glycine, serine, threonine, and isoleucine (see Figure 2). FIGURE 2 | Propanoate metabolism and methylcitrate cycle mapping the sub-set of four enzymes codified by the genes dysregulated in the microarray data corresponding to the glucose-limited chemostat culture with the wild-type E. coli strain (first row, Experiment A, Table 2). The enzymes are indicated by their Enzyme Commission (EC) number, while the names of the genes that encode them are indicated in red. The diagram is a simplified adaptation based on KEGG Pathways (map 00640) after Ontology DAVID 6.8 software processing (FDR =1 × 10 −5 ).
The fourth noticeable point is that the number of up-regulated genes in the so-called "chemostat1-batch" comparison, which is based exclusively on the wild-type strain (Experiment A, Table 2), is as high as the number detected when both mutant strains are compared against the wild-type strain (Experiments B and C, Table 2).
In brief, while the third point above raises the question about the nature of the transcriptional changes implied in the physiology of E. coli when working under increasingly glucoselimited conditions, the fourth point raises the question as to whether the inferences from the expression profile of mutants can really be attributed to the mutations carried out by the strains analyzed, or whether they must be attributed to other causes. In the following, we will present our exploration to determine the extent to which the two questions are separable or interdependent.

Interaction Between Growth Rate and Transcription Profile
Although the use of mutant strains in Experiments C and D pursued the objective of enhancing the visualization of the acetylation-deacetylation regulation in E. coli, the data from Table 2 show that the growth rate leads by itself to substantial changes in the transcriptome, even when no planned mutations were present. In other words, the effects of the mutations overlap with the effects of the growth rate. In statistical terms, it is said that the two variables are "confounded." For this reason, and exploiting the implicit structure of a 2 × 2 factorial experimental design present in the data, we will proceed to determine rigorously which genes have an expression with significant dependence on the growth rate. To this end, we will quantify the so-called "interaction effect" between the growth rate (µ) and the type of mutation carried out by the strain (see Table 3).
Note that according to the concept of super-ratio , a positive interaction effect in Table 3 refers to those genes for which the change in the level of their transcription during the transition from the wild type to the indicated mutant in the glucose-limited chemostat (D = µ = 0.23 h −1 ) was significantly greater than the change in the transcription for the same type of transition but in batch culture (exponential phase, TABLE 3 | Number of genes differentially expressed according to the statistical criteria of interaction between the growth rate (µ) and the type of mutation carried out by the strain (wild-type and cobB or patZ).

Total 50 14
The variables are considered on two levels, high (+) or low (−), such that µ + = batch (exponential-phase, non-limited culture), µ − = chemostat (limited culture); strain + = wildtype, strain − = defective mutant. The values in parentheses indicate the associated false discovery rate (%FDR). For details see Tables S5 and S6. µ max = 0.62 h −1 ). Thus, from Table 3 we conclude that there is a synergy between glucose limitation and some acetylated (or succinylated) transcription factors (or regulators), because the transcription of 12 genes appears as increased in the mutant cobB (i.e., when deacetylase and desuccinylase enzyme CobB is absent).
Interestingly, from the 12 genes with positive interaction between µ and cobB, multiple functions are covered in a quite coherent mode, because all contribute to supporting better survival under glucose limitation. Thus, three genes are associated with improving the influx of nucleotides (codB, nupC, uraA) and two genes are related to the nucleotide balance through the salvage pyrimidine pathway (codA and ndk genes, which codify the enzymes cytosine deaminase and nucleoside kinase, respectively). Another important set of genes is given by the nuoM gene (codifying a protein involved in the generation of proton-motive force, a member of a terminal respiratory chain), the glnL gene (codifying the regulatory protein ntrB, a member of the two-component system that senses the availability of nitrogen), the puuB gene (codifying an oxidase involved in the degradation of the polyamine putrescine), and the mupC gene (which codifies a siderophore involved in the uptake and reduction of ferric ions).
Conversely, a negative interaction denotes that the variation in the expression during the transition from the wild type to the indicated mutant in batch culture (exponential phase, µ max = 0.62 h −1 ) is significantly greater than the one observed in the limited chemostat (D = µ = 0.23 h −1 ). Interestingly, in the case of negative interaction between µ and cobB, 24 out of 38 genes detected (63%) codify non-annotated products (i.e., hypothetical proteins). From the remaining genes, 5 out of 38 (13%) correspond to non-codifying RNA (micC, ryjA, sibd, ffs, and spf ), and 2 out of 38 (5.2%) genes detected correspond to opposing functions (fliR and yjfO, related respectively to the flagellar secretory protein system and to biofilm formation). More importantly, 4 out of 38 (10.5%) genes detected have very related functions, as they belong to a same operon, prpBCDE, which codifies enzymes corresponding to the propanoate cycle. Paradoxically, we have to conclude that the transcription of the set of genes with negative interaction appears as diminished under glucose limitation (or enhanced at the exponential phase), also due to the presence of some acetylated or succinylated transcription factor or regulator. In fact, these dual, opposing effects seem to be mediated by the same acetylation mechanism (this point is further analyzed in section Transcriptome Profiling of CRP Mutants).
Concerning the positive interaction between µ and patZ, only two genes were detected: one is pyrD (codifying the enzyme dihydroorotate dehydrogenase 2) and the other is fhuF (codifying a ferric reductase enzyme). By contrast, 12 genes were detected concerning the negative interaction between µ and patZ ( Table 3). In this case, 9 out of 12 (75%) genes detected still have no annotated function (hypothetical proteins), while 2 out of 12 (17%) genes correspond to non-coding mRNA (ryjA and ffs). In addition, the nrfG gene was detected, which codifies dual enzyme (formyl-dependent nitrite reductase activity, in anaerobic metabolism) and heme lyase from the nrf EFG operon (for the insertion of heme into c552 cytochrome).
The data in Table 3 suggest a low transcriptional responsivity to the deletion of the enzyme PatZ in relation to the observed with the deletion of the enzyme CobB. While the mutant defective in the deacetylase CobB yielded a total of 50 differential genes, the mutant defective in the acetyl-transferase PatZ only accounted for 14 differential genes (see Table 3). The prevalence of CobB upon PatZ effects (ratio = 3.6:1) indicates that there is no symmetrical effect, even though the two enzymes exert opposing actions. This fact suggests two alternatives: (i) it could be possible that other, as yet unidentified acetyl-transferases are present; (ii) more feasibly, the prevailing mechanism of protein acetylation does not seem to be dependent on acetyl-transferase PatZ, but mediated non-enzymatically by acetylphosphate (Kuhn et al., 2014;Wolfe, 2016).

Dissecting the Transcriptional Profiles in the cobB and patZ Strains
Another result arising from our analysis is the verification that transcription profiles kept a very strong correlation between the two mutated strains tested along the more than eight thousand genes examined (R 2 = 0.9854, see Figure 3A). Moreover, the regression over the cloud of points in Figure 3A has a linear trend with a unitary slope (m = 1.0113). This identity line is somewhat unexpected, because the mutant strains were designed to produce opposing effects on the acetylation-deacetylation phenomena, which is not evident from the monotonic, positive relation observed in Figure 3A. The expected relationship probably cannot be observed because it is masked as a consequence of the scant number of differential genes involved. So, the challenge is to identify a small number of altered genes by isolating them from a large series of genes affected by normal, noisy transcription.
The fact that the number of genes compromised by one or the other mutation is actually small in relation to the large proportion of genes that do not vary significantly suggests that data series corresponding to one of them, e.g., the mutant patZ, can also be used as background to determine the differential genes in both mutants if, in addition to itself, the other mutant were also considered under the form of a ratio of expressions (see Figure 3B).
In the case of Figure 3A, the conceptual problem of determining which points (i.e., genes) depart significantly from the slope value = 1 is analytically solved by Q-GDEMAR. Using this method, we can establish the profile of distribution of the ratio between the expressions of cobB and patZ mutants, while making the deconvolution of central data distribution (see column 1, Table 4). In the case of Figure 3B, based on the similarity of this representation with Shewart and Deming's control charts (Tague, 2005), the dysregulated genes can be determined by establishing which points of the graph exceed the control limits (see column 2, Table 4). Importantly, whichever The data came from the chemostat culture mode. The values in parentheses indicate the false discovery rate (%FDR). For details see Table S7 and Table S8. the type of discrimination procedure adopted, given that both data sets of the mutants came from chemostats operated at the same dilution rate, the differential genes are identified herein on the basis of equal influence of the growth rate, and without using information corresponding to the wild-type strain. This point is another important difference with respect to previous analyses done by Castaño-Cerezo et al. (2014). As expected, the analysis of the cobB: patZ ratio by the deconvolution of central data of the distribution renders a much more sensitive detection of genes associated with the mutations tested. In fact, by applying Q-GDEMAR we were able to identify 94 up-regulated genes and 87 down-regulated genes, whereas the second procedure only detected 22 and 11 genes, respectively (see Table 4).
However, by comparing the up-regulated genes from column 1 in Table 4 with the up-regulated genes from column 1 in Table 2, we verify that all 94 genes from Table 4 (based only on the mutants) are actually included in the list of 104 genes previously detected in Table 2 (based only on the wild-type strain). Noticeably, this means that 100% of the genes identified as differentially up-regulated in column 1 of Table 4 cannot be explained by the mutation cobB, as is supposed, but due to the low growth rate generated by the substrate limitation in the culture. Instead, from the 87 down-regulated genes detected in column 1 of Table 4, 96.5% (84 genes) differ from the ones detected in column 1 of Table 2 with the wild-type strain, and hence, they can be reliably attributed to an increased transcription under the mutation patZ.
By applying the second strategy (confidence interval to the fluctuation of the cobB: patZ following the form of representation shown in Figure 3B), we recover those genes that depart significantly from the horizontal reference line (y = 1). In this way, as shown in column 2 of Table 4, we were able to identify 22 genes (FDR = 1.08%), whose transcription is apparently associated with the loss of the gene cobB, while 11 genes (FDR = 2.18%) showed a transcription that is apparently associated with the loss of the gene patZ (see Table S5, Supplementary Materials).
Importantly, the transcriptional effects associated with cobB mutation that could not be previously identified from column 1 of Table 3 can now be identified from column 2 of Table 3. In fact, only 3 out of 22 genes (13.6%) detected in column 2 of Table 3 are shared with the list of up-regulated genes from column 1 of Table 1. That is, following the analysis over the distribution of the second ratio form in Table 3, it has been possible to identify 19 genes whose expression is enhanced by the loss of CobB activity. Interestingly, none of the 11 genes detected as down-regulated in column 2 of Table 3 are common to the list of 87 genes detected as down-regulated in column 1 of Table 3. Consequently, they have to be added to the list of 84 genes whose expression is enhanced by the loss of PatZ, thus bringing the total to 95 genes.

Transcriptome Profiling of CRP Mutants
When we studied the CRP mutants, an additional problem arose, because the mutated genes cannot be cloned at their native position and the genomic insertion had to be carried out into the paaH locus. For this reason, the mutant N was designed, which has to be considered as an additional control besides the wild-type strain (see Materials and Methods). In fact, the change of position in the cloning of native crp gene led to a considerable distortion, as 281 genes ended up being spuriously transcribed (see Table 5).
Note in Table 5 that when the values of WT-N depart significantly from zero, this means that the transcriptional response of the N strain differs from the wild-type strain. Accordingly, when WT-N << 0, this implies an increased transcription of some genes in the mutant N with respect to what is observed in the wild-type strain. By contrast, when WT-N >> 0, this implies a significant attenuation of transcription of some genes in the N strain with respect to the wild-type strain. In spite of the significant number of over-and under-transcribed genes that appeared in Table 5, the validity of the experiments done with the crp mutants is still guaranteed. In fact, we have verified that the signal ratio ( N stationary / N exponential ) kept a high correlation with the homolog ratio (WT stationary /WT exponential ) across the 8,862 genes tested (R 2 = 0.8241, see Figure S1).
Thus, once we had identified the genes involved in Table 5 (see Table S9), appropriate corrections were applied to the list of differential genes corresponding to the CRP mutants. This was done by crossing the cases of "false negatives" from Table 5 against the list of total differentially down-regulated genes, and the cases of "false positives" in Table 5 against the list of total differentially up-regulated genes. In turn, the comparisons were disaggregated according to the phase of culture (exponential or stationary) and the form of crp considered (Q-type or R-type). The net results after these corrections are shown in Table 6.
Although we tested different types of discrimination variables to identify the transcriptional effect of CRP mutants, the best form that we found is the one applied in Table 6. There, the discriminating variable is given by the ratio between the signals from the analyzed mutant with respect to the ones in the wild type for each culture phase. This implies considering the following set of partial ratios: (R/WT) stationary >>1 and (R/WT) stationary <<1, (R/WT) exponential >>1 and (R/WT) exponential <<1, (Q/WT) stationary >>1 and (Q/WT) stationary <<1, (Q/WT) exponential >>1 and (Q/WT) exponential <<1.   The values in parentheses indicate the associated false discovery rate (%FDR). *NDAS, not detectable as significant.
The operation of correction applied to achieve the results shown in Table 6 proved particularly important for the data corresponding to the exponential phase. There, we had to exclude 46.9% of up-regulated genes for crp(Q-type), 15.3% of upregulated genes for crp(R-type), and 13% of down-regulated genes for crp(Q-type). However, only minor corrections were necessary for the remaining cases: 5.6% of down-regulated genes for the crp(R-type) under exponential phase and 1.3% of downregulated genes for crp(Q-type) under stationary phase. No correction was needed for the up-regulated genes for crp(Qtype) under stationary culture phase. The detail of the genes supporting Table 6 is given in Tables S10-S13.
Importantly, based on data in Table 6 we were able to establish which genes are influenced by each type of mutant in each growth phase (see Figure 4).
From Figure 4 we observe several noticeable findings: (i) CRP is shown to act not only at the stationary phase when glucose is exhausted (bottom of Figure 4), but also during the exponential phase when glucose availability is unconstrained (top of Figure 4); (ii) For a given growth culture phase, both deacetylated and acetylated stages of CRP are shown to influence "exclusive" sets of genes (for the genes corresponding to blue and pink areas in Figure 4, see Tables S10-S13); (iii) However, in three out of four cases analyzed, each of the crp mutants (Q-type and R-type) is shown to be able to produce "exclusive" effects in opposing directions (for a given culture phase in the Venn Diagram, compare areas from the left where genes are up-regulated with the homologs at the right where genes are down-regulated). Therefore, from the type of analysis done, we have to conclude that the acetylation stage of CRP by itself does not determine univocally the type of transcriptional effect produced; (iv) Importantly, at the exponential phase we detected two sub-sets where the genes seem be equally responsive to crp(Q-type) and crp(R-type) (see green areas of intersection, Figure 4). The identity of these common genes changes whether the genes are up-or down-regulated (see Tables S10, S11); (v) According to the number of genes affected, the transcriptional effects of CRP deacteylated prevail over CRP acetylated during the exponential phase (the ratio crp(R-type): crp(Q-type) ranged between 7:1 and 2.5:1, whether genes are up-or down-regulated, respectively). This feature is even more manifest at the upregulated genes of the stationary phase, where the same ratio reached a value of 13:1; (vi) However, only the crp(Q-type) was detected to exert inhibitory transcriptional effects in the stationary phase (see Figure 4).

CRP Effects Independent of Deacetylation/Acetylation
The evidence that both mutants, crp(Q-type) and crp(Rtype), shared some common target genes (13.3% in the upregulated genes, 16% in the down-regulated genes, Figure 4) during the exponential phase was somewhat unexpected. A first, simple interpretation of this phenomenon could suggest that transcription of these sets of genes does not depend on the acetylation stage of CRP, but rather on the presence of several conformer complexes between c-AMP and CRP (Heyduk and Lee, 1989;Mukhopadhyay et al., 1999;Tworzydło et al., 2005;Tutar, 2008;Saha et al., 2015).
Supporting the "conformers" hypothesis, it is noticeable that neither the crp gene nor the cyaA gene was responsive to the acetylation stage of the CRP protein. In fact, the microarray probes corresponding to these genes displayed similar signal intensity under all the scenarios tested (data not shown). This point is important, since the CRP protein can influence its own transcription from the crp gene, while the cyaA gene is a direct CRP target that codifies de enzyme adenylate cyclase. This enzyme is required for the biosynthesis of cAMP, a metabolite that modifies allosterically the CRP protein such that CRP can behave effectively as a master regulator (Zheng et al., 2004;Zhou et al., 2014).
Thus, the above observations raise the question as to whether the regulation by acetylation, although conserved from bacteria to eukaryotes, actually has some form of hierarchical order. In this regard, the sub-sets of genes that are influenced by both CRP deacetylated and CRP acetylated (green areas, Figure 4) could represent distinct adaptations along the evolution and/or could have a distinct functional value than those genes preferentially affected by one or the other form (blue and pink areas, Figure 4). In any case, it is clear that crp and cyaA genes, due to their central role, are at the top of the CRP-dependent regulation cascade (see Figure 6). Other possible explanations for the shared genes between acetylated and deacetylated stages of CRP will be presented in the next section.

Ontology Analysis of CRP Mutants
Tables S10-S13 provide exhaustive lists of the differential genes transcribed by the CRP mutants tested, whilst a global summary of this information is shown in Figure 3. By means of ontology analysis, we attempted to determine whether the differential genes identified can be associated with some specific bacterial function in terms of well-standardized, significantly over-represented categories (Rhee et al., 2008;Lin and He, 2012). For the ontology analysis we also adopted as cut-off an FDR value ≤5%. The results of these analyses are summarized in Tables S14-S18.
Some of the well-defined functional groups that we identified are associated with the stationary phase and are driven mainly by the crp(Q-type) mutant. There, within the up-regulated genes, we detected a set of five enzymes related to the metabolism FIGURE 4 | Venn diagram showing the disaggregation of the differential genes corresponding to each of the CRP mutants (Q-type: acetylated form, R-type: deacetylated form) according to the batch culture phase (exponential or stationary). The blue areas in the diagrams correspond to genes that depend exclusively on the Q-type, while the pink areas correspond to the genes that depend exclusively on the R-type. Importantly, the green regions correspond to common genes that can be influenced by both types of CRP mutants.
of Arginine (FDR = 1.8 × 10 −6 , see Table S16), whereas within the down-regulated genes, we identified a set of eleven ribosomal proteins (FDR = 1.2 × 10 −8 ), a set of eleven enzymes related to the pyrimidine metabolism (FDR = 4.9 × 10 −5 ), and a set containing seven H + -ATPases, which are related to the generation of energy in the oxidative phosphorylation (FDR = 1.1 × 10 −4 ) (see Table S17). Interestingly, the same set of five enzymes related with catabolism of arginine detected for the stationary phase of batch culture based on the crp mutants (Table S16) also appeared as significant in the ontology analysis of the stationary phase with the wild-type strain (experiment B in Table 2, ratio >>1).
Although other functional classes containing multiple related genes have also been recognized, including some associated with the crp(R-type) mutant, these functional groups did not satisfy the FDR criterion (see Tables S14, S15). Given that in some of the conditions assayed, we identified 400-500 differential genes (see Figure 4), the low number of genes associated with the ontology classes and the low number of functional categories recognized as significant seems be a challenging problem.
However, this ontology picture can be explained by several factors: (a) CRP is a "master regulator." As such, CRP interacts with 47 local co-regulators and controls the induction of 22 transcription factors (Martínez-Antonio and Collado-Vides, 2003;Martínez-Antonio et al., 2008). This complexity allows CRP to influence up to 283 operons (Ishihama et al., 2016). Hence, the low detection of significant ontological classes could be due to the very high number of pleiotropic effects associated with CRP. In fact, the numerous effects of CRP might appear as "distributed" across multiple functional categories, rather than concentrated. As a consequence, the probability of detecting a given effect will be strongly diminished; (b) Results from the ontology analysis are sensitive to the low number of genes sampled. Although the lists of differential genes sampled contained a large enough number, in some of the experimental groups analyzed around 50% of the genes fall under the category of "hypothetical genes." Since these genes are classified as "unidentified" or "unknown, " they are not accounted for in the computation of the hypergeometric probabilities, thus shortening the list of available genes; (c) For the crp mutants at the exponential phase, the correction applied to the data, based on the (positional) control strain N, proved to be excessively protective (47% of differential genes of the crp(Qtype) mutant were excluded, see section Transcriptome Profiling of CRP Mutants). In fact, the criterion of exclusion used was based on comparing the degree of over-expression (or underexpression) in the N strain with respect to that observed in the wild-type (WT) strain rather than on considering qualitatively if induction (or repression) occurred.
To overcome this problem, alternatively we only excluded the differential genes related to the paaABCDEFGHIJK operon in which the crp gene was cloned. In this way, the number of up-regulated genes at the exponential phase, driven by crp(Qtype) or crp(R-type) indistinctly, increased from 70 to 104. More importantly, from this extended series we were able to verify the significant over-representation of 32 genes associated with the flagellar assembly (FDR = 7.2 × 10 −44 ) together with another 15 genes associated with chemotaxis (FDR = 7.5 × 10 −15 ) (see Table S18). Interestingly, most of the genes in this highly significant set are targets of the products codified by the genes flhC and flhD, which are responsive to the cAMP-CRP complex (see Figures 5, 6).
Reasoning in inverse form, the high significance detected herein for the flagellar and chemotaxis ontologies (FDR = 7.2 × 10 −44 ) offers strong support for the existence of sets of genes with the apparent property of being driven indistinctly by crp(Q-type) or crp(R-type). In fact, the functional coherence of this group and the extremely low value of FDR ruled out the possibility of a methodological artifice or a grouping generated by chance.
However, several questions still emerge: (i) Given that a subset of genes appeared equally transcribed by the two opposing crp mutants, one could wonder about the relevance of the acetylation stage of CRP for the transcription of these genes; (ii) In the contrary, assuming differences in the transcriptional effects associated to each crp mutant: Might it be possible that the CRP acetylation variants acting at different levels of the CRP regulatory network lead to a situation in which the global effect appear "as if " no differences exist? (green areas in Figure 4); (iii) How are explained the "exclusive" differential effects also observed for each of the acetylated CRP stages? (blue and pink areas in Figure 4); (iv) Finally, it should be explained also how a given acetylated (or deacetylated) form of CRP can exert both positive and negative effects within a same growth phase (for a given row, areas with the same color between the left and the right diagrams in Figure 4). To address these questions we have compiled the information concerning the genes involved in Figure 6.
In this regard, it is known that the crp gene is regulated positively by the transcription factor CRA (FruR), but only if the level of fructose 1,6-biphosphate is low. The crp gene has a complex relationship with its negative regulator FIS, since CRP induces the fis gene when FIS is absent, but represses this gene when both are present together (Nasser et al., 2001). In addition, the crp gene has two binding sites for the CRP protein, a stimulatory one located upstream of the promoter, and an inhibitory one located downstream of the promoter (Ishizuka  ) or negative (-|) effects. The chemical intermediaries responsible of the effects among the genes (transcription factors, cAMP, ppGpp) are omitted. The nodes colored in green define a virtual output layer with the flagellar and chemotaxis genes. The yellow node, the operon flhD/C, codifies a master regulator with effects on several flagellar genes as well as on other important extra-flagellar genes (not shown; (Prüss et al., 2003)). The nodes colored in pink correspond to a highly regulated sub-network that include the genes codifying CRP, several sigma factors, and another interdependent regulators. The resulting signal from this complex layer acts as input toward the operon flhD/C and the remaining output genes. The diagram was built by composing the results showed in Figure 5 with the regulatory information given at EcoCyc database. et al., 1994). CRP is not really what is bound, but rather the cAMP-CRP complex (Saha et al., 2015) (see Figure 6).
In turn, the availability of cAMP is a complex function of the cyaA gene activity and the energy status of the cell (levels of ATP, GTP, AMP). However, the cyaA gene is negatively regulated by the cAMP-CRP complex, while ATP can also be diverted toward the synthesis of the alarmone ppGpp mediated by the pyrophosphokinase RelA, but whose activity depends on the nutritional and stress stage of the bacteria (Magnusson et al., 2005;Burgos et al., 2017) (see Figure 6).
On the other hand, the assembly of flagella and the synthesis of components of the chemotaxis system during the exponential growth phase require the simultaneous transcription of a plethora of genes, distributed along non-contiguous operons. However, as is shown in Figure 6, most of these genes are not direct targets of CRP, but are induced by the molecules codified by the genes flhD and fllhC. Importantly, both genes belong to the same operon and are direct targets of the cAMP-CRP complex (see Figure 6).
Note that the transcriptional processes related to the flagellar biogenesis and chemotaxis are mainly mediated by RNA polymerases containing sigma factors σ 28 , σ 54 , and σ 70 . However, the transcription of these sigma factors is regulated by several repressors (ppGpp-DKSA complex, CSGD protein) and the protein CPXR (a repressor of the csgD gene). All these repressors have the notable particularity that they are induced by the cAMP-CRP complex (see Figure 6). Moreover, the IHFA-IHFB complex, which represses the critical genes flhC and flhD, has a dual effect on the genes fliA and rpoF, which codify the sigma factors σ 28 and σ 70 , respectively. In fact, the IHFA-IHFB complex acts as an inducer of CSGD -a sigma factor repressor-, while simultaneously contributing to the induction of these sigma factors. In the same line, the IHFA-IHFB complex also leads to the induction of the factor anti-sigma 28 from the gene flgM (see Figure 6).
In brief, given the complex structure of the regulatory network involved in flagellar assembly, it is highly plausible that if a given acetylated form of CRP were able to influence a gene codifying a repressor negatively, while the other acetylated stage were able to influence the target gene of this repressor positively, then the gene in question could appear up-regulated as being modulated indistinctly by both opposing CRP species. Similarly, a gene could appear as down-regulated indistinctly by both opposing CRP species if one acetylated CRP stage acts positively on a gene codifying a repressor, while the other acetylated CRP stage acts negatively on the target gene of this repressor.
To illustrate this hypothesis, based on the network showed in Figure 6, the following chains of process leading to the synthesis of the FLHD-FLHC complex can be identified: CRP dksa ihfBA fis hns flhD/C CRP relA ihfBA fis hns flhD/C From these set of processes it turns evident that if it were required that the operon flhD/C be maximally transcribed, CRP in Equation (1) might be in a one of the CRP acetylation stages (the able to activate the operon flhD/C), whereas in Equations (2)-(4) CRP could be present in the contrary stage of acetylation (thus inhibiting the genes ompr, dksa, relA). This distribution of CRP variants would provide the maximal induction of the flhD/C operon, thus showing how is feasible that a gene appears from a phenomenological point of view as if they were equally influenced by CRP acetylated and CRP deacetylated even when this was not the real mechanism. This hypothesis to be satisfied requires that both acetylated stages, CRP acetylated and CRP deacetylated , coexist in some proportion. Moreover, both CRP acetylated and CRP deacetylated could have similar affinities for their DNA-binding sites but would differ in their mode to interact with the RNA polymerase. The first requirement does not seem restrictive since it depends only on the balance between the enzyme CobB (a deacetylase) with respect to the concentration of acetyl-phosphate (the main no-enzymatic donor of acetyl moieties, related with the metabolism of acetate). Consequently, the flhD/C operon could also be transcribed in variable extents under distinct biological scenarios according to the increased or decreased balance between these factors.
The fact that CRP acetylated and CRP deacetylated also showed to exert specific effects (see Figure 4) can be explained by assuming that each of CRP variant acts on a given target gene but not on the genes which codify a repressors of these target genes. The case that a given form of CRP can exert both, positive and negative effects on different genes is difficult to justify; it could require a better understanding of the interactions between CRP and other transcriptional regulators with action on these genes.

CONCLUDING REMARKS
The present study aimed to characterize the regulation by the acetylated/deacetylated forms of CRP in Escherichia coli. In this context, we also analyzed their dependence on their most direct regulator (glucose) and some direct effectors (deacetylase CobB, peptidyl N(ε)-Lysine acetyl-transferase PatZ). Although this topic belongs to the ambit of molecular physiology, the clarification of these mechanisms could have important consequences in the field of biotechnological applications (Kuczynska-Wiśnik et al., 2016;Ahsan et al., 2017;Basak et al., 2017;Ishigaki et al., 2017;Venkat et al., 2018), and even for some medical issues (Bernal et al., 2016;Ou et al., 2017).
Our first reflection is that effects of CRP are currently associated with the stationary phase of bacterial cultures. There, the catabolic repression exerted by CRP is relieved, thus enabling the activation of genes for alternative substrate assimilation and a better adaptation to this stage (Deutscher, 2008;Schilling et al., 2015). Apart from the catabolic repression phenomenon, we found few studies examining other effects of CRP in growing bacteria (Gomez-Gomez et al., 1996;González-Flecha and Demple, 1997;Shimada et al., 2013;Baron and Eisenbach, 2017;Zhou et al., 2017). A common feature of these studies is that they focused on a single gene or protein. Our analysis, based on the systems biology approach, lead us to analyse the CRP effects on the genome-wide scale with the intention of not only capturing the trends emerging from the data but also, when feasible, inferring the more general organizing principles behind the CRP regulation. At this genome-wide level, we were able to verify that CRP acts along all the culture phases (i.e., with excess, limitation, or shortage of glucose). However, what was unexpected was our observation that successive transitions toward the stages of increasing glucose limitations appear to be associated with the up-regulation of specific sets of target genes rather than with the loss of genes present when glucose was in excess ( Table 2).
One could conjecture at least two possible explanations: (i) most of the genes are essential to survival, with only a minor fraction responsible for the adaptation; (ii) the genes acquired during the glucose limitation could be an adaptive response to the stress generated by the carbon-source shortage, which would coexist with the pre-existing machinery if stress intensity is not extremely high. However, the first hypothesis must be ruled out, since <10% of the Escherichia coli genome is essential for its survival (Mahadevan and Lovley, 2008). Concerning the second hypothesis, it is well-established that the transcriptional response to "nutrient limitation" differs from the response to "deprivation" (Saldanha et al., 2004). Moreover, in the glucoselimited chemostat, we did not observe any dysregulation in the stress markers (e.g., relA, rpoS, rpoE, rpoN, crp, arcA, arcB, katE, sodA, cpxa, cpxr), nor did the transcription of genes codifying some enzymes of central metabolic pathways appear to be dysregulated (e.g., pfka, aceA, aceE, aceF, aceK, sdhA, sdhB, sdhC, sucD, fumB, pykA, pykF, glpc, pta).
Similar observations have been made for other biological systems. Thus, important variations have been verified in the fluxes that these enzymes catalyze, leading researchers to conclude that regulation under glucose-limiting conditions exists, but must be mainly of a post-transcriptional nature (Daran-Lapujade et al., 2004). In our case, these reports lend support to the possible importance of acetylationdeacetylation as the mechanisms responsible for the observed regulation. In fact, although the effects of glucose limitation overlap with the effects caused by the mutants cobB and patZ, by introducing Experiment A ( Table 2), analyzing one mutant against the other, and applying the factorial analysis, we were ultimately able to identify the genes that really participate in the interaction among these variables, thus offering direct proof that regulation by acetylation/deacetylation is indeed operating (sections Substrate Availability and Genome-Wide Transcription, Interaction Between Growth Rate and Transcription Profile, Dissecting the Transcriptional Profiles in the cobB and patZ Strains). Moreover, from the asymmetry observed in Table 4 between the total number of genes affected by the mutant cobB and the ones in the case of the mutant patZ (ratio 2.6:1), it was inferred that some as yet unidentified acetyl-transferase must exist, or that a non-enzymatic mechanism (e.g., acetyl-phosphate) might be the dominant mechanism with respect to the mediated by the enzyme PatZ. This investigation highlight some aspects that requires further inquiry. This is the case of the propanoate pathway in the wild-type strain under glucose limitation, the observed dysregulation of non-coding RNAs genes under negative interaction between the growth rate and the cobB mutant (genes micC, sibD, ryjA, and spf ), and the dysregulation of genes codifying small RNAs under the negative interaction between growth rate and the patZ mutant (genes ryjA and ffs) (see section Interaction Between Growth Rate and Transcription Profile). Note that the non-coding RNA MicC, by modulating OmpC, can inhibit the bacterial motility, while the small RNA Spf inhibits the formation of the complex succinate:quinone oxide reductase from the operon sdhCDAB. The Ffs small RNA, together with the Ffh peptide, constitute the signal recognition particle (SRP), which directs the insertion of the nascent proteins into the plasma membrane.
As a second comment, we conclude that the acetylation of CRP has proved to be a no ordinary form of regulation. In fact, both CRP acetylated and CRP deacetylated forms can produce positive and negative effects on "exclusive" sets of genes. The number and identity of these genes is conditioned by the glucose availability (section Transcriptome Profiling of CRP Mutants). Specific ontology functional classes appeared to be associated with each of these different sets of genes (section Ontology Analysis of CRP Mutants).
Apart from the exclusive differential effects recognized for each deacetylated/acetylated stage of CRP in a given culture phase, we also identified groups of genes that, being directly or indirectly responsive to CRP, seem not to depend on the CRP acetylation stage. For these genes, susceptible to "shared" effects, we proposed two possible explanations. One is that their regulation does not depend on the acetylation stage of CRP but rather on the conformational stage of the cAMP-CRP complex (section CRP Effects Independent of Deacetylation/Acetylation). Alternatively, these groups of genes might depend on each acetylation stage of CRP, but due to the nested characteristics of the network involved, phenomenologically the regulation occurs as if the observable mechanisms were independent of the acetylation forms (section Ontology Analysis of CRP Mutants). Although in sections Transcriptome Profiling of CRP Mutants, CRP Effects Independent of Deacetylation/Acetylation, Ontology Analysis of CRP Mutants we presented the arguments supporting each of these hypotheses, additional experiments will be required to rule out one or the other. In any case, from the ontology analysis we determined that the assembly of flagella and chemotaxis are the most important functions that, being regulated by the cAMP-CRP, appeared as being equally sensitive to both CRP acetylated and CRP deacetylated (FDR = 7.2 × 10 −44 ) (section Ontology Analysis of CRP Mutants). Our finding of an enhanced transcription of the genes corresponding to the flagellar assembly and chemotaxis during the exponential phase is in agreement with some published reports (Amsler et al., 1993;Sim et al., 2017), in contrast with the generally accepted conception that flagellar activity is associated to the stationary phase.
Recently, it has been suggested that the crp(Q-type) mutant (i.e., the CRP acetylated ) exerts an inhibitory effect on the CRPdriven genes containing the type II promoters, while possibly also simultaneously exerting a positive effect on the CRP-driven genes containing the type I promoters (Davis et al., 2018). However, the hypothesis of the opposed effect of CRP acetylated depending on the type of promoter still remains to be experimentally verified. The most relevant conclusions of this study arose from an E. coli in which the mutated genes crp were cloned in the chromosomal gene paaH and challenged with glucose and acetate. In our investigation we have data concerning to this genetically modified strain with glucose as source of carbon (see section Strains and Conditions Tested). We have considered the data referred to the crp mutants with glucose because we verified that the strain N (the cloning control) exhibits an acceptable transcriptional behavior when compared with the wild-type strain (R 2 = 0.8241, see section Transcriptome Profiling of CRP Mutants). However, the same strain N showed to be no-appropriated when is challenged with acetate. In this case, we observed that the transcriptional response of the strain N differs markedly from the one observed in the wild type strain when tested along 8,862 genes (R 2 = 0.4406, see Figure S1). In any case, beyond of the critical experimental details mentioned, the theory presented by Davis et al. (2018) does not seem to assign any role to the crp(R-type) mutant (i.e., the CRP deacetylated ) for which we have also observed specific, significant effects. Neither does the proposal account for the group of genes that we have herein detected as being equally influenced by both CRP acetylated and CRP deacetylated .
Other researchers, using advanced computational techniques, have built a predictive model of CRP based on the DNA sequences of its binding sites in the E. coli genome (Tsai et al., 2018). This model not only accounts for the induction effects of CRP-responsive genes with promoters of the type I and type II, but also predicts the repression effects of CRP. It could seems then that acetylation/deacetylation of CRP has not any regulatory influence. However, despite the importance of this model, since it is based uniquely on descriptors derived from the sequence of DNA, it cannot deal with questions such as the physiological conditions under which CRP will exert a positive or a negative effect. To be able to address this question it would be required to know the factors beyond of the DNA sequence that determine the binding of CRP to a positive or negative site in a given gene. In turn, this requires to establish if the binding of CRP to a positive or negative sequence in the specific DNA sites studied is a stochastic phenomenon or it depends on some particular form of post-translational modifications in the CRP or in another associated transcriptional regulators (e.g., acetylated/deacetylated; Brown et al., 2017). This question is relevant because we have observed that there is not a simple binary correspondence between acetylated/deacetylated form of CRP and their positive or negative effects. Both forms of CRP can produce both positive as negative effects during the exponential phase as well as at the stationary phase (see Figure 4).
Finally, with a view to future experiments, we have to mention the need to consider the succinylation phenomenon, since the deacetylase enzyme CobB also exerts a desuccinylase action (Colak et al., 2013), while succinylation as a post-translational modification overlaps with many sites of peptidyl N(ε)-Lysine acetylation (Weinert et al., 2013b;Okanishi et al., 2017). Hence, this information could help to clarify or even to modify our comprehension of regulation by acetylation/deacetylation in E. coli.

AUTHOR CONTRIBUTIONS
DG designed and performed the analysis of the data. NT supervised the data analysis. DG and NT wrote the manuscript.

FUNDING
This work was supported by the research projects BIO2014-54411-C2-2-R from the Ministerio de Economía y Competitividad (MINECO, Spain).

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb. 2018.00941/full#supplementary-material Figure S1 | Comparison between the overall transcriptional responses of the N strain (cloning control) with respect to the one in the wild-type strain whether the type of substrate used (glucose or acetate).
Table S1 | Differential genes transcribed in the wild-type strain of E. coli under glucose-limited culture (chemostat, D = 0.3 h −1 ) with respect to the one growing in glucose excess (exponential phase, µ max = 0.62 h −1 in batch). Table S2 | Differential genes transcribed in the wild-type strain of E. coli in stationary phase (batch culture) with respect to the one growing in glucose excess (exponential phase, µ max = 0.62 h −1 in batch). Table S3 | Differential genes transcribed in E. coli carrying defective mutation cobB with respect to the wild-type strain, both growing in chemostat under glucose-limited conditions (D = 0.3 h −1 ). Table S4 | Differential genes transcribed in E. coli carrying defective mutation patZ with respect to the wild-type strain, both growing in chemostat under glucose-limited conditions (D = 0.3 h −1 ). Table S5 | Differential genes transcribed in E. coli according to the interaction in the factorial 2 × 2, where the factors analyzed are the specific growth rate (µ high due to glucose excess in exponential phase or µ low due to the glucose-limited conditions in chemostat) and the mutation cobB (absence or presence). Table S6 | Differential genes transcribed in E. coli according to the interaction in the factorial 2 × 2, where the factors analyzed are the specific growth rate (µ high due to glucose excess in exponential phase or µ low due to the glucose-limited conditions in chemostat) and the mutation patZ (absence or presence). Table S7 | Differential genes using as discrimination variable the ratio between the expression in the two mutants ( cobB: patZ) when applying the confidence limits arising from the Gaussian central region deconvolution (Q-GDEMAR). Table S8 | Differential genes using as discrimination variable the ratio between the expression in the two mutants ( cobB: patZ) when applying the confidence limits on the whole distribution without Gaussian central region deconvolution. Table S9 | Differential genes spuriously transcribed when comparing the cloning-control strain ( N) against the wild-type strain (WT), by applying Q-GDEMAR to the difference (WT-N) as discrimination criterion for the exponential phase and stationary phase of batch culture.
Table S10 | Differential genes transcribed in E. coli carrying crp(Q-type) mutation during exponential growth phase with and without correction by N.
Table S11 | Differential genes transcribed in E. coli carrying crp(R-type) mutation during exponential growth phase with and without correction by N.
Table S12 | Differential genes transcribed in E. coli carrying crp(Q-type) mutation during stationary growth phase with and without correction by N.
Table S13 | Differential genes transcribed in E. coli carrying crp(R-type) mutation during stationary growth phase with and without correction by N.
Table S14 | Ontology analysis for the crp(Q-type) and crp(R-type) mutants in the case of up-regulated genes during the exponential growth phase.
Table S15 | Ontology analysis for the crp(Q-type) and crp(R-type) mutants in the case of down-regulated genes during the exponential growth phase.
Table S16 | Ontology analysis for the crp(Q-type) and crp(R-type) mutants in the case of up-regulated genes during the stationary phase.
Table S17 | Ontology analysis for the crp(Q-type) mutant in the case of down-regulated genes during the stationary phase.
Table S18 | Detail of the genes up-regulated during the exponential growth phase that were induced by both the crp(Q-type) and the crp(R-type) mutants, being enriched in the flagellar assembly and chemotaxis ontology functions.