An updated gene regulatory network reconstruction of multidrug-resistant Pseudomonas aeruginosa CCBH4851.

BACKGROUND: Healthcare-associated infections due to multidrug-resistant (MDR) bacteria such as Pseudomonas aeruginosa are significant public health issues worldwide. A system biology approach can help understand bacterial behaviour and provide novel ways to identify potential therapeutic targets and develop new drugs. Gene regulatory networks (GRN) are examples of in silico representation of interaction between regulatory genes and their targets.
OBJECTIVES: In this work, we update the MDR P. aeruginosa CCBH4851 GRN reconstruction and analyse and discuss its structural properties.
METHODS: We based this study on the gene orthology inference methodology using the reciprocal best hit method. The P. aeruginosa CCBH4851 genome and GRN, published in 2019, and the P. aeruginosa PAO1 GRN, published in 2020, were used for this update reconstruction process.
FINDINGS: Our result is a GRN with a greater number of regulatory genes, target genes, and interactions compared to the previous networks, and its structural properties are consistent with the complexity of biological networks and the biological features of P. aeruginosa. MAIN
CONCLUSIONS: Here, we present the largest and most complete version of P. aeruginosa GRN published to this date, to the best of our knowledge.

Pseudomonas aeruginosa is a ubiquitous and opportunistic pathogen of which infections can affect the lower respiratory tract, skin, urinary tract, eyes, soft tissues, surgical wound, and gastrointestinal system, among others, leading to bacteraemia, endocarditis, and other complications, particularly in health care settings and in immunocompromised patients. , , This Gram-negative bacteria is one of the most difficult to treat, due to its intrinsic resistance, acquisition of resistance through chromosomal gene mutations, and horizontally acquired resistance mechanisms to multiple drugs. Multidrug resistance (MDR) imposes the central difficulty in the selection of appropriate antibiotic treatment and reduces treatment options, especially in nosocomial settings such as healthcare-associated infections (HAI). ,

HAI is a severe public health issue related to high morbidity and mortality rates in hospitalised patients and high healthcare costs. Worldwide, P. aeruginosa is one of the most prevalent agents of HAI.

In Brazil, the Brazilian Health Surveillance Agency ranked P. aeruginosa as the third most common causative agent of HAI in hospitalised patients in adult intensive care units (ICU) and the second in paediatric ICU, being nearly 40% of the reported strains resistant to carbapenems. This class of beta-lactam antibiotics has been widely administered worldwide for treating P. aeruginosa infections and other MDR Gram-negative bacterial infections. Indeed, a significantly higher mortality rate was observed among patients infected with MDR P. aeruginosa clones (44.6%) compared to those infected with non-MDR (24.8%).

The most epidemiologically important mechanism of carbapenem resistance is the production of carbapenemases. Among MDR P. aeruginosa clinical isolates in Brazil, the most prevalent carbapenemase is the São Paulo metallo-β-lactamase (SPM-1). This enzyme is encoded by the gene bla , located on the P. aeruginosa chromosome, and it confers resistance to almost all classes of beta-lactams. The first register of an MDR P. aeruginosa strain carrying the bla gene found in Brazil is from 2003. Widely disseminated in distinct Brazilian geographic regions, SPM-1-producing P. aeruginosa is associated with the clone SP/ST277 and has been isolated from hospital sewage systems, rivers, and microbiota of migratory birds. , The strain P. aeruginosa CCBH4851 belongs to clone SP/ST277, and was involved in an endemic outbreak in Brazil in 2008. This strain is resistant to most antimicrobials of clinical importance, such as aztreonam, amikacin, gentamicin, ceftazidime, cefepime, ciprofloxacin, imipenem, meropenem, and piperacillin-tazobactam, being susceptible only to polymyxin B, and has several mechanisms of mobile genetic elements. ,

To better understand P. aeruginosa’s behaviour, more comprehensive knowledge of gene expression patterns predicted by analysing its gene regulatory network (GRN) is of great value. A GRN consists of a set of transcription factors (TF) that interact selectively and nonlinearly with each other and with other molecules in the cell to regulate mRNA and protein expression levels.

Mathematical modelling and computational simulations are approaches for analysing the GRN and other complex cellular systems influenced by numerous factors. These models allow the construction of biological networks, predict its behaviour under unusual conditions, identify how a disease might develop, and intervene in such development to prohibit cells from reaching undesirable states. In addition, due to their lower cost and high accuracy, such approaches contribute to developing new drugs.

