Host Range Determinants of Pseudomonas savastanoi Pathovars of Woody Hosts Revealed by Comparative Genomics and Cross-Pathogenicity Tests

The study of host range determinants within the Pseudomonas syringae complex is gaining renewed attention due to its widespread distribution in non-agricultural environments, evidence of large variability in intra-pathovar host range, and the emergence of new epidemic diseases. This requires the establishment of appropriate model pathosystems facilitating integration of phenotypic, genomic and evolutionary data. Pseudomonas savastanoi pv. savastanoi is a model pathogen of the olive tree, and here we report a closed genome of strain NCPPB 3335, plus draft genome sequences of three strains isolated from oleander (pv. nerii), ash (pv. fraxini) and broom plants (pv. retacarpa). We then conducted a comparative genomic analysis of these four new genomes plus 16 publicly available genomes, representing 20 strains of these four P. savastanoi pathovars of woody hosts. Despite overlapping host ranges, cross-pathogenicity tests using four plant hosts clearly separated these pathovars and lead to pathovar reassignment of two strains. Critically, these functional assays were pivotal to reconcile phylogeny with host range and to define pathovar-specific genes repertoires. We report a pan-genome of 7,953 ortholog gene families and a total of 45 type III secretion system effector genes, including 24 core genes, four genes exclusive of pv. retacarpa and several genes encoding pathovar-specific truncations. Noticeably, the four pathovars corresponded with well-defined genetic lineages, with core genome phylogeny and hierarchical clustering of effector genes closely correlating with pathogenic specialization. Knot-inducing pathovars encode genes absent in the canker-inducing pv. fraxini, such as those related to indole acetic acid, cytokinins, rhizobitoxine, and a bacteriophytochrome. Other pathovar-exclusive genes encode type I, type II, type IV, and type VI secretion system proteins, the phytotoxine phevamine A, a siderophore, c-di-GMP-related proteins, methyl chemotaxis proteins, and a broad collection of transcriptional regulators and transporters of eight different superfamilies. Our combination of pathogenicity analyses and genomics tools allowed us to correctly assign strains to pathovars and to propose a repertoire of host range-related genes in the P. syringae complex.

The study of host range determinants within the Pseudomonas syringae complex is gaining renewed attention due to its widespread distribution in non-agricultural environments, evidence of large variability in intra-pathovar host range, and the emergence of new epidemic diseases. This requires the establishment of appropriate model pathosystems facilitating integration of phenotypic, genomic and evolutionary data. Pseudomonas savastanoi pv. savastanoi is a model pathogen of the olive tree, and here we report a closed genome of strain NCPPB 3335, plus draft genome sequences of three strains isolated from oleander (pv. nerii), ash (pv. fraxini) and broom plants (pv. retacarpa). We then conducted a comparative genomic analysis of these four new genomes plus 16 publicly available genomes, representing 20 strains of these four P. savastanoi pathovars of woody hosts. Despite overlapping host ranges, cross-pathogenicity tests using four plant hosts clearly separated these pathovars and lead to pathovar reassignment of two strains. Critically, these functional assays were pivotal to reconcile phylogeny with host range and to define pathovar-specific genes repertoires. We report a pan-genome of 7,953 ortholog gene families and a total of 45 type III secretion system effector genes, including 24 core genes, four genes exclusive of pv. retacarpa and several genes encoding pathovar-specific truncations. Noticeably, the four pathovars corresponded with well-defined genetic lineages, with core genome phylogeny and hierarchical clustering of effector genes closely correlating with pathogenic specialization. Knotinducing pathovars encode genes absent in the canker-inducing pv. fraxini, such as those related to indole acetic acid, cytokinins, rhizobitoxine, and a bacteriophytochrome. Other pathovar-exclusive genes encode type I, type II, type IV, and type VI secretion

INTRODUCTION
The Pseudomonas syringae species complex is considered the most relevant phytopathogenic bacterium worldwide, as a prominent research model and as an economically important phytopathogen (Mansfield et al., 2012). It comprises at least 13 phylogroups (PGs) encompassing 15 different Pseudomonas species associated with plants and the water cycle (Berge et al., 2014;Gomila et al., 2017). The complex is also divided at an infrasubspecific level into ca. 65 pathovars defined by characteristic plant host ranges (Bull and Koike, 2015).
The very broad host range of the P. syringae complex contrasts with the host specialization of most individual isolates from agricultural crops. Conversely, environmental isolates show a wider host range (Morris et al., 2013;Berge et al., 2014;Monteil et al., 2016). Moreover, a recent study showed that P. syringae strains form an overlapping continuum of host range potential and that pathovar denominations do not correspond to the P. syringae biology (Morris et al., 2019). Understanding this complexity at the different levels requires the establishment of appropriate pathosystems, facilitating integration of phenotypic, genomic, and evolutionary data. Consequently, the P. syringae complex emerged as one of the most relevant models to study host specificity in bacterial phytopathogens (Vinatzer et al., 2014;Baltrus et al., 2017;Dillon et al., 2019b).
Here, we performed comparative genomics analyses of 20 P. savastanoi strains isolated from olive, oleander, ash, and broom, covering all four well-defined P. savastanoi pathovars of woody hosts. First, we obtained the complete chromosome sequence of Psv NCPPB 3335 and draft genome sequences of Psn Psn23, Psr CECT 4861 and Psf NCPPB 1006. Second, cross-pathogenicity tests allowed correlating the results of the genomics analyses with host specificity. We thus delineated the core and pan-genome of the four P. savastanoi pathovars of woody hosts and identified the set of pathovar-and strain-specific gene complements. Our results provide a collection of putative virulence genes that can now be experimentally tested for their contribution to the definition of host specificity in knot-and canker-forming bacterial pathogens and their evolution into pathovars.

