Systematic Discovery of Archaeal Transcription Factor Functions in Regulatory Networks through Quantitative Phenotyping Analysis

Identification of transcription factor candidates.To prioritize candidate TFs for phenotypic characterization, genes encoding TFs were first detected in the Halobacterium salinarum genome (27) using the Systems Biology Experimental Management System database (SBEAMS) (36), where annotations were cross-referenced with other databases (PFAM, COG [clusters of orthologous groups of proteins], and PSI-BLAST [37–39]) (Fig. 1). Previous ab initio protein structure prediction results (26) and subsequent matches to the Protein Data Bank (PDB) were used to identify possible DNA binding domains in proteins of unknown function. This pipeline resulted in a list of 130 putative TFs, which were included as priors in inference of the H. salinarum EGRIN model (7, 15). To identify candidates of interest for further analysis here, TFs were more stringently defined as those with strong homology (E value ≤ 1 × 10⁻⁷) to a sequence-specific DNA binding domain or structural fold experimentally characterized in other bacterial or archaeal species (13, 26) (Fig. 1). This definition excluded 24 DNA-binding proteins with peripheral roles in transcription (e.g., helicases). Another 18 genes encoding previously published, well-characterized TFs, were also excluded (21, 40–42). Of the 88 remaining TFs, a final subset of 27 were selected for further study based on transcriptional changes during fluctuations in a wide array of environmental conditions (1, 7, 20, 23–25, 43) and functional predictions from EGRIN gene regulatory network models (7, 15) (Fig. 1; see Table S1 in the supplemental material). The TFs in the resultant collection are members of a variety of functional families and contain diverse structural domains (Table S1). These include DNA binding domains known from other archaea and bacteria (winged helix-turn-helix, ribbon-helix-helix) and domains unique to halophilic archaea (e.g., HalX). Some TFs contain ligand binding domains homologous to those of known function in bacteria, such as DtxR family iron-dependent repressors. Other TFs possess domains of novel function specific to halophiles, such as the RosR C-terminal domain (16) (Table S1). Together, the results of these bioinformatic analyses suggest that the panel of selected TFs is representative of the global transcription regulatory landscape of H. salinarum.

• Open in new tab
• Download powerpoint

FIG 1

TF candidate selection pipeline. Genes encoding proteins with a putative DNA binding domain were annotated using sequence databases PFAM (37) and PSI-BLAST (38), structural predictions (26), and protein functions from COG (39). These annotations were stored in the Systems Biology Experiment and Analysis Management System (SBEAMS) database (36), resulting in 130 putative TFs. Transcriptome analysis across 1,495 experimental conditions (43) and GRN network inference models (7, 15) were then used to generate predictions regarding TF functions. Details of GRN predictions and criteria for selection of the final collection of 27 TFs are given in Table S1 in the supplemental material. DBD, DNA binding domain.

Quantification of significant differential growth of TF knockout strains under stress. (i) Experimental design and data.To determine the physiological function for each TF and systematically compare these functions across TFs, knockout mutant strains for each of the 27 selected TFs were grown under five conditions: standard growth, low salinity, paraquat (PQ), hydrogen peroxide, and heat shock (see Materials and Methods). The growth of these knockout mutant strains was compared to that of the isogenic parent control strain, Δura3 strain, a uracil auxotroph used to generate knockout mutants (44) (see Materials and Methods). Growth conditions were chosen based on their relevance to the hypersaline habitat of H. salinarum. Knockout mutants for 10 of these 27 TFs were constructed in previous studies, where the role of each TF was assessed under a single stress condition (Table 1). These strains were included for the following reasons: (i) as a control to validate our methods and (ii) to test possible secondary functions for these previously studied regulators. Strains with deletion of the genes encoding the remaining 17 TFs of interest were constructed in the current study using established genetic methods for H. salinarum (see Materials and Methods) (44) (Table S2). The genotypes of all 27 mutants, regardless of prior publication, were also verified here (see Materials and Methods; Fig. S1 and Table S2). Growth of each knockout and parent control strain was measured under each of the five conditions every 30 min for 48 h, resulting in 210,180 data points (raw data given in Table S3).

