Introduction

Transcription factors (TFs) are key regulators of transcription and thus play crucial roles in mediating multiple biological pathways and events in eukaryotes and prokaryotes. Based on their DNA-binding domains (DBDs), TFs can be classified into different families, such as the LysR, AraC, and LuxR families¹. TFs can activate and repress the expression of a set of genes in response to various environmental and/or specific signal triggers^2,3. Several approaches can be used to study TFs, including chromatin immunoprecipitation sequencing (ChIP-seq)⁴ and ChIP-exo⁵ in in vivo research and DNA affinity purification sequencing (DAP-seq)⁶ and high-throughput systematic evolution of ligands by exponential enrichment (HT-SELEX)⁷ in in vitro research. ChIP-seq is a powerful in vivo technique that enables more accurate genome-wide mapping of TF–DNA interactions than in vitro methods because TFs may interact with other co-regulators in an environment-specific pattern, thereby altering binding preferences in vivo^8,9.

These sequencing approaches have been used to illustrate the biological functions of TFs in several bacterial species, including Mycobacterium tuberculosis, Vibrio cholerae, Salmonella enterica, Escherichia coli, and Clostridium thermocellum^10–15. In addition, several databases and models have been constructed based on experimental validation or machine learning for studying TFs in bacteria, such as CollecTF, PredicTF, PRODORIC, DBTBS, and RegulomePA^16–20. However, the functions and spatiotemporal interactions of most bacterial TFs remain to be clarified.

The gram-negative bacterium Pseudomonas aeruginosa is a major human opportunistic pathogen that grows ubiquitously. P. aeruginosa tends to cause infection in burn victims, patients with cystic fibrosis (CF), and individuals hospitalised for long periods^21,22. It utilises versatile pathways to exert its virulence, such as biofilm formation, quorum sensing (QS), type III and type VI secretion systems (T3SS and T6SS, respectively), motility, siderophore production, oxidative stress resistance, and antibiotic resistance, all of which are under the control of a complicated TF regulatory network (TRN). The metabolic changes regulated by the TRN enhance the pathogen’s adaptability under different infection and survival conditions.

QS, a major virulence pathway in P. aeruginosa, is a bacterial cell–cell communication mechanism involving the expression of many virulence genes. Its regulation occurs hierarchically, comprising the interconnected las, rhl, and pqs systems^23,24. Through its transcriptional regulator, LasR, and in collaboration with its corresponding autoinducer, 3-oxo-C₁₂-HSL (encoded by lasI), the las system triggers the activation of the rhl and pqs systems^23,25. In addition, P. aeruginosa attaches to and colonises surfaces primarily using two forms of motility: (i) swarming and swimming powered by flagella, which are multicellular and individual cell movements, respectively, and (ii) twitching powered by type 4 pili (T4P)²⁶, wherein T4P promote surface attachment by twitching to pull the cell closer to the attachment sites²⁷. The motility of P. aeruginosa is thus the initial step in its pathogenic process and is associated with the upregulation of virulence and the induction of host defence. Overall, the virulence-related pathways, including QS and motility, are controlled by the sophisticated TRN in P. aeruginosa.

Despite the TRN being essential for understanding how bacterial genes are regulated synergistically, only a limited number of binding targets of TFs have been identified in P. aeruginosa over the past decades. For example, ChIP-seq and RNA sequencing (RNA-seq) co-analysis has been used to map a P. aeruginosa genomic regulatory network (PAGnet) of 20 key virulence-related TFs, while HT-SELEX has been used to characterise the DNA-binding specificities of 182 TFs²⁸. DAP-seq has been used to reveal a regulatory network of 55 response regulators (RRs) of two-component systems in P. aeruginosa, 51 of which are from the PAO1 strain²⁹. However, almost half of the TFs remain to be annotated with biological functions, and their downstream targets and upstream regulators remain unknown. To fill this knowledge gap, we performed ChIP-seq of 172 TFs (Table S1) from P. aeruginosa, which revealed a hierarchical and co-association regulatory network among TFs and newly identified virulence-related regulators. We also established a Web-based database to store and search TF-binding patterns together with ChIP-seq and HT-SELEX data, which is available to the public. Our findings and tools will significantly contribute to future mechanistic studies on and a comprehensive understanding of the TRN in P. aeruginosa and other bacterial species.

Results

ChIP-seq of 172 uncharacterised TFs reveals a global transcriptional atlas of P. aeruginosa

