The information for six bacterial genomes used in this work was downloaded from either the NCBI server or RegulonDB: E. coli K-12 MG1655 (NC_000913.3, GCF_000005845.2), B. subtilis 168 (GCF_000009045.1), P. aeruginosa PAO1 (GCF_000006765.1), S. typhimurium LT2 (GCF_000006945.2), S. aureus N315 (GCF_000009645.1), and M. tuberculosis (GCF_000195955.2). For each genome, the FASTA sequence was obtained from the “gbff/gbk” files parsed with an ad hoc program (Supplementary material, ParserGBK.py), to add the appropriate label in the header: NCBI gene ID, local gene ID, gene name, product description, and organism name. Sequences with missing information were annotated as “NODATA.” In addition, the 12,230 genomes of bacteria and 649 archaeal genomes were downloaded from the NCBI RefSeq genome database on May 18, 2021, to infer their GRNs.

The regulatory interactions were obtained from specialized databases [DBTBS for B. subtilis release 5 (Sierro et al., 2008),¹ RegulonDB release 10.9 for E. coli (Santos-Zavaleta et al., 2019a),² M. tuberculosis (Kapopoulou et al., 2011; Sanz et al., 2011), RegulomePA release 1.0 for P. aeruginosa,³ Salmonet release 2.0 for S. typhimurium LT2 (Métris et al., 2017), and for S. aureus N315 (Ravcheev et al., 2011; Poudel et al., 2020)] and posteriorly homogenized, following the same format: First column corresponds to the assigned number by regulatory interaction per organism; second column, TF associated; third column, Target gene; and the other columns indicate the annotations derived from the original networks (Supplementary material). These GRNs are summarized in Table 1.

The protein sequences from each model organism were used as reference to identify the orthologs in an all-vs-all genomes fashion using the program Proteinortho (Lechner et al., 2011) with the following parameters: E-value ≤10⁵, coverage ≥70%, and identity of ≥25%, as previously described for the identification of TFs (Flores-Bautista et al., 2020).

The predictions of Transcription Units (TUs) or operons were obtained using the method described by Moreno-Hagelsieb and Collado-Vides (2002). In brief, the predictions were based on the transcription direction and the intergenic distance (shorter intergenic distances and in the same direction for genes in the same TU).

The reference genomes were used to scan the 12,230 bacterial and 649 archaeal genomes to identify their orthologs and map their interactions considering the following criteria: If the orthologs of the TF and its TG of the model organism were found in a new genome, the interaction was assigned using guilt by association. In a second step, predicted TUs were used to expand the TF-TG interactions as follows: If the first gene of the orthologous TG in an organism corresponded to the first gene in the TU, the other genes belonging to the TU were associated with the same TF. Finally, each network was integrated using all the ortholog assignments with the six reference GRNs. All the network interactions can be inferred by running the script pipeline.sh, provided as Supplementary material and Figure 1.

The GRNs were analyzed by using Cytoscape (Shannon et al., 2003; Otasek et al., 2019) to obtain their degree, clustering coefficient, and other centrality metrics. Hubs were obtained by using networkX from python (Hagberg et al., 2008). In addition, to identify transcriptional co-regulators and modules in a GRN, the CoReg software was used. In brief, CoReg calculates gene similarities based on the number of common neighbors of any two genes in the network (Song et al., 2017).

The GRNs inferred for all the bacterial and archaeal genomes are available through the web server at https://regulatorynetworks.unam.mx/, which is built on HTML5, JQuery, and Php languages, while the data are stored in a MySQL database. For the data display, we use the Cytoscape JS (Franz et al., 2016) framework due its capabilities to represent nodes and edges of the network with determined properties, allowing users to change forms, colors, and layer visualization of the network.

GRNs were compared using two different approaches to establish the reliability of the approach, one based on the ability to recover edges and the second focusing on the ability to recover network motifs (Milo et al., 2002) by comparing the six reference networks with networks for the same genomes generated using our approach.