(ii) Development of a FANOVA model and test for differential growth of knockout mutants under stress.Typically, microbial population growth is modeled using parametric models, such as Gompertz regression (45). The effects of stress and genetics on growth are then quantified by testing for statistically significant differences in the estimated parameters for different conditions (46). However, we previously showed that the impact of genetics and stress perturbations on growth are more accurately captured with a nonparametric Gaussian process (GP) model (35). In contrast to parametric models, GPs have the advantage of learning these relations directly from the data and do not require explicit equations describing the effects of stress and genetics on growth, which are typically unknown in the case of phenotypic discovery. In this study, we again use GPs to model microbial population growth, now in the form of functional analysis of variance (FANOVA) (47). In the FANOVA setting, microbial population growth is divided into a linear combination of different experimental variables: condition, genetic background, and the interaction between the two (see Materials and Methods) (Fig. 2A). Additionally, FANOVA allows for the explicit comparison of different effects in terms of significance, even if the shapes of those effects are quite different. For these reasons, FANOVA is particularly well adapted to modeling and comparing the effects of many different stress and genetic perturbations. This extends our previous GP model (35) by including all data from multiple conditions and genetic backgrounds into a single model, which allows for direct comparison of different effects in a common framework.

• Open in new tab
• Download powerpoint

FIG 2

FANOVA modeling and statistical ranking of TF knockout mutant growth phenotypes across five environmental conditions, standard growth conditions, paraquat stress, peroxide, low salt, and heat shock. (A) Mutants with the largest difference in growth compared to the Δura3 control strain under each condition. Raw data from growth trajectories for individual cultures under standard conditions (thin gray lines) were fit and compared with growth under stress conditions (thin blue lines) using FANOVA (see Materials and Methods). Solid lines indicate the mean of the fit to all replicate trajectories, and shaded regions are the 95% confidence interval of the fit. (B) Functional difference (OD_Δ) of the TF knockout strain relative to the isogenic Δura3 parent strain. The mutants with the top three scoring phenotypes according to ||ODΔ|| are shown. (C) Summary statistical metric ranking (||ODΔ||) for those mutants with strongly different growth trajectories compared to that of the Δura3 strain (||ODΔ||≥ 0.337 across all conditions).

Two metrics were used to assess the significance of differential growth based on our FANOVA model of growth data. The first metric, change in the optical density (OD_Δ) (35), is a function representing the difference in growth levels between the mutant strain [f_m(x)] and parent strain [f_p(x)] over the duration of the growth curve under a specific condition (Fig. 2B and Fig. S2) as follows: ODΔ(t)=fm(t)−fp(t)(1)

Where the OD_Δ curve differs significantly from zero (based on a 95% credible interval), the growth phenotype is considered significantly different from the parent strain for that section of the growth curve. The second metric, ||ODΔ||, represents the overall magnitude of the functional difference between parent and mutant strains under a specific condition and is calculated based on OD_Δ (Fig. 2C and Fig. S3): ||ODΔ||=∑t=t0tn[ODΔ(t)]2(2)

Higher values of ||ODΔ|| indicate an overall larger deviation from parent strain behavior, regardless of positive or negative growth phenotype, and were used as a rank ordering of phenotype severity to prioritize TF mutants for further analysis.

(iii) FANOVA modeling of microbial population growth enables the discovery of new TF knockout phenotypes.On the basis of ||ODΔ|| ranking, the ΔtrmB strain exhibited the strongest growth impairment under standard conditions (Fig. 2A). This result was expected based on previous characterization of TrmB as a global regulator of nutrient metabolism (17, 18, 35, 48, 49). Also as expected from prior work, the ΔrosR mutant with deletion of a key oxidative stress response TF was among the top-ranking quartile of mutants under oxidative stress (16, 35, 50). Consistent with prior Gaussian process model results, the ΔsirR strain exhibited a significant growth impairment under peroxide stress (35); indeed, here the ΔsirR strain exhibited the top-ranking peroxide phenotype of all mutants tested (Fig. 2C). The consistency of these results with those of prior studies demonstrates the validity of the FANOVA model for recapitulating known phenotypes.

