Background

Host specialization is a hallmark of numerous plant pathogens including bacteria, fungi, oomycetes and viruses. Yet, the molecular and evolutionary bases of host specificity are poorly understood. In some cases, pathological convergence is observed for individuals belonging to distant phylogenetic clades. This is the case for Xanthomonas strains responsible for common bacterial blight of bean, spread across four genetic lineages. All the strains from these four lineages converged for pathogenicity on common bean, implying possible gene convergences and/or sharing of a common arsenal of genes conferring the ability to infect common bean.

Results

To search for genes involved in common bean specificity, we used a combination of whole-genome analyses without a priori, including a genome scan based on k-mer search. Analysis of 72 genomes from a collection of Xanthomonas pathovars unveiled 115 genes bearing DNA sequences specific to strains responsible for common bacterial blight, including 20 genes located on a plasmid. Of these 115 genes, 88 were involved in successive events of horizontal gene transfers among the four genetic lineages, and 44 contained nonsynonymous polymorphisms unique to the causal agents of common bacterial blight.

Conclusions

Our study revealed that host specificity of common bacterial blight agents is associated with a combination of horizontal transfers of genes, and highlights the role of plasmids in these horizontal transfers.

In nature, most pathogens are generalists, meaning that they are able to infect multiple hosts, while other pathogens are specialists, meaning that they are highly adapted to a single or few host species [1]. For plant pathogens, adaptation to a specific host plant is a complex process possibly involving diverse molecular determinants and leading to host specificity [2, 3]. Understanding the molecular basis of host specificity can provide new insights into the evolution and ecology of specialist pathogens, and their potential to shift species and to infect new hosts. Bacteria from the genus Xanthomonas infect at least 392 plant species including important crops and ornamentals [4]. Yet, each individual strain is able to infect only one or few plant species. Strains able to cause the same symptoms on the same host range are grouped into pathovars [5]. Although our understanding of the molecular basis of host specificity is still limited, chemotactic sensors, adhesins and type III effectors emerge as key determinants for shaping host specificity in Xanthomonas [6–8]. Chemotactic sensors enable the bacteria to detect attractant or repellent molecules and trigger flagellar motility towards entry sites of the host plant, while adhesins allow the attachment on the host plant surface and biofilm formation, and type III effectors are delivered into the plant cells where they can have different functions including providing pathogen-associated molecular pattern triggered immunity (PTI).

The Xanthomonas axonopodis species complex sensu Vauterin [9, 10] groups more than 30 pathovars infecting a wide range of plants including economically important crops and ornamentals, such as Citrus, Anthurium and Dieffenbachia species, as well as pepper, cassava, cotton, mango, soybean, and common bean. Based on repetitive-sequence-based Polymerase Chain Reaction (rep-PCR) fingerprints, X. axonopodis has been subdivided into six subclusters named 9.1 to 9.6 [11]. More recently, this species complex has been split into the four species X. citri, X. euvesicatoria, X. phaseoli and X. axonopodis [12, 13]. Common bacterial blight of bean (CBB) is the most devastating bacterial disease infecting common bean (Phaseolus vulgaris L.). CBB occurs everywhere where common bean is cultivated and may cause up to 75% yield loss in the most severe cases [14, 15]. Xanthomonas strains responsible for CBB are distributed across four different genetic lineages [16]. The fuscous lineage (fuscans) and the non-fuscous lineages 2 (NF2) and 3 (NF3) belong to X. citri pv. fuscans while the non-fuscous lineage 1 (NF1) belongs to X. phaseoli pv. phaseoli [9, 11, 12]. Pathological convergence between the NF1 and fuscans lineages is associated with horizontal gene transfers (HGT) involving dozens of genes [17]. Horizontal transfer of genes encoding Transcription Activator-Like (TAL) type III effectors was also observed between the four lineages of CBB agents [18]. In particular, all strains from the four genetic lineages display an allele of the tal23A gene, suggesting that this gene is important for Xanthomonas adaptation to common bean.

In order to search for Xanthomonas genes putatively involved in the adaptation leading to common bean specificity in Xanthomonas, we have generated the whole genome sequences of 17 X. citri pv. fuscans and X. phaseoli pv. phaseoli strains. A combination of approaches including a comparison between the phylogeny of genes and the phylogeny of organisms, a parsimony approach to infer gene gains and losses, and a genome-wide search for specific k-mers, was used to search for genes presenting common characteristics unique to strains belonging to the four bean-pathogenic lineages of X. citri pv. fuscans and X. phaseoli pv. phaseoli.

Genome sequencing and phylogeny