In this study, 373 proteins in P. aeruginosa were annotated as TFs with DNA-binding activity based on the existing annotations in the Pseudomonas Genome Database³⁰ (Fig. 1A). To profile the whole TF–DNA binding landscape and construct a comprehensive regulatory network of P. aeruginosa, ChIP-seq experiments were performed on 172 TFs to determine the TF-binding sites (TFBSs) (Fig. 1A). Libraries were constructed using vesicular stomatitis virus glycoprotein tagging, and raw sequence reads were mapped to the PAO1 genome using bowtie2 (Version 2.3.4.1)³¹. Using MACS2³², 81,009 significant binding peaks with a cut-off p-value of 0.001 were identified. Subsequently, the R package ChIPpeakAnno³³ was used to define the peak locations and find the nearest genes. The peak location can be divided into six features: upstream, overlapStart, inside, overlapEnd, downstream, and includeFeature. The percentages of binding peak locations for each TF indicated that more than half of the TFs preferentially bind to the upstream (upstream intergenic regions) and overlapStart (overlapped with the transcription start sites [TSSs] of genes) regions (Fig. 1B). Apart from these regions, some TFs bind to the inside regions (coding regions of genes). Few TFs were found to preferentially bind to the rest of the peak location features.

Figure 1.

To further provide insights into the functional roles of different TFs in gene regulation, we analysed the peak distance to the TSS, which revealed several patterns (Fig. 1C). First, most of the TFs exhibited peak distributions that were concentrated in the proximal region around the TSS. Some TFs, such as ArsR and AraC, showed a broader distribution of peak distance, indicating that these TFs might also be active at positions farther from the TSS. Furthermore, certain TFs, including CRP and DeoR, displayed distinct bimodal or multimodal distributions, while other TFs, such as Fis and AraC, showed smoother and more narrowly focused distributions, suggesting their different binding preferences. The number of peaks of each TF and its DBD types are summarised in a treemap in Fig. 1D. TFs in the LysR and AraC families accounted for almost half of all of the peaks. PA2718, classified in the MerR family, yielded the highest number (2,480) of binding peaks, followed by the TF PA0756 (2,080 binding peaks) from the OmpR family. Furthermore, eight TFs showed no more than 10 peaks, and the TFs PA4074 and PA4806 had only two peaks each. To preliminarily identify the different regulatory functions of TFs via their different binding preferences, we performed Gene Ontology (GO) functional enrichment analysis for each TF and presented the top 10 GO terms (Fig. 1E). The significantly enriched GO terms (BH-adjusted P < 0.05) for all of the 172 TFs analysed by ChIP-seq are listed in Table S2, and these terms suggest divergent functions of the TFs. Of these functional categories, most were related to metabolism, such as transcription, translation, ribosome structure, and GTPase activity, while some of the functional categories were associated with virulence, including T4P, O antigen biosynthesis, and biofilm formation (Fig. 1E).

Hierarchical networks of TFs based on pairwise interactions

TF-coding genes can also serve as targets for other TFs, leading to interactions between the DNA-binding profiles of TFs. To clarify these relationships, lists of binding targets of the TFs obtained using ChIP-seq and HT-SELEX were used to assemble a system-level hierarchical regulatory network. After filtering by TFs that bind to promoter regions, 947 unique pairwise patterns between TFs were identified from a total of 13,375 promoter interactions.

We next defined a statistic index, h, based on the out-degree (O) and in-degree (I) of interactions to measure the direction and hierarchical level of TFs [h = (O-I)/(O+I)] in our established hierarchical regulatory network and thus identify hubs and information-flow bottlenecks³⁴. The h value of each TF can reflect the direction and extent of information flow. A positive h indicates that the TF acts ‘upstream’, whereas a negative value suggests that it acts ‘downstream’. The higher the h value, the stronger the implication of the direction. A density plot of all of the TFs with their h values revealed that they can be approximately divided into three groups: top (h ∈ [0.75, 1]), middle (h ∈ [-0.75, 0.75]), and bottom (h ∈ [-1, 0.75]) levels (Fig. S2A). The top-level TFs are those that tend to control many other TFs, while the bottom-level TFs are those that are more regulated by other TFs than acting as regulators themselves. The middle-level TFs are those that connect the top- and bottom-level TFs. Of the 231 TFs in this hierarchy, 100, 46, and 85 TFs were defined as top-, middle-, and bottom-level TFs, respectively. The complete hierarchical information of all of the TFs is provided in Table S3. To avoid extreme cases and aid visualisation, we summarised the hierarchical regulatory network by removing TFs for which (O+I) was more than 10 (Fig. 2A and Fig. S2B). The network was visualised using Cytoscape software³⁵.

Figure 2.

Additionally, we investigated the factors influencing the hierarchical level of TFs, such as DBD types (e.g., DBD families containing more TFs). TFs from the AraC, GntR, and TetR families tended to cluster at the top level. TFs from the MarR family and two-DBD-type TFs (i.e., TFs with two DBDs) were mostly found in the bottom level. Additionally, among the characterised TFs, 58% were clustered at the bottom level, probably because these TFs directly regulated the downstream genes involved in important phenotypes, such as AlgR, AmrZ, MvfR, FleQ, and VqsM. Overall, the organised hierarchical regulatory network profiling at the three levels showed complex pairwise interactions among TFs in P. aeruginosa.