Surprisingly, however, several novel mutant phenotypes were observed. For instance, the ΔcspD1 mutant showed the strongest growth impairment relative to the parent control strain under PQ stress, and the ΔcopR mutant exhibited the strongest growth impairment under heat shock (Fig. 2A and B). In order to compare the number of significant mutant phenotypes under each condition, a universal cutoff of ||ODΔ|| of 0.337 was chosen (see Materials and Methods). This cutoff allowed rank ordering of the stress conditions themselves in terms of the strength of perturbation to physiology, as measured by the number of significant mutant phenotypes. Peroxide treatment had the strongest effect, with 19 mutants exhibiting significant differences in growth compared to the Δura3 strain (12 mutants grew more slowly than the Δura3 strain, and 7 mutants grew faster than the Δura3 strain [Fig. S2 and S3]). The next strongest effect was low salt (11 mutants, 9 faster and 2 slower than the Δura3 strain), followed by PQ (10 mutants, 2 faster and 8 slower than the Δura3 strain), heat shock and standard conditions (3 mutants in each condition, all with significantly impaired growth phenotypes [Fig. 2C]). We conclude that FANOVA analysis recapitulates known roles of TFs but also suggests novel contributions of TFs to cell physiology.

Phenotype network analysis reveals extensive cross-regulation of stress responses by TFs.Many mutants exhibited significant differential growth phenotypes under multiple conditions tested (Fig. 2). At the rank order ||ODΔ|| cutoff (Fig. 2C), 23 of the 27 mutants studied exhibited significant differential growth under at least one condition (Fig. 3). Network analysis of these phenotypes across conditions enabled the classification of the mutants into three categories. (i) Significant differential growth was detected for 12 mutants relative to the Δura3 strain under two or more stress conditions. These TFs were considered to have “cross-stress” functions (Fig. 3). (ii) Significant differential growth was detected for three mutants (ΔtrmB, ΔphoU, and ΔtroR mutants) under both standard conditions and one or more stress conditions. These TFs were considered to have “growth & stress” functions. (iii) Eight TFs were required for normal growth under only one of the conditions tested here and therefore considered “stress-specific” TFs. However, three of the stress-specific mutants also exhibited growth defects under other stress conditions tested ad hoc in previous studies but not included in the systematic phenomic testing here (e.g., ΔsirR, Δidr1, and Δidr2 mutants are required for metal homeostasis [19, 22]), suggesting that some of these “stress-specific” TFs might also be required for protection across multiple stressors.

• Open in new tab
• Download powerpoint

FIG 3

Phenotype network analysis reveals three major classes of TF mutants and extensive cross-regulation of stress responses by TFs. Node and edge attributes are given in the key at the bottom of the figure. OD_Δ numbers in the edge color legend refer to the maximum (fast growth) or minimum (slow growth) value of the mean of the posterior prediction from the FANOVA model across the entire growth time course. ||ODΔ||numbers for edge thickness refer to the median value of the distribution from ||ODΔ|| boxplots (see Fig. 2C, Fig. S3, and Materials and Methods).

Several mutants with strong growth impairments relative to the parent control strain under one condition exhibited growth improvement under others. For example, although the ΔtrmB mutant grew poorly under standard conditions, it showed significantly improved growth under PQ and osmotic stress conditions (maximum OD_Δ of 0.5 and 1.0, respectively; Fig. 2A and 3). Similarly, the ΔrosR mutant exhibited substantially improved growth under low osmolarity (maximum OD_Δ of 0.8) but strong growth defects under oxidative stress induced by PQ and peroxide. Deletion of the gene encoding a third TF, CspD1, led to increased growth relative to the control strain under standard conditions but impaired growth under oxidative stress (Fig. 2A and 3 and Fig. S2). These observations suggest novel functions for these previously characterized TFs. These opposing phenotypic patterns for individual TFs could result from direct regulation of genes required for growth and/or stress resistance under these conditions. For example, ribosome levels are directly related to growth rate across the tree of life (51, 52), including archaea (25), and H. salinarum RosR directly regulates ribosome biosynthesis genes (50). Alternatively, alteration in growth rate per se could change cellular stress resistance properties; for example, slow growth in wild-type yeast cells is associated with heat shock resistance (53). Together, the classes of TF mutants identified here suggest that a surprisingly high number of TFs are required for growth homeostasis during exposure to multiple stressors, suggesting extensive network interconnectivity.

