Strains used in this study are described in Table 1. Pst DC3000 strains were grown at 28°C in King’s B (KB) medium (King et al., 1954) and on KB medium solidified with 1.5% (w/v) agar. Escherichia coli TOP10 (Invitrogen) strains were grown at 37°C in LB medium and on LB medium solidified with 1.5% (w/v) agar (Bertani, 1951). Kanamycin was supplied at 50 μg/mL and gentamycin was supplied at 10 μg/mL. Plasmid DNA was isolated using Qiagen Miniprep Kit (Qiagen) from overnight cultured E. coli TOP10 cells and subsequently used to transform Pst DC3000 and mutant derivatives by electroporation (Choi et al., 2006).

The Pst DC3000 random barcoded-transposon insertion library was constructed by conjugating the E. coli APA752 donor library containing the barcoded mariner plasmid pKMW3 into Pst DC3000. Mating was performed overnight at a 1 (E. coli) to 4 (Pst DC3000) ratio on nutrient agar (2 g l⁻¹ yeast extract, 5 g l⁻¹ peptone, 5 g l⁻¹ NaCl containing 300 μM diaminopimelic acid). The bacterial mixture was harvested into a single sample, resuspended in liquid KB, plated on KB supplemented with 1.5% (w/v) agar, 50 μg ml⁻¹ kanamycin and 50 μg ml⁻¹ rifampicin plates and incubated for 2 d at 28°C to select for insertion mutants. Approximately 214,000 mutant clones were obtained, each one representing one transposition event. These clones were pooled together in 200 mL of KB media with 50 μg ml⁻¹ kanamycin and 50 μg ml⁻¹ rifampicin and 20% glycerol, aliquoted in 1 mL samples and frozen at −80°C.

Genomic DNA from the Pst DC3000 library was prepared for sequencing as previously described (Wetmore et al., 2015; Price et al., 2018). In brief, genomic DNA was fragmented by ultrasonication to an average size of 300 bp and an Illumina library preparation was performed based on modifications of the NEBNext DNA library preparation kit (New England Biolabs) with PCR enrichment of the transposon insertion sites. Mapping and pool construction were performed using the FEBA scripts MapTnSeq.pl. and DesignRandomPool.pl.,¹ requiring a minimum of 10 reads per barcode to be included in the mapping (Wetmore et al., 2015). Based on the TnSeq mapping data, we used feba/bin/Essentiality.pl. as in Price et al. (2018) to calculate the read density (reads/nucleotides across the entire gene) as well as the insertion density (sites/nucleotides) within the central 10–90% of each gene. We then used the “Essentials” function in feba/bin/comb.R (Price et al., 2018), to predict likely essential genes. This analysis estimates how short a gene could be and still be unlikely to have no insertions by chance (p < 0.02, Poisson distribution). After excluding genes shorter than 350 bp, the read density was normalized by GC content, and the insertion density was normalized so the median gene’s value was 1. Genes were considered essential or nearly essential if both the normalized read density and the normalized insertion density were < 0.2. Because mariner transposons require a TA dinucleotide site for insertion, genes with 0 TA sites were also excluded.

One aliquot (1 mL) of the library glycerol stock was thawed and resuspended to an OD₆₀₀ of approximately 0.3 in 20 mL of KB media supplemented with 50 ug/mL rifampicin and 50 ug/mL kanamycin. The culture was incubated at 30°C until the OD₆₀₀ had doubled, at which time an aliquot was harvested to serve as the time 0 (T0) sample. The remaining cells were pelleted, washed twice with M9 minimal media, and resuspended to a final OD₆₀₀ of 0.02 in 5 mL of M9 minimal media ± 2 mg/mL Casamino acids. Three replicates of each culture condition were incubated at 30C until they achieved an OD₆₀₀ of approximately 1 (close to six doublings), then harvested for BarSeq.