Cross-Pathogenicity Tests
The pathogenicity of P. savastanoi strains was analyzed on O. europaea plants derived from a seed germinated in vitro (originally collected from an "Arbequina" plant), N. oleander plants accession "pink" (single pink flowers) supplied by Viveros Guzmań (Maĺaga, Spain), and F. excelsior and R. sphaerocarpa plants native from Valladolid and supplied by Viveros Fuenteamarga (Valladolid, Spain). For inoculation, bacterial cells grown on LB plates for 48 h at 28°C were collected, washed twice with 10 mM MgCl 2 and adjusted to a final concentration of 10 8 colony-forming unit (CFU)/ml in the same buffer. Olive, oleander, and ash plants were wounded with a scalpel along the stem and 10 6 CFU were placed per wound (Penyalver et al., 2006;Peŕez-Martıńez et al., 2007). For broom plants, 10 7 CFU were injected in each wound using a syringe with needle. For the cross-pathogenicity tests shown in Figure 3 and Supplementary Figure S3 we used two plants per strain. The number of wound sites infected per plant varied between 5 and 10, depending on the size of the plant. For the cross-pathogenicity tests of Psv and Psn strains shown in Figure  4 we used three olive and three oleander plants per strain, with between 10 and 12 wounds sites infected per strain and plant. Wounds were protected with parafilm (Bemis, Neenah, WI, USA) for seven days to maintain humidity. Plants were kept in a greenhouse under natural photoperiod (15 h light/9 h dark) at room temperature (26°C day/18°C night). Symptoms were scored at 90 dpi and captured with a high-resolution digital camera (Nikon DXM 1200; Nikon Corporation, Tokyo, Japan), for olive, oleander, and ash, or with a stereoscopic microscope (Leica MZ FLIII; Leica Microsystems, Wetzlar, Germany), for broom plants. P. savastanoi cells were recovered from olive infected plants as follows. Three inoculated sites excised from two (Figure 3, Supplementary Figure S3) or three ( Figure 4) different plants were homogenized by mechanical disruption, and serial dilutions plated on LB containing Nf and Ch; bacterial populations were estimated by colony counts.

Genome Sequencing and Assembly
Genome sequencing and assembly of the Psv NCPPB 3335 chromosome was performed at BGI Tech Solutions Co., Ltd. (Hong Kong). A library of randomly sheared DNA fragments  (Table 3). c NCPPP 3335 data include both the complete sequence of its chromosome (accession number NZ_CP008742.1) and those of its three native plasmids (pPsv48A, FR820585.2; pPsv48B, FR820586.1 and; pPsv48C, FR820587.2).
(0.5-2.0 kb) was subjected to Illumina GA II (Solexa) sequencing. Reads (coverage, 200×) were qualitatively assessed before assembling with SOAP de novo (Li et al., 2008). Assembly of the NCPPB 3335 chromosome was aided by subtraction of the complete sequences of its three native plasmids (Bardaji et al., 2011). Primer walking and PCR amplification were used to fill the remaining gaps and solve misassembled regions. Sequencing of the Psn Psn23, Psr CECT 4861 and Psf NCPPB 1006 draft genomes was performed at Centro de Investigacioń Biomedica (CIBIR), La Rioja, Spain. Paired-end (100 cycles) sequencing used the Illumina GA ll platform (coverage, 250×). The genomes were assembled using CLC Genomics Workbench. The four genomes were automatically annotated upon submission to GenBank at National Centre for Biotechnology Information 1 (NCBI). The accession numbers of these four genomes are listed in Table 1.

Core-and Pan-Genome Analyses
Core− and pan-genome analyses were performed using BPGA v1.3 (Chaudhari et al., 2016) with assemblies downloaded from the NCBI 1 . In all comparative analyses, Psv NCPPB 3335 was used as reference with the sequence of the complete chromosome and of its three plasmids concatenated in a single file (5′chromosome-pPsv48A-pPsv48B-pPsv48C-3′). Orthologous genes were identified using the USEARCH algorithm (Edgar, 2010) with a threshold of 0.9 (90% blastp identity). To determine whether pan-genomes were opened or closed, we used the medians of the total number of genes found and the curves were then fitted to Heap's law model (Tettelin et al., 2008).

Phylogenetic Analyses
Phylogenetic relationships were predicted by analysis of core genome single nucleotide polymorphisms (SNPs) using Parsnp v1.2 and Gingr v1.2 (Treangen et al., 2014), and 100 bootstrap replicates; trees were visualized and manipulated using MEGA 7 (Kumar et al., 2016). The P. savastanoi phylogeny was also a n a l y z e d b y M L S A u s i n g 4 0 c o n c a t e n a t e d g e n e s (Supplementary Table S1) whose sequences are complete in all 20 genomes (80,902-80,949 nt per strain). Sequence alignment using MUSCLE within MEGA 7 (Kumar et al., 2016) was hand curated to eliminate ambiguities. The tree was constructed using RaxML 7.2.8 within Geneious 8.1.9 (Kearse et al., 2012), with 200 bootstraps. All trees were rooted using P. syringae pv. ciccaronei ICMP 5710 (assembly LJPY01) as outgroup.