The structure of the gene regulatory network corresponds strongly with TF physiological functions.To gain insight into possible regulatory mechanisms underlying such phenotypes, we determined the relationship between growth trajectories by clustering phenotypes according to the OD_Δ metric (Fig. 4) and asked how correlated phenotypes mapped to GRN topology. In particular, the ΔrosR mutant was previously shown to regulate 20 other TFs, 7 of which were included in the strain collection evaluated here (ΔarcR, ΔcspD2, Δhlx2, Δhrg, Δtrh4, ΔVNG0039H, and ΔVNG0194H strains). Hierarchical clustering of OD_Δ phenotypes across the five conditions revealed that these mutants clustered closely together under oxidative stress, including PQ (Fig. 4A) and peroxide conditions (Fig. 4B). Under PQ stress, correlation of the OD_Δ trajectories for these mutants with that of the ΔrosR mutant were significantly enriched for strong positive correlations (ρ^≥0.4) relative to all other pairwise correlations between the 27 mutants (P ≤ 6.90 × 10⁻³ by the hypergeometric distribution test [Fig. 4C]). Under peroxide stress, according to||ODΔ||, knockouts in genes encoding TFs regulated by RosR were also significantly enriched for impaired growth relative to all other mutants (P ≤ 0.034 by the hypergeometric distribution test). Such significant phenotype correlations were not observed under other conditions tested (Fig. 4C and Fig. S4), suggesting that the interconnection between these TFs is coordinated specifically under oxidative stress (50). Together, these results are consistent with the hypothesis that TFs that regulate each other in GRN subnetworks controlling the cellular response to a particular stress are predictive of TF knockout phenotypes and vice versa.

• Open in new tab
• Download powerpoint

FIG 4

TFs that regulate each other have similar OD_Δ phenotype trajectories. (A and B) Heat maps depict hierarchical clustering of OD_Δ trajectories for the 27 TF knockout mutants under paraquat (PQ) (A) and peroxide (B) conditions. TFs transcriptionally regulated by RosR are known (50) and indicated by red text. Colors in the dendrogram represent different clusters. The color scale indicates the mean of the posterior OD_Δ distribution for each mutant across the growth time course (x axis). (C) Each of the seven mutants regulated by RosR are statistically enriched for strongly correlated phenotypes (ρ^≥ 0.4) with the ΔrosR mutant under PQ and standard growth conditions relative to other conditions (red text). Colored squares represent significant correlations (P ≤ 0.001), and white squares represent nonsignificant correlations. See the color scale to the right of the figure for correlation values.

Previous predictions of TF functions from computational GRN models specific to oxidative stress for H. salinarum (15) were compared with phenotypic results obtained here. Of the 27 knockouts studied, 15 were predicted by network analysis to play a role in gene regulation during oxidative stress (15). Of these 15, 14 TF knockouts showed significantly altered growth relative to the Δura3 control in oxidative stress induced by PQ, peroxide, or both (Fig. 2C, Fig. S3, and Table S1). This suggests that the GRN computationally predicted from gene expression data has strong predictive power for the roles of TFs in cell physiology.

