Gene ortholog groups were generated using the PyParanoid analysis pipeline v0.4.1 (Melnyk et al., 2019). Peptide sequences from the following RefSeq genome assemblies were used to construct ortholog groups: GCF_000147055.1 (Dda3937), GCF_003403135.1 (DdiaME23), and GCF_014893095.1 (Ddia6719). From these assemblies, RefSeq gene loci were then matched to their corresponding protein names to allow comparison to the Barcode Sequencing (BarSeq) fitness data. Additionally, Clusters of Orthologous Groups (COG) categories for Dda3937 were downloaded from the IMG database (Chen et al., 2019), GenBank gene names were replaced with their corresponding RefSeq names, and added to this ortholog table, allowing putative COG assignments for orthologous genes in D. dianthicola strains.

Strains used in this study are described in Supplementary Table 1. All bacteria were cultured in Luria-Bertani (LB) medium (10 g tryptone, 5 g yeast extract, and 10 g NaCl per 1 L) (Bertani, 1951) at 28°C, except for pure culture E. coli grown at 37°C. When noted, kanamycin was used at a final concentration of 50 μg/ml. Barcoded transposon libraries were constructed by conjugating the barcoded mariner transposon plasmid pKMW3 from the E. coli WM3064 donor library APA752 (Wetmore et al., 2015) into wild-type Dda3937, DiaME23, and Ddia6719. Recipient Dickeya strains were grown as 3 ml LB liquid cultures overnight and 1.5 ml of each culture was subcultured into 30 ml LB without antibiotics. An entire 1 ml freezer aliquot of the E. coli donor library was thawed and used to inoculate 25 ml LB containing 300 μM diaminopimelic acid (DAP) (Sigma-Aldrich, United States) and kanamycin. All cultures were grown to OD₆₀₀ of approximately 1.0, washed with 10 mM KPO₄, and combined in equal amounts before plating donor-recipient pairs each on 50 LB plates containing 300 μM DAP. Conjugations were incubated at 28°C overnight, and exconjugants were then scraped into 40 ml 10 mM KPO₄. This conjugation mixture was then vortexed, spread onto 200 LB kanamycin plates per strain, and incubated at 28°C for 3 days. All colonies were resuspended in 200 ml LB with kanamycin, diluted back to a starting OD₆₀₀ of 0.2, and grown at 28°C with shaking for 6–8 h, to a concentration of approximately 10⁹ CFU/ml (measured as a 1/10 dilution at OD₆₀₀ = 0.15–0.3). Glycerol was added to the library to a final concentration of 15%, and 1 ml aliquots were frozen at –80°C.

For DNA library preparation, genomic DNA from each library was purified from an entire 1 ml cell pellet using the Monarch Genomic DNA Purification Kit (New England Biolabs, United States). Samples were eluted in 50 μl nuclease-free water. Purified DNA was quantified on a NanoDrop One (Thermo Fischer Scientific, United States), and 500 ng DNA was used as input for the NEBNext Ultra II FS DNA Library Prep kit (New England Biolabs, United States), following the manufacturer’s instructions with modifications as follows. For enzymatic DNA fragmentation, a 12-min incubation time was used. DNA fragments were size selected using AMPure XP magnetic beads (Beckman Coulter, United States) at the recommended ratios 0.4X and 0.2X. We used a modified version of the protocol described in Wetmore et al. (2015), with a two-step PCR used to enrich for transposon insertion sites, based on (Rubin et al., 2021). A custom splinkerette adapter was ligated to fragmented DNA, prepared by annealing oligos:/5Phos/G*ATCGGAAGAGCACACGTCTGGGTTTTTT TTTTCAAAAAAA*A and G*AGATCGGTCTCGGCATTCCC AGACGTGTGCTCTTCCGATC*T (Rubin et al., 2021). Between rounds of PCR and before submitting for sequencing, DNA was cleaned by binding to AMPure XP magnetic beads, using a bead ratio of 0.9X and eluted in 15 μl 0.1X TE buffer for intermediate steps and 30 μl 0.1X TE for sequencing. Finally, the sequencing library was quantified using a Qubit dsDNA HS assay kit (Thermo Fischer Scientific, United States). DNA libraries were submitted for sequencing at the Biotechnology Resource Center (BRC) Genomics Facility at the Cornell Institute of Biotechnology on an Illumina MiSeq to check library quality, followed by sequencing on a NextSeq 500 (Illumina, Inc., United States). All mapping used single-end sequencing for 150 bp fragments.