Prediction of Secretion Systems and Effectors
Codification of structural and regulatory components of the type III (T3SS), type IV (T4SS) and type VI (T6SS)secretion systems was investigated using T346 Hunter (Martıńez-Garcıá et al., 2015a).
To search for known T3SS effectors (T3Es) we used PIFAR (Martıńez-Garcıá et al., 2016). T3Es sequences were manually examined to search for truncations, disruptions, and frameshifts. To identify novel hop genes, we used and ad-hoc pipeline considering two different criteria: the N-terminal sequence features, using EffectiveDB 3 (Eichinger et al., 2015), and presence of potential HrpL boxes 500 nucleotides upstream of the start codon, following previously described criteria (Fouts et al., 2002). Hierarchical clustering of P. savastanoi genomes based on their T3E content was performed using Morpheus 4 with a presence-absence matrix with values of "1" (full-length CDS), "0.5" (truncated/disrupted/frameshifted CDS) and "0" (absence of T3E). Values of "2" and "1.5" were given for T3E paralogs encoded as two full-length CDSs or as one full-length and one truncated CDS, respectively. The matrix was used to measure the distance between strains using the Euclidean distance method. Then, the tree was constructed using the complete linkage hierarchical clustering method.
The gene content of P. savastanoi pathovars was compared using Roary v3.12.0 (Page et al., 2015) as described above. The resulting presence/absence matrix was used to obtain the lists of pathovar-exclusive genes and pathovar-absent genes. Genes were first classified by gene ontology (GO) enrichment analysis using Sma3s.v2 (Muñoz-Meŕida et al., 2014) and then manually curated by blastx analysis (Supplementary Table S2).

General Genome Features
We reported the draft genome sequence of Psv NCPPB 3335 (Rodrıǵuez-Palenzuela et al., 2010) and the complete sequence of its three plasmids (Bardaji et al., 2011). To use this genome as a reference in this and future comparative genomic analyses, we obtained the complete sequence of the Psv NCPPB 3335 chromosome. This yielded a single contig of 6.02 Mb, encoding 5,850 genes and 5,165 protein-coding sequences. Likewise, we generated draft genome sequences of three strains whose pathogenicity in their host of isolation was previously reported, i.e. Psn Psn23, Psf NCPPB 1006 and Psr CECT 4861 (Table 1). Additionally, we analyzed the previously published draft genomes of 19 other strains of the four P. savastanoi pathovars of woody hosts ( Table 1). The relevant characteristics of all these genomes are shown in Table 1.

The Pan-Genome and Core Genome
Six of the genomes correspond to three strains from different collections, independently sequenced by us and others ( Table 1); we therefore selected those of the highest quality for subsequent genomic analyses. The resulting 20 genomes contained 114,601 genes, clustering into a hard pan-genome of 7,953 ortholog groups, of which 3,017 (38%) constitute the hard-core genome, 3,636 (46%) are accessory (encoded in 2-19 genomes), and the remaining 1,300 (16%) are strain-specific (Table 2, Figure 1). Nevertheless, there is a moderate strain-to-strain variation in the number of accessory genes (1,857 to 2,388) and strain-specific genes (14 to 236) ( Figure 1B). We also estimated that only 19 additional ortholog groups are part of the P. savastanoi soft-core genome, i.e. 3,036 groups are present in at least 19 of the 20 genomes (95%).
Using the Heap's law model, we estimated a fitting parameter (g) of 0.15 for all 20 genomes ( Table 2), just above the threshold of g=0 distinguishing open (g > 0) from closed (g < 0) genomes. Therefore, sampling more P. savastanoi genomes will likely not result in a significant increase of the pan-genome size nor in a substantial reduction of the core genome, because their accumulation curves tend towards an asymptote ( Figure 1A). A heat map representing presence and absence of individual genes on each strain is shown in Figure 2.