Validation and characterization of novel TF functions. (i) CopR functions as a regulator of both heat shock and copper overload responses.To support the stress-specific novel functions for individual TFs discovered here, we performed validation experiments. First we tested the observation of secondary functions for the previously characterized TF, CopR (Fig. 3). We detected significantly impaired growth of the ΔcopR mutant relative to the Δura3 parent strain during heat shock (Fig. 2A). The heat shock phenotype of the ΔcopR mutant was the top-ranking mutant under this condition (Fig. 2C). These observations were surprising given that CopR (previously called VNG1179C [22, 54]) was previously characterized as a repressor of P1-type ATPases that export copper during overload (22). To conduct further functional validation of the role for CopR in response to elevated temperature, a wild-type copy of the copR gene was expressed in trans on a plasmid in the ΔcopR deletion background, which returned the growth of this strain to levels indistinguishable from that of the Δura3 parent control (Fig. 5A). This complementation result indicates that the ΔcopR heat shock growth defect was caused by loss of the copR gene alone and not by polar or off-target secondary site genetic effects.

• Open in new tab
• Download powerpoint

FIG 5

Phenotype validation: a novel heat shock function for the copper-responsive regulator CopR. (A) Slow growth under heat shock conditions in the ΔcopR mutant (left graph) is complemented by expression of the copR gene in trans (right graph). (B) Cytoscape gene regulatory network depicting the significant overlap between genes regulated by CopR in response to copper overload (CopR node at the top) and in wild-type cells in response to heat shock (heat node at the bottom) (55) and copper (copper node on the left) (22). Node colors (representing expression levels in wild-type cells under heat shock) and edge line types are indicated in the keys. Gray boxes behind groups of nodes represent arCOG functional categories (61). P values of enrichment were calculated by the hypergeometric distribution test. Genes of unknown function are not shown here for clarity but are given in Fig. S5. (C) Quantitative RT-PCR gene expression of cctA and nirJ genes in the knockout strain compared to the Δura3 control strain. The levels of expression before heat shock (0′) and 30 min after induction of heat shock (30′) are shown. Expression is normalized relative to a control gene whose expression does not change during heat shock (see Materials and Methods). Error bars represent standard errors of the means (SEM) of three biological replicate cultures, each with three technical replicate trials.

To detect potential gene targets of CopR regulation under heat shock, existing transcriptomic and TF-DNA binding data sets were reanalyzed. In the wild-type H. salinarum strain exposed to heat shock, 247 genes were differentially expressed (55, 56). A significant fraction of these heat-responsive genes were also differentially expressed in the ΔcopR mutant grown under copper overload conditions (22) and/or bound by CopR under optimum growth conditions (54) (63 of 247; significance by the hypergeometric distribution test, P ≤ 1.42 × 10⁻⁸; Fig. 5B and Fig. S5). Chaperones and amino acid metabolic functional categories were significantly enriched among heat-induced and CopR-repressed genes (Fig. 5B). In contrast, energy generation, translation, and transcription functions were significantly enriched among heat-repressed and CopR-induced genes (Fig. 5B). Of these 63 heat-responsive CopR-regulated genes, 6 overlapped with the list of 10 genes whose transcripts and proteins were most strongly induced in response to heat shock (56) (hypergeometric test P ≤ 5.0 × 10⁻⁵). These six genes encoded proteins required for cellular repair following heat shock (chaperones CctA and Hsp5) and protection from further damage (metalloprotein NirJ, ferritin DpsA, and anaerobic metabolic genes ArcAC). Together, these CopR-regulated gene functions are consistent with a cellular need to arrest growth to refold and regenerate degraded proteins under elevated temperatures (57).

CopR regulation of gene expression described above was tested under standard and copper overload conditions (22, 54). To validate whether CopR is also specifically required for regulation during elevated temperature, gene expression was measured by quantitative reverse transcription-PCR (qRT-PCR) in the ΔcopR strain immediately before and 30 min after a shift from 42°C to 54°C. cctA expression in the ΔcopR mutant cells was 3- to 4-fold higher under standard conditions than in Δura3 parent control cells and remained elevated upon heat shock (Fig. 5C). This indicates that CopR is required for cctA repression under standard conditions and that relief of CopR repression is necessary for cctA heat induction. In contrast, 7-fold induction of nirJ in the Δura3 strain was abrogated in the ΔcopR mutant, indicating that CopR is an activator of this gene under heat shock (Fig. 5C). Taken together, these validation analyses suggest that CopR functions both as a global regulator of gene expression under heat shock and as a specific regulator of a copper efflux transporter during copper overload (22, 58). The connection between heat shock and copper overload has not been established in archaea; however, because both heat shock and copper overload can induce the accumulation of oxidative radicals, the transcriptional responses to these common types of cellular damage may be linked (4).

