Pantoea ananatis strains (n = 92) used in this study were isolated from diverse sources; weeds, thrips, and onion tissue (foliage, bulbs, and seeds) in Georgia from 1992 to 2019. The metadata for each strain such as the source, year of isolation, and county of origin in Georgia for these strains and their distribution within each category are listed (Table 1).

Among the strains used, 55 strains (59.8%) were isolated from onion foliage or bulb tissue, which constituted the majority of the strains. The remaining 38 strains (40.2%) were isolated from other diverse sources including weeds, and thrips. This is followed by the weeds Richardia scabra L. (8.7%; 8/92), Digitaria spp. (6.5%; 6/92), and Verbena bonariensis (4.3%; 4/92). The strains from various plant sources constituted 8.7% (8/92) of the total strains studied (Table 1; Figure 1B). Strains from thrips (Frankliniella fusca and F. occidentalis) constituted the remaining 13% (12/92). These curated strains were initially identified as P. ananatis by their colony morphology and physiological characteristics such as being: Gram-negative, facultatively anaerobic, positive for indole production, and negative for nitrate reductase and phenylalanine deaminase. Further confirmation was done using a P. ananatis-specific PCR assay as described earlier (Walcott et al., 2002).

Pantoea ananatis inoculum for foliar pathogenicity phenotyping was prepared by transferring single colonies of each bacterial strain from 24 h-old cultures on nutrient agar (NA) medium to nutrient broth (NB). The broth was shaken overnight on a rotary shaker (Thermo Scientific, Gainesville, FL, United States) at 180 rpm. After 12 h of incubation, 1 ml of each bacterial suspension were centrifuged at 5000 × g (Eppendorf, Westbury, NY, United States) for 2 min. The supernatant was discarded, and the pellet was re-suspended in deionized water. Inoculum concentration was adjusted using a spectrophotometer (Eppendorf, Westbury, NY, United States) to an optical density of 0.3 at 600 nm [≈ 1 × 10⁸ colony forming unit (CFU)/ml].

Deletion of the P. ananatis PANS 02-18 pepM_pgb pepM_HiVir genes was conducted by two-step allelic exchange as described by Stice et al. (2020). In brief, approximately 400-bp flanking regions to the targeted genes were directly synthesized as a single joined sequence by Twist bioscience and cloned via BP clonase II using primer-introduced attB1/2 recombination sites into the pR6KT2G Gateway^® compatible sacB-based allelic exchange vector. These deletion constructs were introduced into PANS 02–18 via biparental conjugation with the RHO5 E. coli strain and single crossover events were recovered via gentamicin selection. Second crossover events were recovered via liquid sucrose counter-selection and identified by screening for backbone eviction based on loss of gentamicin resistance and the formation of white colonies on X-gluc. Deletion mutants were identified and confirmed based on PCR and amplicon sequencing using independent primers designed to amplify from genomic regions adjacent to the 400-bp deletion flanking regions.

Pathogenic potential of P. ananatis strains were initially phenotyped on onion scale using a previously described red scale necrosis (RSN) assay (Stice et al., 2018). Red onions (cv. Red Burgundy) were surface sterilized with 70% ethanol and the outermost scale sliced to approximately 3 cm × 4 cm. The resulting scales were set on a sterile petri dish or on sterilized microtube trays, with the bottom covered with sterilized paper towels pre-moistened with distilled water. Each onion scale was then wounded via direct penetration with a sterilized needle and inoculated with 10 μl of approximately 1 × 10⁶ CFU/ml inoculum of P. ananatis. A known onion-pathogenic strain (PNA 97-1) was used as a positive control (Gitaitis and Gay, 1997). Sterile water was used as a negative control. The resulting petri dishes were then laid in an aluminum tray (46 cm × 25 cm × 10 cm) and covered with a plastic lid. These onion scales were then incubated for 5 days in the dark. The area of pigment clearing, and necrotic lesions were measured at 7 days post-inoculation. Strains that did not clear the red anthocyanin pigment or developed necrotic lesions, were declared non-pathogenic. Strains that caused necrosis along with pitting, with a visible zone of pigment clearing were considered pathogenic. Three replications were performed for each strain and in total two experiments were conducted.

Foliar pathogenicity and aggressiveness of P. ananatis strains (n = 92) were determined on A. porrum (cv. King Richard) and A. fistulosum × A. cepa (cv. Guardsman) under controlled greenhouse conditions. P. ananatis strain (PNA 97-1) was used as a positive control for both Allium species (14,31). Seedlings were established in plastic pots (T.O. plastics, Clearwater, MN) with dimension of 9 cm × 9 cm × 9 cm (length × breadth × height) containing a commercial potting mix (Sta-green, Rome, GA, United States). The seedlings were maintained under greenhouse condition at 25–28°C and 70–90% relative humidity with a light:dark cycle of 12:12 h. Osmocote smart release plant food (The Scotts Company, Marysville, OH, United States) was used for periodic fertilization. Bacterial strains were maintained on NA plates and inoculum was generated as described above. Once the primary leaf of each Allium spp. reached 9 cm, seedlings were inoculated using a cut-tip method as described previously (Dutta et al., 2014a). Briefly, a wound was created by cutting the central leaf (2 cm from the apex) with a sterile pair of scissors. Using a micropipette, a 10 μl drop of a bacterial suspension containing 1 × 10⁸ CFU/ml (1 × 10⁶ CFU/leaf) was deposited at the cut-end. Seedlings inoculated with sterile water as described above were used as negative control. Three replications per strain per host were used for one experiment and a total of two independent experiments was conducted. The seedlings were observed daily for symptom development until 5 days post-inoculation (DPI) and were compared with the foliar symptoms displayed by the positive control on each Allium species. The aggressiveness of P. ananatis strains was determined based on the lesion length on each Allium spp. For A. porrum, strains that caused a lesion length of 0.2–0.5 cm were considered less aggressive, 0.5–0.9 cm moderately aggressive, and > 1 cm highly aggressive. For A. fistulosum × A. cepa, strains were considered highly aggressive when a lesion length of > 1.4 cm was observed. Lesion lengths ranging from 0.7 to 1.4 cm were considered as moderately aggressive and strain with lesion length < 0.7 cm was regarded as less aggressive. Bacterial strains that did not display any lesion were considered as non-pathogenic. To confirm if the symptoms were caused by P. ananatis, bacteria were isolated from the region adjoining the symptomatic and healthy tissue on PA-20 semi-selective medium and incubated for 5–7 days at 28°C (Goszczynska et al., 2006). Presumptive colonies were further confirmed using a P. ananatis-specific assay as mentioned above (Goszczynska et al., 2006). Further, strain identity from randomly isolated colonies from A. porrum and A. fistulosum × A. cepa were confirmed by their DNA fingerprints using repetitive extragenic palindrome (rep)-PCR as previously described (Dutta et al., 2014a, b).

