Comparative phylogeography of the plateau zokor (Eospalax baileyi) and its host-associated flea (Neopsylla paranoma) in the Qinghai-Tibet Plateau

Background Specific host-parasite systems often embody a particular co-distribution phenomenon, in which the parasite’s phylogeographic pattern is dependent on its host. In practice, however, both congruent and incongruent phylogeographic patterns between the host and the parasite have been reported. Here, we compared the population genetics of the plateau zokor (Eospalax baileyi), a subterranean rodent, and its host-associated flea species, Neopsylla paranoma, with an aim to determine whether the two animals share a similar phylogeographic pattern. Results We sampled 130 host-parasite pairs from 17 localities in the Qinghai-Tibet Plateau (QTP), China, and sequenced a mitochondrial DNA (mtDNA) segment (~2,500 bp), including the complete COI and COII genes. We also detected 55 zokor and 75 flea haplotypes. AMOVA showed that the percentage of variation among the populations of zokors constituted 97.10%, while the within population variation was only 2.90%; for fleas, the values were 85.68% and 14.32%, respectively. Moreover, the flea Fst (fixation index) values were significantly smaller than in zokor. Although the Fst values between zokors and fleas were significantly and positively correlated (N =105, R =0.439, p =0.000), only a small amount (R2= 0.19) of the flea Fst variations could be explained by the zokor Fst variations. The two animals showed very distinct haplotype network structures from each other while co-phylogenetic analyses were unable to reject the hypothesis of an independence of speciation events. Conclusions Zokors and fleas have very distinct population genetic patterns from each other, likely due to the influence of other sympatrically-distributed vertebrates on the transmission of fleas.


Background
Comparative phylogeography aims to compare geographical patterns of evolutionary subdivision across multiple co-distributed species or species complexes, making it a powerful tool for elucidating the demographic and historical nature of intraspecific evolution [1,2]. Host-parasite systems represent a specific co-distribution phenomenon, in which the parasite's phylogeographic pattern is typically dependent on their host; comparative investigations of host-parasite systems have thus become one of the most promising directions of comparative phylogeography [3,4]. Several studies have shown congruent phylogeographic patterns between parasites and their hosts, such as nematodes (Heligmosomoides polygyrus) and field mice (Apodemus sylvaticus) [5] or mite (Spinturnix myoti) and bat (Myotis punicus) [6]. In contrast, other studies have shown either incongruent (e.g. lice (Polyplax arvicanthis) and rodents (Rhabdomys spp.) [7]) or somewhat incongruent (e.g. nematode (Heligmosomoides spp.) and mice (Apodemus spp.) [8]) phylogeographic patterns between parasites and their hosts. Moreover, many features of the host and its parasite can influence the overall phylogeographic or genealogic congruence, such as the intimacy of the interaction between the two organisms, the scale of vertical transmission of the symbionts through the host generations, and differences in life history traits (population sizes, generation times, and migrating abilities between the hosts and the parasites, etc.) [9][10][11][12].
Fleas belong to the order Siphonaptera; all species are obligatory blood-feeding parasites of endothermic vertebrates, including mammals and birds [13]. The life cycles of fleas consist of four stages: egg, larva, inactive pupa, and adult. All of these stages are highly dependent on the nests or burrows of their hosts. However, as a typical ectoparasite, all four stages can be viewed as free-living phases. Moreover, adult fleas may also leave their hosts and migrate onto new hosts [13,14], which potentially negatively influences the congruence of the phylogeographic pattern between fleas and their hosts. Gomez-Diaz et al. [15] compared the phylogeography of fleas (Xenopsylla gratiosa) and shearwaters (Calonectris spp.) and showed that neither genetic distances among host populations nor their spatial distribution explained the patterns of genetic variability observed in the fleas, likely due to a local adaptation to the sympatric host species or to a parasitic exchange during mixing among breeding populations of the hosts. In addition, Jones and Britten [16] compared the genetic structure between prairie dogs (Cynomys ludovicianus) and fleas (Oropsylla hirsuta); interestingly, the prairie dogs, but not the fleas, showed a significant isolation-by-distance pattern. Moreover, the estimated rates of gene flow among the flea colonies were higher than those among host colonies. Furthermore, these authors argued that other sympatric prairie mammals might be involved in dispersing the fleas, thus resulting in a lack of concordance between the population genetic structures of the host and ectoparasite.
Neopsylla paranoma (Ctenophthalmidae) is a Qinghai-Tibet Plateau (QTP) endemic flea [17]. The main host of N. paranoma is the plateau zokor (Eospalax baileyi), a rodent also endemic to the QTP. During a recent survey (Chen et al., Chinese Journal of Vector Biology and Control, in press), N. paranoma comprised 70% of the total fleas collected from trapped zokors. This species is reported as parasitizing sympatric animals (Himalayan marmots (Marmota himalayana), hamsters (Cricetulus longicaudatus), dipodids (Allactaga sibirica), plateau pikas (Ochotona curzoniae), and some birds), but comprise a very small proportion (<1%) of the total fleas collected from each of these hosts [18]. It should be mentioned that, as a typical subterranean rodent, plateau zokors spend almost all of their life in underground burrows with rather limited dispersal. Moreover, as a solitary species, zokors of the same gender rarely come into contact with each other; only during the mating season do males and females meet in temporary burrows for a very short time period (several minutes) [19]. Hence, unlike in the shearwaters and prairie dogs mentioned above, zokor-associated fleas have limited opportunities to transfer among zokor populations and between zokors and other sympatric animals. As a result, we hypothesized that the phylogeographic patterns of fleas would largely reflect the phylogeographic history of the zokors. Here, we sampled both zokors and fleas, and using mitochondrial DNA sequences, aimed to establish whether there was a high congruence in the phylogeographic patterns between the two species.

