Tandem repeats analysis for the high resolution phylogenetic analysis of Yersinia pestis.

BACKGROUND: Yersinia pestis, the agent of plague, is a young and highly monomorphic species. Three biovars, each one thought to be associated with the last three Y. pestis pandemics, have been defined based on biochemical assays. More recently, DNA based assays, including DNA sequencing, IS typing, DNA arrays, have significantly improved current knowledge on the origin and phylogenetic evolution of Y. pestis. However, these methods suffer either from a lack of resolution or from the difficulty to compare data. Variable number of tandem repeats (VNTRs) provides valuable polymorphic markers for genotyping and performing phylogenetic analyses in a growing number of pathogens and have given promising results for Y. pestis as well.
RESULTS: In this study we have genotyped 180 Y. pestis isolates by multiple locus VNTR analysis (MLVA) using 25 markers. Sixty-one different genotypes were observed. The three biovars were distributed into three main branches, with some exceptions. In particular, the Medievalis phenotype is clearly heterogeneous, resulting from different mutation events in the napA gene. Antiqua strains from Asia appear to hold a central position compared to Antiqua strains from Africa. A subset of 7 markers is proposed for the quick comparison of a new strain with the collection typed here. This can be easily achieved using a Web-based facility, specifically set-up for running such identifications.
CONCLUSION: Tandem-repeat typing may prove to be a powerful complement to the existing phylogenetic tools for Y. pestis. Typing can be achieved quickly at a low cost in terms of consumables, technical expertise and equipment. The resulting data can be easily compared between different laboratories. The number and selection of markers will eventually depend upon the type and aim of investigations.

Background

Within the Y. pestis species, strains are separated into three biovars according to their ability to reduce nitrate and to ferment glycerol [1]. Since Y. pestis was connected to plague by Yersin [2], strains of biovar Antiqua have been generally isolated from Asia and from East and Central Africa, Medievalis was found in Central Asia, and Orientalis worldwide. Y. pestis is thought to have recently evolved from Yersinia pseudotuberculosis, some 1,500 to 20,000 years ago [3]. Based on the biochemical assays, Devignat [1] proposed that Antiqua strains, causing probably the first known pandemic, represent the ancestor. Medievalis strains are suggested to be responsible for the second pandemic whereas the third pandemic was associated exclusively with Orientalis strains. This overall scenario, although not formally established, is supported by the observed higher diversity in the Antiqua and Medievalis biovars as measured by IS typing and the geographic origin of strains [3-5]. "Pestoides" strains are particular Y. pestis isolates found recently in Russia, and which infect unique species of rodents [5,6].

Several molecular methods have been used for genotyping Y. pestis strains, mostly based on pulse-field gel electrophoresis (PFGE), insertion sequence polymorphism and ribotyping. Recently, multiple locus variable-number-of-tandem-repeat (VNTR) analysis (MLVA) was shown to be a promising method for genotyping a number of pathogens [7-15]. When applicable, MLVA is of great interest, because data produced in different laboratories can be easily exchanged and merged. This is especially relevant in the case of pathogens such as Y. pestis, for which strains cannot be easily exchanged for security reasons. In this context, the availability of standardized and easy to set-up typing tools to facilitate the research efforts on this pathogen is important. The complete genomic sequences of Y. pestis CO92 [16], biovar Orientalis, and of strain KIM [17], biovar Medievalis, have been determined, facilitating the identification of tandem repeats and consequently the selection of primers for MLVA [8,15].

Until now only small series of Y. pestis strains were typed by MLVA [8,9] and although the method seemed appropriate to genotype reproducibly and accurately, a large scale study seemed necessary. In the present work, a collection of 180 isolates of Y. pestis, from different geographical origins, and including various Y. pestis reference strains were typed using the 25 VNTRs previously described [8].

Results

MLVA genotyping

