Advanced transcriptomic analysis reveals the role of efflux pumps and media composition in antibiotic responses of Pseudomonas aeruginosa

Abstract Pseudomonas aeruginosa is an opportunistic pathogen and major cause of hospital-acquired infections. The virulence of P. aeruginosa is largely determined by its transcriptional regulatory network (TRN). We used 411 transcription profiles of P. aeruginosa from diverse growth conditions to construct a quantitative TRN by identifying independently modulated sets of genes (called iModulons) and their condition-specific activity levels. The current study focused on the use of iModulons to analyze the biofilm production and antibiotic resistance of P. aeruginosa. Our analysis revealed: (i) 116 iModulons, 81 of which show strong association with known regulators; (ii) novel roles of regulators in modulating antibiotics efflux pumps; (iii) substrate-efflux pump associations; (iv) differential iModulon activity in response to beta-lactam antibiotics in bacteriological and physiological media; (v) differential activation of ‘Cell Division’ iModulon resulting from exposure to different beta-lactam antibiotics and (vi) a role of the PprB iModulon in the stress-induced transition from planktonic to biofilm lifestyle. In light of these results, the construction of an iModulon-based TRN provides a transcriptional regulatory basis for key aspects of P. aeruginosa infection, such as antibiotic stress responses and biofilm formation. Taken together, our results offer a novel mechanistic understanding of P. aeruginosa virulence.


