Phylogenetic Analysis of Entomoparasitic Nematodes, Potential Control Agents of Flea Populations in Natural Foci of Plague

Entomoparasitic nematodes are natural control agents for many insect pests, including fleas that transmit Yersinia pestis, a causative agent of plague, in the natural foci of this extremely dangerous zoonosis. We examined the flea samples from the Volga-Ural natural focus of plague for their infestation with nematodes. Among the six flea species feeding on different rodent hosts (Citellus pygmaeus, Microtus socialis, and Allactaga major), the rate of infestation varied from 0 to 21%. The propagation rate of parasitic nematodes in the haemocoel of infected fleas was very high; in some cases, we observed up to 1,000 juveniles per flea specimen. Our study of morphology, life cycle, and rDNA sequences of these parasites revealed that they belong to three distinct species differing in the host specificity. On SSU and LSU rRNA phylogenies, these species representing three genera (Rubzovinema, Psyllotylenchus, and Spilotylenchus), constitute a monophyletic group close to Allantonema and Parasitylenchus, the type genera of the families Allantonematidae and Parasitylenchidae (Nematoda: Tylenchida). We discuss the SSU-ITS1-5.8S-LSU rDNA phylogeny of the Tylenchida with a special emphasis on the suborder Hexatylina.


Introduction
More than 150 species of fleas feeding on different mammalian hosts, primarily rodents, are vectors of the bacterium Yersinia pestis, a causative agent of plague [1,2]. In natural foci of plague, the dynamics of flea populations are among the main factors controlling the incidence of epizootics that pose a threat to humans inhabiting the areas [3][4][5]. Entomoparasitic nematodes of the order Tylenchida are known to control populations of various insect hosts [6][7][8][9]. The rate of tylenchid infestation in fleas reaches 50-60% in some cases [10,11], when the nematodes cause castration and early death of the flea hosts [9,12,13].
Despite high importance of the Tylenchida as a nematode order harboring entomoparasites and notorious crop pests, their reliable phylogeny is still a challenge. Tylenchid nematodes differ widely in life cycle, parasitic strategies, and the host range that spans plants, fungi, and invertebrates.
Phylogenies obtained from SSU and partial LSU rDNA data often disagree with classifications based on morphology and life cycle [14][15][16][17][18][19][20][21]. Phylogenetic resolution inside the order is far from being clear, which in many respects results from the insufficiency of data available to adequately describe its diversity. As for tylenchid parasites of fleas, only 31 species are described to date [9,[22][23][24][25][26][27][28][29][30][31], with no molecular vouchering. Here we present a study of parasitic nematodes isolated from fleas sampled from different rodent hosts in a natural focus of plague.  major). Three flea species (Citellophilus tesquorum, Neopsylla setosa, and Frontopsylla semura) were sampled on sousliks; two species (Amphipsylla rossica and Ctenophthalmus secundus) were on M. socialis voles; and one species (Mesopsylla hebes) was on jerboas. Fleas were examined for nematode infestation (Table 1). Examination and dissection of fleas were carried out using the dissecting microscope MBS-2 (LOMO, Russia). A half of parasitic nematodes sampled from each flea was preserved for subsequent DNA extraction, and another half was used for morphological analysis. Live fleas infected with nematodes were placed in glass flasks with river sand to obtain free-living forms. Insects were kept in a KBF 720 (E5.2) climate chamber (Binder, Germany) at 26 ∘ C and 80% humidity.