In order to obtain genomic data representative of the diversity of Xanthomonas strains responsible for CBB, we produced whole genome sequences for 17 strains from the four genetic lineages of X. citri pv. fuscans and X. phaseoli pv. phaseoli that affect beans. In addition, we sequenced two strains of X. citri pv. mangiferaeindicae, one strain of X. citri pv. anacardii, three strains of X. oryzae pv. oryzicola, and used 51 other publically available Xanthomonas genomes for a total of 72 whole genome sequences (Table 1). Stenotrophomonas maltophilia strain R551–3 and Xylella fastidiosa strains 9a5c and Temecula1 were used as outgroups for further analyses [19–21]. Annotation revealed from 3209 to 5405 coding sequences (CDS) per Xanthomonas genome (Additional file 1). Among Xanthomonas strains responsible for CBB, chromosome size ranged from 4,957,446 bp in strain CFBP1815 to 5,517,999 bp in strain CFBP6992, with an average GC content ranging from 64.3 to 64.9%. The phylogeny of strains was assessed based on the amino acid sequences of all annotated CDS using CVTree (Fig. 1). The overall topology of this tree was congruent with previous Xanthomonas phylogenies [11, 12, 22]. As described previously, strains responsible for CBB are distributed into four distinct genetic lineages belonging to two different species, X. citri and X. phaseoli [12, 13, 16, 23].

Phylogeny of Xanthomonas strains used in this study. Phylogeny of Xanthomonas strains used in this study with indication of gene gains and losses. The phylogenetic tree is based on whole genome analysis using CVTree [66] with default parameters. Strain aliases are described in Table 1. Stenotrophomonas and Xylella genomes have been used as outgroups. Xanthomonas main phylogenetic groups 1 and 2 [24] and the X. axonopodis species complex [9, 10] are indicated by arrows. Groups 9.2, 9.4, 9.5 and 9.6 [11] are indicated in brackets. Fuscans, NF2, NF3 and NF1 refer to the four genetic lineages of strains responsible for CBB. A parsimony approach was performed to infer gene gains (blue) and losses (red) at levels higher than the pathovar rank, and numbers are displayed at each branch. Red stars highlight cases where gene loss was greater that gene gain. Curved arrows represent horizontal gene transfers (HGT) retrieved by Ks analysis on alignments of 115 candidate genes for bean specificity, with HGT from X. citri pv. fuscans to X. phaseoli pv. phaseoli in green, HGT from X. phaseoli pv. phaseoli to X. citri pv. fuscans in purple, and HGT between X. citri pv. fuscans lineages in red. Numbers in circles correspond to the numbers of candidate genes involved for each HGT. Question marks indicate events for which the origin or end of the HGT was not precise enough to assign any particular lineage

Genome expansion occurred during the evolution of Xanthomonas

To identify the genes shared by different clades of Xanthomonas, we constructed an orthology matrix using OrthoMCL (Additional file 2). Based on this orthology matrix, we performed a parsimony approach to infer gene gains and losses at each branch of the phylogenetic tree (Fig. 1). We did not take into account events occurring on the most distal branches to reduce the bias due to the difference of quality between sequenced genomes. At every branch, one to several hundreds of genes were either gained or lost. A general observation was that gene gains were higher than gene losses, suggesting that genome expansion occurred during the evolution of Xanthomonas (Fig. 1). Only four cases of genome reduction were observed (i) at the origin of the Xylella genus, (ii and iii) along two consecutive branches before and after the split between the X. oryzae species and the X. vasculorum and X. musacearum species, and (iv) at the origin of X. citri pv. malvacearum. The largest gene gain (418) was observed at the origin of Xanthomonas phylogenetic group 2 as defined by Young et al. [24], while the largest gene loss (819) was observed at the origin of the Xylella genus. Few gene losses (9 to 32) were observed before the diversification of each of the four genetic lineages involved in CBB. Of those, the NF1 lineage was the one which gained the most genes (271) followed by the NF2 (225), NF3 (108) and fuscans (83) lineages, respectively.

The pan and core genomes of Xanthomonas reveal extensive horizontal gene transfers between strains pathogenic on common bean