The P. aeruginosa PAO1 strain had its genome sequence published in 2000, providing information regarding genome size, genetic complexity, and ecological versatility. It has been extensively studied since then.

Published in 2011 by Galán-Vasquez et al., the first P. aeruginosa GRN was based on the PAO1 strain (PAO1-2011). Then, in 2019, Medeiros et al. described a GRN reconstruction of CCBH4851 strain (CCBH-2019). Finally, in 2020, Galán-Vasquez et al. published the updated GRN of P. aeruginosa with the PAO1 strain (PAO1-2020), which was much larger than the previous ones, containing new interactions. All works analysed the GRNs main structural properties and regulatory interactions.

This manuscript describes CCBH-2022, an updated GRN of the MDR P. aeruginosa based on the CCBH4851 strain, using as references both CCBH-2019 and PAO1-2020. We characterise regulators, target genes (TGs), transcription factors (TFs), auto-activation interactions, and influential genes of the network.

We analyse the main structural properties of the network, such as degree distribution, clustering coefficient, and relative abundance of network motifs. Finally, we compare the results of our analyses with those from previous GRNs.

MATERIALS AND METHODS

In this work we study the P. aeruginosa CCBH4851 strain, which is deposited at the Culture Collection of Hospital-Acquired Bacteria (CCBH) located at the Laboratório de Pesquisa em Infecção Hospitalar, Instituto Oswaldo Cruz/Fundação Oswaldo Cruz (Fiocruz) (WDCM947; 39 CGEN022/2010). The genome sequence is available in the GenBank database (Accession CP021380.2).

CCBH-2019 and PAO1-2020 models were the bases for the reconstruction of this GRN. CCBH-2022 GRN results from the orthology analysis between the P. aeruginosa PAO1 and CCBH4851 gene sequences. CCBH-2022 model also inherits the orthologs between CCBH4851 and P. aeruginosa PA7 and P. aeruginosa PA14 strains, which were already present in CCBH-2019 GRN. The evolutionary histories of genes and species reconstruction are based critically on the accurate identification of orthologs. Orthology refers to a specific relationship between homologous characters that arose by speciation at their most recent point of origin, , a common ancestor. One of the most common approaches to determining orthology in comparative genomics is the Reciprocal Best Hits (RBH), which relies on BLAST. An RBH occurs when two genes from different genomes find themselves the best scoring match in the opposite genome. , Regulatory interactions between TFs and TGs in the PAO1 GRN were propagated to CCBH-2022 GRN if the TF and the TG formed an RBH. Medeiros et al. designed and implemented an algorithm using the Python programming language to automate and generate a list of RBHs in a tabular format (available as Supplementary data). All the protein sequences from P. aeruginosa CCBH4851 (P1) and P. aeruginosa PAO1 (P2) were considered. BLAST+ was used to query the proteins from P1 against those from P2 (forward results) and P2 against P1 (reverse results). Each P1 query sequence was considered in turns, and its best match from P2 was identified from forwarding results (x). Likewise, each P2’s query sequence was considered from the reverse results, with its best match in P1 (x’). If x = x’, then they are RBH. Local BLASTP searches of each protein set against the other were executed, with the following cut-off parameters: identity ≥ 90%, coverage ≥ 90%, and E-value ≤ 1 e-5, showing the results in tabular format. If the search returned no hits, the gene was considered to have no ortholog within the opposite genome. Manual BLASTP was used to prevent false negatives, aligning these gene sequences with the opposite genome, considering the above parameters. If they still returned no hits but were present in either PAO1-2020 or CCBH-2019 models, the results were evaluated with a literature search to determine if they were accurate and whether they should be part of CCBH-2022.

The final GRN table is available as Supplementary data and is organised into six columns: Regulatory gene, Ortholog of the regulatory gene, Target gene, Ortholog of the target gene, Mode of regulation, and Reference. The first column lists the regulatory genes of P. aeruginosa CCBH4851, while the second column contains orthologs of regulatory genes in the reference strain (PAO1 and PA7 or PA14 from the exclusive interactions in CCBH-2019; the same applies to TG’s orthologs). The third column refers to the target gene in P. aeruginosa CCBH4851, while the fourth column lists orthologs of TGs in the reference strain. Finally, the fifth column describes the mode of regulation, and the sixth column indicates the corresponding data source.