INTRODUCTION
The transcriptional regulatory network (TRN) for any organism controls the expression of genes in response to genetic and environmental perturbations. For pathogenic bacteria, the TRN coordinates important functions for the infection of hosts, such as virulence mechanisms, antibiotic resistance and biofilm formation (1,2). Uncovering the TRN of pathogenic bacteria can be helpful in identifying novel drug targets and the underlying molecular mechanisms of inhibition. Pseudomonas aeruginosa is an important member of the ESKAPEE (Enterococcus faecium, Staphylococcus aureus, Klebsiella pneumoniae, Acinetobacter baumannii, Pseudomonas aeruginosa, Enterobacter spp. and Escherichia coli) group of pathogens that are notoriously multi-drug resistant (3). In addition to its intrinsic antibiotic resistance mechanisms (3), the large genome and genetic plasticity of P. aeruginosa allow it to grow in varied environments, but complicate comprehensive TRN mapping (4).
Previous efforts to characterize TRNs from Pseudomonas species are consolidated on the 'Pseudomonas Genome DB,' which annotated 371 confirmed or putative transcription factors (TFs) in the P. aeruginosa PAO1 genome (5). Efforts have been made to elucidate the TRNs mainly by the characterization of TFs using experimental approaches like transcriptional profiling, chromatin immunoprecipitation, and other similar methods (6,7). Further, Wang et al. validated the TFs in P. aeruginosa involved in virulence using the high-throughput systematic evolution of ligands by exponential enrichment assay (8). These studies describe the mechanistic details of the TF binding sites and regulon annotation. However, these methods are expensive, time-consuming and not derived from the transcriptional profiles. Although these studies are necessary for providing a foundation for understanding the TRN of P. aeruginosa, other approaches harness big data to generate a comprehensive understanding of the TRN. Therefore, in the current study, we use transcriptional profiles of P. aeruginosa to reconstruct the TRN.
Independent component analysis (ICA) is used to separate multivariate signals into additive components. For transcriptomic datasets, it identifies independently modulated sets of genes, called iModulons and quantifies their activity under specific conditions. Since iModulons capture the signals from transcriptional regulators, they can be compared to regulons. Regulons are sets of co-regulated genes defined based on bottom-up approaches using a variety of biomolecular methods, whereas iModulons are defined in a top-down manner using machine learning of entire transcriptomic profiles. Previously, we used ICA to annotate the TRNs of Escherichia coli (9), Staphylococcus aureus (10), Bacillus subtilis (11) and Sulfolobus acidocaldarius (12), which generated valuable hypotheses including putative regulatory interactions, novel associations between regulators and the conditions which may activate them, and specific insights into transcriptomic reallocation during key physiological processes. ICA has also been used to study the effect of adaptive laboratory evolution on the TRN (13,14). Results from all public iModulon analyses are available on our user-friendly online knowledgebase, iModulonDB (https://imodulondb.org/about.html) (15). Previously, we used 364 expression profiles of P. aeruginosa to generate and characterize 104 iModulons, which elucidated the genomic boundaries of biosynthetic gene clusters (BGCs), identified a novel bacteriocin producer, named ribosomally synthesized and post-translationally modified peptides (RiPP), that may play a role in virulence, and uncovered stimulons (clusters of iModulons with correlated activities) specific to iron and sulfur (16). Here, we expanded upon this TRN foundation by using new profiles generated from antibiotic treatment and biofilm growth conditions, since these are highly relevant for P. aeruginosa during infections of human hosts.
In this study, we used 411 high-quality expression profiles (281 publicly available in the SRA database and 130 generated for this study) to reconstruct the TRN of P. aeruginosa by implementing independent component analysis (ICA) using our established pipeline (17). We identified and characterized 116 iModulons comprising a quantitative TRN for P. aeruginosa based on all publicly available, high-quality RNA-sequencing profiles. Our analysis revealed novel insights into pathogenicity and antibiotic responses in this organism. We describe the regulators that modulate antibiotic efflux pumps, the differential activity of beta-lactam antibiotics in bacteriological cell culture media, such as LB and CAMHB, and physiological or eukaryotic cell culture media, and the influence of betalactam antibiotics on the Cell Division iModulon. To facilitate the broad utilization of this informative TRN, all iModulons can be searched, browsed, and studied through iModulonDB.org.

RNA extraction and library preparation
P. aeruginosa (PAO1 and ΔmexB) strains were used in this study. We extracted RNA samples for 45 unique conditions including different media types (M9 + glucose + amino acids (M9glucAA), cation-adjusted Mueller-Hinton Broth, Luria-Bertani (Lennox) Broth (LB), RPMI 1640 + 10% LB), oxidative stress (treatment with paraquat), iron starvation (treatment with DPD), osmotics stress (high NaCl), low pH, various carbon sources (succinate, glycerol, pyruvate, fructose, sucrose, N-acetyl glucosamine), micronutrients (copper, iron, zinc, sodium hypochlorite), antibiotics (gentamicin, azithromycin, piperacillin, ceftazidime, mecillinam, meropenem, tetracycline), and biofilms. For the reference condition M9glucAA, 18 amino acids, with the exception of cysteine and methionine, were added to a final concentration of 50 ug/ml for each amino acid and a glucose concentration of 0.2% (w/v). The concentration of NaCl in our standard LB media is 5 g into 1 l, or 0.5% (w/v) NaCl. All conditions were collected in biological duplicates and untreated controls were also collected for each set to rule out the possibility of the batch effect (Supplementary Table S1).
We followed the protocol outlined by Rajput et al. to perform RNA extraction and library preparation (16). In brief, strains were grown overnight at 37 • C, with rolling, in appropriate media types for the testing condition of choice. Overnight cultures were then diluted to a starting OD 600 of ∼0.01 and grown at 37 • C, with stirring. Once cultures reached the desired OD 600 of 0.4, 2 ml cultures were immediately added to centrifuge tubes containing 4 ml RNAprotect Bacteria Reagent (Qiagen), vortexed for 5 s, and incubated at room temperature for 5 min. Samples were then centrifuged for 10 min at 5000 × g and the supernatant was removed prior to storage at -80 • C until further processing. In conditions involving antibiotic treatment, when the bacterial culture had reached an OD 600 of ∼0.2, antibiotics were added at 2× or 5× their MIC in the appropriate media type and allowed to incubate at 37 • C, with stirring, for an additional hour prior to sample collection.
Total RNA was isolated and purified using a Zymo Research Quick-RNA Fungal/Bacterial Microprep Kit from frozen cell pellets previously harvested using Qiagen RNAprotect Bacteria Reagent according to the manufacturers' protocols. Ribosomal RNA was removed from 1 ug total RNA with the use of a thermostable RNase H (Hybridase) and short DNA oligos complementary to the ribosomal RNA, performed at 65 • C to prevent non-specific degradation of mRNA. The resulting rRNA-subtracted RNA was made into libraries with a KAPA RNA HyperPrep kit incorporating short Y-adapters and barcoded PCR primers. The libraries were quantified with a fluorescent assay (ds-DNA AccuGreen quantitation kit, Biotium) and checked for proper size distribution and average size with a TapeStation (D1000 Tape, Agilent). Library pools were then assembled and a 1× SPRI bead cleanup performed to remove traces of carryover PCR primers. The final library pool was quantified and run on an Illumina instrument (NextSeq, Novaseq).

Fluorescence microscopy
After the collection of RNAseq samples, fluorescence microscopy was performed as previously described (18). 4 l cells were added to a 4 l dye mix containing 100 g/ml FM4-64 (red), 100 g/ml DAPI (blue), and 20 M SYTOX Green (green) in 1× T-base. Samples were transferred to a single-well glass slide containing an agarose pad (1.2% agarose in 20% LB) and imaged on an Applied Precision DV Elite epifluorescence microscope. The exposure times for each wavelength were kept constant for all images.

Data preprocessing
Apart from the in-house generated data, we also downloaded and processed all RNA sequencing data available from NCBI SRA for P. aeruginosa PAO1. In the current study, we used 411 expression profiles, of which 281 is collected from the NCBI-SRA and 130 is generated in the lab. Out of 130 lab-generated expression profiles, 83 include the conditions like osmotic stress, carbon sources (glucose, succinate, pyruvate, glycerol, N-acteyglucosamine), metals (Zn, Cu, Fe), media types (CAMHB, M9 and R10LB) and salt stress, which is used in our previous study. However, in the current study, we added some conditions like antibiotics and biofilms to reconstruct the TRNs. We include antibiotics like gentamicin, azithromycin, piperacillin, ceftazidime, mecillinam, meropenem and tetracyclines in media like CAMHB, R10LB, and M9.
To ensure quality control, data that failed any of the following four FASTQC metrics were discarded: per base sequence quality, per sequence quality scores, per base n content, and adapter content. Samples that contained under 500 000 reads mapped to coding sequences were also discarded. Hierarchical clustering was used to identify samples that did not conform to a typical expression profile.
Manual metadata curation was performed on the data that passed the first four quality control steps. Information including the strain-description, base media, carbon source, treatments, and the temperature was pulled from the literature. Each project was assigned a short unique name, and each condition within a project was also assigned a unique name to identify biological and technical replicates. After curation, samples were discarded if (a) metadata was not available, (b) samples did not have replicates or (c) the Pearson R correlation between replicates was below 0.95. Finally, the log-TPM data within each project was centered on a project-specific reference condition. After quality control, the final compendium contained 411 high-quality expression profiles: 130 generated for this study, plus 281 expression profiles extracted from public databases.

Computing robust components with ICA
To compute the optimal independent components, an extension of ICA was performed on the RNA-seq dataset as described in McConn et al. (23).
Briefly, the scikit-learn (v0.23.2) (24) implementation of FastICA (25) was executed 100 times with random seeds and a convergence tolerance of 10 -7 . The resulting independent components (ICs) were clustered using DBSCAN (26) to identify robust ICs, using an epsilon of 0.1 and a minimum cluster seed size of 50. To account for identical components with opposite signs, the following distance metric was used for computing the distance matrix: where ρ x,y is the Pearson correlation between components x and y. The final robust ICs were defined as the centroids of the cluster.
Since the number of dimensions selected in ICA can alter the results, we applied the above procedure to the dataset multiple times, ranging the number of dimensions from 10 to 360 (i.e. the approximate size of the dataset) with a step size of 10. To identify the optimal dimensionality, we compared the number of ICs with single genes to the number of ICs that were correlated (Pearson R > 0.7) with the ICs in the largest dimension (called 'final components'). We selected the number of dimensions where the number of nonsingle gene ICs was equal to the number of final components in that dimension.

Determination of the gene coefficient threshold
The gene coefficients are determined as described in Sastry et al. (17). Each independent component contains the contributions of each gene to the statistically independent source of variation. Most of these values are near zero for a given component. In order to identify the most significant genes in each component, we iteratively removed genes with the largest absolute value and computed the D'Agostino K 2 test statistic (27) for the resulting distribution. Once the test statistic dropped below a cutoff, we designated the removed genes as significant.
To identify this cutoff, we performed a sensitivity analysis on the concordance between significant genes in each component and known regulons. First, we isolated the 20 genes from each component with the highest absolute gene coefficients. We then compared each gene set against all known regulons using the two-sided Fisher's exact test (FDR < 10 -5 ). For each component with at least one significant enrichment, we selected the regulator with the lowest P-value.
Next, we varied the D'Agostino K 2 test statistic from 50 through 2000 in increments of 50, and computed the F1-score (harmonic average between precision and recall) between each component and its linked regulator. The maximum value of the average F1-score across the components with linked regulators occurred at a test statistic of cutoff of 420 for the P. aeruginosa dataset.

iModulon annotation
The regulator enrichments are determined as described in Sastry et al. (17). The gene annotation pipeline can be found at https://github.com/SBRG/pymodulon/blob/master/docs/ tutorials/creating the gene table.ipynb. Gene annotations were pulled from Pseudomonas genome db (5). Additionally, KEGG (28) and Cluster of Orthologous Groups (COG) information were obtained using the EggNOG mapper (29). Uniprot IDs were obtained using the Uniprot ID mapper (30), and operon information was obtained from Biocyc (31). Gene ontology (GO) annotations were obtained from AmiGO2 (32). The known TRN was obtained from RegPrecise (33) and manually curated from the literature.

Differential activation analysis
The gene coefficients are determined as described in Rajput et al. (16). We calculated the difference in the iModulons activities between two or more conditions by fitting the log-normal distribution to the difference. The statistical significance was checked by calculating the difference between the absolute values of the activities among the iModulons, and further confirmed by the P-values. However, the p-values were adjusted using the Benjamini-Hochberg correction. While calculating the differential activation among the iModulons, the difference with >5 was considered significant.

Prediction of two-components system
Through our iModulon analysis, we identified the TCSs regulating the efflux pumps. It is done in three steps: (i) Identification of the TCSs by using our previously published method (34) and P2CS database (35). (ii) Checking the relationship between the RR gene and the iModulon containing efflux pumps by calculating the PCC between them. (iii) Validating the RR's specificity/relationship by scanning its correlation with all predicted 116 iModulons. Furthermore, we also checked the specificity of the metabolites with respective efflux pumps.

Generating an iModulonDB web page
The P. aeruginosa aeruPRECISE411 iModulonDB (15) web pages were generated using the imodulondb export() function in the Pymodulon package (17). The generated page includes gene information from Pseudomonas Genome DB (https://pseudomonas.com) (5).

Overview of Pseudomonas aeruginosa iModulons
In our previous study, we used 364 transcriptional profiles (281 from literature + 83 lab-generated) of P. aeruginosa, which is non-informative about the antibiotics responses (16). The aeruPRECISE364 profiles included the data from conditions like media types with different nutrient sources, osmotic stress and specific gene-knockouts to identify the TRN associated with biosynthetic gene clusters, and secretion systems, and central metabolism (16). As the antibioticspecific responses are lacking in the previous study, we generated and compiled the antibiotic conditions in the current study.
The present study expands our transcriptome compendium from 364 to 411 high-quality profiles, introducing 47 new conditions including different antibiotic treatments and biofilm growth (Supplementary Table S1, Figure 1A). All transcriptomic profiles were processed using our established pipeline, (17) all samples passed quality control, and all the replicates among the 411 profiles show Pearson's correlation coefficient (PCC) of 0.98. To remove batch effects, all samples were normalized to respective reference conditions (17). We applied our ICA workflow on the expanded compendium to compute an updated TRN composed of 116 independently modulated sets of genes, called iModulons (Supplementary Figure S1b). The new decomposition explains 66.16% of the total variance in the compendium of 411 profiles. Each iModulon was annotated with one of 11 functional categories: antibiotic resistance, amino acid and nucleotide metabolism, biosynthetic gene clusters, carbon source utilization, metal homeostasis, prophages, quorum sensing, stress response, structural components, translational and uncharacterized.
The transcriptional regulators of each iModulon were identified by comparing them with the known regulons. Regulon data was collected from RegPrecise (33) and manually mined from literature (Supplementary Table S1). We used 'regulon recall' and 'iModulon recall' to characterize the association of the 116 iModulons with the known regulons. Regulon recall is the number of shared genes between an iModulon and its associated regulon divided by the total number of genes in the regulon, while the iModulon recall is the number of shared genes divided by the total number of genes in the iModulon. Based on the regulon recall and the iModulon recall and a threshold of 0.6, 116 iModulons are divided into 4 categories: well-matched, regulon subset, regulon discovery, and poorly matched (Supplementary Figure  S1b and Supplementary Notes 1).
We compared the 116 iModulons identified in this study with our previous work (16) which resulted in 104 iModulons ( Figure 1B and Supplementary Table S2). By correlating the gene weights of the iModulons between the two datasets, we could map the previously characterized signals to the new ones. The large overlap between recovered signals indicates that the TRN structure is robust to new data, and the new iModulons show transcriptional response triggered by the new conditions used.
We then sought to use our iModulon structure to predict possible regulators for those systems. Given that regulators very often exhibit feedback, it is common for TFs to be in the iModulons that they also regulate (45)(46)(47). Therefore, any predicted regulators in the efflux-pump containing iModulons should be explored in future studies for their putative role in regulating the pumps and possible associated effects on antibiotic resistance (Supplementary Figures   S2 and S3). Upon searching the literature for these associations, we found that six of the eight interactions were known in past studies, which supports that this approach can uncover meaningful relationships. The final two possible interactions are therefore important as potential new discoveries. Specifically, PA2274 may regulate mexGHI-opmD and PA4878 (brlR) may regulate mexPQ-opmE.
It is also important to know which conditions may activate the TFs and their associated iModulons. We therefore list activating conditions in Table 1 to facilitate further study. We also used the activity levels of the iModulons under various conditions to identify metabolites that may activate the efflux pump containing iModulons, some of which were not known in the literature: DADS (EsrC), metals and organic compound (PprB), and n-acetyl-glucosamine (Glc-NAc) (MexS-1) ( Figure 2B and Table 1). Additional details on each efflux pump can be found in Supplementary Notes 2.
Furthermore, we compared the activities of the identified iModulons containing efflux pumps under different treatments or growth conditions ( Figure 2B). As expected, there was high activity in the CzcR and CueR iModulons in response to copper (CuSO 4 ) supplementation (48). Of interesting note is the high level of activation of the MexS-1 iModulon during growth on GlcNAc, during trace metal supplementation, and during NaOCl treatment (Figure 2B). A few studies suggest the increased concentration of GlcNAc, and metals like Cu, Zn, Fe, etc. are indicators of increased pathogenicity in the host in cystic fibrosis (CF) patients (49,50). Further, the mexEF-oprN efflux pump is known to increase virulence in various conditions including CF (51,52). Our study suggests that the mexEF-oprN efflux pump is specific for GlcNAc, Cu, Fe, and Zn. Thus, we  hypothesize that mexEF-oprN might play a role in increasing virulence in CF conditions. Two-component systems (TCSs) often play a significant role in regulating efflux pumps (53). Recent work from our group (34) and others (47,54) has characterized them and their targets in P. aeruginosa. We, therefore, compared our putative new regulators to the set of known and predicted TCSs. PA2274 may be a response regulator with an associated histidine kinase of PA2275 (Supplementary Notes 2).
Therefore, ICA analysis of the P. aeruginosa transcriptome was helpful in the identification of two regulators (PA4878 (brlR) and PA2274) which potentially regulate the expression of efflux pumps (mexPQ-opmE and mexGHI-opmD, respectively). We also identified substrates that specifically activate each efflux pump (as detailed in the form of a heatmap Figure 2B, Table 1). Further studies guided by these suggested relationships could elucidate the mechanisms underlying cellular transport and ultimately be useful for the design of antimicrobial drugs.

Differential activity of beta-lactam antibiotics between bacteriological and physiological media
Beta-lactam antibiotics are often used clinically against P. aeruginosa infections, typically in conjunction with betalactamase inhibitors, to combat mechanisms of resistance in the pathogen. Recent studies have demonstrated that bacterial sensitivity to antibiotics can be altered depending on the growth medium used (55,56). We determined the minimal inhibitory concentration (MIC) of a variety of antibiotics against an efflux-pump knockout strain of P. aeruginosa ( mexB) in different media types, including the classic susceptibility testing media, cation-adjusted Mueller-Hintron Broth (CAMHB), and the eukaryotic cell culture media, RPMI 1640 supplemented with 10% Luria Broth (R10LB). We found that three of the four beta-lactam antibiotics tested had higher MICs in R10LB compared to CAMHB, whereas two of the three protein synthesis inhibitors had lower MICs in R10LB relative to CAMHB ( Table 2).
Because of the observed differences in MIC, we surveyed the identified iModulons to determine if there were media-specific differences in the transcriptional response to the four beta-lactam antibiotics. We found that activities for certain iModulons containing antibiotic resistancerelated genes, such as efflux pumps, had a differential activity that could potentially be linked to media type (Figure 3A, B). The iModulon SoxR, for example, contains the genes mexGHI and ompD ( Figure 3C) and has upregulated activity during treatment with beta-lactam antibiotics in R10LB and downregulated activity in CAMHB ( Figure  3A), which correlates with the observed differences in MIC. The difference in the MIC value of antibiotics in both media types were also found in our previous studies on Staphylococcus aureus (57)(58)(59). Similarly, there was increased activity of the iModulon CzcR/Fur, involved in divalent cation homeostasis (60), in R10LB compared to CAMHB ( Figure  3A). CzcR/Fur has been shown to mediate antibiotic susceptibility by regulating the expression of OprD, which is a route of entry for carbapenems like meropenem (60). Interestingly, there seemed to be no difference in activities for iModulons CpxR, CueR and MexZ, while activities for the iModulon AmpC were upregulated compared to untreated controls but highly variable ( Figure 3A).
When we compared all iModulon activities of individual beta-lactam treatments between CAMHB and R10LB, we found some media-specific differences. The iModulon TagR1 is upregulated in three of the four beta-lactam treat- ments in R10LB compared to CAMHB ( Figure 3D). The TagR1 iModulon contains genes related to the type VI secretion system, specifically the PpkA protein kinase, which has been shown to mediate bacterial adaptability as well as oxidative stress tolerance and pyocyanin synthesis (61,62). This suggests that TagR1-mediated oxidative stress tolerance may help P. aeruginosa tolerate these antibiotics in the physiological milieu but not in CAMHB, where antibiotic susceptibility is typically tested. Additionally, the iModulon PchR, involved in pyochelin synthesis and iron acquisition (63), has high activity in R10LB conditions compared to the CAMHB samples ( Figure 3D). In comparison, the iModulons Un-9, which encompasses genes in energy metabolism, and Anr-2, which contains genes related to iron acquisition, were both highly activated in the CAMHB conditions relative to the R10LB samples ( Figure 3D). In combination, this data demonstrates that although the antibiotics have the same molecular targets within P. aeruginosa, the transcriptional response of the bacteria to the treatment can vary depending on media type.

Differential activation of the cell division iModulon by betalactam antibiotics
Beta-lactam antibiotics have different affinities for various penicillin-binding proteins and enzymes responsible for coordinating cell wall synthesis and bacterial cell division (64)(65)(66). It has been shown that treatment of Eshcerichia coli with different antibiotics, including cell-wall active compounds, leads to morphologically distinct cytology based on molecular mechanisms of action (18). Mecillinam binds specifically to PBP-2, which is essential for the maintenance of rod-morphology in P. aeruginosa, E. coli, and other Gram-negative bacteria (67). Meropenem targets PBP-2 and PBP-4 (65,68), a non-essential DDendopeptidase (69,70). Piperacillin and ceftazidime have specificity for PBP-3 (65), also known as FtsI, which is a member of the divisome and is important for the septation of daughter cells (71,72). P. aeruginosa treated with these antibiotics have morphological differences based on the PBP of the antibiotic targets ( Figure 4C). Mecillinamtreated cells are smaller and more oblong-shaped, and penicillin or ceftazidime treatment leads to cell filamentation ( Figure 4C). We were interested to see if the introduction of stress on the same biosynthetic pathway, i.e. cell wall, by different molecular means, in this case, beta-lactam antibiotics, would have the same or different transcriptional response in P. aeruginosa.
We identified a single cell division iModulon ( Figure  4A). The cell division iModulon contained the genes murC, murE, and murG, which are important for the synthesis of peptidoglycan (73)(74)(75), as well as divisome components FtsA and FtsQ (76,77) ( Figure 4A). This iModulon also contains PBP3, also known as FtsI ( Figure 4A), which is the specific target for penicillin and ceftazidime (65). The cell division iModulon activity increased in R10LB compared to untreated but not when grown in CAMHB, during beta-lactam treatment, with the exception of mecillinam ( Figure 4B). The untreated samples have similar levels of activity for this cell division iModulon ( Figure 4B). Mecillinam treatment had little effect on the activity in this iModulon and also had a smaller effect on cell morphology. The higher levels of the cell division iModulon activity in the treated samples grown in R10LB may be a major contributor to the observed higher MICs in this media compared to CAMHB. The cells grown in R10LB produce more divisome-related transcripts that likely allow for greater tolerance of cell wall-active antibiotic-induced stress. The difference in MIC might also be due to the activation of efflux pumps, as noted earlier. We also discovered that there was high Cell Division iModulon activity in samples taken from planktonic growth (PRJNA643216) as well as in those grown in high salt conditions (Figure 4b), but low activity in samples taken from P. aeruginosa grown in a biofilm (PR-JNA643216). This finding is confirmatory, since we would expect actively dividing cells, such as those growing planktonically in a bioreactor (78), to have high cell division iModulon activity compared to slower-growing cells, such as those in a biofilm. Identifying the regulatory mechanisms underlying this iModulon's activation should be the topic of further study, as additional drugs which downregulate this response could be administered in conjunction with betalactams to improve their effectiveness.

The PprB iModulon activates under stress conditions and is a signal for the transition from planktonic to biofilm lifestyle
The PprB iModulon consists of the PA1874-PA1877 efflux pump, tight adherence (tad) machinery, and the cupE1-cupE5 (flp fimbriae) chaperone-usher pathway ( Figure 5A-C). We found that the PprB iModulon is highly activated during biofilm growth in M9glucAA and CAMHB conditions as well as during growth with iron (Fe, 5uM), zinc (Zn, 5uM), and n-alkane (equal volumes of the C 8 , C 10 , C 12 , C 14 , and C 16 ) supplementation ( Figure 5D). It has been previously established that PprB is involved in biofilm formation (79,80), including the regulation of tad machinery expression involved in the production of the Type IVb pili (81,82). Additionally, PprB was demonstrated to control fimbriae assembly via CupE, which is specifically expressed during biofilm growth (44). Although the molecules excreted by the efflux pump are unknown, the tad system and cupE1-E5 system both support biofilm formation. The activation of this set of genes under high-Fe conditions supports recent evidence that high-Fe conditions lead to a limited motility phenotype that promotes biofilm formation (83).
The genes PA1874 through PA1877 were included in the PprB iModulon. This operon was previously characterized as a likely type I secretion system (T1SS), important for resistance to aminoglycosides and fluoroquinolones, specifically during growth in biofilms (80,84). Further, these studies observed the increased susceptibility to the aminoglycoside tobramycin during constitutive expression of pprB in biofilm growth, which was thought to be through pprB modulation of membrane permeability (85). We found that the PprB iModulon is activated during planktonic growth with treatment with antibiotics that target protein synthesis, such as gentamicin, azithromycin, and tetracycline, as well as NaOCl and paraquat (pqt), which are both redoxactive compounds (Figures 2C and 5D). This suggests that the PprB iModulon is activated under stress conditions and may be a contributing signal for P. aeruginosa to switch between planktonic and biofilm modes of growth.
We found a significant correlation between the activity levels of the PprB iModulon and the MexS-1 and RpoS-1 iModulons. We found that the MexS-1 iModulon had a strong bimodal activity distribution, and all samples which highly activated the MexS-1 iModulon also activated the PprB iModulon. These samples included the presence of supplemental Fe and Zn, during growth on GlcNAc, as well as during the biofilm mode of growth ( Figure 5E, F). We also found a linear correlation between the activities of the PprB and the RpoS-2 iModulons, the latter of which contains the central stress response sigma factor RpoS. Both iModulons activated during growth in biofilm or osmotic stress (e.g. NaCl), and were highly inactivated during treatment with diallyl disulfide (DADS) ( Figure 5F) Treatment of P. aeruginosa with DADS was shown to down-regulate many genes involved in pathogenicity and virulence, including TCSs, such as PprB (86). These correlations indicate either co-stimulation or hierarchical regulation; in the case of RpoS-2, it suggests that stress correlates with biofilm formation, which is consistent with expectations about the function of biofilm in protecting P. aeruginosa communities. Stewart et al. also show that stress can contribute to antibiotic resistance in biofilms of P. aeruginosa.

DISCUSSION
P. aeruginosa is an important member of the ESKAPEE group of pathogenic bacteria (87). It is an opportunistic pathogen and shows high resistance to diverse classes of antibiotics (88). In the current study, we used ICA on 411 expression profiles of P. aeruginosa to generate 116 iModulons. Our analysis reveals various important aspects of P. aeruginosa transcriptional regulation, such as the prediction of regulators modulating the efflux pumps, specific macromolecules and antibiotics which activate transcription of efflux pumps, differential behavior of beta-lactam antibiotics in bacteriological and physiological media, the impact of beta-lactam antibiotics on the PBPs and cell division, and the role of the PprB iModulon in the transition to biofilm growth under stress conditions.
The expression of efflux pumps is one of the classic mechanisms of antibiotic resistance (89). Apart from antibiotics, efflux pumps transport various compounds like metals and organic compounds out of the cell (90). In this study, we predicted the regulators which may regulate efflux pumps based on correlated transcription (Fig-ure 2 and Table 1). Among the eight iModulons related to efflux pumps captured by our study, we found two potentially novel regulators regulating efflux pumps: PA4878 (brlR) and PA2274 are predicted to regulate mexPQ-opmE and mexGHI-opmD, respectively (Table 1). Further, we confirmed the role of identified regulators in regulating the specific efflux pumps by two step analysis. First, we checked the correlation among the identified response regulator and the efflux pump containing iModulons. Next, we used the activating conditions of the efflux pumps to predict their substrates, which identified new putative relationships.
We found that antibiotics sharing the same molecular target can activate iModulons for efflux pumps or cell division to different levels, and their activities are also dependent on the media type (Figures 3 and 4). Therefore, this work provides a molecular view into transcriptional regulation underpinning the media-specific susceptibility differences observed via MIC between P. aeruginosa treated in the bacteriological media CAMHB and the LB-supplemented eukaryotic cell culture media R10LB (Table 2).
We also identified a single cell division iModulon of P. aeruginosa that had differential activity dependent on media condition, e.g. CAMHB versus R10LB, and mode of growth, e.g. biofilm versus planktonic. Confirming experimental observations of higher MICs during growth in R10LB, we found increased activity of the cell division iModulon during beta-lactam antibiotic treatment in R10LB relative to untreated and to CAMHB conditions. This highlights the ability of iModulons to identify potentially clinically-relevant transcriptional responses from in vitro experimental data sets.
The current study is focused on elucidating the antibiotic-specific behavior of P. aeruginosa through the application of machine learning to a compendium of transcriptomic profiles. The code used to explore the iModulons is available on GitHub (https://github.com/akanksha-r/modulome paeru2.0). To make our data easily accessible to experimental biologists we provide our data in the form of a user-friendly web interface (https://imodulondb.org/dataset.html?organism= p aeruginosa&dataset=precise411). We found that the iModulons are an important tool for the reconstruction of the TRN of P. aeruginosa to address various biological queries, including predicting the regulators and substrates of efflux pumps, examining the differential behavior of beta-lactam antibiotics in different media, exploring the role of biofilm-specific antibiotic resistance systems in the virulence of P. aeruginosa, and elucidating the molecular mechanism of beta-lactam antibiotics on PBPs and cell division. This study demonstrates that iModulon-based TRNs can provide new insights when new transcriptomic data is added; this work can therefore serve as the basis for further studies. With increasing scale, we hope that this approach will culminate in a comprehensive, quantitative, irreducible TRN for this important pathogen.

DATA AVAILABILITY
All the in-house generated sequences were deposited in the NCBI-Sequence Read Archive database (PRJNA717794). The accession number of the deposited reads is provided in Supplementary Table S1. While the X, M and A matrices are available at GitHub (https://github.com/akankshar/modulome paeru2.0). Each gene and iModulon have interactive, searchable dashboards on iModulonDB.org, and data can also be downloaded from there.

CODE AVAILABILITY
The customized code for the ICA analysis, as well as various files including the X, M, A matrices, TRN regulator file, gene, annotated files, gene ontology, and Kegg pathway annotation files, are available on GitHub (https://github.com/ akanksha-r/modulome paeru2.0).