Yersinia pestis Lineages in Mongolia

Background Whole genome sequencing allowed the development of a number of high resolution sequence based typing tools for Yersinia (Y.) pestis. The application of these methods on isolates from most known foci worldwide and in particular from China and the Former Soviet Union has dramatically improved our understanding of the population structure of this species. In the current view, Y. pestis including the non or moderate human pathogen Y. pestis subspecies microtus emerged from Yersinia pseudotuberculosis about 2,600 to 28,600 years ago in central Asia. The majority of central Asia natural foci have been investigated. However these investigations included only few strains from Mongolia. Methodology/Principal Findings Clustered Regularly Interspaced Short Prokaryotic Repeats (CRISPR) analysis and Multiple-locus variable number of tandem repeats (VNTR) analysis (MLVA) with 25 loci was performed on 100 Y. pestis strains, isolated from 37 sampling areas in Mongolia. The resulting data were compared with previously published data from more than 500 plague strains, 130 of which had also been previously genotyped by single nucleotide polymorphism (SNP) analysis. The comparison revealed six main clusters including the three microtus biovars Ulegeica, Altaica, and Xilingolensis. The largest cluster comprises 78 isolates, with unique and new genotypes seen so far in Mongolia only. Typing of selected isolates by key SNPs was used to robustly assign the corresponding clusters to previously defined SNP branches. Conclusions/Significance We show that Mongolia hosts the most recent microtus clade (Ulegeica). Interestingly no representatives of the ancestral Y. pestis subspecies pestis nodes previously identified in North-western China were identified in this study. This observation suggests that the subsequent evolution steps within Y. pestis pestis did not occur in Mongolia. Rather, Mongolia was most likely re-colonized by more recent clades coming back from China contemporary of the black death pandemic, or more recently in the past 600 years.


