Genetic and morphological variation of metacercariae of Microphallus piriformes (Trematoda, Microphallidae): Effects of paraxenia and geographic location.

Host organism offers an environment for a parasite, and this environment is heterogenous within the host, variable among individual as well as between the hosts, and changing during the host's lifetime. This heterogeneity may act as a prerequisite for parasite species divergence. Intraspecific variability related to a certain type of heterogeneity may indicate an initial stage of speciation, and thus poses an evolutionary importance. Here we analyzed genetic and morphologic variation of trematode metacercariae of Microphallus piriformes (Trematoda, Microphallidae). Genetic variability of trematodes was assessed from sequences of cytochrome c oxidase subunit 1 (COI) and internal transcribed spacer region (ITS-1). Morphological variation of metacercarial body shape was for the first time analyzed using geometric morphometrics. Parasites from the White Sea and the Barents Sea coasts demonstrated partial genetic divergence (according to COI sequence analysis) and had significantly different body shape. Neither genetic nor morphological variation of metacercariae was related to intermediate host species. We discuss possible causes of the observed genetic divergence of parasite populations in different geographic regions.

Introduction

Intraspecies variability (due to both genetic polymorphism and phenotypic plasticity) is a fundamental attribute of living organisms, making them able to survive under variable conditions within species ranges (Lewontin, 1957; Selander and Kaufman, 1973; Meyers and Bull, 2002). In turn, local adaptation is an evolutionary important process, acting as a prerequisite for ecological speciation (both allo- and sympatric) (Mayr, 1970; Thibert-Plante and Hendry, 2011; Lenormand, 2012; Nosil, 2012). Although the terminology of ecological speciation and related concepts were developed mainly for free-living organisms, all this reasoning is quite fair for parasites as well (Combes, 2001). Host organism acts as a first-order environment for parasite (Pavlovsky, 1934; Dogiel, 1964). This environment is heterogeneous within a host, variable among individual hosts, and changing during the host's lifetime.

Some parasites, being at the same ontogenetic stage, are able to infect several host species. Such taxonomically different hosts functionally similar for the life cycle of a parasite are called paraxenic, and the phenomenon is called paraxenia (Granovitch, 1996, 1999; 2009, 2016). Evidently, hosts belonging to different species provide different conditions for a parasite (due to their ecological and physiological dissimilarities) contributing to the general variability of parasites' environment. Adaptation to a certain host species represents the parasitic version of local adaptation. This might become a prerequisite for alloxenic speciation (as adaptation to hosts of different species), a parasitic version of ecological speciation (Combes, 2001). It is fair not exclusively for definitive hosts (where sexual reproduction occurs) but for intermediate hosts as well. For example, association with certain intermediate host species may facilitate transmission to certain definitive host species, prompt partially independent parasite's life cycle functioning of these subsets and finally drive their divergence. Paraxenia at the level of intermediate host may cause genetic subdivision of parasite population also due to differential survival of parasites in different host species, which is crucial for supporting of parasites genetic diversity and essential for the whole parasitic system stability (Granovitch, 1999; Granovitch and Mikhailova, 2004). Importantly, first-order environment for a parasite is not continuous but “patchy”, being comprised of many individual hosts. Owing to this, a population of a parasite (all the individuals at all ontogenetic stages within some geographical region, Combes, 2001) is hence subdivided by (a) host individuals (infrapopulations, Margolis et al., 1982),(b) host species harboring the same parasitic stages (“paraxenic hemipopulations” sensu Beklemishev, 1960; Granovitch, 1996, or “xenopopulations” sensu Combes, 2001); (c) host species harboring different ontogenetic stages of a parasite due to complexity of parasite life-cycle (“metaxenic hemipopulations”, Granovitch, 1996).

On the other hand, geographically distant host populations can strongly differentiate due to either local adaptation or neutral genetic shift or both; they will not be evenly matched as an environment for a parasite. Such a host “non-equivalence” within a geographic range is expected to affect polymorphism in a parasite population and might contribute to divergence processes in case of restricted gene flow between parasite populations. Thus, polymorphism and the discrete genetic structure in a parasite population may be related to heterogeneity of both the first-order (set of paraxenic hosts, their sex, age and immune competence etc.; e.g. Granovitch et al., 1987; Granovitch and Gorbushin, 1995; Lively et al., 2004; Koehler and Poulin, 2012; Levakin et al., 2013) and second-order environment (distant geographic regions, contrasting biotopes etc.; e.g. Blair et al., 1997; Semyenova et al., 2006; Sithithaworn et al., 2007; Webster et al., 2007; Krapivin et al., 2018). Study of polymorphism and population genetic structure due to influence of the first and second-order environment can bring new insight in parasitic systems functioning in nature and evolutionary processes occurring there.