MOPS minimal media supplied with 0.2% fructose and 0.3% agar was used for all swimming assays (Neidhardt et al., 1974). Five milliliter of swimming assay medium was added to 35 × 10 mm petri dish and air dried in a biosafety cabinet hood for exactly 30 min before inoculation. Pst DC3000 strains overnight cultures were resuspended and adjusted to an OD₆₀₀ of 0.4 in MOPS liquid media with 0.2% fructose prior to inoculation. One microliter inoculum from Pst DC3000 wild type or derivative strain culture was injected into the agar at the geometrical center of the petri dish. Plates were then transferred to a sealed moisture-enhanced chamber (RH = 100%) and incubated at 25°C for 48 h. For the batch swimming assay, 50 mL of swimming assay medium was added to square plates (86 × 128 mm) instead.

Initial inoculum generated from the RB-TnSeq library was prepared similarly as described above but from a recovery culture where a library aliquot was thawed on ice and added to 50 mL of LB with 50 μg/mL kanamycin. This recovery culture was then incubated at 28°C for 6 h before inoculum preparation. For liquid culture passages, a 1 μL inoculum from the library recovery culture was added to 5 mL of MOPS minimal media with 0.2% fructose and grown for 48 h at 25°C with shaking. At 48 hpi, this culture was spun down and cells were extracted for making new inoculum. Motility assays in both motile and non-motile selection series were set up as described above. For the motile selection series, the whole swimming colony was extracted from the passaging swimming plates with a 6 mm radius biopsy punch (Robins Instruments) at 48 hpi while a sample of cells was extracted from the inoculation point by a 200 μL micropipette for the non-motile selection series (Figure 1A). Bacterial cells were then released from the soft agar by physical homogenizing and resuspended into novel inoculum at OD₆₀₀ of 0.4 in MOPS liquid media with 0.2% fructose.

For the library pre-culture, BarSeq PCR, Illumina sequencing and fitness calculations, the exact protocol was followed as previously described (Helmann et al., 2019). Genomic DNA was purified using Wizard SV Genomic DNA Purification System (Promega), and the 98°C BarSeq PCR as described by Wetmore and colleagues was used to amplify the barcode region of each sample (Wetmore et al., 2015). The PCR for each sample was performed in 50 μL total volume: containing 0.5 μL Q5 High-Fidelity DNA polymerase (New England Biolabs, United States), 10 μL 5X Q5 buffer, 10 μL 5X GC enhancer, 1 μL 10 mM dNTPs, 150–200 ng template gDNA, 2.5 μL common reverse primer (BarSeq_P1), and 2.5 μL of forward primer from one of the 96 indexed forward primers (BarSeq_P2_ITXXX), both at 10 μM (Wetmore et al., 2015). Following the BarSeq PCR, 10 μL of each reaction was pooled, and 200 μL of this DNA pool was subsampled and purified using the DNA Clean and Concentrator Kit (Zymo Research, United States). The final DNA sequencing library was eluted in 30 μL nuclease-free water, quantified on a Nanodrop One (ThermoFisher, United States), and submitted for sequencing at the BRC Genomics Facility at the Cornell Institute of Biotechnology. Prior to sequencing, the quality of each amplicon pool was assessed using a Bioanalyzer. Each sequencing pool was run on a single NextSeq 500 (Illumina, Inc., United States) lane for 75 bp single-end reads.

Sequencing reads were used to calculate genome-wide gene fitness using the FEBA scripts MultiCodes.pl., combineBarSeq.pl., and BarSeqR.pl. (Wetmore et al., 2015). From the raw fastq reads for each sample, individual barcode sequences were identified and counted using MultiCodes.pl., and these counts were combined across samples using combineBarSeq.pl. Using BarSeqR.pl., fitness values for each insertion strain were calculated as the log₂ ratio of barcode abundance following library growth under a given experimental condition divided by the abundance in the time0 sample. The fitness of each gene is the weighted average of strain fitness values based on the “central” transposon insertions only, i.e., those within the central 10–90% of a gene. Barcode counts were summed between replicate time0 samples. For analysis, genes were required to have at least 30 reads per gene in the time0 sample, and 3 reads per individual strain, and gene fitness values were normalized across the chromosome so that the median gene fitness value was 0 (Wetmore et al., 2015).