Morphological Analysis.
Fixation and clarification of nematode preparations were performed using standard techniques described by De Grisse [32]. Material was mounted on slides in a drop of glycerin, bound by a paraffin circlet (http://pest.cabweb.org). Color staining of preparations was not performed. Morphometric analysis was conducted using the light microscope "Leica DM 1000" (Leica, Germany) with an eyepiece micrometer. Pictures of nematodes were taken with the microscope "DFC 425" (Leica, Germany). Published data on morphometrics [23,25,26] were used for comparison.

DNA Extraction, PCR, and
Sequencing. DNA samples were extracted with a Diatom DNA Prep (IsoGen Lab, Russia). rDNA fragments were amplified using an Encyclo PCR kit (Evrogen, Russia) and primers given in Table 2. The amplified rDNA fragments were sequenced using an Applied Biosystems 3500xL DNA analyzer. Sequence reads were assembled with the CAP contig assembly program [33] and proofread with the BioEdit software [34]. For three isolates, almost complete sequences of 18S and 28S rRNA and complete sequences of 5.8 rRNA, internal transcribed spacers ITS1 and ITS2 were assembled. The sequences were submitted to GenBank under accession nos. KF155281-KF155283. For the rest of isolates, partial (750-800 bp) sequences of 18S and 28S rRNA genes were submitted to GenBank under accession nos. KF373731-KF373740.

Phylogenetic Analysis.
The newly obtained rDNA sequences of tylenchid parasites of fleas were aligned with a selected set of other tylenchid sequences obtained from the GenBank. The main selection criterion was to sample representatives of all clades that occur in published SSU and LSU rDNA phylogenies of the Tylenchida [16][17][18][19][20][21]39]. Apart from the D2-D3 LSU rDNA expansion segment commonly used in previous studies, we included all LSU rDNA sequence data available for the Tylenchida, with the exception of Basiria sp. SAN-2005 (accession nos. DQ145619, DQ145667) that in our preliminary analyses (data not shown) demonstrated a disputable affinity to the Tylenchida. For the species Anguina tritici, Globodera pallida, Heterodera glycines, Pratylenchus vulnus, and Radopholus similes the nearly complete rDNA sequences were assembled with appropriate cDNA fragments identified with BLAST [40]. Partial LSU rDNA sequence of Ditylenchus dipsaci was combined with the soil environmental clone NTS 28S 061A 2 b4 (accession no. KC558346), as the clone sequence appeared to represent a close tylenchid relative of D. dipsaci. Chimeric sequences were also created in some cases when closely related partial rDNA sequences were found in the database. All sequences and their accession numbers are listed in Table 3. Cephalobidae and Chambersiellidae were chosen as the outgroup. Alignments were constructed with the MUSCLE program [41] and refined manually using the MEGA 5.0 software package [42]. Three alignments were generated: (1) SSU rDNA, (2) D3 region of LSU rDNA, and (3) concatenated rDNA data including SSU, LSU, 5.8S rDNA, and highly conserved regions of ITS1. After discarding ambiguously aligned positions, the alignments length was 1,723, 592, and 4,930 positions, respectively. Bayesian reconstruction of phylogeny was done with the PhyloBayes software, version 3.2 [43] under the GTR + CAT + DP model [44]. Eight independent runs were performed with 4,000,000 cycles each; the first 3,000,000 cycles were discarded. A consensus tree with Bayesian posterior probabilities was constructed for the remained tree sample. Bayesian reconstruction was also performed using the MrBayes software [45] under the GTR + G8 + I model [46] in two independent runs, each with four Markov chains. The chains were run for 5,000,000 generations, with trees sampling every 1,000th generation. The consensus posterior probabilities were calculated after discarding the first 3,000,000 generations. Partitioning "by genes" was used for the concatenated alignment with all parameters unlinked, except for the topology and branch lengths. In addition, node support was estimated with maximum likelihood bootstrap as implemented in the RAxML software, version 7.2.6 [47], under the GTR + G + I model with 1,000 bootstrap replicates. Alternative topologies were tested using the approximately unbiased (AU) [48] and Kishino and Hasegawa [49] tests implemented in the CONSEL software [50] and the expected likelihood weight test [51] implemented in the TREE-PUZZLE software [52]. TREEVIEW  Adult parasitic females and their progeny were found in the haemocoel of infected fleas. In the infected fleas C. tesquorum, A. rossica, C. secundus, and Mesopsylla hebes, only one generation of parasitic females was observed. Their amount in a flea specimen is determined by the number of free-living infective females that penetrate into the flea larva. We observed 1 to 2 or 1 to 4 adult parasitic females per flea specimen in spring and autumn, respectively. An additional parthenogenetic generation of parasitic females was found in some fleas of N. setosa and F. semura, where        Opperman, unpublished [96].