Genomic DNA was extracted utilizing the E.Z.N.A bacterial DNA kit Omega Bio-Tek (Norcross, GA). A 50 μl of DNA (50 ng/μl) per sample was used for library preparation as per the manufacturer’s instructions at Novogene Bioinformatics Technology Co. Ltd. (Beijing, China). Genomic DNA of each sample was randomly sheared into short fragments of about 200–400 base pairs (bp). The obtained fragments were subjected to library construction using the NEBNext^® DNA Library Prep Kit. After end repairing, dA-tailing, and further ligation with NEBNext adapter, the required fragments (in 200–400 bp size) were PCR enriched by P5 and indexed P7 oligos. The library was subsequently sequenced on Illumina NovaSeq 6,000 platform (Illumina Inc., San Diego, CA, United States). Pair-end sequencing were performed with the read length of PE 150 bp at each end. The raw fastq reads obtained were quality filtered. FastQC was used to assess the raw fastq files. Reads were filtered utilizing Trimmomatic (v. 0.36). The read data were filtered to remove low quality reads/bases and trimmed for reads containing primer/adaptor sequences using Trimmomatic’s ILLUMINACLIP paired end mode with seed mismatches set to two, palindrome clip threshold of 30, and a simple clip threshold of ten (Bolger et al., 2014). Further, all 5′ and 3′ stretches of ambiguous ‘N’ nucleotides were clipped to ensure high quality reads via setting both leading and trailing options to 3. The window size for sliding window was set to 4 with a required threshold of 30.

Trimmed data were re-assessed using FastQC and further used for genome assembly followed by pan-genome analyses. Further, all contigs ≤ 500 bp were removed using Seqtk (1.3). The cleaned reads were assembled using SPAdes (v. 3.15.3; Bankevich et al., 2012). Both the paired and unpaired data were used in assembly at default settings. The scaffolds of the respective 92 P. ananatis strains were annotated using Prokka (v. 1.14.5; Seemann, 2014). The resulting .gff files were used in the downstream pan-genome analysis.

Average nucleotide identity was determined using FastANI (Jain et al., 2018). KEGG gene ontology (GO) assignment was conducted using BioBam BLAST2GO pipeline (Götz et al., 2008; BioBam, 2022). A phylogenetic tree of single nucleotide polymorphisms (SNP’s) of core genome was generated utilizing PanSeq at default settings and RAxML alignment with 10,000 bootstrap replicates as per the previously reported methodologies (Laing et al., 2010; Stamatakis, 2014; Agarwal et al., 2021). The RAxML Boostrap random number used was 9,595, with parsimony random seed of 5,959. For Coinfinder and Roary plots, a tree was generated using FastTree 2.1.11 at default settings utilizing Roary’s core gene alignment (Price et al., 2009).

All annotated genomes passed quality control and were used as inputs in ROARY (v. 3.12.0) at default settings that aided in generating a pan-genome with core and accessory genes. The complete pan-genome matrix in the form of presence and absence variant was used as inputs in SCOARY (v. 1.6.16). The GWAS analysis was conducted using the SCOARY program, which determined association between pangenomes and observed phenotypes (Page et al., 2015; Brynildsrud et al., 2016). This program was operated twice separately on each host plant of interest, once at default parameters, and a second time with a forced maximum value of p 0.05 across all testing parameters. For the gene association/dissociation analysis, complete pangenome of 92 P. ananatis strains were used as inputs for Coinfinder (v. 1.0.1). This program generated both gene-pair associations and dissociations with modification to association significance increased to a value of p of 0.1, and default settings (p = 0.05) for dissociation as previously described (Whelan et al., 2019). Direct comparisons of genetic sequences were performed using the Clustal Omega online server at default settings (Madeira et al., 2022).

A single colony of P. ananatis strain PNA 15–3 (with putative T3SS) and a strain of Pantoea stewartii subsp. indologenes (20GA0713; positive control for T3SS) was suspended and grown overnight in modified Coplin medium (Asselin et al., 2018). Approximately 100 μl of the overnight culture was syringe-infiltrated into the tobacco leaf and the resulting infiltrated area was marked with a black marker. A sterile Coplin lab medium was used as a negative control. The symptom was observed at 48-h post inoculation (hpi) when the image was taken. This experiment was repeated twice.

Phenotyping of 92 P. ananatis strains displayed variability in the level of aggressiveness on both Allium spp. (Table 1). Variations in P. ananatis pathogenicity and aggressiveness on two Allium spp. were considerable (Figures 1A,B; Table 1). Strains screened in this study belonged to different isolation sources (Figure 1C; Table 1). Using the RSN disease phenotyping assay, we observed 61.3% (57/92) of P. ananatis strains displayed typical necrosis of red onion scale whereas 38.7%% (36/92) of strains did not cause necrosis (Figure 1D).

When P. ananatis strains were screened on A. fistulosum × A. cepa, 20.4% (19/92) and 44.1% (41/92) were found to be non-pathogenic and mildly aggressive, whereas 25.8% (24/92) and 9.7% (9/92) of the strains were identified as moderately aggressive and highly aggressive, respectively (Figure 1E). In contrast, on A. porrum, 45.7% (42/92) and 37% (34/92) of the strains were non-pathogenic and mildly aggressive, respectively. Interestingly, a much lower proportion of the strains; 14.1% (13/92) and 3.3% (3/92) identified as moderately aggressive and highly aggressive, respectively on A. porrum (Figure 1F). The percentage of strains that were pathogenic on A. porrum but non-pathogenic on A. fistulosum × A. cepa was only 2.1% (2/92). In contrast, 27% (25/92) of the strains that were pathogenic on A. fistulosum × A. cepa were non-pathogenic on A. porrum. Interestingly, 4.3% (4/92) of strains were highly aggressive on A. porrum but less aggressive on A. fistulosum × A. cepa whereas 9.7% (9/92) of strains were highly aggressive on A. fistulosum × A. cepa but less aggressive on A. porrum. Percentage of strains that were moderately to highly aggressive on both Allium spp. was 4.3% (4/92) whereas 18.8% (17/92) of the strains were non-pathogenic on both hosts tested. All the strains isolated from symptomatic A. porrum or A. fistulosum × A. cepa were identified as P. ananatis by recovery on PA-20 semi-selective medium and a P. ananatis-specific PCR assay as described above. P. ananatis colonies were not recovered from any of the negative control seedlings on PA-20 medium indicating no potential cross-contamination among the inoculated strains.

When comparing pathogenicity results with RSN results the 57 RSN-positive strains, 63.2% (36/57) of strains were pathogenic on both A. porrum and A. fistulosum × A. cepa whereas 0% (0/57) and 26.3% (15/57) of strains were only pathogenic on A. porrum or A. fistulosum × A. cepa, respectively. The remaining RSN-positive strains were non-pathogenic in the leaf tip necrosis assay on both hosts 10.5% (6/57). Among the RSN-negative strains, 41.6% (15/36) strains were non-pathogenic on both hosts, whereas 25% (9/36) were pathogenic on both hosts. Also, 5.5% (2/36) and 27.7% (10/36) of strains were pathogenic on only A. porrum and A. fistulosum × A. cepa, respectively. These results indicate that there is a discrepancy between the RSN phenotype and the foliar necrosis phenotype.