Individuals that are closely related to each other typically share more orthologs than unrelated individuals. Therefore, groups of closely related individuals tend to have a smaller pan genome and a larger core genome than groups of more divergent individuals. As such, the pan and core genomes for the 72 Xanthomonas strains comprised 32,602 and 1144 CDS, respectively, while the pan and core genomes for the 75 strains including the outgroups comprised 34,723 and 816 CDS, respectively (Fig. 2). Similarly, each Rademaker group alone, i.e. 9.2, 9.4, 9.5 and 9.6, had a smaller pan genome (6578, 8222, 9387 and 9437 CDS, respectively) and a larger core genome (3493, 2949, 3056 and 3213 CDS, respectively) than the X. axonopodis species complex, which had pan and core genomes of 19,010 and 2297 CDS, respectively. Strikingly, when grouping strains responsible for CBB belonging to X. citri pv. fuscans and X. phaseoli pv. phaseoli, both the pan and core genomes (10,750 and 3222 CDS, respectively) were larger than the pan and core genomes from groups 9.4 (8222 and 2949 CDS, respectively) or 9.6 (9437 and 3213 CDS, respectively) (Fig. 2). Thus, strains responsible for CBB, although being phylogenetically diverse, had more genes in common than they had with other strains belonging to their respective clades, which was suggestive of extensive HGT among these strains. This result was reminiscent of previous comparative analyses showing that dozens of genes have been horizontally transferred between the fuscans and NF1 lineages [17].

To search for genes potentially involved in the convergence between X. phaseoli pv. phaseoli (i.e. the NF1 lineage) and X. citri pv. fuscans (i.e. the NF2, NF3 and fuscans lineages) to infect common bean, we performed a combination of four different analyses. First, within the OrthoMCL matrix, we searched for CDS specifically present in the genomes of CBB agents and absent from any other Xanthomonas genome, or present in the genomes from all Xanthomonas but not in the genomes of CBB agents. No CDS was retrieved by this analysis. We also searched for CDS specifically present or absent when grouping the NF1 lineage to each of the NF2, NF3, or fuscans lineages. Only one CDS was specifically retrieved in the NF1 and fuscans lineages.

Second, we used the results from the CDS gains and losses approach described above to search for genes shared by all strains from X. phaseoli pv. phaseoli and X. citri pv. fuscans, and gained in the ancestor of one pathovar or the other. This approach unveiled nine CDS shared by all strains responsible for CBB, and gained in either X. phaseoli pv. phaseoli or X. citri pv. fuscans (Table 2). We also searched for CDS shared by the NF1 and each of the NF2, NF3, or fuscans lineage and gained in at least one of these lineages. Four CDS were shared by the NF1 and NF2 lineages, or the NF1 and NF3 lineages, while three CDS were shared by the NF1 and fuscans lineages.

Table 2

Numbers of CDS presenting similarities among the lineages of CBB agents

Lineages studied

Presence/absencea

Gainedb

Monophyleticc

24-mersd

NF1/NF2/NF3/fuscans

0

9

28

108

NF1/NF2

0

4

9

33

NF1/NF3

0

4

5

28

NF1/fuscans

1

3

105

231

Total

1

20

147

400

aCDS specifically present or absent in all the lineages studied compared to all other X. axonopodis strains

bCDS present in all the lineages studied and gained in a least one of these lineages

cCDS monophyletic for the lineages studied

dCDS containing 24-mers specifically present or absent in all the lineages studied compared to all other X.axonopodis strains

Third, we used a phylogenetic approach to search for genes for which strains from X. citri pv. fuscans and X. phaseoli pv. phaseoli formed a monophyletic group. For this, we constructed phylogenetic trees on 3202 CDS present in every X. citri pv. fuscans and X. phaseoli pv. phaseoli strain and in at least one additional strain from group 9.2 or 9.4 and one other strain from group 9.5 or 9.6. The additional strains from groups 9.2, 9.4, 9.5 and 9.6 were located inbetween X. citri pv. fuscans and X. phaseoli pv. phaseoli (Fig. 1). Thus, CDS found as monophyletic for CBB strains could be potential traces of HGT between both pathovars. This approach unveiled 28 CDS for which the four genetic lineages formed a monophyletic group, suggesting that they were horizontally transferred among these lineages (Table 2). Nine CDS were specifically monophyletic for the NF1 and NF2 lineages, five for the NF1 and NF3 lineages, and 105 for the NF1 and fuscans lineages, suggesting that most horizontal transfers occurred among the NF1 and fuscans lineages.

Finally, we used the SkIf tool [25] on the 72 Xanthomonas genomes to search for genes containing short 24-bp sequences (24-mers) specific to strains responsible for CBB, or alternatively genes from strains responsible for CBB lacking 24-mers present in all other strains from the X. axonopodis species complex. In all, we identified 108 CDS containing 24-mers either specifically present or absent from the four lineages (Table 2). Moreover, 33 CDS contained 24-mers specific for the NF1 and NF2 lineages, 28 for the NF1 and NF3 lineages and 231 for the NF1 and fuscans lineages. Similarly to the analysis based on phylogeny, this analysis based on k-mers pointed an overrepresentation of CDS with specific 24-mers shared by the NF1 and fuscans lineages compared to NF2 and NF3 lineages.