Ternary regulatory motifs show flexible relationships among TFs in P. aeruginosa

Apart from their global hierarchical regulatory structure, we investigated the constituent regulatory network motifs of TFs, which revealed small connectivity patterns associated with canonical functions³⁶. Positive and negative feedback loops are prevalent in bacterial regulatory networks, which generally involve combinations of several genes to appropriately respond to stimuli. One of the regulatory patterns was a single auto-regulator that can self-regulate its expression or activity. Through analysis, nine auto-regulator TFs were identified (Fig. 2B). Compared with non-auto- regulators, auto-regulators usually tend to be repressors, which play important roles in maintaining a steady state^37,38. For example, the key T3SS activators PsrA and HrpL can repress their own expression in P. aeruginosa and P. savastanoi, respectively^39–41.

In addition, the most fundamental regulatory pattern was the ternary TF motif. We computed all possible motifs, which are listed in Fig. 2C and Table S4. The basic triangular motifs were monodirectional regulatory structures with five different motifs. The most frequently observed motif (6,535 occurrences) was motif_2, wherein two different TFs co-regulate another TF. The remaining seven types of motifs, which contained toggle switches, occurred much less frequently than the basic motifs. The more complex the relationships among TFs, the lesser their occurrence. In particular, motif_10, motif_12, and motif_13, characterised by two or three mutually regulating TFs, were not found even once (Fig. 2C). The high occurrence of different TFs in diverse motifs revealed a complex and flexible regulatory network in P. aeruginosa. Next, we constructed a basic motif_3 regulatory network based on potential TF–TF interactions to detect the regulatory preferences based on DBD types (Fig. 2D). Basic motif_3 was a typical hierarchical regulatory network in which TF1 regulates TF2, which in turn regulates TF3. We found that the AraC, Fis, and GntR families tended to harbour top-level TFs rather than middle- and bottom-level TFs, while the OmpR family tended to harbour TFs from middle- and bottom-levels more than those from the top level. Furthermore, the LysR and LuxR families revealed non-significant differences among these three TF levels (Fig. 2D). Taken together, our integrated hierarchical regulatory profiling revealed a multifaceted TF–TF connection, reflecting the regulatory preferences of TFs and DBD types in P. aeruginosa.

Clustering of 103 TF-binding motifs

According to their binding sites, a total of 103 binding motifs were determined. The TFs were clustered based on the positional weight matrix similarity of their motifs. A circular phylogenetic tree was constructed using hierarchical clustering of the pairwise similarity matrix (Fig. S3 and Table S5). Each node in the tree represents a TF, colour-highlighted according to its DBD family classification. The clustering analysis revealed four distinct clusters of TFs, one each coloured light green (Cluster 1), light blue (Cluster 2), purple (Cluster 3), and pink (Cluster 4). These clusters indicate groups of TFs with similar binding properties or regulatory functions.

Within each identified cluster, TFs tend to share conserved motif patterns, indicating potential functional relationships. For instance, it was found that the TFs in cluster 1 were enriched in adenine (A) and thymine (T) bases. Further, cluster 4 was more likely to have an inverted repeat (IR) sequence, such as PA0367 and PA1015. This conservation suggests that TFs within the same cluster co-regulate similar sets of genes. Additionally, the clustering uncovered novel associations between TFs from different DBD families, indicating potential new regulatory interactions or evolutionary links that warrant further investigation.

Genome-wide synergistic co-association of TFs in P. aeruginosa

To further investigate the crosstalk among TFs in P. aeruginosa, we computed the co- association scores based on the overlaps between the peaks of all pairs of TFs and ChIP-seq data (Fig. S4A and Table S6). Among a total of 16,175 TF pairs, we found 5,716 significant TF interactions with a co-association score of up to 0.1 based on an elbow statistic (Fig. S4B). The relationships among TFs can be further determined as core clusters if their co-association scores are more than 0.4. The clustering analysis of core networks using Glay clustering⁴² resulted in seven associated modules, each containing a diverse number of TFs with the DBD families highlighted in different colours (Fig. 3A).

Figure 3.

