INTRODUCTION
Schistosomiasis of wildlife is known to be responsible for human cercarial dermatitis (CD) throughout Europe (Chamot et al. Reference Chamot, Tiscani and Rougemont1998). Cercarial dermatitis occurs as the result of an acute allergic reaction caused by prolonged repeated exposure to water infested by non-host-specific schistosomatid cercariae released from freshwater snails (Kolárová, Reference Kolárová2007; Horák and Kolárová, Reference Horák and Kolárová2010). The cercariae penetrate the skin inducing a major inflammatory response presenting in maculo-papulo-vascular eruptions followed by intense itching, which can often lead to secondary infection as a result of scratching, swelling of lymph nodes, erythema and oedema (Lichtenbergova et al. Reference Lichtenbergova, Kolbekova, Kourilova, Kansy, Mikes, Haas, Schramm, Horak, Kolrova and Mountford2008). Experimental studies have also shown that in immuno-deficient or naïve non-specific hosts, the parasites can migrate to lungs or the nervous system causing pathologies beyond those associated with the skin (Kolárová, Reference Kolárová2007; Horák and Kolárová, Reference Horák and Kolárová2010).
Outbreaks of CD across Europe have had significant impact on tourism and local economy, with foci of infection in Germany, Austria, Switzerland, the Netherlands, Czech Republic, France, Poland, Hungary (Chamot et al. Reference Chamot, Tiscani and Rougemont1998). The majority of CD has been attributed to the cercariae of avian schistosomes, however, Majoros et al. (Reference Majoros, Dán and Erdélyi2010) identified natural foci of Schistosoma turkestanicum (syn. Orientobilharzia turkestanicum), a known but neglected agent of CD causing public health concern throughout Asia (Wang et al. Reference Wang, Chen, Zhao, Chen, Zhai, Li and Zhu2009; Tabaripour et al. Reference Tabaripour, Youssefi and Tabaripour2015), utilizing red deer (Cervus elaphus) in the Gemenc Forest region of Hungary as the definitive host. The Gemenc Forest is an area frequented by holiday makers for recreational water sports and fishing as well as commercial hunting, and Majoros et al. (Reference Majoros, Dán and Erdélyi2010) highlighted a long history of CD in the area, which had been assumed to be caused by avian schistosome cercariae.
The discovery of S. turkestanicum in Hungary raised concern of the adaptability of such ‘tropical’ parasites to establish in temperate environments. However, detailed molecular analysis of Hungarian S. turkestanicum indicated that it was a genetically unique lineage distinct from Asian isolates although most closely related to isolates from Iran (Lawton and Majoros, Reference Lawton and Majoros2013). It was hypothesized that the divergence of the Hungarian population coincided with the migration North of deer populations from the Middle East during the course of the last ice age from 110 000 to 10 000 years ago (Lawton and Majoros, Reference Lawton and Majoros2013). However, such studies on Hungarian S. turkestanicum have been undertaken on small sample sizes unable to account for genetic diversity across the parasite population or the impact that infections in individual hosts may have upon such genetic diversity; crucial in understanding the historical events of the establishment of these parasites population.
As highlighted by Brusentsov et al. (Reference Brusentsov, Katokhin, Brusentsova, Shekhovtsov, Borovikov, Goncharenko, Lider, Romashov, Rusinek, Shibitov, Suleymanov, Yevtushenko and Mordvinov2013), detailed studies on the genetic diversity of parasitic organisms can provide a deeper understanding of the genetic structure and phylogenetic history of parasite populations, providing insights into the ecological and evolutionary dynamics of these organisms, often crucial for monitoring and control. Therefore, the aim of this study was to measure genetic diversity of Hungarian S. turkestanicum within and between hosts to disentangle the evolutionary history of this parasite's emergence into Hungary by highlighting molecular characteristics, which could differentiate a long endemic history or recent invasion event. The approach taken in this study was used to provide an insight into the origins and radiation of S. turkestanicum in Hungary in the absence of detailed host provenance information. Sequences of the mitochondrial cytochrome c oxidase (cox1) gene and nuclear ribosomal internal transcribed spacer region (ITS) gene fragment were used to test four distinct but interlinked hypotheses: (1) owing to long divergence times between S. turkestanicum populations from China, Iran and Hungary, all cox1 and ITS haplotypes in the Hungarian S. turkestanicum isolate will be unique to Hungary, and intermediate forms of these haplotypes will be seen between Iranian and the most divergent Hungarian haplotypes; (2) there will be genetic diversity within and between S. turkestanicum subpopulations from different individual deer hosts; (3) the Hungarian S. turkestanicum cox1 and ITS haplotypes will illustrate fine population structure and show demographic stability; (4) due to the long establishment of these populations, positive selection would be driving the variation seen within and between worm populations.
MATERIALS AND METHODS
Ethics statement
Adult S. turkestanicum were collected from deer that were shot by licensed hunters who annually cull the deer within the Gemenc Forests for meat and pelts. Only livers donated by hunters for inspection were used. This material forms part of a long-term monitoring programme on deer health in relation to the prevalence of liver flukes, particularly the invasive species Fascioloides magna.
Sample collection
Adult S. turkestanicum were recovered during necropsy from livers of six individual red deers (C. elaphus) based on dissection, sedimentation and microscopy detection techniques as described by Majoros et al. (Reference Majoros, Dán and Erdélyi2010). It was not possible to check any other organs for blood flukes as only livers had been submitted for a deer health survey by hunters during the annual hunting seasons from 2012 to 2013. The livers were extracted from apparently healthy deer that were shot within the Gemenc Forest Hungary, which encompasses over 95 square miles of gallery forest and inundation areas stretching along the Danube River.
DNA extraction, PCR amplification, sequencing and sequence alignment
From the six deer livers, a total of 56 paired worms were recovered of which only males were used for further analysis to avoid the potential contamination of multiple nuclear genotypes being amplified from ova held within the females. Genomic DNA was extracted from each of the male S. turkestanicum worms individually using the Qiagen DNeasy tissue kit (Qiagen Inc.) as specified by the manufacturer's instructions. Both the mitochondrial marker, the cox1 subunit gene, and the nuclear marker, the ribosomal ITS, were amplified through polymerase chain reaction (PCR) using the primers and methods originally used by Lawton and Majoros (Reference Lawton and Majoros2013). Only 5 µL of each reaction was used to visualize resultant amplicons on a 1% agarose gel stained with GelRed (Bioline). The remaining 20 µL PCR products were sequenced at the Natural History Museum, London, using fluorescent dye terminator sequencing kits (Applied Biosystems), then run on an Applied Biosystems 3730XL DNA analyser.
Resultant DNA sequences were assembled and edited using Bioedit (Hall, Reference Hall1999) and initial BLAST searches performed at NCBI (http://www.ncbi.nlm.nih.gov/) against the GenBank sequence database confirmed all generated sequences to represent S. turkestanicum. Assembled sequences were submitted to GenBank (Cox1 Acc:MF144970-MF145024; ITS Acc: MF145025-MF145080) and two main alignments were constructed using the MUSCLE sequence alignment tool (http://www.ebi.ac.uk) combining all available cox1 and ITS sequences for S. turkestanicum with the data generated from this study. Five unique cox1 sequences from previously published data were included from China (Acc: EU177877, EU221230–31, EU436657, HQ283100) and five from Iran (Acc: AY157200, KC456231–34) and a further 10 sequences generated S. turkestancium from a single deer, also C. elaphus, from the Gemenc Forest (Acc: JX467179-88). This took the total sample size for cox1 up to 75 individual sequences across three localities and a total of seven deers sampled. Although it was possible to obtain the corresponding ITS sequences for the previously published cox1 sequences from the Gemenc Forest (JX56189-98), it was not possible to do this for representations from the other countries. However, four unique sequences from China were added to the alignment (Acc: EU436663-66) and a single available sequence from Iran (Acc: HM803240) were incorporated into the ITS alignment taking the total number of ITS sequences used for analysis to 70 individual worms.
Genetic diversity and phylogenetic analysis of S. turkestanicum populations
Standard population diversity indices, including number of segregating sites (S), number of haplotypes (h), haplotype diversity (Hd), nucleotide diversity (π) and average number of pairwise differences, were estimated for the cox1 and ITS alignments within and between isolates using DnaSP5 (Librado and Rozas, Reference Librado and Rozas2009). To identify the occurrence of any significant genetic differentiation between S. turkestancium populations, analysis of molecular variance (AMOVA) was calculated using Arlequin (Excoffer and Lischer, Reference Excoffer and Lischer2010) for both the cox1 and ITS alignments using country of origin to specify samples for between-group calculations, where genetic distance was considered between haplotypes and their frequency between the populations defined between samples, with statistical significance determined by 1000 permutations. Finally, to provide a phylogenetic insight into the structure of the S. turkestanicum population, the most parsimoniously informative haplotype networks were constructed for both cox1 and ITS alignments in TCS 1.13 (Clement et al. Reference Clement, Posada and Crandall2000).
Historical demography and interhost variation
To infer potential population expansion events, the distribution of pairwise genetic differences (mismatch distribution) analysis was conducted across all isolates and specifically for the Hungarian isolates. This analysis was performed on both cox1 and ITS alignments using DnaSP 5, where distribution of frequency of pairwise differences between sample sequences as a measure of mismatch distribution was performed for 1000 simulations. The raggedness (r) statistic was used to measure the smoothness of the observed mismatch distributions and to test for significant departure from unimodality. For comparison to the r statistic, both Fu's F S and Tajima's D were also calculated using DnaSP5 for both cox1 and ITS alignments to obtain P values to identify any significant departure from neutrality, which could also be indicative of recent population expansion events.
In order to measure genetic diversity of the parasite within and between hosts, nucleotide diversity (π) of the worms from the seven hosts was calculated for cox1 and ITS sequences as described above and a Wilcoxon signed rank test was employed to calculate any significant differences between within-host diversity of cox1 and ITS for populations from all seven hosts. The frequency and distribution of identified cox1 and ITS haplotypes were calculated between hosts in order to highlight any host specificity of parasite genotypes. Finally, to evaluate the presence of any genetic differentiation between parasite populations within individual hosts, AMOVA was again implemented for both the cox1 and ITS sequences with each host being specified as unique samples and thus independent populations. As described previously, genetic distance was considered between haplotypes and their frequency between the populations defined between samples, with statistical significance determined by 1000 permutations.
Detection of signatures of selection
To detect the occurrence of selection at the nucleotide level within and across the different populations of S. turkestanicum the divergence of non-synonymous (Ka) and synonymous (Ks) mutations was estimated in the cox1. This was followed by sliding window analysis with a step size and window size of 3, which was employed to illustrate the contribution of Ka and Ks mutations to nucleotide diversity across the S. turkestanicum cox1 and was also performed using DnaSP. The contribution of both types of mutations to divergence of the cox1 was calculated using MEGA6 as the frequency of total non-synonymous and synonymous in relation to pairwise distance. In order to assess selection, the typical Ka/Ks ratio (ω) was employed where the value of ω indicates the mode of selection with ω = 1 indicating sites under neutrality, ω < 1 indicating the occurrence of purifying selection and ω > 1 signifying the presence of positive selection. Selection across the mitochondrial genome often occurs in a non-random manner within specific portions of genes or gene complexes associated with metabolic requirements imposed from the environment. Thus, a few amino acid substitutions can impose significant changes on the kinetics of reactions, which are catalyzed by mitochondrial proteins. Therefore, changes in physicochemical properties caused in response to amino acid replacements across the S. turkestancium population cox1 phylogeny were evaluated in TREESAAP 3.2 (Woolley et al. Reference Woolley, Johnson, Smith, Crandall and McClellan2003).
RESULTS
Diversity and divergence of cox1 and ITS between geographical isolates
Owing to the over-representation of sampling from Hungary relative to Iran and China, no major results can be presented in the context of the frequency of haplotypes relative to each population from geographical locations. Therefore, the data presented are primarily based on the nuclear divergence of sequences and the sharing of specific haplotypes between localities. As illustrated by Lawton and Majoros (Reference Lawton and Majoros2013), the cox1 showed consistently higher divergence levels than the ITS across all isolates with higher numbers of singleton mutations, a greater number of haplotypes, higher nucleotide diversity and increased average pairwise differences (Supplementary Table S1). Despite the low number of representative sequences from China and Iran compared with Hungary, the genetic diversity was relatively high with each sequence representing a unique haplotype and having comparable haplotype diversity, nucleotide diversity and pairwise difference to that of the Hungarian isolates (Supplementary Table S1). The overall trend of high levels of cox1 divergence was also seen in the Hungarian isolate with 41 sequences from the total 65 being unique haplotypes relative to the 21 sequences from 66 representing unique haplotypes within the ITS. This was also highlighted in the considerably lower haplotype diversity, nucleotide diversity and pairwise differences seen in the ITS relative to that of the cox1 (Supplementary Table S1).
Statistically parsimonious networks were constructed to establish the genealogical relationships between haplotypes of both cox1 and ITS. The cox1 haplotype network shows a clearly defined transition from Asia into Europe with distinct Chinese, Iranian and Hungarian clades (Fig. 1A). Surprisingly, all cox1 haplotypes identified were geographically specific and no haplotypes were shared between countries. The 41 Hungarian cox1 haplotypes also appeared to have become highly radiated with only two haplotypes sharing considerable similarity with those from Iran (Fig. 1A). This pattern of highly structured divergence was not shared in the ITS with no real defined clades due to only a single haplotype being specific to China and 19 haplotypes being specific to Hungary. The majority of haplotypes in Hungary were shared with those identified from Asia with 26 ITS sequences being identical to those from China and further 15 sequences being identical to those generated from Iran (Fig. 1B). The AMOVA for the cox1 gene also indicated that genetic differentiation taking place between the Chinese, Iranian and Hungarian populations showing significant differences in the frequency of pairwise differences between countries (cox1ΦST = 0·41874, P = 0·0001, P < 0·001). This was not reflected in the ITS as the AMOVA did not provide any significant differences between Hungary and the Asiatic isolates (ITSΦST = −0·04994, P = 0·524, P > 0·05).
Statistically most parsimonious networks inferred using the cox1 (A) and ITS (B) from Schistosoma turkestanicum. Topology between networks shows clear incongruence with the geographically distinct lineages seen in cox1 network (A) but lacking in differentiation in the ITS (B). The sizes of the circles are proportional to the number of samples that share that specific haplotype.
Interhost variation in Hungarian S. turkestanicum
In order to assess nucleotide and haplotype diversity within and across hosts in the Hungarian population of S. turkestanicum, nucleotide diversity of the cox1 and the ITS was measured between the worms taken from individual hosts. With the exception of worms sequenced from host 1, overall the cox1 displayed considerably higher nucleotide diversity within hosts than the ITS region, which was shown to be significant (average π overall cox1: 0·04328, ITS: 0·001952, cox1 n: 7, ITS n: 7, Wilcoxon's W: 1, P = 0·02, P < 0·05) (Fig. 2A). The majority of the cox1 haplotypes appeared to be unique to each host with only cox1 haplotypes 12, 14, 25, 32 and 33 being shared between hosts (Fig. 2B). This was contrasted in the ITS with ITS haplotype 2, which is shared with China, appearing in worm populations from all hosts as did ITS haplotype 3, which is shared with Iran; all other ITS haplotypes were host specific (Fig. 2C). The AMOVA analysis was applied by considering each host as an individual population of worms for both the cox1 and ITS. The analysis illustrated there to be genetic differentiation between the host populations (Interhost cox1ΦST = 0·51545, P = 0·0001, P < 0·001; Interhost ITSΦST = 0·148848, P = 0·0001, P < 0·001).
Interhost molecular diversity and genotype frequency of Schistosoma turkestancium from Hungary. Where (A) illustrates the nucleotide diversity of cox1 and ITS of S. turkestanicum from individual deer hosts; (B) illustrates the distribution and frequency of cox1 haplotypes between individual hosts; (C) illustrates the distribution and frequency of ITS haplotypes between individual hosts.
Demographic history and population expansion
The Tajima's D and the Fu's F S showed there to be significant deviation from neutrality in the combined cox1 alignment considering all the isolates as a single population. The Tajima's D was negative at –1·53725 and significant with P = 0·041, P < 0·05, a pattern, which also occurred in the Fu's F S at –42·2694 and highly significant at with P = 0·0000, P < 0·001. This suggests that either selection or historical demographic population processes could have given rise to an excess of low-frequency haplotypes as a result of past population expansion events. The pairwise mismatch analysis for the whole S. turkestanicum cox1 alignment appeared to be bimodal and the comparisons of the experimental distribution of mismatched pairs known as the raggedness index r significantly differed from observed data (r = 0·0046, P = 0·0000, P < 0·001). The bimodal nature of the mismatches along with significant raggedness of the data indicates that the populations are not only distinct but that they have not undergone any major population expansions in their recent histories (Fig. 3AI). This was also true when cox1 sequences from the population in Hungary were analysed without representation of S. turkestanicum from Asia. The Tajima's D was negative at −8·48501, and it was significant with P = 0·02, P < 0·05, and the Fu's F S was −29·3299 and highly significant with P = 0·0000, P < 0·001. The mismatch distribution analysis showed the same characteristics as described previously being significantly differing from the observed data illustrating the bimodal nature of the data and indicating a stable population (r = 0·0057, P = 0·0000, P < 0·001) (Fig. 3BI).
Mismatch distribution curves for cox1 and ITS haplotypes of all Schistosoma turkestanicum isolates (A) and the Hungarian isolate alone (B). Where values of expected frequency of pairwise mismatches under expanding population model are shown with a dashed line and the observed values are indicated with a solid line. In both A and B, the cox1 (I) appears bimodal and considered stable whereas the ITS (A and B II) appears to have undergone a recent expansion event.
However, this was not the case for the ITS alignment with both the full alignment and those sequences representing only Hungary appearing to be evolving neutrally and having gone through a recent population expansion event. The alignment representing all ITS sequences of S. turkestanicum from China, Iran and Hungary Tajima's D was negative at −1·47631 and not significant at P = 0·06, P > 0·05 unlike the Fu's F S that was negative −16·401 but significant at P = 0·0000, P < 0·001. Interestingly, the same pattern was seen in the sequences that represented the Hungarian population with Tajima's D being negative at −1·35104 not significant with P = 0·07, P < 0·05 and Fu's F S also being negative at −14·683 but significant at P = 0·0000, P < 0·001. Yet, in both full and Hungarian datasets, the mismatch distribution for the ITS was unimodal and raggedness index not being significant for both (full dataset r = 0·0482, P = 0·14615, P > 0·05; Hungarian dataset r = 0·0461, P = 0·14028, P > 0·05) (Fig. 3AII and BII).
Signatures of positive selection in cox1
Previous analysis illustrated a deviation from neutrality in the cox1 indicating the occurrence of selection acting upon the cox1 marker. Across the three geographical isolates, there is an indication of positive selection with ω > 1 with diversity associated with non-synonymous mutations being higher than the nucleotide diversity brought about by synonymous mutations. This was a feature that appeared in both the Iranian and Hungarian populations but was absent from the Chinese isolates with ω < 1 (Table 1). The sliding window analysis revealed that substantial contribution of non-synonymous mutations to nucleotide diversity across the cox1 fragment. The greatest peaks in diversity of non-synonymous mutations were seen between 287 and 587 bp with peaks ranging from 0·00872 to 0·20331 (Supplementary Fig. S1). Surprisingly when frequency of synonymous and non-synonymous mutations and their contribution to pairwise distance was assessed across all populations, the synonymous mutations appeared to saturate and plateau at approximately eight mutations resulting in pairwise distances from 0·01 upwards. This was contrasted by the non-synonymous mutations, which appeared to continue to accumulate and contribute to pairwise distances with up to 16 mutations contributing to a pairwise distance of 0·02163 (Supplementary Fig. S1).
Signature of selection in the cytochrome oxidase subunit 1 gene in populations of Schistosoma turkestanicum
Where Ka = divergence based on non-synonymous mutations, Ks = divergence based on synonymous mutations Ka/Ks (ω) indicates the presence or absence of selection, where ω > 1 represents positive selection, ω < 1 represents purifying selection and ω = 1 is indicative of neutrality.
Alterations in physiochemical properties of cox1 was assessed across the phylogeny of the S. turkestanicum phylogeny in TREESAAP 3.2 by using a sliding window analysis to compare the number of observed amino acid replacements relative to those under neutrality. Replacements affecting a total of 20 physiochemical properties were assessed across eight categories with 1–3 being considered to indicate only small effect on amino acid properties within a peptide, 4–5 having a medium effect, and 6–8 representing the most profound alterations in protein biochemistry due to replacement of amino acid (positive selection). TREESAAP employed a z-test to evaluate the significance of replacements and to identify the type of selection acting upon the bases within the window. Of the 20 physiochemical properties to be assessed, 16 appeared to be under some form of selection. However, only two properties had z-scores over ±3·09 with P < 0·001, which has been specified as a minimum score to reduce false-positive detection rate for the sites under positive selection. The properties indicating a radical amino acid replacement due to positive selection were isoelectric point (z = 7·656), which scored a category level 8 and power to be a C-terminal (z = 12·741) with a category level 6 (Table 2). This is illustrative of the phenotypic impacts of the non-synonymous mutations acting on the cox1 of S. turkestanicum.
Physicochemical property changes as a result of amino acid substitutions in the cytochrome c oxidase 1 subunit protein across populations of Schistosoma turkestanicum
Significant TREESAAP global negative (purifying selection) and positive (positive selection) z-scores for amino acid physicochemical property changes across the cox1 protein. Radical categories 5–8 shown with changes where P < 0·001 are most profound and likely, where 5 is of least impact and 8 is most profound change to properties and function of the protein.
P < 0·05*, P < 0·01**, P < 0·001***.
DISCUSSION
The consistent conflict in the cox1 and ITS evolutionary histories is indicative of mito-nuclear discordance in S. turkestanicum. As in Lawton and Majoros (Reference Lawton and Majoros2013), all cox1 sequences generated in the present study were unique to Hungary with the presence of intermediate forms being produced between the Iranian haplotypes (cox1 H6–H9) and those initially generated in the Lawton and Majoros (Reference Lawton and Majoros2013) study (cox1 H10 and H11). Conversely, only 19 ITS haplotypes were specific to Hungary with the majority being identical to either the Iranian ITS sequences (ITS H3) or those from China (ITS H2). This was further supported by the significant differentiation seen in the AMOVA analysis of the cox1 data, but not that of the ITS. Such patterns in discordance between mitochondrial and nuclear markers tend to arise either as a result of interspecies hybridization, populations rejoining after long periods of isolation or in situ discordance as a result of local adaptation (Toews and Brelsford, Reference Toews and Brelsford2012). Although not formally described, mito-nuclear discordance has been shown to occur several times in Schistosoma populations with natural hybrids between Schistosoma mansoni and Schistosoma rodhaini (Morgan et al. Reference Morgan, DeJong, Lwambo, Mungai, Mkoji and Locker2003; Steinauer et al. Reference Steinauer, Hanelt, Mwangi, Maina, Agola, Kinuthia, Mutuku, Mungai, Wilson, Mkoji and Loker2008) and with the identification of a complex introgressive process between species within the Schistosoma haematobium group (Huyse et al. Reference Huyse, Webster, Geldof, Stothard, Diaw, Polman and Rollinson2009; Webster et al. Reference Webster, Webster, Gouvras, Garba, Lamine, Diaw, Seye, Tchuem Tchuenté, Simoonga, Mubila, Mwanga, Lwambo, Kabatereine, Lange, Kariuki, Mkoji, Rollinson and Stothard2013). However, there is no evidence of interspecies hybridization between S. turkestanicum and another Schistosoma species in the present study or any other, and therefore it is suggested that discordance identified in this study arose through a combination of secondary contact of long isolated populations and in situ local adaptation (Paczesniak et al. Reference Paczesniak, Jokela, Larkin and Neiman2013; Pavlova et al. Reference Pavlova, Amos, Joseph, Loynes, Austin, Keogh, Stone, Nicholls and Sunnucks2013; Denton et al. Reference Denton, Kenyon, Greenwald and Lisle2014; Sun et al. Reference Sun, Wang, Zhang, Chapuis, Jiang, Hu, Yang, Ge, Xue and Hong2015). The cox1 haplotype network clearly illustrates the presence of a range of mitochondrial haplotypes giving rise to potentially new Hungarian lineages. This co-existence of divergent mitochondrial lineages and haplotypes in a single sympatric population is a signature of historic isolation and potential secondary contact (Zink and Barrowclough, 2008; Pavlova et al. Reference Pavlova, Amos, Joseph, Loynes, Austin, Keogh, Stone, Nicholls and Sunnucks2013). The pairwise mismatch analysis also indicated that the cox1 haplotypes are stable and have not undergone a recent expansion event contrasted by that of the ITS, which correlates with continuous historical gene flow, again a pattern, which has been described consistently in other long established populations of schistosomes and other fluke species (Morgan et al. Reference Morgan, DeJong, Lwambo, Mungai, Mkoji and Locker2003; Attwood et al. Reference Attwood, Fatih and Upatham2008; Poole et al. Reference Poole, Hellmann, Jensen and Nielsen2010; Zhao et al. Reference Zhao, Jiang, Dong and Nie2012; Brusentsov et al. Reference Brusentsov, Katokhin, Brusentsova, Shekhovtsov, Borovikov, Goncharenko, Lider, Romashov, Rusinek, Shibitov, Suleymanov, Yevtushenko and Mordvinov2013; Lawton and Majoros, Reference Lawton and Majoros2013; Webster et al. Reference Webster, Webster, Gouvras, Garba, Lamine, Diaw, Seye, Tchuem Tchuenté, Simoonga, Mubila, Mwanga, Lwambo, Kabatereine, Lange, Kariuki, Mkoji, Rollinson and Stothard2013; Goulding and Cohen, Reference Goulding and Cohen2014; Sun et al. Reference Sun, Wang, Zhang, Chapuis, Jiang, Hu, Yang, Ge, Xue and Hong2015).
The cox1 sequences of S. turkestanicum show the presence of positive selection within and between isolates. The TREESAAP analysis illustrated a distinct phenotypic change altering the properties of the resultant protein produced by the cox1 gene. Cheviron and Brumfield (Reference Cheviron and Brumfield2009) suggested that mitochondrial haplotypes could be differentially adapted to environmental conditions, which could account for the strong mitochondrial genetic structuring not seen in the nuclear markers, a process documented in a range of vertebrate and invertebrate groups (Paczesniak et al. Reference Paczesniak, Jokela, Larkin and Neiman2013; Pavlova et al. Reference Pavlova, Amos, Joseph, Loynes, Austin, Keogh, Stone, Nicholls and Sunnucks2013; Sun et al. Reference Sun, Wang, Zhang, Chapuis, Jiang, Hu, Yang, Ge, Xue and Hong2014; Tomasco and Lessa, Reference Tomasco and Lessa2011; McDonough et al. Reference McDonough, Šumbera, Mazoch, Ferguson, Phillips and Bryja2015; Morales et al. Reference Morales, Pavlova, Joseph and Sunnucks2015). The TREESAAP analysis revealed that the most important mutation affects isoelectric point of the cox1 protein altering protein folding or molecular charge and the efficiency of the electron transport system (Wikström and Verkhovsky, Reference Wikström and Verkhovsky2011). The cytochrome c oxidases are a multimeric complex, crucial for mitochondrial oxidative phosphorylation (OXPHOS), directly influencing metabolic performance (Wikström and Verkhovsky, Reference Wikström and Verkhovsky2011). The OXPHOS also provides heat to maintain body temperature through proton leakage (Sun et al. Reference Sun, Wang, Zhang, Chapuis, Jiang, Hu, Yang, Ge, Xue and Hong2014) and such adaptive changes have been shown in the cox1 under hypoxic conditions and as a response to tolerate cold temperatures in insects (Katewa and Ballard, Reference Katewa and Ballard2007) and thermal regulation European anchovy (Silva et al. Reference Silva, Lima, Martel and Castilho2014). Therefore, it could be suggested that the differential patterns in mitochondrial haplotypes shown in this present study may be related to the environment where S. turkestanicum infection foci have established, with the differences seen in the Hungarian isolates being an adaptation to the cooler summers and winters of central and Eastern Europe relative to those in Asia. This type of adaptation could allow the larval forms of S. turkestanicum to withstand greater thermal ranges in temperate aquatic systems. However, without sampling from a wide range of environments across Europe and Asia from a variety of hosts, it is not yet possible to assess in detail the role such factors could have on the metabolic evolution of S. turkestanicum.
The potential adaptive responses seen in the Hungarian populations could also be influenced by host as highlighted by Lawton and Majoros (Reference Lawton and Majoros2013). The Hungarian S. turkestanicum exists within a wildlife system unlike that of those sampled from Iran and China, which infect domestic animals. Wang et al. (Reference Wang, Chen, Zhao, Chen, Zhai, Li and Zhu2009) highlighted that S. turkestanicum has a wide host range parasitizing not only domestic cattle but also a wide range of other domestic and wild ungulates, rodents and lagomorphs, and divergence between subpopulations of schistosome species infecting different host species is well documented (Morgan et al. Reference Morgan, DeJong, Lwambo, Mungai, Mkoji and Locker2003; Hanelt et al. Reference Hanelt, Mwangi, Kinuthia, Maina, Agola, Mutuku, Steinauer, Agwanda, Kigo, Mungai, Loker and Mkoji2010; Lu et al. Reference Lu, Wang, Rudge, Donnelly, Fang and Webster2011). In the present study, each individual deer host was shown to harbour its own genetically distinct subset of parasites with some worm haplotypes being host specific for both the cox1 and the ITS. The AMOVA analysis also illustrated that S. turkestanicum did form distinct subpopulations between hosts, which would exist as unique interbreeding populations within a single host. Such patterns of genetic differentiation within parasite populations have been seen between individuals of the same host species within a given locality several times in both S. mansoni (Curtis and Minchella, Reference Curtis and Minchella2000; Sire et al. Reference Sire, Durand, Pointer and Theron2001; Curtis et al. Reference Curtis, Sorensen and Minchella2002; Prugnolle et al. Reference Prugnolle, Theron, Pointer, Jabbour-Zahab, Jarne, Durand and de Meeusm2005; Van den Broeck et al. Reference Van den Broeck, Meurs, Raeymaekers, Boon, Dieye, Volchaert, Polman and Huyse2014) and S. japonicum (Rudge et al. Reference Rudge, Lu, Fang, Wang, Basanez and Webster2009; Lu et al. Reference Lu, Wang, Rudge, Donnelly, Fang and Webster2011). Genetic diversity of parasites within a single host is thought to be predetermined by host mobility and their potential to visit several different infection foci, the age of the host and host immunity (Curtis and Minchella, Reference Curtis and Minchella2000; Curtis et al. Reference Curtis, Sorensen and Minchella2002; Lu et al. Reference Lu, Wang, Rudge, Donnelly, Fang and Webster2011; Van den Broeck et al. Reference Van den Broeck, Meurs, Raeymaekers, Boon, Dieye, Volchaert, Polman and Huyse2014). The high levels of genetic diversity within an individual host could relate to the number of transmission sites visited leading to the host becoming a ‘genetic mixing bowl’ for parasites, thus the offspring are genetically well mixed and parasite variation is randomly distributed between hosts (Curtis and Minchella, Reference Curtis and Minchella2000; Curtis et al. Reference Curtis, Sorensen and Minchella2002; Criscione et al. Reference Criscione, Poulin and Blouin2005). However, mixing and development of diverse infrapopulations is time dependant and Van den Broeck et al. (Reference Van den Broeck, Meurs, Raeymaekers, Boon, Dieye, Volchaert, Polman and Huyse2014) illustrated low genetic heterogeneity of S. mansoni infecting children in Senegal compared with the highly diverse populations found in adults. Age and movement of the host could account for the variation seen in the Hungarian S. turkestanicum as all adults were collected from livers of old large red deer stags, which are known to be wide roaming taking advantage of several wallowing sites all of which being potentially infective foci, along the Danube River and associated water bodies (Majoros et al. Reference Majoros, Dán and Erdélyi2010). It was also suggested that concomitant immunity could play a role, whereby an individual host develops partial immunity to a specific schistosome genotype. Inevitably, this would cause a reduction in the infection success of genetically similar individual schistosomes and increase the recruitment of genetically distinct schistosomes over the life time of the host (Van den Broeck et al. Reference Van den Broeck, Meurs, Raeymaekers, Boon, Dieye, Volchaert, Polman and Huyse2014). The authors of the current study argue that the high frequency of novel host-specific S. turkestnaicum haplotypes also coincides with the concomitant immunity hypothesis. Although it is clear that further work with increased sampling is required to confirm these hypotheses, the patterns in the data presented in the current study are typical of long-term well-established host–parasite interactions as seen in other schistosome host systems (Rudge et al. Reference Rudge, Lu, Fang, Wang, Basanez and Webster2009; Lu et al. Reference Lu, Wang, Rudge, Donnelly, Fang and Webster2011; Van den Broeck et al. Reference Van den Broeck, Meurs, Raeymaekers, Boon, Dieye, Volchaert, Polman and Huyse2014).
The current study further supports Lawton and Majoros (Reference Lawton and Majoros2013) original findings indicating the long establishment of S. turkestanicum in Hungary, with populations of the parasite appearing during the middle of the later part of the last interglacial period of the ice age and became locally adapted to the climatic conditions. Over the thousands of years that followed as deer, and probably other ungulates, moved back and forth from Europe to the Middle East and vice versa, this encouraged gene flow between S. turkestancium populations at nuclear loci. Due to the selective pressure of geographical region-specific metabolic requirements, novel mitochondrial genotypes were unable to establish in a temperate environment and were out-competed by Hungarian-specific mitogenomic adaptations. This would have allowed Asian-type nuclear genes to introgress into the Hungarian population but forced the extinction of Asian-type mitochondrial genotypes. Although in order to verify this hypothesis, extensive sampling of S. turkestancium from the Middle East, and Central and Southern Asia would be required. Nonetheless, such patterns of mito-nuclear discordance have been reported in populations of other ice age animals as a result of prolonged isolation followed by population secondary connection linked to glacial maxima and interglacial periods during the Pleistocene (Gompert et al. Reference Gompert, Forister, Fordyces and Nice2008; Sequeira et al. Reference Sequeira, Sodré, Ferrand, Bernardi, Sampaio, Schneider and Vallinoto2011; Kodandaramaiah et al. Reference Kodandaramaiah, Simonsen, Bromilow, Wahlberg and Sperling2013; Nesi et al. Reference Nesi, Kadjo, Pourrut, Leroy, Shongo, Cruaud and Hassanin2013; Paczesniak et al. Reference Paczesniak, Jokela, Larkin and Neiman2013; Pavlova et al. Reference Pavlova, Amos, Joseph, Loynes, Austin, Keogh, Stone, Nicholls and Sunnucks2013; Denton et al. Reference Denton, Kenyon, Greenwald and Lisle2014; McDonough et al. Reference McDonough, Šumbera, Mazoch, Ferguson, Phillips and Bryja2015; Morales et al. Reference Morales, Pavlova, Joseph and Sunnucks2015). However, the present study illustrates the complexities in disentangling the evolutionary history of an emerging parasite and shows the requirement for detailed understanding of the biological processes that govern variation in mitochondrial and nuclear markers. Such information could be used to better inform study design and illustrates the need for the further development of molecular markers that can be used to address specific epidemiological questions rather than the current dependence on generic markers often used in taxonomic studies. Despite these limitations, molecular studies such as the one presented here are crucial not only for providing initial insights into the biology of S. turkestanicum but also provide techniques and a baseline from which to monitor the spread of an emerging economically important parasite.
SUPPLEMENTARY MATERIAL
The supplementary material for this article can be found at https://doi.org/10.1017/S0031182017000920.
ACKNOWLEDGEMENTS
The authors would like to thank the Molecular Sequencing Facility at the Natural History Museum for performing the DNA sequencing reactions. The authors are also grateful to the Hungarian forestry commission, the rangers of the Gemenc Forest and the hunters for providing material for this study.
FINANCIAL SUPPORT
This research received no specific grant from any funding agency, commercial or not-for-profit sectors.