The 25 loci could be amplified (Figure 1 and Additional file 1) in all 180 isolates (Table 1), with the exception of locus ms06 which failed to yield a PCR amplification product in three strains (corresponding to genotypes 4 and 5, Figure 2) despite numerous attempts. Sixty-one different genotypes were identified in this collection, numbered 1 to 61 (Figure 2). Clustering analysis correctly separates Y. pestis Orientalis isolates comprising genotypes 9 to 51 from Antiqua (genotypes 2, 3 and 5 to 8) and Medievalis (genotypes 4 and 52 to 61). Antiqua strains of essentially African origin (genotypes 2, 3, 5, comprising 6 different strains) and the four Antiqua strains from Russia and Asia (genotype 6 and 7 originating from Russia [18] and genotype 8, Antiqua strains KUMA and Yokohama originating from Manchuria and Japan) are clearly separated. Within biovar Orientalis, a group of isolates among which CO92 (genotype 15) has been assigned to the IS100 typing group O1 using the PCR-based typing method described by Motin et al. [5] (data not shown; results indicated in Figure 2). This demonstrates that the MLVA clustering correlates well with another molecular typing method. Additional PCR-IS100 typing indicates that the other Orientalis strains, mostly from Vietnam, are of the O2a type (Figure 2). This is in agreement with the report from Motin et al. [5] suggesting an association of O2a with South-East Asia. These Orientalis isolates are further separated by MLVA into three main branches comprising genotypes 27 to 51. One representative from each Orientalis genotype (genotypes 9 to 51) was assayed by PCR for the presence of the glpD deletion characterised by Motin et al. All strains yield a PCR product of the size corresponding to the deleted allele (data not shown).

Figure 1

Illustration of the MLVA assay set-up The PCR products of an amplification using primers "ms01" (on the left) or primers "ms04" (on the right) have been loaded on an agarose gel, electrophoresed, and stained with ethidium bromide. Lanes M are a 20 bp ladder molecular weight marker. Samples 1 to 7 correspond to a Yersinia pseudotuberculosis strain (lane 1), Yersinia pestis strains representing genotype 18 in Figure 2 (EV76 DNA from 4 different origins, lanes 2, 4, 5, 7, and a strain from Turkey, lane 6), Yersinia pestis strain KIM (lane 3). The image illustrates how the number of units can be directly deduced by manual reading. The marker names below the gel provide the repeat unit size of the tandem repeat, the expected PCR product size in the CO92 genome, and the corresponding number of units in the CO92 genome.

Table 1

List of the Yersinia pestis isolates used for genotyping

Isolates	Number of isolates	Geographical origin	Biovar
CEB-O-1 to CEB-O-132	132	Vietnam	orientalis
CEB-O-133 and O-134	2	Burma	orientalis
M23, Java9, Java10, TS	4	Indonesia	orientalis
Exu21, Exu38	2	Brazil	orientalis
195P, CEB-O-135 and CEB-O-136	3	India	orientalis
K120	1	SAF	orientalis
EV26, EV76(and 3 replicates), 6/69 (and 2 replicates)	3	Madagascar	orientalis
TW, Yp1255	2		orientalis
GB	1	UK (laboratory)	orientalis
G32, C0-92	2	USA	orientalis
IP523, IP524	2	Senegal	orientalis
Turkey 10/1	1	Turkey	orientalis
CEB-O-137 to CEB-O-139	3	unknown	orientalis
CEB-M-1 to CEB-M-5, PKH4, PKR25, KIM (and 5 replicates)	8	Kurdistan	medievalis
CEB-M-6 and CEB-M-7	2	unknown	medievalis*
KUMA	1	Mandchouria	antiqua
Yokohama	1	Japan	antiqua
5761, 735	2	Russia	antiqua
IP543, 343	2	Congo	antiqua
IP537, A22, 129M22	3	Kenya	antiqua
Turkey 10/5, Turkey 10/3	2	Turkey	antiqua
8786G	1	Georgia	pestoides

*deletion in napA

Figure 2

MLVA genotyping data and cluster analysis Clustering analysis was done using the categorical and UPGMA options. The main biovars are indicated by P (Pestoides), A (Antiqua), M (Medievalis) and O (Orientalis). From left to right, the columns designate the genotype number, the number of isolates with an identical genotype, and the geographic origin of the isolates from each genotype. In three instances, indicated by asterisks (*), the strain origin was unknown. A number of classical strains available for this investigation are indicated in brackets. The Medievalis strain lacking the napA locus is indicated by **. O1 or O2a refers to the nomenclature defined by Motin et al. [5].

Most Medievalis strains are also clustered into one major branch (genotypes 52 to 61) with the exception of the strain representing genotype 4. The Pestoides isolate from Georgia (genotype 1) is (weakly) grouped with the African Antiqua strains. Among Antiqua isolates, the KUMA and Yokohama isolates show the identical MLVA genotype 8, and possess a seemingly specific ms09 allele.