Post-sequencing, 1,594,092,228 raw reads were obtained and after stringent quality-filtering and trimming nearly 87% of the total reads (1,413,144,772 quality reads) were retained. The FastQC results indicated the sequence quality “passed,” as the majority of per-nucleotide and sequence qualities achieved high scores with no issues reported. For example: a good score can be ascertained with an average quality score of 30–40, with an exponentially increasing quality score distribution. Sequences that failed to incorporate into the final pangenome, or showed signs of contamination, were removed from the study. All P. ananatis sequences used in this study were submitted to NCBI (bio project PRJNA825576). Their corresponding accession numbers are listed in the Supplementary Table 1. The ANI matrix indicates that the strains investigated and utilized for genome analysis were indeed P. ananatis as the scores were more than 95% (Figure 2). The smallest ANI value was between PNA 18-6S vs. PANS 19-17 with an ANI score of 96.2% while the highest value was observed with PNA 99-6 vs. PNA 99-7 with a score of 99.4% (Figure 2).

An overview of the final pangenome shows a core genome (occurs in 99% or more genomes, N > = 91) of 2,914 genes, a soft-core (occurs in 95–99% of genomes, N = 87 to 91) of 687 genes, a shell genome of 1833 genes (occurs in 15 to 95% of the genomes, N = 14–87), and a cloud genome (occurs in 0 to 15% of genomes, N = 0–14) of 9,196 genes for a total of 14,630 genes (Figure 3A). Details of the number of core and accessory genes contributed by each strain are shown in Figure 3B. A visual representation of the total presence and absence of genes within the pangenome where genomic differences in accessory components as well as the homogeneity of the core genome across the strains can be observed (Figure 3C). The phylogenetic tree produced that was used as input for both the gene-pair coincidence (GPC) analysis and the ROARY plots script is shown in Supplementary Figure 1.

The assignment of GO terms resulted in 19,323 annotations (Supplementary Figure 2). Among the genes annotated and assigned to biological processes (BP) within the pangenome, 3,828 are dedicated to cellular processes, 3,109 to metabolic processes, 930 to localization, 761 to biological regulation, 726 to the regulation of biological processes, 514 to the response to stimulus, 191 for signaling, 152 for the interspecies interaction between organisms, 101 for locomotion, 65 to viral processes, 46 for the negative regulation of biological processes, 43 for detoxification, 42 to developmental processes, 35 for positive regulation of biological processes, 36 for reproduction, 12 for nitrogen utilization, 5 for carbon utilization, 4 for multicellular organismal process, and one for immune system process (Supplementary Figure 2).

Among the genes that are assigned to molecular functions (MF), 3,182 are for catalytic activity, 2,713 for binding activity, 676 for transporter activity, 263 ATP-dependent activity, 252 with transcription regulator activity, 89 with molecular transductor activity, 82 with structural molecule activity, 37 with small molecule sensor activity, 34 for antioxidant activity, 30 with toxin activity, 23 with translocation regulation activity, 15 with molecular function activity, 12 for cytoskeletal motor activity, 6 for molecular carrier activity, and finally one assigned with nutrient reservoir activity. Among the genes assigned to cellular components, 2,894 are assigned as a cellular anatomical entity, 210 are protein-containing complexes, and 2 are virion components (Supplementary Figure 2). Further differentiation of these groups are available in the Supplementary Figure 3.

To determine the relationship between phylogeny of bacterial strains and their pathogenicity on Allium hosts, a phenotypic tree based on SNPs of core genes was constructed using RAxML and PanSeq (Figure 4) and is visually represented using the Interactive Tree of Life online tool (Letunic and Bork, 2016). When assessing the phylogenetic tree in its totality, it is difficult to determine a precise pattern except for strains from the same year of isolation tend to group together. This potentially indicates that these strains in the same group are genetically closely related. Despite this lack of obvious pattern in the overview of the phylogenetic tree, there are several clades where the terminal taxa are sorted based on their pathogenicity. Overall, these results provide support that changes in pathogenicity are not the result of strain lineage, but rather an expansive accessory genome.

A pangenome utilizing ROARY was built, and the strength of gene association to the pathogenic phenotype on seedlings (A. porrum and A. fistulosum × A. cepa) was calculated using SCOARY. A total of 836 genes were found associated with RSN phenotype in A. fistulosum × A. cepa (p < =0.05).

Further, Benjamini-Hochberg correction (BHC) and Bonferroni correction were utilized in an attempt to avoid false positive associations (p < =0.05) and as a result only 50 genes were found significantly associated with the A. fistulosum × A. cepa foliar pathogenic phenotype. However, due to the GWAS results for A. porrum (discussed below) we found that the statistical association is not reflective of the biological role for some of the laboratory-based verified genes of interest. As such we included potential false positive associations in our analysis, but only focused on those that are shared between multiple analyses or those that seem to be similar to known gene clusters of interest. The GWAS results are listed in Supplementary File 1 in order of their statistical association. Within this total set of 836 genes, we found two divergent copies of phosphoenolpyruvate mutase (pepM) genes, annotated as “phosphonopyruvate hydrolase.” The first gene hvrA/pepM belongs to the previously described ‘HiVir’ cluster. Among the top twenty significantly associated genes with corrected p-values in P. ananatis affecting A. cepa × A. fistulosum, the HiVir cluster genes were found to be prominent (Supplementary file 1).

The second pepM gene belonged to a separate gene cluster (pgb) previously described by Polidore et al. (2021), which ranked at 451 based on naïve significance value. Five phosphatase genes, one gene related to chemotaxis, three related to virulence-region associated virB, and several genes in the previously described alt gene cluster were also associated with the pathogenicity (Jain et al., 2018; Stice et al., 2018, 2020). Two copies of the fliC gene, which encodes the flagellin monomer, were also found to be related with P. ananatis’ pathogenicity on A. fistulosum × A. cepa (Macnab, 2003; He et al., 2012). Flagellar motility has been previously observed to be important for onion leaf virulence (Weller-Stuart et al., 2017). Among the associated genes, we also screened for genes (annotated or hypothetical) to assess if they occur in clusters. Within the top 50 significantly associated genes we found at least five hypothetical gene clusters (group_4714–4719, group_3715–3726, group_5180–5182, group_5653–5660, group_5704–5778) as well as the HiVir gene cluster (Supplementary File 1).

Using A. porrum pathogenic strains, a total of 243 genes were found associated (naïve significance of p < = 0.05) with the pathogenic phenotype. However, none of the predicted genes were associated with the phenotype when the Bonferroni correction, and BHC were applied. When the naïve p-values were selectively screened for previously described genes known for pathogenicity and virulence in onion (HiVir genes, alt), we observed the HiVir cluster to be significantly associated with the phenotype; however, it ranked lower (rank: 91–101) compared to other annotated or hypothetical genes. The loss of statistically significant genes post-correction is likely due to the overall “potential resistance phenotype” this particular cultivar of A. porrum seems to display in our dataset. In following sections, we describe the biological relevance of the HiVir cluster for foliar pathogenicity in A. porrum and as such we believe that the naïve p-values are sufficient for further investigation despite the potential for false positives.