Together, these four analyses unveiled respectively 0, 9, 28, or 109 CDS presenting features unique to CBB agents. The analysis based on presence/absence seemed to be too stringent for unveiling any CDS, while the analysis based on k-mers was the most sensitive, suggesting that SkIf was an appropriate tool for finding common traits shared by phylogenetically distant strains. Most of these CDS were found redundantly by two or more analyses, for a total of 115 non-redundant CDS (Table 3). The most represented functions encoded by these 115 predicted CDS were hypothetical proteins (26 CDS), followed by membrane-related proteins (10 CDS), two-component system proteins (six CDS), putative secreted proteins (five CDS), reductases (five CDS), RNA-related proteins (five CDS), Type III secretion system-related proteins (five CDS), TonB-dependent proteins (four CDS), Type IV secretion system-related proteins (three CDS), Type VI secretion system-related proteins (three CDS), DNA-related proteins (three CDS) and transcription regulators (three CDS) (Table 3).

cgenes with 100% identity over 95% of their length according to Aritua et al. [17]

We hypothesized that the CDS potentially involved in the specific adaptation to common bean should bear nonsynonymous polymorphisms specific to Xanthomonas strains pathogenic on common bean. Analysis of the alignments for the 115 candidate CDS highlighted 44 CDS with nonsynonymous sites retrieved exclusively in X. citri pv. fuscans and X. phaseoli pv. phaseoli (Table 3). More than one third of these CDS (16/44) encoded hypothetical proteins. Among the other CDS, three encoded type IV secretion system proteins TrbI TraG and TraF, two encoded putative secreted proteins, two encoded type III secretion system proteins XopA and XopAD, two encoded DNA topoisomerases, and others encoded proteins of various functions (Table 3).

Specificity to common bean is associated with successive waves of horizontal gene transfers

Strain CFBP6546 from the NF1 lineage was used as reference for further analyses. Its genome contained one chromosome and three extrachromosomal plasmids formerly described as plasmid a, plasmid b and plasmid c in strain 4834-R [26]. Most candidate genes (95/115) were located on the chromosome, while 20 were located on the plasmid a (Fig. 3). This corresponds to a density of one candidate gene per 50.9 kbp in the chromosome, and one per 3.5 kbp for the plasmid a, while none were retrieved in plasmids b or c. Interestingly, all the CDS found in plasmid a contained specific nonsynonymous sites (Table 3). Thus, plasmid a appeared as an important vector of genes involved in the adaptation to common bean. Another observation was that the chromosome contained regions with various 24-mers shared by the NF1 lineage and any of the fuscans, NF2 or NF3 lineages (in green, blue or black in Fig. 3, respectively). This suggests that the regions shared by the NF1 and the other lineages diverged since the split between the NF2, NF3 and fuscans lineages. By contrast to what was observed for the chromosome, all specific 24-mers found in plasmid a were simultaneously shared by the four genetic lineages of strains responsible for CBB (in purple in Fig. 3), indicating that these regions have been shared between the four lineages recently enough to still have 100% identity between each other. Together, these results suggest that 24-mers retrieved in the chromosome correspond to more ancient HGT events than those retrieved in plasmid a.

Fig. 3

Mapping of the 24-mers specific for strains pathogenic on common bean. The innermost rings represent the reference chromosome or plasmids with associated coordinates. Colored lines represent 24-mers specifically retrieved in X. citri pv. fuscans and X. phaseoli pv. phaseoli strains (purple), or in the NF1 plus fuscans lineages (green), or in the NF1 plus NF2 lineages (blue), or in the NF1 plus NF3 lineages (black). Red numbers correspond to the identifiers of the 115 genes listed in Table 3