(ii) A novel oxidative stress response function for CspD1, a conserved cold shock family protein.The ΔcspD1 strain exhibited the strongest growth defect of all 27 mutants under oxidative stress induced by PQ and the sixth strongest phenotype under peroxide (Fig. 2C and 3). The CspD1 sequence showed strong similarity to the cold shock family of proteins (E value of 2.30 × 10⁻²¹; Table S1). Cold shock domain (CSD) proteins are broadly conserved from bacteria to humans but serve diverse functions, including RNA protection, inhibition of DNA replication during oxidative stress, and transcriptional regulation under various stress conditions (59, 60). Thus, we further investigated the role of CspD1 in response to multiple stresses. We found that, although the cspD1 gene was significantly repressed in response to high temperature (Fig. 5B) (56), the growth of the ΔcspD1 mutant strain was indistinguishable from that of the Δura3 control strain under heat and cold shock conditions (Fig. S6). In addition, cspD1 expression was induced under copper overload conditions in a CopR-dependent manner (Fig. 5B) (22); however, the ΔcspD1 mutant did not exhibit impaired growth under copper overload conditions (22). These data suggest either that CspD1 does not play a role in the temperature shock and copper overload responses or that CspD1 is functionally redundant under these conditions with other, as yet unknown, TFs.

To further test the role of CspD1 in oxidative stress, we expressed the cspD1 gene in trans on a plasmid in the ΔcspD1 knockout background. The oxidative stress phenotype was complemented in this strain, confirming that deletion of cspD1 was responsible for the observed phenotypes (Fig. 6A). In previously published transcriptomic data, the expression of the cspD1 gene was strongly correlated with fluctuations in oxygen levels (1), induced 30 min after an increase in oxygen (cross-correlation ρ^ = 0.730; Fig. 6B). Such expression follows a pattern similar to that observed for rosR, which encodes a known global regulator of the oxidative stress response (16, 50). In the ΔcspD1 strain exposed to fluctuating oxygen levels over time (see Materials and Methods), the transcription of 132 of the 660 known oxygen-responsive genes in H. salinarum (1) exhibited significantly altered expression (Fig. 6C and Table S4). Of these 132 genes, 89 are induced in a CspD1-dependent manner during the transition from anaerobic to aerobic conditions. According to arCOG (clusters of orthologous genes for archaea) ontology (61), these aerobic genes are significantly enriched for functions crucial for cell growth (e.g., translation; P ≤ 5.68 × 10⁻²⁰ by the hypergeometric distribution test; Table 2). In contrast, CspD1 is required to repress 43 genes under anaerobic conditions (Fig. 6C and Table 2). Of the 132 genes requiring CspD1 for appropriate expression under oxygen fluctuation, 106 are also differentially expressed under oxidative stress induced by PQ (Fig. 6D) (15).

• Open in new tab
• Download powerpoint

FIG 6

Phenotype validation: a novel oxidative stress function for the cold shock family protein CspD1. (A) Slow growth under oxidative stress conditions in the ΔcspD1 mutant (left graph) is complemented by expression of the cspD1 gene in trans (right graph). (B) Line plot depicting the expression of the cspD1 gene (left axis) during fluctuations in oxygen concentrations (gray line, right axis). cspD1 expression is compared to that of the gene encoding known oxidative stress regulator RosR (left axis). (C) Line plot depicting the expression of genes requiring CspD1 for appropriate dynamic expression in response to oxygen. Each line represents the mean expression value of the groups of genes indicated in the legend. Expression data and annotations for individual genes are given in Table S4. WT, wild type. (D) Expression profiles for 106 of the 132 CspD1-dependent genes differentially expressed under PQ conditions in wild-type cells (15). Thin lines represent the expression of individual genes. Thick lines show the mean of induced or repressed genes. (E) Overlap between EGRIN predictions for CspD1 target gene regulatory influences (7) and differentially expressed genes (DEG) in the ΔcspD1 background. (F) Similar to panel E except that predictions under oxidative stress conditions are shown (15).