The trematode Microphallus piriformes (Odhner, 1905) (Plagiorchiida, Microphallidae) from the North Atlantic intertidal has a dixenic life cycle. Five periwinkle species of the Littorina Ferussac, 1822 genus (L. saxatilis (Olivi, 1792), L. arcana Hannaford Ellis, 1978, L. compressa Jeffreys, 1865 – “saxatilis” group; L. obtusata (Linnaeus, 1758), L. fabalis Turton, 1825 – “obtusata” group) serve as paraxenic intermediate hosts; several seabird species (gulls, e.g. Larus argentatus Pontoppidan, 1763, eiders, most commonly Somateria mollissima (Linnaeus, 1758), oystercatchers, e.g. Haematopus ostralegus Linnaeus, 1758) are definitive hosts (Galaktionov et al., 2012). M. piriformes, like other species from “pygmeus” group, does not have a free-living stage: metacercariae are formed within daughter sporocysts. Metacercariae within mature daughter sporocysts look like adult worms which have not started egg production yet. They have almost completely developed somatic organs, but do not possess a mature reproductive system. This accelerates maturation and potentiates reproduction efficacy in case of definitive host infection (Galaktionov, 1993; Galaktionov and Dobrovolskij, 2003). Infective metacercariae are remarkably variable in their body size and shape among hosts individuals (our own observations). In addition, this trematode species has a high level of genetic diversity (Galaktionov et al., 2004, 2005, 2008). We consider M. piriformes as a convenient model to study genetic and morphological variation related to paraxenia and geographical distribution, since this parasite is widely distributed along the coasts of the Northern Europe and infects several intermediate hosts.

Material and methods

Sampling

Sporocysts with metacercariae of Microphallus piriformes were sampled during summer expeditions of 2017 and 2018 from wild populations at 4 geographic locations (Fig. 1; Table 1). In total 72 samples were collected.

Fig. 1

The map of the study region (image: TerraMetrics, map data: Google). Sample collection sites (Tromsø city, Kiberg settlement, Dalnie Zelentsy settlement, Sredny Island) are shown.

Table 1

Information on sampling sites. Only M. pygmeaus which was used for genetic analysis as an outgroup was collected in the Yakovleva inlet.

Collection sites	Water area	Year	Coordinates
Korga-Islet in the Levina inlet, vicinities of Sredniy Island	Kandakaksha Bay,White Sea	2018	66°18′04.09″N33°27′27.02″E
Yakovleva inlet, vicinities of Sredniy Island	Kandakaksha Bay,White Sea	2018	66°18′50.01″N33°50′20.06″E
Vicinities of Dalnie Zeletsy	Barents Sea	2017	69°06′56″N36°04′10″E
Vicinities of Kiberg	Varanger Fjord, Barents Sea	2017	70°17′07″N30°59′53″E
Vicinities of Tromsø	Norwegian Sea	2018	69°40′58″N18°56′34″E

Identification and photography

Adult individuals of littorinid snails L. arcana, L. compressa, L. fabalis, L. obtusata and L. saxatilis were collected from intertidal area and transported to the laboratory. For sample preparation snails were dissected under MBI-10 binocular microscope to identify species, sex and possible trematode infection. Infected hepatic tissues were separated and soaked in distilled water, where parasites were carefully washed out from host tissues using glass pipettes. Washed parasites were used for (1) temporary preparations for identification and photographing; (2) fixation in 96% ethanol, final DNA isolation and molecular analysis.

To complete relax and flattening of metacercariae on temporary preparation, glass slides were heated in a drop of distilled water for 2 min at 70 °C using heating table (Galaktionov, 1980, 1993) (Supplementary Fig. 1). Photographic pictures were taken with MBI-10 binocular microscope coupled with MFU photo adapter and Canon EOS 1200D camera; camera settings remained unaltered throughout the study.

Molecular markers analysis