Sequence data were analyzed using the scripts from the FEBA package v1.3.1 (Wetmore et al., 2015). MapTnSeq.pl was used to identify the barcode and location in the genome of each read with identifiable transposon sequence from both MiSeq and NextSeq reads, based on the “model_pKMW3.2” sequence. DesignRandomPool.pl was used to assemble the mutant pool using barcodes seen in a single location 10 or more times. All TnSeq mapping and BarSeq fitness calculation code is available at http://bitbucket.org/berkeleylab/feba/ (Wetmore et al., 2015). Mapping scripts were run on a Cornell University BioHPC Cloud 40-core Linux (CentOS 7.6). server with 256 GB RAM.

Using the output from MapTnSeq.pl, gene essentiality predictions were made using https://bitbucket.org/berkeleylab/feba/src/master/bin/Essentiality.pl and the function “Essentials” from https://bitbucket.org/berkeleylab/feba/src/master/lib/comb.R (Wetmore et al., 2015). Using the median insertion density and the median length of genes > 100 bp, this method calculates how short a gene can be and still be unlikely to have no insertions by chance (P < 0.02, Poisson distribution); genes shorter than this threshold are then excluded (Price et al., 2018). For the Dickeya strains examined here, the minimum gene length for a gene to be predicted as essential for growth in LB was 175 bp (Dda3937) or 150 bp (DdiaME23 and Ddia6719). Protein-coding genes are then considered to be essential or nearly essential if there are no fitness values and the normalized central insertion density score and normalized read density score as computed by the FEBA package were < 0.2 (Price et al., 2018).

For a given BarSeq experiment, a single transposon library freezer aliquot was thawed and recovered in 25 ml LB containing kanamycin at 28°C until OD₆₀₀ = 0.5–0.7, approximately 6–8 h. At this point, two 1 ml cell pellets were frozen as time0 measurements, and the remaining culture was washed in 10 mM KPO₄ and used to inoculate experimental samples.

All in vitro cultures were grown in 1 ml volumes in 24-well plates. In each well, 50 μl starter culture at 0.3 OD₆₀₀ was added to 950 μl medium containing kanamycin. Media tested were LB, Potato Dextrose Broth (PDB)(pH 5) (Sigma-Aldrich, United States), and M9 minimal medium (M9) as described in M9 Minimal Medium (Standard) (2010) but containing 0.4% glycerol instead of 0.4% glucose. Plates were incubated at 28°C with shaking at 200 rpm. After 1 day (LB and PDB) or 2 days (M9), each 1 ml sample was pelleted and frozen prior to genomic DNA extraction.

Prior to inoculation, all tubers were rinsed and then surface sterilized by submerging in 70% ethanol for 10 min, followed by two washes with distilled water. Inoculum was standardized to 10⁹ CFU/ml by measuring a 1/10 culture dilution at OD₆₀₀ = 0.3 (corresponding to 10⁸ CFU/ml), and 10 μl was inoculated in two replicate stab wounds created by pushing a 200 μl pipet tip roughly 3 mm into each tuber. Six replicate tubers were used for each bacterial strain and potato cultivar. Inoculated tubers were stored in plastic bags at 28°C. Two days after inoculation, ∼2 cm length cores were taken at each site of inoculation using a 1 cm diameter cork borer. Duplicate cores from each tuber were pooled in 8 ml 10 mM KPO₄ and shaken at 200 rpm at 28°C for 10 min. For each sample, 2 ml of bacterial suspension were pelleted and frozen prior to DNA extraction.

Genomic DNA was purified from cell pellets using the Monarch Genomic DNA Purification Kit (New England Biolabs, United States). Purified DNA samples were eluted in 30 μl nuclease-free water and quantified on a Nanodrop One (Thermo Fischer Scientific, United States). After gDNA extraction, the 98°C BarSeq PCR as described in Wetmore et al. (2015) was used to specifically amplify the barcode region of each sample. 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 (46–49 samples per pool), 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, 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). Scripts to calculate gene fitness values were run on a Cornell University BioHPC Cloud 40-core Linux (CentOS 7.6). Server with 256 GB RAM. 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 for a given transposon library 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 (Wetmore et al., 2015). Gene fitness values were normalized across the chromosome so that the median gene fitness value was 0. All experiments described here passed previously described quality control metrics (Wetmore et al., 2015; Price et al., 2018).