Next, we summarised the potential functional crosstalk among the middle-scale co- associated five TFs (PA2417, PA2718, PA0756, PmiR, and AgtR) in cluster 3. Among a total of 1,241 target genes, 372 were co-bound by at least three TFs, and 35 were co- bound by all five TFs, indicating that these TFs have a high level of collaboration in regulatory patterns (Fig. 3B). There were smaller intersections among different combinations of TFs, highlighting potential co-regulation or shared binding sites. Binding targets co-bound by more than four TFs are summarised in Fig. 3C. The dense network and connectivity indicate a complex regulatory landscape, where these TFs may orchestrate a coordinated regulation of gene expression. For example, peak visualization of PA2504 revealed distinct binding patterns of these five TFs compared with the control group (Fig. 3D). PA2504 is the sole partner of the PppH RNA hydrolase related to several important cellular factors⁴³. In addition, four of the five TFs (barring PA2718) can bind to the promoter region of pqsH, which is involved in the pqs pathway (Fig. 3E). Notably, PmiR has previously been verified to bind to the pqsH promoter region and thereby regulate virulence⁴⁴. PA2504 and pqsH are highlighted in blue and green, respectively, in Fig. 3C.

Twenty-four master regulators are associated with virulence

After elucidating the binding targets of 172 TFs, we investigated their potential functions in the virulence-related pathways and pathogenesis of P. aeruginosa. For this purpose, we performed an enrichment analysis using a hypergeometric test to define the master regulators (BH-adjusted P < 0.05) that might play important roles in regulating specific pathways. We primarily focused on TF-binding profiles involving the promoter regions of genes involved in nine pathways, namely QS, motility, biofilm production, antibiotic resistance, T6SS, T3SS, reactive oxygen species (ROS) resistance, pyocyanin production, and siderophore production. Accordingly, 134 of the 172 TFs were found to bind to at least one target gene associated with the abovementioned virulent pathways, and 24 were identified as master regulators of genes involved in six pathways, namely siderophore production, biofilm production, pyocyanin production, QS, motility, and ROS resistance pathways (Fig. 4A and Table S7).

Figure 4.

The intersection of master regulators across these virulent pathways revealed that six TFs (PA0167, PyrR, PA0707, PA0877, PA1504, and PA1484) only played an important role in regulating the motility pathway, and four TFs (PA1380, PA0815, PA5428, and PA3973) jointly regulated the motility and other pathways, including biofilm production, pyocyanin production, and QS pathways (Fig. 4B). Specifically, PA0815 co-regulated motility, biofilm production, and QS pathways. PA1380 jointly regulated motility and pyocyanin production pathways. PA3973 and PA5428 jointly regulated motility and biofilm production pathways. The regulatory network of these four master regulators and their target genes is summarised in Fig. S5, represented by lines connecting the TFs to the targets involved in specific pathways.

Predicted functions of regulators characterised in two metabolomic pathways

P. aeruginosa uses different global regulators associated with metabolism that ensure its survival and enhance its adaptability under fluctuating environmental conditions. For instance, it uses different regulatory mechanisms to respond to nutrient changes, oxygen limitation, or CF-associated lung infection conditions^45–48.

To define the core TFs that regulate metabolism, we performed master regulator analysis in metabolic pathways (p.adjust < 0.05) based on the TFBSs located in the promoter regions revealed by ChIP-seq for all TFs. We identified 14 TFs potentially involved in two pathways, namely the tricarboxylic acid (TCA) cycle and ribosome function pathways (Fig. S6A-B). The radial plot highlights the master regulators within the TCA cycle and ribosome function pathways, and their joint regulation with their targets is presented in the regulatory network below. PA2957 and PA5438 were found to play putative roles in the TCA cycle by potentially binding to the promoter of certain genes, including icd, sucA, and sdhC, indicating that they coordinate the expression of genes crucial for central metabolism and energy production (Fig. S6B). Furthermore, certain TFs, such as PA0611, PA0403 (PyrR), and PA1283, were found to be involved in regulating a broad array of ribosomal genes, such as rpsB, rpsA, and rplL. Notably, TF PyrR was found to play a putative role in regulating motility in P. aeruginosa. Overall, these regulators, which are related to the core metabolic pathways that enhance the adaptation of P. aeruginosa to different infection conditions and contribute to its pathogenicity, could offer potential drug targets for tackling P. aeruginosa infection in the future.

A global summary of transcription regulatory networks and functions in P. aeruginosa

After studying the TFBSs of most of the TFs in P. aeruginosa both in vivo and in vitro, we organised these datasets to yield a comprehensive connection of virulence-related master regulators. The QS and motility pathways consisting of 18 and 20 TFs, respectively, were the key pathways in the virulent landscape and showed a close collaboration with the other five pathways (Fig. 5A). Four TFs (ExsA, CprR, AlpR, and PA2449) were involved in T3SS, a direct host infection mechanism used by P. aeruginosa. Additionally, six TFs (PA4436, PA0929, PA3771, PA2534, PA3699, and PA5218) were involved in siderophore production, as identified using HT-SELEX, and showed correlations with only one TF, PA2534, involved in the ROS resistance pathway. Four TFs (CprR, PA0708, PA1520, and PA2479) were characterised in the stringent response and persister cell pathways and were correlated with other virulence-related pathways, namely QS, ROS resistance, T3SS, and antibiotic resistance pathways.

