Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

Identification of genetic relationships and subspecies signatures in Xylella fastidiosa

BMC Genomics volume 20, Article number: 239 (2019) Cite this article

Abstract

Background

The phytopathogenic bacterium Xylella fastidiosa was thought to be restricted to the Americas where it infects and kills numerous hosts. Its detection worldwide has been blooming since 2013 in Europe and Asia. Genetically diverse, this species is divided into six subspecies but genetic traits governing this classification are poorly understood.

Results

SkIf (Specific k-mers Identification) was designed and exploited for comparative genomics on a dataset of 46 X. fastidiosa genomes, including seven newly sequenced individuals. It was helpful to quickly check the synonymy between strains from different collections. SkIf identified specific SNPs within 16S rRNA sequences that can be employed for predicting the distribution of Xylella through data mining. Applied to inter- and intra-subspecies analyses, it identified specific k-mers in genes affiliated to differential gene ontologies. Chemotaxis-related genes more prevalently possess specific k-mers in genomes from subspecies fastidiosa, morus and sandyi taken as a whole group. In the subspecies pauca increased abundance of specific k-mers was found in genes associated with the bacterial cell wall/envelope/plasma membrane. Most often, the k-mer specificity occurred in core genes with non-synonymous SNPs in their sequences in genomes of the other subspecies, suggesting putative impact in the protein functions. The presence of two integrative and conjugative elements (ICEs) was identified, one chromosomic and an entire plasmid in a single strain of X. fastidiosa subsp. pauca. Finally, a revised taxonomy of X. fastidiosa into three major clades defined by the subspecies pauca (clade I), multiplex (clade II) and the combination of fastidiosa, morus and sandyi (clade III) was strongly supported by k-mers specifically associated with these subspecies.

Conclusions

SkIf is a robust and rapid software, freely available, that can be dedicated to the comparison of sequence datasets and is applicable to any field of research. Applied to X. fastidiosa, an emerging pathogen in Europe, it provided an important resource to mine for identifying genetic markers of subspecies to optimize the strategies attempted to limit the pathogen dissemination in novel areas.

Background