We also found significantly associated genes (n = 123), which were shared between the two hosts, with 48 of the top 50 associated genes (statistically significant post-correction) in A. fistulosum × A. cepa occurred in both GWAS results (Supplementary File 1). Some of the known genes that were shared between the A. fistulosum × A. cepa and A. porrum include the entire HiVir gene cluster, pemK_2 (mRNA interferase), soj (sporulation initiation inhibitor protein), parM (Plasmid segregation protein), umuD (protein UmuD), tibC (glycosyltransferase), ycaD (uncharacterized MFS transporter), dadA (D-amino acid dehydrogenase 1), frbC (2-phosphonomethylmalate synthase), amiD (N-acetylmuramoyl-L-alanine amidase), and rfbB (dTDP-glucose 4,6-dehydratase). The genes that constitute the thiosulfinate tolerance cluster (alt) only appeared in A. fistulosum × A. cepa association with the phenotype with the following annotations; xerC (tyrosine recombinase XerC), altA/nemA (N-ethylmaleimide reductase), gor (glutathione reductase), altJ/osmC (peroxiredoxin OsmC), and altD/trxA (thioredoxin; Supplementary File 1). Genes that are members of the larger OVRA region, but not alt-specific genes in were also associated with the pathogenicity phenotype, and they include rbsC (ribose import permease protein RbsC), rbsB (ribose import binding protein), rbsA (ribose import ATP-binding), and altD/trxA (thioredoxin; Supplementary File 1). These results indicate that despite the lack of significance post-correction for the A. porrum GWAS results, there are still clusters of biological relevance for other Allium spp. might be important. Due to the proven role of the HiVir cluster in foliar pathogenicity, the remainder of this manuscript will focus on genes that are shared among multiple analyses, independent of the corrected significance values.

Gene-pair association of the P. ananatis pan-genome resulted in a total of 165 genes separated into 39 individual groups (Table 2; Supplementary File 2). Of the 165 genes, 45 genes are shared with the genes that are predicted based on GWAS for pathogenic phenotype on A. fistulosum × A. cepa and only two genes are shared with the genes that are predicted via GWAS for A. porrum pathogenicity (Table 2; Supplementary File 2). An overview of the associative Coinfinder output can be seen in Figure 5. Of the groups that also occurred on the GWAS analysis, 9/10 are saturated with genes that only associate with the pathogenic phenotype for A. fistulosum × A. cepa, and only one group is saturated with genes that associate with the pathogenic phenotype on A. porrum. These results indicate that in our pangenome there is a stronger associative pressure on genes that are specific to one host or the other, and there is no evidence for gene association between genes that associate for both hosts.

Gene-pair dissociation of the P. ananatis pangenome resulted in 255 genes separated into 50 groups of dissociated genes (Table 3; Supplementary File 2). Of the 255 genes, 22 are shared with the genes associated with the pathogenic phenotype on A. fistulosum × A. cepa as predicted by GWAS, whereas only three genes are shared with the pathogenic phenotype on A. porrum. Components with dissociating gene-pairs that also occur on the GWAS output are summarized in Table 3. A full summary of the genes, their coincidence values, and their groups can be found in Supplementary File 2. An overview of the dissociative Coinfinder output can be seen in Figure 6. Among these groups, 4/6 show a dissociative relationship between genes that associate with the A. fistulosum × A. cepa disease phenotype and genes that do not associate with the pathogenic phenotype on either host. Finally, there is one group where there is a dissociative relationship between genes associated with the pathogenic phenotype in A. fistulosum × A. cepa or A. porrum, but not both (Table 3). These results indicate that there might be selective pressure leading to the evolution/acquisition of host-specific virulence factors rather than those that may be generally useful for virulence like the HiVir cluster.

Five strains were chosen based on their species-specific pathogenicity on Allium hosts. The strains PNA 15-3, PANS 99-14 were pathogenic on A. porrum, but non-pathogenic on A. fistulosum × A. cepa and were all RSN-negative. The strains PNA 07-10, PNA 07-1, and PNA 05-1 were all pathogenic on A. fistulosum × A. cepa and RSN-positive, but non-pathogenic on A. porrum (Table 1). Gene presence and absence were compared manually among these strains. All strains possessed the HiVir cluster except for the strain PNA 15-3. The absence of this cluster is the likely cause for the strain’s inability to cause foliar lesions on A. fistulosum × A. cepa, and necrosis on onion scale.

The strain PNA 15-3; however, carries genes that indicate the presence of a type III secretion system, a virulence pathway that uncommon in P. ananatis (Supplementary Table 2). When comparing this type-III secretion system to those found in Kirzinger et al. (2015), it appears to show similarities with PSI-1b. Attempts to align this sequence to the type III secretion system found in P. stewartii subsp. indologenes indicated low sequence similarity. When PNA 15-3 was inoculated into tobacco leaf panels, no hypersensitive response was observed (Supplementary Figure 4). We found 43 genes proximal to each other surrounding the stcC gene in PNA 15-3. A total of 35 genes in the cluster were annotated as hypothetical proteins. The other eight genes were annotated as: oleC (olefin beta-lactone synthetase), gacA (response regulator GacA), mxiA (protein MxiA), hrcN (type III secretion ATP synthase HrcN), spaP (surface presentation of antigens protein SpaP), spaQ (surface presentation of antigens protein SpaQ), yscU (yop proteins translocation protein U), sctC (type 3 secretion system secretin), and dctD (C4-dicarboxylate transport transcriptional regulatory protein DctD). Further, we utilized NCBI database nucleotide BLAST to query the type III secretion system sequence at default values for the P. anantis taxid in the WGS database. We observed an 88% similarity with over 98% query coverage for PANS 99-23, PANS 99-26, PANS 200-1, PNA 86-1, UMFG54 (JACAFO010000015.1), NRRL B-14773 (JACEUA010000002.1), and DE0584 (VDNR01000019.1). Using NCBI database nucleotide blast at default values for the P. anantis taxid in the nr nucleotide collection, we observed an 88% identity with LCFJ-001 (CP066803.1) and FDAARGOS_680 (CP054912.1) chromosomal sequences.

