Identification of genetic relationships and subspecies signatures in Xylella fastidiosa

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. Electronic supplementary material The online version of this article (10.1186/s12864-019-5565-9) contains supplementary material, which is available to authorized users.


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.
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][34][35]. But the blast analysis suggested its possible presence, as partial 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). 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 LongXy-lefa#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).
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.

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  Codes of strains having one copy of 16S rRNA Codes of strains having two copies of 16S rRNA 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, peptidylproline activity and chemotaxis. 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).
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   (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. 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 (r 2 = 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).

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.

Bacterial strains and growth conditions
The seven strains of X. fastidiosa (

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 Tru-Seq 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.

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.
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.
Authors' 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.
Ethics approval and consent to participate Not applicable.

Consent for publication
All the authors read and approved the manuscript.