Xylella fastidiosa is a species of plant pathogenic bacteria endemic in the Americas, but listed as quarantine pests elsewhere (https://gd.eppo.int/taxon/XYLEFA/categorization). However, since 2013, various cases of emergences have been reported in Europe (Italy, France, Germany and Spain) on large ranges of host plants including olive trees, grapevine, and ornamentals [1,2,3,4,5]. In Italy, assuming that X. fastidiosa started spreading in 2010, a recent model approach suggested that it will progress through olive orchards to infect the northernmost recorded orchards within 43.5 years [6]. In France, bacterial introduction was estimated between 1985 and 2001, depending on the modeled scenarios [7, 8]. Several records of X. fastidiosa in imported materials (i.e. mostly coffee plants) were also reported over the same period in Europe [9,10,11,12].

X. fastidiosa is a genetically diverse species that is currently divided into six subspecies (subsp. fastidiosa, pauca, multiplex, sandyi, morus, tashke), the four-first being the most damaging and now being found in numerous countries worldwide. But the genetic diversity of the genus Xylella is undoubtedly underestimated. Yet another species, X. taiwanensis, was recently proposed for the strains causing leaf scorch on nashi pear tree, a disease that was reported more than 25 years ago in Taiwan and initially thought to be caused by a X. fastidiosa strain [13]. Recombination is known to drive X. fastidiosa evolution and adaptation to novel hosts [11, 14,15,16]. For example, the subspecies morus has been proposed for grouping strains issued from large events of intersubspecific recombination that were associated with a host shift [16]. The recent outbreaks and interception of imported, contaminated materials in Europe as well as investigations in South America also revealed the existence of previously unknown Sequence Types of several X. fastidiosa subspecies [2, 9, 11, 17, 18].

Because management and regulations of X. fastidiosa outbreaks in France depends on the subspecies of X. fastidiosa, it is of major importance to precisely define these subspecies, understand the robustness of these groupings and their meaning in terms of specific or shared genetic material. One way to resolve such a series of interrogations is the achievement of comparative genomics to identify similarities and specificities between groups of individuals. Yet, exploring big datasets is not trivial and requires dedicated bioinformatic tools to be cost- and time-effective. Various applications use k-mers mostly to analyze sequence reads to improve the quality of genome, transcriptome and metagenome assembly [19,20,21,22,23,24,25,26]. K-mer are all the possible substrings of length k that are contained in a nucleotide character string. K-mer-based methods can also be employed on whole genome sequences to taxonomically assign organisms [27, 28]. Moreover, several tools were developed to calculate pairwise relationships, like the average nucleotide identities using blast (ANIb) or MUMmer (ANIm) algorithms, and the tetranucleotide frequency correlation coefficients (TETRA), which can be accessed online through JSpecies [29] or with workstation installation of python3 pyani module [30].

Here, we developed SkIf (Specific k-mers Identification) and applied it to gain a better understanding of X. fastidiosa clustering in subspecies through the detection of genomic regions specifically associated with X. fastidiosa subspecies. We also used this tool to identify specific k-mers within 16S rRNA gene and assess the occurrences of X. fastidiosa in the SILVA database as a first attempt to mine large databases to evaluate the worldwide dispersion of X. fastidiosa subspecies.

Results

The genome sequence dataset

The dataset used in this study gathered 47 Xylella genomes sequences, including 46 X. fastidiosa and one X. taiwanensis specimen (Table 1). The X. fastidiosa subspecies tashke could not be included as no strain or genome sequence are available. In some analyses, the three strains belonging to X. fastidiosa subsp. sandyi were separated into two groups, containing either the original strain Ann-1 (sandyi) or the more recently discovered relatives CO33 and CFBP 8356 (sandyi-like) both belonging to the unusual sandyi ST72 [17]. The CFBP 8073 strain, described as an atypical X. fastidiosa subsp. fastidiosa strain [11] was either included or not for analyses of this subspecies. The X. fastidiosa genome sequences involve 39 publicly available ones and seven newly individuals. The strains sequenced in this work were selected based on their country of isolation, genetic diversity and host range, inferred from their belonging to the subspecies fastidiosa (CFBP 7969, CFBP 7970, CFBP 8071, CFBP 8082, and CFBP 8351), sandyi (CFBP 8356) and multiplex (CFBP 8078). Genome sequence characteristics are described in Table 2.

Table 1 List of the 47 Xylella genome sequences used in this study
Full size table
Table 2 List of Xylella fastidiosa strains sequenced for this study and genome properties
Full size table

Evaluation of strain synonymy using k-mers

CFBP 7970, the X. fastidiosa and X. fastidiosa subsp. fastidiosa type strain [31], has various synonymous names (ATCC 35879, DSM 10026, LMG17159, http://www.straininfo.net/strains/901514 or http://www.bacterio.net/xylella.html) in other collections, as a result of strain exchanges between the American (ATCC), German (DSMZ), Belgian (LMG), and French (CIRM-CFBP) collections [32]. As no genome sequences were available at the beginning of this study for any of these strains, CFBP 7970 was included in our dataset and its genome was sequenced. Later on, the genome sequences of the strains ATCC 35879 and DSM 10026 were released. Although the genome sequences of the three strains were very similar (99.83–99.97% ANIb), they were not strictly identical (Additional file 1). The use of SkIf identified DNA fragments that were specific to two genome sequences, but absent in the third one, yielding in 95, 192, and 594 k-mers specifically present in the pairs ATCC 35879/DSM 10026, CFBP 7970/DSM 10026, ATCC 35879/CFBP 7970, respectively (Additional file 2). The absence of some mers in genome sequence could be due to sequencing artefacts (e.g. sequencing technology employed, average coverage and assembly methods) that resulted in specific SNPs (Table 3). However, 16 mers (ranging from 29 to 9,178 nt and totalizing 36,845 bp in size) were detected in a single contig (41,458 bp) in CFBP 7970 and into several DSM 10026 contigs but were absent from ATCC 35879 genome sequence (Additional file 2). These sequences shared high identity levels with a plasmid found in multiple Xylella subspecies, known as pXF-De Donno (subsp. pauca De Donno strain), pXF-RIV5 (subsp. multiplex RIV5 strain), pXF-FAS01 (subsp. fastidiosa M23 strain) or present but unnamed elsewhere (like in subsp. pauca CoDiRO strain and subsp. multiplex Dixon strain) [2, 33,32,35]. But the blast analysis suggested its possible presence, as partial matches (3.8 -6 kb in total) with high but not perfect identity levels (95–99%) were found when searched against the genome sequence of ATCC 35879 (Additional file 2; Additional file 3).

Table 3 Properties of genome sequences of strains CFBP 7970, DSM 10026 and ATCC 35879
Full size table

Data mining of the 16S rRNA SILVA database to assign occurrences of Xylella sp.

The availability of the 47 Xylella genome sequences renders possible the analysis of the allelic diversity of the 16S rRNA marker gene. This housekeeping marker is widely used in bacterial phylogeny and taxonomy studies for various reasons including its vertical inheritance and ubiquity in prokaryotes. It is also commonly used to survey microbial communities and as such is a marker to survey largely the environment. A total of 74 16S rRNA gene sequences was retrieved from our dataset, these either being present in one (n = 20) or two copies (n = 27) in X. fastidiosa and X. taiwanensis (Table 4). We detected 19 SNPs (over 1547 nucleotides, 1.22%) specific to X. taiwanensis PLS 229 16S rRNA (Additional file 4).

Table 4 Repertoire of 16S rRNA gene sequences in 47 genomes of Xylella sp
Full size table

Specific mers were searched for within the Xylella genus (i.e., the in-group included the 74 Xylella 16S rRNA copies; the out-group included all the SSU sequences from the Silva database other than Xylella-tagged) and the X. fastidiosa species (i.e., the genome sequence of X. taiwanensis PLS 229 strain was included in the out-group). Five long-mers (referred to as LongXyl#1 to #5) specific to Xylella genus and four long-mers (referred to as LongXylefa#1 to #4) specific to the X. fastidiosa species were identified. LongXyl (23-43 nt) and LongXylefa (23-31 nt) mers located between positions 212 and 866, and positions 202 and 1013, respectively, which include the V3-V4 hypervariable regions widely used in community profiling approaches (Additional file 5). Specific signatures obtained from eight nucleotide positions in the 16S rRNA alignment discriminated alleles from X. fastidiosa subsp. fastidiosa, multiplex, morus, sandyi and pauca (Table 5).

Table 5 Specific signatures in 16S rRNA nucleotide sequences to discriminate X. fastidiosa subspecies
Full size table

The occurrence of these specific mers in 16S rRNA nucleotide sequences was investigated within the SILVA rRNA database (Additional file 4). A large proportion of the sequences retrieved from the SILVA database (n = 118/195) covered less than half of the total gene length, and only a minority (n = 70/195) covered at least three fourth of the length. After nucleotide alignment, 53 of these sequences including the LongXyl and LongXylefa specific signatures were retained (Additional file 4). Based on their genetic signatures, 51 sequences were assigned to subsp. multiplex (n = 32), fastidiosa (n = 11), pauca (n = 5), sandyi (n = 2) and morus (n = 1), while two sequences (EU560720.1 and EU560722.1) could not be assigned because they did not cover most of the eight discriminant nucleotide positions (Additional file 5). The validity of this presumptive taxonomic affiliation was consolidated by the description of the sample, which were in adequacy with the current host range of these subspecies. Indeed, samples from subsp. fastidiosa mainly come from alfalfa or grapevine, while those from subsp. multiplex were collected on various oak species, those from subsp. sandyi were isolated from oleander, those from subsp. pauca from periwinkle, coffee or olive trees, and the one from subsp. morus came from mulberry.

Identification of allelic variants specific to each X. fastidiosa subspecies

Beyond focusing on a single gene (16S rRNA), we applied SkIf to a whole genome-based analysis. Seven groups of strains were defined: fastidiosa (two groups), pauca, multiplex, morus and sandyi (two groups). The two fastidiosa groups differed by the presence/absence of the strain CFBP 8073, while one sandyi group included only the original member Ann-1, and the sandyi-like group included only CFBP 8356 and CO33 strains. Specific mers were searched for in each group against all the others.

Overall, long-mers were identified all along the genomes (Fig. 1), mainly matching coding sequences (71–80% depending on the subspecies), with one to several long-mers in the identified CDS (Table 6A; Additional file 6). Gene set enrichment analysis was performed by comparing the predicted functions associated with these specific CDSs to the overall predicted proteomes. Fischer’s exact test revealed multiple GO terms over- (mostly) or under- (rarely) represented for the seven groups. Overall, only one GO term (catalytic activity) was always found enriched, except for the subspecies morus, while 10 other GO terms were found enriched in all groups except in subspecies morus and sandyi (Table 7; Additional file 7; Fig. 2). Several GO terms were only identified in a single subspecies (Table 8), suggesting that associated mechanisms might be key markers of X. fastidiosa subspecies evolution. As for subspecies pauca it concerned 175 GO terms, including 20 terms associated with the bacterial cell wall/envelope/plasma membrane and 16 related to nucleotide metabolic/biosynthetic process, especially for purine, as well as 4 terms under-represented dealing with viral or symbiont processes. As for subspecies fastidiosa (without CFBP 8073) it concerned six GO terms related to DNA modification and vitamin process. The subspecies multiplex specific GO terms deal with metabolic process, catalytic activity and conformation of DNA and organelle organization. The subspecies morus had only one ontology enriched, associated with DNA replication.

Fig. 1
figure 1

Distribution of k-mers along the X. fastidiosa genome sequences. Frequency of core (black) and specific (colored) k-mers mapped onto the genome of reference (mentioned into brackets) of each subspecies. a subsp. fastidiosa with or without CFBP 8073 strain. b subsp. sandyi. c subsp. morus. d subsp. multiplex. e subsp. pauca

Full size image
Table 6 Main features related to the specific mers identified in X. fastidiosa subspecies
Full size table
Table 7 Main Gene Ontologies (GO) identified as enriched in almost all the subspecies for the CDS harboring specific mers
Full size table
Fig. 2
figure 2

Overlap between gene ontologies differentially represented (over or under) in genes harboring specific k-mers. a Relationships between six groups: FAS and FAS2, subsp. fastidiosa without or with CFBP 8073, respectively; MUL, subsp. multiplex; SAN and SAN2, subsp. sandyi and sandyi-like, respectively; PAU, subsp. pauca. Note: the subsp. morus is not indicated on the Venn diagram as it was identified only one GO term, specific to it. b Relationships between three groups (pauca, multiplex and the third one resulting from the grouping of subsp. fastidiosa, sandyi and morus)

Full size image
Table 8 Selected differentially represented Gene Ontologies of CDS with specific mers in X. fastidiosa subspecies or subclades
Full size table

Reconstruction of the parental origin of the subspecies morus

The subspecies morus was proposed to group strains pathogenic on Morus that derived from largescale intersubspecific homologous recombination events between ancestors from at least subspecies fastidiosa and multiplex. This assumption is based on the analysis of seven housekeeping genes [16]. We challenged it with whole-genome sequence datasets to further understand the contribution of the subspecies fastidiosa, multiplex and others in the parenthood of the subspecies morus. We used SkIf to identify mers specific of the morus group (i.e. two strains, Mul-MD and Mul0034) plus one of each of the other groups (Table 6B; Additional file 6). The highest level of specific mers was found for the combination morus x fastidiosa x sandyi: the highest mer cumulated size represented 5% of the Mul0034 genome in size (which does not mean that all the 5% are unique to Mul0034). The morus x multiplex (3.6%) and morus x pauca (< 0.1%) relationships were lower. To illustrate these findings, the specific mers were mapped onto the Mul0034 genome and were found to be distributed all along the sequence (Fig. 3). These results on closest relationships between genomes of subspecies fastidiosa, sandyi, and morus are coherent with the ANIb values (Additional file 1). Enrichment tests identified shared GO terms for the various combinations indicated (Table 6B). For the combination morus x multiplex, one GO term linked to the amino acid biosynthetic process was specifically recorded. At the level of fastidiosa x sandyi x sandyi-like x morus, 28 GO terms were specifically identified, associated with various processes like cellular component or protein complex disassembly, peptidyl-proline activity and chemotaxis.

Fig. 3
figure 3

Distribution of k-mers specific to the X. fastidiosa subspecies morus and others. a core k-merome of X. fastidiosa species. b specific subsp. morus. c, d, e specific subsp. morus + sandyi and/or sandyi-like. f specific subsp. morus + multiplex. g, h specific subsp. morus + fastidiosa (with/without CFBP 8073 strain). i specific subsp. morus + fastidiosa (with CFBP 8073) + subsp. sandyi + subsp. sandyi-like. j specific subsp. morus + pauca. Frequency of k-mers are mapped onto the genome of reference for subsp. morus (Mul0034)

Full size image

A focus on these 28 categories was performed. Due to the redundancy within the GO hierarchical nomenclature, the 28 GOs were reduced to six. All the CDS harboring specific long-mers in Mul0034 genome were retrieved, as well as their closest homologs in the in- (subsp. fastidiosa, sandyi, sandyi-like, morus) and the out- (subsp. multiplex or pauca) groups. This corresponded to 26 CDS found each in a single copy in all the 46 Xf genomes, indicating that they belong to the Xf core genome. For each gene, the sequences were aligned together with the specific long-mers (80 long-mers in total), against the sequence in Mul0034 used as a reference. First, perfect identity in long-mer sequence was conserved among all the genomes of subsp. fastidiosa, sandyi, sandyi-like, and morus. In contrast, SNPs were always found in the alignment of multiplex and pauca sequences. In comparison with multiplex, 46/80 long-mers had only synonymous SNPs and 34 had non-synonymous SNPs in the gene sequences. Considering the 26 CDS, 18 harbored non-synonymous SNPs. In comparison with pauca, half long-mers had only synonymous SNPs and half had non-synonymous SNPs in the gene sequences. Considering the 26 CDS, 21 harbored non-synonymous SNPs (Additional file 8).

Genetic diversity within the subspecies pauca

ANIb values clearly showed genetic heterogeneity among strains of the subspecies pauca. Three lineages were differentiated: subclade I.1 included seven citrus strains (9a5c, U24D, Fb7, CVC0251, CVC0256, J1a12 and 11,399) and two coffee strains (3124 and 32); subclade I.2 included two coffee (6c and COF0324) and one Prunus (Pr8x) strains, and subclade I.3 included the ST53 strains CoDiRO, COF0407, OLS0478 and OLS0479. Two strains, CFBP 8072 and Hib4 were isolated, as they appeared outside (Additional file 1).

A search for specific mers within these three subclades (Table 6C, Additional file 6) and the identification of associated GO terms (Table 8; Additional file 7) were performed. At the level of subspecies pauca it concerned 175 GO terms, including 16 terms related to nucleotide metabolic/biosynthetic process, especially for purine and 20 terms associated with the bacterial cell wall/envelope/plasma membrane. A focus on these 20 categories was performed. Due to redundancy within the GO hierarchical nomenclature, the 20 GOs were reduced to eight. All the CDS harboring specific long-mers in 9a5c genome were retrieved, as well as their closest homologs in the in- (subsp.pauca) and the out- (subsp. multiplex or fastidiosa, sandyi, sandyi-like, and morus) groups. This correspond to 105 CDS found each in a single copy in almost all the 46 Xf genomes. For each gene, the sequences were aligned together with the specific long-mers (746 k-mers in total), against the sequence in 9a5c used as a reference. First, perfect identity in long-mer sequence was conserved among all the genomes of subsp. pauca, except for 6 long-mers with small variants in a few pauca strains. In contrast, SNPs were always found in the alignment with multiplex and fastidiosa, sandyi, sandyi-like, morus sequences. In comparison with multiplex, 390/746 long-mers had only synonymous SNPs, 354 had non-synonymous SNPs and 2 long-mers were not found. Considering the 105 CDS, 93 harbored non-synonymous SNPs. In comparison with fastidiosa, sandyi, sandyi-like, and morus, 368 long-mers had only synonymous SNPs and 378 had non-synonymous SNPs in the gene sequences. Considering the 105 CDS, 92 harbored non-synonymous SNPs (Additional file 8). The search for enriched GO was also performed at the subclade level. For the subclade I.1, only one term was found (catalytic activity, GO:0003824). In all other studied cases (subclades or individual strains within subsp. pauca), all the GO terms identified were under-represented. This might be explained by a recent evolution in these strains/subclades, rendering the corresponding SNPs in these functional gene ontologies less frequent. Two regions harboring specific mers in 13 pauca strains (subclades I.2 and I.3, plus Hib4) and absent in the others (subclade I.3 plus CFBP 8072) were identified (Fig. 4). These include genes encoding various enzymes (endonuclease, hydrogenase, hydrolase, integrase/recombinase, methyltransferase, peptidase, polyketide synthase, reductase, terminase, topoisomerase) and genes associated with bacteriophages (Additional file 6). For subclade I.3, 10 GO terms were specific, dealing with transport, recombination, and organelle part. CFBP 8072 has two specific GO terms nucleoside triphosphate biosynthetic process and monocarboxylic acid biosynthetic process. As for the Hib4 genome, 12 terms were found as unique, associated with iron-sulfur complex, or transport activity.

Fig. 4
figure 4

Distribution of k-mers specific to the X. fastidiosa subspecies pauca and its subclades. a core k-merome of X. fastidiosa species. b specific subsp. pauca. c specific of subclade I.2, I.3 and strain Hib4 from subsp. pauca. Frequency of k-mers are mapped onto the genome of reference for subsp. morus (Mul0034)

Full size image

Identification of chromosome and plasmid specific islands unique in X. fastidiosa Hib4

The k-mer approach identified three large genomic regions that were specific to the strain Hib4. A fragment of 34,148 bp (long-mer2422 in Additional file 6) appeared to be chromosomic. It contained 34 genes coding for 12 hypothetical/conserved proteins, 8 conjugal transfer proteins (including TraG), 3 membrane proteins, 2 methyltransferases, and one acriflavine resistance protein B, DEAD/DEAH box helicase, DSBA oxidoreductase, hemolysin secretion protein D, integrating conjugative element protein pill (pfgi-1), lytic transglycosylase, multidrug transporter, RAQPRD family plasmid, and superoxide dismutase (Additional file 6). The screen (blastn) of the Nucleotide collection (nr/nt) database revealed that it shared high identities (> 90%) with sequences of Cupriavidus sp., Comamonas testosterone, Pseudomonas aeruginosa, Klebsiella pneumoniae and Bordetella petrii, but the largest fragments cover no more than 60% of the X. fastidiosa long-mer (Additional file 9).

Two large regions of 32,804 bp (long-mer4538 in Additional file 6) and 16,015 bp (long-mer4596) localized onto the plasmid pXF64-HB. Together with smaller specific long-mers (long-mer4536 to 4596), they accounted for 60,224 bp over the 64,251 bp total size of this plasmid https://www.ncbi.nlm.nih.gov/nuccore/NZ_CP009886.1). It contained 39 genes, including genes coding hypothetical proteins (23), conjugal transfer proteins (7; TraH, I, J, K, N, Q, U, W), and one DNA topoisomerase, endonuclease, helicase, lytic transglycosylase, membrane protein, mobilization protein, protein mobD, relaxase, and TrbA (Additional file 6). The plasmid could have been acquired from a strain of Paraburkholderia hospita (93% identity over 86% length of the plasmid), P. aromaticivorans (86% identity over 83% length) or even Burkholderia vietnamiensis (81% identity over 76% length) or Xanthomonas euvesicatoria (80% identity over 72% length) (Additional file 9).

Robust whole genome-based X. fastidiosa clustering with shared k-mers

After looking at specific k-mers in whole genome sequences using SkIf, we employed a complementary approach to draw a robust image of the genetic relationships among individuals, based on shared k-mers. Simka [23] provided a distance matrix that was transformed in a similarity matrix corresponding to the percent of shared k-mers to assess strain relationships (Additional file 10). The k-mer-based dendrogram showed a general distribution into three major clades, represented by the subspecies pauca (clade I), multiplex (clade II), and the union of subspecies fastidiosa¸ sandyi and morus (clade III; Fig. 5). It is congruent with the one obtained with ANIb illustrated by with a strong linear regression (r2 = 0.9945; Fig. 6). The current clustering of X. fastidiosa in five subspecies should be restricted to three subspecies, a proposal that is supported by ANIb and shared k-mers values (99.00% and 0.86, respectively) for clade III (Fig. 6, Additional files 1 and 10). This mostly differ from the view obtained with a MLSA scheme (7 genes) by the repositioning of the subspecies morus (Fig. 5). We finally mapped the key points resulting from SkIf (specific k-mers) analysis on the dendrogram (shared k-mers) to illustrate how X. fastidiosa genetic diversity can be associated with particular traits (Fig. 5).

Fig. 5
figure 5

Phylogenetic representation of X. fastidiosa using k-mers, ANIb and MLSA schemes. All the representations were constructed using the 46 X. fastidiosa, with addition of the X. taiwanensis genome sequence. a Whole genome-based dendrogram built with distance matrixes obtained after running simka (shared k-mers of 22 nucleotides) or ANIb (1020 nt) algorithms. Some specificities and similarities in enriched gene ontologies or identification of plasmid and chromosomic sequences specific to X. fastidiosa are highlighted at nodes or subclades. b Maximum-Likelihood (ML) tree constructed with 1000 replicates for bootstrap values using the concatenated sequences (4161 bp) of seven housekeeping genes from a MultiLocus Sequence Analysis (MLSA) scheme. Key features related to X. fastidiosa subspecies obtained through the combination of specific k-mer identified and gene ontologies enrichment tests are indicated at nodes

Full size image
Fig. 6
figure 6

Inter- and intrasubspecies comparisons of ANIb and shared k-mers values. a Boxplot of the ANIb values calculated from our genome dataset. b Boxplot of the shared k-mer values. c Dot plot of the ANIb and shared k-mer mean values. Linear regression and its corresponding r2 is indicated. For intrasubpecies comparisons, the number of plotted values corresponds to [(number of genome) 2 - number of genome]. For intersubspecies comparisons, it corresponds to [(2 * number of genome subspecies A * number of genome subspecies B)]. Number of genomes: fastidiosa (13), morus (2), sandyi (3), multiplex (10), pauca (18)

Full size image

Discussion

While tools based on k-mers are mainly used to improve genome assembly [36, 37], SkIf (https://sourcesup.renater.fr/wiki/skif/) was developed to quickly extract information from genomic datasets from already assembled genomes. It allows to decipher genomic fragments associated with traits shared by a group of sequences of interest. This strategy is applicable to any scientific questions requesting the comparison of user-defined groups of sequences.

Because management of X. fastidiosa outbreaks in France depends on the subspecies of X fastidiosa, it is of major importance to precisely define these subspecies, understand the robustness of these groupings and their meaning in terms of shared and specific genetic material. In order to detect genomic regions specifically associated with a group of organisms (i.e. a subspecies) we applied SkIf to gain a better understanding of X. fastidiosa clusterings in subspecies. This tool was also used to mine large databases as a first step to evaluate worldwide dispersion of X. fastidiosa in natural settings.

The phylogeny provided by shared k-mers was highly similar to the one based on ANIb, a reference method for analyzing phylogeny of bacteria [38, 39]. However, phylogenies were much more quickly constructed using shared k-mers than were ANI calculations in JSpecies. Here, k-mers of 22 nt were used while ANIb and TETRA are calculated from k-mers of 1020 and 4 nt, respectively, and ANIm values result from the maximal unique match decomposition of two genomes [40,41,42].

The current grouping of X. fastidiosa in five subspecies is inappropriate and is not supported by genomic data. This is obvious regarding the phylogenies reconstructed using shared-k-mers and ANIb (Fig. 5) and it is coherent with a previous proposal [43]. Three-well demarcated genomic clusters were retrieved in phylogenetic trees reconstructed from 46 genome sequences. X. fastidiosa subsp. fastidiosa embraced, in addition to the classical subsp. fastidiosa strains, the more recently proposed sandyi and morus subspecies. Mean ANIb values of 99% are found within the former subspecies, while ANIb value with the two later are below 98% (Fig. 6). The two-other subspecies, multiplex and pauca, were well supported, even if subsp. pauca showed a clear divergence between the lineage of strains isolated from citrus and coffee in Brazil vs. the lineage of other strains isolated from coffee from Central America and olive. Indeed, mean ANIb value of 99.44% was calculated within the multiplex subspecies, while ANIb values below than 97% were found with the other two subspecies. Concerning pauca, mean ANIb value of 98.48% was calculated for the 18 genome sequences included in this subspecies while ANIb values of less than 97% were obtained with the two-other species. This mean ANIb value of only 98.66% within pauca clearly illustrate the largest diversity found in this subspecies in comparison to the fastidiosa and multiplex ones. This grouping in a subsp. fastidiosa sensu largo also matches with an enrichment in 28 GO terms associated with various processes like cellular component or protein complex disassembly, peptidyl-proline activity and chemotaxis. Interestingly, GO-enriched associated CDSs, which harbor k-mers specific to this clade, have homologs in all other X. fastidiosa genomes, but most often these homologs present non-synonymous SNPs, avoiding perfect matches with the k-mers. More importantly, this suggests diversity in protein sequences with putative impact on their functions. Referring to the definition of the species and a threshold value of ANI at 95% [39] values calculated here on 46 genomes sequences indicate that X. fastidiosa with its diversity forms a unique species.

Grouping the subspecies morus within a subspecies fastidiosa sensu largo is coherent with the results of the analysis made to uncover the origin of the subsp. morus. This subspecies was proposed to group strains pathogenic on mulberry trees that derived from recombination events between ancestors of the subspecies fastidiosa and multiplex [16]. The use of SkIf showed that the cumulated size of the mers uniquely shared within genomes of the clade III (subsp. morus, fastidiosa, sandyi and relatives) represents 5% of the Mul0034 genome, while those uniquely shared between subsp. morus and multiplex count for 3.7%. That showed a closest proximity of morus with fastidiosa and sandyi strains than with multiplex. But because some mers are uniquely associated with the two morus genomes (Additional file 6), some genetic material of an unknown origin has been introduced in morus subspecies genome during evolutionary history.

The evolutionary history of X. fastidiosa is also driven by the acquisition of genetic material from heterologous origin. Recently the genome sequence of Hib4 strain was released. This strain presents three large genomic fragments that are unique within X. fastidiosa. One of these regions is chromosomic, the two others locate on a large plasmid (~64kbp) that is not found in any other X. fastidiosa genome sequence but share strong homology with plasmids from Burkholderia hospita DSM17164 [44], P. aromaticivorans BN5 [45], Burkholderia vietnamensis G4 [46], or Xanthomonas euvesicatoria LMG930 [47] strains (Additional file 9). It should be mentioned that so far it is the only case of a plasmid that is not distributed in various strains within X. fastidiosa and that originates from a non-Xylella strain. Thus, it is tempting to hypothesize on how the acquisition could have occurred, while not easy as these species were isolated from various natural environments including water, soil and plants. Interestingly, X. euvesicatoria LMG932 was isolated from Capsicum frutescensi in Brazil [48], the country of origin of X. fastidiosa Hib4 strain isolated from Hibiscus fragilis. Several strains of Burkholderia vietnamensis were isolated from Coffee plants in Mexico [49], while other were isolated from a Brazilian cystic fibrosis patient [50]. These findings illustrate the presence of putative plasmid donor either in the country (Brazil) were Hib4 was isolated and on a host (coffee) in a country (Mexico) were X. fastidiosa is known to occur. Another particularity of Hib4 is that its harbors two specific regions, only shared with strains of the subsp. pauca subclades I.1 and I.2 (Fig. 4). These genomic regions presumably result from a bacteriophage origin. They could have been acquired specifically by a common ancestor of subclades I.1, I.2 and Hib4 or they could have been lost during evolution, accentuating the degree of divergence with other pauca strains.

Microbial collections exchange strains, but mistakes during collection curation cannot be totally excluded, engendering distribution of mislabeled strains, as already shown for X. fastidiosa [51]. SkIf proved useful to survey microbial collections for synonyms. Here, following the release of the genome sequence of the X. fastidiosa type strain from three origins (CFBP 7970 in the present study, ATCC 35879, DSM 10026), SkIf was used to check the relevance of their synonymy. While not strictly identical due to sequencing and assembly biases, the most striking feature was the putative absence of a large fragment in ATCC 35879, corresponding to a plasmid carrying a complete type IV secretion system [34]. Yet, based on our analysis, it is possible that the plasmid could be present in ATCC 35879, but could have been partially lost during the read filtering and assembly process, or even during strain cultivation. The synonymy between the three strains would therefore be valid. An alternative scenario might be that the plasmid-like sequences have been integrated into the chromosome of ATCC 35879 whilst the plasmid was lost, and in this case, only CFBP 7970 and ATCC 35879 are indeed synonymous. The definite answer will come from an analysis of the raw reads of ATCC 35879 to check for the presence/absence of the plasmid (raw reads are currently not available in SRA) or from a plasmid extraction from the specimen stored at ATCC.

Occurrences of X. fastidiosa found in Silva rRNA database were assigned to subspecies that are coherent with sample designation. Specific mers of the genus Xylella and of the various subspecies of X. fastidiosa were retrieved in the V3-V4 region of the 16S rRNA encoding genes. One use of these tools is to taxonomically assign at the subspecies level the occurrences of X. fastidiosa from large database. The sample description and especially the name of the plant species of isolation allow to validate the assignation provided by specific mers. It should however be noticed, that some plant species like almond, olive tree, oleander, coffee tree, or citrus may be host of several X. fastidiosa subspecies (www.pubmlst.org/xfastidiosa) [52] and in consequence this a posteriori validation will not always be possible. Long read sequencing will generate more full length 16 s rRNA gene sequences which will facilitate subspecies discrimination. Another tempting use of these tools could be to survey large metagenome database for occurrences of these markers. This is however currently not feasible due to an astonishingly too long time required to download data or incapacity to browse those databases using our markers. Another use will be to design primers and if required probes for PCR detection-identification of X. fastidiosa in plant material.

Conclusions

Skif is a freely available, bioinformatic tool dedicated to the identification of specific mers. Although the results presented here were applied in the context of the emerging plant pathogen Xylella fastidiosa in Europe, this software is useful to answer many other questions beyond this scope. It is adapted to mine various group of sequences (gene, protein, genome, metagenome databases) defined by the user to identify specific or shared features. In the context of X. fastidiosa it allowed to i) refine the current grouping in subspecies that are not supported by genomic data; ii) trace the origin of the subspecies morus, a plasmid from Hib4 strain and the extent of synonymy among specimen representing the same initial strain in microbial collections; and iii) design markers that are specific to each subspecies of X. fastidiosa.

Methods

Bacterial strains and growth conditions

The seven strains of X. fastidiosa (Table 2) used in this study were provided by the French Collection of Plant-Associated Bacteria (CIRM-CFBP; http://www6.inra.fr/cirm_eng/CFBP-Plant-Associated-Bacteria). Strains were grown on B-CYE [53] medium up to 8 weeks at 28 °C. Experiments with X. fastidiosa living cells were carried out under quarantine at IRHS, Centre INRA, Beaucouzé, France under the agreement no. 2013119–0002 from the Prefecture de la Région Pays de la Loire, France.

Genomic DNA extraction

For genome sequencing, bacterial material was harvested on agar plates and suspended in 4.5 ml of sterile, ultrapure water. Genomic DNA was extracted with the NucleoSpin Tissue kit (Macherey-Nagel), following the manufacturer’s recommendations. DNA was recovered in 100 μl of elution buffer (5 mM Tris/HCl, pH 8.5) with final concentration ranging from 3 to 12 μg. Quality and quantity of extracted genomic DNA were checked by depositing an aliquot on agarose gel combined to the use of a nanodrop (Thermo Scientific).

Library preparation

Genomic DNA solutions were homogenized at 20 ng/μl in 55 μl of resuspension buffer to prepare libraries of sonicated, purified, blunted, and adenylated DNA fragments of 350 bp, following the instructions of the Illumina TruSeq DNA PCR-Free Sample Preparation Guide – Low Sample (LS) Protocol (Catalog #FC-121-9006DOC, Part #150361887 Rev. B, November 2013). Adapters were ligated using the Illumina TruSeq DNA Free PCR LT kit. Libraries were individually quantified and then mixed in a single, equimolar pool (40 nM) also quantified by qPCR following the recommendations of the Library Quantification kit (Kapa Biosystems).

Genome sequencing, assembling, and annotation

For sequencing, diluted libraries (4 nM) were denatured as described (Illumina Preparing DNA libraries for Sequencing on the MiSeq protocol), resulting in 20pM denatured DNA. The final DNA concentration used for sequencing was 12pM in a 600 μl volume containing 1% of PhiX control. The sample was deposited in a V3 cartridge. The seven X. fastidiosa genomes were sequenced with the Illumina MiSeq v3 600 cycles technology at the ANAN plateform, SFR QuaSav, Angers, Fr. Genome assembly was performed using a combination of Velvet [54], SOAPdenovo and SOAPGapCloser [55] assemblers. Structural and functional annotations were conducted with Eugene-PP algorithm [56], using a concatenation of the Swissprot database and the publicly available X. fastidiosa 9a5c [57] and Temecula1 [58] genomes.

Definition of the acronyms

In this study, we used the following definitions for these five key terms: i) mer: a sequence within a nucleotide character string; ii) k-mers: all the possible substrings of length k that are contained in a nucleotide character string; iii) long-mer: a result of the concatenation of overlapping and/or consecutive mers; iv) specific mers: sequences that are exclusively found for members of the group of interest (in-group) while small variants (i.e. with a few indels or SNPs) can be found in some members of the out-groups without being strictly identical; and (v) shared mers: sequences that are found in all the members of different groups of interests or that are common to two individuals in the case of pairwise comparisons.