Both fitness scores and an absolute t-like statistic were used to prioritize genes of interest. This t-value was calculated as the fitness score divided by the square root of the maximum variance calculated in two ways (Wetmore et al., 2015). High fitness value genes were selected by criteria fitness score >2 or <−2 and absolute t-score > 5. Data analysis was conducted in R v.4.2.3 and data visualization were carried out with ggplot2. Principal component analysis was performed on the gene fitness matrix across all 14 passages excluding time0 with R package prcomp.

Standard marker exchange mutagenesis was used to sequentially produce a Pst DC3000 strain with each gene of interest deleted (Kvitko and Collmer, 2011). Deletion constructs for PSPTO_1042 (pIY32) and PSPTO_1980 (pIY33) were each made using synthetic linear DNA fragments from Twist Bioscience that contained between 800 to 1,000 bp flanking both ends of each gene joined together with EcoRI and HindIII sites added to the 5′ and 3’ends, respectively (Supplementary Table S5). Each deletion construct retained the first and last six codons of the original gene to minimize polar effects on downstream genes. The DNA fragments were incorporated into pK18MobSacB by restriction enzyme digestion and ligation with T4-ligase (Thermo Fisher Scientific). Deletion construct for PSPTO_0406 pZB72 was created by PCR amplification of 0.7 and 0.8-kb flanks to PSPTO_0406 with oSWC4709/4711 and oSWC4712/4713, respectively. Then, gel purified PCR fragments were joined by a second PCR amplification with primers oSWC4710/4714 containing XbaI restriction sites. The approximately 1.4 kb product was gel purified, digested with XbaI and ligated with pK18mobsacB cut with the same restriction enzyme. The deletion construct for PSPTO_4229 (pIY26) was created in a similar fashion. Deletion constructs were confirmed by restriction analysis and sequencing before transformed into wild-type Pst DC3000 by electroporation (Choi et al., 2006) and selection for kanamycin resistant merodiploids. Sucrose was used to select for recombinants that had subsequently eliminated the sacB-containing deletion construct plasmid backbone and confirmed to have lost kanamycin resistance. The deletion was confirmed by sequencing PCR products amplified with primers that anneal to sequences flanking the deleted locus.

Overexpression plasmids were constructed with Gateway cloning. An expression cassette of each target gene was PCR amplified with CACC motif on 5′ end and incorporated in pENTR/D/SD backbones using TOPO cloning (Thermo Fisher Scientific). Then these DNA fragments were incorporated into pBS46 (Swingle et al., 2008) with Gateway LR Clonase II enzyme mix (Thermo Fisher Scientific). After whole plasmid sequencing (Plasmidsaurus) confirmation, these plasmids were then used to transform corresponding Pst DC3000 deletion mutant strains by electroporation (Choi et al., 2006). The transformants were then selected by Gentamycin resistance.

To screen for Pst DC3000 genes important in swimming motility, we constructed a library of random barcoded insertion mutants using a barcoded mariner transposon (Wetmore et al., 2015). High-throughput sequencing was used to map insertion locations in the genome as well as their corresponding 20 nucleotide barcode sequences. We identified 125,581 barcodes that were seen 10 or more times and mapped to the genome at 73,967 different locations. Of these, 73,138 insertions (at 43,173 distinct locations) were located in the central 10–90% of coding sequences and could therefore be used for fitness analysis (4,691 of 5,619 protein-coding genes). Overall, this mutant library consists of a median of 11 insertion strains per protein-coding gene. We also used these mapping data to predict 403 (8%) putative essential genes on the chromosome, based on protein-coding genes containing at least one TA dinucleotide site and at least 350 bp in length, while having no or very low abundance of transposon insertions (Supplementary Table S1). Enriched functional category annotations among these genes include translation, gene expression, cell wall biogenesis, lipid biogenesis and cofactor metabolic process (Supplementary Table S1).