Nucleotide synonymous substitution rates at silent sites (Ks) is an estimation of neutral evolution because it does not take into account the nonsynonymous sites that can be under selection pressure. Therefore, Ks can be used as an approximation of the time of divergence between genes or taxa, with higher Ks value meaning longer time of divergence between two sequences [27, 28]. For each of the 115 candidate genes found in CBB agents, we performed multiple alignments. We could not perform Ks analysis on 15 genes that were lacking outgroups (Table 3), therefore we tested only 100 out of 115 genes. Among these 100 genes, 18 were recombinants according to RDP software analysis (Table 3). For these 18 recombinants, Ks values were independently calculated on both sides of the breakpoints. We calculated pairwise Ks values for different combinations of strains including X. citri pv. fuscans, X. phaseoli pv. phaseoli and closely related strains including X. citri pv. anacardii, X. citri pv. aurantifolii, X. phaseoli pv. syngonii, X. phaseoli pv. dieffienbachiae and X. phaseoli pv. manihotis (Fig. 1, Additional file 3). We then used these Ks as relative time divergence estimations to infer if a HGT occurred between NF1, NF2, NF3 and/or fuscans lineages, as well as the direction of this HGT. For example, for gene m00100580b the mean Ks value between strains from the NF1 and fuscans lineages was 7.40e-03 +/− 2.85e-10. This value was lower than the Ks values when comparing the fuscans lineage to its closest relatives from the NF2 or NF3 lineages (Ks = 4.58e-02 +/− 1.32e-09 or 4.08e-02 +/− 1.54e-09, respectively), or when comparing the NF1 lineage to its closest relative X. phaseoli pv. manihotis (Ks = 5.33e-01 +/− 0.00). These results indicate that m00100580b was more similar between the NF1 and fuscans lineages than between these lineages and their closest relatives, meaning that m00100580b was horizontally transferred between the ancestors of the NF1 and fuscans lineages. Moreover, the Ks value between the NF1 lineage and the NF2 or NF3 lineages (Ks = 4.36e-02 +/− 1.04e-09 or 4.01e-02 +/− 0.00, respectively) was lower than between the NF1 lineage and X. phaseoli pv. manihotis (Ks = 5.33e-01 +/− 0.00). Therefore, m00100580b was closer between NF1 strains and other strains from group 9.6 than from it’s closest relatives, meaning that the horizontal transfer was directed from the fuscans lineage to the NF1 lineage. This analysis confirmed HGT for 88 out of 100 genes tested (Fig. 1, Table 3, Additional file 3). The vast majority of HGT was directed from X. citri pv. fuscans to X. phaseoli pv. phaseoli, while only four HGT occurred from X. phaseoli pv. phaseoli to X. citri pv. fuscans. In particular, 55 HGT events were detected from the fuscans lineage to the NF1 lineage. In addition to having been transferred between distant lineages, 16 and 1 genes were also transferred between the fuscans and the NF2 lineages, or the fuscans and the NF3 lineages, respectively (Fig. 1). Moreover, eight genes had Ks = 0.00 +/− 0.00 between the NF2 and NF3 lineages, and nine between the NF2, NF3 and fuscans lineages, suggesting that HGT events also occurred between these lineages (Fig. 1, Additional file 3). Together, our results show that several more or less important waves of HGT occurred between the ancestors of phylogenetically distant strains responsible for CBB.

Finally, GC content is often used as a mean to detect HGT from foreign origin [29]. Out of the 115 candidate CDS, only two CDS (m00105390 and m00200160) presented an atypical GC content (α < 0.05) within the genome of strain CFBP6546 (Table 3). This result was not unexpected, as all strains from the NF1, NF2, NF3 and fuscans lineages have a similar GC content around 64% (Table 1), therefore HGT between these strains was not expected to result in a shift of GC content.

We performed a comparative genomics analysis to detect genes putatively involved in Xanthomonas specificity to common bean. For this, we generated the whole genome sequence from 17 strains representing the diversity of the four genetic lineages belonging to X. citri pv. fuscans and X. phaseoli pv. phaseoli. We used a combination of comparative genomics approaches that led to the discovery of 115 genes bearing features unique to CBB agents. Out of these 115 genes, 108 were retrieved using the SkIf tool based on specific 24-mer search [25]. Previous analyses based on identity percentage unveiled 63 genes sharing 100% identity over at least 95% of their length among strains from the NF1 and fuscans linages [17]. Only nine of these genes were retrieved within our list of 115 genes (Table 3). This difference can be explained by the fact that we discarded most of the genes shared only by the NF1 and fuscans lineages and retained genes similar in all four genetic lineages of CBB agents (Table 2). On the other hand, we unveiled numerous genes that did not share 100% identity among the NF1 and fuscans lineages for their whole length, but instead shared small specific sequences of 24 nucleotides or more. Whether these similarities lie within functionally important domains of the encoded proteins remains to be studied.

