The third approach, implemented here, is to consider a family of transcription factors (TFs) and then identify binding motifs for as many TFs as possible. This is mainly motivated by the desire to analyze the structure of protein–DNA interactions and the co-evolution of TFs and the motifs they recognize (Desai et al., 2009; Huang et al., 2009; Camas et al., 2010; Leyn et al., 2011; Ravcheev et al., 2012). An important issue in such studies is to connect TFs to the cognate TF binding sites (TFBSs) identified by phylogenetic footprinting and other computational techniques (Conlan et al., 2005; Wels et al., 2006; Liu et al., 2008). This problem is either solved experimentally or addressed computationally, for instance for regulons controlled by local TF from specific protein families (Rigali et al., 2004; Francke et al., 2008; Sahota and Stormo, 2010; Ahn et al., 2012; Kazakov et al., 2013). Phylogenetic profiling of TF genes and motifs upstream of candidate regulon members is an alternative bioinformatics approach for assigning TFs to putative regulons (Rodionov and Gelfand, 2005). The comparative analysis of ligand-binding domains in TFs also helps identify ligand specificity determinants and propose models of functional diversification within large and functionally heterogeneous families of TFs (Kazanov et al., 2013).

Here we study the LacI family of bacterial transcription factors (LacI-TFs). The namesake of the family, the lactose repressor LacI of *E. coli,* has been the model object for the analysis of bacterial transcriptional regulation since the classical papers of Jacob and Monod (1959, 1961). The family was established by the analysis of similarity of protein sequences, and simultaneously the similarity of DNA motifs recognized by the family members was noted (Weickert and Adhya, 1992). At the same time, it was observed that its DNA-binding domains of LacIfamily regulators are similar to the helix-turn-helix domains of other TFs (Nguyen and Saier, 1995), whereas the ligand-binding domain of LacI-TFs is homologous to the periplasmic proteins of ABC-transporters (Mauzy and Hermodson, 1992; Fukami-Kobayashi et al., 2003). Interestingly, this domain was also seen in combination with another DNA-binding domain, winged helixturn-helix of the GntR family (Franco et al., 2006). While early observations on a limited number of sequences suggested that the history of the family involved a series of duplications at the early stage and a low level of duplications later on (Nguyen and Saier, 1995), more data, derived from bacterial genomes becoming available, demonstrated that the duplications were occurring throughout the history of the family (Fukami-Kobayashi et al., 2003). Given its large size and high level of structural similarity, yielding reliable multiple alignments, this family was widely used as a model for algorithms for identification of functionally important residues both uniformly conserved and specificitydetermining (Mirny and Gelfand, 2002; Fukami-Kobayashi et al., 2003; Kalinina et al., 2004; Pei et al., 2006; Tungtur et al., 2007; Parente and Swint-Kruse, 2013). Some of these predictions were further tested in experiment (Meinhardt and Swint-Kruse, 2008; Camas et al., 2010; Tungtur et al., 2010, 2011).

The structural similarity of the DNA motifs recognized by the LacI-TFs was used to identify regulons in *Lactobacillus plantarum* (Francke et al., 2008) and *Dickeya dadantii* [*Erwinia chrysanthemi*] (Van Gijsegem et al., 2008). The structural aspects of interactions of the LacI-TFs with DNA and ligands and between subunits in dimers have been summarized in a recent review (Swint-Kruse and Matthews, 2009).

Here we report results of a large-scale, manual comparative genomics analysis of LacI-TFs aimed at identification of their binding sites, motifs, and regulons in bacterial genomes. By analyzing the genomic and metabolic context of the reconstructed regulons and by combining this analysis with information gathered from literature, we inferred the biological roles and molecular effectors for a large number of the studied LacI-TFs. As result, we made a number of observations on the distribution of orthologous LacI-TFs in genomes, the statistics of the binding sites' arrangement in regulatory regions, and the number and functional characteristics of the regulated genes. By combining the functional annotations with phylogenetic analysis, we proposed evolutionary models of functional diversification for a number of LacI-TF groups. The obtained reference dataset of 1281 regulons in 272 genomes was deposited in the RegPrecise database (Novichkov et al., 2013).

# **MATERIALS AND METHODS**

The genomes were downloaded from the MicrobesOnline database (Dehal et al., 2010). TFs from the LacI family were identified by similarity searches and domain predictions in the Pfam database (Finn et al., 2014). LacI-TFs consist of two characteristic domains, an N-terminal HTH DNA-binding domain (PF00356) and a C-terminal effector-binding domain, which is homologous to periplasmic binding proteins of sugar ABC transporters (PF00532, PF13377, or PF13407). Gene orthology was defined by the bidirectional best-hit criterion implemented in the GenomeExplorer software (Mironov et al., 2000) and validated by phylogenetic trees from the MicrobesOnline database (Dehal et al., 2010). Genes were considered as orthologs if they (i) formed a mono- or paraphyletic branch of the phylogenetic tree and (ii) demonstrated conserved chromosomal gene context. For genes in the reconstructed regulons, correspondence to both of these criteria was sufficient for these genes to be considered as orthologs. For the studied LacI-family regulators, an additional criterion of orthology was used. Thus, orthologous groups of LacI-TFs should (i) form a mono- or paraphyletic group on the phylogenetic tree; (ii) have a conserved gene context; (iii) have highly similar TFBS motifs; and (iv) have the same effector specificity (known or predicted based on the regulon content).

For the regulon reconstruction, we used a previously established comparative genomics approach (Rodionov, 2007) implemented in the RegPredict Web server (Novichkov et al., 2010). This approach includes prediction of putatively regulated genes, inference of TFBSs, construction of positional weight matrices (PWMs) for TFBS motifs, and a further search for additional regulon members on the basis of predicted TFBSs in gene promoter regions. Overall, three main strategies were used for the reconstruction of regulons: (1) construction of PWMs on the basis of known TFBSs, for regulons being previously analyzed in model organisms; (2) prediction of novel TFBS motifs in promoter regions of regulated genes, for regulons with only regulated genes but not TFBSs known; and (3) prediction of putatively co-regulated genes followed by the inference of putative TFBS motifs in their promoter regions and attribution of a candidate TF to a putative regulon. Presumably, regulated genes were predicted by the analysis of conserved gene neighborhoods around an analyzed LacI-TF gene. Data about known regulated genes and LacI-TFBS motifs were extracted from the literature and from the RegTransBase (Cipriano et al., 2013), RegulonDB (Salgado et al., 2013), DBTBS (Sierro et al., 2008), and CoryneRegNet (Pauling et al., 2012) databases.