As a complementary analysis, a minimum spanning tree analysis was performed (Figure 3). This kind of analysis is applicable to categorical data sets (see also [19]). The creation of hypothetical types (open circles) further minimizes the summed distance of all branches of the tree. The numbers indicated in the circles correspond to the genotype numbers listed in Figure 2. The Orientalis strains (genotypes 9 to 51) are closely related, and are grouped into a single complex. The O2a strains from Vietnam constitute a well-defined subgroup. The backbone of the O2a cluster is made by genotypes 28, 40, 46, 48, which altogether represent 97 isolates. The two Antiqua strains KUMA and Yokohama from Asia (genotype 8) are clearly located outside of the Orientalis group in this analysis, in a position intermediate between the Medievalis group (genotypes 52 to 61) and the Orientalis group. All strains classified as Medievalis strains based on the nitrate reductase assay fall within this group, with one exception, the strain from Kenya defining genotype 4 which is very close to Antiqua genotype 5 (two strains from Kenya and Congo). The six genotypes representing Antiqua strains (African strains, genotypes 2, 3, 5; Asian strains, genotypes 6 to 8) are very loosely connected to each other, suggesting a very high diversity of this biovar.

Figure 3

Minimum spanning tree analysis A minimum spanning tree was constructed using the genotyping data provided in Additional file 1. The parameters used for the construction of the tree authorized the introduction of hypothetical (missing) links which appear as open circles. A number of reference strains are indicated in brackets.

At least two independent types of Medievalis strains

The strain representing genotype 4 is nitrate-reductase negative, and for this reason has been phenotypically assigned to the biovar Medievalis. However, its position in Figures 2 and 3, being very close to genotype 5 (Antiqua strains) and very distinct from the Medievalis cluster, and its geographic origin from Kenya prompted us to investigate the origin of nitrate-reductase deficiency in more detail. The Medievalis strain KIM is nitrate-reductase deficient because of a single point mutation in the napA gene [17]. We have analyzed all Medievalis strains from genotypes 52 to 61 and they showed the same point mutation (data not shown). In contrast the napA gene in the strain representing genotype 4 has been inactivated by a deletion as seen by the absence of amplification of the gene in spite of the use of different primers from this locus, and the absence of an hybridisation signal in a Southern blot experiment (data not shown).

A simple MLVA assay comprising 7 markers

We have evaluated the possibility to define a smaller set of markers for routine typing and comparison of new strains with the data from the present collection. Table 2 lists the main characteristics of the 25 loci. Their name is indicated in the first column. The second column indicates when relevant the name of the corresponding marker in Klevytska et al [9]. The third and fourth columns indicate the repeat size and the number of alleles in our collection of strains and that of Klevytska et al. [9]. The fifth and sixth columns indicate the size range observed in our collection of strains, and the seventh column contains the polymorphism index of each marker. Seven markers (ms01, ms04, ms06, ms07, ms46, ms62, ms70) have a polymorphism index value above 0.6, seven others have an index value between 0.4 and 0.5, the last eleven have polymorphism index below 0.4. Fifty-seven genotypes are resolved with the 7 most polymorphic loci, instead of 61 when using the full set of 25 markers (analysis not shown).

Table 2

Main characteristics and behavior of the 25 tandem repeat loci

locus	alias	repeat size (bp)	nb of alleles*	allele size range (units)	allele size range (bp)	polymorphism index*
ms01		18	6	6–11	192–282	0.61
ms04	M58	17	5(6)	5–9	179–247	0.71 (0.78)
ms05	M59	17	3(4)	10–12	274–308	0.45 (0.68)
ms06		60	10	3–11	306–786	0.73
ms07	M37	10	4(5)	7–10	164–194	0.65 (0.72)
ms09		18	7	5–36	178–736	0.48
ms15		15	2	9–10	222–237	0.29
ms20	M51	15	3(2)	7.5–9	230–253	0.42 (0.49)
ms21		18	4	6–11	224–314	0.18
ms35		15	2	8–9	204–219	0.35
ms38		16	4	5–8	185–233	0.42
ms40		17	3	7–9	214–248	0.38
ms41		17	3	5–7	183–217	0.38
ms44		17	2	7–8	233–250	0.23
ms45	M42	12	3(2)	6–8	149–173	0.37 (0.18)
ms46		17	9	6.5–14.5	243–380	0.71
ms51		18	2	6–7	182–200	0.12
ms54		22	2	6–7	259–281	0.09
ms56		16	4	6–11	204–284	0.35
ms62	M34	9	13(11)	5–20	222–357	0.88 (0.81)
ms69		16	3	5–7	163–195	0.47
ms70		9	7	3–11	119–191	0.62
ms71		14	4	4–7	143–185	0.45
ms73		18	3	3–6	171–225	0.3
ms74		15	4	5–8	175–220	0.41