Figure 5.

To investigate the phylogenetic classification of the binding targets of the identified TFs on a genome-wide scale, we derived a global atlas revealing the relationships between the protein functions of target genes and the family categories of TFs (Fig. 5B). The Clusters of Orthologous Groups of proteins (COGs) database was used to classify proteins based on the orthology concept⁴⁹. Based on the annotated COGs of P. aeruginosa, we found that most of the TF targets were involved in metabolism, followed by poorly characterised categories. The P. aeruginosa PAO1 genome is 6.3 Mbp long, and approximately 40% of it remains uncharacterised but is believed to encode well-conserved hypothetical proteins (HPs), which might have indispensable and similar functions. Given that most of the TF targets are poorly characterised, it is crucial to study these HPs in the future to provide more insights into P. aeruginosa pathogenicity. Altogether, our results provide a resource of TFs and their putative target genes in P. aeruginosa and, consequently, a basis for understanding the regulatory network of virulence and orthological functions.

A Web-based database of all TFBSs in P. aeruginosa

To enable a quick search of the TFBSs in P. aeruginosa, we developed a Web-based database containing our ChIP-seq and HT-SELEX results at https://jiadhuang0417.shinyapps.io/PATF_Net/, which is available to the public. It has four parts: homepage, data page, upload, and contact. The homepage briefly describes the background information and schematic workflow of ChIP-seq and HT-SELEX. The data page provides two types of data: a network plot of TFs and targets, and detailed TF–target information in a table format. The network plot illustrates the top 50 peaks for ease of visualisation. The code was generated using R, and the figures were created using BioRender (https://www.biorender.com/).

Overall, this platform provides several services, including individual TF–target network visualisation, detailed peak annotations in a table format, TF–target interaction visualisation, data download, and data upload for those who want to share their results on this online platform. The primary purpose of PATF_Net is to store, search, and mine valuable information on P. aeruginosa TFs for researchers investigating P. aeruginosa infection in the future.

Conservation and evolution of TFs in P. aeruginosa

Differences between species are mainly due to differences in genes, whereas essential genes are usually well conserved, especially in subspecies. To investigate the conservation of TFs in P. aeruginosa, we performed a pan-genome analysis of 100 strains containing model strains and clinical strains using the Roary software⁵⁰. A total of 21,432 genes were identified across the 100 P. aeruginosa strains, of which only 4,544 genes were classified as core genes (existing in more than 99% of strains), revealing that most of the genes were cloud genes (existing in less than 15% of strains) (Fig. 6A). However, when focusing on the TFs among these 100 strains, almost 89% of TFs (331 of 373) were found to be well-conserved (Fig. 6B), which implied that these recognised TFs from PAO1 play essential roles in biological activities in all of the P. aeruginosa subspecies.

Figure 6.

The bar plot provides detailed information on the rest of the non-core TFs in PAO1, in addition to several important TFs, such as LasR and ExsA, that are lacking in some strains (Fig. 7C). As is known, the lasR mutant of P. aeruginosa is associated with CF lung disease progression and is considered a marker of an early CF-adaptive phenotype⁵¹. ExsA is considered an important activator of the T3SS pathway in P. aeruginosa, which is characterised by an acute infection phenotype. Additionally, VqsM might have a relatively special function in PAO1 among P. aeruginosa subspecies in addition to playing an important role in the QS and antibiotic resistance pathways. Taken together, the core- and non-core TFs were both vital and worth studying for understanding the metabolism and virulence-related regulatory mechanisms in P. aeruginosa.

To further investigate the conservation of TFs in P. aeruginosa, we re-analysed the ChIP-seq data of PhoB and RpoN in the PA14 strain^52,53 and compared them with those in the PAO1 strain. Fig. S7A shows genome-wide binding peaks of PhoB and RpoN throughout the genomes of PAO1 and PA14. We found that in PA14 and PAO1, PhoB had 725 and 63 peaks, while RpoN yielded 363 and 1,238 peaks, respectively. Next, we used MEME⁵⁴ to compare the motifs of PhoB and RpoN in PAO1 and PA14 based on their binding sequences from ChIP-seq (Fig. S7B). PhoB and RpoN shared a well-conserved motif in the two strains even though they showed different peak distributions in the genomes. They had different numbers of peaks, and of the 690 and 345 respective annotation targets of PhoB and RpoN in PA14, 114 and 66 targets, respectively, were unique to PA14. The rest of the targets were homologous between PAO1 and PA14. We next performed target annotation and observed that the genes regulated by these two TFs were significantly (Fisher’s test) overrepresented in the two tested strains (Fig. S7C). For example, RpoN has been found to play an important role in regulating T6SS in both PAO1 and PA14^53,55. Overall, the conservation of the representative TFs PhoB and RpoN in PAO1 and PA14 indicates that conserved TFs may share similar functions in different species. This knowledge is expected to guide drug development for different infectious strains.