The strains PNA 07-10, PNA 07-1, and PNA 05-1 shared several genes that do not occur in PNA 15-3, or PANS 99-14. These genes include the alt cluster, a cluster of 12 genes with three annotations (argT_3: lysine/arginine/ornithine-binding periplasmic protein, group_2282: ureidoglycolate lyase, and dapL: LL-diaminopimelate aminotransferase), a cluster of 12 genes (virB_2: virulence regulon transcriptional activator VirB, uspA_2: universal stress protein A, galE_2: UDP-glucose 4-epimerase, ybjJ_2: inner membrane protein YbjJ, nudK_3: GDP-mannose pyrophosphatase NudK, group_5271: phosphorylated carbohydrates phosphatase, mtnP: S-methyl-5′-thioadenosine phosphorylase, gph_2: phosphoglycolate phosphatase, arnB_2: UDP-4-amino-4-deoxy-L-arabinose—oxoglutarate aminotransferase, arnB_3: UDP-4-amino-4-deoxy-L-arabinose—oxoglutarate aminotransferase, perA: GDP-perosamine synthase, and iolG_5: inositol 2-dehydrogenase/D-chiro-inositol 3-dehydrogenase). These strains share 10 more genes in common, with 7 annotations (bepF: efflux pump periplasmic linker BepF, group_5443: adaptive-response sensory-kinase SasA, phoP_3: alkaline phosphatase synthesis transcriptional regulatory protein PhoP, parA_2: plasmid partition protein A, yedK_1: SOS response-associated protein YedK, ppaC: putative manganese-dependent inorganic pyrophosphatase, and crcB_2: putative fluoride ion transporter CrcB). Of the 60 genes shared between these strains, none of them are associated with the A. porrum disease phenotype and 58 genes are associated with the A. fistulosum × A. cepa disease phenotype. The two genes that were not associated with the A. fistulosum × A. cepa disease phenotype are annotated as hypothetical genes. Of these 60 genes, only one appears in the gene pair dissociation component 30 as the hypothetical gene “tar” dissociating with group_397 (Supplementary Table 2).

Apparent gene clusters shared by the A. porrum pathogenic PNA 15-3 and PANS 99-14 include a moderate gene cluster of 7 hypothetical genes and caf1M (chaperone protein Caf1M), another gene cluster of 7 hypothetical genes as well as 4 annotated genes (amiD_3: N-acetylmuramoyl-L-alanine amidase AmiD, yraI_2:putative fimbrial chaperone YraI, htrE_2: outer membrane usher protein HtrE, fimC: chaperone protein FimC). Half of a third cluster is shared between the two strains, with one hypothetical, with one hypothetical protein and three annotated proteins (group_5495: HTH-type transcriptional regulator PgrR, iolS_2: aldo-keto reductase IolS, ywrO_2: general stress protein 14). The final cluster consists of 10 genes without annotations, and group_7087 (replicative DNA helicase). Of the 39 shared genes, none appeared to be associated with the foliar pathogenic phenotype on A. porrum whereas 28 genes did appear to be associated with the foliar pathogenic phenotype for A. fistulosum × A. cepa. Furthermore, of these genes, none appeared in any gene-pair dissociation output. However, 19 of these genes appeared in the gene-pair association output within components 9 (N = 1), 11 (N = 5), 31 (N = 3), 37 (N = 2), 38 (N = 3), and 39 (N = 4; Supplementary Table 2).

The HiVir gene cluster found in the RSN-positive, A fistulosum × A. cepa pathogenic strains PNA 07-10, PNA 07-1, and PNA 05-1 did not contain any single nucleotide polymorphism (SNP) compared to that of the RSN-positive, wild type P. ananatis PNA 97-1 (Figure 7). However, three unique SNPs that resulted in missense mutations were identified in the RSN-negative, A. porrum pathogenic strain PANS 99-14. These mutations included alanine (A) to valine (V) change in amino acid position 7 in hvrA (pepM) gene, glutamine (Q) to lysine (K) change in amino acid position 352 in hvrB gene and, lysine (K) to arginine (R) change in amino acid position 11 in hvrH gene. These variant SNPs could be associated with disruption of the pantaphos pathway and loss of necrosis-associated phenotypes (Figure 7).

Comparative genome analysis focusing on only strains that were pathogenic on A. porrum, but non-pathogenic on A. fistulosum × A. cepa (PNA 15-3, PANS 99-14) vs. strains (PANS 99-11, PANS 99-12, PNA 06-4, PNA 99-9) that were pathogenic and highly aggressive on both hosts identified another pepM gene, which appeared to be a member of a secondary phosphonate biosynthetic cluster (Table 4). In 8/92 of P. ananatis strains (PANS 02-12, PANS 99-31, PANS 200-2, PANS 2-5, PANS 2-7, PANS 2-8, PANS 99-11, and PANS 99-12) there was a putative phosphonate biosynthetic cluster with 14 genes and a total length of approximately 19,000 bp. This cluster showed high sequence similarity to the pgb-cluster mentioned in Polidore et al. (2021), which was not responsible for generating onion-bulb rot symptoms. In our annotations, the left-flank of the cluster begins with a prophage integrase intS and is followed by the pepM phosphoenolpyruvate mutase (5′–3′ 900 bp long). The following cpdA is a 3′,5′-cyclic adenosine monophosphate phosphodiesterase cpdA (3′–5′ 771 bp). The third component of the cluster is fabG encoding 3-oxoacyl-[acyl-carrier-protein] reductase FabG3 (765 bp 3′–5′). The fourth gene is a phosphonopyruvate decarboxylase, aepY (1,176 bp, 3′–5′). The fifth gene of the cluster is phnW encoding for 2-aminoethylphosphonate-pyruvate transaminase (1,095, 5′–3′). The sixth component is asnB1 encoding putative asparagine synthetase [glutamine-hydrolyzing] (1,758 bp, 3′–5′). Following asnB2 is spsI1, encoding for Bifunctional IPC transferase and DIPP synthase (771 bp 5′–3′). The gene asd1 follows spsl1 that encodes aspartate-semialdehyde dehydrogenase (1,056 bp 5′–3′). The ninth component of the cluster is the MFS 1 transporter (1,227 bp, 5′–3′). The tenth component is glyA1, a serine hydroxymethyltransferase (1,359 bp, 5′–3′). The CDP-alcohol phosphatidyltransferase is the eleventh gene (600 bp 5′–3′). The twelfth gene is spsI2, a glucose-1-phosphate adenylyl/thymidylyltransferase Bifunctional IPC transferase and DIPP synthase (720 bp, 5′–3′). The thirteenth gene is aspC1, aspartate aminotransferase (1,173 bp, 5′–3′). The fourteenth gene is UDP-2,3-diacylglucosamine diphosphatase (762 bp, 3′–5′). Following the UDP-2,3-diacylglucosamine diphosphatase is another transposase. An interesting observation of this phosphonate cluster is the inclusion of phosphonopyruvate decarboxylase directly within the set of genes. This characteristic is unique when comparing it to the HiVir. The phosphonopyruvate decarboxylase has been described to play a critical role in the generation of phosphonates via the stabilization of the PEP mutase reaction (Supplementary Table 3; Kirzinger et al., 2015).