Candidate motifs in upstream regions of regulated operons were identified by the Discover Profile tool in RegPredict (Novichkov et al., 2010). A search for palindromic DNA motifs of 14- to 24-bp length was carried out within putative promoter regions from −400 to +100 bp relative to the translational gene start. Motifs were further manually validated by phylogenetic footprinting, that is, analysis of conserved islands in multiple alignments of DNA fragments (Shelton et al., 1997). The constructed PWMs were further used to search for additional regulon members using the Run Profile tool in RegPredict. The lowest score observed in the training set of known and/or predicted TFBSs was used as the threshold for a site search in genomes. To eliminate false-positive TFBS predictions, the consistency check approach (Ravcheev et al., 2007; Rodionov, 2007) and/or functional relatedness of candidate target operons were used. In this approach, an operon can be considered as regulated when its upstream region contains putative TFBSs with a score higher than the threshold, and such sites can be found in a number of related genomes. Operons were defined as groups of genes satisfying the following criteria: same direction of transcription, intergenic distance up to 100 bp, absence of internal TF-binding sites, and conservation of the locus structure in a number of related genomes. All predicted TFs, motifs, and sites are available in the RegPrecise database (http://regprecise*.*lbl*.*gov/) (Novichkov et al., 2013), where they are publicly available within the TF family collections of regulons.

Functional gene annotations were extracted from the literature and uploaded from the SEED (Disz et al., 2010), UniProt (Magrane and Consortium, 2011), MicrobesOnline (Dehal et al., 2010), and KEGG (Kanehisa et al., 2012) databases. Known functional annotations for a particular gene were expanded to all orthologous genes. For prediction of gene functions, both the comparative genomics and context-based methods were used (reviewed in Osterman and Overbeek, 2003; Overbeek et al., 2005; Rodionov, 2007). Multiple alignments of protein and DNA sequences were built by MUSCLE (Edgar, 2004). Phylogenetic trees were constructed using a maximum-likelihood algorithm implemented in PhyML 3.0 (Guindon et al., 2010) and visualized via Dendroscope (Huson et al., 2007) and iTOL (Letunic and Bork, 2011). Sequence logos for DNA motifs were drawn with WebLogo (Crooks et al., 2004).

# **RESULTS**

The comparative genomics workflow for regulon reconstruction implemented in the RegPredict Web server (Novichkov et al., 2010) and the RegPrecise database (Novichkov et al., 2013) includes three steps: (i) selection of a taxonomic group of related bacteria; (ii) selection of a subset of diverse genomes that represent a given group; and (iii) reconstruction of regulons in the selected genomes. For the analysis of LacI-TF regulons, we selected a set of 344 representative genomes from 39 taxonomic groups from 7 bacterial phyla (Table S1 in Supplementary Material). Among the analyzed lineages, there are 19 taxonomic groups of Proteobacteria (183 genomes), 9 groups of Firmicutes (72 genomes), and 7 groups of Actinobacteria (57 genomes). The Bacteroides, Chloroflexi, Deinococcus-Thermus, and Thermotogae phyla are each represented by a single taxonomic group and have 32 genomes in total.

# **REPERTOIRE OF LacI-TF GENES IN BACTERIAL GENOMES**

To estimate the abundance of LacI-family TFs (LacI-TFs) in the studied genomes, we collected primary LacI-TF sets using a similarity search and the existing prokaryotic TF compilations. In total, 2572 proteins were found unevenly distributed in most (309/344; 90%) of the studied genomes, whereas 10% of the genomes do not encode putative LacI-TFs (Table S1 in Supplementary Material). The largest average numbers of LacI-TFs per genome were found in several lineages of the Actinobacteria phylum including Streptomycetaceae and Bifidobacteriaceae (from 32 and to 17 regulators), in two lineages of Proteobacteria–Rhizobiales and Enterobacteriales (15 regulators in each group), and in two lineages of Firmicutes–Bacillales and Enterococcaceae (12 regulators in each group). The remaining taxonomic groups possess less than 10 LacI-TFs per genome on average. Noteworthily, the Methylophilales, Neisseriales, Nitrosomonadales, Oceanospirillales, Magnetospirillum/Rhodosprillum, and Desulfovibrionales groups completely lack LacI-TFs in their genomes. The absence of LacI-TFs in these taxonomic groups of Proteobacteria can be related to (i) relatively small proportion of sugar catabolic genes in their genomes (as LacI-TFs mostly control sugar catabolism, see below); (ii) increasing usage of TFs from other families to compensate the contraction of the LacI-TF pool.

#### **STATISTICS OF RECONSTRUCTED REGULONS AND REGULOGS**

The entire set of identified LacI-TFs was broken into taxonomic group-specific orthologous groups that were subjected to further comparative genomics analysis using the RegPredict Web server (Novichkov et al., 2010). Normally, an orthologous group contained no more than one TF per genome. However, in some cases TFs formed by recent, mainly genome-specific, duplication were assigned to the same orthologous group (Table S1 in Supplementary Material). By analyzing orthologous groups of regulators in each taxonomic group, candidate motifs and binding sites were predicted for 1303 LacI-TFs (50% of all putative LacI-TFs) in 272 bacterial genomes (80% of studied genomes). The main outcome of this analysis is an annotated regulog, which is defined as a set of genome-specific regulons controlled by orthologous TFs. Overall, we inferred 1281 LacI-TF regulons that constitute 322 populated regulogs unevenly distributed across 39 studied taxonomic groups of genomes (Tables S1, S2 in Supplementary Material). The reconstructed regulons included 7465 candidate sites, 6076 operons, and 13,558 genes.

The taxonomical distribution of the reconstructed LacI-TF regulons is highly uneven but generally follows the distribution of all LacI-family TFs (**Figure 1**): 57% of regulons were from Proteobacteria, about 30% from Firmicutes, 7% from Actinobacteria, about 1–2% from each of Thermotogales, Bacteroides, Chloroflexi, and Deinococcus/Thermus. Yet, compared to the genomic distribution of all putative LacI-TFs (Table

S1 in Supplementary Material), the Actinobacteria phylum is underrepresented. Based on the phylogenetic analysis of LacI-TF proteins, regulators from the reconstructed regulogs were merged into larger orthologous groups that were consistent with the taxonomy, had regulated orthologous genes, and had similar binding motifs. TFs were assigned to an orthologous group if they formed a mono- or paraphyletic branch (see below) in the phylogenetic tree (Figure S1 in Supplementary Material). As in most gene families containing multiple paralogs resulting from frequent duplications, losses, and horizontal transfers, the resolution of orthology in some cases was difficult and required arbitrary decisions that are supported by the genomic context and/or functional attributes of the reconstructed regulons.

As a result, the studied LacI-TFs were classified into 190 orthologous groups characterized by conserved DNA motifs and regulated pathways (Table S2 in Supplementary Material). Twothirds of the obtained orthologous groups (125/190) contain TFs from a single regulog, which is defined as a set of orthologous regulons in a group of closely related genomes. Thirty-seven orthologous groups include two regulogs, whereas the remaining 28 groups were assigned to three or more regulogs (**Figure 2A**). The total number of regulons (and corresponding TFs) per orthologous group of LacI-TFs varies between 1 and 59, with the average being 6.7 (**Figure 2B**). The maximal number of groups was observed for the groups including two TFs (37 groups). Orthologous groups containing up to 5 and between 6 and 10 regulons constitute 60 and 20% of all groups, respectively. The most populated groups of LacI-TF orthologs were found for the global catabolite control regulator CcpA in Firmicutes (59 regulons, 6 regulogs), the ribose repressor RbsR in Proteobacteria (52 regulons, 10 regulogs), the maltose repressor MalR in Firmicutes (42 regulons, 6 regulogs), as well as sugar catabolism regulators FruR (40 regulons, 6 regulogs), GntR (37 regulons, 7 regulogs), and GalR (27 regulons, 5 regulogs) in γ-Proteobacteria.

#### **GLOBAL AND LOCAL REGULONS**

The reconstructed LacI-TF regulons demonstrate drastic differences in the numbers of predicted target genes and operons. The majority of regulons (1198/1288, 93%) include 20 or fewer genes (**Figure 3A**), and further, three-fourths of these regulons contain between 2 and 7 genes, whereas 26 regulons have only 1 target gene. With respect to the number of regulated operons (**Figure 3B**), the largest portion of the studied LacI-TFs regulates one (31%) or two (38%) operons. We divided all reconstructed LacI-TF regulons into two main categories depending on their size (number of regulated genes and operons) and functional diversity (number of regulated pathways). A total of 125 regulons (12 regulogs) were classified as global, since each of them (i) contained more than 15 target genes that were arranged in at least 7 operons and (ii) controlled multiple metabolic pathways. The remaining LacI-TF regulons (1163/1288, 90%) were classified as local, each having a smaller number of targets and controlling a single metabolic pathway, usually a particular carbohydrate catabolic pathway.

Almost one-half of the identified global regulons (59 regulons, 6 regulogs) are operated by orthologs of the *B. subtilis* catabolite control regulator CcpA in all analyzed taxonomic groups of Firmicutes. CcpA regulons were previously described in detail for bacteria from the Bacillaceae (Sonenshein, 2007;

Fujita, 2009; Leyn et al., 2013), Staphylococcaceae (Seidl et al., 2009; Ravcheev et al., 2011), Lactobacillales (Mahr et al., 2000; Zheng et al., 2011, 2012; Zotta et al., 2012; Ravcheev et al., 2013b), and Clostridiaceae (Antunes et al., 2012) lineages. Another large group of global regulons (31 regulons, 3 regulogs) are operated by orthologs of the *E. coli* purine repressor PurR in three related taxonomic groups of γ-Proteobacteria— Enterobacteriales, Pasteurellales, and Vibrionales. PurR is a global transcriptional regulator of *E. coli*, controlling biosynthesis of purines, some steps of biosynthesis of pyrimidines, polyamine metabolism, and nitrogen assimilation (Ravcheev et al., 2002; Cho et al., 2011). The remaining three global regulogs— FruR in Enterobacterales, PckR in Rhizobiales, and GapR in Rhodobacterales (totaling 35 regulons)—that control the central and periphery carbohydrate catabolic pathways in diverse groups of Proteobacteria are described in more detail below.

Autoregulation was observed for 72% (943 TFs) of the studied LacI-TFs. The autoregulation of TFs was more typical for local regulons, as expected. Thus, less than one-half of the studied global regulators (55 TFs) demonstrated autoregulation, whereas more than three-fourths of the analyzed local regulators (888 TFs) were autoregulated.

#### **BINDING SITES AND MOTIFS**

Binding motifs of the considered LacI-TFs are palindromes formed by highly conserved inverted repeats, which is consistent with previous studies (Francke et al., 2008; Camas et al., 2010; Milk et al., 2010). The distance between the repeats is usually constant for a given orthologous group of TFs, although in rare cases there is some flexibility. The overwhelming majority of LacI-TF binding motifs (1251/1303 TFs) are even palindromes (16, 18, 20, or 22 nt long), whereas non-canonical palindromes (17, 19, or 21 nt long) were found for only 4% of regulators (e.g., LacR, GalR, and EbgR). The characteristic feature of even palindromes is the presence of a consensus CG pair in the center of the palindrome (1173 TFs). Nonetheless, in some even palindromes the central pair can be either different (49 TFs) or degenerate (29 TFs).

More than 75% of identified LacI-TFBSs are located within the area between 140 and 30 bp upstream of the start codon (**Figure 4A**). Less than 1% of sites are localized within coding regions, including experimentally demonstrated *E. coli* sites of LacI deep within the *lacZ* gene (Lewis, 2005) and PurR in the *purB* gene (He et al., 1992). Approximately 7% of sites are localized far upstream (*>*200 bp), and while some of them might in fact regulate divergently transcribed genes, experimental examples of such localization are known, e.g., the PurR site upstream of *prsA*, again in *E. coli* (He et al., 1993).

About 20% of regulated operons are preceded by more than one binding site. Thus, tandem sites were observed upstream of 1118 operons. The multiplicity of sites within regulatory regions of divergently transcribed genes (i.e., divergons) increases with the divergon length. Groups of three and four adjacent sites were found for 116 and 5 regulated operons, respectively.

The histogram of intersite distances for pairs of sites localized in the same intergenic region has pronounced peaks at 13, 22, and 32 bp (**Figure 4B**). The latter two lengths are multiples of the DNA helix step and hence are clearly indicative of cooperative

binding of two TF dimers. Several LacI-TF regulogs, e.g., RafR from Enterobacteriales and ScrR from Burkholderiales, have only double sites at the distance of 21–22 bp, suggesting that the cooperative binding in this case is obligatory. Overlapping sites, situated at a distance of about 13 bp, were observed upstream of 110 operons. Such an arrangement may be functional, as has been demonstrated for GntR-binding sites upstream of *gntKU* of *E. coli* (Tsunedomi et al., 2003a) and it is conserved for more than one-half of the operons regulated by GntR and its orthologs in γ-Proteobacteria.

A trivial explanation—that these observations are an artifact does not seem plausible for two reasons. Firstly, such multiple sites are conserved and even preferred in some orthologous groups. Secondly, the artifact hypothesis does not explain the preferred distance of 13 bp. In the available tertiary structures of LacI-TFs complexed with DNA (Schumacher et al., 1994; Barbier et al., 1997; Kalodimos et al., 2002), seven base-pairs nearest to the site center form contacts between bases and side residues. Hence, the site-overlap region strongly overlaps with the zone of specific TF-DNA contacts. While tetramer binding has been suggested, it is difficult to reconcile this with the structural data, as dimers bound at these distances cannot interact. It is possible that there exists a specific binding mode involving partial unwinding of the DNA strands.

#### **REGULATED METABOLIC PATHWAYS AND EFFECTORS**

By assessing the functional content of the reconstructed regulons, we tentatively predicted possible biological functions and effectors for 190 orthologous groups of LacI-TFs. As result, metabolic pathways were predicted for 182 groups of LacI-TFs. These include 54 groups that were only assigned to the general category of sugar metabolism, and for them, a specific sugar catabolic pathway remained unknown (Table S2 in Supplementary Material). We compared the predicted regulon functions with previous results of experimental studies available for 24 selected LacI-TFs. The previously established functions of these regulators are in good agreement with the target-regulated pathways that were predicted in this work. Based on the metabolic pathway reconstruction and the knowledge of pathway metabolites, a range of possible molecular effectors was suggested for 108 groups of LacI-TFs (Table S2 in Supplementary Material). Of these, effectors were previously known for regulators from 21 LacI-TF groups.

As expected, the overwhelming majority of the studied LacI-TF orthologous groups control carbohydrate metabolism (176/182; 96% of groups with assigned pathways). At that, most of the orthologous groups containing local LacI-TFs are assigned to specific carbohydrate utilization pathways. In contrast, five orthologous groups of local regulators including AdeR, HpxR, and UriR control the nucleoside utilization pathways, whereas a local regulator NtdR controls the neotrehalosadiamine biosynthesis. Most global regulators from the LacI family also are involved in the control of carbohydrate metabolism (FruR, CcpA, PckR, GapR), whereas PurR regulates several key metabolic pathways including purine and pyrimidine biosynthesis.

In agreement with the observed tendency of LacI-TF to control sugar metabolism, we report that carbohydrates constitute the largest class of effectors for these regulators (103/107; 96% of orthologous groups with assigned effectors). The majority of carbohydrate effectors assigned to the LacI-TF groups are monosaccharides and their derivatives (26/45), including hexoses (e.g., glucose, galactose, mannose), pentoses (e.g., ribose, xylose), sugar phosphates (e.g., fructose-1-phosphate, allose-6-phosphate), sugar acids (e.g., gluconate, galacturonate), sugar alcohols (e.g., ribitol), and amino sugars (e.g., N-acetylglucosamine). The second-largest category of carbohydrate effectors contains various oligosaccharides (14/26), including common disaccharides cellobiose, maltose, sucrose, and trehalose or their phosphorylated derivatives (e.g., cellobiose-6-phosphate, sucrose-6-phosphate). Finally, non-carbohydrate effectors of LacI-TFs include nucleobases (e.g., guanine and hypoxanthine are co-repressors of PurR in *E. coli* Schumacher et al., 1994), nucleosides (e.g., cytidine and adenosine are inducers of CytR in *E. coli* Barbier et al., 1997), and proteins (e.g., phosphoprotein HPr-Ser46-P is a co-repressor of CcpA in *B. subtilis* Schumacher et al., 2011).

Functional analysis of reconstructed LacI-TF regulons revealed that many sugar utilization pathways are regulated by two or more non-orthologous regulators. These include catabolic pathways for at least 10 distinct types of carbohydrates that are controlled by more than 90 orthologous groups of LacI-TFs (**Table 1**; Table S1 and Figure S1 in Supplementary Material). The observed large numbers of non-orthologous regulators for the glucoside and galactoside catabolic pathways correlate with structural diversity of glucose- and galactose-containing oligosaccharides that can be **Table 1 | Sugar utilization pathways controlled by non-orthologous LacI-TFs.**

utilized by bacteria in diverse natural habitats. On the other hand, the diversity of regulators for several other sugars including ribose and sucrose suggest a high frequency of convergent evolutionary events when the same ligand specificity has evolved independently in different branches of the LacI family.

For example, the sucrose utilization pathway is regulated by LacI-TFs from at least 11 orthologous groups from the phyla of Proteobacteria (11 lineages) and Firmicutes (5 lineages), as well as a single lineage of Actinobacteria (**Figure 5**). Analysis of the respective sucrose regulons revealed multiple distinct combinations of sucrose uptake transporters including permeases (*scrT*, *cscB*, *sut1*), phosphotransferase systems (PTSs) (*scrA*) and porins (*scrY*, *scrO*, *omp*), and sucrose catabolic enzymes including phosphorylases (*scrP*), hydrolases (*scrB*), and fructokinases (*scrK*). The effectors sucrose and sucrose-6-phosphate were assigned to the respective groups of ScrR regulators based on the type of regulated sucrose-specific transporters, i.e., a permease or a PTS, respectively. Interestingly, some families such as Bacillales and Enterobacteriales contain non-orthologous ScrR regulators with different effectors. As expected, DNA motifs of TFs controlling the sucrose catabolic pathway are well conserved within orthologous groups but are clearly different between non-orthologous regulators (**Figure 5**).

#### **GLOBAL REGULONS FOR CENTRAL SUGAR METABOLISM**

The LacI family contains a number of global regulators for central carbohydrate metabolic pathways, in particular the previously known regulators CcpA in Firmicutes and FruR in Enterobacteria. Here, we report a comparative genomics reconstruction of orthologous FruR regulons in γ-Proteobacteria, while the reconstructions of CcpA regulons in different lineages of Firmicutes has been reported previously (Ravcheev et al., 2011, 2013b; Antunes et al., 2012; Leyn et al., 2013) and is available in

shown by square brackets.

the RegPrecise database. Further, we report identification of three novel non-orthologous regulons (named PckR, GapR, and GluR) that control central carbohydrate metabolism in three lineages of α-Proteobacteria, namely, Rhizobiales, Rhodobacterales, and Caulobacterales. Below we provide functional analysis of regulons for each of these four LacI-family regulators in Proteobacteria.

The fructose repressor FruR (also known as the catabolism repressor and activator Cra) is a global regulator of central metabolism in *E. coli* (Saier and Ramseier, 1996). FruR/Cra coordinates the carbon flow by repressing glycolytic genes involved in the Embden–Meyerhof, Entner–Doudoroff, and pentose– phosphate pathways and by activating gluconeogenesis genes. The comparative genomics reconstruction of orthologous FruR regulons in γ-Proteobacteria revealed that the regulon size correlates with the taxonomy of studied groups and with the phylogeny of the FruR proteins (**Figure 6**). In Vibrionales and Pseudomonadales, the fructose utilization operon *fruBKA* and the *fruR* gene are the only members of the reconstructed FruR regulons; therefore, FruR operates as a local regulator in these lineages. In the Enterobacteriales, the FruR regulon is expanded to cover genes of the central glycolytic pathways, a part of the TCA cycle, and several fermentation and respiration pathways. Further regulon expansion to the genes of glyoxylate bypass is observed in the *Escherichia*, *Salmonella*, *Citrobacter*, *Enterobacter*, and *Klebsiella* species. Finally, in the closely related *Escherichia* and *Salmonella* species, the FruR/Cra regulon is expanded to include the Entner–Doudoroff pathway genes. Similar trends in the evolution of global regulons in closely-related bacterial species

circles show the numbers of genomes with correspondent regulation.

were previously demonstrated for the PhoP regulon in enterobacteria (Perez and Groisman, 2009). On the other hand, in Pasteurellales, the regulon is degrading: FruR is absent in the *Haemophilus* spp. and *Actinobacillus pleuropneumoniae*, whereas in *Pasteurella multocida*, *Mannheimia succiniciproducens, A. succinogenes*, and *A. aphrophilus* the *fruR* gene is seemingly intact, but no candidate FruR-binding sites could be detected.

The hypothetical TF PckR (SMc02975 in *Sinorhizobium meliloti*) was previously annotated as a putative regulator of the phosphoenolpyruvate carboxykinase *pckA* (EMBL accession number AF004316.1); however, it had not yet been studied experimentally. Orthologs of PckR were found in 10 out of 15 analyzed genomes from the Rhizobiales order including species from the Rhizobiaceae, Brucellaceae, Phyllobacteriaceae, and Xanthobacteraceae families. By using the comparative genomics approach, we identified the putative PckR binding motif and reconstructed the PckR regulons in each of these genomes (**Figure 7**). Furthermore, by analyzing the binding-site position within promoter regions, we predicted a negative or positive mode of PckR regulation. As result, PckR was predicted to function as a dual transcriptional regulator that represses glycolytic genes from the Embden–Meyerhof and Entner–Doudoroff pathways (*glk*, *fba*, *pykA*, *zwf-pgl-edd*, *eda*) and activates genes from the gluconeogenesis and TCA cycle (*pckA*, *mdh-sucCDAB*, *sdhABCD*).

In the Xanthobacteraceae family, the reconstructed PckR regulon includes a minimal number of genes (*pckA*, *pckR*, *edd*, and *hpr*), whereas it is significantly expanded in the other three families of Rhizobiales and contains from 14 to 28 genes per genome. The most conserved members of the reconstructed PckR regulons include the *pckA* and *edd* genes (in 10 and 9 genomes, respectively), the *zwf-pgl* and *mdh-sucCDAB* operons (in 8 genomes), and the *fba*, *glk*, and *eda* genes (in 7 genomes). In most of the analyzed genomes, the *pckR* genes are not clustered with their target genes, which is a common feature of many global regulators in bacteria. However, in the genome of *Xanthobacter autotrophicus*, *pckR* is divergently transcribed with the 6-phosphogluconate dehydrogenase gene *gnd*, which is preceded by a putative PckRbinding site.

Orthologs of PckR and their cognate DNA motifs were found in Rhizobiales but not in other lineages, suggesting that the PckR regulon was introduced relatively recently in the evolution of α-Proteobacteria. PckR from Rhizobiales can be considered as a partial functional replacement of the Enterobacterial Cra/FruR (see above) and the *Shewanella* HexR regulators (Leyn et al., 2011) that both play a pleiotropic role modulating the direction of carbon flow through different carbohydrate metabolic pathways. The molecular effector of PckR is unknown. A plausible hypothesis is that PckR dissociates from its DNA sites in response to an intermediate of the central glycolytic pathway in rhizobia.

A novel global transcriptional regulon for carbohydrate metabolism genes named GapR (RSP\_1663 in *Rhodobacter capsulatus*) was identified in all 13 studied genomes from the Rhodobacteraceae family (**Figure 7**). GapR was predicted to recognize a 20-bp palindromic DNA consensus, which is distinct from the PckR-binding consensus. The reconstructed GapR regulons contain from 7 to 18 genes per genome organized in 5–13 operons. GapR regulates glycolytic genes involved in the Embden–Meyerhof and Entner–Doudoroff pathways (*zwf-pglpgi*, *edd-eda*, *fba*, *gapA*, *gapB*, *eno*, *pykA)*, gluconeogenesis (*pckA*, *pycA*), fructose utilization (*scrK*), pentose phosphate pathway (*tal*), and the TCA cycle (*mdh*). A similar DNA motif was identified upstream of the *gapR* genes in five genomes, suggesting their autoregulation. Similarly to PckR in Rhizobiales, the *gapR* genes do not cluster on the chromosome with their target genes in most genomes. However, *Roseobacter* possess two copies of *gapR*. One of these paralogs is divergently transcribed with the *gapB* gene, which is preceded by a putative GapR-binding site. The molecular mechanism and effector for GapR regulators remain to be elucidated.

Although PckR and GapR regulate overlapping sets of genes from the central carbohydrate metabolism in two distinct lineages of α-Proteobacteria, these regulators are not orthologous to each other (Figure S1 in Supplementary Material) and recognize different DNA motifs (**Figure 7**). Another non-orthologous regulator from the LacI family named GluR (CC2053 in *Caulobacter crescentus*) was predicted to control the central carbohydrate metabolism in the Caulobacteraceae family (**Figure 7**). GluR recognizes a conserved 20-bp palindromic consensus, different from the PckR- and GapR-binding motifs above. The reconstructed GluR regulon in the *Caulobacter* spp. is composed of the glycolytic genes *zwf-pgl-edd-glk*, *pykA,* and *gnl-gfo*, as well as the gluconeogenic gene *ppdK* encoding pyruvate-phosphate dikinase and the *sucABCD* operon involved in the TCA cycle. The predicted regulatory gene *gluR* is located immediately downstream of the target operon *zwf-pgl-edd-glk* and is preceded by a conserved GluR-binding site in all three analyzed *Caulobacter* genomes.

Orthologs of GluR were not found in other α-Proteobacteria, although bacteria from the Caulobacterales lineage have three paralogs (BglR1-3) that are predicted to control the cognate operons involved in the β-glucoside utilization (Figure S1 in Supplementary Material). The molecular effector of GluR is not known. Based on its close similarity to the predicted BglR repressors that possibly respond to β-glucosides and/or glucose, we propose that GluR dissociates from its DNA sites in response to glucose. In confirmation of our hypothesis, it has been demonstrated that glucose induces expression of the Entner–Doudoroff pathway genes and that the *edd* and *glk* genes are essential for the glucose utilization in *C. crescentus* (Hottes et al., 2004).

# **DISCUSSION**

We used the comparative genomics reconstruction of regulons for analysis of the LacI family of bacterial TFs. This choice was based on the following features. First, the LacI family is large, varied, and broadly distributed in bacteria. Second, proteins from this family have a rigid domain structure and a highly conserved structure of TFBS motifs. This study resulted in a detailed reconstruction of 1281 LacI-TF regulons in 272 bacterial genomes. Most (∼90%) of the analyzed TF-LacI regulons are local, i.e., they control a small number of genes and operons that are involved in only one metabolic pathway. However, some LacI regulons are global, controlling tens to hundreds of genes involved in multiple metabolic pathways. In addition to the reconstruction of previously known global regulons, such as FruR/Cra, PurR,

and CcpA, we identified three novel regulators for the central carbohydrate metabolism in α-Proteobacteria, PckR, GluR, and GapR, and reconstructed their corresponding regulons. For two of these global regulons, FruR and GluR, we reconstructed their possible evolutionary histories. Both these TFs likely originated from local regulators during a process of gradual regulon expansion.

A large-scale phylogenetic analysis of LacI-TFs reveals numerous examples of various evolutionary processes for regulators and their regulons including divergent evolution (diversification of TF functions and binding specificities after duplication), convergent evolution (appearance of the same function in distantly related branches of a phylogenetic tree), and formation of paraphyletic groups (origin of novel functions and specificities, noncharacteristic for a given branch of TFs). Below we discuss these evolutionary processes in more detail and provide examples of functional diversification within the LacI-TF family.

The LacI-TF phylogenetic tree demonstrates that some orthologous groups of TFs form branches consistent with the taxonomy, with TFs regulating orthologous genes, and with recognizing similar motifs, but those branches included an internal clade with demonstrated differences in the motif structure, effector specificity, or regulon content, i.e., these TFs formed so-called paraphyletic groups. The most interesting example of paraphyly is the branch of ribose repressors RbsR in β- and γ-Proteobacteria that have the purine repressor PurR as an excluded clade. Phylogenetic analysis revealed the presence of PurR orthologs in only three bacterial orders, Enterobacteriales, Pasteurellales, and Vibrionales. The closest PurR paralogs in these groups are ribose repressors (RbsRs) (Figure S1 in Supplementary Material). Most probably, PurR was originated by duplication of RbsR in the common ancestor of Enterobacteriales, Pasteurellales, and Vibrionales. The Pseudomonadales and β-Proteobacteria with a single RbsR repressor seemingly feature the ancestral state. The RbsR orthologs from the above three orders retain the ligand specificity but have a slightly modified DNA binding motif, compared to RbsR of Pseudomonadales and β-Proteobacteria (**Figure 8**). On the contrary, PurR retained the motif but changed the ligand and the default state, as it binds DNA in the presence of its ligand, whereas RbsR, like the majority of the LacI-TFs, binds DNA in the absence of the ligand.

This example illustrates a possible way of formation of paraphyletic branches by duplication of an ancestral gene with subsequent conversion of one copy of the copies, resulting in the origin of a novel function. In the case of RbsR and PurR, duplication and conversion are rather ancient, and we observe only the result of such evolution but cannot observe the intermediate states of the evolutionary process. On the contrary, for the α-glucoside utilization regulon AglR and the trehalose regulon ThuR, such an intermediate state is clearly observable. AglR regulators from Rhizobiales and Rhodobacteriales form a paraphyletic branch with ThuR regulators from Rhizobiales as an excluded clade (Figure S1 in Supplementary Material). Both AglR and ThuR are local regulators, each controlling expression of the regulator genes and one other operon, divergently transcribed with the latter. The regulated operons contain homologous genes for kinases and ABC transporters, but non-homologous genes for hydrolases. The TFBS motif for ThuR (natcnAAAnCGnTTTngatt) is different from the one for AglR (nnntcAAAGCGCTTTgannn). Thus, during diversification, ThuR changed both the ligand and motif specificity. For the paraphyletic group CelR in γ-Proteobacteria, the AscG regulator in Enterobacteriales is an excluded clade. In

the case of the CelR and AscG regulators, their TFBS motifs and sets of regulated genes are drastically changed after duplication, but the effector specificity (cellobiose-6P) and the target metabolic pathway (cellobiose utilization) are retained.

The phylogenetic tree for the analyzed LacI-TFs (Figure S1 in Supplementary Material) is patchy, complicating the reconstruction of the evolutionary history. However, in some cases we can observe two or more adjacent branches, each corresponding to one orthologous group. A natural explanation is that these groups appeared as a result of initial duplication of a TF gene followed by diversification of copies. The LacI-TF tree contains multiple examples of monophyletic taxon-specific branches that consist of proteins from a single bacterial lineage. One such branch includes the *Bacteroides* UxaR, UxuR, and KdgR regulators that control the catabolic pathways for galacturonate, glucuronate, and 2-keto-3 deoxygluconate, respectively. Two other monophyletic branches include (i) kojibiose regulator KojR and the unknown sugar regulator Caur\_3448 from *Chloroflexus* and (ii) ribose regulator RbsR and uridine regulator UriR from *Corynebacteria* spp. Previously, a similar situation was observed for the ROK family of sugarspecific regulators in the deeply branched Thermotogales lineage (Kazanov et al., 2013).

Other examples of divergent evolution of regulator specificity are often demonstrated by adjacent branches in the phylogenetic tree. Thus, idonate repressor IdnR in Enterobacteriales and gluconate repressor GntR in multiple lineages of γ-Proteobacteria are the closest paralogs (Figure S1 in Supplementary Material). In the IdnR regulons, the *idnK* and *idnT* genes are the closest homologs of the GntR-regulated genes *gntK* and *gntU*, respectively. Thus, duplication affected not only a TF gene, but also some of the regulated genes. The ability of GntR to recognize IdnR-binding sites upstream of the *idnK* and *idnDOTR* operons in *E. coli* (Tsunedomi et al., 2003b) also confirms a recent duplication of these regulators. Structural similarity of sugar effectors for GntR and IdnR also points to a recent duplication and further specialization of IdnR in Enterobacteriales.

Another scenario of TF diversification is duplication of an ancestral TF gene followed by acquisition of novel regulated genes and, accordingly, new effector specificity. An example is provided by a branch containing fructose (FruR) and sucrose (ScrR) orthologous groups in γ-Proteobacteria. Sucrose is a fructosecontaining disaccharide. Because the FruR- and ScrR-regulated genes are functionally and structurally different, the regulator was most probably duplicated alone, and then one copy acquired new regulated genes. The acquisition of a novel regulatory function was coupled with the changes in the cognate TFBS motif. Here, divergent TFs have structurally similar effectors, fructose-1,6-biphosphate for FruR and sucrose-6-phosphate for ScrR.

Based on the analysis of paraphyletic branches and adjacent monophyletic branches, three main types of the origin of TF with novel functions can be described. The first type is duplication of both the TF gene and a regulated gene or operon followed by diversification, as in the case of ThuR in Rhizobiales. During diversification, the TF and regulated genes change their specificities, some regulated genes may be lost, and novel genes may be included in the regulon. The second type is duplication of a regulator followed by acquisition of novel regulated genes, as for PurR in γ-Proteobacteria or for AscG in Enterobacteriales. The third type is rare, with only one example—the SCO5692 regulon in Streptomycetaceae. In this case, novel specificities for TFs originated without duplication, probably resulting from the loss of regulated genes.

In the process of acquisition of a new function, three characteristics of a TF regulon can be changed: (i) a set of regulated genes, (ii) effector specificity, and (iii) a TFBS motif structure. In most cases of TFs with a novel function, we observed the change of at least two of these characteristics. Change of all three characteristics is rarer and is usually observed in TFs from deeply branched lineages of bacteria such as Thermotogales. Most probably, in these cases change of all three characteristics is a result of a long evolutionary history. Change of only one of these characteristics is observed for only recently duplicated TFs, for which the diversification process is not yet complete.

In summary, the obtained extensive dataset for the LacI-TF family provides numerous examples of various evolutionary processes for regulators and their regulons. These data are publicly available in the RegPrecise database within the LacI family collection, which will enable further detailed analysis of signature residues in both DNA- and ligand-binding domains of regulators and establishment of the correlations between these residues and specificities toward the DNA motifs and molecular effectors they recognize.

# **AUTHOR CONTRIBUTIONS**

Dmitry A. Rodionov, Pavel S. Novichkov, and Mikhail S. Gelfand conceived and designed the research project. Dmitry A. Ravcheev, Dmitry A. Rodionov, and Mikhail S. Gelfand wrote the manuscript. Dmitry A. Ravcheev, Dmitry A. Rodionov, Matvei S. Khoroshkin, Olga N. Laikova, Olga V. Tsoy, Natalia V. Sernova, and Svetlana A. Petrova performed comparative genomics analysis to reconstruct TF regulons. Dmitry A. Rodionov provided the quality control of annotated regulons in the RegPrecise database. Matvei S. Khoroshkin analyzed statistical properties of TF regulons. All authors read and approved the final manuscript.

# **ACKNOWLEDGMENTS**

This work was supported by the Russian Foundation for Basic Research (14-04-00870, 14-04-91154) and by the Russian Academy of Sciences via the programs "Molecular and Cellular Biology" and "Living Nature." This research was also partially supported by the Genomic Science Program (GSP), Office of Biological and Environmental Research (OBER), and U.S. Department of Energy (DOE) and is a contribution of the Pacific Northwest National Laboratory (PNNL) Foundational Scientific Focus Area. Preliminary reconstruction of regulons was done by students of the Spring 2013 Bioinformatics course at the Faculty of Bioengineering and Bioinformatics in the Lomonosov Moscow State University as a part of regular coursework.

# **SUPPLEMENTARY MATERIAL**

The Supplementary Material for this article can be found online at: http://www.frontiersin.org/journal/10.3389/fmicb. 2014.00294/abstract

# **REFERENCES**

in proteobacterial genomes. *Nucleic Acids Res.* 29, 774–782. doi: 10.1093/nar/ 29.3.774

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 06 April 2014; paper pending published: 12 May 2014; accepted: 28 May 2014; published online: 11 June 2014.*

*Citation: Ravcheev DA, Khoroshkin MS, Laikova ON, Tsoy OV, Sernova NV, Petrova SA, Rakhmaninova AB, Novichkov PS, Gelfand MS and Rodionov DA (2014) Comparative genomics and evolution of regulons of the LacI-family transcription factors. Front. Microbiol. 5:294. doi: 10.3389/fmicb.2014.00294*

*This article was submitted to Microbial Physiology and Metabolism, a section of the journal Frontiers in Microbiology.*

*Copyright © 2014 Ravcheev, Khoroshkin, Laikova, Tsoy, Sernova, Petrova, Rakhmaninova, Novichkov, Gelfand and Rodionov. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Connecting lignin-degradation pathway with pre-treatment inhibitor sensitivity of *Cupriavidus necator*

*Wei Wang1 \*†, Shihui Yang2 \*†, Glendon B. Hunsinger 2, Philip T. Pienkos <sup>2</sup> and David K. Johnson1*

*<sup>1</sup> National Renewable Energy Laboratory, Biosciences Center, Golden, CO, USA*

*<sup>2</sup> National Renewable Energy Laboratory, National Bioenergy Center, Golden, CO, USA*

#### *Edited by:*

*Katherine M. Pappas, University of Athens, Greece*

#### *Reviewed by:*

*Harold J. Schreier, University of Maryland Baltimore County, USA Claudio Avignone-Rossa, University of Surrey, UK*

#### *\*Correspondence:*

*Wei Wang and Shihui Yang, 15013 Denver West Parkway, Golden, 80401 CO, USA e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
*†These authors have contributed equally to this work.*

To produce lignocellulosic biofuels economically, the complete release of monomers from the plant cell wall components, cellulose, hemicellulose, and lignin, through pre-treatment and hydrolysis (both enzymatic and chemical), and the efficient utilization of these monomers as carbon sources, is crucial. In addition, the identification and development of robust microbial biofuel production strains that can tolerate the toxic compounds generated during pre-treatment and hydrolysis is also essential. In this work, *Cupriavidus necator* was selected due to its capabilities for utilizing lignin monomers and producing polyhydroxylbutyrate (PHB), a bioplastic as well as an advanced biofuel intermediate. We characterized the growth kinetics of *C. necator* in pre-treated corn stover slurry as well as individually in the pre-sence of 11 potentially toxic compounds in the saccharified slurry. We found that *C. necator* was sensitive to the saccharified slurry produced from dilute acid pre-treated corn stover. Five out of 11 compounds within the slurry were characterized as toxic to *C. necator,* namely ammonium acetate, furfural, hydroxymethylfurfural (HMF), benzoic acid, and p-coumaric acid. Aldehydes (e.g., furfural and HMF) were more toxic than the acetate and the lignin degradation products benzoic acid and p-coumaric acid; furfural was identified as the most toxic compound. Although toxic to *C. necator* at high concentration, ammonium acetate, benzoic acid, and p-coumaric acid could be utilized by *C. necator* with a stimulating effect on *C. necator* growth. Consequently, the lignin degradation pathway of *C. necator* was reconstructed based on genomic information and literature. The efficient conversion of intermediate catechol to downstream products of cis,cis-muconate or 2-hydroxymuconate-6-semialdehyde may help improve the robustness of *C. necator* to benzoic acid and p-coumaric acid as well as improve PHB productivity.

**Keywords:** *Cupriavidus necator***, pre-treatment inhibitor, saccharified slurry, deacetylation, lignin degradation, biofuel, polyhydroxylbutyrate (PHB), genomics**

## **INTRODUCTION**

Lignocellulosic biomass is considered a renewable and sustainable source for energy production. The great environmental, energy security, and economic benefits of biofuels have driven research on biomass throughout the world, resulting in several pilot and demonstration scale projects, based primarily on cellulosic ethanol. One of the leading routes for cellulosic ethanol production is based on the deconstruction of biomass to a monomeric sugar solution produced by chemical pre-treatment and enzymatic saccharification termed a saccharified slurry (hereafter referred to as slurry) followed by fermentation of the sugars to ethanol. However, the physico-chemical properties of ethanol limit its further penetration into the current petroleum-based transportation fuel infrastructure (Serrano-Ruiz and Dumesic, 2011). The lower energy density of ethanol due to its higher oxygen content compared to conventional petroleum hydrocarbons, is a significant disadvantage to using bioethanol as a fuel alternative.

Current transportation fuel infrastructure has been established for hydrocarbon-based petroleum. Many advanced biofuels under development exhibit advantages such as higher energy densities, longer carbon chains and lower oxygen numbers. For example, there are some fuels derived from fatty acids with carbon numbers in the range from C12 to C22 that can be upgraded to high-energy-density hydrocarbons via hydrotreating. Therefore, development of hydrocarbon fuels, which are compatible with the current fuel infrastructure as replacements or blendstocks for gasoline, jet, and diesel, has recently received great attention.

The primary routes leading to hydrocarbon synthesis from biomass are thermochemical, biochemical, and hybrid approaches which employ elements of both. Bacteria, yeasts, or fungi can naturally synthesize fatty acids, isoprenoids, and polyalkanoates for energy storage. Although these compounds are mostly exploited in pharmaceutical, nutritional, and packaging sectors, they also have great potential for production of hydrocarbon fuels (Zhang et al., 2011). However, several obstacles need to be solved for future hydrocarbon production at the commercial scale.

Moreover, the efficient utilization of lignin, a major cell component, would greatly improve the economics of advanced biofuel production. Due to its highly recalcitrant aromatic structure, lignin is seen largely as an impediment to biochemical processes due to its reinforcement of the cellulose/hemicellulose matrix complicating the saccharification process (Zeng et al., 2014). Owing to the high energy state of its aromatic ring structures, lignin could have great value as a substrate for energy production, but so far very limited opportunities for microbial transformation of lignin to high value-added fuels or fuel intermediates have been identified. Some bacteria are known to degrade lignin monomers via the β-ketoadipate pathway (Harwood and Parales, 1996). Currently lignin is almost untouched in lignocellulosic ethanol fermentations and remains behind in the fermentation residue after distillation. These residues are mostly burned to provide power in current cellulosic ethanol production processes, which places a low value on this material. To reach the goal of utilizing lignin as a carbon source for advanced biofuel production, the composition of lignin degradation compounds as well as the toxic effects of these compounds on hydrocarbon-producing microorganisms should be investigated.

In this study we are seeking a biological understanding of the impact of model inhibitors, which include lignin degradation products as well as furans and acetate, on hydrocarbon-producing microorganisms. *Cupriavidus necator* (Syn. *Alcaligenes eutrophus, Ralstonia eutropha*) has been extensively studied for production of polyhydroxybutyrate (PHB), which consists of a C4 repeating unit that can be thermally depolymerized and then decarboxylated to propene (Fischer et al., 2011; Pilath et al., 2013), an intermediate, which can be upgraded to hydrocarbon fuels via commercial oligomerization technologies. *C. necator* has been reported to be able to utilize lignin monomers as a carbon source (Pérez-Pantoja et al., 2008).

The genome sequence of *C. necator* H16 has been published with *in silico* genome modeling and a developed genetics system (Pohlmann et al., 2006; Park et al., 2011; Brigham et al., 2012). In addition, several transcriptomic studies have recently been reported (Peplinski et al., 2010; Brigham et al., 2012), and the genome sequences for a number of other *Cupriavidus* spp. are also now available (Amadou et al., 2008; Pérez-Pantoja et al., 2008; Janssen et al., 2010; Lykidis et al., 2010; Poehlein et al., 2011; Cserhati et al., 2012; Hong et al., 2012; Van Houdt et al., 2012; Li et al., 2013). This information will facilitate future comparative genomics and systems biology studies to develop *C. necator* H16 as a robust and metabolically diverse hydrocarbon-intermediate production strain. Genomics is applied in this study to explore the metabolic pathways related to lignin utilization and response to toxic compounds in slurries, which will provide perspectives for strain metabolic engineering toward future economic hydrocarbon production using lignin.

# **MATERIALS AND METHODS STRAINS AND MEDIA**

The strain used in this study is a glucose-utilizing mutant of *C. necator* H16 (wild-type H16 is not able to metabolize glucose) (Orita et al., 2012), *C. necator* 11599, which was purchased from NCIMB culture collection. It is routinely cultured in LB at 37◦C. A minimal medium recipe was selected for this study (Cavalheiro et al., 2009). Specifically, the defined minimal medium for *C. necator* (per liter, pH 6.8) was: 10 g glucose, 1.0 g (NH4)2SO4, 1.5 g KH2PO4, 9 g Na2HPO4 12H2O, 0.2 g MgSO4 7H2O, 1.0 mL trace element solution. The Trace Element Solution (per liter): 10 g FeSO4.7H2O, 2.25 g ZnSO4.7H2O, 0.5 g MnSO4.5H2O, 2 g CaCl2.2H2O, and 1 g CuSO4.5H2O, 0.23 g Na2B4O7.10H2O, 0.1 g (NH4)6MO7O24, 10 mL 35% HC1.

#### **PRODUCTION OF SACCHARIFIED SLURRY AND MOCK MEDIA**

A deacetylated saccharified slurry, which was produced from the modified sulfuric acid pre-treatment and enzymatic hydrolysis of corn stover including an added deacetylation step before pre-treatment, was used in this study (Chen et al., 2012). The composition of the mock sugar media simulating the saccharified slurry is summarized in **Table 1**. The composition is based on the composition of the saccharified slurry in fermentation media at the level of 20% total solids.

### **GROWTH OF** *C. NECATOR* **ON SACCHARIFIED SLURRY**

*Cupriavidus necator* was first grown in 5 mL of LB in 125 mL baffled flasks, cultured at 200 rpm, and 37◦C. After 1 day, a 10% inoculum was added to 50 mL of fermentation media in a 250 mL flask and incubated in a shaker at 37◦C and 180 rpm for 4 days. The fermentation media contained either mock sugar slurry as shown in **Table 1** or saccharified slurry supplemented with tryptone (10 g/L) and yeast extract (5 g/L) as nutrients. Mock slurry was added at a level to achieve the same sugar concentrations (e.g., the glucose concentration in the 2X-diluted mock medium was 50 g/L). All experiments were run in duplicate.

### **PHB ANALYSIS**

The PHB content of the bacterial cells was determined by a quantitative method that used HPLC analysis to measure the crotonic acid formed by acid-catalyzed depolymerization of PHB (Karr et al., 1983). Cell mass samples were freeze-dried before analysis. PHB-containing dried bacterial cells (15–50 mg) were then digested in 96% H2SO4 (1 mL) at 90◦C for 1 h. The reaction

vials were then cooled on ice, after which, ice-cold 0.01N H2SO4 (4 mL) was added followed by rapid mixing. The samples were further diluted 20- to 150-fold with 0.01N H2SO4 before analysis by HPLC.

The concentration of crotonic acid was measured at 210 nm using an HPLC equipped with a photodiode array detector (Agilent 1100, Agilent Technologies, Palo Alto, CA). A Rezex RFQ Fast Acids column (100 × 7*.*8 mm, 8μm particle size, Phenomenex, Torrance, CA) and Cation H+ guard column (BioRad Laboratories, CA) operated at 85◦C were used to separate the crotonic acid present in the reaction solutions. The eluent was 0.01N H2SO4 at a flow rate of 1.0 mL min−1. Samples and crotonic acid standards were filtered through 0.45μm pore size nylon membrane syringe filters (Pall Corp., NY) prior to injection onto the column. The HPLC was controlled and data were analyzed using Agilent ChemStation software (Rev.B.03.02).

### **CHARACTERIZATION OF POTENTIAL INHIBITORS**

The selection of chemicals for sensitivity assays was based on analyses of saccharified slurry by GC-MS and LC-MS analyses performed at NREL as well as ICP-MS analyses carried out by Huffman Labs (Golden, CO). The compounds identified with high concentrations in the slurry and their potential derivatives were selected for further investigation. The compounds included ammonium (added to neutralize pre-treated corn stover) with two common anions of acetate (released by hydrolysis of hemicelluloses) and sulfate (from the sulfuric acid pre-treatment); sugar degradation products furfural and HMF; lignin monomers vanillin, coumaric acid, ferulic acid, and 4-hydroxybenzaldehyde as well as benzoic acid, a common intermediate from lignin monomer (coumarate and cinnamate) degradation (**Figure 5**). We also included products from the oxidation of lignin monomers, vanillic acid and 4-hydroxybenzoic acid, in a similar concentration range to their aldehyde forms (**Table 2**). The 1X concentrations for some chemicals (e.g., benzoic acid) in the slurry used for testing are based on concentrations used in a previous study.

### **BIOSCREEN C HIGH THROUGHPUT TOXICITY ASSAY**

The high-throughput Bioscreen C assay was carried out as reported previously (Franden et al., 2009, 2013). Briefly, *C. necator* cells were revived from overnight LB culture with OD600nm adjusted to 3.0 using minimal medium. This cell suspension was used as seed culture to inoculate Bioscreen C plates containing 290μL minimal medium per well at an initial OD600nm of 0.1. Growth was then monitored using the Bioscreen C instrument (GrowthCurves USA, NJ) with three technical replicates. The experiments were repeated at least two times.

*Cupriavidus necator* grown in the absence of potential inhibitory compounds was used as the control, and inhibition studies utilized cultures of *C. necator* challenged with different concentrations of each compound ranging from 0.1- to 10-fold (0.1X to 10X, **Table 2**). Stock solutions of compounds at 10X concentrations were prepared by dissolving the compounds to be tested in the minimal medium. These stock solutions were then diluted in minimal medium for testing at lower concentrations. For certain compounds with low aqueous solubility, incubation at 55–60◦C in a hot dry bath for several hours was needed for complete dissolution. The pH of the stock solutions was adjusted to the desired point of 6.8 using ammonium hydroxide (NH4OH) or sulfuric acid (H2SO4) and then filter sterilized before using.

# **GENOMIC INVESTIGATION FOR LIGNIN DEGRADATION PATHWAY RECONSTRUCTION**

To identify the enzymes related to lignin degradation, the protein sequences of *C. necator* H16 were extracted and reannotated functionally. Briefly, 6626 protein sequences were downloaded from NCBI (Genbank#: AM260479) and imported into CLC Genomics Workbench (V5.5) as the reference protein sequences for Blast search. In addition, the protein sequences were also imported into Blast2GO for the functional annotation (Gotz et al., 2008). The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were then extracted as well as the information of euKaryotic Orthologous Groups (KOG), enzyme code, and the reaction substrate(s) and product(s). The potential homologous gene(s) in *C. necator* H16 were identified by re-iterated BlastP searches. The information on protein product and conserved domains were examined and the pathway was reconstructed with the enzyme and pathway information from a literature search.

### *Gas chromatographic analysis*

Analysis of samples was performed on an Agilent 7890 GC equipped with a 5975 MS (Agilent Technologies, Palo Alto, CA). Sample compounds were separated using a 30 m × 0.25 mm × 0.25 mm DB-FFAP column (Agilent). HP MSD Chemstation software (Agilent) equipped with NIST database Rev. D.03.00 was used to determine the identity of the unknown compounds found within the samples.

Each sample was placed on an auto-sampler (Agilent) and injected at a volume of 1μL into the GC-MS (Agilent). The GC-MS method consisted of a front inlet temperature of 250◦C, MS transfer line temperature of 280◦C, and a scan range from 35 to 550 m/z. A starting temperature of 50◦C was held for 1.5 min and then ramped at 11◦C/min to a temperature of 165◦C with no hold time, then continued at a ramped rate of 35◦C/min to 250◦C and held for 9.617 min. The method resulted in a run time of 28 min for each sample.

# **RESULTS AND DISCUSSIONS**

# **GROWTH KINETICS OF** *C. NECATOR* **IN THE SACCHARIFIED SLURRY**

*Cupriavidus necator* did not grow in the original deacetylated saccharified slurry (**Figure 1A**) even in the 2-fold diluted saccharified slurry (**Figure 1B**) though it was able to grow in media containing the mock slurry. The sugars provided by the mock slurry allowed *C. necator* for superior growth over the LB alone in both 1-fold concentration and 2-fold dilution conditions with fast growth and high final OD600nm value. This complete inhibition of growth in real slury suggests that *C. necator* is not tolerant to the combination of all toxic compounds present in the saccharified slurry.

### **KINETICS OF** *C. NECATOR* **IN THE SACCHARIFIED SLURRY AFTER ACTIVATED CARBON (AC) TREATMENT**

It is well-known that acid pre-treatment of corn stover can generate various kinds of potential toxic inhibitors including organic

*The concentration (mM) following the chemical abbreviation is the concentration in the slurry calculated from analytical results of GC-MS, LC-MS, or ICP-MS.*

acids and aldehydes, and can also add to the toxicity through the build up of inorganic salts from the sulfuric acid used for pre-treatment and ammonium hydroxide added for neutralization (e.g., organic and inorganic acids of levulinic acid, vanillic acid, hydroxybenzoic acid, sulfuric acid, and aldehydes of HMF

**Table 3 | The compositions of deacetylated saccharified slurry before and after activated carbon (AC) treatment.**

*ND, Not detected; the detection limitation for furfural is 2.6 mM, HMF is 2.8 mM.*

and furfural). Furfural is generally considered the most potent inhibitor to various microbial catalysts such as *E. coli*, *Z. mobilis*, and yeast (Zaldivar et al., 1999; Liu et al., 2004, 2005, 2008; Gorsich et al., 2006; Allen et al., 2010; Bowman et al., 2010; He et al., 2012; Huang et al., 2012; Franden et al., 2013). The maximum amount of furfural was 2 g/L in the slurries in this study. Shake flask growth experiments showed that furfural greatly inhibited cell viability when the concentration was above 2 g/L (data not shown).

A chemical treatment process using activated carbon (AC) was applied in this study in order to improve growth by removing potential inhibitors such as furfural from the slurry. AC was added to the slurry at 0.1 g/mL loading and incubated for 2 h at 130 rpm, 24◦C. The compositions of saccharified slurry before and after AC-treatment indicated that all the furfural and about 30% of the acetate was removed by the AC treatment (**Table 3**).

The impact of AC treatment on the toxicity of diluted saccharified slurry (4-fold dilution), which provided about 25 g/L of glucose to the growth media, was tested. As shown in **Figure 2A**, untreated saccharified slurry was very toxic to *C. necator,* while AC-treated saccharified slurry was much less inhibitory to growth, although it was still inhibitory compared to the

mock slurry control. It is noteworthy that the consumption of glucose from AC-treated slurry was far slower than that from mock slurry and at 24 h, only about 25% of the total glucose was removed from the culture medium compared to 85% in the mock slurry (**Figure 2B**). This indicates that inhibitors remained in AC-treated slurry although furfural was completely removed. The remaining inhibitors are either unidentified residual components in the slurry (other than furfural) or

components released from pre-treated corn stover during saccharification. It is noteworthy that besides the removal of all furfural by AC-treatment, our preliminary data indicated that a significant fraction of the other potentially toxic compounds such as acetic acid (**Table 3**) and lignin degradation products (e.g., vanillin, 4-hydroxybenzaldehyde, and p-coumaric acid) were also removed by AC-treatment. The impact of these individual toxic compounds in the saccharified slurry will be evaluated in the following section.

The cell yield from the AC-treated slurry (4-fold diluted, 4X) was about one third that of the culture grown on the mock slurry (4X), with similar amounts of PHB accumulated in the cells from both growth conditions 96 h post-inoculation (**Figure 2C**). However, when 2-fold diluted slurry was used, increased inhibition led to a lower cell mass and accordingly a lower PHB yield even with the carbon treated slurry (**Figure 2C**). All these results indicate that the toxicity of the biomass-derived slurry had a major impact on PHB production. Moreover, the low cell viability on 2x diluted slurry indicates that the inhibitory effect on cells is likely due to the combined contribution of several inhibitors besides furfural. Although the efficacy of the AC-treatment is obvious in term of cell growth, it only partially mitigated the inhibition. Identification of the other inhibitors in the AC-treated slurry is needed to help eliminate or mitigate the observed inhibition.

#### **INHIBITOR SENSITIVITY INVESTIGATION USING THE HIGH-THROUGHPUT BIOSCREEN C ASSAY**

To systematically explore the impact of inhibitors in the saccharified slurry on *C. necator*, a high-throughput Bioscreen C growth assay was performed using growth medium augmented with the potentially toxic compounds that had been identified at concentration ranges discussed above (**Table 2**). Briefly, growth curves were generated by subtracting OD readings from the test wells in the Bioscreen C from the background wells containing blank medium only. Typical growth curves with *C.*

**Table 4 | The growth rates and responses of** *C. necator* **to various furfural concentrations (mM) supplemented in the minimal medium.**

**Table 5 | The IC50 values and the ratio of the IC50 value to Conc. value (the concentration of the toxic compound identified in the saccharified slurry, mM) for** *C. necator.*

**Table 6 | Response values of** *C. necator* **to 11 potentially toxic compounds identified in the saccharified slurry at different concentrations (mM) supplemented in the minimal medium.**

*The values in bold font indicate stimulation of growth by the compound at the corresponding concentration.*

**concentrations, 5 days post inoculation.** R1, R2, and R3 are three technical replicate wells in the Bioscreen C honeycomb plate.

*necator* grown with varied concentrations of furfural are shown in **Figure 3**. The high-throughput nature of this assay allowed for utilization of both biological replicates on multiple plates and technical replicates on the same plate. Generally speaking, correlation between growth curves in both technical and biological replicates was quite high (data not shown), emphasizing the power of this method. The results were consistent with previous shake flask experiments showing that furfural greatly

inhibited cell growth when the concentration was above 2 g/L (ca. 20 mM).

Growth rates (μ in terms of h−1) for each growth curve were then calculated as described previously (Franden et al., 2009, 2013). The response values, given as the percentage of the growth rate compared with the control in the absence of supplemented inhibitor, were then calculated for each concentration (**Table 4**). The response values were then used to determine the

concentration of inhibitor (IC50) that resulted in 50% growth compared to the control. The application of IC50 in this study was used to determine the top inhibitors in the slurries by comparing the IC50 of different compounds. This approach has frequently been used as a general toxicity indicator for potential inhibitors (Franden et al., 2013). IC50 values for 4-hydroxybenzaldehyde, vanillin, ferulic acid, 4-hydroxybenzoic acid, and vanillic acid were above the highest values tested, and the IC50 values for the remaining compounds were listed in **Table 5**. The growth of *C. necator* was inhibited by 5 out of the 11 compounds tested. The aldehydes, furfural and HMF, were the most toxic compounds in the slurry. The IC50 value of 9 mM for furfural was lower than the concentration detected in the saccharified slurry (12.6 mM or 1.2 g/L). Although the IC50 for benzoic acid had the lowest value, its concentration in the slurry was only 0.037 mM and less than 10% of the IC50 value making it the least toxic component in the slurry. The lignin degradation products were minimally toxic in the relevant concentration range except for p-coumaric acid with an IC50 value about 4 times higher than the concentration in the saccharified slurry. Ammonium acetate was found to be more toxic than ammonium sulfate (**Table 5**).

A few key points, however, must be made that aren't captured in **Table 4**. One major observation is that some of the compounds found in saccharified slurry stimulate growth by as much as 50% (**Table 6**). Growth stimulation with some of these compounds (e.g., ammonium sulfate, ammonium acetate, benzoic acid, and coumaric acid) was seen at low concentrations, but these compounds became inhibitory at higher concentrations. This indicated that *C. necator* might utilize these compounds as supplemental carbon or nitrogen sources at low levels. The utilization of acetate by *C. necator* in shake flask experiments (data not shown) supports this hypothesis that the higher response values could be an indicator of growth stimulation by utilizing the supplemented substrates (**Figure 2B**), which is consistent with previous report that addition of acetate as a supplementary substrate improved the cell growth and PHB production in *C. necator* DSMZ 545 (Sharifzadeh Baei et al., 2009). The combination of growth stimulation by some compounds and inhibition by others complicated analysis of the impact of saccharified slurry on growth and productivities and should be further explored in the future.

An additional point concerning the bacterial toxicity profiles is that they can provide additional information that could be used to understand microbial physiology and to propose genetic targets for metabolic engineering. For example, as the concentration of benzoic acid increased and the response changed from growth stimulation to inhibition, the color of the culture medium also became darker (**Figure 4**). This color change was not seen

in un-inoculated control wells. The availability of the *C. necator* genome sequence will facilitate our understanding of this phenomenon using a genomics approach that will be explained in the following section.

#### *C. NECATOR* **LIGNIN-DEGRADATION PATHWAY RECONSTRUCTION**

*Cupriavidus necator* is capable of utilizing lignin monomers and recently a *Cupriavidus* sp. (*C. basilensis* B-8) has been reported to be able to utilize kraft lignin (Shi et al., 2013). With that information in hand, we then set out to reconstruct the lignin degradation pathway of *C. necator* based on genome information and literature reports to help understand the bottleneck for further development of *C. necator* as a lignocellulosic biofuel-production strain (**Figure 5**).

Although the enzymatic reaction from lignin to lignin monomers and most reactions from lignin monomers to the metabolic intermediates are not yet fully revealed, homologous genes encoding the enzyme BenA, B, C, D were identified in *C. necator* that would convert benzoate to catechol. When the concentration of the supplemented benzoate is lower than the concentration characterized in the slurry, catechol will be further converted into cis,cis-muconate or 2-hydroxymuconate-6 semialdehyde through two different metabolic pathways leading to the accumulation of the key intermediate of acetyl-CoA for PHB production (**Figure 5**). Previous study already indicated that the accumulation of catechol in *Pseudomonas mendocina* can lead to its conversion of catechol into 1, 2-benzoquinone (Parulekar and Mavinkurve, 2006), which caused the culture medium to turn dark and at the same time inhibited the cellular growth since benzoquinone is a toxic agent against microorganisms such as *E. coli*, *Pseudomonas fluorescens,* and *Erwinia amylovora* (Beckman and Siedow, 1985). Our GC/MS study also indicated the correlation between the catechol appearance and benzoic acid disappearance when 1 g/L benzoic acid was supplemented into the minimum medium (**Figure 6**). Although further experimental data are needed to quantify the disappearance of benzoic acid and the appearance of catechol and 1,2-benzoqionone, our current study suggest that at high benzoate concentrations, *C. necator* can not convert catechol completely to cis,cis-muconate or 2-hydroxymuconate-6-semialdehyde and catechol began to accumulate which potentially caused the formation of 1, 2 benzoquinone and then turned medium into dark brown color and inhibited cell growth (**Figures 4**–**6**).

In addition, homologous genes were also identified in *C. necator* to carry out the reaction from p-coumarate to 4 hydrobenzoate and then feed into TCA cycle (**Figure 5**), which indicates that *C. necator* can utilize lignin monomer coumarate as carbon source for cell growth and may explain the stimulus effect of coumarate on *C. necator* growth in the lower concentration range. When high concentration of p-coumarate was supplemented into the medium, high concentration of catechol may be produced and the bottleneck reaction from catechol to cis,cis-muconate or 2-hydroxymuconate-6-semialdehyde (**Figure 5**) caused the inhibitory effect as discussed above. However, it is possible that other metabolic pathways that we have not been covered will also involved in and even play a key role on lignin monomer utilization and pre-treatment inhibitor sensitivity.

This study is an attempt and shows the feasibility to connect genomic information to explain microbial physiological phenomena. Although it could provide genetic targets for metabolic engineering to improve strain robustness (e.g., overexpression catechol dioxygenase to drive catechol into TCA instead of catechol accumulation), other experimental approaches such as systems biology study are needed to completely understand the mechanism of pre-treatment inhibitor sensitivity and genetic studies are especially required to confirm the hypothesis generated by bioinformatics study.

#### **AUTHOR CONTRIBUTIONS**

Wei Wang, Shihui Yang, Philip T. Pienkos, and David K. Johnson designed the experiments. Wei Wang carried out the flask assay and AC treatment. Shihui Yang carried out the Bioscreen C assay and genomic study. Glendon B. Hunsinger performed the GC/MS. Shihui Yang and Wei Wang did the spectrum scanning. Wei Wang, Shihui Yang, Philip T. Pienkos, and David K. Johnson analyzed the data, and wrote the manuscript.

# **ACKNOWLEDGMENTS**

This work was supported by the US Department of Energy, Bioenergy Technology Office (BETO) under contract number DE-AC36-08-GO28308 to NREL. NREL is a national laboratory of the US Department of Energy, Office of Energy Efficiency and Renewable Energy, operated by the Alliance for Sustainable Energy, LLC. We thank Michael Guarnieri and Thieny Trinh for preliminary toxicity experiments with *C. necator*, Ashutosh Mittal and William E. Michener for PHB and GC/MS analyses.

### **SUPPLEMENTARY MATERIAL**

The Supplementary Material for this article can be found online at: http://www*.*frontiersin*.*org/journal/10*.*3389/fmicb*.*2014*.* 00247/abstract

# **REFERENCES**

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 09 January 2014; accepted: 06 May 2014; published online: 27 May 2014. Citation: Wang W, Yang S, Hunsinger GB, Pienkos PT and Johnson DK (2014) Connecting lignin-degradation pathway with pre-treatment inhibitor sensitivity of Cupriavidus necator. Front. Microbiol. 5:247. doi: 10.3389/fmicb.2014.00247*

*This article was submitted to Microbial Physiology and Metabolism, a section of the journal Frontiers in Microbiology.*

*Copyright © 2014 Wang, Yang, Hunsinger, Pienkos and Johnson. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Elucidation of *Zymomonas mobilis* physiology and stress responses by quantitative proteomics and transcriptomics

# *Shihui Yang1,2,3\*, Chongle Pan4,5, Gregory B. Hurst 5, Lezlee Dice1,2, Brian H. Davison1,2 and Steven D. Brown1,2\**

*<sup>1</sup> Biosciences Division, Oak Ridge National Laboratory, Oak Ridge, TN, USA*

*<sup>5</sup> Chemical Sciences Division, Oak Ridge National Laboratory, Oak Ridge, TN, USA*

#### *Edited by:*

*Katherine M. Pappas, University of Athens, Greece*

#### *Reviewed by:*

*Harold J. Schreier, University of Maryland Baltimore County, USA Claudio Avignone-Rossa, University of Surrey, UK*

#### *\*Correspondence:*

*Shihui Yang, National Bioenergy Center, National Renewable Energy Laboratory, 15013 Denver West Parkway, Golden, CO 80401, USA e-mail<EMAIL_ADDRESS>Steven D. Brown, BioEnergy Science Center, Oak Ridge National Laboratory, 1 Bethel Valley Road, Oak Ridge, TN 37831, USA e-mail<EMAIL_ADDRESS>
*Zymomonas mobilis* is an excellent ethanologenic bacterium. Biomass pretreatment and saccharification provides access to simple sugars, but also produces inhibitors such as acetate and furfural. Our previous work has identified and confirmed the genetic change of a 1.5-kb deletion in the sodium acetate tolerant *Z. mobilis* mutant (AcR) leading to constitutively elevated expression of a sodium proton antiporter encoding gene *nhaA,* which contributes to the sodium acetate tolerance of AcR mutant. In this study, we further investigated the responses of AcR and wild-type ZM4 to sodium acetate stress in minimum media using both transcriptomics and a metabolic labeling approach for quantitative proteomics the first time. Proteomic measurements at two time points identified about eight hundreds proteins, or about half of the predicted proteome. Extracellular metabolite analysis indicated AcR overcame the acetate stress quicker than ZM4 with a concomitant earlier ethanol production in AcR mutant, although the final ethanol yields and cell densities were similar between two strains. Transcriptomic samples were analyzed for four time points and revealed that the response of *Z. mobilis* to sodium acetate stress is dynamic, complex, and involved about one-fifth of the total predicted genes from all different functional categories. The modest correlations between proteomic and transcriptomic data may suggest the involvement of posttranscriptional control. In addition, the transcriptomic data of forty-four microarrays from four experiments for ZM4 and AcR under different conditions were combined to identify strain-specific, media-responsive, growth phase-dependent, and treatment-responsive gene expression profiles. Together this study indicates that minimal medium has the most dramatic effect on gene expression compared to rich medium followed by growth phase, inhibitor, and strain background. Genes involved in protein biosynthesis, glycolysis and fermentation as well as ATP synthesis and stress response play key roles in *Z. mobilis* metabolism with consistently strong expression levels under different conditions.

**Keywords:** *Zymomonas mobilis***, microarray, proteomics and metabolomics, acetate, pretreatment inhibitor, stress responses, quantitative proteomics, systems biology**

### **BACKGROUND**

Yeast strains are among the leading current generation industrial biocatalyst microorganisms for fuel production (Hahn-Hagerdal et al., 2006). However, engineered bacteria such as *Zymomonas mobilis*, *E. coli, Bacillus subtilis* are also being developed and deployed to address commercial biofuel catalyst requirements (Dien et al., 2003; Inui et al., 2004; Romero et al., 2007; Alper and Stephanopoulos, 2009). *Z. mobilis* is an ethanologenic bacterium with many desirable industrial characteristics such as high-specific productivity and yield, high ethanol tolerance, and wide pH range (Gunasekaran and Raj, 1999; Panesar et al., 2006; Rogers et al., 2007). Recently, transformation efficiency has been improved by modifying the DNA restriction-modification systems (Kerr et al., 2010), and the inhibitor tolerance genes have been identified to improve the pretreatment inhibitor tolerance using genes from *Z. mobilis* (Yang et al., 2010a,b) or from *Deinococcus radiodurans* (Zhang et al., 2010). The genome sequences for strains ZM4, NCIMB 11163, 10988, 29291 and 29292 have been determined (Seo et al., 2005; Kouvelis et al., 2009, 2011; Pappas et al., 2011; Desiniotis et al., 2012), and the ZM4 genome annotation was improved recently (Yang et al., 2009a). Genome-scale *in silico* metabolic modeling analysis have been reported (Lee et al., 2010; Widiastuti et al., 2011; Rutkis et al., 2013) and recombinant strains have been engineered to express and secret cellulase (Linger et al., 2010) or ferment hexoses and pentose sugars such as xylose and arabinose (Zhang et al., 1995; Deanda et al., 1996).

A core challenge in next-generation biomass-based cellulosic biofuel is the recalcitrance of biomass to breakdown into sugars (Himmel et al., 2007; Alper and Stephanopoulos, 2009).

*<sup>2</sup> BioEnergy Science Center, Oak Ridge National Laboratory, Oak Ridge, TN, USA*

*<sup>3</sup> National Bioenergy Center, National Renewable Energy Laboratory, Golden, CO, USA*

*<sup>4</sup> Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN, USA*

Biomass pretreatment regimes are required to release the sugars, which can create inhibitors such as furfural, acetate, and vanillin (Almeida et al., 2007; Pienkos, 2009). The existence of pretreatment inhibitors increases production costs due to lower production rates and decreased yields. The development and deployment of robust inhibitor-tolerant biocatalysts for efficient fermentation of biomass to biofuel will be a critical component for successful production of biofuels at industrial-scale quantities to meet sustainability and energy security challenges associated with fossil fuels (Almeida et al., 2007). Several recent transcriptomics studies have used microarray approach to characterize mutant strain or *Z. mobilis* stress responses (Yang et al., 2009b, 2013; Hayashi et al., 2011; He et al., 2012a,b; Jeon et al., 2012; Skerker et al., 2013).

Acetic acid is an important inhibitor produced by the de-acetylation of hemicelluloses during biomass pretreatment. Unlike another major inhibitor furfural which is volatile and can be converted into a less toxic product of furfural alcohol (Liu et al., 2005; Heer and Sauer, 2008; Franden et al., 2009; Agrawal and Chen, 2011), acetate is stable during fermentation and poses a constitutive stress on the growth and ethanol production of *Z. mobilis* (Yang et al., 2010a,b). A *Z. mobilis* mutant strain, designated AcR, generated by chemical mutagenesis and selection, is able to produce ethanol efficiently in the presence of 20 g/L sodium acetate (NaAc), while the parent ZM4 is inhibited significantly above 12 g/L (Joachimstahl et al., 1998). Through comparative genome sequencing and next-generation sequencing (NGS)-based genome resequencing, we characterized the AcR mutant and identified a 1.5-kb deletion in strain AcR, which likely truncated the promoter of the *nhaA* gene encoding a sodium proton antiporter. We have carried out genetics study to confirm the association of 1.5-kb deletion in AcR mutant with its sodium acetate tolerance phenotype, we further performed microarray study to identify the differentially expressed genes between wildtype ZM4 and AcR mutant, and identified that *nhaA* gene is consistently upregulated in AcR mutant background compered to wild-type ZM4 under different conditions of NaCl and NaAc stress (Yang et al., 2010b).

Although we confirmed the 1.5-kb deletion in AcR mutant background leading to *nhaA* gene overexpression for enhanced sodium acetate tolerance phenotype (Yang et al., 2010b), we haven't systematically explored the global transcriptional profile difference between AcR mutant and wild-type ZM4 especially in the condition of minimal medium (MM), which is potentially more relevant to industrial fermentation conditions and more stressful to *Z. mobilis* than rich media (RM) we used for previous study (Yang et al., 2010b). In addition, there is no systems biology for this important industrial strain in minimal medium yet. To further explore the acetate-tolerance differences between ZM4 and AcR in MM condition, comprehensive microarraybased transcriptomic profiles and <sup>14</sup>*/*15N-labelled quantitative proteomic data were generated for ZM4 and AcR in MM with a large number of *Z. mobilis* proteins detected and quantified, which will be useful for future biocatalyst development. In addition, data collected in this and previous studies were used to hypothesize condition-responsive (strain, media, growth phase, or treatment) genes.

# **RESULTS**

# **PHYSIOLOGICAL RESPONSE OF** *Z. MOBILIS* **TO SODIUM ACETATE IN MINIMAL MEDIUM (MM)**

The growth of *Z. mobilis* wild-type ZM4 and AcR mutant in MM supplemented with 0, 12, or 16 g/L NaAc as well as 8.65 or 11.4 g/L NaCl with same corresponding Na+ molar concentrations as NaAc was assessed using a Bioscreen C instrument (Growth Curves USA, NJ) under anaerobic conditions to determine the effect of NaCl and NaAc on *Z. mobilis* growth and to decide on an appropriate NaAc concentration for subsequent systems biology studies. *Z. mobilis* grew more slowly in MM and attained lower cell densities (Additional File 1) compared to that in RM conditions (Yang et al., 2010b). Consistent with earlier RM results, wild-type ZM4 growth was arrested when NaAc was added to minimum medium at 16 g/L, and differences were observed between ZM4 and AcR with NaAc supplemented at concentrations of 12 g/L (Additional File 1). The concentration of 10 g/L was chosen for NaAc treatment in subsequent systems biology studies.

Similar to previous reports in RM (Joachimstahl et al., 1998; Yang et al., 2010b), mutant strain AcR also outperformed wildtype strain ZM4 in MM supplemented with NaAc (**Figure 1**, Additional File 1). Strain AcR overcame the acetate stress and reached stationary phase after about 130 h post-inoculation while ZM4 reached stationary phase after approximately 166 h. Both strains achieved similar final cell densities based upon OD600nm readings. Strain ZM4 had a longer lag phase and began to consume glucose much later than AcR, with most of the glucose in AcR cultures consumed at 148 h post-inoculation while cultures of ZM4 took 166 h to use the majority of the glucose in MM (**Figure 1**). As glucose was consumed earlier by AcR cultures, there was also a concomitant earlier onset for ethanol production compared to ZM4. Both strains produced similar yield of ethanol by the end of the experiment with acetate concentration kept steady during the whole experiment (**Figure 1**).

### *Z. MOBILIS* **WILD-TYPE ZM4 AND AcR MUTANT PROTEOMIC PROFILING DIFFERENCES IN MM**

Quantitative proteomics was used to compare the proteomic differences between *Z. mobilis* wild-type ZM4 and AcR mutant cells at 148 and 190 h post-inoculation. We identified 705 chromosomal encoded proteins and 6 plasmid derived proteins for the 148 h comparison, and 728 chromosomal proteins and 7 plasmid proteins for the 190 h comparison. We thus identified 638 chromosomal proteins and 4 plasmid proteins in common for the two time points. Altogether, this study identified 795 unique chromosomal proteins and 9 unique plasmid proteins, which is about 46% of the predicted *Z. mobilis* ZM4 proteome based upon a reannotation of the genome (Yang et al., 2009a) (Additional Files 2, 3A). The PI and MW distributions for the identified proteins were similar to that of the predicted theoretical PI and MW distributions for all proteins (Additional File 4). This study represents, to our knowledge, the largest quantitative measurement of the *Z. mobilis* proteome published to date.

The 148 h time point comparison between ZM4 and AcR strains revealed 120 proteins that had at least 1.5-fold differences in abundance and were significantly different based upon

the ProRata likelihood algorithm for quantitative shotgun proteomics and a confidence interval of 90% as described previously (Pan et al., 2006). This included 55 upregulated proteins and 65 downregulated proteins in AcR compared to ZM4, of which 11 were upregulated and 25 were downregulated more than 2-fold (Additional File 2). There were 107 proteins with a significant difference and at least a 1.5-fold change in levels for stationary phase comparison (190 h) using the same criteria. In strain AcR, 47 proteins were upregulated and 60 downregulated compared to ZM4, of which 10 were upregulated and 26 were downregulated and had changes of at least 2-fold (Additional File 2). Approximately half of the proteins that were identified as being different for one time point were also statistically different in the other time point (Additional Files 2, 3B–C).

The interactions among the down-regulated or up-regulated proteins with at least a 1.5-fold change were also analyzed for their previously documented interactions using STRINGS database (Jensen et al., 2009). The 14 proteins consistently upregulated in AcR for both time points compared to ZM4 had fewer interactions (Additional Files 2, 3D). Among the 39 proteins downregulated in AcR at both the 148 and 190 h time points post-inoculation, stress-responsive proteins such as catalase (ZMO0918), glutathione synthetase (ZMO1913), glutaredoxin 2 (ZMO0070), and a glutaredoxin-related protein (ZMO1873) were downregulated and most of them have been shown to interact with one another (Additional Files 2, 3E). Other proteins downregulated and connected included preprotein translocase (ZMO1896-8), cofactor synthesis of biotin synthase (ZMO0094), and amino acid biosynthesis such as shikimate 5-dehydrogenase (AroB, ZMO0041), threonine synthase (ThrC, ZMO1891), and 3-isopropylmalate dehydratase large subunit (LeuC, ZMO0105) (Additional Files 2, 3E).

Although our intention was to identify differentially expressed proteins more confidently by combining datasets from both time points than merely the second time point of stationary at 190 h to exclude the growth effect, the first comparison however will be potentially problematical due to the growth difference between AcR mutant and wild-type ZM4 at 148 h time point (**Figure 1**). Therefore, although comparison for second time point of stationary phase at 190 h is reasonably accurate, caution should be taken to interpret proteomic comparison result of second time point due to potential growth effect impact. To overcome this problem, samples from same growth phase between AcR mutant and ZM4 wild-type were used for our transcriptomic comparisons as discussed below.

### *Z. MOBILIS* **WILD-TYPE ZM4 AND AcR MUTANT TRANSCRIPTOMIC PROFILING DIFFERENCES IN MM**

Transcriptomic profiles for ZM4 and AcR in the presence of NaAc in MM were examined using NimbleGen high density expression arrays, essentially as reported previously (Yang et al., 2010b). The expression profiles generated in this study have been deposited in the GEO database (accession number GSE25443). About seventeen hundred genes were identified to be significantly differentially expressed using ANOVA modeling with strain (ZM4 and AcR) and time points post-inoculation as variables in MM (Additional File 5), which covered nearly all of the reannotated *Z. mobilis* ZM4 genes (Yang et al., 2009a). Hence dynamic gene expression changes were observed and for NaAc responses, 474 genes were significantly differentially expressed with at least a 2-fold change in expression values (Additional File 5).

Nine differentially expressed genes from different functional categories with a broad range of expression ratios were chosen for real-time quantitative PCR (RT-qPCR) validation (Additional File 6). RT-qPCR results indicated a high degree of concordance between microarray and RT-qPCR data with R-squared correlation coefficient values of 0.88, 0.87, 0.71, or 0.78 for the time points of 130, 148, 166, or 190 h respectively (Additional File 7).

Four genes were upregulated and 30 genes downregulated significantly with at least 2-fold changes in AcR when comparing AcR expression profiles to ZM4 profiles for all the time points (Additional File 5). Similar to the transcriptomics study in RM as reported previously (Yang et al., 2010b), the sodium/proton antiporter gene (ZMO0119) was also upregulated and ZMO0117 was downregulated in AcR in MM conditions. The other three genes that were consistently upregulated in AcR were Lasparaginase (ZMO1683), beta-fructofuranosidase (ZMO0375), and aldose 1-epimerase (ZMO0889). The 29 remaining genes that were consistently downregulated encoded mostly ribosomal proteins as well as proteins related to chemotaxis and flagellar biosynthesis, electron transport, fatty acid biosynthesis and protein export (Additional File 5).

### **CORRELATIONS BETWEEN PROTEOMIC AND TRANSCRIPTOMIC DATA OF ZM4 AND AcR IN MM**

We also examined the correlations between the gene or protein expression ratios at two time points of 148 and 190 h post-inoculation. There were 632 common proteins identified in both the 148 and 190 h time points. When the relative log2 ratios (AcR/ZM4) were compared, a moderate correlation (*R*<sup>2</sup> <sup>=</sup> <sup>0</sup>*.*47) was observed between the different sampling times (**Figure 2A**). Similarly, the correlation of transcriptomic profiles (AcR/ZM4) between 148 and 190 h was moderately high (*R*<sup>2</sup> <sup>=</sup> <sup>0</sup>*.*59) for the 530 common genes that showed statistically significant differential expression, but had not undergone any fold-change filtering (**Figure 2B**). When transcriptomic and proteomic data from the same time point (either 148 or 190 h) were compared, there are 460 gene-protein pairs for the 148 h time point that correlated poorly (*R*<sup>2</sup> <sup>=</sup> <sup>0</sup>*.*09) (**Figure 2C**), and the correlation between the 302 geneprotein pairs at 190 h post-inoculation is also relatively low (*R*<sup>2</sup> <sup>=</sup> <sup>0</sup>*.*17)( **Figure 2D**).

We have reported that the correlation between transcriptomic and proteomic data increases when only the significant genes or proteins were included for comparison, and the correlation even increases more when the ratios of transcriptomics or proteomics comparison increases from our recent *Z. mobilis* ethanol stress experiment of high-density microarray and shot-gun proteomics (Yang et al., 2013). Similarly, in this study, when only the genes or proteins with a significant and a ≥1.5-fold change were used for comparison purposes, the correlation between transcriptomic and proteomic data increased to *<sup>R</sup>*<sup>2</sup> <sup>=</sup> <sup>0</sup>*.*30 (*<sup>N</sup>* <sup>=</sup> 35) and *<sup>R</sup>*<sup>2</sup> <sup>=</sup> <sup>0</sup>*.*25 (*<sup>N</sup>* <sup>=</sup> 32) for 148 and 190 h comparisons, respectively (**Figure 2E**). Eighteen common gene-protein pairs were identified between proteomics and transcriptomics data in both 148 and 190 h post-inoculation with at least a 1.5-fold significant change, and all were downregulated in AcR (**Table 1**). The correlation for the proteomic study between 148 and 190 h went up to *<sup>R</sup>*<sup>2</sup> <sup>=</sup> <sup>0</sup>*.*79 (*<sup>N</sup>* <sup>=</sup> 18) when significance and fold change were considered (**Table 1**, **Figure 2F**). However, the correlation for the transcriptomic study between 148 and 190 h were decreased to 0.01 (*N* = 18) when significance and fold change criteria were applied (**Table 1**, **Figure 2F**). This result may suggest that change at transcriptional level is rapid which is consistent with our transcriptomic study in *Clostridium thermocellum* that genes were differentially expression even after several minutes post ethanol or furfural shock (Yang et al., 2012; Wilson et al., 2013). Moreover, this study also suggested that the change at protein level is relatively steady with potential post-transcriptional modification existing.

### **METABOLIC PATHWAY ANALYSES OF PROTEOMIC AND TRANSCRIPTOMIC DATA OF ZM4 AND AcR IN MM**

The PathwayTools Omics Viewer (Karp et al., 2002, 2011) was used to further examine transcriptomics and proteomics data and their relationships. Genes in the fatty acid and hopanoid biosynthetic pathways are upregulated in ZM4 compared to AcR (Additional File 8A). Hopanoids have been found in a variety of bacteria including *Z. mobilis* and are reported to protect against the toxic effects of ethanol (Flesch and Rohmer, 1989; Hermans et al., 1991; Horbach et al., 1991; Shigeri et al., 1991; Welander et al., 2009). For example, they have also been reported to play a role in *Rhodopseudomonas palustris* membrane integrity and pH homeostasis (Welander et al., 2009). However, there is debate of their role in *Z. mobilis* (Moreau et al., 1997) and we didn't identify obvious differential gene expression when *Z. mobilis* encountered ethanol stress (Yang et al., 2013). Further study is required to elucidate their role in *Z. mobilis*, and their relative contribution, if any, to sodium acetate sensitivity or tolerance to the strains used in this experiment.

In other pathway comparisons of AcR to ZM4 at 148 and 190 h, the strains had similar patterns of protein detection, with 120 or 107 differentially expressed proteins with ≥1.5-fold change at each time point respectively were used for analyses (Additional Files 8B,C). We identified 509 chromosomal gene/protein pairs that have no significant changes between AcR and ZM4 in MM with NaAc in both 148 and 190 h post-inoculation (Additional File 2), the majority of which are involved in the ED and mixedacid fermentation pathways, and tRNA charging (Additional File 8A–C). When the log2 transformed peptide hits of all the proteins identified in 148 and 190 h post-inoculation were used for pathways analyses, proteins involved in ED pathway and ethanol production are among the most abundant proteins (Additional Files 2, 8D). The abundance of ED pathway and ethanol production enzymes and the relatively stability of protein levels between ZM4 and AcR indicates that this pathway for carbon and electron flow are core elements in the physiology of both strains.

Recently, a number of studies have added important details to our understanding of *Z. mobilis* physiology (Kouvelis et al., 2009; Yang et al., 2009a, 2010a,b; Kerr et al., 2010; Linger et al., 2010; Zhang et al., 2010; Widiastuti et al., 2011). However, many genes encode proteins for which we now have proteomics evidence of their expression, their functions remain unknown. The largest quantitative measurement of the *Z. mobilis* proteome obtained in this study may help elucidate *Z. mobilis* genes, proteins, their functions and regulation.

### **ZM4 AND AcR TRANSCRIPTOMIC PROFILES IN DIFFERENT CONDITIONS**

In this study, we also combined *Z. mobilis* microarray data collected from previous studies with expression data collected in this experiment to provide greater insights in *Z. mobilis* physiology and gene regulation. Forty-four microarrays were used in the analysis (Additional File 9) to compare the gene expression differences between the variables of strain (AcR and ZM4), growth phase (exponential and stationary), media (RM and MM), and treatment (NaCl, NaAc, and control of RM only). Considering the growth differences between ZM4 and AcR in this

experiment (**Figure 1**) and the fact that *Z. mobilis* has dramatic transcriptional differences between exponential and stationary phase (Additional Files 5, 10) (Yang et al., 2009b), the relative large number of genes significantly differentially expressed at time point 148 h post-inoculation between AcR and ZM4 as well as the low transcriptomic correlation between 148 and 190 h post-inoculation is likely associated with the growth phase differences. To ensure compatible comparisons were being made among different time points in RM or MM under different conditions, the correlations among all the forty-four microarrays were investigated using the JMP Genomics 4.0 (SAS, NC) and based upon a hierarchical clustering and growth curve analyses (**Figure 3**), time points were assigned to either exponential or stationary phase for further statistical analyses (**Figure 3**, Additional File 9). Differences between the variables of strain, growth phase, media, and treatment were then analyzed by ANOVA, as described previously (Yang et al., 2009b, 2010b). The ANOVA identified that nearly every gene in the *Z. mobilis* genome was differentially expressed under one or more of the many conditions tested in this global analysis (Additional File 11), leading to an opportunity to identify the strain, media, growth phase or treatment responsive genes, as discussed below.

**Table 1 | Gene-protein pairs identified in both transcriptomic and proteomic studies at both 148 and 190 h post-inoculation, and those only identified in either 148 or 190 h post-inoculation.**

*(Continued)*

#### **Table 1 | Continued**

*A148 and A190: array result comparing AcR to ZM4 at time point 148 and 190 h post-inoculation respectively; P148 and P190: proteomics result comparing AcR to ZM4 at time point 148 and 190 h post-inoculation respectively. All the gene-protein pairs are statistically significant with at least 1.5-fold difference between AcR and ZM4 in both transcriptomic and proteomic studies. Negative number indicates gene/protein was downregulated comparing AcR to ZM4 and positive number indicates upregulation.*

#### *Strain-specific genes*

Gene expression profiles for ZM4 and AcR were compared and when all the conditions (Additional File 9) were taken into account, only five genes were significantly differentially expressed with at a least 2-fold change between strains. The significant differentially expressed strain specific genes included two genes with increased expression and three genes downregulated in strain AcR (Additional Files 11, 12). We have previously reported that *nhaA* (ZMO0119) is upregulated and ZMO0117 is downregulated in AcR compared to ZM4 in exponential and stationary growth conditions for RM (Yang et al., 2010b). In this study, which considered more conditions, we identified that that another gene (hypothetical protein ZMO1787, 82aa) was consistently upregulated in AcR under the conditions tested and two AcR genes encoding a predicted permease (ZMO0055, 263 aa) and a conserved hypothetical protein (ZMO0025, 234 aa) were consistently downregulated Additional File 11).

#### *Media-responsive genes*

Gene expression differences between MM and RM in the presence of NaAc were compared (Additional File 9). When the strain responses of both AcR and ZM4 in both exponential and stationary phases were combined to examine media differences, 232 genes were identified as being upregulated and 247 downregulated in MM compared to RM (Additional File 13). Almost half of upregulated genes encode hypothetical proteins with unknown functions, and many of the remaining upregulated genes encode proteins involved in stress sensing and responses. These included two-component signal transduction genes (7), transcriptional regulators (11), sigma factors (*rpoD*, *rpoH*), flagellar and chemotaxis genes, and other stress responsive genes such as an putative operon containing heat-inducible transcription repressor *hrcA* (ZMO0015) and chaperone gene *grpE* (ZMO0016). Other genes included cold shock gene (ZMO0925), DnaJ-class chaperone (ZMO1069) and chaperone *hspD* (ZMO0989), three glutaredoxin-related genes, three thioredoxin-like genes, and genes associated with DNA repair such as UvrABC system gene *uvrB* (ZMO0362), *mutM* (ZMO1187), *mutS* (ZMO1907) and ZMO1426 (Additional File 13).

There are more previously documented interactions among the downregulated genes in MM relative to RM compared to the upregulated genes (Additional File 14). These MM-downregulated genes are mostly related to the central carbon metabolism, which included genes such as 6-phosphogluconolactonase gene *pgl* (ZMO1478), phosphoglycerate mutase gene *pgmA* (ZMO1240), gluconolactonase gene *gnl* (ZMO1649), pyruvate decarboxylase gene *pdC* (ZMO1360), and pyruvate dehydrogenase gene *pdhA* and *pdhB* (ZMO1606, ZMO1605) (Additional Files 13, 14B). Genes involved in amino acid biosynthesis, encoding ribosomal proteins and genes for nucleotide biosynthesis, as well as genes related to energy metabolism such as electron transport system and ATP generation were also mainly down-regulated in MM (Additional Files 13,14B).

### *Treatment-responsive genes*

The following conditions were used to examine treatmentspecific responses: ZM4 grown in RM; both ZM4 and AcR growth in RM with NaCl, RM with NaAc, or MM with NaAc (Additional File 9). These conditions were used to make the following comparisons: treatment of NaCl or NaAc in RM vs. RM for strain ZM4; and as the treatment of NaAc vs. NaCl in RM for both ZM4 and AcR strains (Additional File 13).

The NaCl-responsive genes in the ZM4 wild-type background included 88 downregulated genes and 47 upregulated genes (Additional Files 13, 15A,B). More genes were responsive in a ZM4 background grown in RM supplemented with NaAc, including 159 downregulated and 103 upregulated genes. This was approximately twice the number as responsive to the less severe NaCl treatment (Additional Files 13, 15C,D). Differentially expressed genes shared between NaCl and NaAc in RM compared to RM included 27 upregulated and 66 downregulated genes (Additional Files 13, 15E,F). Genes upregulated in both NaCl and NaAc treatment related to cysteine synthesis, with 6 genes clustered together (*cysD, I, J, K, N*, ZMO0006 and ZMO0055) and tryptophan synthesis (*trpB, F*), and most of the downregulated genes encoded hypothetical proteins with unknown function with the exception of several related to flagellar biosynthesis and sucrose metabolism (Additional Files 13, 15).

When the responses of *Z. mobilis* in RM with NaAc were compared to that of NaCl, there were 37 upregulated and 41 downregulated genes (Additional Files 13, 15I,J). Genes related to ribosomal proteins and amino acid biosynthesis were downregulated and upregulated genes related to the Entner–Doudoroff (ED) pathway (ZMO1518 and *pgi*), energy metabolism [e.g., electron transport (ZMO0021, ZMO1851, ZMO1885) and ATP synthesis gene *atpC* (ZMO0242)], cell wall formation (ZMO1724), regulator gene *zrp* (ZMO0372), and several transporter related genes such as signal peptidase I gene *lepB* (ZMO1710), Sec-independent protein translocase tatA/E homolog gene (ZMO1220), TolQ biopolymer transport gene (ZMO0161) and biopolymer transport gene *exbD* (ZMO1715), Fe2<sup>+</sup> transport system gene *feoB* (ZMO1541), and potassium transport system gene *kup* (ZMO1209).

#### *Growth phase-dependent genes*

When strain profiles in the presence of NaAc were considered together and a ≥2-fold change was applied, 331 genes were significantly differentially expressed in a phase-dependent manner for RM (Additional Files 11, 12). When the RM NaCl profiles were analyzed, 661 genes were identified (Additional Files 11, 12). Similar patterns were identified for the comparisons when variable of strain was taken into consideration separately (Additional Files 11, 12). In the presence of NaAc, the stationary phase upregulated genes were less correlated than the downregulated ones (Additional Files 16). The genes downregulated in stationary phase are related to ribosomal proteins, chemotaxis and flagellar systems, amino acid and nucleotide biosynthesis, and electron transport (Additional Files 11, 16A). Except for those encoding hypothetical proteins, genes upregulated in stationary phase are related to stress responses such as the catalase gene (ZMO0918), glutaredoxin gene (ZMO0070), thioredoxin gene (ZMO1705), ferredoxin gene *fdxN* (ZMO1818), organic hydroperoxide resistance gene (ZMO0693), ATP-dependent Clp protease gene *clpB* (ZMO1424), RNA polymerase sigma-32 factor gene *rpoH* (ZMO0749), and integration host factor gene *ihfB* (ZMO1801) (**Additional Files 11, 16B**).

# **DISCUSSION**

## **GENES CONTRIBUTING TO SODIUM ACETATE TOLERANCE PHENOTYPE OF AcR**

AcR mutant had advantages under NaAc stress in both RM (Yang et al., 2010b) and MM conditions, with shorter lag phase leading to earlier glucose consumption and ethanol production than ZM4 cells (**Figure 1**). We have reported previously that the truncation within gene ZMO0117 in AcR causes the consistent upregulation of the sodium proton antiporter gene *nhaA* (ZMO0119), which is the determining factor of AcR mutant for sodium acetate tolerance in RM (Yang et al., 2010b). In this study, the *nhaA* gene also had higher expression levels in strain AcR compared to ZM4 in MM (Additional File 11). ZM4 *nhaA* expression in unamended RM was higher (∼2-fold) in stationary phase compared to exponential phase (Additional File 11) and differential expression was also measured for other conditions and strain comparisons (Additional File 13). Acetate accumulates in stationary phase *Z. mobilis* fermentations (Yang et al., 2009b) and greater expression of *nhaA* gene may ameliorate cellular function and activity under such stressful conditions.

In addition, compared to ZM4, the ZMO0117 gene in AcR showed a similar profile of down-regulation in MM (Additional File 11) to that of in RM (Yang et al., 2010b). Proteomic data also showed that the peptide hits for ZMO0117 were less abundant in AcR compared to ZM4 at 148 h post-inoculation (Additional File 2). However, the NhaA protein was not detected in any sample in this experiment (Additional File 2) and yet differential expression was observed for this gene, which demonstrates the value of conducting integrated experiments and may point to a technical challenge in detecting hydrophobic proteins like NhaA that have many transmembrane spanning domains (Poetsch and Wolters, 2008; Gilmore and Washburn, 2010).

As discussed above, this study also identified other AcR specific genes. For example, ZMO0055 was significantly downregulated in AcR compared to ZM4 in all conditions with NaAc and responded to other conditions too (Additional File 11). ZMO0055 is predicted to belong to the sulfite exporter TauE/SafE family protein (pfam01925), a family of integral membrane proteins involved in transporting anions across the cytoplasmic membrane during taurine metabolism as an exporter of sulfoacetate (Weinitschke et al., 2007). As discussed above, ZMO0055 is closely related to CysC, D, I, N (ZMO0003, 4, 5, 8) based on protein analysis using String pre-computed database (Additional File 15), which were upregulated in both NaCl and NaAc treatment. In addition, recent work by Skerker et al. (2013) also indicated that cysteine synthases, which are required for L-cysteine biosynthesis from sulfate, are related to hydrolysate tolerance. They proposed that the increased demand for sulfite and cysteine is due to the increased needs of glutathione, which is formed from glutamate and cysteine, during *Z. mobilis* growth in stressful condition (Skerker et al., 2013). Therefore, the downregulation of sulfite exporter potentially could help Z. *mobilis* retain sulfite for enhanced sodium acetate tolerance phenotype. However, further detailed study is needed to understand the role and contribution of these AcR specific genes, which will help reveal the complete difference between AcR and ZM4 and at the same time help assign the function for these AcR-specific genes with hypothetical function (ZMO1787 and ZMO0025).

#### **RELATIONSHIP AMONG MEDIA, INHIBITOR, GROWTH PHASE, AND GENE EXPRESSION DYNAMICS**

The gene expression intensity mean value of the significantly expressed genes (Additional File 17) were further used for hierarchical clustering to compare each gene at 14 different conditions, the result indicated that minimum medium had the most dramatic effect on gene expression compared to those in RM, followed by growth phase, inhibitors (NaAc*>*NaCl), and strain (Additional File 18A). Gene expression profiles in MM were clustered together and separated from the rest in RM (Additional File 18A). *Z. mobilis* growth in MM has a longer lag phase and generates a lower final cell density compared to that in RM (Additional File 1) (Yang et al., 2010b), which indicates that one or more factors are limiting *Z. mobilis* growth in MM, and MM is possibly a more stressful environment for *Z. mobilis* cells which is further supported by the upregulation of stress-responsive genes in MM compared to that in RM (Additional Files 13, 14).

The second level of clusters was between exponential and stationary growth phases, consistent with an earlier global profiling report for strain ZM4 that growth phase tends to have more dramatic effect on gene expression than a single stressor (Yang et al., 2009b). The greatest number of differentially expressed genes between exponential and stationary phases is observed for strain ZM in RM, followed by strain ZM in RM with NaCl, RM with NaAc and then ZM4 in MM (Additional File 12). In RM, *Z. mobilis* growth is inhibited by NaCl but NaAc is more inhibitory (Yang et al., 2010b), and a similar trend was observed in MM (Additional File 1). In this experiment, the inhibitor (NaAc) was added into the medium before the fermentation allowing the strains to respond and adapt to the conditions for several generations before samples were taken for systems biology studies. Hence, when transcriptomic and proteomic profiles between strains are compared under the same growth condition, the differences likely represent differences in homeostasis that allow AcR to function better than ZM4. The two strains were in different growth phases for the first proteomics time point, which confounds strain comparisons for tolerance but does still allow gene-protein relationships to be examined. In general, the large datasets reported in this study will allow others to investigate aspects such as codon bias and RNA secondary structures.

In addition, the microarray results for each gene at different conditions can help us understand the microbial physiology and identify consistently strong or inducible promoters for metabolic engineering application. For example, in this study, 52 genes and 6 other genetic features were among the top 2.5% of all genetic features with consistently strong expression intensity. Except for a dozen of genes encoding hypothetical proteins, they involve in protein biosynthesis, glycolysis and fermentation as well as ATP synthesis and stress response indicating their key role on *Z. mobilis* metabolism (Additional File 17). Several papers have been published recently to study the stress responses of *Z. mobilis* to different inhibitors (e.g., ethanol, furfural) using systems biology approaches (Yang et al., 2009b, 2010b, 2013; He et al., 2012a,b; Jeon et al., 2012; Skerker et al., 2013), we are currently working on several transcriptomics studies using microarray and next-generation sequencing based strand-specific RNA-Seq technique. Combining all these systems biology datasets, we will revisit this topic to investigate the impact of different variables of media, carbon source (glucose or xylose), different inhibitors, and growth phase on *Z. mobilis.* The transcriptional profiles will be further compared to help identify the condition specific promoters more confidently, and the discrepancy between transcriptomics and proteomics at different conditions will help understand the post-transcriptional regulation mechanism.

# **CONCLUSIONS**

Our study has provided the global profiling of the model ethanogenic bacterium *Z. mobilis* at both transcriptomic and proteomic levels in minimal medium (MM) for the first time. The results indicated that AcR had similar advantage in MM as in RM. AcR mutant can overcame acetate stress earlier with shorter lag phase, earlier glucose utilization and ethanol production than wild-type ZM4, and sodium proton antiporter gene (*nhaA*) also plays an important role in sodium acetate tolerance in MM as that of in RM we reported previously (Yang et al., 2010b). In addition, this is also the first attempt that we are aware of to combine massive transcriptomic data for condition-specific gene identification. The proteomic and transcriptomic data generated in this study and the one we have reported will provide massive datasets for future metabolic modeling and strain improvement.

# **MATERIALS AND METHODS**

#### **STRAINS AND GROWTH CONDITIONS**

Wild-type *Z. mobilis* ZM4 was obtained from the American Type Culture Collection (ATCC 31821). *Z. mobilis* acetate tolerant strain AcR has been described previously (Joachimstahl et al., 1998). ZM4 and AcR were cultured in RM (Glucose, 20.0 g/L; Yeast Extract, 10.0 g/L; KH2PO4, 2.0 g/L, pH5.0) (Yang et al., 2009b) for routine strain maintenance. MM was similar to the one reported to isolate auxotrophic *Z. mobilis* mutants (Goodman et al., 1982): 20 g Glucose, 1 g KH2PO4,1gK2HPO4, 0.5 g NaCl, 1 g (NH4)2SO4 was dissolved in 986.5 mL distilled H2O with pH adjusted to pH 5.0 before autoclave. The sterilized MM broth was then cooled to below 55◦C before 500 mg MgSO4.7H2O (2.5 mL 200 g/L stock), 25 mg Na2MoO4.2H2O (1 mL 25 g/L stock), and 10 mL Vitamin (ATCC MD-VS) were added to make a 1-L MM.

Growth assays were used to identify the effects of sodium acetate (NaAc) and sodium chloride (NaCl) on the growth of *Z. mobilis* in MM for subsequent fermentation experiments using a Bioscreen C instrument (GrowthCurves USA, NJ). Fermentations were conducted in 7.5-L BioFlo110 bioreactors (New Brunswick Scientific, NJ) fitted with agitation, pH, and temperature probes and controls, and bacterial growth was monitored turbidometrically by measuring optical density at 600nm with a model 8453 spectrophotometer (Hewlett-Packard, CA.) as described previously (Yang et al., 2009b), except that the fermentation volume was 4 L. Samples were harvested during fermentation at different time points (**Figure 1**) as described previously (Yang et al., 2009b).

#### **HPLC**

HPLC analysis was used for the measurements of the concentration of glucose, acetate, and ethanol in 0.2µm-filtered samples taken at different time points during fermentation (**Figure 1**) and analyzed as described previously (Yang et al., 2009b).

#### **PROTEOME SAMPLE PREPARATION**

Duplicate mixtures of microbial cells that were metabolically labeled with either 15N ammonium sulfate for the wild-type *Z. mobilis* ZM4 culture or 14N ammonium sulfate for the AcR mutant culture were prepared by mixing equal weights of cell pellets from duplicate cultures for each strain. Cell mixtures were lysed by sonication in ice-cold 50 mM Tris-HCl (pH 7.5) buffer, and unbroken cells were removed by centrifugation at 5000 × g for 10 min. Protein concentration for each sample was determined with the RC DC™ protein assay (Bio-Rad Lab, CA). The two fractions from each cell mixture were digested using the following protocol. The proteins were denatured and reduced with 6 M guanidine and 10 mM dithiothreitol (DTT) (Sigma-Aldrich, MO) at 60◦C for 1 h. The denatured proteome fractions were diluted 6-fold with 50 mM Tris/10 mM CaCl2 (pH 7.6), and sequencing grade trypsin was added at the ratio of 1:100 (wt:wt). The first digestion was run overnight at 37◦C and, after adding additional trypsin, the second digestion was run for 5 h at 37◦C. The samples were then reduced with 20 mM DTT for 1 h at 60◦C and were desalted using Sep-Pak Plus C-18 solid-phase extraction (Waters Co, MA).

## **QUANTITATIVE PROTEOMICS MEASUREMENT**

The protein digests were examined with LC-MS/MS using twelvestep split-phase MudPIT (MacCoss et al., 2002; McDonald et al., 2002) in duplicate. The samples were loaded via a pressure cell (New Objective, MA) onto a 250-*u*m-I.D. fused silica front column fritted into an M-520 filter union (Upchurch Scientific, WA). The column packing consisted of 2 cm strong cation exchange resin Luna® and 2 cm C18 reverse-phase resin Aqua (Phenomonex, CA). A 100-*u*m-I.D. PicoFrit column (New Objective, MA) was packed with 15 cm C18 reverse-phase resin. The front column was connected with the PicoFrit column and then placed in-line with a Dionex Ultimate quaternary HPLC. Two-dimensional LC separation was performed with twelve salt pulses, each of which was followed by a 2-h reverse-phase gradient. MS/MS analysis was performed on an LTQ linear ion trap instrument (ThermoFinnigan, CA) with dynamic exclusion enabled. Each full scan (400–1700 *m/z*) was followed by three data-dependent MS/MS scans at 35% normalized collision energy. The full scans were averaged from five microscans and the MS/MS scans were averaged from two microscans.

#### **QUANTITATIVE PROTEOMICS DATA ANALYSIS**

All MS/MS scans were searched in two iterations against the FASTA database containing all annotated *Z. mobilis* proteins using the SEQUEST program (Eng et al., 1994). In the first iteration, the molecular masses of amino acids containing 14N were used, and, in the second iteration, the masses of amino acids containing 15N were used. The peptide identifications from the two iterations were merged. The DTASelect program (Tabb et al., 2002) was used to filter the peptide identifications and to assemble the peptides into proteins using the following parameters: retaining the duplicate MS/MS spectra for each peptide sequence (DTASelect option: -t 0), fully tryptic peptides only, with a delCN of at least 0.08 and cross-correlation scores (Xcorrs) of at least 1.8 (for parent ion charge state, *z* = +1), 2.5 (*z* = +2), or 3.5 (*z* = +3). Selected ion chromatogram extraction, peptide abundance ratio estimation and protein abundance ratio estimation were completed with the ProRata program as described previously (Pan et al., 2006, 2008).

#### **MICROARRAY ANALYSIS AND qRT-PCR VALIDATION**