Ks comparisons, showed that a majority of these genes were involved in HGT between X. citri pv. fuscans and X. phaseoli pv. phaseoli. Therefore, HGT was the predominant force leading to similarities between the genomes of X. citri pv. fuscans and X. phaseoli pv. phaseoli. Finding HGT events within these genes validated our approach. In particular, SkIf was an interesting tool because in addition to being more sensitive than gene gain and loss or phylogenetic approaches, it was not based on gene alignments, thus less sensitive to annotation and/or sequencing biases. HGT events occurred at different moments of the evolution of Xanthomonas strains having common bean as host. The vast majority of these HGT were directed from X. citri pv. fuscans to X. phaseoli pv. phaseoli. This strongly suggest that X. citri pv. fuscans was originally pathogenic on common bean, and that X. phaseoli pv. phaseoli subsequently acquired the ability to cause CBB on bean due to successive acquisitions of novel genes and/or novel alleles coming from the three X. citri pv. fuscans lineages. This can be compared to our knowledge on the origin of the lineages and their genetic diversity. The causal agent of CBB was first isolated and identified by Smith in 1897 [30] as a yellow pigmented strain, later shown as belonging to the NF1 lineage. Burkholder later isolated the first fuscous strains from beans grown in Switzerland in 1924 [31]. However there are no data to document a putative pre-existence of one, the other, or both types of strains prior to their first identifications. The genetic diversity of the yellow and fuscous strains was revealed by various methods. Amplified or restriction fragment length polymorphism analyses [16, 32–34], amplified polymorphic DNA fragments [32, 35], pulse field gel electrophoresis [33], and multilocus sequence analysis [23] all revealed that both types of strains are more or less equivalent in terms of genetic diversity. This suggests that diversification of both lineages occurred around the same time, and thus that the ancestors of these two lineages may have coexisted. As a consequence, X. citri pv. fuscans may be the descendant of the original CBB agent that had transferred determinants useful for adaptation on common bean to the ancestor of X. phaseoli pv. phaseoli. Therefore, the ancestor of X. phaseoli pv. phaseoli appears as a recombinant that emerged as a new common bean pathogen through the acquisition of novel genes and alleles. In quarantine areas such as Europe, Turkey, Barhain, Azerbaijan and Israel, seed lots are routinely tested using a method from the International Seed Testing Association involving isolation of bacterial strains, pathogenicity tests, and specific PCR assays [36, 37]. Our results could serve for improving PCR-based monitoring of CBB agents by designing PCR primers on genes presenting sequences unique to CBB agents and potentially important for common bean specificity. Such primers could potentially detect novel HGT of these genes in strains unrelated to X. phaseoli pv. phaseoli and X. citri pv. fuscans, thus allowing to forecast new threads potentially dangerous for common bean production.

Very diverse functions were retrieved among the proteins encoded by the 115 candidate genes. Interestingly, 44 genes contained nonsynonymous polymorphisms specific for strains responsible for CBB, suggesting that they may play an important role in common bean specificity. Although 17 of these 44 genes encoded hypothetical proteins, it appears that most other genes encoded proteins involved in pathogenicity or in the interaction with the plant environment. The type IV secretion system was particularly represented with genes encoding TraF, TraG and TraI, and is involved in the translocation of macromolecules such as proteins important for pathogenicity, or DNA for mediating HGT [38, 39]. Thus, sharing similar proteins of the type IV secretion system could favour HGT among strains responsible for CBB. Indeed, strains from the NF1, NF2, NF3 and fuscans lineages have been found in La Réunion Island in 2000 (Table 1), indicating that sympatry exists among all these lineages, rendering further HGT events possible [18]. Three genes encoding proteins related to the type III secretion system were retrieved. XopA [40] and HpaH [41] are two proteins that may be involved in the structure of the type III secretion system, while XopAD is a type III effector of unknown function consisting of multiple semi-conserved 42 amino acids SKW repeats [42, 43]. The type III secretion system is pivotal for the virulence of most Gram negative plant pathogenic bacteria, and repertoires of type III effectors have been described as potentially important factors for host specificity and host adaptation in Xanthomonas [7, 44, 45] and other genera such as Pseudomonas [46]. Moreover, our analysis pointed out one diguanilate cyclase and one cGMP specific phosphodiesterase, two proteins involved in the metabolism of cyclic di-GMP that may play a role in biofilm formation [47] and pathogenicity [48]. One TonB-dependent transporter was also retrieved. TonB-dependent transporter are outer membrane receptors involved in molecule uptake such as iron siderophore complexes or nutrients and may play a role in host specificity [3, 49]. Other proteins putatively involved in pathogenicity were retrieved, such as ThrC, a threonine synthase involved in the virulence of X. oryzae pv. oryzicola in rice [50], XpsD, an outer membrane protein from the type II secretion system that is putatively involved in the secretion of cell wall degradative enzymes during infection [51], or IcmF, a protein of the type VI secretion system, which is involved in the interaction with other bacteria and may participate in pathogenicity [52, 53]. One flavine reductase and one xanthine dehydrogenase, two proteins putatively involved in oxydoreduction pathways were retrieved, and may be involved in the response to stress during the interaction with common bean. In addition to genes putatively involved in pathogenicity or virulence, our analysis unveiled genes involved in more general metabolism pathways, such as PhnB involved in the biosynthesis of tryptophan [54], or two DNA topoisomerases involved in the relaxation of the supercoiled DNA molecule during transcription, replication or recombination [55]. On one hand, our analysis unveiled only one or few genes within a given function, while the functions retrieved correspond to pathways often involving dozens of genes. This suggested that slight modifications within a given pathway would be sufficient to impact host specificity. On the other hand, the genes retrieved here encompass almost all the stages of host plant colonization by the bacteria, from the ability to mobilize trophic resources for multiplication to the interaction with other microorganisms, biofilm formation, response to oxidative stress, and inhibition of plant defences. Therefore, the ability to infect a particular plant seems to require not just one or a few adaptative determinants but an arsenal of factors allowing a global adaptation to a specific niche including the plant and, as a consequence a fine tuning and coordination of the activity of these determinants.