TFs play a pivotal role in regulating gene expression, and studying their conservation and evolution could provide insights into how gene expression patterns are conserved or modified across different species and help to predict gene regulatory networks in different species. For example, the conservation of a TF across multiple species suggests that its targets and regulatory interactions are likely to be conserved as well. Therefore, to study the conservation and evolution of P. aeruginosa TFs across other species, we conducted phylogenetic analyses of the 373 TFs of P. aeruginosa across other bacteria, fungi, archaea, plants, and animals. We found that three P. aeruginosa species showed the most conserved TFs, followed by P. savastanoi, a phytopathogen (Fig. 6D). Compared with bacteria, the conservation values of the TFs were lower in fungi, archaea, plants, and animals. However, several orthologues of TFs were still noted among these organisms. For instance, PA2032 had more than 25% and approximately 30% identity with aminoadipate aminotransferase in humans and mice, respectively. Altogether, the high-level evolutionary conservation of TFs in P. aeruginosa strongly suggests that the regulatory mechanisms are common to a wide range of bacteria.

Discussion

Bacterial pathogens use many mechanisms to establish infection and cause diseases in human hosts, making them a major public health concern worldwide. These bacteria can secrete a wide range of molecular particles that recognise and bind to host cell targets, thus damaging or preventing host responses⁵⁶. Pathogenic bacteria use strategies to adapt to environmental stressors such as nutrient limitation or antibiotic treatment⁵⁷. The underlying mechanism of these strategies involves sensing external signals and then responding to different environmental conditions. TFs are key molecules involved in responding to host or environmental signals by precisely modulating transcription levels; hence, their pivotal roles in infection conditions cannot be overstated. Nevertheless, our comprehension of bacterial TFs remains considerably limited, even in extensively examined model organisms such as P. aeruginosa. To address this knowledge gap and shed more light on TFs in P. aeruginosa, we conducted ChIP-seq experiments aimed at elucidating the TFBSs for 172 TFs with relatively unexplored biological functions.

The basic function of TFs is to recognise specific DNA sequences (i.e., motifs) at the promoter region and then upregulate or downregulate the transcription of the target gene⁵⁸. Apart from promoter regions, TFs have also been identified to directly regulate the transcription of coding regions in bacteria⁵⁹. ChIP-seq is a valuable technique to obtain high-resolution DNA–protein interaction mapping on a genome- wide scale to study the potential binding targets of TFs. The technique is highly sensitive and can detect low-abundance DNA-binding events, allowing the detection of weak TF–DNA interactions. In vitro assays, such as DAP-seq or HT-SELEX, offer superior scalability to ChIP-seq but present their own challenges^29,60. It is difficult to investigate the crosstalk and interactions between TFs by studying proteins individually. Thus, we devised a strategy involving TF-overexpressing plasmids. This strategy can enhance the visualisation of more subtle binding peaks, particularly when investigating co-associations and TFs lacking well-defined signals. Notably, the literature underscores the collaborative nature of TFs, which often exhibit interactions with other auxiliary proteins at multiple promoter sites⁶¹. As shown in Figures 3 and S3, we present a TF-binding motifs clustering tree and a TF co-association regulatory network based on ChIP-seq data and core clusters of TFs with high co-association scores. Overall, our above-described analysis has illustrated a regulatory atlas to study the biological functions of P. aeruginosa TFs in specific contexts.

P. aeruginosa uses TFs to control biological functions, including metabolism and virulence. Given the public health burden of P. aeruginosa infection, many studies have investigated the roles of TFs in this pathogen. These studies have characterised the regulatory mechanisms of individual TFs, such as FleQ, AlgR, LasR, VqsM, and PvrA^62–66. Additionally, a TRN comprising 690 genes and 1,020 regulatory interactions with six biological modules and main motifs was established in 2011 using published data⁶⁷. Furthermore, we developed a PAGnet containing 20 key virulence-associated TFs based on ChIP-seq and RNA-seq data⁶⁸. A previous study used machine learning to identify modules of a group of genes related to known TFs from published data⁶⁹. Another study used DAP-seq to determine a regulatory network of 55 RRs in P. aeruginosa²⁹. To further identify and characterise unknown TFs, we used HT-SELEX to describe the binding specificities of 182 TFs²⁸. Compared with the abovementioned studies, the present work profiled a comprehensive regulatory network of P. aeruginosa, combining existing and new data with experimental verification.