Phylogeny of P. savastanoi Pathovars
We analyzed the phylogenetic relationships among the 20 P. savastanoi strains ( Table 1) using SNPs from a core genome alignment. This tree ( Figure 2) shows a very similar topology to a MLSA-based tree constructed using 40 concatenated genes (Supplementary Table S1, Supplementary Figure S2).
We could identify five well-differentiated clades (Figure 2), two of which correspond to pathovars Psr and Psf. However, pathovar Psv appears not to be monophyletic being distributed into two different clades. Two strains of Psv (NCPPB 3335-pABC and ICMP 1411) clustered closely to the Psr strains (Psv-I cluster), well-separated from the remaining Psv strains (Psv-II cluster). This last group also includes two oleander strains (ICMP 13786 and 0485_9), whereas the remaining Psn strains clustered together in a separate clade.
Psr CECT 4861 induced well-defined knots at most inoculation sites only on broom plants ( Figure 3A, Supplementary Figure S3). Ash and olive plants showed tissue swelling at most inoculation sites (80% and 60%, respectively), while infected oleander plants remained symptomless ( Figure  3A, Supplementary Figure S3). These results agree with the narrow host range reported for P. savastanoi isolates from broom (Álvarez et al., 1998).
On ash plants, Psf strains NCPPB 1006 and CFBP 5062 induced excrescences (necrotic swellings) at most inoculation points ( Figure 3A, Supplementary Figure S3). Excrescences also appeared in up to a third of infected sites on olive plants, supporting the reported ability of Psf to also induce symptoms in this host (Janse, 1981;Janse, 1982;Iacobellis et al., 1998). Inoculated oleander plants were symptomless.
All four Psv strains induced visible symptoms on olive and ash plants, although they were variable and strain dependent ( Figure 3A, Supplementary Figure S3). For instance, Psv NCPPB 3335 and Psv PseNe107 generated knots on olive plants, which were more necrotic on wounds inoculated with the later, while Psv CFBP 1670 and Psv DAPP-PG722 primarily induced swelling of the tissues surrounding the wounds. This agrees with the reported strain-and cultivar-dependent severity of symptoms in inoculated olive plants (Penyalver et al., 2006;Peŕez-Martıńez et al., 2007). Ash plants primarily developed Analyses of pv. nerii and pv. savastanoi considering the oleander strains 0485_9 and ICMP 13786 as belonging to pv. nerii, as they were originally classified. The other analyses were done including these two strains in pv. savastanoi, as indicated by the cross-pathogenicity tests.  (Janse, 1981). Nevertheless, both Psn CFBP 5067 and Psn23 appear to induce systemic infections because inoculated oleanders showed secondary symptoms, such as flower bud swellings, leaf curling, and tiny knots on the leaves (see Supplementary Figure S4 for CFBP 5067). Broom plants infected with Psn showed no symptoms ( Figure 3A, Supplementary Figure S3).
All the strains reached the highest populations densities at 90 days post-inoculation (dpi) in their host of isolation, with A B FIGURE 1 | (A) Core-and pan-genome of P. savastanoi strains. Boxes indicate changes in number of gene families to the number of genes added sequentially, with median values denoted with a horizontal black line and standard deviation with vertical bars. The pie chart shows the frequency distribution of ortholog groups of genes. (B) Strain-specific P. savastanoi genes. Number of genes in the core-genome (center), accessory genes (petals) and strain-specific genes (in brackets), calculated using BPGA. Psv, Psn, Psf, and Psr, P. savastanoi pathovars savastanoi, nerii, fraxini, and retacarpa, respectively. Asterisks indicate oleander strains included in Psv (see Figure 4). medium to low populations in the other hosts. No viable counts were recovered from symptomless broom plants infected with Psf NCPPB 1006 ( Figure 3B).
In summary, our results agree with the host range previously reported for the four well-established P. savastanoi pathovars of woody hosts. Additionally, we also found that some Psv strains can occasionally induce the formation of slightly swollen lesions on broom, and that strain Psr CECT 4861 also induces the formation of slight symptoms at some inoculation sites not only on broom, but also on olive and ash ( Figure 3A, Supplementary Figure S3).

Phylogeny of P. savastanoi Strains Isolated From Olive and Oleander Reflects Differential Pathogenicity in These Two Hosts
Psv strains did not cluster in a monophyletic lineage and two oleander strains clustered with four olive strains (Figure 2,  Supplementary Figure S2). To confirm their pathovar assignation, we therefore examined the pathogenicity in olive and oleander of the six olive and oleander strains not included in the cross-pathogenicity tests shown in Figure 3.
Psn strains ICMP 16944 and ICMP 13781 induced the largest knots in olive plants, but their bacterial populations were the lowest ( Figure 4). Additionally, they were the only strains inducing knots or swellings on oleander plants, were they reached very high populations. Conversely, oleander strains ICMP 13786 and 0485_9 induced symptoms in olive plants, were they reached the highest bacterial populations, but not in oleander. Thus, these two strains resemble the typical behavior of Psv strains and, although isolated from oleander, they belong to pathovar savastanoi.

Secretion Systems
The distribution of the T3SS, T4SS, and T6SS determinants was investigated using the web-based tool T346Hunter (Martıńez-Garcıá et al., 2015a).
As in Psv NCPPB 3335 (Rodrıǵuez-Palenzuela et al., 2010), all the genomes encoded two different T3SS clusters: the tripartite pathogenicity island, the most prevalent T3SS within the P. syringae complex (Baltrus et al., 2017;Dillon et al., 2019b), and an additional T3SS resembling that found in Rhizobium species ( Figure 5) (Joardar et al., 2005;Martıńez-Garcıá et al., 2015b). Although most of the corresponding deduced products showed 100% identity among genomes, we also found amino acid variations characteristic of each of the five phylogenetic groups ( Table 4).
The complete set of T4SS-A genes (virB1-virB11 and virD4) carried by plasmid pPsv48B (Bardaji et al., 2011) was found in 12 genomes ( Figure 5). Additionally, all 20 strains carried incomplete sets of T4SS-B genes (tra and/or trb genes) and T4SS-C genes (tfc genes) ( Figure 5). The T4SS-C is involved in propagation of genomic islands (Juhas, 2015), and has not been previously described in the P. syringae complex. As described for Psv NCPPB 3335 (Rodrıǵuez-Palenzuela et al., 2010), P. savastanoi strains encode two different T6SS clusters, one highly similar to that found in P. syringae pv. phaseolicola 1448A and a second homologous to the T6SS hcp1 cluster of P. syringae pv. tomato DC3000. Considering the composition of both T6SS clusters and additional chromosomal genes, all 20 genomes contains 12 of the 13 genes constituting the T6SS core (Boyer et al., 2009), only lacking that encoding the forkheadassociated (FHA) domain-containing protein VCA0112 ( Figure 5).  Figure 1; PT, pathotype strain.