*the values reported in [9] are indicated in brackets.

Discussion

In the present report, an MLVA typing assay comprising 25 markers has been applied to a collection of Y. pestis isolates of various origins, but with a strong bias towards Orientalis strains from South-East Asia. One hundred and eighty strains or isolates (Table 1) have been genotyped and 61 different MLVA genotypes were identified (see Additional file 1). Clustering analysis and minimum spanning tree analysis suggest relations between the different strains and biovars which are in excellent agreement with current knowledge. In spite of the very limited number of Antiqua and Medievalis strains which could be investigated here, the data obtained suggest the existence of two groups of Antiqua strains. The first group from Russia and Asia represented by genotypes 6, 7 (from Russia) and 8 (KUMA, Yokohama) holds an intermediate position between the Medievalis and the Orientalis group. The second group comprises the African Antiqua strains. IS100 typing distinguishes the KUMA and Yokohama isolates which were typed as A2 and A1b type, respectively [5]. They are very similar in DNA microarrays studies [20].

Medievalis and Orientalis strains derive from Antiqua strains by the loss of metabolic functions, respectively the capacity to reduce nitrate, and to metabolize glycerol. Whereas all Orientalis strains investigated so far (and including this report) are derived from a single ancestor carrying a simple deletion in the glpD gene [5], we report here that the Medievalis phenotype can be associated with at least two independent mutation events in the napA gene. This underlines the fact that the initial biochemical tests should be complemented, or replaced, by direct molecular analyses of the glpD and napA genes.

It is tempting to speculate about possible scenarios suggested by Figure 3 and current knowledge. Pestoides strains originate from Central Asia (reviewed in [6]), and are proposed to be an outgroup of the Y. pestis group. Their genetic composition is relatively distinct from Y. pestis, in terms of plasmids and chromosome, as assayed by DNA array analysis for instance [20]. In Figure 3, the Pestoides strain studied is very distantly related to an hypothetical missing link close to the center of the figure. Hypothetical missing links (open circles) are created by the Minimum Spanning Tree software if they result in a reduction of the total tree length. Many missing links are suggested by the software in this area. They enable the connection of the three branches made of the African Antiqua group, the Medievalis group, and the Orientalis group. The strains closest to this central position are the Antiqua strains from Asia. This and the position of the Pestoides strain suggest that all three biovars may originate from Asia. Indeed, much more strains from these regions will be needed to test this hypothesis. For instance, strain Nicholinsk51 was shown by [5] to share some features of Orientalis strains while being an Antiqua strain from Asia. The authors favored the hypothesis that this strain was a revertant from the Orientalis phenotype. We would rather predict that Nicholinsk51 will be placed by MLVA analysis between the central group of Figure 3 and the Orientalis group.

In total, in our study and that of Klevytska et al. [9], 61 VNTRs were characterized, and still more exist in the Y. pestis genome [15] that can be tested for the selection of an optimal set for MLVA if necessary. Different questions can be addressed by different sets of tandem repeat loci. Global phylogenetic investigations will best be done using loci with a low or moderate mutation rate. Forensics, or local outbreaks investigations, may use loci with a higher mutation rate and particularly simple sequence repeats such as the tetranucleotide described by Adair et al. [21]. Mutation rate of a tandem repeat locus and phylogenetic value is very poorly predicted from sequence analysis [15,22] so that it has to be experimentally measured by typing collections of isolates, as done here (Table 2 and Additional file 1).

Conclusion