The interactions between transcription factor proteins and the genes they regulate in an organism define a directed graph. For the computational analysis, the structure of GRN can be represented as a directed graph, formed by a set of vertices (or nodes) connected by a set of directed edges (or links). Basic network measurements are related to vertex connectivity, the occurrence of cycles, and the distances between pairs of nodes, among other possibilities.

The degree of vertices is the most elementary characterisation of node i, and k(i) is defined as its number of edges. In directed networks, there are incoming (k-in degree) edges and outgoing (k-out degree) edges. The degree distribution can follow a functional form P(k) = Ak -γ, called power-law distribution, where P(k) is the likelihood that a randomly chosen node from the network has k direct interactions, A is a constant that ensures that the P(k) values add up to one, and γ is the degree exponent. , , , , According to Albert, this function indicates high diversity of node degrees, with the P(k) value decaying as a power law that is free of a characteristic scale, resulting in the absence of a typical node in the network that could be used to characterise the rest of the nodes. Most real networks with structural information available exhibit this scale-free behaviour, deviating from a Poisson distribution expected in a classical random network. ,

Studies have shown a scale-free structure in cellular metabolic networks, , protein interaction networks, including in cancer, , transcription regulatory networks, and GRN. , , , Following the literature, , , , , there are some qualitative and quantitative characteristics to ensure that a network is scale-free: the power-law distribution appears as a straight line on a log-log plot; the γ value usually is in the range 2<γ<3; and the presence of high-degree nodes, called hubs, the most highly connected nodes, with most nodes clustered around them. The hubs demonstrate the absence of a uniform connectivity distribution in the network, presenting the 80-20 rule (also referred to as the Pareto principle), with small-degree nodes being the most abundant. However, the frequency of high-degree nodes decreases slowly. Hubs are fundamental for determining therapeutic targets against an infectious agent. Scale-free networks are heterogeneous, so random node disruptions generally do not lead to a significant loss of connectivity. However, the loss of the hubs causes the breakdown of the network into isolated clusters. Some studies validate these general conclusions for cellular networks. , ,

In the GRN, determining the vertices with the highest k-out degrees is a method for identifying a hub, The degree threshold is the exact number of interactions that characterise a hub, and this criterion differs from one study to another. The degree threshold adopted in this work was the average number of connections of all nodes having at least two edges, resulting in a cut-off value of 16 connections.

Motifs are connectivity patterns, a small set of recurring regulation patterns from which the networks are built , that are associated with specific functions. A triangle, i.e., three fully connected vertices, is the simplest type of motif. These genes are a regulator, X, which regulates Y, and gene Z, which is regulated by both X and Y. Triangles can be closed (three connections within the set) or open (two edges). This 3-genes motif is the feedforward loop (FFL) and the most common in GRN, appearing in gene systems in bacteria and other organisms, , with the possibility of either activation or repression in each of the three regulatory interactions. The coherent type-1 FFL and the incoherent type-2 FFL occur more frequently in transcriptional networks.

The clustering coefficient is the probability that two genes with a common neighbour in a graph are also interconnected. This measure has two popular definitions: the local and global clustering coefficient. The local clustering coefficient of vertex i, Ci, is defined as Ci = 2ei / ki (ki - 1), where ei is the number of edges connecting node i with degree ki, and ki (ki - 1) / 2 is the maximum number of edges in the neighbourhood of node i. In GRNs, the local clustering coefficient C(i) is interpreted as the interaction between genes forming the regulatory groups. The clustering coefficient of a network, C, is calculated by the average of Ci over all vertices. , Not considering the directionality of the edges, the global clustering coefficient is the ratio between the number of closed triangles and the total number of triangles (open or closed) in the network. C(k) represents the mean clustering coefficient over the vertices with degree k. Some biological networks tend to present high clustering coefficient values, e.g., in the protein-protein interaction network of S. cerevisiae, ≈ 0.18. ,

The network density measure is the number of edges of the network over the maximum possible number of edges, measuring the interconnectivity between vertices, and is strongly correlated to the potential to generate gene expression heterogeneity. The network diameter is the path length between the two most distant nodes. The average path length is the measure that indicates the distances between pairs of vertices (the average of the shortest path length over all pairs of nodes in the network).