Swimming in soft agar media is generally adaptive because it helps cells access nutrients and reduce competition with other clonemates. We predict that when a colony develops from a library of pooled mutants, motile cells will move from the inoculation point toward the periphery and leave the less motile mutants enriched in the center of the colony. Sampling the diversity in an entire colony after it forms should show a bias for fully motile mutants and a depletion of non-motile mutants. This bias should be reversed in the center of the colony near the inoculation point, where there should be an enrichment for non-motile and less-motile mutants. In contrast, growth in liquid with shaking moves the bacteria by external forces and should reduce or eliminate the adaptive benefits that flagella provide, while still maintaining selection for genes needed for growth.

We passaged our barcoded transposon insertion mutant library in swimming medium and sampled the entire colony or only the center at the inoculation point to select for the motile and non-motile populations, respectively. As a baseline control, we grew the library in liquid media to identify genes needed for growth and survival in the medium. In all cases we used MOPs minimal medium with 0.2% fructose (m/v) that was either prepared as liquid or semi-solidified with 0.3% agar (m/v) (Neidhardt et al., 1974). All three selection series were started from the same culture of the barcoded transposon insertion mutant library (Time 0 sample) and replicated in triplicate. For the liquid selection series, we inoculated MOPs liquid medium and incubated for 48 h at 25°C with shaking before sampling for DNA extraction and to inoculate a second passage. For the motile and non-motile selection series, the inoculum was injected into the medium and the colony was allowed to grow for 48 h at 25°C in a humid chamber. We then collected the cells either from the entire colony (motile selection series) or only from the area in the vicinity of the inoculation point (non-motile selection series) (Figure 1A and Methods). Cells recovered from the motile and non-motile selection series were then each used for DNA extraction and as inocula for subsequent passages in swimming medium. This isolation and re-inoculation cycle was repeated for a total of 7 passages for the motile and non-motile selection series.

We followed the formula developed by Wetmore and colleagues for calculating fitness scores (Wetmore et al., 2015). Mutant fitness scores were calculated as the weighted average of the log2 of the ratio of barcode abundance after each passage to that of barcode abundance in the time 0 samples, then the gene fitness scores as the weighted average of fitness scores for all strains with transposon insertions in each gene (Wetmore et al., 2015). Since fitness scores are calculated based on the abundance of typically loss-of-function mutants, positive-value fitness scores are associated with genes whose function are antagonistic for fitness in the specified environment and negative fitness scores are conversely associated with genes that contribute a beneficial function under the imposed selection. We focused our analysis on genes with absolute fitness score values above 1 and considered genes with absolute fitness score values above 2 and absolute t-scores above 5 to have strong influence on fitness and survival. We found 30 genes from the liquid series, 202 genes from the motile selection series and 153 genes from the non-motile selection series with absolute fitness scores above the strong influence threshold (i.e., fitness score absolute values of greater than 2) (Table 2; Supplementary Table S2). We conducted Gene Ontology (GO) annotation enrichment analysis on the 30 genes identified in the liquid series and as expected found amino acid metabolism and organic acid biosynthesis to be the most strongly enriched functional categories (Supplementary Table S2), suggesting these genes are necessary for growth under the nutritional restriction of MOPs minimal liquid media.

We assessed the efficacy of the non-motile selection for mutants with defects in motility by testing the swimming phenotype of randomly sampled individual strains from the populations collected after each passage. For comparison, we determined the swimming phenotypes of individuals sampled from the motile selection series passages where the imposed selection favored motile cells. The motile selection series showed little to no change in motility across the seven passages (Figure 1B). In contrast, the proportion of mutants with motility defects steadily increased in non-motile selection series (Figure 1B). The effect of selection for low motility was also evident in the landscape of gene fitness scores where a clear divergence emerged in the principal component analysis between the non-motile selection series compared to the motile selection series (Figure 1C).