Once set-up, MLVA is a very powerful and reproducible genotyping method and it is hoped that this simple molecular tool will help unravel the molecular phylogeny of Yersinia pestis when being applied to a larger number of isolates. In comparison, MLST analysis [3] proved to be almost non informative within Y. pestis. DNA array analysis [20,23] shows a very low resolution with only a few different genotypes identified so far. IS typing by Southern blotting has a very high discriminatory power, but the resulting data is not easily comparable between different laboratories. PCR-IS typing developed by Motin et al. [5] provides exchangeable data but with a much reduced resolution as compared to classical IS typing. MLVA typing may thus turn out to be the method of choice for Y. pestis, once more isolates will have been typed, common genotype databases put together, and reference collections of markers selected. As one step in this direction, the data has been made accessible on our Web site . This includes not only the full dataset which can be recovered from Additional file 1, but also the possibility to run queries with new MLVA data. This may be of use for investigators lacking the specialized tools or expertise required to run MLVA clustering analyses. A very satisfying typing resolution (but not a robust phylogenetic analysis) can readily be achieved by PCR amplification of only 7 loci.

Methods

Bacterial strains and isolates

Most of the strains are from the collection maintained by the French ministry of defense at Centre d'Etudes du Bouchet (CEB) and others came from French medical military institutions [24]. The Y. pestis strains were isolated mostly in Asia between 1964 and 1979, some in Africa (Kenya, Congo), Kurdistan and Madagascar. Additional reference strains and DNA isolates of different biovars, as identified by the source laboratories, are from the Institute of Microbiology Federal Armed Forces, Munich (Germany), or were kindly provided by Prof. F. Allerberger (Vienna, Austria), Prof. H. Tschäpe (Wernigerode, Germany), Pr H. Mollaret, Dr E. Carniel (Paris, France) and Dr. A. Rakin (Munich, Germany) (Table 1). Thermolysates were prepared by heating a bacterial suspension in water for 30 min at 95°C.

Biovar assignment by molecular tests

The presence of an intact glycerol-3-phosphate dehydrogenase (glpD) gene was tested by PCR as described [5]. Sixty three isolates were typed by PCR for the presence of four IS100 insertions (vlm04 and vlm05, vlm06, vlm25, vlm28) as described by Motin et al. [5]. The single nucleotide change in napA resulting in nitrate reductase deficiency in strain KIM [17] was used in a SNP typing assay. Primers napAFor (5' GCGCTAAAAGAGAAAGGCCCGA 3') and either SNPnapKIM (5' AGAGCACGAAGGCATCGGCTTA 3') or SNPnapCO92 (5' AGAGCACGAAGGCATCGGCTTC 3') were used in two PCR reactions at 58°C annealing temperature producing a 230 bp amplicon in either reaction.

Minisatellite PCR amplification and genotyping

PCR reactions and analyses were performed as described [8] using 25 polymorphic markers (Figure 1). Markers ypms04, ypms05, ypms07, ypms20, ypms45 and ypms62 correspond respectively to marker M58, M59, M37, M51, M42 and M34 in Klevytska et al. [9] (Table 2).

Data management and analyses

Gel images were analyzed using the bionumerics software package version 3.5 (Applied-Maths, Sint-Martens-Latem, Belgium) as previously described [10]. The number of motifs in each allele was deduced from the amplicon size. The resulting data were analyzed with bionumerics as a character dataset. Clustering analysis was done using the categorical parameter and the Ward coefficient. The minimum spanning tree was constructed with the following options: (a) in case of equivalent solutions in terms of calculated distances, the priority rule used was to select the tree with the highest number of branches connecting genotypes differing at only one locus ("Highest number of single locus variants" option).; (b) the creation of hypothetical types (missing links) reducing the total length of the tree was allowed. Polymorphism indexes for each locus were calculated as 1 minus the sum of the squares of the frequency of each allele within the different genotypes identified. The data produced was made accessible from a Web page as previously described [10] taking advantage of the BNServer application (Applied-Maths, Sint-Martens-Latem, Belgium).

Authors' contributions

FAM did most of the typing work and CP and GV did the error checking analysis. CP performed the biovar assignment by molecular tests and the PCR-IS typing. FR was in charge of the CEB strain collection and prepared the DNA samples. HN provided DNA from a number of additional reference strains. GV initiated and managed the project, and was in charge of the Bionumerics database and clustering analyses. CP, HN, and GV wrote the report. All authors read and approved the final manuscript.

Additional File 1

this file contains the genotyping data for each of the 61 different genotypes. It is a tab-delineated text file.

Click here for file