Several genes are connected in the GRN. When the nodes interact through a direct or an indirect link (intermediate connections), they are considered part of a connected component. These associations are the concept of network connectivity, and for this analysis in the present work, network interactions were considered undirected.

Analysing the structural characteristics (connected components, hubs, and motifs) can help determine the best approach to disturb a network to promote a desired phenotype in the cell.

For CCBH-2022 structural analyses, the R programming language and RStudio were used. Scales, dplyr, tibble, readr and igraph packages were used for data manipulation and plotting the structural analyses. , , The igraph library was used to compute most properties described above: the in and out degrees, centrality, clustering coefficients, feed-forward loop motifs, connectivity, cycles, paths, and hierarchical levels analyses.

The illustrations of the GRN, the hub’s network, and the connectivity analysis were made in Cytoscape. All figures are presented with higher resolution in the Supplementary data.

The codes for the structural analysis in R and for finding RBH in python, implemented by Medeiros et al., and the CCBH-2022 file in CSV format are available as Supplementary data in our Github repository (https://github.com/FioSysBio/CCBH2022).

RESULTS

CCBH-2022 consists of 5452 regulatory interactions among 3186 gene products, of which 218 were identified as regulatory genes and 2968 as target genes. Of the 218 regulatory proteins, 87 are TFs, 19 are sigma factors (SF), and 13 are RNAs. Of these 13 RNAs, 11 are SF as well. The tables containing their relations are presented in the Supplementary data.

Given the 6577 predicted protein-coding genes of P. aeruginosa CCBH4851, the model organism in this study, the current network represents roughly 50% of the genome, against 16.52% from CCBH-2019.

Specific regulatory genes and their interactions were kept as described in CCBH-2019, such as the ones resulting from the P. aeruginosa PA7 and P. aeruginosa PA14 orthology, and in dedicated biological databases and scientific literature, e.g., IHF (integration host factor). This bacterial DNA-bending protein, essential in gene expression regulation, is absent in the CCBH4851 genome. However, Delic-Attree et al. demonstrated that P. aeruginosa contains the IHF protein composed of the products of the himA and himD genes. These genes act in combination as a TF for several TGs, and all were listed as regulatory genes in CCBH-2019. Consequently, equivalent annotations to the previous CCBH4851 GRN were maintained.

CCBH-2022 has 5452 edges, and these interactions were classified into activation (“+”), repression (“-”), dual (“d”, when, depending on the conditions, the regulatory gene act as an activator or a repressor), and unknown (“?”), as described in biological databases and scientific literature. An illustration of CCBH-2022 is presented in Fig. 1.

Fig.1:

visualisation of CCBH-2022. Yellow circles indicate regulatory genes, light blue circles indicate target genes (TGs), black lines indicate an unknown mode of regulation, green lines indicate activation, and red lines indicate repression. Purple lines indicate a dual-mode of regulation. A: the gene regulatory networks (GRNs) large highly connected network component; B: all regulatory and TGs with no connections with A.

Regarding the structural measurements of the updated network, the summarised statistical results are presented in Table I. It contains the standard measures (the number of nodes and edges, number of autoregulatory motifs, network diameter, and average path length), the number of feed-forward motifs, and clustering coefficients. Also, Table I presents a comparison with data from PAO1-2011, CCBH-2019, PAO1-2020, and CCBH-2022.

TABLE I

Comparison of structural statistic measures between PAO1-2011, CCBH-2019, PAO1-2020, CCBH-2022

	PAO1-2011	CCBH-2019	PAO1-2020	CCBH-2022
Vertices	690	1046	3009	3186
Edges	1020	1576	5040	5452
Regulatory genes	76	138	173	218
Target genes	593	908	2709	2968
Positive regulation	779	772	3851	3829
Negative regulation	218	454	390	649
Dual regulation	11	13	10	19
Unknown regulation	12	337	789	955
Autoregulation (total)	29	72	50	91
Positive autoregulation	16	21	24	29
Negative autoregulation	13	39	15	46
Unknown autoregulation	-	12	11	17
Feed-forward loop motifs (total) a	137	208	702	968
Coherent type I feed-forward loop motifs a	82	79	226	239
Incoherent type II feed-forward loop motifs a	3	4	8	10
Density	2.12e-03	1.44e-03	6.07e-04	5.99e-04
Diameter	9	12	12	12
Average shortest path length	04.08	4.80	04.01	4.67
Global clustering coefficient	2.28e-02	3.2e-02	3.03e-03	4.42e-03
Local clustering coefficient	2.5e-01	1.92e-01	1.63e-01	1.87e-01

a: number of feed-forward loop motifs determined using the igraph package.

Since both CCBH-2022 and PAO1-2020 contain significant updates from their previous counterparts, the comparison between CCBH-2022 and PAO1-2020 is most relevant. CCBH-2022 had a density of 5.99e-04, slightly lower than the density of PAO1-2020 (6.07e-04) but showed the same order of magnitude. The diameter was 12, the same as CCBH-2019 and PAO1-2020 and higher than PAO1-2011, which was 9. The average shortest path distance was 4.67, higher than PAO1-2020 (4.01) but slightly lower than CCBH-2019 (4.80). Similar to the previous GRN, CCBH-2022 was disconnected, showing one large connected component (3102 genes) and 20 small connected components.

The degree distributions of the four networks can be seen in Fig. 2A-D, with A and B being the incoming and C and D the outcoming degree distribution. Fig. 2B, D is on a log-log axis, and the straight line is consistent with a power-law distribution. For k-in, the estimated value for γ was 2.79, within the range 2<γ<3, consistent with a power law distribution. For PAO1-2020, the corresponding value was 2.67, 2.89 for CCBH-2019 and 2.71 for PAO1-2011.

Fig. 2:

graphical representation of structural measurements of CCBH-2022 (red) compared to the previously published networks: PAO1-2011 (purple), CCBH-2019 (orange), and PAO1-2020 (green). (A-B) incoming degree distribution of the four gene regulatory networks (GRNs); (C-D) outgoing distribution of the four GRNs. The distributions are plotted on a linear (A, C) and on a logarithmic scale (B, D); (E) local clustering coefficient distribution; (F) clustering coefficient by degree.

The distribution of local clustering coefficients can be seen in Fig. 2E. CCBH-2022 had a global clustering coefficient equal to 4.42e-03, higher than PAO1-2020 (3.03e-03). The scatter plot in Fig. 2F shows the correlation between the local clustering coefficient C(i) and the degree k(i).

The most frequent mode of regulation in CCBH-2022 is activation, occurring in 70,2% of the total interactions in the network, followed by roughly 12% of repression mode and 17.8% of dual or unknown mode. Autoregulation occurs when a gene regulates its expression, and the prevalence in CCBH-2022 is of negative autoregulatory motifs.

The most abundant motif in all four networks was the coherent type I FFL, with 239 in CCBH-2022 (PAO1-2011: 82; CCBH-2019: 79; PAO1-2020: 226). In addition, there were 10 incoherent type II FFL motifs in CCBH-2022 (PAO1-2011: 3; CCBH-2019: 4; PAO1-2020: 8).

Table II shows the 30 most influential hubs in CCBH-2022 and PAO1-2020.

TABLE II

The 30 most influential hubs of CCBH-2022 and PAO1-2020

CCBH-2020			PAO1-2020
Gene	Total number of connections (k-out)	Function	Gene	Total number of connections (k-out)
rpoD	740	Control of expression of housekeeping genes 70	rpoD	749
rpoN	650	Nitrogen metabolism, adhesion, quorum sensing (QS), biofilm formation 71	rpoN	658
algU	353	Positive regulation of response to oxidative stress 72	algU	357
sigX	298	Positive regulation of cell growth 73	sigX	319
rpoS	278	QS, Biofilm, virulence, antibiotic resistance 74	fliA	281
fliA	270	Adhesion, flagellin biosynthesis 75	rpoS	271
rpoH	184	Heat-shock response 76	rpoH	194
gacA	121	Monolayer and biofilm formation 77	gacA	128
algR	119	Cell motility, biofilm formation 78	algR	122
amrZ	109	Cell motility, biofilm formation 79	amrZ	115
lasR	106	QS, regulation of elastin catabolic process 80	lasR	95
fleQ	92	Regulation of mucin adhesion and flagellar expression 81	pvdS	91
fur	88	Control of expression of siderophores and exotoxin A 82	sphR	90
pvdS	87	Iron metabolism, pyoverdine, virulence 83 , 84	fleQ	85
sphR	74	Sphingosine catabolic process 85	fur	69
mvfR	65	QS, regulation of lyase activity, control production of virulence factors 86	vqsM	65
vqsM	61	QS, control production of virulence factors 87	mvfR	62
anr	58	Regulation of oxidoreductase activity 88	pchR	57
rhlR	56	QS, regulation of lipid biosynthetic and proteolysis 89 , 90	anr	53
mexT	53	Antibiotic efflux pump 91	mexT	51
pchR	47	Regulation of pyochelin siderophore, ferripyochelin receptor synthesis 92	argR	46
argR	46	Controls arginine uptake and metabolism 93	fecI	44
gbdR	44	Regulation of cellular amino acid metabolic process 94	gbdR	42
pmrA	43	Antibiotic efflux pump 95 , 96	rhlR	40
fecI	41	Regulation of iron ion transport 97	phoB	40
soxR	40	Antibiotic efflux pump 98	pmrA	40
phoB	40	Cell motility, regulation of cellular response to phosphate starvation 99 , 100	soxR	39
vfr	37	QS, exotoxin A regulator, cell motility 101	dnr	34
dnr	34	Regulation of nitrogen compound metabolic process 102	himA	30
rsaL	34	QS, biofilm formation, regulation of virulence factors 103	himD	30

An analysis was performed to determine whether the hubs are interconnected through direct interactions (Fig. 3).

Fig.3:

connectivity relationships among the 30 most influential hubs of CCBH-2022. Yellow circles indicate regulatory genes considered hubs, light blue circles indicate target genes, black lines indicate an unknown mode of regulation, green lines indicate activation, and red lines indicate repression. Purple lines indicate a dual-mode of regulation.

DISCUSSION

The reconstruction and analysis of the P. aeruginosa GRN contribute to a better understanding of its antibiotic resistance mechanisms. It also contributes to a greater knowledge of related cellular behaviours, such as adaptation and pathogenicity, mainly based on an MDR strain such as CCBH4851.

In this work, we have good coverage of roughly 50% of the genome on this updated network. The genome of reference strain PAO1 has 6.2Mbp, and PAO1-2020 has a coverage of 50% as well, with 5040 interactions and 3006 genes. However, considering that the CCBH4851 genome has 6.8Mbp and has 5452 edges, and 3186 nodes, we can affirm that, to the best of our knowledge, this study presents the largest GRN of P. aeruginosa that has been assembled to date.

On the structural aspects, the charts in Fig. 2 and data in Table I make clear that CCBH-2022 represents a substantial improvement in terms of network completeness and complexity when compared with the previous P. aeruginosa GRNs since it includes more nodes, edges, and network motifs, and when comparing clustering coefficients (Fig. 2E-F). For the in silico approach, the network structural analysis is essential to understand the network architecture and performance.

The structural differences between CCBH-2022 and PA01-2020 results from additional information available due to the new version of PAO1 and recent experimental work on characterising the complete closed genome of P. aeruginosa CCBH4851.

The structural measures of CCBH-2022, such as node degree distribution and clustering coefficient, are consistent with a qualitative description of a scale-free network type. Indeed, the degree distribution followed the power-law distribution (Fig. 2B, D): a small number of nodes had many connections (the hubs) and many nodes had few connections.

The local clustering coefficient and node degree correlation (Fig. 2F) showed that nodes with lower degrees had greater local clustering coefficients than nodes with higher degrees. These characteristics are representative of several biological processes, e.g., RNA binding. ,

CCBH-2022 showed a lower density value than PAO1-2020. The density of both GRNs was low due to the dynamic and structural flexibility of the networks, a characteristic typical of natural phenomena-based networks, and because the nodes were not all interconnected. However, CCBH-2022 density was lower probably because it has 20 small connected components disconnected from the larger one (Fig. 1), while PAO1-2020 had 12 separated components. The variation in the number of connected components is plausible due to their size difference and the biological information about interactions available for the reconstruction.

All the previous P. aeruginosa GRNs are disconnected graphs, showing one large connected component and a separated few small connected components, and there may be several reasons for this disconnection in specific points. According to Medeiros et al., interactions among all genes are not expected since some genes in an organism are independent of each other, compartmentalised or global, constitutive or growth phase-dependent, and are triggered in different growth phases, thus resulting in a disconnected network, which corroborates with the observed low density. The reason can also be from loss of existing interactions or a gain of interactions still not fully described from additional strain-specific blocks of genes acquired by horizontal gene transfer. The large number of connected components found in CCBH-2022 results from connectivity parameters and the global clustering coefficient. Both structural measures are affected by the same biological behaviours.

CCBH-2019 presented more negative regulations than PAO1-2011, a trend that continued between CCBH-2022 and PAO1-2020. Also, the most frequent regulatory activity in CCBH-2022 is activation, but ~50% of the autoregulation was negative, which may be a consequence of the increase in negative autoregulation in the overall network interactions compared to the previous ones. Negative auto-regulation in biological systems is commonly observed. The Escherichia coli GRN exhibited the same pattern, with negative autoregulation prevailing concurrently with the positive regulation in the overall network. The continuity of biological processes is ensured by positive autoregulation. For example, quorum sensing, biofilm formation, secretion of toxins, virulence, and resistance factors production, once initiated, must reach a final stage to have the expected effect. In CCBH-2022, genes involved in these processes, such as lasR, rhlR, pvdS, algU, dnr and anr, have positive autoregulation (and are amongst the 30 principal hubs).

Negative cycles are also crucial for life-sustaining cyclic processes such as metabolic processes and cellular homeostasis. In CCBH-2022, genes involved in arginine metabolism (iscR, desT, lexA, hutC, and mvat) showed a predominance of negative mode of autoregulation. Negative autoregulation is associated with cellular stability. It rapidly responds to variations in concentrations of proteins, toxins, and (or) metabolites to avoid undesired effects such as the energy cost of unneeded synthesis. In CCBH-2022, algZ (transcriptional activator of AlgD, involved in alginate production), lexA (involved in the SOS response), metR (involved in swarming motility and methionine synthesis), , ptxR (affects exotoxin A production) and rsaL (quorum-sensing repressor) presented negative autoregulatory interactions. Autoregulation is common among genes positioned upstream in GRN with crucial developmental functions. ,

The FFL motifs are essential for the modulation of cellular processes according to environmental conditions. CCBH-2022 has 968 FFL motifs, which are patterns of structural structures, while PAO1-2020 has 702. There are 239 coherent type I FFL motifs in CCBH-2022, an abundant presence. According to Mangan and Alon these motifs act as sign-sensitive delay elements, i.e., a circuit that responds rapidly to step-like stimuli in one direction (ON to OFF) and as a delay to steps in the opposite direction (OFF to ON); the temporary removal of the stimulus ceases transcription, so the activation of expression requires a persistent signal to carry on. The incoherent type II FFL motif was less represented but also found in all the GRNs, with a total of 10 in CCBH-2022. Contrastingly with the coherent FFL, the type II FFL acts as a sign-sensitive accelerator, i.e., a circuit that responds rapidly to step-like stimuli in one direction but not in the other.

One last characteristic revealed by the structural analysis was the presence of hubs. The hub’s network (Fig. 3) shows the connection among their interactions; they are all interconnected and belong to the largest connected component of the GRN (Fig. 1A). This connectivity reflects the importance of the influential genes. The hubs can be considered the basis of the GRN. They are crucial in searching for potential drug targets for developing new drugs, as in direct interaction with their specific targets or for an indirect interaction with the subsequent process regulation triggered by them. CCBH-2022 hubs are mainly associated with efflux pump mechanisms (mexT, pmrA, soxR), , , alginate biosynthesis (algU, algR, rpoN), and biofilm formation (rpoN, rpoS, gacA, amrZ).

Table II shows the 30 hubs of PAO1-2020. They are very similar to CCBH-2022 hubs, with some changes in the k-out connections. However, two CCBH-2022 hubs were not in the 30 most influential hubs of PAO1-2020: vfr, a global virulence factor regulator that directly regulates 37 genes, and rsaL, associated with bacterial tolerance to antibiotics, including ciprofloxacin and carbenicillin which directly regulates 34 genes. In PAO1-2020, vfr directly regulates only 12 target genes, while rsaL regulates 19, being one of these an exclusive PAO1-2020 interaction. There are 25 exclusive vfr interactions in CCBH-2022 and 16 exclusive rsaL interactions compared with PAO1-2020. These distinctions can be explained by the fact that the P. aeruginosa CCBH4851 strain is more virulent and multidrug-resistant and also because CCBH-2022 is 7.6% larger in the number of regulatory interactions (5452) than PAO1-2020 (5040), and 20,7% larger in the number of regulatory genes (218) than PAO1-2020 (173). The tables containing these exclusive interactions are in the Supplementary data. These facts strongly indicate that the operation of the main network hubs is not identical. The functioning of CCBH4851 is different, probably due to the greater influence of these two critical genes associated with multi-drug resistance and antibiotic tolerance mechanisms.

The Vfr gene’s role in regulating virulence factor production is related to the production of exotoxin A, a toxin that modifies specific target proteins within mammalian cells and induces necrosis in different tissues and organs in MDR P. aeruginosa infections. , The Vfr gene also regulates the las and consequently the rhl quorum-sensing system, two systems that together control the expression of several genes associated with virulence factor production, including alkaline protease, exotoxin A, pyocyanin, and rhamnolipid, as well as critical genes such as rpoS (the 5th most influential gene in CCBH-2022). The signal receptor (R gene) is one of the essential components of the las and rhl QS systems. It is necessary for coding the transcriptional activator protein (R protein). The lasR and rhlR genes are among the 20 principal hubs. In interactions present only in CCBH-2022, vfr regulates genes associated with virulence and alginate production.

The rsaL gene has an important role in P. aeruginosa as a global regulator of quorum sensing, virulence, and biofilm formation. , Fan et al. showed that the mutation of rsaL increased bacterial tolerance to ciprofloxacin and carbenicillin. In interactions present only in CCBH-2022, rsaL regulates mostly genes of the phz1 and phz2 clusters, showing the control that rsaL has on phenazine expression in P. aeruginosa, driving the production of phenazine-1-carboxylic acid (PCA) which is further converted in the virulence factor pyocyanin. , The pyocyanin production contributes to bacterial tolerance to ciprofloxacin and carbenicillin.

Pseudomonas aeruginosa evades antimicrobial activity during treatment and exerts antimicrobial resistance by mainly intrinsic resistance mechanisms. Examples of resistance mechanisms are multi-drug efflux pumps, biofilm synthesis, enzymatic inactivation/degradation, drug permeability restriction, production of beta-lactamases, acquired resistance by a mutation in drug targets, and acquisition of resistance genes via horizontal gene transfer.

There is a directed regulatory connection from alginate biosynthesis to iron metabolism and some antibiotic resistance mechanisms. The algU, algR, rpoN, pvdS, and fecI genes are related to these processes , and are among the most influential hubs.

Pseudomonas aeruginosa has multiple efflux pump systems that prevent the antimicrobial agents from accumulating in adequate concentration to cause an effect in the cell, extruding the drug out. Efflux pump systems are associated with resistance to beta-lactams, fluoroquinolones, tetracycline, chloramphenicol, macrolides, and aminoglycosides. Differential expression or mutations of efflux system genes are also contributing factors for carbapenem and aminoglycoside resistance. The mexT, pmrA, soxR genes, related to multidrug antibiotic efflux pumps, are also amongst the most influential hubs.

The fleQ gene is also among the hubs and affects psl (polysaccharide synthesis locus) genes and regulates the efflux pump genes, mexA, mexE, and oprH, by brlR. , The psl cluster comprises 15 exopolysaccharide biosynthesis-related genes organised in tandem that are important for biofilm formation.

The mexT and soxR genes positively regulate an efflux pump system, and several virulence factors, , and pmrA regulate efflux pumps and the polymyxin B and colistin resistance. , ,

Efflux pumps also help biofilm formation. Biofilms are also related to protection from the host immune system and antibiotic penetration and tolerance, preventing them from entering the microbial population and inhibiting its action as a first-line defense mechanism. , , The rpoN, rpoS, gacA, algR and amrZ hubs participate in the regulation of P. aeruginosa biofilm.

This system biology approach to characterise the MDR P. aeruginosa CCBH4851 regulatory network may lead to the development of strategies to disrupt the connectivity of these essential processes, thus, possibly decreasing the pathogenicity and suppressing the resistance of this bacterium.

In conclusion - This manuscript reports the reconstruction and structural analysis of the largest P. aeruginosa regulatory network available in the literature to date. This work can give new insights into identifying novel candidate antibiotic targets and contributes to an increase in our understanding of the behaviour of this bacterium.

This network’s dynamic model construction is one of our future studies, intending to help researchers working on experimental drug design and screening. The goal is to predict the dynamic behaviour better and improve the understanding of P. aeruginosa, allowing the simulation of normal and stress conditions to discover potential therapeutic targets and help develop new drugs against P. aeruginosa bacterial infection.