Distribution of T3SS Effectors
The P. savastanoi pan-genome encodes 43 known T3Es ( Figure  6), including truncated genes, distributed in 37 T3E families from the Hop database 6 and 33 of the new merged families (Dillon et al., 2019a). T3E truncations were analyzed and assigned to groups according to their extent (Figure 7). All the strains contained the same truncations of hopAZ1, hopAA1, and hopM1. Conversely, T3E allelic variants encoded by specific pathovars, and therefore possibly involved in host specificity, include those only found in all Psn (avrPto1, hopAT1, and hopA2) or Psr strains (hopAF1-2 and hopAO1). Additionally, while all Psv and Psn strains encode a full-length coding sequence (CDS) for avrRpm2, Psr, and Psf strains either encode a truncated version of this gene or do not contain it. Specific truncations affecting single strains were also identified for hopA2, hopAF1-1, hopAO1, and hopW1 (Figure 7).
We found two novel T3E candidates not homologous to known T3Es (Table 5) but preceded by a hrp box promoter, characteristic of genes regulated by HrpL (Fouts et al., 2002). Additionally, their deduced products meet at least two of the three features of N-terminal secretion signals found in Hop 6 http://www.pseudomonas-syringae.org/ A B FIGURE 3 | Cross pathogenicity tests of P. savastanoi strains. (A) Symptoms generated by P. savastanoi strains at 90 days post-inoculation. Control (-), negative control, plants inoculated with 10 mM MgCl 2 . (B) Populations of P. savastanoi strains in different hosts. Bars represent mean populations at 90 days post-inoculation from three biological replicates, with standard error; different letters indicate means that are significantly different using ANOVA test followed by the Bonferroni t-test (p < 0.05). Psv ICMP 4352 was obtained from the CFBP collection as strain CFBP 1670. Psv, Psn, Psf and Psr, P. savastanoi pathovars savastanoi, nerii, fraxini, and retacarpa, respectively.
proteins . The gene coding for WP_032074452 was found in all 20 P. savastanoi genomes, whereas that for WP_057446723 was exclusively found in Psr genomes.
Hierarchical clustering based on T3E content identified three main clusters composed of i) all Psf genomes (Psf cluster), ii) all Psr genomes (Psr cluster) and iii) all Psv and Psn strains (Psv-Psn cluster) ( Figure 6). The Psv-Psn cluster was divided into two sub-clusters, one comprising all eight Psv strains and the other the four Psn strains ( Figure 6). Thus, T3E clustering of P. savastanoi genomes is similar to the core genome phylogeny ( Figure 2) indicating that the four pathovars contain characteristic sets of T3Es, which are thus likely contributing to define their pathogenicity profile.
We found 21 core P. savastanoi T3E genes, encoded as fulllength ORFs or truncated versions in all the genomes. Three additional T3Es (hopAF1-2, hopAO2, and hopD1) were found in 19 of the 20 genomes analyzed (95%) and were included in the soft-core effectorome. Four T3E genes are candidates to contribute to host range definition: the Psr-exclusive avrRps4, hopF1, and hopZ5 genes, and the hopG1 gene, only absent in all three Psr strains ( Figure 6). Known T3E encoded in single pathovars but not in all strains, thus likely not having a relevant role in host specificity, are: hopBK1 (Psv strains), hopZ1, hopH3, hopZ4, hopAW1, avrA1, avrB2, and hopAV1 (Psf strains).