The EGRIN model based on gene expression under a wide array of conditions accurately predicted a significant fraction of the 132 CspD1-dependent genes (42%; P ≤ 6.33 × 10⁻³² by the hypergeometric distribution test) (7) (Fig. 6E). These EGRIN-predicted genes were also significantly enriched for functions involved in ribosome biogenesis (P ≤ 1.8 × 10⁻²⁶). Similarly, EGRIN predictions based solely on gene expression under oxidative stress also predicted that CspD1 regulates a significant fraction of the 132 CspD1-dependent genes (P ≤ 3.83 × 10⁻¹⁵) (15) (Fig. 6F). Together, these results provide validation of GRN functional predictions and phenotype analysis, implicating CspD1 in the regulation of functions critical to growth under oxidative stress and fluctuating oxygen conditions.

Conclusions and perspectives.Data and analyses presented here enabled the discovery of physiological roles for 17 previously uncharacterized TFs in the archaeal species H. salinarum. New physiological roles for previously characterized TFs were also revealed (e.g., CspD1 and CopR). This demonstrates the power of our combined high-throughput growth analysis and quantitative growth modeling for discovering unknown gene functions. The functions of a large fraction of genes and pathways remain unknown even in well-characterized model microorganisms. Recently, other high-throughput, genome-wide forward and reverse genetic approaches have made great strides in gene functional discovery, such as population genomics in Saccharomyces cerevisiae (62), clustered regularly interspaced short palindromic repeat interference (CRISPRi) in Bacillus subtilis (63), and genome-wide knockout collections in Escherichia coli (33, 64), as well as nontraditional model organisms such as methanogenic archaea (65). Previously we showed that Gaussian process (GP) methods for phenotype discovery are generally applicable across diverse species (35). Additionally, GPs have been shown to be useful in other domains of biological research, such as modeling genome-wide expression data (66–69). Here we extend the use of GP methods in the context of functional ANOVA to compare the growth of a large collection of strains (Fig. 2), demonstrating the promise of these methods for genome-wide functional discovery across a variety of species.

Here we directly test the predictions of computationally inferred, global GRN models such as EGRIN (7, 15). Because these and other statistical GRN inference models were inferred directly from gene expression data, the direct impact of the GRN activity on cellular physiology remained unclear. Here we show that models such as EGRIN predict not only gene expression but also the phenotypic impact of such expression (Fig. 2 and 6). In previous work, we also demonstrated the accuracy of EGRIN in predicting TF functions. For example, hypotheses generated from EGRIN predictions enabled the discovery of RosR, an archaeon-specific, novel master regulator of oxidative stress response (7, 15, 16, 50). Similarly, EGRIN model predictions regarding the cross-regulation of phosphate metabolism and methanogenesis pathways were validated by gene knockout studies in methanogens (30). Here we also observed extensive cross-regulation: each TF was important for resistance of multiple stressors, and multiple TFs played a role in surviving each stressor (Fig. 3). Such cross-regulation of gene expression has also been observed in bacteria as a means to integrate environmental cues that cooccur (e.g., heat and singlet oxygen [8, 70]) or that induce functionally related response pathways (e.g., metal homeostasis and oxidative stress [71]). Together with these previous studies, the broader investigation of 27 TF knockout phenotypes reported here demonstrates that GRN models such as EGRIN are effective tools for generating accurate hypotheses regarding TF functions. The combination of GRN modeling and phenomic validation reveals the direct impact of a complex web of regulatory interactions on cell physiology (Fig. 4 and 5). This work supports an emerging general principle that cross-regulation between TFs within the GRN enables a coordinated response to a variety of environmental stimuli.