General information
A total of 130 zokor-flea pairs were collected from 17 locations (Table 1, including the GPS coordinates; Figure 1). In 12 locations, 7~12 pairs were collected. In the other 5 locations, only 1~4 pairs were collected because of a low parasitism rate for the fleas. The concatenated mitochondrial DNA (mtDNA) sequences of the zokors were 2,502 bp in length, successively consisting of complete sequences of six genes: tRNA-Cys, tRNA-Tyr, COI (cytochrome c oxidase subunit I), tRNA-Ser, tRNA-Asp, and COII (cytochrome c oxidase subunit II). The main segments were the COI and COII genes, which have 1,545 bp and 684 bp, respectively. The concatenated mtDNA sequences of the fleas were 2,407 bp in length, successively consisting of complete sequences of five genes: tRNA-Cys, tRNA-Tyr, COI, tRNA-Leu, and COII. Moreover, the main segments were also COI and COII genes, which have 1,536 bp and 681 bp, respectively.
From the 130 zokor sequences, a total of 287 variable sites and 55 haplotypes (Zh1~Zh55) were detected. No haplotype was shared by two or more populations. For the 130 flea sequences, a total of 157 variable sites and 75 haplotypes (Fh1~Fh75) were detected. Of these, there were two haplotypes shared by two (Fh64, shared by QL1 and QL2) and three (Fh25, shared by GN2, GN3, and XH) populations, respectively. The distributions of the haplotypes are listed in Table 1.

Population structure
AMOVA (Analyses of Molecular Variance, based on all 17 locations) showed that, for zokors, the percentage of variation between populations constituted 97.10%, while the within population variation constituted only 2.90%. For fleas, the percentage of variation among populations and within populations constituted 85.68% and 14.32%, respectively. All these AMOVA tests were highly significant (p = 0.000).
The Fst (fixation index) values and geographic distances among the 15 zokor/flea populations are listed in Table 2. A Mantel test showed that the correlation of (R 2 ) the zokor Fst with the geographic distance was 15.84% (P = 0.000), while the flea Fst values were also significantly correlated with geographic distance (p = 0.000). However, the R 2 value was 46.20%, much higher than that of the zokors. A Wilcoxon Signed Ranks Test revealed that the flea Fst values were significantly smaller than zokor Fst (N =105, Z = −7.712, p =0.000). In addition, a Spearman correlation analysis showed that the Fst values of zokors and fleas were also significantly and positively correlated (N =105, R =0.439, p =0.000).

Network and co-phylogenetic results
Both zokors and fleas had very distinct network features. The 55 zokor haplotypes formed into six distinct clades (Z1~Z6, Figure 2), corresponding to six geographic regions ( Figure 1). However, the 75 flea haplotypes formed only three clades: F1 included four populations while F2 included only one population and F3 included the other 12 populations related with the majority of the haplotypes (61 of 75) (Figures 1 and 3).