We found 202 genes that strongly influenced fitness and survival in the motile selection series (absolute values above 2). We compared this list of genes to those found in the liquid series control to help categorize the roles of genes in the motile selection series. We found that 29 of the 30 genes identified with high fitness scores (absolute values above 2) in the liquid series also had similar high-absolute-value scores in the motile selection series, indicating that the nutritional requirements are the same in both conditions. We examined the GO annotation of the remaining 173 genes and found motility was the most strongly enriched category of annotations (Fisher’s exact test, FDR-adjusted p < 0.001, Supplementary Table S2). Consistent with our hypothesis that motility is beneficial for growth in soft agar swimming media, we found homologs of 34 out of 47 genes in the P. aeruginosa flagellar synthesis and regulatory pathways (Dasgupta et al., 2003) showed strongly negative fitness scores across the passages (Figure 2; Supplementary Tables S2, S3), with 44 out of these 47 genes conserved in Pst DC3000. The most surprising exception was that motA (PSPTO_4953) and motB (PSPTO_4952), which stood out as having strongly positive fitness scores, indicating that these flagellar motor proteins are paradoxically antagonistic for growth in swimming media.

These results suggest that thresholding genetic fitness scores can indeed identify genes with a strong impact on motility. However, results from this thresholding approach do not compare the relative magnitude of each gene’s impact on swimming motility and thus does not help identify the most likely novel motility regulators. We therefore needed to find a quantitative approach to rank genes in their likelihood of encoding motility regulators regardless of their current annotation.

Non-motile selection should enrich for and increase the relative fitness of mutants with insertions in genes that contribute positively to motility and cause these genes to have increased fitness scores, in contrast to the motile selection. Consistent with this, we observed that known motility genes like cheY (PSPTO_1980) displayed much increased fitness scores in the non-motile selection series compared to the motile selection series. We therefore compared gene fitness scores throughout the non-motile selection series to those throughout the motile selection series, with the expectation that the magnitude of the difference in fitness scores should correlate with the genes’ contribution to motility. We found 600 genes with fitness scores that differed between the two series across passages (t-test, FDR-adjusted p < 0.05), and again found chemotaxis and motility as the most enriched functional category (Supplementary Table S4). This approach successfully identified major components of the chemotaxis regulation cheAYVRW and 10 MCPs but failed to identify flagellar genes from class I to III (Supplementary Table S4) (Dasgupta et al., 2003).

We used the fitness score differences between the two selection series as a quantitative metric for ranking genes according to their likelihood of functioning in motility and conducted principal component analysis (PCA) on the resulting profiles of 600 differentially fitness genes (DFG). We found two groups of motility-altering genes that were strongly separated by PC1, each including known motility-related genes such as chemotaxis response regulator cheY and flagellar stator motAB, respectively (Figure 3A). Clustering results from K-means analysis (X = 4) similarly identified the two groups of genes containing cheY and motAB, respectively, as unique clusters (Supplementary Figure S1), suggesting there is functional similarity among genes within each group. We thus used PC1 value of −0.5 and 0.5 as thresholds to identify individual genes in each of the two groups as being potentially important for swimming motility (vertical lines, Figure 3A). We found 18 genes in the che group that presumably contribute to swimming motility and 13 genes in the mot group that appeared antagonistic to swimming motility (Table 3). Genes in the che group displayed negative fitness scores in the motile selection series, with a complete reversal to positive trend in the non-motile selection series (Figure 3B, left panel) while the mot group genes displayed a positive trend of fitness scores in the motile selection series and shifted to a neutral, slightly negative trend in the non-motile selection series (Figure 3B, middle panel). Genes from both che and mot groups displayed stronger divergence in their fitness scores between motile and non-motile series than the rest of DFGs (Figures 3B,C). These results are consistent with our hypothesis that genes affecting swimming motility would have fitness scores that differ most between the two series.

We next investigated individual genes from the predicted motility-contributing che group because their mutant phenotypes should be easily identifiable in swimming assays. We chose PSPTO_0406 (dipA), PSPTO_1042 (chrR) and PSPTO_4229 (hypothetical protein) as our candidate motility-contributing genes because their fitness score differences between non-motile and motile selection series strongly suggest they positively contribute to swimming motility (Figure 4A).