Identification of Pathovar-Specific Genes
We identified several variable regions (VR) longer than 3.9 kb that were exclusively present or absent in distinct monophyletic groups (Figure 8). In general, the patterns of presence/absence of these VRs are congruent with the organismal phylogeny ( Figure  2) and suggest their acquisition before pathovar differentiation followed by further vertical inheritance. Six VRs were of interest for being chromosomal and exclusively absent from all Psf strains (Figure 8, Supplementary Figure S1). VR3, VR4, VR5, and VR6 are enriched in hypothetical proteins (HPs) whereas VR1 and VR2 contain genes encoding putative virulence factors. VR1 encodes the iaaMH operon, responsible for the biosynthesis of indole-3-acetic acid (IAA) via indole-3-acetamide, and two genes involved in rhizobitoxine biosynthesis (rtxA and rtxC), a phytotoxin modulating plant-microbe interactions by inhibition of ethylene synthesis (Sugawara et al., 2006). In turn, VR2 encodes the type II secretion system proteins GspH and GspI and a pectate lyase (Figure 8, Supplementary Figure S1). With the exception of VR2, the other five VRs were predicted as genomic islands, being also often bordered by mobile genetic elements.
The core and pan-genome profiles of the 20 strains were also analyzed at the pathovar level using BPGA ( Table 2). Psf showed the largest pathovar soft-core genome size, with a nearly closed pan-genome (g=0.065), suggesting that the gene content of Psf strains is more homogeneous than that of the other three pathovars. We also conducted this analysis considering the original assignation of strains 0485_9 and ICMP 13786 to pv. nerii ( Table 2), which completely prevented the identification of exclusively present/absent genes in both Psn and Psv. This illustrates the confounding effects that might originate from an improper identification of the strains under analysis, and reinforces the need for appropriate phenotypic analyses to support conclusions derived exclusively from genomic data. We explored the gene repertoire specific to pathovars using two complementary bioinformatics tools, Roary (Page et al., 2015) (Supplementary Table S2) and PIFAR ( Figure 5) (Martıńez-Garcıá et al., 2016). A total of 123 pathovarexclusive genes and 69 pathovar-absent genes were identified and classified into eight categories, among which HPs and metabolism-related genes were the most common ( Table 2,  Supplementary Table S2). Many of the pathovar-exclusive and pathovar-absent genes are likely carried by plasmids and other mobile genetic elements (Figure 8, Supplementary Figure S1), because we found a fair number of genes characteristic of these elements, such as genes involved in DNA replication, recombination, mutation, and repair, and toxin/antitoxin systems (Bardaji et al., 2019). Although around 24% of the pathovar-exclusive/absent genes encoded HPs, other pathovarexclusive genes encode proteins with a possible role in virulence and/or host specificity, such as 11 transcription factors, three cdi-GMP-related proteins, two chemotaxis-related proteins, 18 transporters belonging to eight superfamilies, the Psn-exclusive secretion systems chaperones PapD (type I secretion system; T1SS) and ShcF (T3SS), four Psr-exclusive T3Es (see also Figure  6 and Table 5) and the ferric-pseudobactin receptor PupB, also exclusive of Psr strains.

DISCUSSION
Genomics-based methods offer unprecedented resolution to analyze strain diversity and evolutionary relationships within the P. syringae complex, besides identifying candidate genes involved in pathogenicity and host range definition Baltrus et al., 2011;Bartoli et al., 2015;Monteil et al., 2016;Baltrus et al., 2017;Dillon et al., 2019a;Dillon et al., 2019b). However, understanding the biological meaning of the observed genomic differences depends on the availability of experimental data confirming specific phenotypes, such as virulence and host range (Baltrus and Orth, 2018). Specifically, testing more than one hundred P. syringae strains on different hosts showed that the host range is an overlapping continuum, with closely related strains exhibiting from extremely narrow to very broad host ranges (Morris et al., 2019). Our genomic analyses of 20 strains representing all four well-established P. savastanoi pathovars of woody hosts were combined with cross-pathogenicity tests, supporting the identification of candidate genes involved in host range definition. In fact, our combination of bioinformatics and functional analyses was critical to guarantee a correct pathovar assignation, ensuring the correct interpretation of the comparative data.
The P. savastanoi Pathovars Pan-Genome The P. savastanoi pan-genome (7,953 ortholog families for 20 genomes) ( Table 2) is approximately 10 and 3.5 times smaller than those of the P. syringae complex (391 genomes) and PG3 (143 genomes), respectively (Dillon et al., 2019b). Additionally, the pan-genomes of the P. syringae complex and PG3 are still open (Dillon et al., 2019b), whereas the pan-genome for each P. savastanoi pathovar is nearly closed ( Table 2). The large homogeneity among P. savastanoi genomes results from a low content of pathovar-and strain-specific genes (16% strainspecific genes vs. ca. 60% in the P. syringae complex; Figure  1), probably reflecting close phylogenetic relationships, a large repertoire of shared niche-specific genes and, perhaps, a natural inability to acquire foreign DNA (Peŕez-Martıńez et al., 2007). However, it is possible that the analyzed genomes may not represent the true diversity of the species and its four pathovars of woody hosts. Thus, the analysis of additional P. savastanoi genomes from strains isolated from woody host, as well as from the surface of herbaceous plant hosts or from nonagricultural environments, could contribute to the genetic diversity of this pathogen, resulting in more open pangenomes and reducing the number of host range determinants identified in this study. In fact, other authors have reported the existence of further genotypical and phenotypical variability among P. savastanoi strains, such as the ability to produce levan or an atypical host range (Marchi et al., 2005;Sisto et al., 2007;Cinelli et al., 2014;Moretti et al., 2017;Morris et al., 2019).