Incongruent genetic structure
Our results showed that the zokors and fleas had very distinct population genetic structures from each other. First, the AMOVA results revealed that the variation between the flea populations (85.68%) was considerably lower than that among zokor populations (97.10%). Second, flea Fst values were significantly smaller compared to the zokor Fst (p =0.000). Third, each zokor haplotype was uniquely distributed within a single population; for fleas, however, two haplotypes were shared between two or more populations. Fourth, the co-phylogenetic analyses were unable to reject the hypothesis of independence of speciation events between the two animals.
The results of the network analyses also showed a distinct geographic haplotype distribution pattern between the two species. The JZ2 population, which was an independent clade in zokors, showed a closer relation with its neighboring populations, JZ1 and BM. Also, the 12 populations (DT, GD, GH, HL, GN2, HY, HZ, MY, QL1, QL2, XH, and ZK), which were grouped into four separate clades in zokors, were grouped into a single clade in fleas (Table 1 and Figure 1). These results indicate that fleas have a lower level of genetic differentiation than zokors.

Effects of geographic barriers
The QTP is the largest and highest plateau on Earth. The QTP began a severe uplift during the Pliocene, and the uplift continued until the Quaternary glaciations. As a result, the plateau began to topographically diversify due to the intricate and continual development of new mountains and watercourses [20], which likely became geographical barriers that helped create various genetic structures for endemic animals [21]. It should be mentioned that with different dispersion abilities, each species of animal living on the QTP tends to have very distinct genetic structures. For example, there is a reported absence in distinct genetic structures among populations of a snow finch (Pyrgilauda ruficollis), which is a strong flier and tend to move seasonally. In contrast, a high level of genetic structure was reported for the ground tit (Pseudopodoces humilis; a resident bird species) and a lizard (Phrynocephalus vlangalii), because the two species tend not to disperse (reviewed in [22]). Like other subterranean rodents (such as the tiny tuco-tuco (Ctenomys minutus) [23]), geographic barriers seem to have played an important role in shaping the phylogeographic structure of zokors: the Yellow River separated Z1 and Z2 from Z3~Z6; the extremely high altitude at and around the JZ2 group isolated the haplotypes in this region into a basal clade Z1, whereas the Huangshui Basin separated Z4, Z5, and Z6 from each other. However, it is noteworthy that distance, rather than physical barriers, seems to have played an important role in shaping the geographic patterns of fleas: the R 2 of the IBD effects in fleas (46.20%) was much higher than that of the zokors (15.84%). The most evident example of the different effects of geographic barriers on population genetics occurred between the XH and GN2 groups. The geographic distance between the two populations was only 57.87 km, but the Fst between them was 0.99; as such, they were grouped into two different clades Z3 and Z6. For fleas however, the Fst between the two populations was only 0.14, and subsequently, they were grouped into a monophyletic clade F3. Taken together, these results suggest that the QTP-associated geographic barriers had a strong isolation effect in the genetic variation of zokors; in contrast, this effect was much weaker on fleas.
The third-part flea hosts?
As mentioned above, the differences in life history traits, such as population sizes, generation times, and migrating abilities, between the hosts and the parasites could directly or indirectly influence the interactions between the two organisms [10][11][12]. Although little is known about the ecology of the flea (N. paranoma), much more is known about zokors. Zokors spend most of their lives underground and have a very limited migration ability (only about 66.9 meters per year) [24]. As such, genetic interactions between the populations resulting from zokor migrations would also be very limited. In the   [19]. As a result, the migration of fleas via zokor contact should also be very limited. In other words, even if the fleas in our study have very distinct life history traits compared to the zokors that they live on, the phylogeographic features of fleas should still largely depend on the zokors, provided no other vertebrates were involved in this host-parasite system. Unexpectedly, however, we found that the fleas had lower inter-population genetic variation levels than zokors, suggesting that some other host animals have also mediated their migration and subsequently had a significant influence on their phylogeographic patterns. The degree of association between fleas and hosts varies among different host/flea species [25]. Zhang [18] reported that Himalayan marmots, hamsters, dipodids, plateau pikas, and some birds, all of which are good migrators, could also be infected by N. paranoma. It is puzzling, however, how zokors might exchange their parasites with above-ground animals. Pikas are the most closely related sympatric animals to plateau zokors. In many cases, pikas occupy the burrows of zokors [26], which might provide opportunities of flea exchange between the two species. Moreover, as a key species in the QTP [27], plateau pikas also share habitats with most of the other species in the QTP, as well. We hypothezed that the plateau pikas first serve as a mediator of fleas between zokors and above-ground animals (including birds), and then these above-ground animals help to disperse the fleas among different regions. It should be mentioned that there is a haplotype (Fh25) that is shared by GN2, GN3, and XH, which were separated by the Yellow River. Since a large majority of mammalian animals seem unable to migrate over the Yellow River, we suggest that birds might play a more important role for cross-river flea transmission.