Of the 92 tested strains, 35 do not contain a complete alt gene cluster while the remaining 57 do possess the entire gene cluster. Of the 92 tested strains, 22 lacked a complete HiVir cluster and the remaining 70 strains possessed the entire gene cluster. Of the strains that produced foliar lesions, 45 had both alt and HiVir clusters, whereas 20 strains had only the HiVir gene cluster. Of the 92 tested strains, 3 strains only had alt whereas 7 strains lacked both gene clusters. Some foliar lesions were formed by the 7 strains (PANS 99-22, PANS 99-26, PANS 200-1, PANS 99-36, PNA 98-3, PNA 11-1, and PNA 15-3) that lacked both clusters, however the lesions varied in size and consistency between replicates. The strains that lacked both gene clusters were also RSN-negative. Of the 7 strains, showed some degree of foliar lesions on A. cepa × A. fistulosum; however, PNA 15-3 showed moderately aggressive lesion length on A. porrum but did not produce any foliar lesion on A. cepa × A. fistulosum. Of the 92 tested strains, only 8 strains contained the pgb cluster, while 84 strains lacked it (Figure 8).

When comparing the pgb cluster against the HiVir cluster, only the annotated pepM and MFS transporter are found to be common features (Table 4). In both clusters, the phosphoenolpyruvate mutase occurs first, with the MFS transporter being at the center of the cluster (5′–3′: 9th gene in pgb and 5′–3′: 9th in HiVir). There are no other shared annotated genes between the clusters. However, sequence alignment using Clustal Omega revealed 48.3% similarity between pepM from HiVir and pgb cluster. Similarly, the MFS transporters from HiVir and pgb clusters displayed 47.6% sequence similarity.

Based on the RSN assay the wild-type strain (PANS 02-18) and the single pepM mutant strain in the pgb cluster (ΔpepM_pgb) produced considerably large necrotic areas on red-onion scale compared to the single pepM mutant strain in the HiVir cluster (ΔpepM_HiVir; Figure 9A). Based on the seedling pathogenicity assay, the single pepM mutant strain in the HiVir cluster (ΔpepM_HiVir) had significantly lower necrotic lesion length on both Allium hosts compared to the wild-type strain (PANS 02-18) and the single pepM mutant strain in the pgb cluster (ΔpepM_pgb; Figures 9B,C). In both hosts, the deletion of pepM gene in the pgb cluster did not significantly affect the foliar lesion length compared to the wild-type strain (Figures 9D,E). While the double mutant strain where pepM genes were deleted in both the HiVir and the pgb clusters (ΔpepM_HiVirΔpepM_pgb) appears to have displayed a higher average lesion length than that of the single mutant strain (ΔpepM_HiVir) on A. porrum (Figures 9B,D). This result is surprising; however, we found no statistical significance between the mean lesion length of the ΔpepM_HiVirΔpepM_pgb vs. the ΔpepM_HiVir strains (Figure 9D). In A. fistulosum × A. cepa, the lesion lengths did not differ significantly between the double mutant strain and the single mutant strain (ΔpepM_HiVir).

Phenotyping of 92 P. ananatis strains displayed variability in the level of aggressiveness in both Allium spp. A considerably higher percentage of strains were either non-pathogenic or less aggressive on A. porrum (82.6%) than on A. fistulosum × A. cepa (65.1%). Also, a considerable percentage of strains were moderately or highly aggressive on A. fistulosum × A. cepa (34.9%) compared to A. porrum (17.4%). Interestingly, when a subset of less-aggressive strains (on both Allium spp.) was previously assayed on onion seedlings they were moderately-to-highly aggressive on onion (Stice et al., 2018; Agarwal et al., 2021). These observations potentially indicate that both A. porrum and A. fistulosum × A. cepa are inherently less susceptible to P. ananatis compared with the typical bulb onion. Also, when A. porrum and A. fistulosum × A. cepa were compared to each other, the former tends to show less severe symptoms compared to the later. However, we acknowledge that only one cultivar of each Allium spp. was evaluated, and it is possible that other cultivars or varieties of these hosts might show a range of susceptibility to P. ananatis. One aspect of these observations might be in part explained by the genetic nature of the A. fistulosum × A. cepa. cv. Guardsman itself. Being a hybrid between bunching onions and the typical bulb onion, it may be reasonable to expect A. fistulosum × A. cepa to be more susceptible to P. ananatis strains that were collected from symptomatic A. cepa tissues. However, without deeper genetic investigation of this hybrid we cannot predict for certain if any susceptibility-related genes or phenotypes were inherited. The P. ananatis culture collection used in this study also favors areas where onions are grown such as the Vidalia region and Tift County in Georgia, United States, which may provide a bias for aggressive bacterial strains on cultivars that are hybridized with A. cepa. Despite this, we recovered some strains that were more aggressive on A. porrum than A. fistulosum × A. cepa. Examples of these are PANS 99-11, PANS 99-12, and PNA 06-4 (Table 1).

Of the tested strains 56 showed a positive reaction to the RSN assay, while other strains (n = 36) did not. Most of the strains that were pathogenic on A. fistulosum × A. cepa were also able to cause necrosis on red onion scale. Based on the previous reports, HiVir is important for RSN-positive phenotype and foliar necrosis in onion, and it is likely that foliar lesions on A. fistulosum × A. cepa is also governed by the same gene cluster. Similarly, in A. porrum, a trend between RSN-positive phenotype and foliar pathogenicity was observed with majority of the strains. Further mutational analysis also indicated that HiVir is important for foliar pathogenicity on both A. fistulosum × A. cepa and A. porrum. However, we identified several strains that did not follow this trend. For example, while possessing a complete HiVir cluster and being RSN-positive, the PANS 19-8, and PANS 19-10 strains were unable to cause foliar lesions on A. porrum or A. fistulosum × A. cepa. This indicates that either mutation in their nucleotides/SNPs or other/alternative pathogenicity factors might be involved with this group of strains. Despite the importance of the HiVir cluster, there were 10 strains that seemed to cause foliar lesions despite lacking a complete HiVir cluster (PANS 99-22, PANS 99-26, PNA 18-9S, PANS 200-1, PANS 99-36, PNA 03-2, PNA 11-1, PNA 18-8S, PNA 98-3, and PNA 15-3). These observations indicate that other pathogenicity factors might be involved in these strains. Consistent with prior observations by Polidore et al. (2021), we also observed that the pgb cluster is not important for foliar pathogenicity in these Allium spp. The PANS 02-18 ΔpepM_pgb and wild-type strains displayed similar foliar pathogenicity phenotype in A. porrum or A. fistulosum × A. cepa.

In this study, we generated a pan-genome of 92 P. ananatis strains where we identified a conserved core genome of 2,914 genes, with a larger accessory genome of 9,196 genes. Earlier pan-genome reports identified similar values of core genes, with varying numbers of P. ananatis strains used for the analysis (De Maayer et al., 2014; Sheibani-Tezerji et al., 2015; Stice et al., 2018; Agarwal et al., 2021). In this work, however, we used a larger set of accessory genes compared to previously observed pan-genome study by Agarwal et al. (2021). We observed 6,808 cloud genes compared to 9,196 cloud genes in our current study. This discrepancy is like due to the use of larger number of diverse strains in this study compared to Agarwal et al. (2021) and is a reasonable increase for an open pangenome (Costa et al., 2020). Due to the cosmopolitan nature of P. ananatis and the extensive host range, it is entirely plausible that the bacterium would have a sizable pangenome when comparing populations from diverse hosts (De Maayer et al., 2014). It may be prudent to utilize a larger collection of strains from non-allium hosts for further pan-genomic assessments, where accessory genes may aid in resolving these strains further down beyond species taxonomic classification. In further attempts to discern the low statistical strength of the gene-associations, we noticed that several genes with high homology were assigned separate annotation tags. It is possible that these tags could artificially inflate the pangenome, leading to a weakened statistical association (Supplementary Figure 4). Even if the pangenome was artificially inflated, the GWAS results indicate significant association of genes that are known virulence factors (such as the HiVir cluster), as well as putative virulence factors that were found previously to be associated with foliar pathogenicity in A. cepa (Agarwal et al., 2021).

