Comparative Phylogenetic Studies on Schistosoma japonicum and Its Snail Intermediate Host Oncomelania hupensis: Origins, Dispersal and Coevolution

Background Schistosoma japonicum causes major public health problems in China and the Philippines; this parasite, which is transmitted by freshwater snails of the species Oncomelania hupensis, causes the disease intestinal schistosomiasis in humans and cattle. Researchers working on Schistosoma in Africa have described the relationship between the parasites and their snail intermediate hosts as coevolved or even as an evolutionary arms race. In the present study this hypothesis of coevolution is evaluated for S. japonicum and O. hupensis. The origins and radiation of the snails and the parasite across China, and the taxonomic validity of the sub-species of O. hupensis, are also assessed. Methodology/Principal Findings The findings provide no evidence for coevolution between S. japonicum and O. hupensis, and the phylogeographical analysis suggests a heterochronous radiation of the parasites and snails in response to different palaeogeographical and climatic triggers. The results are consistent with a hypothesis of East to West colonisation of China by Oncomelania with a re-invasion of Japan by O. hupensis from China. The Taiwan population of S. japonicum appears to be recently established in comparison with mainland Chinese populations. Conclusions/Significance The snail and parasite populations of the western mountain region of China (Yunnan and Sichuan) appear to have been isolated from Southeast Asian populations since the Pleistocene; this has implications for road and rail links being constructed in the region, which will breach biogeographical barriers between China and Southeast Asia. The results also have implications for the spread of S. japonicum. In the absence of coevolution, the parasite may more readily colonise new snail populations to which it is not locally adapted, or even new intermediate host species; this can facilitate its dispersal into new areas. Additional work is required to assess further the risk of spread of S. japonicum.

Pleistocene; this has implications for road and rail links being constructed in the region, which will breach biogeographical barriers between China and Southeast Asia. The results also have implications for the spread of S. japonicum. In the absence of coevolution, the parasite may more readily colonise new snail populations to which it is not locally adapted, or even new intermediate host species; this can facilitate its dispersal into new areas. Additional work is required to assess further the risk of spread of S. japonicum.

Background
In China schistosomiasis in humans is caused by infection with the parasitic blood-fluke Schistosoma japonicum (Trematoda: Digenea). Schistosomiasis causes major public health problems in China and the Philippines [1]. Despite over 45 years of integrated control efforts, approximately one million people, and more than 1.7 million bovines and other mammals, are currently infected in mainland China [2]. The disease develops after direct penetration of the skin by free-swimming parasite larvae (cercariae) released from the snail intermediate host.
Exposure to cercariae of S. japonicum occurs through contact with wet vegetation or walking through rice fields or other flooded areas inhabited by infective snails. In the case of S. japonicum, only sub-species of the snail Oncomelania hupensis (Gastropoda: Pomatiopsidae) are able to act as intermediate host [3]. The snails become infected by similarly mobile and penetrative larvae (miracidia) passed in the feces of definitive hosts, which include a wide variety of mammals (up to 28 genera and 7 orders [1,4,5]).
Much of the transmission in mainland China occurs across the generally flat and marshy lake-land areas around the middle to lower Yangtze river, where Oncomelania hupensis is the snail intermediate host. In the highland areas of Southwest China (Sichuan and Yunnan) O. hupensis robertsoni is the intermediate host. Fewer people are infected in the highland areas; however, recrudescence is most marked in this region and in Sichuan between 1989 and 2004, the disease re-emerged in 8 counties (prevalence: 1.4% to 18.3% in humans, and up to 22.3% in cattle and 0.16% in snails) [6]. A third subspecies is found in China, O. hupensis tangi, which is endemic to hilly or coastal areas of Fujian, Guangdong and Guangxi provinces, although S. japonicum transmission appears to have been interrupted in these areas. In the Philippines, a further sub-species is responsible for transmission, namely O. hupensis quadrasi, and in Sulawesi O. hupensis lindoensis is responsible. Other subspecies exist in Taiwan and in Japan (where the only other species of Oncomelania, O. minima, is also found) but these are not, or are no longer, relevant to human health [7].
Clearly, S. japonicum continues to pose a serious public health problem and inter-disciplinary research is required to understand the patterns of transmission and persistence of the disease. Phylogeographical studies can shed light on the problem, particularly in determining where the disease comes from, explaining current distributions of the intermediate host and predicting future epidemiology. Comparative phylogenetics can help detect patterns of hostparasite coevolution and indicate any potential for regional adaptation. Despite the potential of such approaches, relatively little work has been done in this area for S. japonicum. The present study was therefore performed in order to apply current phylodynamic techniques to the estimation of sources and tracts of dispersal for this parasite and to test for the signatures of hostparasite coevolution during the evolution of S. japonicum and its intermediate hosts (the latter necessarily having the defining influence on the distribution of the parasite). Strictly, the techniques used here test for phylogenetic congruence rather than directly for coevolution. The absence of phylogenetic congruence would make long-term coevolution unlikely; thus providing an indirect test for the latter. The study aimed to utilise the largest and most geographically extensive data set available (Fig 1).
A geographical strain complex has been revealed within Schistosoma japonicum using partial DNA sequences of mitochondrial genes [22,23]. In 2005 microsatellites were used to assess variation among S. japonicum populations from 8 geographical locations in 7 endemic provinces across mainland China. The populations were found to fall into two main clades, a The locations of the sample sites in China and Japan. The coloured spots indicate areas from which snail and worm populations were sampled (red and green respectively). Most of the red spots lie on top of a green spot, indicating that the snails and the worms were sampled from the same locality. The coloured regions define biogeographical areas, within which there are no significant barriers to snail dispersal (e.g., no high mountains). Philippine samples (worms and snails) were included, but omitted from the map to increase resolution. Abbreviations are listed in Table 1. Plotted using the OpenStreetMap package in R. Sichuan/Yunnan clade and a lake/marshland clade of the middle and lower Yangtze river drainage [24]. The finding was later corroborated [25] using partial mitochondrial gene sequences (cytb-nd4L-nd4, cox3, nad1 nad4, nad5 and 16S-12S), with the additional observation that the lake/marshland clade included a highly diverse lower Yangtze sub-clade [25,26]. Molecular phylogenetic studies of Oncomelania hupensis are relatively few in number, but a  [13]. A study with similarly wide geographical coverage to that used here has been performed [27]; this employed DNA sequence variation in O. hupensis to indicate that isolation by distance was more significant in shaping genetic divergence than isolation by environment; however, this study which covered 29 localities, was population genetic (not comparative phylogenetic) and did not include the parasite.