Conclusions
Despite the unique subterranean lifestyle of the zokor, the population genetic patterns of the fleas were distinct from their hosts. Moreover, fleas had a lower level of genetic variation between populations compared to zokors, suggesting other sympatrically distributed vertebrates serve as a mediator for flea transmission.

Sampling
During April and May, 2012, individual zokors were sampled using ground arrows [28] from the east of Qinghai Province, China. The geographical information at each sampling site was recorded using an Etrex GPS unit (Garmin, Taiwan). Each of the trapped zokors was immediately put into a plastic bag filled with diethyl ether for ten minutes. The dead fleas were collected and stored in 75% alcohol. Muscle samples from each zokor were collected and fixed in 95% alcohol. The taxonomy of the fleas was identified under light microscopy. All animal work in this study was conducted with ethical approval from the Ethics Committee, Northwest Institute of Plateau Biology, Chinese Academy of Sciences. No specific permissions were required for the locations/activities in this study. The field studies did not involve endangered or protected species.

DNA extracting, amplifying and sequencing
Total DNA of each zokor sample was extracted using standard methods for animal tissue [29]. The sequences of the three mitochondrial DNA (mtDNA) segments were amplified by polymerase chain reaction (PCR) using three pairs of primers (Table 3) designed with reference to the plateau zokor mitochondrial genome (accession No. JN540033.1). Total DNA from each flea was extracted using a spin column kit (DNeasy tissue kit; Qiagen, Germany), according to the manufacturer's protocol. The sequences of three target mtDNA segments were amplified using three primer pairs (Table 3). Two were designed according to the mitochondrial genome of a mecopteran species (Boreus elegans, NC_015119.1), which was viewed  as the closest sister group of the fleas [30], and the other one is according to Whiting [30]. PCR amplifications were performed in total reaction volumes of 50 μL, containing 10mM Tris-HCl (pH 8.3), 1.5 mM MgCl 2 , 50 mM KCl, 100 μM of each dNTP, 0.25 μMol each primer (synthesized by Sangon Ltd., China), 0.5 μl of template DNA, and 2.5U Taq DNA polymerase (Takara). The reaction mixtures were denatured at 95°C for 5min and subjected to 35 cycles of 45 s at 95°C, 1 min at 50~55°C (depending on the primer pair), 1.5 min at 72°C, and a final extension step of 7 min at 72°C. PCR products were purified using a CASpure PCR Purification Kit, following the manufacturer's recommended protocol (Casarray, Shanghai, China). Purified DNA products were sequenced in both directions with the PCR primers on an ABI 373 automated sequencer. Both strands of each segment were sequenced using forward and reverse primers. Sequences were recorded in both strands with an overlap of at least 30%. The sequences of each segment were checked by eye and aligned using CLUSTAL W [31] and refined manually. Finally, the three segments of each species were combined for further analyses.

Genetic variation calculation
Basic polymorphism data (variable sites, haplotype distribution, etc.) were determined using DnaSP (version 5) [32]. The Arlequin (version 3.11) [33] software was used to calculate AMOVA (Analyses of Molecular Variance) and Fst values among populations. We then used ArcGIS software (version 9.3) to calculate the geographic distances among the sampling locations, and used Arlequin to carry out Mantel test [34] to examine isolation-by-distance effects of the two species (1000 permutations were used). We used a non-parametric 2 Related Samples Test (Wilcoxon Signed Ranks Test) in SPSS software (version 20.0) to compare the zokor/flea Fst values. We also used a Spearman correlation in SPSS to test the relationships between Fst values of zokors and fleas. It should be mentioned that, in two locations (GD and ZK) there was only one zokor-flea pair each, which could induce unrealistic Fst estimates such as negative values. Hence, for the analyses relating Fst values, these locations were excluded.

Network and co-phylogenetic analysis
Roehl data files were generated with DnaSP. The haplotype networks of zokors and fleas were constructed with the program NETWORK (version 4.6.1.2). The Median Joining method [35] was used to calculate network structure. Then the drawn networks were saved as bitmap and were introduced into CorelDraw (version X6) to prepare the final feature figures.
The optimized versions of Pierre Legendre's Parafit (i.e. AxParafit) [36] was used to test the co-phylogenetic congruence between zokors and fleas. First, the between group mean distance (Kimura 2-parameter model) matrices based on the 17 populations and 15 populations (GD and ZK excluded) were calculated with MEGA (version 5.2) [37]. The global-fit congruence and individual host-parasite interaction analysis were tested with 100,000 permutations using zokor distance matrix, flea distance matrix, and the zokor-flea association file. Each individual zokor-flea interaction is determined to be significant if either its ParaFit 1 or Parafit 2 p-value <0.05.