We identified 244 genes that were significantly associated with foliar pathogenicity in A. porrum whereas 836 total genes were associated with foliar pathogenicity in A. fistulosum × A. cepa, with 50 genes displaying significant agreement between naïve, Bonferroni, and BHC corrections for foliar pathogenicity in A. fistulosum × A. cepa. Among the genes associated with virulence in two hosts, 123 genes were shared. For both hosts, the HiVir cluster was found within the top-100 significantly associated genes. The occurrence of this cluster grants some additional credibility to other genes that show stronger significance with phenotypic association. Further investigation is required to determine the validity of these genes that are truly associated with the pathogenicity phenotype, or are errors from random sampling of the genome. In addition to this, most of the genes that occurred in the GWAS output apart from the HiVir cluster are annotated as hypothetical and would require further characterization to determine their relevance. Among these 123 genes, 48 of the top 50 statistically significant genes associated with the pathogenic phenotype in A. fistulosum × A. cepa were also associated with the pathogenic phenotype in A. porrum. While this alone is not enough to state their relevancy, it lends some credibility that these hypothetical genes may be useful as general virulence factors for the Alliums spp. and should undergo downstream mutational analysis to assess their functions.

Some of the genes with non-hypothetical annotations that were shared between the A. fistulosum × A. cepa and A. porrum include the entire HiVir gene cluster (some listed as hypothetical), pemK_2 (mRNA interferase), soj (sporulation initiation inhibitor protein), parM (Plasmid segregation protein), umuD (UmuD, translesion DNA polymerase subunit), tibC (glycosyltransferase), ycaD (uncharacterized MFS transporter), dadA (D-amino acid dehydrogenase 1), frbC (2-phosphonomethylmalate synthase), amiD (N-acetylmuramoyl-L-alanine amidase), and rfbB (dTDP-glucose 4,6-dehydratase). Upon manual investigation of the local pemK_2 region there seems to be a repeating pattern of seven genes, three flanking to the left, and four on the right. Utilizing BLAST for these gene sequences against the total gene sequences available for our P. ananatis strains shows that the pemK_2 region appears frequently throughout the entirety of the pangenome (Supplementary Table 5). The pemK gene is a known factor in toxin/antitoxin systems that are vital for bacterial competition and function (Lee et al., 2012; Klimina et al., 2013; Poluektova et al., 2017). Unfortunately, there is little information within the literature pertaining to the potential diversity and utility of P. ananatis toxin/antitoxin systems, including pemK. Without functional analysis, it is not possible to determine whether it is an Allium-specific virulence factor as opposed to a coincidental gene cluster, or if the annotation provided is correct. We also found an antitoxin gene, higB that was associated with A. fistulosum × A. cepa pathogenicity. Despite the lack of information, the region may be a valuable target for a toxin/antitoxin system within our P. ananatis strains from Georgia. Another annotated gene of interest includes tibA, an adhesin/invasion autotransporter. The sporulation initiation inhibitor protein, soj, is noted as possess a “centromere-like function involved in forespore chromosome partitioning inhibition of Spo0A activation” in Bacillus subtilis (Kunst et al., 1997). P. ananatis is not a spore forming bacteria; however, the inclusion of this gene with other genes that are similarly annotated for DNA manipulation may indicate that there is a requirement for maintaining genetic stability. The genes that constitute the “thiosulfinate tolerance; alt cluster,” only appeared in the A. fistulosum × A. cepa GWAS output with their Uniprot annotations of xerC (tyrosine recombinase XerC), altA/nemA (N-ethylmaleimide reductase), gor (glutathione reductase), altJ/osmC (peroxiredoxin OsmC), and altD/trxA (thioredoxin). Non-alt members of the OVRA region include rbsC (ribose import permease protein RbsC), rbsB (ribose import binding protein), and rbsA (ribose import ATP-binding).

Overall, GWAS was able to determine genes associated with foliar necrosis in A. fistulosum × A. cepa and A. porrum hosts. Follow-up experiments will test the validity of these hypothetical and annotated gene clusters for their relevance in pathogenicity and virulence in A. fistulosum × A. cepa and A. porrum hosts.

By comparing strains that were pathogenic on only A. fistulosum × A. cepa or A. porrum we hoped to significantly reduce the background noise (non-relevant genes from accessory) that may potentially result from the extensive P. ananatis pan-genome. Here we found several genes that belonged to strains that were only pathogenic to one host or the other, as well as a few interesting gene clusters.

One of the gene clusters of interest appears to be saturated with genes with annotated function of a typical type III secretion system. When using nucleotide blast in NCBI against the nr and whole genome sequence database, there were several other P. ananatis strains that shared a high consensus to the sequence (N = 11; in NCBI). Most of the strains in NCBI with an annotated type III secretion-system were isolated in Georgia, USA (N = 7; in NCBI) with two strains from onions and five strains from weeds. The remaining four strains with potential type III secretion-system were not isolated from in Georgia, USA. Further work is required to assess if the annotations are correct and investigate their role in onion pathogenicity.

Among the clusters shared by the A. fistulosum × A. cepa pathogenic strains PNA 07-10, PNA 07-1, and PNA 05-1, we found the alt cluster, two larger gene clusters (N = 12 genes, N = 22 genes), and two small gene clusters (N = 7 genes, N = 3 genes). Of the total 60 shared genes, none of them are on the A. porrum GWAS output, and 58 appeared in the A. fistulosum × A. cepa GWAS output. Only one gene appeared in the gene pair dissociation component 30 as the hypothetical gene tar that dissociated with group_397 (Supplementary File 4). Among the clusters shared by the A. porrum pathogenic strains; PNA 15-3 and PANS 99-14, we observed three gene clusters with 7, 8, and 11 genes each. Of these 39 genes, 28 were identified through GWAS as associated with the disease phenotype for A. fistulosum × A. cepa, but not associated with the disease phenotype in A. porrum (Supplementary Table 2). None of these genes appeared in the gene-pair dissociation output. However, 19 of these genes were found within components 9 (N = 1), 11 (N = 5), 31 (N = 3 genes), 37 (N = 2 genes), 38 (N = 3 genes), and 39 (N = 4 genes) in the gene-pair association analysis (Supplementary Table 2).