Aims of the study
Although there have been past population phylogenetic studies of both Oncomelania and S. japonicum none has used as geographically diverse, taxon and character rich, data set as the present study. Also, this is the first quantitative comparison of the phylogenetic history of the snails and the parasites. Relatively few phylogenetic studies have considered the origins of these taxa, although Davis [28] proposed a Gondwanan origin for the Pomatiopsidae (including Oncomelania and proto-S. japonicum), with rafting to mainland Asia (via the Indian Craton after the break up of Gondwana) and colonisation of Southeast Asia and China (Tibet/Yunnan) via the northern-India-Myanmar, Brahmaputra-Irrawaddy river corridor in an West to East direction, around 18 Ma. In contrast, an East to West hypothesis has been proposed for the Chinese Pomatiopsidae taxa [29], with precursors of Oncomelania colonising originally from Australasia and via the Philippines, along island chains created by low sea levels and by tectonic movements (rafting). After reaching Japan, Proto-Oncomelania gives rise to Oncomelania hupensis in mainland China; the latter then recolonises Japan, the Philippines and Sulawesi (replacing antecedent forms). In a recent paper [30], the East to West hypothesis was tested and the radiation across China dated at 15-5 Ma (by molecular clock); however, the radiation of S. japonicum did not appear to be isochronous with that of the present day intermediate hosts [31,32]. Nevertheless, the history of host utilization in Schistosoma has been regarded as an evolutionary battle with snail defences [11,[33][34][35][36][37], with the schistosome under significant pressure to evolve counter measures to snail immune responses, or to track snail divergence in an evolutionary arms race (the 'Red Queen hypothesis' [38]). Schistosomiasis researchers postulating coevolution have evidenced this by citing restriction of the parasites to certain snail lineages [11], high levels of genetic variation in naturally exposed snail populations [34], and evidence for selection in schistosomes exposed to snails previously selected for resistance [36]. The latter study [36] did use genetic (microsatellite) variation to demonstrate higher rates of parasite divergence in resistant laboratory lines of snails; however, this was over a very short time-scale and the resistant populations had not evolved under natural conditions. The study also used the Schistosoma mansoni:Biomphalaria glabrata system (which involves far higher infection rates than seen with S. japonicum-see Discussion). Consequently, further studies are needed into the role of coevolution in the evolutionary history of Schistosoma.
The present investigation was performed in view of the lack of studies on S. japonicum comparing the parasite and snail intermediate host evolutionary histories, the alternative hypotheses regarding the origins of these taxa and the colonisation of mainland China, and the recent availability of additional DNA sequence data (improving geographical coverage). It is important to address these questions because they can help explain the current distribution of the parasite within the range of the snails, which is relevant to snail control strategies and the potential for range expansion, and to help assess the likely impacts of environmental manipulation such as the South-to-North-Water-Transfer project in China and the construction of road links. The Greater Mekong Subregion (GMS) Chengdu-Kunming corridor and GMS North-South Corridor from Yunnan, through Laos and to Thailand (NSEC1), both involve tunnels and rapid links through the mountains of Sichuan and Yunnan to Laos [39]. Consequently, the work will not only improve our understanding of host-parasite coevolution, but also shed light on the impacts of development in the region.

Methods Sampling
Three new samples of Oncomelania were taken and sequenced for this study; these were O. hupensis from Jiangshan in Jiangxi Province, China (cox1), and from Xinjin in Chengdu, Sichuan Province, China (cox1 and rrnL), and O. hupensis nosophora from Kanagawa Kiyokawa in Honshu, Japan (cox1 and rrnL)-these loci all belonging to the mitochondrial genome, with partial DNA sequences of each locus obtained using the popular Folmer [40] and Palumbi [41] primers, to allow comparisons with data from earlier studies. The three samples were identified on the basis of conchology, morphology and ecological habit (following [9]). The identifications were performed by YS and NH, two greatly experienced medical malacologists. The DNA sequencing followed procedures detailed elsewhere [30]. The remainder of the data were obtained from GenBank to give a total of 265 cox1 sequences and 70 rrnL sequences, apparently orthologous with the new sequence data above. Sequence data for Tricula hortensis from Sichuan China, also a pomatiopsid snail [42], were included for use as an outgroup. In addition to the data for Oncomelania, 14 complete cox1 and cox3 DNA sequences were obtained from the GenBank for samples from 14 Schistosoma japonicum transmission localities across China. A corresponding Schistosoma mekongi sequence was obtained as an outgroup taxon for the S. japonicum analyses. Full details of the data used and sampling range are given in Fig 1 and Tables 1, 2 and 3. The data set included all available relevant sequence data in the GenBank at the time of searching; however, some data (10 sequences) were excluded because of uncertain origin (e.g., DQ112269, DQ112270) and taxonomy (e.g., EU780224). The cox1 and cox3 loci were chosen following earlier work, which indicated that, for Schistosoma, cox3 was the locus of choice in terms of consistent phylogenetic signal and sufficient number of phylogenetically informative characters per site; this was followed by nad4, however, the present data showed a haplotype diversity <1 at this locus and cox1 was used because it gave high support values and "correct" topology in the same study [43]. Many of the localities listed in Tables 1 and 2 were not published either in GenBank or in the papers where the data were first presented. In these cases the localities were found by accessing field work reports cited in the paper presenting the sequence data (if present), matching published sample codes with associated accession numbers, reference to museum specimen accession numbers, contacting original researchers (or the local officials who accompanied them or those who arranged their collections) and referring to the personal observations of the present authors. In some cases, incomplete (e.g., to County level only) location data was given, or ambiguously transliterated place names, in these cases location records were completed and transliterations resolved. The locations in Tables 1 and 2 have been changed (where necessary) to use the closest place name (village etc) appearing on Google Maps; this is for the convenience of future authors.
After submission of this manuscript, new data for Philippine and Japanese Schistosoma japonicum were added to the GenBank by unpublished authors of another laboratory. These data are included in Table 3, but were not included in the hypothesis testing; nevertheless, a phylogeny was estimated using this extended data set in order to assess the impact of the new data on the conclusions of this study (if any).

Initial handling of data and selection of partitioning scheme and substitution models
Sequence data were extracted from GenBank using Biopython 1.61 Bio.SeqIO [44]. The sequences were aligned using MUSCLE 3.8.31 [45], with default settings. To reduce computing time in subsequent analyses, the alignments were inspected in SeaView 4.4.2 [46] and the first 240 and last 60 bps of the complete cox1 gene for the S. japonicum (worm) alignment were excluded (the ingroup showed no variation in these regions). Taxa with identical sequences at a locus (gene) were then excluded, leaving one representative of each haplotype: for the worms identical sequences only occurred within villages or townships, but for the snails identical sequences were found up to county level. In data sets where sequences for two loci were concatenated, data were only excluded if sequences were identical between two taxa at both loci. The final data set for the worms comprised a concatenated cox1+cox3 sequence, and had 16 taxa (19 in the extended data set) and 1994 characters (after removal of 2 identical sequences). As cox1 and rrnL data were not both available for all snail taxa, there were two data sets for the snails. The cox1 data set comprised 146 taxa and 599 characters (after removal of 129 identical sequences) and a second data set was made using pyfasta 0.5.2 to select all cox1 sequences for which there was a corresponding rrnL sequence; these sequences were then concatenated to produce a cox1+rrnL "both loci" alignment of 51 taxa and 1110 characters. The reading frame of the protein coding loci was determined using ExPASy Translate [47]. In addition to the alignments for individual sequences, population level data sets were produced for the snails because these were expected to be easier to visualise and detect dispersal tracts. To achieve this, the geographical range of the samples was divided into 10 biogeographical regions (see Fig 1) such that within each region there was no barrier to dispersal of Oncomelania; thus there would be only isolation by distance and no major ecological (e.g., lack of suitable habitat) or physical (e.g., mountain ranges or ridges) barriers; however it should be noted that the Tai Hu Plain (THP) region is likely to be interrupted by industrialised belts where there are no Oncomelania habitats. In contrast, each region is separated by mountains or similar areas of highland or sea. Consensus sequences were produced for the individuals within each region and a population sequence alignment estimated. The population data set included 13 ingroup taxa (JAP has two island populations, and the Philippines (PHL) was also included).
Phylogenetic analysis was performed using two fundamentally different approaches; the non-parametric Maximum Parsimony (MP) approach and the parametric Maximum Loci sampled are cox1 and cox3 (extracted from longer contigs, including complete mitochondrial genomes). Data are: GenBank accession numbers (GACs), collection localities (with coordinates) and source references; unpub, indicates unpublished data with locations drawn from isolate codes, specimen voucher numbers or details in the GenBank. ADL indicates adult worms from laboratory lines based on field-collected infected snails, FCC stands for cercariae from snails collected in the field and used directly to infect laboratory hosts. doi:10.1371/journal.pntd.0003935.t003 Likelihood approach. Two contrasting methods were used so that resilience of the hypothesis to changes in method and associated assumptions could be used to reveal weakly supported clades. The rationale was that any clade that was represented in phylogenies found by both methods (and well supported) would be a reliable reconstruction of phylogenetic history for these taxa. In addition, the strategy helps reveal artefacts associated with a specific class of methods.

Phylogenetic analysis by Maximum Likelihood (ML) with RAxML
RAxML 7.4.8 [48] was chosen to implement the ML analysis because, among currently available programs, RAxML has shown best performance in terms of inference times and final likelihood values, and has a rapid boot-strapping algorithm which allows clade support to be estimated in reasonable time frames, even when estimating null distributions. A series of test runs were used to determine optimum values or states for the settings of the RAxML analyses (details published elsewhere [30]). The apparently optimum partitioning strategy and evolutionary model for each partition was determined using PartitionFinder 1.0.1 [49], under a BIC criterion and models restricted to those implementable in RaxML. The chosen models for the snails were: cox1, GTR+G codons CP 123 ; cox1+rrnL, GTR+G codons CP 123 and rrnL separately partitioned; populations, GTR+G codons CP 111 cox1 and rrnL partitioned separately. For the worms the models were, for cox1+cox3, GTR+G codons CP 112 with cox1 and cox3 codons in the same partition (i.e., there were two only partitions). RAxML was run with 100000 bootstrap replicates, using every fifth tree found by bootstrapping as a starting tree for a series of 20,000 full ML searches. The CAT approximation in RAxML was not used. Three main runs were performed with different random number seeds. After checking that the independent runs led to trees of the same topology and very similar levels of support (±1%), the bootstrap trees for all three runs were pooled and a 100% (strict) majority rule consensus tree reported. For the hypothesis testing, a data set was constructed that included all snail taxa for which there was a worm sample taken at the same locality. These data included ten taxa and 1994 characters. The model and partition scheme for the worms was the same as for the full worm data set, and for the snails a two partition model was again chosen: GTR+G (rrnL, cox1codon1, cox1codon2) and cox1codon3.
Phylogenetic analysis by Maximum Parsimony (MP) with POY POY 5.1.1 [50] (parallelised build) was used to implement a Maximum Parsimony approach. The use of MP afforded an analysis free of the assumptions underlying ML methods, and the use of POY with its dynamic homology approach (where characters (transformation series) are inferred during the process of phylogenetic reconstruction) frees the analysis of any effects particular to the alignment inferred by MUSCLE [51]. A sensitivity analysis was used to choose the weighting (gap cost etc) and partition schemes for each data set, protocol published elsewhere [52].

Hypothesis testing
In order to test whether the evolution of the DNA sequences was consistent with a particular hypothesis, such as coevolution of Oncomelania hupensis and Schistosoma japonicum, the Incongruence Length Difference test [53] (ILD) was used as implemented in PAUP Ã 4.0b10 [54], with 5000 replicates. The test employed a sub-set of the data which included only polymorphic sites (for reasons published elsewhere [55]). The test, which randomly exchanges segments between the snail and worm data partitions, should give ML outcomes which are not significantly different from those of the original data if the snails and worms evolved to a common history. The ILD has been shown to be rather conservative when used as a test of topological congruence if phylogenies with trees of similar topologies are compared, with the opposite effect observed if the trees compared differ markedly in structure (e.g., internal branch length differences) [56,57]. Noise in the phylogenetic signal can also lead to type-I errors in the ILD test [58]. Consequently, the hypotheses were also tested using Monte Carlo simulation in the manner of the SOWH test [59]. In the case of the test for coevolution, the test statistic is the likelihood ratio of the phylogeny inferred for the worm data (unconstrained) and the same phylogeny inferred with the evolutionary history constrained to that of the snails (represented by the ML tree estimated for the snail data). A null distribution of the test statistic was calculated by simulating many data sets using Seq-Gen 1.3.3-1 [60] and the ML parameters of the substitution model inferred for the worm data, but constrained under the null hypothesis (the ML tree for the original worm data constrained by the snail ML tree). For each simulacrum the ML tree was estimated both unconstrained and constrained by the snail ML tree, and a likelihood ratio computed. The null distribution then being a distribution for the amount of discord expected to occur when the worm phylogeny had evolved according to the same history as the snails. If the likelihood ratio for the original data falls above the 95 th percentile of the null distribution, the hypothesis that both worms and snails evolved to the same (i.e., the worms') phylogenetic history can be rejected at P<0.05. The replicates were performed using RAxML (with settings as for the original worm data, but bootstrapping set to terminate according to a convergence criterion based on the extended majority rule consensus trees), and they continued until the null distribution appeared to have stabilised, as judged by a plateau of t-values with increasing replicate number.

Visualisation of phylogenies and phylogeographies
In cases where the topologies resulting from phylogenetic analysis by ML and MP did not agree for a particular data set, a strict consensus tree was generated from the two trees so that discordant clades were represented by polytomies. Consequently, the resolved clades in the final trees for each data set were those supported by both methods. In order to visualise the phylogeographies of both snails and worms, the phylogenies were mapped in Marble Virtual Globe 1.8.3 using the Phylo2GE R script. In addition, phylogenies were compared (topologically) using the compare2Trees algorithm [61], which scores each pair of branches, between the trees, according to the common taxon partitions they define, with the branches then paired so as to maximise the overall score; this process yields a score (S) value for a pair of trees.

Phylogenetic reconstructions
To enable RAxML and reduce computing time 129 identical and/or ambiguous cox1 sequences were excluded from the original snail data set (i.e., the sequence alignment for Oncomelania) to leave 146 taxa. Only one of the 129 exclusions involved identical sequences at the inter-regional scale. The haplotype diversity by region was roughly inversely proportional to the sampling intensity (i.e., number of sequences per region), except for YEB and PHL, and to a lesser degree SCB, which had small sample sizes and low haplotype diversities (Table 4).
Replicate phylogenetic analyses run for the snail data, with different random seeds and only the cox1 data, failed to result in a common tree (S<0.88), and there was poor agreement between the results of the ML and MP searches (S<0.55); the phylogeny also contained many unresolved clades. Consequently, the cox1 data were considered to lack phylogenetic signal and were not considered further in this study. In contrast, the phylogenies estimated for both loci (cox1 and rrnL concatenated) showed good agreement among replicate runs (S>0.92) and between ML and MP (S>0.75).
The strict (100%) consensus tree between the ML and MP searches is given in Fig 2. Considering O. hupensis, the only biogeographical regions characterised as monophyletic clades are JAP and YEB. Nevertheless, the Sichuan populations (SCB and SAV), which lie in an isolated mountainous area, form an unresolved near-monophyletic clade, except for three SAV individuals which form a further unresolved clade at the base of the O. hupensis clade (this could be a result of long branch attraction due to saturation or a lack of apomorphies in younger clades [62], with slippage of long branches leading to SAV taxa, towards the root of the tree). The YEB populations, also mountainous but separated from the Sichuan populations by the Hengduan Range, formed a distinct basal clade at the same level as the three extraneous SAV samples. Thus, the basal clades of the phylogeny are occupied solely by O. h. robertsoni (and O. minima). The remaining major clade, which appears derived from O. h. robertsoni contains all the other O. hupensis samples, including the non-Chinese taxa O. h. quadrasi (Philippines), which is basal to the clade, and O. h. nosophora (Japan)-suggesting that these taxa did not diverge in situ. Although O. h. nosophora is monophyletic and O. h. robertsoni is near monophyletic (it's apparent polyphyly might be explained by long branch attraction), O. h. hupensis is polyphyletic and includes O. h. tangi and O. h. formosana (of FCP and Taiwan, respectively); this suggests that the latter two may be populations of O. h. hupensis. The GCP taxon lies at the base of the Chinese O. h. hupensis clade and this gives some support to the case for O. h. guangxiensis. The populations of the lake and marshland region (DLB, FCP and THP) form an unresolved clade, suggesting that there are few barriers to migration (gene-flow) between them. In view of this, the population phylogeny (where data for individuals is pooled to provide consensus sequences for each region) provides a representative summary of the phylogeny (Fig 3). The population phylogeny agrees with that for both loci, based on individuals, except that it shows GCP as forming a sub-clade, of the Lake and Marshland, Taiwan and Japan clade, with PLB and DLB.
The phylogeny estimated for the worms (Fig 4) differed from that of the snails in certain key features. For example, some DLB and THP populations are basal to the western mountain clades (YEB, SCB and SAV). As in the snail phylogeny, the Taiwan population falls within a Table 4. Numbers of individuals sequenced (N) and haplotype diversity (H) by region for the Cox1 and rrnL data sets (snails) and the Cox1+Cox3 data set (worms).

Cox1 rrnL Cox1+Cox3
Region clade comprised solely of Lake and Marshland mainland Chinese taxa; however, this clade excludes DLB and so contains only THP, PLB and TAI, also the Taiwan sample clusters with PLB forming a clade derived from the Chikou THP samples.

Hypothesis testing
The snail and corresponding worm phylogeny (i.e., that estimated in exclusion of taxa not held in common; Fig 5) showed a low level of correspondence (S = 0.46). Consequently, it was necessary to test the null hypothesis that the snails and worms had evolved to a common evolutionary history, i.e., that of the snails. Initially the ILD test was used to detect significant conflict in the phylogenetic signals of the snail and worm data. The test was significant (P = 0.00004) suggesting that the two data partitions are the result of different evolutionary processes. To test further the null hypothesis, a SOWH test was performed; this resulted in a likelihood ratio for the observed data (unconstrained / constrained by the null hypothesis tree) of 75.52 and a 95 th percentile for the null distribution of 4.26. Consequently, the null hypothesis can be rejected in the light of these data (P<0.0001). In view of these findings it appears highly unlikely that the evolutionary radiation of Schistosoma japonicum across China was shaped or driven by that of the snail intermediate hosts.

Phylogeography
Geospatially projected phylogenies (Fig 6 and Files A and B in S1 Dataset) assist in phylogeographical interpretation and in the present study they reveal clear differences between the radiations of the snails and the worms. The map for the snails (Fig 6A) suggests initial colonisation of the valleys of the western mountains (SAV, SCB, YEB) by a proto-Oncomelania hupensis robertsoni; thereafter these lineages, established in their respective valleys and basins, appear to have stabilised after initial divergence (no further cladogenesis), and remained so throughout most of the history of O. hupensis. The mountain clade appears to have given rise to the Lake and Marshland and East-coast clades (THP, FCP) of Chinese O. h. hupensis, with radiations back into Japan (as O. h. nosophora) and to Taiwan and the Philippines more recently in the history of this species. A second, slightly more recent, radiation from DLB, west to GCP and eastwards to PLB, then occurs. In contrast the worms show a history where the Lake and  Marshland (eastern) populations appear to involve cladogenic events that occur throughout the history of Chinese S. japonicum, whereas the western mountain taxa are stable following initial establishment in their respective valleys. Unlike the snails, the worms show a most recent colonisation event of Taiwan that is associated with the PLB region rather than THP/FCP. After submission of this manuscript, new data for Philippine and Japanese Schistosoma japonicum were added to the GenBank. In order to determine the position of these additional populations in the phylogeography and their congruence with the snail phylogenies, an extended data set was subjected to phylogenetic analysis. File C in S1 Dataset and S1 Fig give the resulting kml and an image showing the extended phylogeny projected onto a globe, respectively.

Taxonomy and phylogenetics
Of the four Chinese Oncomelania hupensis sub-species described by Davis[9], only two are supported in the present study. Oncomelania h. robertsoni is not polyphyletic and all    [11]; therefore it is possible that all three sub-species in the large derived clade indicated by phylogenetic analyses for both loci and for snail populations (Figs 2 and 3) are in fact O. h. hupensis. The phylogeny estimated for the worms, in which the lowland populations form a clade distinct from those of the highland populations (SAV, SCB and YEB), is consistent with morphological, host-utilisation, and maturation rate observations that suggested independent lineages or sub-species for the highland and lowland S. japonicum in China [63]. Nevertheless, O. h. hupensis appears paraphyletic with two DLB (and one THP) populations forming separate clades near the base of the tree. Divergence of some Hunan/Hubei (DLB) populations away from those of other Lake and Marshland regions could result from the particularly long history of intensive control efforts [4] in these provinces, with slippage of these clades towards the outgroup owing to long branch attraction coupled with a lack of apomorphies among younger clades [64]. It is also interesting to note that the zoophilic strain [65] from Taiwan is not genetically distinct from those strains capable of infecting humans, although it is one of the two most derived members of the lake and marshland clade.

Hypothesis testing
The lack of concordance between the snail and worm phylogenies found in the present study could result from heterochronous evolution of the host and parasite in response to different palaeo-geographical or climatic environments. In the closely related taxon, Schistosoma mekongi, the radiation of the parasite (dated 2.1-1.0 Ma) was shown to be independent of that of its intermediate host . Divergence events among the snails were considered to be a concerted response to the final Indosinian orogeny around 5 Ma, with S. mekongi colonising the snails across its present range much later (early Pleistocene) [31,32]. Using a molecular clock the introduction of O. hupensis across mainland China has been dated to the early Miocene (c.a., 22 Ma), with high rates of cladogenesis 8-2 Ma and linked to the exceptionally warm and humid climate in the region at that time and tectonic upheaval in Japan [30]. The divergence of the Schistosoma japonicum clade has been dated at 4.6 Ma [32]; this, implies that the radiation of O. hupensis occurred before that of S. japonicum. If the radiation of the snails and worms is heterochronous there is no opportunity for coevolution; the implication is also that ancestral intermediate hosts differed from those of the present, which again makes coevolution unlikely.
Coevolution might be expected in Schistosoma species such as S. mansoni, which infect snail populations at relatively high prevalence and achieve high rates of cercariogenesis [66], but seems unlikely in S. japonicum because of its lower prevalence in the intermediate host populations. The prevalence of natural infections in China ranges from 0.038% (Jiangsu in 2011) to 7.8% (Anhui in 2013) [67]. In cases where the snails experience a low probability of becoming infected, they are under little pressure to invest resources in defence [68]. In contrast, prevalences as high as 75.7% have been reported for natural populations of S. mansoni in Biomphalaria glabrata in Brazil [69]. Factors such as the generalised nature of gastropod immune systems, and evidence for frequent host switches in the parasites' evolutionary histories, also make a "Red Queen" scenario unlikely for these schistosomes [32,68,70,71]. Consequently, molecular clock dating for other members of the S. sinensium clade and their intermediate hosts (also close relatives of Oncomelania) and the low prevalence of infection in O. hupensis suggest that the lack of evidence for coevolution found in this study is to be expected.
It could be argued that comparison of naturally infected (and infective) snails, rather than snails from the local populations transmitting the parasite (but not necessarily themselves infected) would be more likely to reveal signs of coevolution, i.e., there could be sub-populations or sub-strains of snails that have been coevolving with the worms, rather than the general snail populations sampled in this study. The existence of such sub-populations is questionable, as is the existence of resistant and susceptible lineages of snails in schistosomiasis. It has been shown that any snail taken from a natural population of B. glabrata will become infected if exposed to enough miracidia from a natural S. mansoni population [72]. Consequently, a "resistant" snail line is merely a sub-population selected to be discordant with the epitopes expressed by a particular Schistosoma line. Even authors working on the B. glabrata-S. mansoni association, have observed a complete reversal in resistance phenotype after a few laboratory generations and note that genotypic responses would only occur in associations where prevalences are high [36]. In view of this, it is unlikely that resistant sub-strains of O. hupensis exist and infection probably occurs more at random across the general snail population (perhaps influenced by ecology and spatial coincidence). Nevertheless, it would be interesting to repeat the present study using naturally infected snails from the localities studied and seek to detect any signs of coevolution.

Phylogeography
As mentioned above, earlier studies date the introduction of O. hupensis to mainland China at around 22 Ma, at which time the region was significantly less mountainous, and three general clades appear to have been rapidly established; these span China, with the O. h. robertsoni clade in the far West, the Dongting Lake Basin (DLB) clade in the center of the Lake and Marshland region, and the Tai Hu Plain (THP) clade near the East coast of mainland China (Fig 6A). O. h. robertsoni appears to have diverged little since its initial colonisation of the mountain valleys in which it is found today, and its lineages have probably been isolated therein since the second major uplift of the Himalaya about 7 Ma [73]. The other two clades appear to have undergone two successive, more recent, cladogenic events, which in the case of the DLB clade, gave rise to the GCP and PLB populations (to the West and East, respectively; Fig 6A node 3). The THP clade gave rise to FCP, Taiwan, Japanese and Philippine populations of O. h. hupensis (Fig 6A  node 2). Such a recolonisation of Japan by mainland Chinese O. hupensis is consistent with the "East to West" hypothesis [29]. The divergence of the western clades is likely to have occurred around 8 Ma when the Himalayan uplift altered global climate and triggered increasing aridity in the region, which would have fragmented existing Oncomelania populations [30]. Interestingly, the Taiwan population is also included in the lake and marshland clades even though Taiwan has been separated from the mainland since Pleistocene [74,75]and the more recent divergences within O. hupensis occurred before 2 Ma. It is possible that tsunami events could have exchanged snails between the mainland and Taiwan. Indeed, the March 2011 Pacific tsunami demonstrated that large aggregates of material may cross even oceanic distances in less than 15 months and that freshwater pools on these may harbour viable communities of exotic aquatic organisms (including molluscs) [76]. Oncomelania is also capable of aestivating out of water for several months. In addition, intermittent land bridges occurred, linking Taiwan during the Quaternary [77]. The relative lack of genetic variation among the Taiwan populations also suggests a recent colonisation of the island (or extinction of long established lineages followed by recent recolonisation). The Guangxi Plain (GCP) populations formed a clade distinct from other O. hupensis and may have been isolated from the other O. hupensis taxa since the late Miocene/late Pliocene by uplift along the margins of the Youjiang Basin (Jiangnan range) [78,79].
Although the radiation of S. japonicum is described above as occurring around 4 Ma later than that of the snails, the western (mountain) clades of the parasite still show the same initial divergence and then absence of cladogenesis as do the snails. The ancestral hosts of S. sinensium group parasites appear to be rodents (especially Rattus and its sister group) and it is therefore likely that S. japonicum radiated in China in concert with the Pliocene radiation of Rattus which began in Southeast Asia [80]. After colonising the valleys of Sichuan and Yunnan in rodents radiating into China from Southeast Asia, the parasites would become isolated in these valleys by cooling and increasing aridity in the Pleistocene [81]; thus suppressing further cladogenesis. In contrast, the lake and marshland clades undergo a series of cladogenic events spread from around the early Pliocene towards the Recent. Initially DLB and THP clades are established, together with a second THP clade that is derived from the O. h. robertsoni clades ( Fig  6B node 1), this is followed by progressive diversification of one branch of the second THP clade (Fig 6B node 2), whilst the DLB-associated THP clade remains stable. The second THP clade most recently gives rise to an ancestral form that diversifies into a PLB and a Taiwan clade (Fig 6B node 3). The possibility of a long-distance dispersal from Sichuan in the western mountains to THP near the East coast of China (Fig 6B arrow) is an intriguing one; however, the possibility of misidentification or laboratory error concerning Genbank deposited sequences must also be considered where highly inconsistent relationships are found The introduction needs to be dated in future work and might be related to traffic down the Yangtze river (human activities) and the extensive cladogenesis and spread of the Sichuan strain in the coastal lowland areas is an unexpected event, which might relate to better development of the parasite in naïve human hosts (or presence of a more dynamic host population than in the mountain regions).
Analysis of the extended data set (S1 Fig) shows the Japanese S. japonicum arising from the same basal THP lineage as the Taiwan population. In contrast the Philippine clade arises from the basal YEB clade along with taxa from Sichuan. Although this relationship appears to mirror that of the snail phylogeny, it is a relatively recent divergence in the parasites and a relatively earlier one in the snails; thus it does not increase the degree of phylogenetic congruence between the hosts and parasites.
The possibility of extensive radiation and dispersal, after long-distance introduction of a Sichuan strain of S. japonicum to the coastal region, is important in view of the fact that the mountains of Yunnan and Sichuan appear to have formed a barrier to dispersal of O. h. robertsoni transmitted S. japonicum for perhaps millions of years. The observation is particularly relevant in the context of the South-to-North-Water-Transfer project (Eastern Route) which will transfer water from endemic areas in Jiangsu province to areas in Shandong Province, towards the Yellow river, where O. hupensis has yet to be found but where conditions appear to be favorable for transmission [82]. Oncomelania (and the associated schistosomiasis) are most widespread in the lake and marshland areas of the middle and lower Yangtze river drainage; the distribution of snail and parasite in the mountainous areas of Sichuan in more patchy, and in Yunnan they are found only in a restricted area around Dali [83]. As the GMS road projects will breach the mountain ranges between Sichuan and eastern China and Sichuan and Yunnan (and further South into Laos and Thailand), it is important to understand the colonisation history of S. japonicum.

Conclusions
The present work has led to the rejection of the hypothesis of coevolution for Schistosoma japonicum and Oncomelania hupensis on the basis of the samples available. The finding is consistent with models regarding the relative timing of the radiations of the two groups proposed in earlier studies, and with observations of a low prevalence of infection in these snails. Nevertheless, it is still possible that the parasites show some adaptations to a snail population through which they have been cycling for some time, but the analysis does suggest that the worms are not highly evolved/restricted to a particular sub-species or strain of snail. Consequently, host-switching or acquisition is more likely than would be implied by a Red Queen scenario. The findings also suggest that O. h. formosana (of Taiwan) and O. h. tangi (of Fujian) might be O. h. hupensis and not distinct sub-species. The phylogeographical reconstructions suggest that at least one long-distance dispersal event occurred across China between the western mountain populations of S. japonicum and the East coast region. The event appears to have triggered extensive cladogenesis and dispersal (including to Taiwan) on the East coast. Consequently, further work is necessary to confirm and to date this long-distance dispersal and to detect any further such events, so that their origins and driving forces can be determined. Further work is also required with a richer data set as support for some of the clades indicated in the analyses was less than 90%. The findings are particularly important in view of the infrastructure development plans which will breach the mountain barriers between Sichuan, Yunnan and Southeast Asia. The results also have implications for the spread of S. japonicum as, in the absence of coevolution, the parasite may more readily colonise new snail populations to which it is not locally adapted, or even new intermediate host species, and this can facilitate its dispersal into new areas. The work also lends support to the East-West hypothesis for the origin and dispersal of Oncomelania and the Pomatiopsidae.