First, based on the orthologous annotations made with Proteinortho, we created a GRN for each of the reference bacteria using a naming convention that ensures that genes that are orthologous among them share the same name in each of the six GRNs used as reference. Then, we created networks for each reference organism based on the regulations transferred from the other five GRNs, again using the consensus gene names. GRNs with consensus gene names were then compared, following two procedures implemented in LoTo (Martin et al., 2017). We employed binary classification metrics to evaluate the similarities between pairs of GRNs as follows: Edges present in both compared networks are considered true positives (TPs), genes only present in one of the networks are false negatives (FNs) if they are only in the reference network, True Negatives (TNs) are the edges absent in both compared networks, and false positives (FPs) if they appear only in the network we compared with the reference. This edge-based approach is used to compare predicted GRNs versus reference networks, and it indicates overall network similarity (The DREAM5 Consortium et al., 2012). The second approach relies on the presence or absence of the motifs defined by Milo et al. (2002) that have been related to functional patterns in GRNs. Instead of considering TF-gene interactions, in this second approach, we considered TP motifs present in both compared networks, FN motifs are only found in the reference network, and FPs are only present in the network compared against the reference GRN.

LoTo calculates several metrics, but here we only focused on the most employed ones:

and

To establish a baseline and determine whether the results from our approach are significant versus what can be expected by chance, we also created a protocol to determine the expectancy of a transferred TF-gene regulation by chance. We randomized the names TFs for the whole inferred networks 10,000 times to calculate expected TP, FP, TN, and FN values by comparing these randomized networks against their reference counterparts. This protocol ensures comparisons of random networks with the same characteristics, e.g., edges, TFs, and genes, against their actual reference. We then employed a G-test as implemented in SciPy (Virtanen et al., 2020) to determine whether the observed number of edges considered TP, FP, TN, and FN can be from the same distribution as that observed for the predicted networks without randomization.

In order to evaluate and expand the GRNs of the six model organisms, the number of TFs, TGs, and their interactions was determined. To do this, we downloaded six GRNs, and their interactions were displayed by using Cytoscape. In this work, we considered TFs as those proteins that activate or repress gene expression but do not belong to the transcriptional basal machinery; therefore, sigma factors, antiterminators, terminators, and sensor proteins, among other proteins, were excluded from the resulting data set (Martínez-Núñez et al., 2013). Table 2 shows the number of interactions associated with each organism. The most studied bacterial species, E. coli K-12, has 3,616 interactions based on experimental evidence, followed by S. typhimurium LT2 (2,969 interactions) and B. subtilis with 2738TF-TG interactions, whereas the GRN of S. aureus contains the smallest number of interactions, with 709. This difference could be a consequence of the experimental evidence accumulated over the years and the number of experiments carried out and performed with each organism; i.e., there is a bias inherent to the experimental analysis towards specific organisms. For instance, in a recent collection of 668 experimentally characterized TFs in bacteria and archaea organisms (Flores-Bautista et al., 2020), 33.5% was associated with E. coli K12, 23% with different strains of M. tuberculosis, and 19% with B. subtilis 168; i.e., 76% of the complete collection is concentrated in few organisms; in contrast, 24% of the collection is distributed among 78 different prokaryotes. This contrast in the information is also evident in more general databases, such as UniProtKB/Swiss-Prot, where E. coli K-12 is the bacterial organism with more proteins deposited and curated manually in the database.⁴

To determine the number of interactions shared between the six model organisms, we first used the program Proteinortho to assign orthology relationships between all proteins in the proteome of each bacterium. Once orthologous proteins were determined, we inferred regulatory interactions between organisms based on the presence of an orthologous TF and an orthologous target of that TF in the model GRN. In the second step, the interactions were expanded by using the TU assignments, as described in Materials and Methods. This comparison showed that E. coli and S. typhimurium LT2 share a high number of interactions, because of their phylogenetic closeness. In contrast, the actinobacterium M. tuberculosis is the organism with the lowest number of shared interactions with the other bacterial models as a consequence of its phylogenetic distance; only 12% (in average) of its interactions are shared with other bacteria (see Supplementary material).