While the inclusion of the alt cluster in the A. fistulosum × A. cepa strains PNA 07-10, PNA 07-1, and PNA 05-1 is unsurprising, as they were isolated from symptomatic onions, its absence from PNA 15-3 is unexpected. The strain PNA 15-3 was isolated from symptomatic onion bulbs, and we would expect the presence of alt cluster as it aids in colonization of the onion bulb (Stice et al., 2020). The strain PANS 99-14 was isolated from an asymptomatic Digitaria spp. and may not need to rely on an alt cluster to survive in this environment. Despite this, both strains were able to generate a lesion on the A. porrum foliar tissue, and both strains failed to produce a positive result in the red-onion scale necrosis assay. These results indicate that some of the shared genes between these strains may aid in increased fitness in the A. porrum foliar environment that is not present in the A. cepa bulb tissue, or the A. fiustulosum × A. cepa foliar tissue. The large percentage of these shared genes (19/39) occur in the gene-pair association output, seems to suggest that these genes occur together at a higher frequency than others. Much like how the alt cluster provides protection to P. ananatis in the thiosulfinate-rich in bulb, it is possible that some of these clusters may provide protection to bacteria in the diverse Allium foliar environments.

We also aligned the HiVir gene clusters of A. fistulosum × A. cepa pathogenic (PNA 05-1, PNA 07-1 and PNA 07-10) and A. porrum pathogenic strains (PNA 15-3, and PANS 99-14) against each other and against the wild type P. ananatis strain PNA 97-1. No single nucleotide polymorphism (SNP) leading to missense mutation was identified in the RSN positive, A. fistulosum × A. cepa pathogenic strains (PNA 05-1, PNA 07-1 and PNA 07-10 strains). However, several missense mutations were present in the hvrA, hvrB and hvrH genes of the RSN negative, A. porrum pathogenic strain (PANS 99-14). TThese mutations were found only in hvr genes of PANS 99-14 but not in the hvr genes of the RSN positive strain (PNA 97-1). According to Metcalf and Van Der Donk (2009) and Polidore et al. (2021), hvrA and hvrB genes encode enzymes that are essential for the proposed phosphonate-toxin ‘pantaphos’ biosynthesis pathway. It is thus possible that the production of phosphonate toxin is compromised by these mutations. However, functional analysis needs to be conducted to confirm the impact of these mutations. In the case of RSN negative PNA 15-3, the strain lacked HiVir cluster but was still able to cause foliar lesions on A. porrum. It is possible that the pathogenicity of A. porrum in RSN negative PNA 15-3 and PANS 99-14 strains might possible be mediated by the genes other the genes in the HiVir cluster. Our results are limited by the number of strains used in the direct comparison; however, they do seem to indicate that there could possibly other Allium-specific virulence factors other than this phosphonate gene cluster (HiVir). These results are further visualized by aligning two figures (Figures 4, 8), where the pattern of the presence and absence of HiVir and alt gene clusters alone are not sufficient to map host pathogenicity (Supplementary Figure 5).

In this work, we utilized gene-pair coincidence as a supporting methodology to predict genes in P. ananatis are potentially relevant in Allium spp. pathogenicity. Hypothetically, genes that are important for survival in pathogenic bacteria should associate throughout a pangenome as their co-occurrence is beneficial for survival. Likewise, gene combinations that compromise survival in specific environments should dissociate with each other as natural selection selects against non-optimized populations. In this work, we hypothesized that the utilization of gene-pair coincidence should provide a phenotype-independent method of validation for the predicted genes from the pangenome via the phenotype-dependent GWAS methodology.

Our gene-pair association analysis generated 39 networks with a total of 165 individual genes, where two common genes were associated with the pathogenic phenotype in A. porrum and 45 common genes were associated with the pathogenic phenotype in A. fistulosum × A. cepa. In gene networks with genes that associated with the pathogenic phenotype (networks: 1, 3, 5, 11, 15, 16, 22, 24, 31, 32), the entire component is found in the GWAS result, indicating that these pathogenicity-associating genes are evolutionarily are potentially associated with each other via evolutionary process in our pangenome. Of these associated gene pair networks, none of the genes are associated with the disease phenotype for both hosts, only one or the other. Component 24 is unique in that only genes with the pathogenic phenotype in A. porrum were present. It is possible that these genes provide a unique advantage to overcome A. porrum host resistance. Unfortunately, both genes are annotated as hypothetical. The genes found in components 1, 3, 5, 11, 15, 16, 22, 24, 31, and 32 may be more useful for bacterial survival in A. fistulosum × A. cepa as opposed to being general virulence factors. Gene-pair association has provided support for further investigation of several potential gene groups that could potentially be harder to distinguish utilizing GWAS alone. Surprisingly, neither the HiVir gene cluster nor the alt cluster appear in the associative network. We would assume that both the gene clusters should co-occur together as they are relevant in bacterial pathogenicity and virulence in A. cepa and are quite prevalent in the pangenome. It could possibly due to limited number of genomes used in this study as compared to the number of strains/genomes needed to close the pan-genome. Another explanation could be the way the ROARY/Coinfinder organizes gene information and it is possible that the known virulence factors were omitted from pairwise analysis due to the lack of orthologous gene families. To determine if the issue was caused by noise because of genes not present in the gene clusters, we conducted analysis after removing individual genes from the ROARY csv file. However, this only strengthened resulting p-values, but did not improve the overall results. As such a larger sample size of strains with a comprehensive accessory genome could be included to better support GPC and GWAS results.

Our gene-pair dissociation analysis generated 50 gene-pair networks with a total of 255 genes, where only three genes were associated with the pathogenic phenotype in A. porrum and 22 genes associated with the pathogenic phenotype in A. fistulosum × A. cepa. Here, dissociation is dominated by networks where only a fraction of the dissociated genes also associated with the pathogenic phenotype in either of the Allium hosts. Dissociation networks 10, 17, 23, and 35 contained genes that associated with the pathogenic phenotype in A. fistulosum × A. cepa but dissociated with the genes that did not associate with the pathogenic phenotype in A. fistulosum × A. cepa. Dissociation network 37 is particularly interesting in that it showed dissociation between one hypothetical gene that associated with the pathogenic phenotype in A. porrum, and two genes that associated with the pathogenic phenotype in A. fistulosum × A. cepa. Whether these hypothetical genes have contrasting functions are yet to be evaluated; however, it is worth investigating if these genes play roles in host-pathogen-environment interactions. Dissociation network 21 is the only network where two genes that associated with the pathogenic phenotype in A. fistulosum × A. cepa dissociated with two genes that did not associate with the pathogenic phenotype in either hosts. These results are expected as genes that associate with the foliar pathogenic phenotype should dissociate with genes that do not associate with the same phenotype. Here we did not observe the HiVir or the alt genes dissociated with other genes, indicating that these gene clusters do not compete with each other in P. ananatis’ pangenome.

Again, these observations provide some level of confidence that the genes being predicted in the GWAS output are playing a role that enables them to be associated with the foliar pathogenicity phenotype on both Allium hosts. The diversity of GPC and their occurrence in the phenotype-dependent analysis enforce the assumption made previously that there could be several mechanisms of causing disease in Allium species other than phosphonate-based toxins.