We focused on genes having fitness values > 1 or < –1 and absolute t-like test statistic > 4. This t score is an estimate of the reliability of the fitness measurement for a gene, and is equal to the fitness value divided by the square root of the maximum variance calculated in two ways (Wetmore et al., 2015). With these cutoffs, we also calculated gene fitness values comparing replicate time0 samples (Price et al., 2018; Liu et al., 2021). Across 6 (each Dda3937 and DdiaME23) and 2 (Ddia6719) replicate time0 samples, 0 gene fitness values had fitness > 1 or < –1 and absolute t > 4. Data were analyzed in R v4.0.3 (R Core Team, 2017) and visualized using the package ggplot2 v3.3.5 (Wickham, 2016). The principal components analysis was performed on the gene fitness matrix for each Dickeya strain using the R function prcomp, which performs centered singular value decomposition.

All raw Illumina reads used for mapping and fitness assays have been deposited in the Sequence Read Archive under BioProject accession number PRJNA692477. Individual sample accession numbers are listed in Supplementary Table 2. Annotated scripts used for computational analysis are available at http://github.com/tylerhelmann/dickeya-barseq-2021. Experimental fitness values are publicly available on the interactive Fitness Browser at http://fit.genomics.lbl.gov.

To enable direct comparison of gene fitness measurements between strains, we constructed a database of homologous gene families using the PyParanoid analysis pipeline (Melnyk et al., 2019). Based on clustering of all predicted protein sequences from Dda3937, DdiaME23, and Ddia6719, 3,821 total homolog groups were identified, representing 88.1% of the total input sequences. Of these, 3,310 groups contained single-copy genes in all three strains. For each group, gene loci, protein identifiers, and gene descriptions are listed in Supplementary Table 3. This table also contains COG assignments matched from the Dda3937 IMG genomic annotation (Chen et al., 2019).

To measure contributions of individual genes to fitness, we constructed barcoded transposon mutant libraries in the Dickeya strains using a barcoded mariner E. coli donor library (Wetmore et al., 2015). These libraries ranged in size from 334,893 to 541,278 mapped genomic insertional strains, with 37–62 median strains per hit gene (Table 1). Of the three strains tested, only one gene in the Ddia6719 genome did not contain any TA dinucleotide sites and was therefore inaccessible to the mariner transposon. Mapped insertions were evenly distributed across the chromosome of each strain (Supplementary Figure 1).

Based on analysis of the TnSeq mapping data, essential genes were predicted using the FEBA RB-TnSeq analysis pipeline (Wetmore et al., 2015; Price et al., 2018). We identified 374–426 genes per strain that are likely to encode essential or near-essential genes for growth in LB (Supplementary Table 4). Using the ortholog group assignments for these genes, 316 of these predicted essential genes (74–84%) are shared among all three strains (Supplementary Figure 2). Most predicted essential genes are in the functional categories of “Translation, ribosomal structure, and biogenesis,” “Cell wall/membrane/envelope biogenesis,” “Coenzyme transport and metabolism,” “Energy production and conversions,” and “Replication, recombination, and repair” (Supplementary Table 5).

To generate genome-wide gene fitness values for the barcoded transposon libraries, each strain was grown in the rich media LB and Potato Dextrose Broth (PDB) as well as M9 minimal medium supplemented with 0.4% glycerol (Supplementary Figure 3). Strain fitness values were calculated as a log2 ratio of barcode abundance following sample growth with barcode abundance measured in the time0 duplicate samples. Gene fitness is the weighted average of individual strain fitness values (Wetmore et al., 2015). For fitness calculations, insertions in the first and last 10% of coding regions were excluded, with insertions in the remaining 80% of the gene considered “central.” While 91–92% of genes in all strains contained centrally mapped insertions, not all genes were used in fitness calculations due to low read or insertion abundance. We focused our analysis on genes with fitness values > 1 or < –1, and absolute t-score > 4 (Supplementary Table 6). Across all conditions, we calculated fitness values for 3,705 (Dda3937), 3,761 (DdiaME23), and 3,528 (Ddia6719) genes, representing 88, 90, and 86% of the total genes in each strain, respectively.