We confirmed the functions and provided potential molecular mechanisms of a number of previously identified TFs to explain their phenotypes. PA0797 is known to regulate the pqs system and pyocyanin production⁷⁰. In the present study, it was also found to bind to the pqsH promoter region and its motif was visualised. PA5428 was found to bind to the promoter regions of aceA and glcB genes⁷¹, which was also demonstrated in our ChIP-seq results. PA4381 (CloR) was found to be associated with polymyxin resistance in a previous study⁷² and to be possibly related to ROS resistance in the present study. Furthermore, PA5032 plays a putative role in biofilm regulation and also forms an operon with PA5033, an HP associated with biofilm formation⁷³.

Our study is the first to illustrate the hierarchical regulatory network and ternary motifs of and genome-wide co-association relationships among TFs in P. aeruginosa, contributing to a better understanding of the fundamental traits of pathogenicity of this species. Traditional bacterial eradication strategies involving antibiotics often focus on pathogen extermination, yet this approach can inadvertently foster selective growth advantages and precipitate the emergence of drug-resistant strains. Therefore, prioritising TF-targeted drugs could potentially alleviate selective growth pressure and enhance future preventative measures against P. aeruginosa infection. Overall, we gathered the binding site data of most of the uncharacterised TFs in P. aeruginosa and established a Web-based database for storing and searching this information. The TF regulatory mechanisms clarified in this database can help to decode the pathogenic mechanisms and guide the development of new drugs targeting virulence-related TFs to combat P. aeruginosa infection. The conservation and evolutionary patterns of TFs also suggest that these TFs perform similar functions in diverse species, thus underscoring the broader applicability of our research to studies investigating TFs in other microbes.

Methods

Strains, plasmids, and primers

The bacterial strains, plasmids, and primers used in the present study are listed in Table S8. The P. aeruginosa PAO1 strain, and its derivatives were grown at 37 °C in LB (Luria-Bertani) broth with shaking at 220 rpm or on LB agar plates. Antibiotics were used for Escherichia coli at the following concentrations: kanamycin at 50 μg/ ml and carbenicillin at 60 μg/ml.

ChIP-seq and analysis