Interestingly, 19 of these 44 candidate genes were retrieved on plasmid a, suggesting that this plasmid played a major role for pathological convergence of CBB agents. Plasmid a carries an additional type III effector gene encoding an effector from the Transcription Activator-Like (TAL) family that was horizontally transferred between the NF1, NF2, NF3 and fuscans lineages [18]. Plasmids are genetic elements that favour HGT, but transfers of whole plasmids often induce a fitness cost for the bacteria [56]. More generally, horizontally transferred genes tend to be lost if not providing selective advantages for recipient strains [57]. Interestingly, the nine candidate genes retrieved by the gain and loss approach were all located on plasmid a (Table 3). The maintenance of these novel genes in the four genetic lineages of CBB agents is a testament to the importance of these genes for the bacteria. Except for two genes encoding proteins involved in the type IV secretion system, the other seven genes encoded proteins of unknown function. It would be interesting to perform functional characterization of these genes, and further analyse their implication in common bean specificity. Analysing the expression patterns during infection would be a natural extension of this study, and a first step towards functional validation of these genes.

Together, our results indicate that consecutive waves of HGT occurred between phylogenetically distant Xanthomonas strains able to cause the same symptoms on the same host plant: common bean. These HGT led to specific combinations of genes only retrieved in strains responsible for CBB, which provided new insights into the evolution of these bacteria towards infecting common bean. Mining for candidate genes for host specificity could be generalized to other polyphyletic pathovars such as pathovars euvesicatoria, vesicatoria, perforans, and gardneri forming a group of strains pathogenic on pepper and tomato [8]. Such analyses could both give new information on the molecular bases of host specificity, and provide new tools for enhancing epidemiological surveillance of strains pathogenic on a given host, or detecting recombinant strains presenting a high potential of emergence through the acquisition of novel genes.

Bacterial genomes and strains