Introduction
Yersinia (Y.) pestis subspecies pestis is the causative agent of human plague. Cases are annually registered by the WHO and nowadays mostly occur in Asia, Africa, and America [1]. Three major pandemics affecting geographic regions previously devoid of established foci are known to Western history, and Y. pestis spread to all continents except Australia and Antarctica [2,3]. The zoonotic plague disease can be transmitted from natural host reservoirs, mostly rodents, via various vectors to other mammals including humans. It is therefore a multi-host and multi-vector pathogen [4].
In the past ten years the use of modern molecular genetics to investigate isolates recovered from most natural foci as well as remains from victims of past pandemics has dramatically increased our understanding of the population structure, origin and spread of this major pathogen. The current view is that Y. pestis can be divided in biovars or ecotypes, grouped into subspecies pestis and subspecies microtus. Subspecies microtus comprises a number of biovars mostly harmless for humans. Microtus was initially investigated by microbiologists from the Former Soviet Union (FSU) under the name pestoides and subsequently under the different phenotype-based biovar designations Caucasica, Ulegeica, Altaica, Hissarica, and Talassica [4]. More recently two additional biovar designations were defined to cover Chinese microtus lineages, namely Xilingolensis and Qinghaiensis [5]. Whole genome sequencing and large scale SNP analysis has provided a robust branching order of the main clades within Y. pestis. The Caucasica biovar recovered so far only from nearby foci in Georgia, Armenia, Azerbaijan and Russia, corresponds to branch 0.PE2 in the SNP typing nomenclature proposed by Achtman and colleagues [4,6,7]. Together with 0.PE7, 0.PE2 branched out most ancestrally from the linear tree leading from Y. pseudotuberculosis to Y. pestis subspecies pestis biovar Orientalis. The strains defining the 0.PE7 clade were first identified as peculiar by Li Microtus and strains from the 0 and 1 branches so far investigated by MLVA25 and by SNP analysis are shown [5,6]. Three Ulegeica, two Hissarica and nine Altaica strains not investigated by SNP analysis are also included. For completion, Table 1 gives further information about assignment of biovar, genotype, and origin. Colors reflect MLVA clustering as suggested by Li et al. [5]. The SNP branch assignment of each strain as defined by Morelli et al. is indicated (column Morelli2010) together with the strain ID and biovar designation [6]. The results of CRISPR analysis according to Cui et al. are shown in column group [18]. Bootstrap support values are indicated. The figure shows the satisfying terminal branches clustering achieved by MLVA but the sometimes incorrect and usually low bootstrap values of deep branching nodes illustrating the complementarity of the two methods. doi:10.1371/journal.pone.0030624.g001 Figure 2. MLVA clustering and SNP branch assignment of 68 previously published Y. pestis pestis branches 1 and 2. Sixty-eight strains from the 1 and 2 branches previously investigated by both MLVA25 and SNP analysis are displayed [5,6]. For completion, Table 1 gives further information about assignment of biovar, genotype, and origin. Colors reflect MLVA clustering as suggested by Li et al. [5]. The SNP branch assignment of each strain as defined by Morelli et al. is indicated (column Morelli2010) together with the strain ID and biovar designation [6]. Bootstrap support values are indicated for each node. The results of CRISPR analysis according to Cui et al. are given in column group [18]. * This strain shows a Medievalis phenotype due to a different mutation in the napA gene compared to the mutation causing the Medievalis phenotype in the Medievalis biovar, as demonstrated by Pourcel et al. [13]. doi:10.1371/journal.pone.0030624.g002  Table 1 gives further information about assignment of biovar, genotype, and origin. Basic correlation and grouping of genotypes is similar compared to previously published Fig. 2 in Morelli et al. [6]. doi:10.1371/journal.pone.0030624.g003 et al. [5]. Only two strains corresponding to this clade have been reported so far (Figure 2 in [5]). C1961001 and C1962002 were recovered from Xinghai district in Qinghai province and importantly C1962002 was isolated from a human patient according to published information [5]. This would qualify clade 0.PE7 as a subspecies pestis biovar rather than microtus. The next branch 0.PE3 is represented by the unique Angola microtus strain, the geographic origin of which is uncertain [6]. It is followed by the branch leading to both 0.PE4 (microtus biovars Xilingolensis and Qinghaiensis) and 0.PE1 (represented by pestoides A, B, C, D with no correspondence provided in terms of microtus biovar designation [6,7]). The rest of the ancestral 0 branch is currently populated predominantly by strains originating from China focus B in the Xinjiang province which define three nodes 0.ANT1, 0.ANT2, and 0.ANT3 [6,7] potentially pathogenic for humans. The investigation of human remains associated with the Black Death demonstrated that the associated Y. pestis strains were almost coincident with the 3.ANT node [8][9][10] indicating that branches 1 (Orientalis biovar and Antiqua strains from Africa) and 2 (Antiqua strains from Tibet, Manchuria and Medievalis biovar) are less than 700 years old. The finding of many 0.ANT branches in China suggests that the Black death Y. pestis evolved in or near western China, and spread via a number of radiations to Southeast Asia, Africa, Europe, South and North America, leading to country-specific lineages [5,6].
Up to now, several hundred Y. pestis strains from the majority of known foci all over the world were analyzed and typed using MLVA based on VNTR loci selected from a collection of more than 60 loci shown to be polymorphic within Y. pestis [3,5,7,[11][12][13][14][15][16][17]. A significant fraction of these strains has also been typed by Clustered Regularly Interspaced Short Prokaryotic Repeats (CRISPR) [18] analysis and large-scale SNP typing [6]. Regarding Y. pestis, the comparability of those methods particularly MLVA and SNPs, and hence the applicability of progressive hierarchical resolving assays using nucleic acids (PHRANA) as earlier described for Bacillus anthracis has not been investigated so far [19,20].
Numerous Chinese and FSU isolates were amongst the investigated strains, but only four Mongolian isolates from two foci have Table 1. Overview of Y. pestis subspecies, biovar, genotype, and natural foci as suggested by different authors [4,5,6,7,18], and as deduced in this study.  [6] has not the same meaning as intermedium defined by Li et al. [5] which refers to Rhamnose positive Y. pestis pestis isolates. # prefix refers to foci as described by Anisimov et al. [4]. Numbers without # refer to Mongolian foci as shown in Figure 4. doi:10.1371/journal.pone.0030624.t001 been analyzed. Mongolia is a place of numerous highly active plague foci [4,21]. Y. pestis can be isolated in almost any province of Mongolia, and human plague is recorded since 1897 there, but was present for a much longer time in Siberian marmots [22]. In particular, Western Mongolia is an exceptional region in terms of Y. pestis diversity, as the Altaica and Ulegeica microtus biovars as well as Y. pestis subspecies pestis coexist in a relatively limited geographic space [4]. The present work was carried out to characterize 100 Mongolian Y. pestis strains from 37 different natural sampling places applying recent molecular analysis tools. CRISPR analysis and MLVA with 25 loci were used as a quick first line classification assay. Resulting data were compared to strains previously characterized by CRISPR, MLVA and SNP analysis [5,6,18]. MLVA clusters containing no strain previously typed by SNP analysis were assigned to SNP nodes by typing key SNPs on selected strains.

Direct comparison and aggregation of published MLVA and SNPs clustering
The work by Morelli et al. [6] was used to evaluate the relevance of CRISPR and MLVA cluster analysis carried out by Li et al., and Cui et al. [5,18], and to link the different clusters. One hundred and thirty-one strains (subsequently called linking strains) were investigated by both Li et al. and Morelli et al. [5,6]. MLVA clustering of this common set of strains using data from Li et al. [5] is shown in Figures 1 and 2. For each strain the SNP branch determined by Morelli et al. [6] is indicated. For instance, the 0.PE4a and 0.PE4b branches correspond to the Qinghaiensis microtus biovar, whereas the 0.PE4c and 0.PE4d branches correspond to the Xilingolensis microtus biovar (Figure 1 and 3, Table 1). The correspondence between 0.PE1 strains (pestoides A, B, C, D) and Altaica could be deduced by comparing MLVA data from Achtman et al. [7] and Li et al. [5], taking advantage of loci included in both assays ( Figure 1, Table 1). Figures 1 and 3 also include Ulegeica (from Mongolia) and Hissarica (from Uzbekistan) strains investigated only by CRISPR analysis and MLVA [5,18]. MLVA clustering suggests with moderate support that the Hissarica biovar is closest to the 0.PE1 and 0.PE4 microtus branches, but SNP typing will be required to confirm this assumption given the long MLVA branch leading to the Hissarica strains ( Figure 1). The interest of combining the MLVA discriminatory power, clustering efficiency and low cost with the phylogenetic robustness of SNPs illustrated here is in agreement with similar findings obtained for Bacillus anthracis [19,20]. Also, a  Table S1. Genotypes are explained in Table 1  MLVA clustering of the Mongolian isolates and tentative assignment of MLVA clusters to SNP branches One hundred Mongolian Y. pestis strains were analyzed by CRISPR and MLVA analysis and compared to previously published data of 366 [18], and more than 500 strains [5], respectively.
The 25 VNTR loci could be amplified in 96 of the 100 isolates. Sixty-five different MLVA25 genotypes are identified. Fifty-four of these are new compared to the current MLVA25 data [5]. The 100 isolates fall within six main clusters. Three clusters are Y. pestis subspecies microtus (11 isolates), the three others are Y. pestis subspecies pestis (89 isolates). The 11 microtus subspecies isolates fall into either Altaica (4 isolates from foci 7 and 8a), Xilingolensis (3 isolates from foci 23 and 33), or Ulegeica (4 isolates from foci 8, 10, 15) ( Figure 4). The remaining 89 strains belong to the biovar Antiqua, and are distributed to all known Mongolian foci ( Figure 4). Figures 5  and 6 show the resulting assignment for the 100 isolates together with previously investigated isolates from Li et al. [5].
Three of these clusters can be confidently assigned to a SNP branch owing to the co-clustering of at least one linking strain in the cluster: 0.PE1 (Altaica), 0.PE4 (Qinghaiensis/Xilingolensis), 2.ANT3 ( Figure 5 and 6, Table 2). These three clusters comprise 4, 3, 6 isolates corresponding to 3, 2, and 2 MLVA25 genotypes respectively. The four 0.PE1 (Altaica) isolates are very closely related to previously investigated strains from Mountain-Altai focus #36 in Anisimov et al. [4] (Figures 1, 4 and 5). One isolate from focus 8a shows the same MLVA25 genotype as two previously reported strains #2131 and I-3446 in Li et al. [5]. The three 0.PE4 (Qinghaiensis/Xilingolensis) isolates are closely related to previously published Xilingolensis (0.PE4c/d) strains isolated from the L focus in China (Figures 1, 4, Table 1). The MLVA25 genotype is identical or almost identical to one strain from the same focus called KP in Li et al. [5]. The H focus as defined by Zhou et al. [24], and Li et al. [25] is located in Manchuria, China, south of the L focus which hosts the Xilingolensis microtus biovar (Figure 4).
The three other clusters (87 strains) do not contain a linking strain which would allow a robust SNP branch assignment. One cluster corresponds to the Ulegeica microtus biovar (four isolates and MLVA25 genotypes) not included in Morelli et al. [6] (Figure 5/ green box). This assignment is deduced from the co-clustering with three Ulegeica strains investigated by Li et al. [5]. All three previously investigated Ulegeica strains originate from Mongolia. Strains of the two other MLVA-clusters correspond to Y. pestis subspecies pestis isolates ( Figure 5 and 6/brown and turquoise boxes). The smaller cluster comprises five Y. pestis pestis Antiqua isolates (4 MLVA25 genotypes) from foci 6, 8a, 30, 32. They are most closely related to isolates from Tuva focus #37 in Russia ( Figure 4 and 5/brown) immediately adjacent to the western border of Mongolia [4,5]. The third unassigned cluster is by far the most numerous and frequent in Mongolia (78 isolates, 41 MLVA25 genotypes; Figure 6/turquoise). It is closest to a small group of seven isolates shown in supplementary Figure 2 in Li et al. [5]. These seven isolates (5 MLVA25 genotypes), C1976008, C1976001, C1989002, C1961006, C1972002, C1972001, C1972003 (Figure 6), were collected in Akesai, Gansu province and Wulan, Qinghai province, China [5]. MLVA25 clustering tends to link these two clusters to 3.ANTa or 0.ANT branches. The relative MLVA genotype diversity of the two groups is shown in Figure 7. The diversity observed in Mongolia is much larger than the diversity in China, but this may be due to the larger number of available Mongolian strains. The two groups are clearly resolved by MLVA suggesting low level of strain circulation and cross contamination between the Mongolian and Chinese foci.

SNP typing of selected strains
The assignment of three clusters to the 0.PE1, 0.PE4 and 2.ANT3 SNP-defined branches as suggested by MLVA25 clustering with linking strains could be confirmed for selected strains by typing relevant SNPs ( Table 2). The remaining three MLVA clusters were positioned on the SNP tree by typing a few selected strains from each cluster for key SNPs according to Morelli et al. [6] (Table 2).
Ulegeica strain MNG 2972 representative of the Ulegeica cluster could be assigned to branch III-VI by analyzing all 56 SNPs in this branch [6]. Eleven SNPs, s85, s90, s463, s846, s849, s940, s951, s1099, s1221, s1248, and s1351 showed the derived genotype, the other 45 showed the ancestral state. This enables the precise positioning of the Ulegeica branching node in between III-VI and indicates that Ulegeica is the most recent microtus branch characterized so far. We propose to call the Ulegeica clade 0.PE8 in agreement with the published SNP branch nomenclature (the Hissarica biovar could be assigned as 0.PE9, if additional SNP analyses will confirm this assumption) ( Table 1).
Eight selected strains from the remaining two MLVA25 clusters, MNG 649, MNG 3054, MNG 2986, MNG 2853, MNG 2645, MNG 3088, MNG 3143, MNG 2881, revealed a derived genotype for the tested SNPs for all branches connecting nodes 0.PE3a and 3.ANTa, and an ancestral genotype for the tested SNP for the branches XI-3.ANTa (as well as XIII-XI, 1.IN2a-XIII, XII-XI) and VIII-3.ANTa (as well as VIII-2.ANT3a and 2.ANT3a-2.ANT2a). This demonstrates that the two clusters are branching out within the 0.ANT3a, VIII and XI nodes ( Table 2). Further SNP typing and whole genome draft sequencing of a few selected strains will allow to determine the exact positioning of the two Mongolian clades.
Progressive hierarchical resolving assays using nucleic acids (PHRANA) In this study, MLVA25 analysis allows a classification of the 100 Mongolian strains into 65 genotypes defining six clusters, corresponding or closely related to known biovars, or recently described lineages (Table 1, 2). Partial SNP typing was then applied to selected strains in order to anchor more precisely clusters devoid of linking strains. This approach is a PHRANA approach in which MLVA, rather than SNP typing as initially proposed, is used as a first line assay. The CRISPR typing yields 14 different genotypes. The CRISPRs alleles in Y. pestis pestis are exceptional by the fact that most spacers (except for the oldest a1-a6) originate from the chromosome (or the plasmids) as compared for instance to Y. pseudotuberculosis spacers [26]. The new spacers identified in this work also originate from the chromosome. Very interestingly, Ulegeica is unique among the microtus lineages in that its recent spacers a82-a85-a86 originate from the chromosome. In this respect, Ulegeica is close to Y. pestis pestis.

Conclusions
The present investigation illustrates and confirms the large variety of the Mongolian microtus biovars Ulegeica, Altaica, and Xilingolensis present in a close proximity ( Figure 5, 6, and 7). It suggests that western Mongolian foci or the adjacent Siberian foci are likely places of emergence of Ulegeica, the most recent clade. Xilingolensis would have spread throughout Mongolia, to focus L in Manchuria (Figure 4). Qinghaiensis is found further south in central China (focus M) [5]. The Mongolian microtus Ulegeica clade is shown to be the most recent microtus branch along the linear tree leading from Y. pseudotuberculosis to Y. pestis subspecies pestis biovars Intermedium, Antiqua, Orientalis, Medievalis. Ulegeica contributes to the filling of a large gap.
The absence in Mongolia of Y. pestis pestis lineages branching along the III-0.ANT3a segment is consequently surprising and might indicate that the presence of Y. pestis pestis in Mongolia is the result of a secondary introduction of strains from China [6] perhaps in the last hundred years as human infections are reported since 1897 in this country. The Gansu province south of Mongolia, in which closest neighbors from the most frequent Mongolian Y. pestis pestis are present, is a likely source. Alternatively, Mongolian 0.ANT representatives might have been replaced by the more recently emerged 3.ANT lineage, and become extinct or at least very rare in Mongolia. More detailed whole genome sequencing and SNP analysis will be required in order to test these two hypotheses and precisely deduce the direction of the dissemination of the 3.ANT lineage across Mongolia and China. In addition, a systematic  MLVA typing of Y. pestis strain collections may enable the identification of other rare clades as previously illustrated [5,6].

Strains and DNA
In this study only strains isolated from wildlife animals or their parasites were investigated to focus on genotypes occurring in nature and with a clear geographic assignment, in contrast to strains recovered from patients who may have travelled recently (Table S1). The investigated plague-strains were conserved in glycerol stocks. For this work, they were recovered on Hottinger's agar at 28uC for 24 h and subcultured on Columbia blood agar at 28uC for 24 h. Thermolysates were prepared by heating a bacterial suspension for 30 min at 95uC. The 100 Y. pestis strains investigated in this study were collected between 1960 and 2007 from 37 sampling sites in Mongolia (Figure 4) distributed over 13 aimags (provinces). They were isolated from various parasitic plague-vectors, such as Oropsylla silantiewi, but also from lice or ticks (Table S1). Parasites were collected from mammalian host species, such as Marmota sibirica. All strains revealed both Y. pestis specific virulence plasmids pMT1, and pPCP1 when investigated by previously published real-time PCR [27].

MLVA markers and PCR amplification
Twenty-five VNTR markers were applied [5]. Three loci were co-amplified in a single multiplex PCR and the resulting products were analysed on a CEQ8000 capillary electrophoresis machine (Beckman-Coulter, Marseille, France) essentially as described [14,28]. The resulting data were analyzed and merged with the previously generated MLVA database including more than 500 Y. pestis strains the majority of which from Central Asia using BioNumerics software package version 6.5 (Applied-Maths, Sint-Martens-Latem, Belgium) [5]. The tree was rooted using two Y. pseudotuberculosis isolates as an outgroup. MLVA data corresponding to pestoides strains A, B, C and D in Achtman et al. [7] were kindly provided by Dr. Paul Keim.

SNP typing
At least one SNP was selected for each relevant branch to determine branching of the Mongolian plague strains within the previously published SNP minimum spanning tree [6]. Each SNP was amplified by conventional PCR, sequenced and analyzed by alignment with Y. pestis type strain CO92 (Table 2).

Ethics Statement
The bacterial strains in this study were obtained from nonvertebrate vectors, collected from various mammal-species (Table  S1). Mammals were trapped in one-door live traps, as previously described [34].