Motility and biofilm formation are inversely controlled by c-di-GMP levels in the cell, with low c-di-GMP levels favoring biofilm dispersal and motility (Christen et al., 2010; Matsuyama et al., 2016; Khan et al., 2023; Kreiling and Thormann, 2023). Cellular c-di-GMP levels are managed by sets of phosphodiesterase (PDEs) and diguanylate cyclase (DGCs) enzymes that degrade and synthesize c-d-GMP, respectively (Khan et al., 2023; Kreiling and Thormann, 2023). Among the genes identified as contributing to motility, PDE-encoding PSPTO_0406 (dipA) had the largest response to selection for low motility (Figure 4A). We generated a genomic dipA knock-out mutant as well as dipA overexpression plasmid to assess its impact on swimming motility. Both dipA⁺ overexpression strain and ΔdipA mutant were partially impaired swimming motility compared to the wild-type Pst DC3000 (Figures 4B,C). The relative levels of c-di-GMP can be compared using Congo red, which binds exopolysaccharide produced while c-di-GMP levels are elevated (Bordi et al., 2010). We grew our strains of interest on Tryptone medium supplied with Congo Red and found neither the dipA deletion strain or dipA⁺ overexpression strain was morphologically distinguishable from the wild type (Figure 4D) (Bordi et al., 2010). This suggests dipA is unlikely to be the dominant PDE in Pst DC3000 and that the global c-di-GMP level in wild-type Pst DC3000 might already be below detection threshold of the Congo Red assay. Consistent with this result, we saw no difference between these strains using a P _cdrA :gfp reporter (data not shown) (Rybtke et al., 2012; Cerna-Vargas et al., 2019).

PSPTO_1042 (chrR) together with PSPTO_1043 encode an ChrR/RpoE-like anti-sigma/sigma factor pair that responds to oxidative stress (Butcher et al., 2017). PSPTO_1043 was proposed to also regulate motility based on its regulation of DGC encoding PSPTO_2591 (Butcher et al., 2017). If true, deletion of PSPTO_1042 chrR should de-inhibit PSPTO_1043 activity and lead to strong PSPTO_2591 expression and elevated c-di-GMP level. Deletion of PSPTO_1042 (chrR) indeed resulted in strong suppression of swimming motility and the plasmid-based expression fully restored motility level back to wild type (Figures 4B,C). However, we found no visible difference with the ΔchrR mutant growing on Congo Red (Figure 4D).

PSPTO_4229 encodes a hypothetical protein of 273 amino acids with an HD-related output domain (HDOD) (Lee et al., 2016) spanning from amino acids 23 to 210. Deletion of PSPTO_4229 resulted in moderately inhibited swimming motility while over-expression almost abolished swimming completely (Figures 4B,C). Neither deletion nor overexpression of PSPTO_4229 resulted in visible difference in the Congo Red assay (Figure 4D). We attempted to identify protein residues or domains necessary for its activity by screening a plasmid expression library for PSPTO_4229 mutants that no longer inhibit motility. We isolated 40 individual PSPTO_4229 mutant expression constructs that failed to inhibit motility when overexpressed and found 30 mutants with frame-shift or nonsense mutations that are predicted to encode truncated proteins. We found 10 mutants predicted to encode full-length protein with between 2 and 6 single amino acid substitutions each and 3 nonsense mutants predicted to encode more than the first 190 amino acids, among which was a mutant that fully preserved the wild-type sequence until amino acid 193 (Supplementary Figure S2). These results suggest that the C-terminal 80 amino acids contain sequences necessary for structure or function of the PSPTO_4229 protein. We also used the ScanNet tool (Tubiana et al., 2022) to scan PSPTO_4229 and found amino acids 84 and 85 as having a high probability of participating in protein–protein interactions, which was close to a mutation hotspot identified in our fail-to-inhibit mutant screen (Figures 5A,B). We cloned a PSPTO_4229 derivative with R84A R85A amino acid substitutions to test its impact on motility. We found replacing these animo acids eliminated motility interference from PSPTO_4229 overexpression (Figures 5C,D).