In order to infer new interactions among the six bacterial genomes, they were compared and their interactions were assigned based on the presence of the TF-TG orthologous pairs. In this regard, Table 1 shows the number of new assignments and their proportion per organism. From this analysis, we found between 776 and 2,499 new interactions, with S. typhimurium LT2 and P. aeruginosa the organisms determined to have more new interactions inferred. These larger numbers for S. typhimurium LT2 and P. aeruginosa are probably a consequence of their phylogenetic closeness with E. coli K-12 (Fukushima et al., 2002) in comparison to the other organisms used as models. It is important that some regulatory interactions were found in more than one organism; therefore, the sum of the rows may not correspond to the total number of new interactions, as is the case for the regulator PhoB (NP_414933.1) of E. coli K-12, which regulates the cytochrome bd-I ubiquinol oxidase subunit (NP_415262.1), as inferred from the interactions previously described in the B. subtilis and M. tuberculosis networks.

Regarding the reliability of interactions predicted by our approach, we compared networks with only TF-TG interactions derived from homology relationships for each of the six species with the respective reference GRNs. The comparisons were made by considering this to be a binary classification problem, and thus, edges (and graphlets) in both the reference network and the predicted GRN are TPs, edges only in the reference are FNs, and edges only in the predicted network are FPs. These results, shown in Table 2 for single edges and in Table 3 for graphlets, indicate a varying range of values depending on the compared bacteria. For recall (R), the rate of recovered TF-TG interactions ranged from 0.028 for M. tuberculosis to 0.52 for S. enterica, whereas precision (P), which indicates the likelihood that the existence of an edge is correctly predicted, ranged from 0.20 for P. aeruginosa to 0.69 for S. enterica. These results are significantly different from those expected by chance, as shown by the very low p-values obtained with the G-test. When the same metrics for the presence and absence of graphlets were used (Table 3), we found a similar trend for each model GRN but with lower values for each metric. Lower values for the metrics calculated with graphlets are expected, since a single edge that differs between two networks often affects various graphlets.

Based on the expanded networks, we identified new TF-TG interactions described in Table 4 that must be exhaustively analyzed. In this regard, we found an increase in the number of targets, TFs, nodes, and interactions for all the bacterial and archaeal extended networks. For instance, for M. tuberculosis H37Rv, there was an increase of 776 new interactions (305 new TGs and 31 new TFs), whereas for B. subtilis, 1821 new interactions (36 new TFs and 553 new TGs) were identified. Therefore, we performed a literature search to find evidence to support our predictions. Based on these searches, and considering the 1,676 new interactions for E. coli K-12 (56 new TFs and 570 new TGs), we identified that 179 of these interactions have been described in the literature (Supplementary material); however, they are not deposited in RegulonDB. In particular, we found that the interaction of SoxS and ompW in the GRN of E. coli and inferred from S. enterica has been experimentally described. In E. coli, ompW is regulated by three TFs, as described in RegulonDB; however, we found that it could be also regulated by SoxS (Zhang et al., 2020) in a negative fashion.

We also found a new interaction, where CpxR could be regulating tar gene. This TF together CsgD has been described in bacterial adhesion, and belongs to the stationary-phase response (Santos-Zavaleta et al., 2019b). Experimental evidence suggests that CpxR and CsgD repress the transcription of fliA, flgM, and tar (Dudin et al., 2014), in addition to bglg and bglb (Mattéotti et al., 2011), and PdhR and lipA (Kaleta et al., 2010). These regulatory interactions identified by our orthologs inferences have not been documented in RegulonDB.

The GRNs inferred in all the bacterial and archaeal genomes are available through a web server whose interface is shown in Figure 2. The GRN of user-selected organisms are shown in an embedded interactive display that through a very intuitive mouse-based interface allows the user to select subnetworks and different types of regulatory interactions. Edge and node colors can also be redefined, as well as the layout used in the network visualization, depending on their properties. Additionally, the user can display and visualize information related to Genes (name, protein ID, initial and end coordinates, and strand), and edges among nodes representing genes, including information about whether this is a new or known edge and the organism from which it was derived. Additionally, if information is available, by clicking on the node name or protein identifier, you can access the NCBI/Uniprot page related to the gene of interest.

Entire GRNs or used defined subnetworks can be downloaded in standard format for further inspection with tools such as Cytoscape, that in addition, connect our tool to the whole array of apps already available for this visualization tool. For more information and a more detailed description of both the input and output files, see the website https://regulatorynetworks.unam.mx/ or http://132.247.46.6/, help section, where an example is provided.