Heteroderidae
Heterodera glycines     up to 16 specimens per flea were observed. As in other entomoparasitic nematodes, the propagation rate depends on the host age. Thus, in young fleas up to 10 juveniles was found per flea specimen, whereas up to 1,000 juveniles of different stages were contained in some old fleas ( Figure 2). After the 2nd molt the number of juveniles is maximal, and 3rd stage juveniles massively migrate to the rectal section of the flea intestine for exit to the environment. In some cases, the observed infestation level was so high that nematodes penetrated distal segments of the flea legs, from where they have no way to the environment.

Morphological Analysis of Entomoparasitic Stages in Nematode Isolates and Their Taxonomic Identification.
Analysis of morphology of entomoparasitic stages suggests that the studied nematode isolates from three distinct groups. A single generation of parasitic females was observed in the first two groups and an additional parthenogenetic generation-in the third group. According to morphometric data on adult parasitic females (Tables 4-6), the first two groups belong to the genera Rubzovinema or Spilotylenchus and the third group to the genus Psyllotylenchus. Photographs of parasitic females of Rubzovinema sp., Spilotylenchus sp., and Psyllotylenchus sp. are depicted in Figure 3. Figure 4 shows their distribution among flea samples studied. According to morphometric evidence, parasitic females and juveniles of the genera Rubzovinema and Spilotylenchus are very similar. However, in the first two groups of isolates we found characters bearing discriminative and identificational value. In particular, the oesophageal glands in juveniles III of the first group are poorly developed. This is a distinctive feature of the genus Rubzovinema, where males and females have shortened oesophageal glands located close to the nerve ring. In the second group of isolates, oesophageal glands are well developed and elongated, which is characteristic of the genus Spilotylenchus. In the first group, the stylet possesses a heavily sclerotized distal spear with a length of approximately half the total stylet length and has a stem with a weaker sclerotization and widening to the base. This stylet structure is characteristic of the genus Rubzovinema, and stylet length (18.5 (14-22) m) is in accordance with morphometrics given in the description of this genus [26]. In the genus Spilotylenchus, the stylet varies in shape but always possesses a shortened conical distal spear. In the second group of isolates, the stylet structure was similar to that of Spilotylenchus. Also, the vulval lips of the first group are more protruded than in Spilotylenchus. Other features, including the morphometrics, vary widely in both genera, which hampers taxonomic identification. Nevertheless, based on distinctive traits, we identified the first and second group of isolates as Rubzovinema sp. and Spilotylenchus sp., respectively.
In the genus Rubzovinema, the single species described to date is Rubzovinema ceratophylla [26]. This species is known to parasitize exclusively the flea Citellophilus tesquorum that feeds on sousliks. The specimens of Rubzovinema studied in this work were isolated from five flea species, C. tesquorum, Neopsylla setosa, Frontopsylla semura, Amphipsylla rossica, and Ctenophthalmus secundus, of which the latter two were sampled on mouse-like rodents. Also, the parasitic females of Rubzovinema sp. differed from R. ceratophylla by morphology; they have a shorter tail and more protruded vulval lips. A morphometric comparison of Rubzovinema sp. and R. ceratophylla is given in Table 4.
The parasitic females of Spilotylenchus sp. were isolated from the flea Mesopsylla hebes associated with jerboas. The females were not identified to the species level because of a small number of available specimens and the lack of a freeliving stage. A morphometric comparison of Spilotylenchus sp. and the morphologically closest species Spilotylenchus maisonabei [23] is given in Table 5.
In the genus Psyllotylenchus, descriptions of most species are fragmentary and incomplete, which precluded the species identification of the Psyllotylenchus isolates from the fleas N. setosa and F. semura feeding on sousliks. A morphometric comparison of Psyllotylenchus sp. and the type species of this genus, Psyllotylenchus viviparous [25], is given in Table 6.
The 18S and 28S rDNA sequences of Rubzovinema sp. specimens from A. rossica and C. secundus were 100% identical, which indicates that the isolates belong to the same species. The sequences of Rubzovinema sp. ex C. tesquorum, Rubzovinema sp. ex N. setosa, and Rubzovinema sp. ex F. semura diverged from one another and from the gene sequences of Rubzovinema sp. ex A. rossica and Rubzovinema sp. ex C. secundus by 0.4-0.7%, which corresponds to the levels of intraspecific variation [14,[114][115][116][117][118][119]. The 18S and 28S rDNA sequences of Psyllotylenchus sp. ex N. setosa and Psyllotylenchus sp. ex F. semura were 100% identical, indicating that they belong to the same species. The 18S and 28S rDNA sequences of Rubzovinema sp. and Psyllotylenchus sp. diverge by 1.2% and 1.9%, respectively. Those of Spilotylenchus sp. ex M. hebes were found to be more divergent. The degree of divergence of the 18S rDNA sequence of Spilotylenchus sp. ex M. hebes from those of either Rubzovinema sp. or Psyllotylenchus sp. was 2.4%; the D3 expansion segment of 28S rDNA diverged by 13.1% and 12.0%, respectively. The observed divergence rate of rDNA sequences agrees well with published evidence on entomoparasitic nematodes [14,[114][115][116][117][118]. Thus, intraspecific divergence of 18S rDNA in Deladenus siricidicola is 1% [120], of D2 and D3 expansion segments in the phytoparasite Bursaphelenchus xylophilus is from 0% to 0.6%, and the interspecific variation between the     phytoparasites B. xylophilus and Bursaphelenchus mucronatus is from 1.7% to 3.7%. The spacers ITS1 and ITS2 are generally more diverged; the intra-and interspecific variation for these species is from 0 to 3.1% and 11.2 to 13.4%, respectively [121][122][123]. Molecular vouchering is proved to efficiently complement morphological species identification in nematodes [73,122,[124][125][126][127][128]. Combining the rDNA and morphological data confirms the species identity within each of the three studied groups of isolates.

Phylogenetic Analysis.
In phylogenetic analyses of rDNA we used a dataset with extensive species and gene sampling (SSU-ITS1-5.8S-LSU) compared to earlier published tylenchid phylogenies, most of which were based on SSU rDNA or D2-D3 expansion segments [17, 19-21, 39, 129]. The SSU-ITS1-5.8S-LSU rDNA tree topology ( Figure 5) is highly similar to other published phylogenies of tylenchids. In this tree, tylenchomorphs are represented by the sister groups Aphelenchidae and Tylenchida. Most of the tylenchid clades occur in published trees but often contradict classifications based on morphology, as it was also noted by other authors [17, 19-21, 39, 129] The nonredundant rDNA data on the Hexatylina in Gen-Bank mostly represents the D2-D3 expansion segments of LSU rDNA. To maximize species sampling of the Hexatylina, we chose the D3 expansion segment as the molecular marker. The phylogenetic tree with the Anguinoidea as an outgroup is shown in Figure 6. In this tree, the suborder Hexatylina consists of two well-supported clades, in accordance with previously published D2-D3 rDNA phylogenies [19,20,39] (Table 7). All three tests do not discriminate between the alternative placement of the flea parasites as closest to the Allantonema, Parasitylenchus, or Deladenus branches; however, its positioning outside this grouping is not rejected only by a less conservative Shimodaira-Hasegawa test [50].
The second major clade of the Tylenchida includes representatives of the classic infraorder Anguinata, with a wellsupported monophyletic origin, except for a few species. They belong outside the second clade and may initially have been wrongly identified.
The genus Howardula is paraphyletic in published rDNA and mitochondrial COI phylogenies [71]. Such characters of Howardula as the degeneration of oesophagus, tail shape, and the absence of stylet in males seem to have evolved independently by convergence. The paraphyletic genus Deladenus is more closely related to either ancestral forms of the Hexatylina or forms typical to the Anguinata. The infraorder Anguinata includes soil-dwelling nematodes, mostly mycetophagous or parasitizing various parts of plants. However, an unidentified entomoparasitic nematode was also grouped within the Anguinoidea [39]. The life cycle of Deladenus spp. is an irregular alternation of free-living and entomoparasitic forms. The nematode D. siricidicola is able of producing an unlimited number of free-living generations in the absence of the host larvae of siricid 0.1 Primarily entomoparasitic, most with free-living mycetophagous or plant-parasitic generation Obligate parasites of plant roots Feeding on algae, mosses, and fungi; parasites of plants Mostly mycetophagous Bacterial feeding Figure 5: Phylogenetic tree of Tylenchida, inferred from SSU-ITS1-5.8S-LSU rDNA sequences. Topology was inferred using the PhyloBayes software (maxdiff = 0.36). Node support values are shown as follows: the first two values are Bayesian posterior probability assessed using the PhyloBayes and MrBayes software, respectively, and the third is bootstrap support assessed by the ML method. Thick lines lead to the nodes, in which at least one support value of posterior probability is 0.95 and higher. Names of clades (framed) are mainly given by type genera included in them (with the exception of Iotonchioidea). Formal taxonomic position (family by [8]) is shown on the right to the color bar. Colors indicate the ecologies (see the legend). Names of the species of Hexatylina that have a mycetophagous stage in their life cycle are shown in blue. The three robust major branches of Tylenchida are marked by gradient.    Figure 4) and on (2) the relationships within the monophyletic branch that includes the studied group of nematodes parasitizing fleas (designated by asterisk).

Topology
Rank pine-killing wood wasps [132]. Like in Anguinata, the freeliving forms of Deladenus spp. are fungal feeding. Such characters of Deladenus asthe mycetophagy, enlargement of subventral glands in entomoparasitic females versus their reduction in free-living forms, the hypertrophy of dorsal glands, and stylet reduction in free-living forms seem to be symplesiomorphic. Resemblance with the Anguinata is also typical of other mycetophagous free-living forms: Hexatylus (Neotylenchidae), Rubzovinema (Neotylenchidae), Prothallonema (Sphaerularioidae) Helionema (Hexatylina dubia), and Paurodontidae. For the latter, the entomoparasitic stage is expected but has never been observed. The relationship between the Hexatylina and Anguinata was earlier hypothesized based on morphology [7,8,130,133,134]. On rDNA phylogenies of tylenchids, the monophyly of the Hexatylina + Anguinata is either supported [19] or not rejected [20]. In the SSU-ITS1-5.8S-LSUrDNA tree obtained in this study, the monophyly of the Hexatylina + Anguinata has the Bayesian posterior probability of 0.91, but the maximum-likelihood bootstrap support is low; the AU and Kishino and Hasegawa tests did not discriminate between alternative hypotheses. According to our SSU-ITS1-5.8S-LSU rDNA phylogeny ( Figure 5), the major robust branches of the Tylenchida are incongruent with morphology-based classifications suggesting three rather than four suborders (the rank is adopted from morphological systems of tylenchids). Among them, the Hexatylina and Anguinata (both are monophyletic) are likely to be sister groups. The third emerged suborder includes representatives of three classic suborders: Tylenchina, Hoplolaimina, and Criconematina, among which only the latter does not contradict morphology-based classifications.
Considering ecological traits coded in Figure 5, the mycetophagy and/or facultative ectophytoparasitism are likely to be ancestral in the Tylenchida. Sedentary phytoparasites (root-knot species of Meloidogyne, the false root-knot genus Nacobbus, and cyst-forming Heterodera and Globodera) and other obligate endoparasites of plants evolved several times from free-living or facultative sedentary forms, as it was previously hypothesized in accordance with the concept of evolutionary trend to endoparasitism in phytonematodes [135]. Similarly, obligate endoparasites of insects from the Hexatylina are likely to have evolved from mycetophagous forms, with some species retaining the ancestral mycetophageous stage in the life cycle (e.g., species of the paraphyletic genus Deladenus and flea nematodes of the genus Rubzovinema). An interesting specific case in the Hexatylina is the genus Fergusobia that includes plant parasites associated with insects [68,70], which may have transited to plant parasitism via entomoparasitism [39].

Ribosomal DNA Phylogeny of the Flea Nematodes and
Their Classification. The nematodes of fleas do not group with the families known as their relatives in morphologybased systems, as these families do not form monophyletic groups in the tree. However, they do group with both type genera of the families Parasitylenchidae and Allantonematidae (Parasitylenchus and Allantonema, resp.). This grouping is preceded by a successive divergence of Deladenus durus and Deladenus siricidicola ( Figure 5). As mentioned above, the pronounced free-living form in Deladenus seems to be ancestral to this group.
Only 31 tylenchid species that parasitize in fleas have been described to date. They differ by morphology, life cycle, and the host specificity, and belong to the five genera: Spilotylenchus (8 species), Psyllotylenchus (20 species), Incurvinema (1 species) Kurochkinitylenchus (1 species), and Rubzovinema (1 species). According to the classification of Siddiqi [8], the genera Spilotylenchus and Psyllotylenchus belong to the family Parasitylenchidae, whereas the genus Rubzovinema is a member of the Neotylenchidae. The two families represent two superfamilies, Iotonchioidea and Sphaerularioidea, respectively. All rDNA phylogenies published to date suggest that these superfamilies are paraphyletic [19,20,39], which is also inferred in our study with an extensive gene and taxon sampling.
A high degree of rDNA similarity in the three studied species suggests a closer relationship of these species than that assumed by the accepted system of classification. Earlier, Slobodyanyuk proposed to unite all known flea parasites into one family, the Spilotylenchidae. Its four subfamilies, Spilotylenchinae, Rubzovinematinae, Psyllotylenchinae, and Kurochkinitylenchinae, are discriminated based on the life cycle features [28]. In Spilotylenchinae and Rubzovinematinae, the entomoparasitic stage is represented by parasitic females of one heterosexual generation. In Psyllotylenchinae, in addition to the heterosexual generation, a parthenogenetic generation occurs in the flea haemocoel. In Kurochkinitylenchinae, two heterosexual generations exist in the haemocoel: the first generation produces parasitic females and the second generation produces both females and males [28]. Siddiqi also considered the unification of all flea tylenchids into one family but observed the need for further evidence in support [8].
Our results strongly suggest the inclusion of the three genera, Rubzovinema, Psyllotylenchus, and Spilotylenchus, in one family, the Spilotylenchidae [28]. The ribosomal DNA genetic distance within the family Spilotylenchidae is much smaller than that of certain tylenchid genera, for example, Meloidogyne (Figure 4) or Pratylenchus [19,84].

Host Specificity of Flea
Nematodes. The majority of tylenchid nematodes are monoxenous or oligoxenous; in particular, flea parasites were thought to be strictly host specific. Earlier papers suggested the lack of strict host specificity in Psyllotylenchus pawlowskyi and Psyllotylenchus viviparous [13,25]. However, later these species were found to be heterogeneous and sustained revision [9,[27][28][29]. Spilotylenchus pawlowskyi and Spilotylenchus caspius were referred to as single-host parasites of the flea Coptopsylla lamellifer [27,136]. Kurochkinitylenchus laevicepsi and Spilotylenchus ivashkini also share the same flea host, Nosopsyllus laeviceps [28,29]. Before our study, the genus Rubzovinema was known to contain a single species, Rubzovinema ceratophylla, which parasitizes exclusively the flea Citellophilus tesquorum.
We found that at least two out of the three studied species are not single-host parasites. Psyllotylenchus sp. was shown to parasitize two flea species feeding on sousliks, Frontopsylla semura and Neopsylla setosa. Rubzovinema sp. was found on five flea species feeding on different rodent hosts: C. tesquorum, F. semura, N. setosa (all sampled from sousliks), Ctenophtalamus secundus, and Amphipsylla rossica (all sampled from voles). A. rossica, F. semura, and C. tesquorum belong to different families of the superfamily Ceratophylloidea (Leptopsyllidae and Ceratophyllidae), whereas C. secundus and N. setosa belong to the superfamily Hystrichopsylloidea. Unlike the host-specific R. ceratophylla, the studied Rubzovinema sp. parasitizes taxonomically distant fleas feeding on different rodents. Thus, the common opinion that flea nematodes are strictly host specific should be revisited.
As the two species of Rubzovinema demonstrate, even closely related parasites may exhibit different host range size. Among other known examples are the entomoparasitic nematodes of the genus Howardula parasitizing various beetles and flies [71,137,138], many phytonematodes [8], sibling species of parasitoid flies [128], and herbivorous insects [139]. The host range of parasites is an indicator of their evolutionary strategy in the ecosystem. Multihost parasites can be considered ecological generalists, in contrast to specialists that coevolve with a particular host. Generalists and specialists play different roles in the ecosystem [140], where they keep in balance, taking advantages and disadvantages of the two strategies. The advantages of generalization are yet to be explained by evolutionary biologists, whereas advantages of specialization are obvious, and it is generally accepted that evolution favors specialism [141,142]. In the flea parasites, this trend is demonstrated by a greater species diversity of ecological specialists, the genera Spilotylenchus and Psyllotylenchus.
Nevertheless, the generalist Rubzovinema sp. was most abundant in the studied samples, which indicates that extending the host range may be evolutionarily successful. Besides the need to combat the immune response of several hosts, which is a requirement to widen the hosts range [143], the free-living stage of Rubzovinema sp. is to adapt to diverse microbioclimatic conditions of complex environments of rodent habitats. Multihost parasites pay a cost of adapting to alternative conditions [141,144] compensated by stable survival of the species. Considering the spatial and temporal dynamics of flea populations feeding on a particular rodent host (one or two flea species usually dominate over a sampling season), multihost nematode parasites gain an advantage of their relative independence of population waves of either flea hosts or their rodent hosts. A higher infestation rate observed for Rubzovinema sp., compared to the two other studied species, may be an indicator of a greater ecological plasticity of this multihost parasite.

Entomoparasitic Nematodes in Natural Foci of Plague.
In natural foci of plague, the epizootic dynamics are influenced by numerous climatic and biotic factors. The spatial and temporal population dynamics of the plague agent, Y. pestis, affect the population dynamics of the flea vectors and their mammalian hosts. Members of the transmission route of the plague agent also closely interact with other living organisms. For example, parasites of fleas that in turn feed on rodents are hyperparasites that play the role of high-level control agents on the ecosystem level, the role that entomoparasitic nematodes share with the bacterial plague agent. Highlevel control agents render the epidemiological state of a natural focus of disease less predictable. On the one hand, a lower density of the flea vector population reduces the plaque transmission rate; on the other, its growth causes an exponential decay of the host rodent population [145] below its epidemiological threshold, above which there is a threat of spillover of plague infection into human population [145]. Hypothetically, nematode-induced decrease of flea population is able to increase the number of rodents above the threshold and thus trigger an epidemic. The dual effect of high-level control agents is well exemplified by cases, when during plague episodes the extermination of rodents by humans causes the return of infection through stimulating the migration of fleas, the plaque vectors [5].
The studied entomoparasitic nematodes possess high potential as control agents of the flea vectors of plague owing to their high propagation rate within the flea host ( Figure 2) and high infestation level (up to 21% observed in this study and from 50 to 60%, as estimated by other authors [10,11]). One of the studied nematode species, Rubzovinema sp., is a multihost parasite. Host-specific parasites reach the optimal level of pathogenicity by maintaining the tradeoff between pathogenicity and transmissibility. Adding of a new host to a multihost system makes the model more complicated [141]. The multihost parasite Rubzovinema sp. is expected to exhibit different levels of pathogenicity with respect to different flea hosts which, in turn, play different roles in the transmission of plague. Epizootics cause sporadic mortality in local populations of all members involved in the interaction with the plague agent, and their survival is contingent on migrations within a metapopulation. It is the case when the Cope's law [139,146] governs the extinction of specialists on a shorter time scale rather than a geological period, and evolution may favor the ecological generalists, such as Rubzovinema sp.
Some authors surmised the involvement of entomoparasitic nematodes in the transmission of the plague agent [4], as it was observed that biofilms of Yersinia pestis adhere to cuticle receptors of Caenorhabditis elegans [147][148][149]. In this perspective, nematodes parasitizing fleas in natural foci of plague take on greater importance, as they may provide for the transmission route that does not include a mammal [4]. Further studies will clarify the role of flea nematodes in the transmission of plague infection.