The chromatin immunoprecipitation (ChIP) procedures were modified from a previous study⁷⁴. Briefly, for the VSVG-tagged, the ORF was amplified by PCR from the PAO1 genome and cloned into pAK1900 plasmid by HindIII/BamHI for the overexpressed TFs through HindIII site by using ClonEx- press MultiS One Step Cloning Kit (Vazyme, China). Wild-type P. aeruginosa containing empty pAK1900 or pAK1900-TF-VSV-G was cultured in LB medium supplemented at 37°C with shaking until mid-log phase (OD600 = 0.6), then we cross-linked the samples with 1% formaldehyde for 10 min. Subsequently, cross-linking was stopped by the addition of 125 mM glycine. Samples were centrifuged and washed thrice with a Tris buffer (20 mM Tris-HCl, pH 7.5, 150 mM NaCl). The resulting pellets were resuspended in 500 μl IP buffer (50 mM HEPES–KOH [pH 7.5], 150 mM NaCl, 1 mM EDTA, 1 % Triton X- 100, 0.1% sodium deoxycholate, 0.1% SDS, and mini-protease inhibitor cocktail (Roche), and then s The DNA was broken into pieces (100 to 300 bp) with an ultrasonic processor. Insoluble cellular debris was removed by centrifugation at 4 °C, and the supernatant was used as the input sample in IP experiments. We next added 25 μl agarose-conjugated anti-VSV antibodies (Sigma) in the IP buffer. Washing, reverse cross-linked, and purification of the ChIP DNA were conducted. We used the NEXTflex ChIP-Seq kit (Bio Scientific) to construct the DNA fragment library, and agarose gel was used to cut DNA fragments between 150 and 250 bp. After sequencing, the raw reads were trimmed by Trim Galore⁷⁵ with default parameters. Then, the filtered reads were aligned to the PAO1 genome (NC_002516) using bowtie2. Only the unique reads after alignment will be used to perform peakcalling with MACS2 (P < 0.001). MEME-ChIP was used to identify consensus motifs with all peaks. TF target genes were then annotated by the R package ChIPpeakAnno.

Identification of virulence-related master regulators

We first generated gene lists associated with nine pathways: QS, motility, biofilm production, antibiotic resistance, T6SS, T3SS, ROS resistance, pyocyanin, and siderophores. Then, we calculated the number of overlap genes between target genes of TF and genes involved in the above eight pathways. The statistical significance of the master regulator was identified using the Hypergeometric test (BH-adjusted P < 0.05). R package ggplot2⁷⁶ was used to visualize the result.

GO analysis

GO enrichment analysis of TFs target genes was identified using R package clusterProfile⁷⁷. The enriched GO terms with BH-adjusted P < 0.05 were defined as significantly enriched.

Pan-genome analysis

The information on the model and clinical P. aeruginosa strains^78,79 can be found in Table S9. Briefly, the genome sequence of all strains were annotated by prokka⁸⁰ and the outputs were then used as input files for roary.

Statistical analysis

Student’s t-test was performed to analyze the RT-qPCR and lux production results using Microsoft Office Excel 2019. *P < 0.05, **P < 0.01, and ***P < 0.001, and results represent means ± SD. All experiments were repeated at least twice.

Data availability

The ChIP-seq data was uploaded to National Center for Biotechnology Information Gene Expression Omnibus database under accession GSE241603 and GSE271817. Analysis codes have been deposited at https://github.com/dengxinb2315/PS-PATRnet-code.

Acknowledgements

This study was funded by the Health and Medical Research Fund (20190942), Shenzhen Science and Technology Fund (JCYJ20210324134000002), National Natural Science Foundation of China grants (32172358), the General Research Funds of Hong Kong (11103221, 11101722, and 11102223), and the Theme-based Research Scheme (T11-104/22-R). The funders had no role in study design, data collection, interpretation, or the decision to submit the work for publication.

Additional information

Contributions

X.D. designed the project and obtained funding, material support and study supervision. X.D., J.H. wrote the manuscript and generated the figures and tables. J.H. performed bioinformatic analysis. J.H., Y.S., F.C., and J.L performed the experiments.

S.L. helped to write the code. Y.S. and L.H. helped to edit the manuscript. All authors reviewed the manuscript.

Supplementary information

Fig. S1. RNA-seq for all TFs and ChIPed TFs. Previous RNA-seq data shows the expression level of all TFs and ChIPed TFs. The circle size indicates the log2- transformation raw counts, and the color highlighted for ChIPed TFs represented different DBD families.

Fig. S2. Distribution of Hierarchy Height h. h = (O-I)/(O+I), O presents the extent they regulate other TFs, and I indicates the level other TFs regulate them. (A). TFs can be divided into almost equal three levels depending on index h: Bottom level (-1.0 < h < 0.75), Middle level (-0.75 < h < 0.75), and Top level (0.75 < h < 1.0). (B). The number of TFs resided in three levels after removing TFs with O+I <10.

Fig. S3 Circular phylogenetic tree of TF binding motifs. Circular phylogenetic tree representing the clustering of TFs based on their motif similarity. The tree was constructed using hierarchical clustering with the Ward.D2 method applied to the pairwise similarity matrix of motifs. The tree were visualized using ggtree⁸¹. Each node represents a TF, color-coded by its DBD family classification as indicated in the legend on the right. The motif logos are aligned with the respective TFs, with their orientations adjusted to maintain legibility around the circular layout.

Fig. S4. Co-association network of TFs. (A). Heatmap reveals a full co-association pattern of all TFs. (B). The density of co-association score for all TF pairs. Determining co-association score as 0.1 (dashed in red) of significant TF co- associations based on an elbow statistic.

Fig. S5. Co-regulation of virulence-related master regulators. Graph diagram of interactions involving target genes of 4 TFs, including PA1380, PA0815, PA5428, and PA3973.

Fig. S6. Regulators involved in TCA and Ribosome pathways. (A). Radar plots show the putative regulator in TCA and ribosome pathways. (B). Graph diagram of interactions involving target genes of 14 TFs, including PA0611, PA5218, PA5431, PA1759, PA0403 (PyrR), PA 1283, PA4381, PA5255 (AlgQ), PA4784, PA0475, PA0893, PA0448, PA2957 and PA5438.

Fig. S7. Conservation comparison of PhoB and RpoN. (A). The coverage of PhoB and RpoN peak regions over PAO1 and PA14 chromosomes. Each line shows the location and log2 Fold Enrichment of peaks signal in the chromosome. (B). The binding motifs of PhoB and RpoN were analyzed via MEME-ChIP. All peaks were used to define the binding motif. (C). Comparison of overlapped targets enriched in PAO1 and PA14 of PhoB and RpoN. The Fisher test made the significance of the overlap ratio.

Table S1. List of 172 ChIPed TFs.

Table S2. Enriched GO terms of 172 ChIPed TFs.

Table S3. Hierarchical regulatory network.

Table S4. Ternary regulatory motifs.

Table S5. PWM pairwise similarity scores.

Table S6. Co-association network.

Table S7. Virulence-related master regulators.

Table S8. Strains and primers used in this study.

Table S9. List of strains for pan-genome analysis.