Correlation Between Phylogeny and Pathogenicity
Seven of the strains analyzed here were repeatedly shown to cause distinct syndromes in their host of isolation (Psv NCPPB 3335, DAPP-PG722, and ICMP 1411;Psn Psn23;Psf NCPPB 1006 and CFBP 5062;Psr CECT 4861). However, pathogenicity data were not available for the remaining 13 strains, which were assigned to a P. savastanoi pathovar (Tables 1, 3) based on their host of isolation and, in some cases, on biochemical tests.
The cross-pathogenicity tests for nine P. savastanoi strains in olive, oleander, ash, and broom plants (Figure 3, Supplementary Figure S3) agree with the host range reported for the four wellestablished P. savastanoi pathovars of woody hosts. This supports using their genomes for comparative genomics analysis to identify candidate genes involved in host range definition. Psv, Psn, and Psf strains showed differences in virulence and in the ability to grow and survive in plant tissues, agreeing with previously reported cultivar-dependent The Pfam database (https://pfam.xfam.org/). c hrp value, counter variable used to identify hrp boxes in the promoter regions, 500 nucleotides before the start of the gene. d Prediction of N-terminal secretion signal in Hop proteins. A, high percentage of serine in the 50 N-terminal amino acids; B, isoleucine, leucine, valine, alanine, or proline in the third or fourth N-terminal amino acids; C, absence of acidic residues within the 12 N-terminal amino acids . Asterisks indicate oleander strains included in pathovar Psv (see Figure 4). e Psv, Psn, Psf and Psr, P. savastanoi pathovars savastanoi, nerii, fraxini and retacarpa, respectively; + = presence, -= absence of a protein ortholog in the genomes of all strains included in the corresponding pathovar analyzed in this study Table 1. virulence variability (Penyalver et al., 2000;Moretti et al., 2017). Additionally, we show that symptom severity is not always correlated with population levels in planta. In general, P. savastanoi strains showed the highest population level in their host of isolation ( Figure 3B), regardless of their ability to also induce visible symptoms on other hosts. Nevertheless, maintenance of reduced bacterial density in non-host plants could favor the appearance of better-adapted variants and the development of emerging diseases in new hosts.
Agreeing with previous data (Gomila et al., 2017;Dillon et al., 2019b), Psf, Psn, Psr, and Psv cluster in a discrete genomic branch in the SNP ML tree (Figure 2) with non-tumorigenic pathovars as the phylogenetically closest relatives (Dillon et al., 2019b). This suggests that Psf, Psn, Psr, and Psv originate from a common ancestor, which likely acquired the ability to induce tumors in plants. These four pathovars also corresponded with well-defined genetic lineages indicating distinct evolutionary events leading to host specialization, although it is necessary to make two considerations. First, the Psv genomes are distributed in two separated clades, clustering with either Psr or Psn genomes (Figure 2). Although intriguing, the occurrence of different genetic lineages for a given pathovar is not new nor surprising (Gironde and Manceau, 2012;Martıń-Sanz et al., 2012;Martıń-Sanz et al., 2013;Hulin et al., 2018), because pathovars are defined for causing distinctive disease syndromes and not necessarily for their genetic relatedness (Young et al., 1992). Nevertheless, olive is a frequent host for the four pathovars and, since strains of P. savastanoi can cause knots in diverse species (Caballo-Ponce et al., 2017a), it is also possible that the two Psv genetic lineages show differing host ranges. In fact, pathogenicity assays on 15 herbaceous and woody hosts showed that the host range of Psv PseNe107 (Psv-II) is broader than that of Psv NCPPB 3335 (Psv-I) (Morris et al., 2019). A second relevant consideration is the clustering with Psv strains of two putative Psn strains (Figure 2). These two strains were more like Psv than Psn for several characteristics, such as their virulence gene repertoire ( Figure 5) or their pattern of T3E genes ( Figure 6). Nevertheless, the ability of strains 0485_9Ĕ and ICMP 13786 to cause knots in olive but not in oleander, and their better survival in olive than in oleander tissues (Figure 4) warrants their unequivocal classification as Psv strains (Janse, 1982).

T3SS and Effectors
All 20 genomes encode a canonical T-PAI T3SS cluster ( Figure  5, Table 4). Since it is essential for knot formation and the induction of a visible hypersensitive response (HR) by Psv and Psn (Caballo-Ponce et al., 2017a), and given its high sequence FIGURE 8 | Genomic comparison of Pseudomonas savastanoi pathovars. Comparison made using the GView Pangenome analysis tool. Black arrows, selected variable regions (VR) that could be pathovar-specific (see Supplementary Figure S1). Strain abbreviations as in Figure 1, indicating the phylogenetic grouping from Figure 2 to the right of the legend. and functional conservation ( Table 4) (Büttner and He, 2009), we could expect a similar role for Psf and Psr. The number of T3Es in P. syringae complex strains varies from four to nearly 50 (Dillon et al., 2019a), and the average in woody hosts strains is 29 (Nowell et al., 2016). Similarly, P. savastanoi strains encode 28 to 35 T3E genes ( Figure 6). As expected, the only three core T3E genes of the P. syringae complex primary PGs (avrE1, hopAA1, and hopAJ2) and the broadly distributed hopM1 (Dillon et al., 2019b), are also P. savastanoi core T3Es. P. savastanoi strains also encoded several T3E genes associated with pathogenicity on woody hosts, i.e. the core members hopAO1, hopBL1, and hopBL2 (Matas et al., 2014;Nowell et al., 2016;Castañeda-Ojeda et al., 2017), as well as hopAY1, hopH3, hopZ5, and hopAF1-1 ( Figure 6). Additionally, all 20 genomes encode hopAN1 and hopJ1, present in all 11 P. syringae complex PGs, but removed as T3Es or reclassified as T3SS helpers (Dillon et al., 2019b).
Truncations of T3E genes are common in P. syringae (Baltrus et al., 2011;Dillon et al., 2019a), and we considered them here ( Figure 7) because some of their products are translocated into plant cells and interfere with plant defenses (Zumaquero et al., 2010;Matas et al., 2014). T3Es are considered the primary determinants of host specificity in the P. syringae complex at both the pathovar-species and race-cultivar levels (Collmer et al., 2000;Fouts et al., 2003;Cunnac et al., 2009;Lindeberg et al., 2009;Mansfield, 2009). We therefore searched for exclusive T3E genes present/absent in specific pathovars. All three Psr genomes encode avrRps4, hopF1, hopZ5, and the putative novel T3E gene B7R56_RS15315. Additionally, they lack hopG1, a gene found in the remaining 17 strains. These four known T3E genes were shown to contribute to the virulence of P. syringae strains but also function as avirulence factors in soybean and Arabidopsis (avrRps4), bean (hopF1), Arabidopsis (hopZ5) and cherry trees (hopG1) (Hinsch and Staskawicz, 1996;Tsiamis et al., 2000;Jayaraman et al., 2017;Hulin et al., 2018). Thus, avrRps4, hopF1, and hopZ5 might contribute to Psr virulence in broom plants and/or restrict the infection of ash, oleander, and olive. Furthermore, Psr strains lack hopG1 and hopAF1-2, or encode a truncated HopAF1-2 protein lacking its catalytic domain ( Figure 7). Therefore, the absence of these T3Es could be related to their recognition by broom resistance genes.
Although we did not find Psn-exclusive T3E genes, Psn strains encode a specific 3′ truncation of avrPto1 (Figure 7). Considering that point mutations affecting the C-terminal region of AvrPto abolish avirulence in tobacco (Shan et al., 2000), the specific Psn truncation might prevent avirulence in oleander. Other Psn-exclusive allelic variants possibly involved in host range definition are those of hopAT1 and hopA2 (Figure 7). Additionally, Psf and Psr strains do not encode a chimera of avrRpm2, with an N-terminal region homologous to that of the HopF4 family (Lo et al., 2017).

Pathovar-Specific Genes
Using three complementary strategies, we identified a repertoire of pathovar-specific genes potentially associated with host specificity (Figures 5, 8, and Supplementary Table S2). GView identified several genomic regions exclusively present/absent in all strains from monophyletic branches of the P. savastanoi phylogeny, most of which correspond to single pathovars ( Figure 8). Since this strategy is not suitable for the identification of individual genes encoded by specific pathovars, we used the additional bioinformatics tools Roary (Supplementary Table S2) and PIFAR ( Figure 5). Nevertheless, many of the pathovar-exclusive and pathovar-absent genes identified by Roary were also detected using GView or PIFAR, increasing the reliability of the results obtained.
As expected, all non-tumorigenic Psf strains lack the iaaMH operon ( Figure 5), required for the biosynthesis of IAA and knot formation by Psv and Psn strains (Comai and Kosuge, 1983;Iacobellis et al., 1998;Aragón et al., 2014). Tumorigenic P. savastanoi pathovars, encode a copy of this operon in VR1, a genomic island also encoding the rhizobitoxine operon and two glutamine biosynthesis-related genes that likely supply glutamine for rhizobitoxine biosynthesis (Sugawara et al., 2006). In Bradyrhizobium elkanii, this enol-ether amino acid toxin enhances nodulation by inhibiting ethylene biosynthesis in host roots. Remarkably, production of IAA in B. elkanii is restricted to strains also synthetizing rhizobitoxine. Therefore, rhizobitoxine might also be involved in knot formation.
All pathovars, except Psn, contained genes for the biosynthesis of the phytotoxin phevamine A. This molecule, composed of a modified spermidine, L-phenylalanine, and L-valine, contributes to the virulence of P. syringae DC3000 on Arabidopsis suppressing immune responses induced by flagellin (O'Neill et al., 2018). We did not find genes for the biosynthesis of other phytotoxins involved in virulence in woody hosts (Hulin et al., 2018).
Gene ptz, involved in cytokinins production, was neither found in Psf strains nor in the genomes of Psn ICMP 13781 and the oleander strains here reassigned into Psv (0485_9 and ICMP 13786) ( Figure 5). Gene ptz is plasmid encoded in most Psv and Psn strains, and is absent in some Psn strains (Peŕez-Martıńez et al., 2008). In addition, failure to produce cytokinins is correlated to symptoms attenuation, both on olive and oleander (Surico et al., 1985). However, Psn ICMP 13781 consistently induced tumors in olive as large as those produced by other Psn strains carrying gene ptz (Figure 4), suggesting that induction of full symptoms by Psn on olive is also dependent on other virulence determinants, as reported for Psv strains (Caballo-Ponce et al., 2017a).
In conclusion, besides T3SS genes whose role in host specificity has already been demonstrated in the P. syringae complex, our results suggest the existence of other putative virulence factors that could also contribute to host range definition, such as proteins secreted by the T2SS and T6SS, genes involved in the biosynthesis of the phytotoxins rhizobitoxine and phevamine A, a large number of transcriptional regulators, proteins involved in the metabolism of c-di-GMP, signal-transduction proteins and transporters of eight different superfamilies.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
JM, ST, CM, PR-P, and CR planned and designed the research and analyzed and interpreted the data. AM-P, AP and EC-P performed experiments. AM-P, AP, JM, PR-P, and CR carried out bioinformatics analyses. AM-P, AP, JM, and CR designed and prepared Figures and Tables. AM-P, AP, EC-P, JM, and CR wrote the manuscript.