DNA isolation was performed using minispin collumns for genomic DNA extraction (ExtractDNA Blood, Evrogen) following recommendations of the manufacturer (http://evrogen.ru/kit-user-anuals/ExtractDNA_Blood.pdf). DNA samples were stored at −20 °C until use.

Two molecular markers: fragments of the mitochondrial cytochrome oxidase subunit I (COI) and internal transcribed spacer (ITS-1) were used for the analysis. The primer sequences used for amplification are in Table 2.

Table 2

Sequences of primers used for amplification of molecular markers.

Direction	Marker	Sequence 5’ -> 3′	Fragment length	Reference
Fw	COI	TTTTTTGGGCATCCTGAGGTTTAT	368 bp	Bowles et al. (1992)
Rev	COI	TAAAGAAAGAACATAATGAAAATG		Bowles et al. (1992)
Fw	ITS-1	ACACCGCCCGTCGCTACTA	347 bp	Galaktionov et al. (2012)
Rev	ITS-1	TGGACGAAACTGCGCGCTTC		Authors

Ready PCR mixture (iQТМ SYBR® Green Supermix, Bio-Rad) was used for amplification following manufacturer's recommendations. PCR was performed using Veriti 96-Well Thermal Cycler (ThermoFisher Scientific). Amplification method was initial denaturation at 95 °C for 3 min, in cycle denaturation at 95 °C for 30 s, primers annealing at 56 °C for 30 seс, in cycle elongation at 72 °C for 30 seс, final elongation at 72 °C for 4 min; 30 cycles.

Presence of PCR products of the expected length was assessed by 1%-agarose gel electrophoresis in Tris-acetate-EDTA buffer, pH 8.3 (TAE, Evrogen); SYBR green (50x SYBR® Green I for PCR, Evrogen) was added to samples while in-gel loading for visualization; electropherogram was visualized with BioRad ChemiDoc MP at wavelength 497 nm.

Amplicons were sequenced by Sanger in both directions using 3500xL Genetic Analyzer (ThermoFisher Scientific) and BigDyeТМ Terminator v3.1 Cycle Sequencing Kit (ThermoFisher Scientific) with technical support of Alexey Masharsky. Sequencing results were analyzed using ChromasPro v.1.7.4 software (Technelysium Pty Ltd); sequence alignment was performed with MEGA X v.10.0.4 software (Kumar et al., 2018).

Phylogenetic analysis

We used Bayesian inference for phylogenetic reconstruction in MrBayes v.3.2 software (Ronquist et al., 2012). Nucleotide substitution model was GTR + G + I (Lanave et al., 1984). Optimal nucleotide substitution model was assorted using PAUP4 (Swofford, 2002) and MrModeltest v.2.3 (Nylander, 2004) software packages. Analysis was performed as two independent runs, five chains in each (four heated and one cold; the first 25% samples from the cold chain were discarded) for 15 000 000 generations with a sample frequency of 1000, print frequency of 1000 and diagnostics calculated every 1000 generations. The convergence between two runs was tested by comparison of statistical parameters in the Tracer Software (Rambaut et al., 2018). 65 COI sequences of 72 were used for analysis, as we failed to get sequences of acceptable quality for other 7 samples. We also used sequences of Microphallus pygmaeus (Levinsen, 1881), Microphallus triangulatus Galaktionov, 1984 and Microphallus similis (Jägerskiöld, 1900) as outgroups in our analysis. Results of phylogenetic analysis were visualized using FigTree v.1.4.3 software (Rambaut, 2009).

M. piriformes haplotypes network was constructed according to TCS algorithm (Clement et al., 2000) using PopART software (Leigh and Bryant, 2015).

Geometric morphometrics

Tps-software package (Rohlf, 2004; The Stony Brook Morphometrics) was used for geometric morphometrics analysis. Two-dimensional images were compiled to tps format using tpsUtil v.1.74, the coordinate locations of landmarks and semilandmarks were digitised in tpsDig2 v.2.3. We used nine landmarks to model metacercarial body shape: two along middle body line at anterior and posterior poles, two at edges of most wide part of posterior body compartment, two at edges of most narrow body part (at the transition between anterior locomotory and posterior generative body parts), one at the gut bifurcation, and two on the body edges at the level of the gut bifurcation. Additionally, 22 semilandmarks forming four curves were used to describe the contour of the body (Fig. 2). Images of 46 clones (clone = metacercariae from different daughter sporocysts within 1 host individual) were analyzed (10 metacercariae from each clone).

Fig. 2

Scheme of landmarks and semilandmarks on the metacercarial body.

To minimize the effect of allometric variability only fully formed metacercariae with all anatomical details visible were used in the analysis. These metacercariae are characterized by completely formed type II marita glands (the final stage of maturation, no morphogenic changes are known after this stage up to maritogony process starting in the final host; Galaktionov, 1991, 1993).

Statistical analysis

A geometric morphometric analysis was carried out in R (R v.3.5.3; R Core Team, 2019) using geomorph (Adams et al., 2019) and abind (Plate, Heiberger, 2016) packages. Shape variables were derived from a Generalized Procrustes Analysis (GPA; Gower, 1975), which superimposes and aligns landmark configurations (through translation, rotation and isometric scaling) removing the information not related to shape. Aligned procrustes coordinates were subjected to Principal Component Analysis (PCA) to visualize shape variation of metacercariae. Proportion of variance explained by each principal component was computed from eigenvalues.

Non-parametric Multivariate Analysis of Variance (perMANOVA; Anderson, 2001) on the matrix of procrustes coordinates was used (1) to compare the shape of metacercariae from L. saxatilis and L. obtusata from the White and Barents Seas; and (2) to test for possible allometric relationships of shape and size of metacercariae. The design included the effect of the host individual nested within a particular species to account for non-independence of the observations belonging to the same clone of metacercariae. We tested significance of shape variations using randomized residual permutation procedure (RRPP; Collyer and Adams, 2018, Collyer and Adams, 2019; RRPP package v.0.4.2). Post-hoc pairwise comparisons were held using perMANOVA; p-values were corrected using Holm-Bonferroni procedure (Holm, 1979). Morphological disparity was assessed using morphol.disparity () function from the geomorph package for groups of meracercariae from different host species and sampling locations. Relative morphological disparity was computed as a percentage of the sum of within-group disparities.

To compare centroid size of metacercariae from different hosts and sites, we applied a General Linear Mixed Model with Gaussian error distribution (GLMM; Pinheiro, Bates, 2000) using lme4 package (Bates et al., 2015). The model included one fixed predictor (group of metacercariae defined by host species and site) and one random effect (host individual). The model assumptions were checked visually on the plots of residuals. Effect significance was tested with Likelihood Ratio Test (LRT). Intraclass Correlation Coefficient (ICC) was used to assess the size variability of metacercariae from an individual snail. Tukey's test was used for post hoc comparisons (Quinn, Keough, 2002).

The results of analysis were visualized using ggplot2 (Wickham, 2016) and cowplot (Wilke, 2019) packages. Tables were produced using knitr (Xie, 2019) and huxtable (Hugh-Jones, 2019) packages.

Results

Genetic polymorphism in populations of M. piriformes. ITS-1

All analyzed ITS-1 sequences were identical. This supports, on one hand, the unity of M. piriformes as a species and, on the other, it justifies the affiliation of all the samples to the same species. P-distance between sequences of M. piriformes and closely related species M. pygmaeus was 3.75% (12 SNP, 1 indel).

Genetic polymorphism in populations of M. piriformes. COI

COI gene proved to be rather polymorphic in M. piriformes: among 64 samples sequenced we have found 29 different haplotypes of COI fragment; three were observed most often. The Bayesian inference based on COI-sequences (Fig. 3) implies a limited gene flow between populations of the White and Barents Seas, as there are several location-specific clades on the tree, e.g., clades including samples exclusively from ether Korga-islet (White Sea) or Kiberg (Barents Sea).

Fig. 3

Bayesian inference based on COI sequence (369 bp); 15000000 generation; GTR + I + G substitution model; A posteriori probabilities are indicated by node shapes; sample name includes parasite species (pir – M. piriformes, pyg – M. pygmaeus, tri – M. triangulatus, sim – M. similis), sample number, geographic region and location (WSk – White Sea, Korga-Islet; WSy – White Sea, Yakovleva; DZe – Barents Sea, Dalnie Zelentsy; Kib - Barents Sea, Kiberg; Tro - Norwegian Sea, Tromsø), host species (sax – L. saxatilis; arc – L. arcana; comp – L. compressa; obt – L. obtusata; fab – L. fabalis). Identical haplotypes and the FST-value of differentiation between populations (Weir, and Cockerham, 1984) are shown in Supplementary Table 1. Branch color reflects geographic region.

Strong effect of geographic location on haplotypes distribution is obvious also in the haplotypes network (Fig. 4A): although most often haplotypes were detected in all three examined regions, there were also many location-specific haplotypes, including several singletons, clustering according geographic region. As evidenced by both the Bayesian inference (Fig. 3) and haplotypes network (Fig. 4B) there was no detectable population subdivision based on host species: haplotypes of M. piriformes parasitizing hosts of different species were interspersed after both analyses.

Fig. 4

Haplotype networks, COI sequence (369 bp); TCS algorithm; dashes correspond to mutations. A: color reflects sampling location. B: color reflects host species.

Intrapopulational polymorphism of metacercarial body shape

Parasite samples from L. saxatilis and L. obtusata hosts only, collected at the White and Barents (Kiberg and Dalnie Zelentsy) Sea shores, were used for analysis: there were five groups in comparison (L. saxatilis was the only host in Dalnie Zelentsy populations). Metacercarial body shape variation related to sampling site and host species is visible on the PCA ordination plot (Fig. 5); PerMANOVA accounting for intraclonal variability confirmed the effects of the factors “location” and “host species” (p < 0.01) on metacercarial body shape. Significant differences were found between all combinations of locations (except samples from Dalnie Zelentsy vs Kiberg, L. saxatilis as a host), but there were no differences in body shape between metacercariae from different hosts within the same location; this result is also visible on the general PCA-ordination plot (Fig. 6, Supplementary Table 2).

Fig. 5

PCA-ordination of individual PC1 can be interpreted as a deepness of a “waist” between locomotory and generative body parts; PC2 can be interpreted as a width of locomotory body part. B: Pairwise post-hoc comparison; significant value are shown as bold (considering Holmes correction for multiple comparison); host species: sax – L. saxatilis; obt – L. obtusata; sampling site: Kib - Barents Sea, Kiberg; Kor – White Sea, Korga-Islet; Zel – Barents Sea, Dalnie Zelentsy.

Fig. 6

A: PCA-ordination of individual Ellipses show 95%-confidence intervals for the five groups in comparison. B: Mean body shapes of . Transformation grids reflect deviation from the overall mean body shape.

within the same hemipopulation (hemipopulation is the totality of parasite individuals at the same stage of the life cycle within all host individuals within a given area (Beklemishev, 1960)) and within the same infrapopulation at different sampling sites (separately for different host species) was also assessed. This analysis revealed that metacercarial body shape varied most strongly within the White Sea hemipopulation; this was fair for the both host species. The same tendency was observed in case of infrapopulations. The metacercarial body shape variability was comparatively low within the Barents Sea hemi- and infrapopulations (both Kiberg and Dalnie Zelentsy) (Fig. 7, Supplementary Tables 3 and 4).

Fig. 7

Variability of metacercarial body shape within hemipopulations and infrapopulations of A: Absolute and relative morphological disparity (MD) of metacercariae within hosts of the same species. B: Distribution of morphological disparity (MD) within individual snails grouped by host species and sampling location.

Interactions of metacercarial body shape and size

Theoretically, differences in body shape might originate from allometry – i.e. age-related changes in proportions of the whole body and particular organs during ontogenesis (Klingenberg, 1996, 2016). To avoid allometry effects, only metacercariae with completely formed anatomical structures (see methods) were included into analysis. Nevertheless, the question of changes in body shape proportions after maturation (when metacercariae are ready to infect a final host and when transition to adult worm occurs) is still vague.

With respect to this, we evaluated interactions between body shape and size (using centroid size as a proxy) in M. piriformes metacercariae and found no significant allometric effect (perMANOVA). However, metacercarial body size varied significantly between the hosts species and locations (LRT = 39.5, df = 4, p < 0.01). Notably, the size of metacercariae from different hosts within the same location did not significantly differ (Fig. 8, Supplementary Table 5). High intraclass correlation coefficient (ICC = 0.71) implies that metacercariae of the same infrapopulation are of very similar size; this fits well with data on synchrony of microphallid development within daughter sporocysts (Galaktionov, 1993).

Fig. 8

Body size of Mean centroid size and 95% confidence interval obtained via bootstrap.

Discussion

In this study we evaluated genetic and morphological intraspecific variability in metacercariae of Microphallus piriformes at several levels of variability: between populations (effects of different geographic locations), between paraxenic hemipopulations (effects of intermediate host species), between clones (effect of a host individual), and between individuals (intraclonal variability).

Geography-related genetic polymorphism within M. piriformes populations

Our results indicate some gene flow limitations between populations of different geographic regions (strong between the White and Barents Sea populations; moderate between the two populations from the Barents Sea). Both intermediate (periwinkles) and final (birds) hosts of M. piriformes are widely distributed along the Northern Atlantic sea coasts. Hypothetically, gene flow between spatially disconnected populations might be maintained by bird migrations. Consistently, isolates of M. piriformes from the Vaygach Island (Barents Sea) and Iceland coast (straight distance exceeds 3000 km) demonstrated almost identical UP-PCR fingerprinting patterns (Galaktionov et al., 2008), while UP-PCR fingerprinting did not reveal any geography-related subdivision in any other species of the “pygmaeus” group (Galaktionov et al., 2005). This might be a consequence of a bird migration-maintained gene flow. Shown in our study prevalence of limited number of haplotypes shared by all studied populations agrees well with such a scenario. On the other hand, some heterogeneity of M. piriformes populations related to sampling site was also recorded in previous studies. In particular, genetic differences between several populations of M. piriformes from the White Sea were detected by RAPD approach (Khalturin et al., 2000). Similarly, genetic divergence between populations from the White and Barents seas (but not between the White and Norwegian seas) was demonstrated with UP-PCR approach (Galaktionov et al., 2004). Thus, our results agree well with the earlier data. Some irregularity in the patterns of population genetic differentiation, in our opinion, might be explained by migratory activity of final host birds funneled to specific routes. Gulls, eiders and waders regularly migrate along the North Atlantics shore (Bakken et al., 2003; Bianki, 1967, Bianki, 1989; Coulson et al., 1984; Dementiev and Vuchetich, 1947; Haftorn, 1971; Kokhanov and Skokova, 1967; Krasnov et al., 2015; Mineyev, 2003; Strom, 2006; Tatarinkova, 1970; Anker-Nilssen et al., 2000; Noskov et al., 2016). The schemes of bird migratory routes, reconstructed from published data, are presented in the Supplementary Fig. 2. Obviously, not all the populations of M. piriformes can be connected via bird migration, but only those lying along migratory routes. Locations of bird settlement sites are another factor, exacerbating genetic breach between parasite populations: wintering sites at the coasts of the White and Barents seas were recorded for gulls and eiders (Koryakin and Kondratiev, 1983; Kohanov and Shkliarevitch, 1985; Shkliarevitch, 1979; Сrаmр and Simmons, 1983; Scott and Rose, 1996). Final host settlement provides an opportunity for a year-round autonomous sustainability of local parasite populations, and as a consequence, their divergence from each other. Short adult longevity (no longer than 10 days) also strongly limits transmission capacity via bird migrations (Belopol'skaya, 1983; Galaktionov, 1993; Galaktionov and Dobrovolskij, 2003). Thus, spatially limited and transient connection of trematode populations by hosts migrations cause partial differentiation of the former along the geographical range. It is worth noting that the formation of such partial subdivision in the genetic structure of a parasite species can act as a prerequisite for allopatric speciation (Mayr, 1970; Huyse et al., 2005).

Evidence of partial isolation between different geographic regions was detected by the mitochondrial marker (COI) only, not the nuclear one (ITS-1). Mitochondrial genomes are expected to be more strongly affected by the factors discussed below. Although adaptive neutrality and clonality of mitochondrial genomes were questioned during the last decade, these are still expected to evolve faster than nuclear ones due to haploidy and maternal germline bottleneck, and thus are better indicators of shallow evolutionary events (rev. in Zink and Barrowclough, 2008; Galtier et al., 2009). Mitochondrial marker COI is widely used (often in combination with other molecular markers for more reliable results) in studies of genetic variation and genetic structure of trematode groups (Blasco-Costa et al., 2016). Often, these studies reveal the presence of several genetic lineages and a high haplotype diversity in fluke populations. Sometimes, genetic structure of a parasite species is consistent with the geographical distance between its populations (e.g. Blair et al., 1997; Semyenova et al., 2006; Stothard et al., 2009). Nevertheless, in other cases, such clear dependence can't be reliably identified (e.g. Lopatkin et al., 2010; Tatonova et al., 2013) for different reasons (influence of autogenetic processes, hosts migrations, insufficient sampling, etc.). In many studies, including ours, the authors conclude that genetic structure of trematode species and genetic lineages they form are significantly affected by features of parasites biology and life cycle and movements of their hosts (especially definitive) (Criscione and Blouin, 2004; Huyse et al., 2005; Keeney et al., 2009; Blasco-Costa et al., 2011).

was not observed in this study. Importantly, 3 species (Littorina arcana, L. compressa, L. fabalis) were relatively rarely detected as hosts of M. piriformes due to their low densities. Therefore, no informative conclusions on genetic specialization of the parasites to these hosts may be done. Most of data were gathered for parasites obtained from L. saxatilis and L. obtusata hosts, and will be discussed in detail below.

All the five species of paraxenic hosts of M. piriformes belong to the subgenus Littorina (Neritrema Recluz, 1869) (Littorinidae, Caenogastropoda). They are phylogenetically closely related and comparatively recently diverged from each other (estimated divergence time does not exceed 5 Ma) (Reid, 1996, 2012). These species differ in several sets of morphological, physiological and ecological characters. Thus, geographic ranges of Littorina arcana and L. compressa are narrower and more patchy compared to other species; they were not recorded from several regions where L. saxatilis, L. fabalis and L. obtusata are present, e.g. the White and North Seas, the Iberian Peninsula coast, etc. (Reid, 1996). There are also differences in the microbiotopic distribution, e.g. L. compressa, L. fabalis and L. obtusata inhabit the fucoid area only, while L. saxatilis and L. arcana can live in the upper part of the intertidal where fucoids are absent; L. fabalis occupies lower intertidal and upper subtidal area and more often is associated with Fucus serratus, in contrast to L. obtusata usually preferring F. vesiculosus and Ascophyllum nodosum (Reid, 1996; Sergievsky et al., 1997; Iastchenko and Granovitch, 2002; Granovitch et al., 2004, 2013). Littorina spp. also exhibit functional differences in tolerance to low salinity, desiccation and other stressors, as well as physiological plasticity (e.g. Todd, 1964; Clarke et al., 2000; Sokolova et al., 2000; Maltseva et al., 2016; Muraeva et al., 2017; Maltseva et al., 2020). Altogether, these differences might form preconditions for host-related speciation among paraxenic hemipopulations of parasites. It is known that subdivision of trematode populations at the genetic level may emerge in connection with paraxenia after expansion of host range via either new host species invasion (usually phylogenetically closely related) or cospeciation (e.g. LoVerde et al., 1985; Theron and Combes, 1988; Semyenova et al., 1995; Blair et al., 1997; Galaktionov et al., 2008, Koehler et al., 2011).

Nevertheless, the most probable divergence scenario of the microphallids from the “pygmaeus” species group assumes final host species and geographic isolation as drivers, and not the intermediate host species (Galaktionov et al., 2012). Despite a number of differences between the species of intermediate hosts, the possibility of host species-driven evolution of M. piriformes seems unlikely for several reasons:

The common ancestor of M. piriformes and other “pygmaeus” species expanded to the North Atlantic from the North Pacific between 3.5 Myr and 2.4 Myr BP (together with the Littorina (Neritrema) species and gull hosts) and infected several available Neritrema species (ancestors of the “obtusata” and “saxatilis” groups had already been diverged by then, Reid, 1996; Reid et al., 2012). The northern communication between Arctic-Atlantic and Pacific via Bering Strait closed again, and this promoted speciation of M. piriformes in the North Atlantic and M. calidris in the North Pacific; the branch leading to other “pygmaeus” species specified due to switch to seaducks species as final hosts (Galaktionov et al., 2012). Speciation of L. compressa, L. fabalis and L. arcana occurred in the North Atlantic after the Bering Landbridge emergence (Reid, 1996), and this did not cause partitioning of M. piriformes, which is able to successfully infect all those species. Accordingly, no evidence of the intermediate host-related subdivision in modern populations of M. piriformes was revealed in this study, while there was an evidence for geography-related subdivision.

It is also worth mentioning, that some Neritrema species are very similar at the both genetic and functional levels, e.g. L. saxatilis and L. arcana with indirect evidence of regular hybridization (Mikhailova et al., 2008, 2009; Granovitch et al., 2013; Maltseva et al., 2016, 2020), therefore hardly providing any background for parasite differentiation.

Additionally, some facts suggest that this trematode species is rather sensitive to host's physiological properties, which in turn is a prerequisite for host-driven speciation. In comparison to other species from “pygmaeus” group, M. piriformes is characterized by relatively narrow range of intermediate hosts: it can infect only snails from Littorina (Neritrema) subgenera. In contrast, e.g. M. pygmaeus can also heavily infect L. littorea (Linnaeus, 1758) (Werding, 1969; Pohley, 1976; Lauckner, 1984), while M. pseudopygmaeus Galaktionov, 2009 is known for its ability to infect a set of hosts of several families and orders (Galaktionov, 2009; Galaktionov et al., 2012). Such narrow specialization of M. piriformes reduces the likelihood of the formation of subdivision in connection with the species of intermediate host.

Affirming this, previous studies detected paraxenia-related genetic variability of M. piriformes. Trematodes from either L. saxatilis or L. obtusata (the hosts significantly different at both ecological and functional levels) collected from the same geographic location demonstrated distinct RAPD patterns depending on host species (Khalturin et al., 2000). Our findings are not necessarily inconsistent with this data due to the lesser sensitivity of our approach (only two markers were used in this study, so our conclusion that there is no intermediate host species related subdivision may not be approved by further whole-genome scans). Instead, together these results suggest that geography-related differences are more profound than those related to the intermediate host species.

Metacercarial body shape variability

Data on metacercariae (and also other ontogenetic stages) body shape variability in M. piriformes are scarce. Significant variability of the metacercarial body shape was shown for M. pseudopygmaeus which is closely phylogenetically related to M. piriformes, yet, there was no clear relations to different host species (represented by a rather diverse set of L. saxatilis, Onoba aculeus, Margarites helicinus, Solariella varicosa, Epheria vincta, Lacuna neritoides, Cryptonatica clausa) with exception of metacercariae from M. groenlandicus, which were relatively distinct in their body shape. Noteworthy, the same study revealed that the body shape of metacercariae from M. groenlandicus tended to correlate with geographic location (Dalnezelenetskaya bay vs Kolguev Island; Galaktionov, 1993). Still, the author stressed that the variability revealed was not constrained to discrete variants but rather continuous with all kinds of intermediate states. This is congruent with our results: there was a metacercarial body shape continuum, and those from either different host species or geographic locations tended to form certain variants more often with no break between these variants (Fig. 5, Fig. 6). The physiological and genetic background of this body shape variability is completely unclear.

Quite interestingly, not only the variability in body shape per se was revealed in our study, but also the differences between geographic locations in the degree of morphological disparity: it was maximal within the White Sea hemi- and infrapopulations. Our genetic data do not allow to suspect a higher genetic diversity of the White sea population (compared to the other two from the Barents Sea) as a possible explanation for an increased level of inter-hemipopulational variability. The higher degree of intraclonal shape variability (based on identical or very similar genetic background) might be interpreted similar to fluctuating asymmetry. Generally, there are two fundamental properties of any organism development restricting its phenotypic variation: canalization (among individuals) and developmental stability (within individuals) (Waddington, 1942; Schmalhausen, 1949; Debat and David, 2001; Willmore et al., 2007). It is well known that under pressure of environmental stressors the accuracy of the molecular and physiological background for both canalization and developmental stability may be compromised, which in turn, would cause an increase in phenotypic variation within and among individuals (e.g. Scharloo, 1991; Parsons, 1997; Hoffmann and Hercus, 2000; Lazic et al., 2015). The increase in both types of phenotypic variation was detected in the White Sea population which might indicate some kind of disturbance. Importantly, high level of phenotypic variation was accompanied with decreasing of the average body size, which might be due to the same reason. The organism size is a function of the growth rate and duration, which depend on temperature, nutrition and oxygen availability in ectotherms (rev. in Shingleton, 2011). The White Sea has a relatively low salinity (24–28 ppm), high seasonal salinity variation in the surface layer, long ice cover period (up to 8 month) (Naumov and Fediakov, 1993). These conditions are likely not optimal for periwinkles, the M. piriformes intermediate hosts Littorina saxatilis and L. obtusata (L. compressa and L. arcana are absent in this area), as evidenced by the sizes of these snails from the Barents and White Sea populations (Granovitch et al., 2004). Temperature together with oxygen availability (which solubility in water depends on temperature) form a latitudinal gradient shaping the size increase at high latitudes demonstrated in different groups of marine invertebrates, e.g. benthic amphipod crustaceans (Chapelle and Peck, 1999), benthic gastropods of the Turridae family (McClain and Rex, 2001) and oceanic nematodes (Soetaert et al., 2002). Body size reducing effects of hypoxia on growth were also described in not aquatic insects, e.g. Drosophila (rev. in Shingleton, 2011). The size differences in the Littorina snails between the Barents and White Sea populations are fully in line with such reasoning, as the White Sea population is the most southern among the studied. Suboptimal functioning of hosts at the White Sea is expected to affect parasite, since it raises basal metabolic costs (rev. in Sokolova, 2013). M. piriformes metacercarial body size forms a trend collinear with latitude increase (Fig. 8). This trend also corresponds well to the hosts’ size distribution on one hand, and to degree of metacercarial shape disparity within hemi- and infrapopulations on the other. Altogether, our results imply that fundamental environment-development interaction patterns known for free-living organisms are also legitimate for parasites, and that phenotypic variation in a parasite may act as an indicator of host functioning.

Conclusion

In our study we applied both genetic and morphometric approach to evaluate the degree of populations subdivision of the dixenic trematode M. piriformes. Both methods showed consistent results: there were signs of gene flow restriction between parasite populations from remote locations where genetic exchange is limited due to both absence of regular bird migration and short parasite lifetime within the bird. Additionally, a strongly increased level of intra- and interclonal metacerarian body shape disparity in the White Sea region together with a reduced mean body size were revealed using geometric morphometrics. We suggest that this may be result of a restricted nutrient supply, development destabilization and decanalization emerging from suboptimal environmental conditions for the intermediate hosts in that region. We failed to find any differences related to the intermediate host species, which contrasts with some earlier data (Khalturin et al., 2000). This implies, that inter-region differences are more pronounced than those between paraxenic hemipopulations; more sensitive methods are needed to detect the latter.