Introduction
In the marine communities, the digenetic trematodes (Digenea) are the most prevalent among parasitic organisms (Cribb et al. Reference Cribb, Cutmore and Bray2021). Particular interest to this group arises from their impressive diversity, complex life cycles, and high abundance of cryptic species (Blasco-Costa and Poulin Reference Blasco-Costa and Poulin2017; Pérez-Ponce de León and Poulin Reference Pérez-Ponce De León and Poulin2018). Investigations on complexes of cryptic or pseudocryptic (distinguishable morphologically but only following molecular distinction) species are promising not only for the elucidation of the actual biodiversity, but also for understanding speciation processes. The most suspicious to be cryptic complexes are the digeneans with broad host specificity, wide geographical range, and high morphological variability. Many representatives of the superfamily Hemiuroidea from the boreal and arctic seas match these criteria, and indeed, those studied by molecular methods turned out to be pseudocryptic (Bouguerche et al. Reference Bouguerche, Huston, Karlsbakk, Ahmed and Holovachov2024; Krupenko et al. Reference Krupenko, Kremnev, Gonchar, Uryadova, Miroliubov, Krapivin, Skobkina, Gubler and Knyazeva2022a). For others, however, molecular data are sparse.
Previously, based on both molecular and morphological observations, we found out that digeneans from the genus Lecithaster in the White Sea likely belong to the species L. salmonis Yamaguti, Reference Yamaguti1934 (Krupenko et al. Reference Krupenko, Kremnev, Skobkina, Gonchar, Uryadova and Miroliubov2022b). This result suggested a wide distribution of the species, including the Arctic and Pacific Northwest (seas of the Far East). Based on infection rates data, fish of the subfamily Salmoninae were detected as the most favourable definitive hosts of L. salmonis. The life cycle of this species was also elucidated: naticid gastropods Cryptonatica affinis (Gmelin, 1791) serve as its first intermediate hosts and planktonic copepods as the second intermediate hosts (Krupenko et al. Reference Krupenko, Kremnev, Skobkina, Gonchar, Uryadova and Miroliubov2022b). Interestingly, in the isolates from the White Sea, we observed minor variation in 28S rDNA and ITS2 fragments, which was consistent neither with the definitive host species nor with the morphological features. This could be a sign of cryptic speciation.
In the present paper, we explore genetic diversity of L. salmonis in seas of the Arctic and the Far East through mitochondrial (cox1 and nad5) and nuclear (28S, ITS1, ITS2) markers to estimate the geographic range of L. salmonis and to check whether it represents a complex of cryptic or pseudocryptic species. Next, we study the morphology of L. salmonis sexual adults (maritae) to find out whether there are differences between the genetically distinct groups. Finally, we describe cercariae and sporocysts of L. salmonis from a new intermediate host, Nucella freycinetii (Deshayes, 1839) and emphasize the importance of the data on intramolluscan life cycle stages of the Digenea for species recognition.
Material and methods
Sampling
Samples were collected in 2019–2023 from the northern European seas (White Sea, Barents Sea, Pechora Sea), the East Siberian Sea, and from the Far East (Sea of Okhotsk, Sea of Japan) (Table 1). The lists of the examined definitive and first intermediate hosts (fishes and gastropods) are shown in Supplementary Tables S1 and S2, respectively. The number of fish species hosting L. salmonis, of fish species examined, was 11 of 17 in the White Sea; 2 of 8 in the Barents Sea; 2 of 5 in the Pechora Sea; 2 of 4 in the East Siberian Sea; 3 of 16 in the Sea of Okhotsk; and 2 of 24 in the Sea of Japan (Table S1). For several fish hosts, we do not have precise infection intensity data, so the estimates in the table are based on the number of preserved worms. We primarily expected to find Lecithaster sporocysts and cercariae in Cryptonatica and other naticid gastropods, in Nucella (Muricidae), and in Odostomia (Pyramidellidae), considering the previous studies (Ching Reference Ching1960; Hunninen and Cable Reference Hunninen and Cable1943; Køie Reference Køie1989; Krupenko et al. Reference Krupenko, Kremnev, Skobkina, Gonchar, Uryadova and Miroliubov2022b). However, we also inspected all other available gastropods, representatives of the families Margaritidae, Trochidae, Solariellidae, Colloniidae, Tegulidae, Capulidae, Littorinidae, Velutinidae, Buccinidae, Nassariidae, Columbellidae, Cylichnidae, Scaphandridae, Lottiidae, and Batillariidae (Table S2). Thus, infection by Lecithaster-alike sporocysts and cercariae was searched for in 15 gastropod species from the White Sea, 13 species from the Barents Sea, 7 species from the Sea of Okhotsk, and 16 species from the Sea of Japan.
Geographical locations of sampling sites
Specimens from our previous study on Lecithaster salmonis in the White Sea (Krupenko et al. Reference Krupenko, Kremnev, Skobkina, Gonchar, Uryadova and Miroliubov2022b) were also analysed, including intramolluscan stages, maritae, and two metacercariae from planktonic copepods.
Live maritae were heat-killed prior to fixation in 96% ethanol. Seventeeen maritae specimens were obtained from frozen fish. Sporocysts with cercariae putatively belonging to Lecithaster were fixed in 96% ethanol. Cercariae were also fixed in 2.5% glutaraldehyde in sea water for morphological studies.
Molecular analysis
A complete list of isolates involved in the molecular analysis is presented in Table 2. To extract DNA, we used fragments of maritae (piece from the posterior end of the body, behind the vitellarium, or sometimes from the anterior extremity) and sporocysts fixed with 96% ethanol. The samples were dried completely in 1.5-ml tubes, incubated in 200 μl of 5% solution of Chelex® 100 resin (Bio-Rad, Hercules, California, USA) with 0.2 mg/ml proteinase K (Evrogen, Moscow, Russia) at 56°C for 3–4 h, then heated for 8 min at 90°C and centrifuged for 10 min at 16,000g. The supernatant with DNA was transferred into a new tube and stored at −20°C.
Isolates, their origin, and GenBank accession numbers for mtDNA sequences. OTUs recovered from mtDNA analysis are provided according to the Results section (see below).
* Isolates of Atopkin et al. (Reference Atopkin, Nakao, Besprozvannykh, Ha, Nguyen and Sasaki2020).
We amplified fragments of the cox1 and nad5 mitochondrial genes, fragments of 28S rDNA (either D2, or D1–D3 domains), fragments containing complete ITS1 (with partial 18S and 5.8S rDNA) and ITS2 (with partial 5.8S and 28S rDNA), with primers and conditions listed in Table 3. The draft assembly of the L. salmonis mitochondrial genome (our unpublished data) was used to design primers for the fragments of cox1 and nad5 mitochondrial genes, in Geneious Prime 2023.2.1 software (https://www.geneious.com).
Primers and PCR temperature profiles used in this study
The PCR mixture contained 4 μL of ScreenMix-HS (Evrogen), 0.5 μL of each primer (10 pmol μL−1), 2 μL of DNA, and 13 μL of PCR-grade water (Evrogen). PCRs were run on a BioRad T100 thermal cycler (Bio-rad Laboratories Inc., USA). PCR products were stained with ethidium bromide and visualized through electrophoresis in a 1% agarose gel. Sequencing was performed with the PCR primers on an AB3500xL genetic analyser (Applied Biosystems, MA, USA). Geneious Prime 2023.2.1 was used to analyse chromatograms and to build alignments. The relevant data for comparison were obtained from GenBank. Annotated sequences of the ribosomal operon (KR703279) and mitochondrial genome (KR703278) of Brachycladium goliath (van Beneden, 1858) Fraija-Fernández, Aznar, Raga, Gibson & Fernández, 2014 (Briscoe et al. Reference Briscoe, Bray, Brabec and Littlewood2016) were used to map our sequences. Positions with double peaks were checked twice via independent PCR runs and sequencing.
Pairwise genetic distances (as the number of base differences per site) were calculated in MEGA 11 (Tamura et al. Reference Tamura, Stecher and Kumar2021). Haplotype networks were constructed in PopART 1.7 (Leigh and Bryant Reference Leigh and Bryant2015) with the TCS network algorithm (Clement et al. Reference Clement, Snell, Walke, Posada and Crandall2002). Species delimitation was conducted on the basis of cox1 and nad5 gene fragments with a barcode gap approach using Assemble Species by Automatic Partitioning (ASAP) web interface with varying parameters (Puillandre et al. Reference Puillandre, Brouillet and Achaz2021).
Phylogenetic analysis was run using the cox1 dataset. The alignment comprised 83 newly obtained sequences (identical ones were removed) and sequences MH626006–8 of L. salmonis by Atopkin et al. (Reference Atopkin, Nakao, Besprozvannykh, Ha, Nguyen and Sasaki2020). The outgroups were representatives of other Lecithaster species, L. confusus Odhner, Reference Odhner1905 (MH625996–600), L. sayori Yamaguti, 1938 (MH626004) and Lecithaster sp. (MH626005), and Hysterolecithoides epinepheli Yamaguti, Reference Yamaguti1934 (all from Atopkin et al. Reference Atopkin, Nakao, Besprozvannykh, Ha, Nguyen and Sasaki2020). Substitution model selected by bModelTest (Bouckaert and Drummond Reference Bouckaert and Drummond2017) was TN93+G+I. The Bayesian inference (BI) analysis was conducted using Monte Carlo Markov Chain (MCMC) analysis available in Bayesian Evolutionary Analysis by Sampling Trees (BEAST2) (Bouckaert et al. Reference Bouckaert, Vaughan, Barido-Sottani, Duchêne, Fourment, Gavryushkina, Heled, Jones, Kühnert, De Maio, Matschiner, Mendes, Müller, Ogilvie, Du Plessis, Popinga, Rambaut, Rasmussen, Siveroni, Suchard, Wu, Xie, Zhang, Stadler and Drummond2019) on XSEDE at the CIPRES Science Gateway (https://www.phylo.org). Two independent runs of MCMC were performed, each with 10,000,000 generations and sampling every 1,000 generations. The trace files were checked for convergence with Tracer v1.7 (Rambaut et al. Reference Rambaut, Drummond, Xie, Baele and Suchard2018). The log files were combined using LogCombiner, discarding the first 10% as burn-in. Trees were summarized with TreeAnnotator using the maximum clade credibility tree option and with node heights as mean heights. The Maximum likelihood (ML) analysis was run in IQ-TREE 1.6.12 (Nguyen et al. Reference Nguyen, Schmidt, von Haeseler and Minh2015) with UFBoot2 (ultrafast bootstrap) (Hoang et al. Reference Hoang, Chernomor, von Haeseler, Minh and Vinh2018) with 1,000 replicates.
Morphology
Cercariae morphology was studied on temporary mounts of live specimens or those fixed with 2.5% glutaraldehyde in sea water. Maritae and sporocysts voucher specimens were stained with acetocarmine (Sigma Aldrich, Germany), destained in 0.1 M HCl in 70% ethanol, dehydrated in a graded alcohol series, clarified in xylol, and mounted in BioMount medium (Bio Optica, Italy). Voucher specimens are deposited in the Helminths collection of the Zoological Institute of the Russian Academy of Sciences (ZISP), section Trematoda (non-type material). The whole mounts were studied under Leica DM 500 or Leica DM 2500 compound microscopes (Leica Microsystems, Germany) in bright field or with differential interference contrast (DIC). Photographs were taken using a Nikon DS Fi3 (Nikon, Tokyo, Japan) or with a smartphone camera. Measurements were made by using Fiji software (Schindelin et al. Reference Schindelin, Arganda-Carreras, Frise, Kaynig, Longair, Pietzsch, Preibisch, Rueden, Saalfeld, Schmid, Tinevez, White, Hartenstein, Eliceiri, Tomancak and Cardona2012). All measurements are in micrometres. In the cercariae description, means are followed by ranges of values in parentheses.
We compared maritae of genetically distinct groups (based on mtDNA) to determine whether they differ in morphological features. To minimize host-dependent morphological variation, we used only the maritae from the Salmoninae fish (Oncorhynchus keta, O. gorbuscha, and Salvelinus malma). Characteristics of the reproductive system (except egg size) are highly variable (Krupenko et al. Reference Krupenko, Kremnev, Skobkina, Gonchar, Uryadova and Miroliubov2022b), and from the present data on cox1 (see below), we can say that the whole range of variation was observed within one genetically distinct group. Thus, in this study, only the following somatic characters were analysed: body length and body maximal width, forebody length, length and width of oral and ventral suckers, length and width of pharynx, and length of preoral lobe. The egg size was analysed separately. Multidimensional statistical methods were used to analyse the morphological data (Quinn and Keough Reference Quinn and Keough2002). All calculations were performed using the R environment (R Core Team Reference Team2023); graphs were created with the ‘ggplot2’ package (Wickham Reference Wickham2016), if not stated otherwise. P-values lower than 0.05 were considered significant. The p-values in multiple tests were adjusted with the Holm-Bonferroni correction. The morphometric data were explored with principal components analysis of the correlation matrix with the help of vegan R package (Oksanen et al. Reference Oksanen, Simpson, Blanchet, Kindt, Legendre, Minchin, O’Hara, Solymos, Stevens, Szoecs, Wagner, Barbour, Bedward, Bolker, Borcard, Carvalho, Chirico, De Caceres, Durand, Evangelista, FitzJohn, Friendly, Furneaux, Hannigan, Hill, Lahti, McGlinn, Ouellette, Ribeiro Cunha, Smith, Stier, Ter Braak and Weedon2022) and visualized with ggfortify package (Tang et al. Reference Tang, Horikoshi and Li2016). Before the analysis, the raw data were double-centred to remove the effect of specimen body size. To test whether the distinct genetic lineages showed differences in morphology, we used permutational multidimensional analysis of variance (PERMANOVA) on the matrix of euclidean distances with 1,000 permutations from the vegan R package. The same analysis was used to assess the effect of the host species on the maritae morphology. One-factor linear model with one predictor was fitted to test the effect of a genetic lineage on the mean egg length of a marita specimen.
Results
Overview
Maritae consistent with the diagnosis of Lecithaster salmonis were recovered from 17 of 55 examined fish species (Table S1). The highest intensities were in fishes from the subfamily Salmoninae. Maritae supposedly belonging to other species of the genus Lecithaster were not found. In the Far East, L. salmonis maritae were recovered only from Salmoninae fish. Of 38 examined species of gastropods from all sampled localities, putative life cycle stages of Lecithaster were found only in Cryptonatica affinis from the White and Barents Seas, and in Nucella freycinetii from the Sea of Okhotsk (Table S2). The number of isolates involved in the molecular analysis from each host-locality combination is shown in Supplementary Tables S1 and S2. Twelve isolates from the White Sea are from our previous work on L. salmonis (Krupenko et al. Reference Krupenko, Kremnev, Skobkina, Gonchar, Uryadova and Miroliubov2022b), including two metacercariae from planktonic copepods, Metridia longa, and Pseudocalanus sp. We use the same isolate IDs as in the previous work.
cox1 gene
Fragment of the cox1 gene was sequenced for 118 isolates, trimmed to 594 bp, and submitted to GenBank under accession numbers PQ165659–776. Sequences MH626006–8 of L. salmonis by Atopkin et al. (Reference Atopkin, Nakao, Besprozvannykh, Ha, Nguyen and Sasaki2020) were also added into the alignment, though they did not match our fragment completely: they were shorter by 81 bp at the 5’-end. The results of the cox1 analysis are summarized as a haplotype network. It comprises 79 haplotypes (Figure 1). In the ASAP analysis, the best score suggested that the dataset was most likely partitioned into six groups with a threshold distance of 11.9%.
Haplotype network of Lecithaster cf. salmonis isolates (N=121) based on partial cox1 gene sequences. Circle size represents the haplotype frequency. Black dots indicate missing haplotypes. Number of hatch marks corresponds to the number of substitutions between haplotypes. Coloured background indicates OTUs derived from ASAP. Isolates taken from intermediate hosts labelled in purple and from definitive hosts labelled in black. Asterisks mark the isolates obtained from the single specimen of chum salmon from the Sea of Japan. Abbreviations: Clu_pal, Clupea pallasii; Cor_lav, Coregonus lavaretus; Cry_aff, Cryptonatica affinis; Ele_naw, Eleginus nawaga; Gad_mor, Gadus morhua; Hyp_jap, Hypomesus japonicus; Lim_lim, Limanda limanda; Met_lon, Metridia longa; Myo_qua, Myoxocephalus quadricornis; Myo_sco, Myoxocephalus scorpius; Nuc_fre, Nucella freycineti; Onc_gor, Oncorhynchus gorbuscha; Onc_ket, Oncorhynchus keta; Onc_kis, Oncorhynchus kisutch; Osm_mor, Osmerus mordax; Pse_sp, Pseudocalanus sp.; Sal_mal, Salvelinus malma; Sal_sal, Salmo salar.
These six haplogroups are marked as operational taxonomic units (OTU) 1–6 on the haplotype network. Only one of them, OTU1, comprised isolates from all the hosts within a life cycle: Cryptonatica affinis as the first intermediate, planktonic copepods as the second intermediate, and various fish species as the definitive hosts. Four of the haplogroups (OTU2, OTU3, OTU5, OTU6) included only the maritae specimens; of these, OTU5 was represented by a single isolate. Noteworthy, some maritae which belonged to the different OTUs originated from the same definitive host specimen. The most remarkable case was a single specimen of the chum salmon Oncorhynchus keta from the Sea of Japan, which was parasitised by maritae from four different OTUs (marked with an asterisk in Figure 1). OTU4 was represented only by the isolates from the first intermediate host Nucella freycinetii.
The OTU1 was encountered in five seas: White, Barents, Pechora, Sea of Okhotsk, and Sea of Japan (Figure 2). OTU2 was found in the the White, Pechora, and East Siberian Seas. OTU3 was found in the Sea of Japan and the Sea of Okhotsk. All specimens of OTU4 originated from the Sea of Okhotsk. OTU5 and OTU6 were found only in the Sea of Japan (the single specimen of chum salmon). Three sequences by Atopkin et al. (Reference Atopkin, Nakao, Besprozvannykh, Ha, Nguyen and Sasaki2020) fell into the OTU6, and they originated from a nearby region, the Pacific coast of Japan.
Geographic distribution of the OTUs derived from cox1 gene analysis. Circle size represents the number of isolates derived from each region. The sectors of the pie chart are proportional to the ratio of the OTUs in the region.
In the haplotype network, OTU1 had a distinct central haplotype and formed a star-like pattern (Figure 1). OTU6 showed a similar pattern, whereas OTU4 and OTU5 had a limited number of isolates, and no pattern was identified. The OTU2 had many unique haplotypes and lacked a dominant one.
The barcode gap was not distinct for the OTUs derived from ASAP: the maximal intragroup pairwise distance was in OTU2 (0.0171±0.0053, 10 substitutions), and the minimal intergroup pairwise distance was between OTU1 and OTU3 (0.0155±0.052, 9 substitutions).
The results of BI analysis of cox1 data are shown in Figure 3. The position of the sequences from the OTU2 is unresolved within the clade OTU2+OTU4, but all the other nodes are well-supported. The results of ML analysis do not fully match BI: sister relationships were not supported between OTU2 and OTU4, and between OTU1 and OTU3.
Bayesian phylogenetic tree of Lecithaster cf. salmonis inferred from cox1 gene data. Posterior probabilities are shown at nodes, followed by bootstrap support from the tree built for the same dataset with ML method. Support values lower than 0.9 (BI) and 75 (ML) are not shown. Scale bar shows the substitution rate. Branches are collapsed within OTUs which appear monophyletic and supported. GenBank accession numbers for sequences from OTU1, OTU3, OTU4, and OTU6 are in Table 2.
nad5 gene
To see if the OTUs derived from the cox1 analysis would be consistent with the results based on other mitochondrial markers, we sequenced a 603 bp-long fragment of the nad5 gene for 21 isolates (PQ231533–53). We picked three to five isolates from each OTU (except OTU5 with a single isolate), trying to embrace the distant haplotypes within the cox1–haplogroups and thus avoid artificial grouping in the nad5 analysis.
The resulting haplotype network comprised 21 haplotypes, so each one was unique (Figure 4). In the network, six OTUs are easily defined. ASAP run with the K80 substitution model gave the best partition into six groups (score 1.50) – the same as those derived from cox1 analysis. However, when run with other parameters (JC69 or p-distance), ASAP gave the same best score (1.50) to partitions with five and six groups. In the five-group partition, OTU2 and OTU5 were together. The barcode gap was evident in the six-group partition: the maximal intragroup pairwise distance was in OTU1 (0.0152±0.0050, 9 substitutions), and the minimal intergroup pairwise distance was between OTU2 and OTU5 (0.0221±0.064, 13 substitutions). In the five-group partition, the barcode gap was obscure: 14 substitutions were minimal intergroup distance (between OTU1 and OTU3), and maximal distance within group was 13 substitutions between isolates D35.2 and D39.6 (belonging to OTU2 and OTU5, respectively, in the cox1 analysis).
Haplotype network of Lecithaster cf. salmonis isolates (N=21) based on partial nad5 gene sequences. Circle size represents the haplotype frequency. Black dots indicate missing haplotypes. Number of hatch marks corresponds to the number of substitutions between haplotypes. Coloured background indicates OTUs derived from ASAP. Isolates taken from the first intermediate host labelled in purple and from the definitive host labelled in black. Asterisks mark the isolates obtained from the single specimen of chum salmon from the Sea of Japan. Abbreviations for host species are the same as in Figure 1.
There was no close enough outgroup with sequenced nad5 gene to run BI or ML analyses.
Nuclear markers
We sequenced three fragments of ribosomal DNA for several isolates of each OTU derived from the mtDNA analyses. A long fragment of 28S rDNA (D1–D3) was sequenced for ten isolates, and for one more isolate, only a short one (D2) was obtained (PQ165648–58). Several previously published sequences of our isolates (OM850373, OM850382, and OM850385–6 from Krupenko et al. Reference Krupenko, Kremnev, Skobkina, Gonchar, Uryadova and Miroliubov2022b) and MH625980–1 of Atopkin et al. (Reference Atopkin, Nakao, Besprozvannykh, Ha, Nguyen and Sasaki2020) were also included into the analysis. The alignment of long fragments did not yield any extra differences compared to the short one, except for the probably low-quality reads at the 5’- and 3’-ends. Thus, our alignment was trimmed to 425 bp. Five polymorphic sites were present in this fragment (Figure 5). Sites with Y had double peaks indicative of heterozygosity in isolates from OTU1, OTU3, OTU4, and OTU6. Generally, both inter- and intragroup variability was evident in this fragment, and demarcation of the OTUs discovered in the mtDNA analyses was not clear.
Polymorphic sites within alignment of three rDNA fragments (ITS1, ITS2, and 28S). OTUs from cox1 and nad5 analysis are marked by a grey outline. Newly generated sequences are in bold. Asterisks mark isolates from Atopkin et al. (Reference Atopkin, Nakao, Besprozvannykh, Ha, Nguyen and Sasaki2020).
For the same set of isolates, a fragment containing ITS2 flanked with the partial 5.8S and 28S rDNA was analysed. The dataset comprised newly obtained sequences (PQ165626–35) and ones from the previous publications (MH625993–4 of Atopkin et al. Reference Atopkin, Nakao, Besprozvannykh, Ha, Nguyen and Sasaki2020; OM850400, OM850407, OM850409, and OM850412–13 of Krupenko et al. Reference Krupenko, Kremnev, Skobkina, Gonchar, Uryadova and Miroliubov2022b). The sequences were trimmed, so the alignment was 556 bp long. Sequences MH625993–4 were shorter by 73 bp at the 3’-end. Ten polymorphic sites were present (Figure 5). In five of these, double peaks indicated possible heterozygosity. OTU2, OTU4, and OTU5 were identical in the analysed fragment. The other three OTUs differed from those and from each other, and also showed intragroup variation.
The ITS1 region flanked with the partial 18S and 5.8S rDNA was successfully sequenced only for 12 isolates from the same set (PQ165636–47). Yet, all six OTUs were covered. After trimming the alignment to the shortest sequence, it was 656 bp long and included 11 polymorphic sites. Two sequences from OTU1 and one sequence from OTU3 were identical. Intragroup variation was present within OTU1, OTU2, and OTU3. Double peak was present in one sequence (Y in the position 203 in D35.2). One position in D4.9 was ambiguous (M).
Morphometry of maritae
For the morphometric analysis, we took maritae from four genetically distinct groups: OTU1, OTU2, OTU3, and OTU6. In this analysis, we had OTU3 and OTU6 only represented by specimens from O. keta, OTU1 mostly represented by specimens from O. gorbuscha, with one specimen from O. keta, and OTU2 from O. gorbuscha and S. malma. Specimens from non-salmonine fishes were not taken into the morphometric analysis to reduce the effect of host species on the results. Only heat-killed maritae were used. Measurements are in Table 4.
Morphometric data for genetically distinct groups from mtDNA analysis. Means with 95% confidence intervals for normal distribution are given, followed by the range of values in parentheses. Ten eggs measured from each specimen
In the PCA analysis, the first two components explained about 90% of the variation (Figure 6a). The first component had positive loadings for the oral sucker length and width, pharynx length and width, and preoral lobe length, and negative loading for the body length. The second component was positively correlated with the forebody length and negatively with the body width. Two distinct groups were visible in the biplot: (1) OTU3 and OTU6, (2) OTU1 and OTU2. Maritae belonging to the first group seemed to have a shorter preoral lobe and smaller oral sucker and pharynx compared to body length, and larger body width to forebody length ratio.
Result of comparing morphometric data from four OTUs. (a) Biplot of a principal components analysis. The ellipses represent the 95% confidence intervals of the groups’ centroids. (b) P-values from post-hoc comparisons after PERMANOVA on the matrix of euclidean distances between isolates, values lower than 0.05 are in bold. (c) Boxplot of mean egg length of maritae belonging to four OTUs. Abbreviations: Onc_gor, Oncorhynchus gorbuscha; Onc_ket, Oncorhynchus keta; Sal_mal, Salvelinus malma.
PERMANOVA generally supported this grouping. There was a significant effect of the OTU on morphological features of the maritae (F3=8.7, p=0.001, R2=0.52). Post hoc comparisons showed that morphology of OTU2 significantly differed from that of OTU3 and OTU6, and morphology of OTU1 from that of OTU3, while other groups were not significantly different (Figure 6b). Host species (O. gorbuscha or S. malma) did not significantly affect the morphological features of OTU2 maritae (F1=3.9, p=0.2).
Egg length was analysed for ten maritae specimens from OTU1, six from OTU2, eight from OTU3, and eight from OTU6 (Table 4). Ten eggs were measured from each specimen. We did not find the effect of OTU on mean egg length (F3 = 1.2, p=0.33) (Figure 6c).
Lecithaster cf. salmonis OTU1 intramolluscan stages
Previously, we described sporocysts and cercariae of L. salmonis from naticid gastropod Cryptonatica affinis (Krupenko et al. Reference Krupenko, Kremnev, Skobkina, Gonchar, Uryadova and Miroliubov2022b). Newly obtained molecular data on cox1 are based on isolates from the previous paper as well as new ones from C. affinis, and all of them belong to the OTU1.
Host: Cryptonatica affinis (Gmelin, 1791)
Known localities: Keret Archipelago (Kandalaksha Bay, White Sea)
New localities: Velikaya Salma Strait (Kandalaksha Bay, White Sea); Dalniye Zelentsy (Kola Peninsula, Barents Sea).
Prevalence: 16.1% (N = 56) in Velikaya Salma Strait; 1.6% (N = 247) at Dalniye Zelentsy.
Sites: Reproductive gland, digestive gland, albumen gland, kidney.
Synonymized names: Cercaria saccocaudata Chubrik, 1966
Lecithaster cf. salmonis OTU4 intramolluscan stages (Fig. 7)
Host: Nucella freycinetii (Deshayes, 1839)
Locality: Nagaev Bay (Taui Bay, Sea of Okhotsk).
Prevalence: 3% (N = 165).
Sites: Reproductive and digestive glands.
Sporocysts very long (Figure 7d), easily torn. Four whole specimens fixed with ethanol 3834–8907 long, with 173–203 maximal width. Parts of sporocysts up to 11,550, so perhaps actually longer. Birth pore terminal at anterior body end.
Intramolluscan life cycle stages of Lecithaster cf. salmonis OTU4. (a–c) Schemes of pre-infective (a), infective (b) cercariae and invasion process (c). (d) Sporocyst stained with acetocarmine, arrowhead marks anterior end. (e–g) Cercarial developmental stages (e, f) and pre-infective cercaria (g) in DIC. (h) Live pre-infective and infective cercariae. (i) Delivery tube eversion in live cercaria under cover glass pressure. Abbreviations: ap, aperture; c, caecum; cb, cercarial body; cc, caudal cyst; dt, delivery tube; d1–d2, caudal cyst diameters; eg, excretory granules; ev, excretory vesicle; gs, granulated structure; fa, filamentous appendage; ic, inner cyst layer; oc, outer cyst layer; os, oral sucker; ph, pharynx; vs, ventral sucker.
Cercariae of cystophorous type, with tail comprising caudal cyst and delivery tube. Measurements of pre-infective cercariae based on 30 glutaraldehyde-fixed specimens from two different host individuals. Pre-infective cercariae with filamentous appendage 31 (25–38) long, attached near cyst aperture (Figure 7a, g). During development, appendage forming as dorsal outgrowth on caudal cyst (Figure 7e, f). Body of pre-infective cercariae 126 (114–142) long, 31 (28–35) wide. Oral sucker 16 (15–17) × 18 (15–20), pharynx 8 (7–10) × 10 (8–12), ventral sucker 19 (18–21) × 19 (17–21). Caeca primordia reaching middle of forebody (Figure 7a). Excretory vesicle with enlarged distal part 14 (11–18) long. Excretory ducts containing large spherical granules (Figure 7a, g). Excretory formula not determined. Within caudal cyst of pre-infective cercaria, granulated structure present near aperture (Figure 7a, g).
Measurements of infective cercariae based on 16 live specimens from one host individual. Caudal cyst comprising two layers (Figure 7b, h). Inner cyst firm, crescent-shaped, with one end roundish and another pointed, 80 (74–83) × 79 (74–84) (measured as shown on Figure 7h), 53 (50–57) thick. Outer cyst thin, probably inflated in seawater, 94 (88–102) in diameter, 71 (68–78) thick.
Delivery tube measured in live specimens when everted under coverslip pressure (N = 12), 446 (423–483) long. Two projections one opposite to another at 96 (85–108) from base (Figure 7c, i). Delivery tube maximal width 26 (23–28), proximal to projections.
Discussion
Studying intraspecific diversity and cryptic species complexes is essential for understanding the speciation processes in parasitic organisms. In the present paper, we access this problem by exploring the geographical distribution and genetic and morphological variation in a fish digenean Lecithaster salmonis. This species was first described by Yamaguti (Reference Yamaguti1934) from Oncorhynchus keta in Japan (specific locality unknown). Later, its molecular characterization was performed through 28S rDNA, ITS2, and cox1 gene (Atopkin et al. Reference Atopkin, Nakao, Besprozvannykh, Ha, Nguyen and Sasaki2020). Based on the similarity in ribosomal markers and morphological characteristics, it was hypothesized that L. salmonis is common in the White Sea, and the other species of the Lecithaster are absent there (Krupenko et al. Reference Krupenko, Kremnev, Skobkina, Gonchar, Uryadova and Miroliubov2022b). Now, we showed that digeneans recognized as L. salmonis are also found in Barents, Pechora, East Siberian Seas, Sea of Okhotsk, and Sea of Japan. The highest infection rates were observed in the fishes from the subfamily Salmoninae.
Morphological distinction of L. salmonis from other species of the same genus can rely on the egg size: 22–24 from the original description (Yamaguti Reference Yamaguti1934), 19–24 from molecular-based research (Atopkin et al. Reference Atopkin, Nakao, Besprozvannykh, Ha, Nguyen and Sasaki2020; Krupenko et al. Reference Krupenko, Kremnev, Skobkina, Gonchar, Uryadova and Miroliubov2022b; present data). Characteristics of vitellarium and male duct do not allow species identification. Considering this, the several previous records of Lecithaster spp. from the Arctic and Pacific Northwest may actually pertain to L. salmonis. Specifically, Zhukov (Reference Zhukov1960, Reference Zhukov1963) regarded L. salmonis as a junior synonym of L. gibbosus; however, his specimens from the Bering Sea had eggs 20–22 long, and the highest prevalence was registered in Salvelinus malma. From the Sea of Japan and at the Shikotan Island, Zhukov (Reference Zhukov1960) described Lecithaster with egg length 21–25, and some of his isolates could also belong to L. salmonis rather than L. gibbosus (egg length 25–27 in the original description by Odhner (Reference Odhner1905)). The existence of L. gibbosus in the Pacific is actually doubtful as its first intermediate host is Brachystomia eulimoides (Hanley, 1844), found only in the North East Atlantic (Høisæter Reference Høisæter2014). High prevalence and intensity of Lecithaster in S. alpinus from the Kara Sea (Issaitschikov Reference Issaitschikov1933) may also suggest that it was actually L. salmonis. The same is true for the data of Shulman and Shulman-Albova (Reference Shulman and Shulman-Albova1953) on Salmo salar from the White Sea.
The present study confirmed previous data on how to tell apart L. salmonis from the other species of the genus; however, new doubts arose about the integrity of this species. Our results are controversial in answering the question whether L. salmonis represents a single species or is a complex of cryptic or pseudocryptic species. Mitochondrial markers (cox1 and nad5) show that there are likely six lineages within our sample (OTU1–6), but the barcode gap for them is not clear. BI and ML analyses of cox1 mtDNA do not match each other well. Both analyses support monophyly within three lineages, OTU1, OTU3, and OTU6; BI also gives evidence that OTU4 is monophyletic. OTU5 was represented by a single isolate; it forms a well-supported clade together with OTU6. The position of isolates from OTU2 is unresolved, likely due to high variation within this lineage. Probably, there is not enough data yet for a reasonable conclusion on the phylogeny within L. salmonis, and the closest available outgroups differ too much from it.
Nuclear ribosomal markers give an even more complex picture. They demonstrate substantial variability, but it does not match groups derived from the cox1 and nad5 analyses. Different OTUs may include isolates which are identical in all three analysed markers (D4.8 and D10.3 from OTU1 and D40.7 from OTU3). At the same time, each OTU had intragroup variability (except OTU5 represented by a single specimen). Ambiguities with double peaks were abundant in the analysed ribosomal fragments, and they may indicate heterozygotes. Thus, cross-fertilization between individuals from different OTUs may be hypothesized. However, there is an alternative explanation of the double peaks, an intra-individual variation of rDNA which was described in various organisms (e.g., Bik et al. Reference Bik, Fournier, Sung, Bergeron and Thomas2013; Carranza et al. Reference Carranza, Baguñà and Riutort1999; Keller et al. Reference Keller, Chintauan-Marquier, Veltsos and Nichols2006; Simon and Weiß Reference Simon and Weiss2008).
Specialization of the OTUs to certain different species of the definitive hosts can be discussed as a basis for the divergence. However, our data show that it is unlikely. OTU1 and OTU2 were encountered in a variety of fish species from different taxa. We found OTU6 only in O. keta, but isolates of Atopkin et al. (Reference Atopkin, Nakao, Besprozvannykh, Ha, Nguyen and Sasaki2020) originate from Osmeridae fish. In OTU3, we also have specimens only from O. keta, but the sample size is too small to state it as a specific definitive host.
Analysis of maritae morphometry was not helpful in delimiting the OTUs detected with mitochondrial gene data. The egg size did not show any significant differences between the OTUs. Results from PCA and PERMANOVA based on morphometry of somatic characters actually contradicted genetic data. Maritae of OTU3 and OTU6 do not differ significantly in morphology, although these groups clearly diverged in the cox1 and nad5 analyses. Simultaneously, OTU1 and OTU3 are genetically close (both in mitochondrial and ribosomal markers) but demonstrate significant differences in the morphometric characters. Here, we should note that the morphology may be affected not only by the genetic divergence but by the other factors, like ‘crowding effect’ or ‘host-induced variation’ (Stunkard Reference Stunkard1957). We suppose the effect of infection intensity is negligible for our sample: all the specimens in the morphometric analysis originate from hosts infected with dozens of L. salmonis maritae. The host-induced variation cannot be ruled out; it occurs widely in the Digenea (Blasco-Costa and Locke Reference Blasco-Costa and Locke2017; Cribb et al. Reference Cribb, Bray, Justine, Reimer, Sasal, Shirakashi and Cutmore2022; Hildebrand et al. Reference Hildebrand, Adamczyk, Laskowski and Zaleśny2015). To reduce the effect of the host species on morphology, we took into the analysis only the isolates from fish belonging to the subfamily Salmoninae, most favourable hosts for L. salmonis. However, that does not mean that all salmonines are equally favourable. Our sample is not sufficient to make a wide-scale analysis of fish species’ effect on maritae morphology. Only for the OTU2, such analysis was performed (three specimens from each O. gorbuscha and S. malma) and did not show significant differences. Still, similarity of OTU3 and OTU6 maritae may arise from the fact that all specimens were obtained from O. keta (moreover, a single specimen), and the vast majority of isolates from other OTUs were obtained from different salmonine species.
Finally, a set of characters which seems the most promising to determine species boundaries in our sample is linked with the cercariae morphology and identity of the first intermediate host. These data were collected only for OTU1 and OTU4, and they suit well to tell one from another. The shape of the caudal cyst and everted delivery tube and the length of the filamentous appendage are the most differential characteristics of cercariae. Also, OTU1 and OTU4 utilize representatives of two different orders of gastropods as the first intermediate hosts: Cryptonatica affinis (Littorinimorpha) and Nucella freycinetii (Neogastropoda), respectively. Such broad specificity is extremely rare for the Digeneans. The most impressive example is Microphallus pseudopygmaeus Galaktionov, 1980, which utilizes representatives of distant families of the Littorinimorpha (Littorinidae, Hydrobiidae, Rissoidae, Naticidae, Falsicingulidae) as well as Trochida (families Margaritidae, Solariellidae) (Galaktionov et al. Reference Galaktionov, Blasco-Costa and Olson2012). However, there is no apparent morphological divergence between M. pseudopygmaeus from different intermediate hosts like the one we see in L. salmonis. We regard differences in cercariae structure together with host species substantial enough to claim that OTU1 and OTU4 belong to different species.
Data on intramolluscan stages are essential to understand if the other OTUs deserve recognition as separate species. Previously, Ching (Reference Ching1960) described lecithasterid Cercaria A and Cercaria B in the Pacific Northeast from two other Nucella species: N. emarginata (Deshayes, 1839) and N. lamellosa (Gmelin, 1791), respectively. Of these, the former strongly resembles cercaria of OTU4. The second one may be the larva of another OTU within L. salmonis. Despite that OTU2 and OTU4 were quite close in our molecular analysis, it is highly unlikely that any species of the Nucella is a host for OTU2 because these gastropods are absent in the White Sea (except for the part close to the Barents Sea) and East Siberian Sea.
The host distribution is one of the most important factors that constrain the range of parasites. So, to discuss the distribution of L. salmonis OTUs, we need to overview the data on the specificity for different hosts in the life cycle. The specificity for the definitive host is wide (as discussed above). The role of the second intermediate host can be performed by copepods from different genera (present data on OTU1). Experimental data on different species of Lecithaster show that successful infections can be performed with quite random copepod species (Ching Reference Ching1960; Hunninen and Cable Reference Hunninen and Cable1943; Køie Reference Køie1989). However, some restrictions are evident: Boyce (Reference Boyce1969) showed that Centropages abdominalis Sato, 1913 and Pseudocalanus minutus (Krøyer, 1845) are preferred over Calanus pacificus Brodsky, 1948 and Neocalanus plumchrus (Marukawa, 1921). As for the first intermediate host, the specificity is probably much narrower (based on data on OTU1 and OTU4). That is why we expect the first intermediate host to confine the distribution of L. salmonis OTUs.
The most widely distributed across the studied area is OTU1, encountered in all the sampled regions except the East Siberian Sea, with the majority of isolates originating from the northern European seas. Utilization of C. affinis gives OTU1 an opportunity for wide distribution through the Arctic, North Pacific, and, possibly, North Atlantic. Isolates belonging to OTU2 were sampled in the White, Pechora, and East Siberian Seas. Absence of OTU2 isolates in the Barents Sea may be either due to the insufficient sample size from this area or to some environmental restrictions. The White, Pechora, and East Siberian Seas share a combination of two characteristics: low temperature and low salinity (Charkin et al. Reference Charkin, Kosobokova, Ershova, Syomin, Kolbasova, PYu, Leusov, Dudarev, Gulenko, Yaroshchuk, Startsev, Fayman, Krasikov, Zverev, Bessonova, Ulyantsev, Elovsky, Yurikova, Kobyakov, Zimina, Gerasimova, Tishchenko and Didov2024; Chernov and Tolstikov Reference Chernov and Tolstikov2020; Sukhotin et al. Reference Sukhotin, Denisenko and Galaktionov2019). These conditions may be favourable for OTU2 in some way, possibly for the distribution of its first intermediate host which is yet unknown. The presence of OTU3–6 only in the Far East may also be linked with the distribution of their first intermediate hosts. Numerous gastropods like Nucella freycinetii, the first intermediate host of OTU4, occur all through the Far East, from the Bering Sea to the Sea of Japan, but do not pass to the north from the Bering Strait (Kantor and Sysoev Reference Kantor and Sysoev2006). In general, geographical range is uneven in six OTUs recovered from mtDNA analysis, but their divergence is not linked to the present geographical isolation. The representatives of different OTUs are encountered in sympatry: OTU1 and OTU2 in the White and Pechora Seas (even together in one O. gorbuscha specimen); OTU1, OTU3, and OTU4 in the Sea of Okhotsk; OTU1, OTU3, OTU5, and OTU6 in the Sea of Japan (notably, in a single O. keta specimen).
The criteria for species delimitation are widely discussed and tested in the research on the digenean diversity (Bray et al. Reference Bray, Cutmore and Cribb2022; Cribb et al. Reference Cribb, Bray, Justine, Reimer, Sasal, Shirakashi and Cutmore2022; Cutmore et al. Reference Cutmore, Corner and Cribb2023; Pérez-Ponce de León et al. Reference Pérez-Ponce De León, Solórzano-García, Huston, Mendoza-Garfias, Cabañas-Granillo, Cutmore and Cribb2024). The most favoured combination includes molecular data, morphological differences, and distribution in the definitive hosts. The case of L. salmonis is instructive in highlighting contradictions between these characters. We regard this situation as a sign of recent or ongoing speciation process. We also suggest referring to the genetic lineages of L. salmonis defined through mtDNA as ‘Lecithaster cf. salmonis’ followed by an OTU number. Involvement of data on intramolluscan stages seems essential for such ambiguous cases. Also, an important thing to do in future is to figure out which OTU matches the type material of Yamaguti and must be referred to as L. salmonis sensu stricto.
The results of cox1 mtDNA analysis together with geographical distribution data may give some clues to the history of the speciation process within L. salmonis. Although the sample from the Far East was smaller than from the northern European seas, most OTUs were recovered from the former region. This makes us hypothesize that L. salmonis originates from the Pacific, which would resemble many other marine species with circumpolar or amphi-boreal distribution (Laakkonen et al. Reference Laakkonen, Hardman, Strelkov and Väinölä2021). The star-like pattern with prominent central haplotype in OTU1 suggests bottleneck occurrence. It could be a result of invasion of an ancestor from the Pacific into the Arctic and subsequent isolation during glaciation accompanied with sea level lowering and closure of the Bering Strait. Present occurrence of OTU1 in the Pacific would be a result of recent reinvasion in this case. Dominant haplotype and star-like pattern is also present in OTU6. Yet, we hesitate to go into the further speculations based on scarce data on this group and unclear phylogenetic relationships between lineages within L. salmonis.
Supplementary material
The supplementary material for this article can be found at http://doi.org/10.1017/S0022149X24000890.
Acknowledgements
We are grateful to colleagues and institutions who helped with or facilitated material collection, including Dr. Kira Regel, Dr. Aleksei Miroliubov, the Educational and Research Station ‘Belomorskaia’ of Saint Petersburg University (SPbU), the White Sea Biological Station (WSBS) ‘Kartesh’ of the Zoological Institute of the Russian Academy of Sciences (ZIN RAS), the N. A. Pertsov White Sea Biological Station of Lomonosov Moscow State University, Biological Research Station of Murmansk Marine Biological Institute of RAS in Dalnie Zelentsy, and Marine Biological Station ‘Vostok’ of Institute of Marine Biology of RAS. During research, we used the equipment of the resource centres ‘Molecular and Cell Technologies’ and ‘Centre for Culture Collection of Microorganisms’ of the Research Park of SPbU. We thank Dr. Anna Romanovich for excellent sequencing. We are also grateful to Dr. Elizaveta Skalon for consulting on the phylogeny reconstructions and to Dr. Isabel Blasco-Costa for an enthusiastic discussion on the earlier version of this manuscript after the ‘Trematodes-2024’ conference in Brisbane.
Financial support
The study was performed with financial support of Russian Science Foundation, project No 23-24-00376, https://rscf.ru/project/23-24-00376/.
Competing interest
The authors declare none.
Ethical standard
This study was conducted in compliance with all institutional, national, and international guidelines on the care and use of animals.