All strains used in this study are listed in Table 1. The strains used for genome sequencing were provided by the CIRM-CFBP (International Center for Microbial Ressources - French Collection for Plant-associated Bacteria, https://www6.inra.fr/cirm_eng/). Genome sequencing was performed using the following procedure. Genomic DNA was prepared from overnight liquid cultures of bacteria previously grown on 10% TSA medium (tryptone at 1.7 g/L, soybean peptone at 0.3 g/L, glucose at 0.25 g/L, NaCl at 0.5 g/L, K2HPO4 at 0.5 g/L, agar at 15 g/L, pH 7.2) for 2 days at 28 °C. DNA was extracted and purified by the method of Klotz and Zimm [58]. Illumina sequencing was performed by Genoscope (20 strains, paired end reads of 300/500 bp) or GATC Biotech (three strains, with combined paired-end reads of ca. 250 bp and 3 kb mate-pair reads). Genome assembly was performed using a combination of SOAPdenovo (version 2.04) [59], SOAPGapCloser (version 1.12) [60] and Velvet (version 1.2.02) [61] assemblers. Sequenced genomes were estimated to be > 93% complete and < 3% contaminated (Additional file 4) using CheckM (version 1.0) [62]. The pathogenicity of all CBB strains was confirmed on common bean plants from cultivar Flavert as described in Ruh et al. [18]. The seeds from cultivar Flavert were kindly provided by Vilmorin (La Ménitré, France) and are available at the bean collection of the CIAT (Center for Tropical Agriculture, Colombia, http://genebank.ciat.cgiar.org/genebank/main.do).

Annotation and phylogenomic analyses

Structural and functional annotation of whole genome assemblies was performed using the automated pipeline Eugene-PP (version 1.2) [63], using SWISS-PROT as protein database and training protein database (http://www.uniprot.org/). Additional functional annotation of all predicted CDS was performed with InterProScan (version 4) [64]. A presence/absence matrix of ortholog groups was constructed using OrthoMCL (version 2.0) on amino acid sequences from all predicted CDS at an inflation index of 1.5 [65]. This matrix was then used for defining core and pan genomes. Phylogenetic trees were constructed using CVTree (version 4.2) [66] using the aminoacid sequences of all predicted CDS from the 75 genomes used in this study. CDS gains and losses were analysed using the Most Parcimonious Reconstruction function from the APE package (version 3.2) [67] to search for the most parsimonious succession of events explaining the repertoire of ortholog groups at each node of the phylogenetic tree.

A phylogenetic approach was used to search for genes for which strains from X. citri pv. fuscans and X. phaseoli pv. phaseoli form a monophyletic group. For this we selected 3202 CDS using an R script to search the orthology matrix for genes that were present in all X. citri pv. fuscans and X. phaseoli pv. phaseoli strains and in at least one strain from Rep-PCR group 9.2 or 9.4 [11] plus one another strain from Rep-PCR group 9.5 or 9.6, in order to avoid getting trees were X.citri pv. fuscans and X. phaseoli pv. phaseoli appear monophyletic due to a lack of correspondig genes in the strains inbetween. CDS were aligned using MAFFT (version 7) with L-INS-I strategy [68]. Neighbour-joining trees were constructed using APE (version 3.2) under the Kimura 80 model [67]. CDS monophyletic for all X. citri pv. fuscans and X. phaseoli pv. phaseoli strains, or alternatively for the NF1 and another lineage (i.e. NF2, NF3, or fuscans), were retrieved using the APE package (version 3.2) [67].

A k-mer-based approach was used to search for genes containing short specific sequences present in all strains from X. citri pv. fuscans and X. phaseoli pv. phaseoli but absent in other strains. For this, we used SkIf (version 1.0) [25] with a k-mer size of 24 (or 24-mer), and using X. citri pv. fuscans strain 4834-R genome as reference. The same approach was used to search for genes containing 24-mers absent in strains belonging to X. citri pv. fuscans and X. phaseoli pv. phaseoli but conserved in all other strains from the X. axonopodis species complex, using X. citri pv. anacardii strain CFBP2913 genome as reference.

Recombination and HGT analyses

The 115 genes presenting specific traits of adaptation to common bean were aligned using MAFFT (version 7) with L-INS-I strategy [68]. Intragenic recombination events were then searched using a suite of programs implemented in RDP (version 4.16) [69], RDP [70], Geneconv [71], MaxChi [72], Chimaera [73], Bootscan [74] and 3seq [75]. Default parameters were used for each method except for Bootscan (window = 150, step = 20, neighbor joining trees, 200 replicates, 95% cut-off, J&N model with Ti:Tv = 2, coefficient of variation = 2). Ks was calculated using DNAsp (version 5) [76]. For each gene, the occurrence and dating of HGT events were estimated by comparing Ks values from 28 different pairwise combinations listed in Additional file 3. For example, NF1 and fuscans strains belong to phylogenetically distant strains, thus if the Ks between strains from genetic lineages NF1 and fuscans was lower than the mean Ks between other lineages, it was indicative of recent HGT between the ancestors of NF1 and fuscans. Direction of events were assessed by comparing the Ks values for outgroups belonging to Rep-PCR groups 9.4 and 9.6 (Fig. 1). For recombinants, separate analyses were performed for each region on both sides of the recombination point.

Acknowledgements

The authors thank Sébastien Carrère for his help on genome annotations and Laurent D. Noël for sharing the genome sequence of strain CFBP6865. The authors also thank the French Network on Xanthomonads (FNX) (https://www.reseau-xantho.org/) for recurrent scientific exchanges on Xanthomonas. Authors benefited from interactions promoted by COST Action FA 1208 (https://www.cost-sustain.org). We thank the CIRM-CFBP (Beaucouzé, France) for strain preservation and supply.

Funding

The research leading to these results has received grants from the Genoscope (3X 154/AP2006-2007, XANTHOMICS 18/AP2009-2010) and the French National Research Agency (XANTHOMIX ANR-2010-GENM-013-02). LSG was funded by a postdoctoral grant from the XANTHOMIX ANR project, and another postdoctoral grant from Angers-Loire Metropole, France. MR was funded by a PhD grant from Angers-Loire Metropole, France. This work was supported by France Génomique National infrastructure, funded as part of “Investissement d’avenir” program managed by Agence Nationale pour la Recherche (contrat ANR-10-INBS-0009).

Availability of data and materials

The datasets (genome sequences) generated through various projects and used in this study have been deposited in GenBank under accession numbers listed in Table 1.

Authors’ contributions

MAJ, LSG and NWGC designed the study. This study was partly lead by LSG during her post-doc and is part of the doctoral research of MR supervised by MAJ and NWGC. MAJ, NWGC, LG and RK obtained grants to fund the study. VB managed the genome sequencing at Genoscope. NWGC, LSG, MAJ, MB, MR, SB, AD, and VB contributed to data analyses and interpretation. NWGC, MR and MAJ wrote the manuscript with inputs from all co-authors. All authors read and approved the final version of the manuscript.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.