Identification of shared or specific k-mers

The percent of shared k-mers between two genome sequences were calculated from the distance matrix built using Simka [23]. Parameters were selected as follow: “-kmer-size 22”, “-abundance-min 1”. SkIf (v1.2) was developed in C ++ and sequence reading was done using Bio++ bpp-seq library [59]. To identify genomic regions that are specific to a group of sequences of interest, SkIf construct an abundance matrix of all mers of sequences. This matrix is used to identify the mers present in all the sequences of the group of interest (in-group) and absent in all the other sequences. Parameters were selected as follow: “-k 22”, “-a dna”, “-g = in-group list”. Then, it maps the specific mers of the in-group to the reference genome sequence of the group and provides their precise locations. By comparing mer length and the positions of the various occurrences, SkIf concatenates the overlapping mers into long-mer using the script “getLongestKmersNC.pl” (option: -k 22; available with Skif). A list of located mers or long-mers specific to the group of interest was hence obtained. Finally, we developed a wrapper for accessing this process in a user-friendly Galaxy tool (https://iris.angers.inra.fr/galaxypub-cfbp). Hence, SkIf allows to extract all specific mers of a dataset. The optimal size of the mer was fixed to 22 nt to optimize the ratio of in-group to out-group specific sequences, after a comparison of a range from 18 to 26 nt was done (data not shown).

Analyses of genome and nucleotide sequences and phylogeny

For the analyses of the seven housekeeping genes used in MLSA-MLST scheme designed for X. fastidiosa (https://pubmlst.org/xfastidiosa/info/primers.shtml) and the 16S rRNA gene and the synonymy of SNPs, nucleotide sequences were aligned using the Geneious suite, with the default parameters of the ‘MUSCLE Alignment’ and the ‘Map to Reference’ options [60]. Maximum-Likelihood (ML) tree was constructed with 1000 replicates for bootstrap values using the concatenated sequences (4161 bp) of seven housekeeping genes from the MLSA scheme. ANIb values were calculated using Pyani [30]. Similarity matrix (based on ANIb or shared k-mers) were transformed into distance matrix (1-ANIbs*100 or 1-shared k-mers*100) in the dist format of R using as.dist and clustered using Ward’s method [61] for hierarchical clustering. Conversion of the distance matrix into dendrograms relied on as.phylo function from R ape package [62]. Blastn (v2.8.0+) analyses were run against the nucleotide collection (nt; 46,977,437 sequences) [63].

Enrichment tests and Venn diagram representations

Enrichment analyses with a Fisher’s Exact Test were performed with Blast2GO v4.1 [64], using the Gene Ontology functional annotations to compare gene lists carrying specific k-mers against all the gene of the reference genomes to identify statistically significant enrichment in biological processes or molecular functions. Venn diagrams were built using jvenn (http://jvenn.toulouse.inra.fr/app/index.html [65].

Development of a galaxy-based website and user guidelines

The SkIf pipeline is free for use online (https://iris.angers.inra.fr/galaxypub-cfbp). Required input files are a zip file with all the fasta files for the in-group genome sequences; a zip file with all the fasta files for the outgroup genome sequences; the length of the k, and the identifier of the reference genome sequence from the in-group. Output files are text files with the list of the k-mers and long-mers specific to the in-group if existing. A wiki page describing SkIf is accessible at https://sourcesup.renater.fr/wiki/skif.

Mining of 16S rRNA database

For the genus analysis, the in-group included the 74 Xylella 16S rRNA copies and the out-group included all the small subunit rRNA gene (SSU) sequences from the Silva database (https://www.arb-silva.de/; release 128) [66] other than Xylella-tagged. For the species analysis, X. taiwanensis PLS 229 strain was included in the out-group. All sequences affiliated to X. fastidiosa were included in the in-group, while all the other, non-X. fastidiosa, were included in the out-group. Ambiguous sequences (e.g. double assignation to Xylella and Xanthomonas) were excluded. SkIf software was used as described above to identify specific k-mers (−k 22) in the in-group and to concatenate consecutive ones in long-mers.

Abbreviations

ANIb:

Average nucleotide identities using blast

ANIm:

Average nucleotide identities using MUMmer

CDS:

Coding sequence

GO:

Gene ontology

ICE:

Integrative and conjugative elements

rRNA:

Ribosomal RNA

SkIf:

Specific k-mers Identification

SNP:

Single nucleotide polymorphism

TETRA:

Tetranucleotide frequency correlation coefficients

References

  1. Saponari M, Boscia D, Nigro F, Martelli GP. Identification of DNA sequences related to Xylella fastidiosa in oleander, almond and olive trees exhibiting leaf scorch symptoms in Apulia (southern Italy). J Plant Pathol. 2013;95(3):668.

    Google Scholar 

  2. Denancé N, Legendre B, Briand M, Olivier V, de Boisséson C, Poliakoff F, et al. Several subspecies and sequence types are associated to the emergence of Xylella fastidiosa in natural settings in France. Plant Pathol. 2017;66:1054–64.

    Article  CAS  Google Scholar 

  3. European and Mediterranean Plant Protection Organization (EPPO). First report of Xylella fastidiosa subsp. fastidiosa on Nerium oleander in Germany. EPPO Reporting Service. 2016;7:133.

    Google Scholar 

  4. European and Mediterranean Plant Protection Organization (EPPO). First report of Xylella fastidiosa in Spain. EPPO Reporting Service. 2016;11:213.

    Google Scholar 

  5. European and Mediterranean Plant Protection Organization (EPPO). Xylella fastidiosa detected in mainland Spain and update for Baleares. EPPO Reporting Service. 2017;7:133.

    Google Scholar 

  6. Strona G, Carstens CJ, Beck PSA. Network analysis reveals why Xylella fastidiosa will persist in Europe. Sci Rep. 2017;7(71):1–8.

    CAS  Google Scholar 

  7. Soubeyrand S, de Jerphanion P, Martin O, Saussac M, Manceau C, Hendrikx P, et al. Inferring pathogen dynamics from temporal count data: the emergence of Xylella fastidiosa in France is probably not recent. New Phytol. 2018;219:824–36.

    Article  PubMed  PubMed Central  Google Scholar 

  8. Denancé N, Cesbron S, Briand M, Rieux A, Jacques M-A. 2017. Is Xylella fastidiosa really emerging in France ? In: Costa J, Koebnik R, eds. 1st Annual Conference of the EuroXanth – COST Action Integrating Science on Xanthomonadaceae for integrated plant disease management in Europe. Dec. 13–15, Coimbra, Portugal: EuroXanth, 7.

  9. Bergsma-Vlami M, van de Bilt JLJ, Tjou-Tam-Sin NNA, Helderman CM, Gorkink-Smits PPMA, Landman NM, et al. Assessment of the genetic diversity of Xylella fastidiosa in imported ornamental Coffea arabica plants. Plant Pathol. 2017;66:1065–74.

    Article  CAS  Google Scholar 

  10. European and Mediterranean Plant Protection Organization (EPPO). Xylella fastidiosa detected in Coffea spp. plants imported into Switzerland. EPPO Reporting Service. 2015;10:181.

    Google Scholar 

  11. Jacques MA, Denancé N, Legendre B, Morel E, Briand M, Mississipi S, et al. New coffee-infecting Xylella fastidiosa variants derived via homologous recombination. Appl Environ Microbiol. 2016;82(5):1556–68.

    Article  CAS  PubMed Central  Google Scholar 

  12. Bergsma-Vlami M, van de Bilt JLJ, Tjou-Tam-Sin NNA, van de Vossenberg BTLH, Westenberg M. Xylella fastidiosa in Coffea arabica ornamental plants imported from Costa Rica and Honduras in the Netherlands. J Plant Pathol. 2015;97:395.

    Google Scholar 

  13. Su CC, Deng WL, Jan FJ, Chang CJ, Huang H, Shih HT, et al. Xylella taiwanensis sp. nov., causing pear leaf scorch disease. Int J Syst Evol Microbiol. 2016;66(11):4766–71.

    Article  CAS  PubMed  Google Scholar 

  14. Nunney L, Yuan X, Bromley RE, Stouthamer R. Detecting genetic introgression: high levels of intersubspecific recombination found in Xylella fastidiosa in Brazil. Appl Environ Microbiol. 2012;78(13):4702–14.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Nunney L, Ortiz B, Russell SA, Ruiz-Sanchez R, Stouthamer R. The complex biogeography of the plant pathogen Xylella fastidiosa: genetic evidence of introductions and subspecific introgression in Central America. PLoS One. 2014;9:e112463.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  16. Nunney L, Schuenzel EL, Scally M, Bromley RE, Stouthamer R. Large-scale intersubspecific recombination in the plant-pathogenic bacterium Xylella fastidiosa is associated with the host shift to mulberry. Appl Environ Microbiol. 2014;80(10):3025–33.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  17. Loconsole G, Saponari M, Boscia D, D’Attoma G, Morelli M, Martelli GP, et al. Intercepted isolates of Xylella fastidiosa in Europe reveal novel genetic diversity. Eur J Plant Pathol. 2016;146(1):85–94.

    Article  CAS  Google Scholar 

  18. Coletta-Filho HD, Francisco CS, Lopes JRS, Muller C, Almeida RPP. Homologous recombination and Xylella fastidiosa host–pathogen associations in South America. Phytopathol. 2017;107:305–12.

    Article  Google Scholar 

  19. Rizk G, Lavenier D, Chikhi R. DSK: k-mer counting with very low memory usage. Bioinformatics. 2013;29(5):652–3.

    Article  CAS  PubMed  Google Scholar 

  20. Ghandi M, Lee D, Mohammad-Noori M, Beer MA. Enhanced regulatory sequence prediction using gapped k-mer features. PLoS Comput Biol. 2014;10(7):e1003711.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  21. Melsted P, Halldórsson BV. KmerStream: streaming algorithms for k -mer abundance estimation. Bioinformatics. 2014;30(24):3541–7.

    Article  CAS  PubMed  Google Scholar 

  22. Abo RP, Ducar M, Garcia EP, Thorner AR, Rojas-Rudilla V, Lin L, et al. BreaKmer: detection of structural variation in targeted massively parallel sequencing data using kmers. Nucl. Acids Res. 2015;43(3):e19.

    Article  PubMed  CAS  Google Scholar 

  23. Plaza Onate F, Batto J-M, Juste C, Fadlallah J, Fougeroux C, Gouas D, et al. Quality control of microbiota metagenomics by k-mer analysis. BMC Genomics. 2015;16:183.

  24. Mapleson D, Garcia Accinelli G, Kettleborough G, Wright J, Clavijo BJ. KAT: a K-mer analysis toolkit to quality control NGS datasets and genome assemblies. Bioinformatics. 2017;33(4):574–6.

    CAS  PubMed  Google Scholar 

  25. Marinier E, Zaheer R, Berry C, Weedmark KA, Domaratzki M, Mabon P, et al. Neptune: a bioinformatics tool for rapid discovery of genomic variation in bacterial populations. Nucl. Acids Res. 2017;45(18):e159.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Pandey P, Bender MA, Johnson R, Patro R. Squeakr: an exact and approximate k-mer counting system. Bioinformatics. 2018;34(4):568–75.

    Article  CAS  PubMed  Google Scholar 

  27. Hasman H, Saputra D, Sicheritz-Ponten T, Lund O, Svendsen CA, Frimodt-Møller N, et al. Rapid whole-genome sequencing for detection and characterization of microorganisms directly from clinical samples. J Clin Microbiol. 2014;52:139–46.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Chattaway MA, Schaefer U, Tewolde R, Dallman TJ, Jenkins C. Identification of Escherichia coli and Shigella species from whole-genome sequences. J Clin Microbiol. 2017;55(2):616–23.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Richter M, Rosselló-Móra R, Glöckner FO, Peplies J. JSpeciesWS: a web server for prokaryotic species circumscription based on pairwise genome comparison. Bioinformatics. 2016;32(6):929–31.

    Article  CAS  PubMed  Google Scholar 

  30. Pritchard L, Glover RH, Humphris S, Elphinstone JG, Toth IK. Genomics and taxonomy in diagnostics for food security: soft-rotting enterobacterial plant pathogens. Anal Methods. 2016;8:12–24.

    Article  Google Scholar 

  31. Schaad NW, Postnikova E, Lacy G, Fatmi M, Chang CJ. Xylella fastidiosa subspecies: X. fastidiosa subsp. piercei subsp. nov., X. fastidiosa subsp. multiplex subsp. nov., and X. fastidiosa subsp. pauca subsp. nov. Syst Appl Microbiol. 2004;27:290–300 (Erratum, 27:763.).

    Article  CAS  PubMed  Google Scholar 

  32. Verslyppe B, De Smet W, De Baets B, De Vos P, Dawyndt D. StrainInfo introduces electronic passports for microorganisms. Syst Appl Microbiol. 2014;37:42–50.

    Article  PubMed  Google Scholar 

  33. Giampetruzzi A, Chiumenti M, Saponari M, Donvito G, Italiano A, Loconsole G, et al. Draft genome sequence of the Xylella fastidiosa CoDiRO strain. Genome Announc. 2015;3:e01538–14.

    PubMed  PubMed Central  Google Scholar 

  34. Rogers EE, Stenger DC. A conjugative 38 kb plasmid is present in multiple subspecies of Xylella fastidiosa. PLoS One. 2012;7:e52131.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Giampetruzzi A, Saponari M, Almeida RPP, Essakhi S, Boscia D, Loconsole G, et al. Complete genome sequence of the olive-infecting strain Xylella fastidiosa subsp. pauca De Donno. Genome Announc. 2017;6(5):e00569–17.

    Google Scholar 

  36. Al-Okaily A. HGA: denovo genome assembly method for bacterial genomes using high coverage short sequencing reads. BMC Genomics. 2016;17:193.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  37. Powers JG, Weigman VJ, Shu J, Pufky JM, Cox D, Hurban P. Efficient and accurate whole genome assembly and methylome profiling of E coli. BMC Genomics. 2013;14:675.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Konstantinidis KT, Tiedje JM. Genomic insights that advance the species definition for prokaryotes. Proc Natl Acad Sci U S A. 2005;102:2567–72.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Richter M, Rossello-Mora R. Shifting the genomic gold standard for the prokaryotic species definition. Proc Natl Acad Sci U S A. 2009;106:19126–31.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Delcher AL, Kasif S, Fleischmann RD, Peterson J, White O, Salzberg SL. Alignment of whole genomes. Nucleic Acids Res. 1999;27:2369–76.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Teeling H, Waldmann J, Lombardot T, Bauer M, Glöckner FO. TETRA: a web-service and a stand-alone program for the analysis and comparison of tetranucleotide usage patterns in DNA sequences. BMC Bioinformatics. 2004;26(5):163.

    Article  CAS  Google Scholar 

  42. Goris J, Konstantinidis KT, Klappenbach JA, Coenye T, Vandamme P, Tiedje JM. DNA-DNA hybridization values and their relationship to whole-genome sequence similarities. Int J Syst Evol Microbiol. 2007;57:81–91.

    Article  CAS  PubMed  Google Scholar 

  43. Marcelletti S, Scortichini M. Genome-wide comparison and taxonomic relatedness of multiple Xylella fastidiosa strains reveal the occurrence of three subspecies and a new Xylella species. Arch Microbiol. 2016;198(8):803–12.

    Article  CAS  PubMed  Google Scholar 

  44. Goris J, Dejonghe W, Falsen E, De Clerck E, Geeraerts B, Willems A, et al. Diversity of Transconjugants that acquired plasmid pJP4 or pEMT1 after inoculation of a donor strain in the A- and B-horizon of an agricultural soil and description of Burkholderia hospita sp. nov. and Burkholderia terricola sp. nov. Syst Appl Mircobiol. 2002;25(3):340–52.

    Article  CAS  Google Scholar 

  45. Lee Y, Jeon CO. Paraburkholderia aromaticivorans sp. nov., an aromatic hydrocarbon-degrading bacterium, isolated from gasoline-contaminated soil. Int J Syst Evol Microbiol. 2018;68(4):1251–7.

    Article  PubMed  Google Scholar 

  46. Nelson MJK, Montgomery SO, O’Neill EJ, Pritchard PH. Aerobic metabolism of trichloroethylene by a bacterial isolate. Appl Environ Microbiol. 1986;52(2):383–4.

    CAS  PubMed  PubMed Central  Google Scholar 

  47. Jones JB, Lacy GH, Bouzar H, Stall RE, Schaad NW. Reclassification of the Xanthomonads Associated with Bacterial Spot Disease of Tomato and Pepper. Syst Appl Microbiol. 2004;27:755–62 [Erratum: Syst. Appl. Microbiol. 2006;29:85–86.

    Article  CAS  PubMed  Google Scholar 

  48. Albuquerque P, Caridade CMR, Rodrigues AS, Marcal ARS, Cruz J, Cruz L, et al. Evolutionary and experimental assessment of novel markers for detection of Xanthomonas euvesicatoria in plant samples. PLoS ONE. 2012;7(5):e37836.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Estrada-De Los Santos P, Bustillos-Cristales R, Caballero-Mellado J. Burkholderia, a genus rich in plant-associated nitrogen fixers with wide environmental and geographic distribution. Appl Environ Microbiol. 2001;67(6):2790–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Carvalho GM, Carvalho AP, Folescu TW, Higa L, Teixeira LM, Plotkowski MC, et al. Transient isolation of Burkholderia multivorans and Burkholderia cenocepacia from a Brazilian cystic fibrosis patient chronically colonized with Burkholderia vietnamiensis. J Cyst Fibros. 2005;4(4):267–70.

    Article  PubMed  Google Scholar 

  51. Nunney L, Elfekih S, Stouthamer R. The importance of multilocus sequence typing: cautionary tales from the bacterium Xylella fastidiosa. Phytopathol. 2012;102:456–60.

    Article  CAS  Google Scholar 

  52. Yuan X, Morano L, Bromley R, Spring-Pearson S, Stouthamer R, Nunney L. Multilocus sequence typing of Xylella fastidiosa causing Pierce's disease and oleander leaf scorch in the United States. Phytopathol. 2010;100:601–11.

    Article  CAS  Google Scholar 

  53. Wells JM, Raju BC, Nyland G, Lowe SK. Medium for isolation and growth of bacteria associated with plum leaf scald and phony peach diseases. Appl Environ Microbiol. 1981;42:357–63.

    CAS  PubMed  PubMed Central  Google Scholar 

  54. Zerbino DR, Birney E. Velvet: algorithms for de novo short read assembly using de Bruijn graphs. Genome Res. 2008;18:821–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Luo R, Liu B, Xie P, Li Z, Huang W, Yuan J, et al. SOAPdenovo2: an empirically improved memory-efficient short-read de novo assembler. Gigascience. 20112;1:18.

  56. Sallet E, Gouzy J, Schiex T. EuGene-PP: a next generation automated annotation pipeline for prokaryotic genomes. Bioinformatics. 2014;30:2659–61.

    Article  CAS  PubMed  Google Scholar 

  57. Simpson AJ, Reinach FC, Arruda P, Abreu FA, Acencio M, Alvarenga R, et al. The genome sequence of the plant pathogen Xylella fastidiosa. Nature. 2000;406:151–9.

    Article  CAS  PubMed  Google Scholar 

  58. Van Sluys MA, de Oliveira MC, Monteiro-Vitorello CB, Miyaki CY, Furlan LR, Camargo LE, et al. Comparative analyses of the complete genomes sequences of Pierce’s disease and citrus variegated chlorosis strains of Xylella fastidiosa. J Bact. 2003;185:1018–26.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  59. Guéguen L, Gaillard S, Boussau B, Gouy M, Groussin M, Rochette NC, et al. Bio++: efficient extensible libraries and tools for computational molecular evolution. Mol Biol Evol. 2013;30(8):1745–50.

    Article  PubMed  CAS  Google Scholar 

  60. Kearse M, Moir R, Wilson A, Stones-Havas S, Cheung M, Sturrock S, et al. Geneious basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics. 2012;28:1647–9.

    Article  PubMed  PubMed Central  Google Scholar 

  61. Ward JH Jr. Hierarchical grouping to optimize an objective function. J Am Stat Assoc. 1963;58:236–44.

    Article  Google Scholar 

  62. Paradis E, Claude J, Strimmer K. APE: analyses of phylogenetics and evolution in R language. Bioinformatics. 2004;20(2):289–90.

    Article  CAS  PubMed  Google Scholar 

  63. Zhang Z, Schwartz S, Wagner L, Miller W. A greedy algorithm for aligning DNA sequences. J Comput Biol. 2000;7:203–14.

    Article  CAS  PubMed  Google Scholar 

  64. Götz S, García-Gómez JM, Terol J, Williams TD, Nagaraj SH, Nueda MJ, et al. High-throughput functional annotation and data mining with the Blast2GO suite. Nucl Acids Res. 2008;36(10):3420–35.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  65. Bardou P, Mariette J, Escudié F, Djemiel C, Klopp C. jvenn: an interactive Venn diagram viewer. BMC Bioinformatics. 2014;15:293.

    Article  PubMed  PubMed Central  Google Scholar 

  66. Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucl. Acids Res. 2013;41(D1):D590–6.

    Article  CAS  PubMed  Google Scholar 

  67. Zhang S, Flores-Cruz Z, Kumar D, Chakrabarty P, Hopkins DL, Gabriel DW. The Xylella fastidiosa biocontrol strain EB92-1 genome is very similar and syntenic to Pierce's disease strains. J Bacteriol. 2011;193(19):5576–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Schreiber HL, Koirala M, Lara A, Ojeda M, Dowd SE, Bextine B, et al. Unraveling the first Xylella fastidiosa subsp. fastidiosa genome from Texas. Southwest. Entomol. 2010;35(3):479–83.

    Google Scholar 

  69. Chen J, Xie G, Han S, Chertkov O, Sims D, Civerolo EL. Whole-genome sequences of two Xylella fastidiosa strains (M12 and M23) causing almond leaf scorch disease in California. J Bacteriol. 2010;192:4534.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Chen J, Wu F, Zheng Z, Deng X, Burbank LP, Stenger DC. Draft genome sequence of Xylella fastidiosa subsp. fastidiosa strain Stag's leap. Genome Announc. 2016;4:e00240–16.

    PubMed  PubMed Central  Google Scholar 

  71. Van Horn C, Chang CJ, Chen J. De Novo whole-genome sequence of Xylella fastidiosa subsp. multiplex strain BB01 isolated from a blueberry in Georgia, USA. Genome Announc. 2017;5(6):e01598–16.

    PubMed  PubMed Central  Google Scholar 

  72. Bhattacharyya A, Stilwagen S, Ivanova N, D'Souza M, Bernal A, Lykidis A, et al. Whole-genome comparative analysis of three phytopathogenic Xylella fastidiosa strains. Proc Natl Acad Sci U S A. 2002;99:12403–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Chen J, Huang H, Chang C-J, Stenger DC. Draft genome sequence of Xylella fastidiosa subsp. multiplex strain Griffin-1 from Quercus rubra in Georgia. Genome Announc. 2013;1:e00756–13.

    PubMed  PubMed Central  Google Scholar 

  74. Guan W, Shao J, Davis RE, Zhao T, Huang Q. Genome sequence of a Xylella fastidiosa strain causing sycamore leaf scorch disease in Virginia. Genome Announc. 2014;2:e00773–14.

    PubMed  PubMed Central  Google Scholar 

  75. Schuenzel EL, Scally M, Stouthamer R, Nunney L. A multigene phylogenetic study of clonal diversity and divergence in north American strains of the plant pathogen Xylella fastidiosa. Appl Environ Microbiol. 2005;71:3832–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  76. Giampetruzzi A, Loconsole G, Boscia D, Calzolari A, Chiumenti M, Martelli GP, et al. Draft genome sequence of CO33, a coffee-infecting isolate of Xylella fastidiosa. Genome Announc. 2015;3:e01472–15.

    PubMed  PubMed Central  Google Scholar 

  77. Guan W, Shao J, Zhao T, Huang Q. Genome sequence of a Xylella fastidiosa strain causing mulberry leaf scorch disease in Maryland. Genome Announc. 2014;2:e00916–3.

    PubMed  PubMed Central  Google Scholar 

  78. Alencar VC, Barbosa D, Santos DS, Oliveira ACF, de Oliveira RC, Nunes LR. Genomic sequencing of two coffee-infecting strains of Xylella fastidiosa isolated from Brazil. Genome Announc. 2014;2:e01190–13.

    Article  PubMed  PubMed Central  Google Scholar 

  79. Niza B, Merfa MV, Alencar VC, Menegidio FB, Nunes LR, Machado MA, et al. Draft genome sequence of 11399, a transformable citrus-pathogenic strain of Xylella fastidiosa. Genome Announc. 2016;13:e01124–16.

    Google Scholar 

  80. Su CC, Deng WL, Jan FJ, Chang CJ, Huang H, Chen J. Draft genome sequence of Xylella fastidiosa pear leaf scorch strain in Taiwan. Genome Announc. 2014;2:e00166–14.

    PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We thank Muriel Bahut (ANAN technical facility, SFR QUASAV, Angers, FR) for genome sequencing, CIRM-CFBP (Beaucouzé, INRA, France; http://www6.inra.fr/cirm_eng/CFBP-Plant-Associated-Bacteria) for strain preservation and supply, and the CATI BBRIC for the galaxy tools allowing assembling the reads and annotating the sequences, and using Blast2GO. We acknowledge Charles Manceau (Anses, Angers, FR) for his contribution while applying for funding. We thank Matthieu Barret for fruitful discussions and critical reading of the manuscript.

Funding

ND salary was funded by the regional program “Objectif Végétal, Research, Education and Innovation in Pays de la Loire”, project SapAlien 2015–2017, supported by the French Region Pays de la Loire, Angers Loire Métropole and the European Regional Development Fund. This work received support from the European Union’s Horizon 2020 research and innovation program under grant agreement 635646 POnTE (Pest Organisms Threatening Europe). The present work reflects only the authors’ view and the EU funding agency is not responsible for any use that may be made of the information it contains.

Availability of data and materials

Strains of Xylella fastidiosa deposited at the CIRM-CFBP (Angers, FR) are available upon request (https://www6.inra.fr/cirm_eng/CFBP-Plant-Associated-Bacteria). The SkIf pipeline is free to use online (https://iris.angers.inra.fr/galaxypub-cfbp). All the sequence files related to the k-mer analyses (raw data, alignments, blast, GO enrichment tests) are submitted together with the manuscript as Additional files. Genome sequences were deposited at NCBI under the following accessions numbers: PHFQ00000000 (CFBP 7969), PHFR00000000 (CFBP 7970), PHFP00000000 (CFBP8071), PHFS00000000 (CFBP 8078), PHFT00000000 (CFBP 8082), PHFU00000000 (CFBP 8351), and PHFV00000000 (CFBP 8356). Related data and reannotation of the public genomes (fasta, genbank and gff3 formats) are available upon request.

Author information

Author notes

  1. Nicolas Denancé and Martial Briand are contributed equally.

Authors and Affiliations

  1. IRHS, INRA, AGROCAMPUS-Ouest, Université d’Angers, SFR 4207 QUASAV, 42 rue Georges Morel, 49071, Beaucouzé cedex, France

    Nicolas Denancé, Martial Briand, Romain Gaborieau, Sylvain Gaillard & Marie-Agnès Jacques

Authors
  1. Nicolas Denancé
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Martial Briand
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Romain Gaborieau
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Sylvain Gaillard
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Marie-Agnès Jacques
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

ND designed the in silico analyses. MB, RG, and SG designed the bioinformatics tools. ND and MB performed the in silico analyses and interpreted the data. MAJ conceived the study, applied for funding, and interpreted the data. ND and MAJ wrote the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Marie-Agnès Jacques.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

All the authors read and approved the manuscript.

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.

Additional files

Additional file 1:

Pairwise comparison of 47 Xylella sp. genomes using average nucleotide identity based on blast (ANIb). (DOCX 63 kb)

Additional file 2:

K-mers resulted from the comparison of synonymous strains CFBP 7970, DSM 10026 and ATCC 35879. (ZIP 45 kb)

Additional file 3:

Blast analysis of the known X. fastidiosa plasmid sequences against the genome sequence of ATCC 35879. (DOCX 30 kb)

Additional file 4:

Raw data of analyses of the 16S rRNA gene repertoire in Xylella. (ZIP 12 kb)

Additional file 5:

X. fastidiosa 16S rRNA sequences from Silva database carrying the five long-mers and taxonomically assigned to a subspecies with the SNP-based code. (DOCX 21 kb)

Additional file 6:

Raw data of the k-mers identified to be specific of the different X. fastidiosa subspecies or combinations of several subspecies. (ZIP 3012 kb)

Additional file 7:

Raw data of the gene ontologies enrichments tests with Blast2GO. (ZIP 22422 kb)

Additional file 8:

Analysis of (non-)synonymous SNPs in homologs to CDS harboring specific k-mers for selected enriched GOs. (ZIP 549 kb)

Additional file 9:

Blast analysis of the three large fragments specific to X. fastidiosa subsp. pauca Hib4 strain. (DOCX 34 kb)

Additional file 10:

Pairwise comparison of 47 Xylella sp. genomes using the occurrence of shared k-mers of length 22 bp. (DOCX 68 kb)

Rights and permissions

Open Access This 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.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Denancé, N., Briand, M., Gaborieau, R. et al. Identification of genetic relationships and subspecies signatures in Xylella fastidiosa. BMC Genomics 20, 239 (2019). https://doi.org/10.1186/s12864-019-5565-9

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-019-5565-9

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us