Principal component analysis showed gene fitness values of the three tuber conditions overlapped (Figure 1), and so these samples were jointly considered as a single “Tuber” condition for some subsequent analyses.

As the libraries were constructed on LB medium, relatively few mutations deleterious in LB were maintained in the populations (Figure 2). Figure 3 presents these data split by COG category. Limited mutations in genes categorized as “cell wall/membrane/envelope biogenesis” (mdoGH) and “cell cycle control, cell division, chromosome partitioning” (ftsX) were present in the mapped populations but generally detrimental in LB for all three strains (Supplementary Figure 4). Even in LB, some variation was apparent between strains, such as disruptions in the gene encoding the cell division protein ZapB which decreased competitive fitness in Dda3937 but not DdiaME23 or Ddia6719 (Supplementary Figure 4).

The rich medium PDB provided a very different gene fitness profile than LB, likely due to nutritional differences and its slightly acidic nature (pH 5). For all three strains, more genes were detrimental (fit > 1) and beneficial (fit < −1) in PDB relative to LB (Figure 2). Genes in diverse metabolic categories contributed to competitive fitness, including “amino acid transport and metabolism,” “carbohydrate transport and metabolism,” “cell wall/membrane/envelope biogenesis,” “coenzyme transport and metabolism,” “inorganic ion transport and metabolism,” “nucleotide transport and metabolism,” “signal transduction mechanisms,” “transcription,” and “translation, ribosomal structure, and biogenesis” (Figure 3). For example, in all three strains oligopeptidase A and the low affinity potassium transporter Kup were specifically important in PDB for growth (Supplementary Figure 5). Disruptions in the two-component system RtsAB were specifically beneficial for Dda3937 in PDB, as were disruptions in the zinc uptake transcriptional repressor Zur (Supplementary Figure 5). Though LB and PDB are both complex rich media, pH and/or specific available nutrients differed enough to clearly separate the gene fitness profiles for Dda3937 and Ddia6719, though not for DdiaME23 (Figure 1).

In the minimal medium M9 containing glycerol as a carbon source, important genes included categories such as “amino acid transport and metabolism,” “carbohydrate transport and metabolism,” “coenzyme transport and metabolism,” and “nucleotide transport and metabolism.” While many amino acids were limiting in both M9 and tuber samples, arginine biosynthetic genes (argCEFGH) were uniquely important in M9, suggesting the presence of available arginine in tubers (Figure 4). Conversely, mutations in many “cell motility” genes had a large positive effect, though this effect was often limited to Dda3937 (Figure 5). This is indicative of the limited energy available in minimal medium, and the high energy cost of motility.

To calculate genome-wide gene fitness values in an ecologically and economically relevant condition, we inoculated the transposon libraries into tubers of three potato cultivars: “Atlantic,” “Dark Red Norland,” and “Upstate Abundance.” As each transposon library contains over 300,000 unique strains, we inoculated approximately 10⁷ cells into each tuber (10 μl of a 10⁹ CFU/ml solution). After 2 days incubation at high humidity, we recovered cells by streaming for barcode sequencing and calculation of gene fitness values. Many genes involved in amino acid biosynthesis that were important for growth in M9 were also important in tubers (leuAC, thrC, serB), highlighting potentially limiting factors for growth during potato soft rot (Figure 6).

The pectin degradation protein 2-dehydro-3-deoxy-D-gluconate 5-dehydrogenase KduD was specifically important for growth in tubers (Figure 7). Interestingly, we identified several putative DNA-binding or helix-turn-helix transcriptional regulators where mutant strains had strain-specific increased fitness in tubers (Figure 7). Insertions in the Ddia6719 helix-turn-helix transcriptional regulator HGI48_RS01985 increased fitness in tubers, while insertions in the paralog HGI48_RS02000 had no effect on fitness. There is no ortholog for this gene in Dda3937, and the ortholog in DdiaME23 had no disruption phenotype in any condition tested (Figure 7).