INTRODUCTION
Didymozoids (Didymozoidea) are digeneans of oceanic pelagic fish with world-wide distribution in tropical to subtropical areas, considered one of the most taxonomically complex digenean families (Pozdnyakov and Gibson, Reference Pozdnyakov, Gibson, Bray, Gibson and Jones2008). In general, their life cycle follows typical digenean propagation where small pelagic fish and cephalopods act as intermediate and pelagic fish are final hosts (Lester and Newman, Reference Lester and Newman1986; Pozdnyakov, Reference Pozdnyakov1996). Apart from that, there are reasons to assume that they have acquired a three- or four-hosts type development: definitive hosts are orally infected and juvenile individuals follow complex migration from their digestive tract to the parasitation site, where two individuals enclose in a capsule, representing true maritae as they mature. However, at the initial stage of adult morphogenesis, young didymozoids have been recorded in the intestines of a wide number of small fish, suggesting that those juveniles are transported between crustaceans (first intermediate hosts) to large pelagic fish (definite hosts) by small fish. The fourth host, at least for some genera, may have emerged when an additional predator has been acquired into the chain of predators along which the parasite is propagated (see Galaktionov and Dobrovolskij, Reference Galaktionov and Dobrovolskij2010). The exact migration route of a didymozoid juvenile from final host's digestive tract to its predilected sites of parasitation is unknown given the sampling and experimental limitation related to large pelagic fish. Similarly, the mechanisms these hermaphroditic or gonochoristic species acquire to find their cospecific are still not understood.
Although data on the occurrence and effect of didymozoids in different hosts have been largely reported, few studies have focussed on the didymozoid phylogeny. It was suggested that didymozoids inhabiting Pacific (Thunnus orinetalis) and Atlantic bluefin reflect the habitat selection in their phylogenetic relationships (Mladineo et al. Reference Mladineo, Bott, Nowak and Block2010) that consequently resulted in the remarkable number of sites they colonize in hosts (Mladineo, Reference Mladineo2006).
In the Atlantic bluefin tuna (Thunnus thynnus, Scombridae), most didymozoid species are encapsulated in pairs in a connective-tissue capsule (Di Maio and Mladineo, Reference Di Maio and Mladineo2008) in skin, muscle and external mucoses, digestive tract, peripheral nervous tissues and kidney, causing different degrees of pathology (Mladineo et al. Reference Mladineo, Žilić and Čanković2008). Didymosulcus katsuwonicola (syn. Didymocystis wedli) is the most abundant didymozoid species infecting farmed bluefin tuna, known to elicit inflammatory and necrotic changes in gill tissues in particular rearing conditions (Mladineo, Reference Mladineo2006; Mladineo et al. Reference Mladineo, Šegvić and Petrić2011). Although nothing much is known of its epidemiological patterns in the final hosts, developed molecular tools might help to shed light on its life style. Therefore our goal was, using mitochondrial cytochrome c oxidase subunit I locus, to asses: (i) genetic structure of two D. katsuwonicola populations from bluefin tuna in the Adriatic Sea, isolated temporally one from the other (2008 and 2013) and (ii) origins of individuals sharing a single cyst from the parasite population sampled in 2013 by analysing their haplotype diversity.
MATERIALS AND METHODS
Parasite and fish sampling
Atlantic bluefun tuna were reared in a facility at WS part of Island of Brač (43°05′67″N, 15°28′63″E) Adriatic Sea. Tuna gill arches (n = 25) that showed the presence of didymozoid cysts were sampled during the harvest in February 2013 and processed as described in Mladineo et al. (Reference Mladineo, Šegvić and Petrić2011) to isolate encysted individuals of D. katsuwonicola (n = 56). Each isolated individual was inspected under the microscope to determine its maturity stage (juvenile/adult, based on the presence of eggs) and fixed in absolute ethyl alcohol for genetic characterization. Special care was taken to avoid mixing individuals that originated from a single cyst.
Additionally, sequences of D. katsuwonicola (n = 25) collected from the same site in 2008 were downloaded from National Center for Biotechnology Information (NCBI) for further population analysis (GenBank accession numbers: KF379712–KF379736).
DNA isolation, cytochrome c oxidase subunit 1 gene (COXI) fragment amplification, sequencing and alignment
The same methodology as described in Mladineo et al. (Reference Mladineo, Bott, Nowak and Block2010) was employed to obtain a 745 bp partial mitochondrial COX1 gene of D. katsuwonicola. Obtained GenBank accession numbers are KF712395–KF712450.
Genetic diversity and population structure of D. katsuwonicola
Nucleotide sequences of COX1 gene of D. katsuwonicola (in total n = 81) were compiled and aligned with Clustal X 1.83 (Thompson et al. Reference Thompson, Gibson, Plewniak, Jeanmougin and Higgins1997), using default parameters and further verified by GBlocks. Molecular diversity was measured using Dnasp 5.0 (Librado and Rózas, Reference Librado and Rózas2009) and Arlequin 3.5 (Excoffier and Lischer, Reference Excoffier and Lischer2010). Values for the number of haplotypes (H), polymorphic sites (S), haplotype diversity (h; Nei, Reference Nei1987), nucleotide diversity (p; Nei, Reference Nei1987) and the average numbers of pairwise nucleotide differences (k; Tajima, Reference Tajima1983) were estimated. Pairwise and overall distances among haplotype sequences were calculated in MEGA 5 (Tamura et al. Reference Tamura, Dudley, Nei and Kumar2007). A substitution model and the gamma distribution shape parameter for the rate of heterogeneity among sites were determined using Modeltest 3.07 (Posada and Crandall, Reference Posada and Crandall1998) based on the Hierarchical Likelihood Ratio Tests. For COX1 data, the Hasegawa, Kishino and Yano (HKY) model (Tamura and Nei, Reference Tamura and Nei1993) of evolution with the gamma shape parameter was selected for the analysis of molecular variance (AMOVA) and phylogenetic analysis. Populations were first constructed based on temporal isolation: population 2008 (sequenced deposited in GenBank KF379712–KF379736, previous study, n = 25) and population 2013 (present study, n = 56). Further, to determine haplotype relations and origins of individuals within a single cyst, didymozoids sampled in 2013 were divided randomly in two groups so that each individual belonged to different population (each population n = 28).
For the evaluation of hypothesized patterns of spatial genetic structure, a hierarchical AMOVA was used to partition variance components attributable to population variance and to individuals within the populations, where 10 000 permutations were performed to test the significance of pairwise population comparison. Pairwise genetic differentiation between populations was estimated using the fixation index F ST and statistical significances were tested with 10 000 permutations. Both AMOVA and F ST calculations were performed in Arlequin 3.5. For COX1 data, genetic distances between haplotypes were corrected using the HKY model of nucleotide substitution with a gamma shape parameter (G = 0·01). The null hypothesis of population panmixia was also tested in Arlequin 3.5 using an exact test of the differentiation of haplotypes among populations.
Tajima's D (Tajima, Reference Tajima1989) and Fu's Fs (Fu, Reference Fu1997) statistics were calculated to verify the null hypothesis of selective neutrality in relation to mtDNA sequences, which would be expected with population expansion. Mismatch distributions (Harpending, Reference Harpending1994) were constructed using Arlequin 3.5. The shapes of the mismatch distributions were used to deduce whether a population has undergone sudden population expansion (Rogers and Harpending, Reference Rogers and Harpending1992). The fit between the observed and expected distribution was tested using the Harpending raggedness index (HRI; Harpending, Reference Harpending1994) and sum of squared deviation (SSD) for the estimated stepwise expansion models (Schneider and Excoffier, Reference Schneider and Excoffier1999), implemented in Arlequin 3.5. Significance was assessed on the parameters with permutation tests under the null hypothesis that sudden population expansion cannot be rejected.
Graphic presentation of genetic diversity and population structure was inferred activating R statistics (R.2.15.2) command in the Arelquin 3.5.
Phylogenetic analysis
Phylogenetic relationships and distributions of didymozoid haplotypes in the Adriatic were reconstructed with median-joining networks of mutations (Bandelt et al. Reference Bandelt, Forster and Rohl1999) using Network v 4.5.1.6 software (available at: http://www.fluxus-engineering.com/sharenet.htm). Bayesian inference (BI; Larget and Simon, Reference Larget and Simon1999) analysis was performed in MRBayes v3.1.2 (Huelsenbeck and Ronquist, Reference Huelsenbeck and Ronquist2001) using a HKY+G (Tamura and Nei, Reference Tamura and Nei1993) evolutionary model of nucleotide substitution with a gamma parameter selected in jModelTest 0.1.1 (Posada, Reference Posada2008) using datasets from 2008 and 2013. Four incremented ‘heated’ Markov chains were carried through 2 500 000 generations, while sampling every hundredth generation, and 6250 samples were discarded as burnin. Markov chain Monte Carlo parameters were calculated with default properties. A consensus tree with the 50% majority rule was constructed from the tree output file produced in the BI analysis and visualized using FigTree (http://tree.bio.ed.ac.uk/software/figtree/).
RESULTS
On average, sampled tuna (2013, n = 25) had 1·63 cysts, ranging from one to maximum seven cysts counted on the individual tuna. From total number of 37 isolated cysts, most of encysted individuals were mature specimens (n = 70). Only in one cyst, a single juvenile without sibling occupant was observed, while in another cyst, the pair consisted of one mature and one immature specimen.
Genetic structure of temporally distant populations (2008 and 2013) of D. katsuwonicola
Sequence divergence (Tamura and Nei distance) among COX1 haplotypes sampled in 2008 and compared with 2013 data ranged from 0 to 2·8%, with an average of 0·3%.
In total, 50 variable sites were observed and 39 haplotypes were detected for the COX1 fragment. Sequence divergence (Tamura and Nei distance) among COX1 haplotypes ranged from 0 to 2·3%, with an average of 0·4%. Among 50 polymorphic sites observed in COX1 fragment, 35 were singleton variable sites and 15 were parsimony informative. Among 39 haplotypes defined, most (31, 79·5%) were unique and represented by a single individual. Only one haplotype was shared among two designated populations (H1), which was also the most abundant haplotype present in two populations. Relative frequency of H1 in 2008 and 2013 didymozoid population was 76 and 25%, respectively. Genetic diversity indices for each population are summarized in Table 1, indicating high levels of haplotype diversity and low nucleotide diversity. As well, COX1 displayed high value of average number of nucleotide differences k.
Table 1. Descriptive statistics of genetic diversity of populations of D. katsuwonicola of T. thynnus caught in the Adriatic Sea in the 2008 and 2013 based on COX1 sequence data
N, sample size; H, number of haplotypes; S, number of polymorphic loci; h, haplotype diversity (±s.d.); π, nucleotide diversity (±s.d.); k, mean pairwise difference (±s.d.).
Genetic structure of individual D. katsuwonicola within a single cyst
Unambiguous sequences for the mitochondrial COX1 of D. katsuwonicola population sampled in 2013 were obtained and deposited in GenBank (KF712395–KF712450). The length of the COX1 in the analysis was 727 bp (n = 56). In total, 41 variable sites were observed and 34 haplotypes were detected for the COX1 fragment. Sequence divergence (Tamura and Nei distance) among COX1 haplotypes ranged from 0 to 2·3%, with an average of 0·4%. Among 41 polymorphic sites observed in COX1 fragment, 27 were singleton variable sites and 14 were parsimony informative. Among 34 haplotypes defined, most (28, 79·4%) were unique and represented by a single individual. Only six haplotypes were shared among the two designated populations (H1, H4, H10, H25, H27 and H34) and H1 was the most abundant haplotype present in two populations (Supplementary Table S1). Genetic diversity indices for each population are summarized in Table 2, indicating high levels of haplotype diversity and low nucleotide diversity. As well, COX1 displayed high value of average number of nucleotide differences k.
Table 2. Descriptive statistics of genetic diversity of populations of D. katsuwonicola of T. thynnus caught in the Adriatic Sea in 2013 based on COX1 sequence data
N, sample size; H, number of haplotypes; S, number of polymorphic loci; h, haplotype diversity (±s.d.); π, nucleotide diversity (±s.d.); k, mean pairwise difference (±s.d.).
Population genetic structure
Genetic differentiation among populations (2008 vs 2013) was assessed using F ST pairwise comparison. Global F ST values were very low, showing no significant genetic structure in the investigated populations (F ST (COX1) = 0·056 72, P = 0·0000). Same pattern of genetic structure was confirmed by AMOVA, which attributed 5·67% of the genetic variation to variability among populations, and 94·33% variation within populations (Table 3). Over all samples, non-differentiation exact P-values were not significant (P = 0·263 90 ± 0·040 65), not rejecting that the temporal population of D. katsuwonicola in Adriatic T. thynnus is panmictic.
Table 3. AMOVA for temporary distant (2008 and 2013) populations of D. katsuwonicola, based on mitochondrial COX1 sequences data
When we have randomly divided didymozoids inhabiting a single cyst collected in 2013 in two groups, similarly, F ST pairwise comparison showed low values (F ST (COX1) = −0·013 53, P = 0·0000) and AMOVA attributed −1·35% of the genetic variation to variability among populations, and 101·35% variation within populations (Table 4). Over all samples, non-differentiation exact P-values were not significant (P = 0·963 96 ± 0·0142), not rejecting that the population of D. katsuwonicola in Adriatic T. thynnus is panmictic.
Table 4. AMOVA for D. katsuwonicola populations, based on mitochondrial COX1 sequences data
Phylogenetic and network analyses
The topology of trees built from COX1 locus using BI was mostly unresolved in its larger part, while inferred sister clades formed groups with posterior probability ranging from 53 to 82%. Those clades included exclusively didymozoids collected in 2013 that belonged to same haplotypes (Fig. 1). This situation was reflected in the haplotype network that showed a star-like phylogeny with most of the unique haplotypes closely related to the common central haplotype (H1) (Fig. 2). The later was composed of 19 and 14 D. katsuwonicola collected in 2008 and 2013, respectively. If parasite underwent expansion, the central common haplotype was likely the ancestral one, from which stepwise nucleotide substitutions could explain all other unique haplotypes. The most divergent haplotype (H16) that developed after 668 mutations belongs to didymozoid collected in 2013 (Dk26 in Fig. 1).
Fig. 1. Unrooted phylogenetic tree inferred by Bayesian analysis of mitochondrial cytochrome c oxidase subunit 1 (COX1) gene locus fragments of D. katsuwonicola isolated from gills of the Atlantic bluefin tuna (T. thynnus) in 2008 (taxa DAd) and 2013 (taxa Dk). Posterior probability values are shown in different colours, from red (53%) to blue (82%), as represented on the scale bar (left).
Fig. 2. Phylogenetic network of 39 haplotypes from D. katsuwonicola infecting gills of the Atlantic bluefin tuna (T. thynnus). Network lines mark relationships between individual haplotypes (circles). Circle sizes correspond to numbers of sequences that belong to certain haplotypes. Colours of circles represent haplotypes collected in different years: 2008 (black) and 2013 (white), respectively. Light grey nodes represent hypothetical haplotypes necessary for emergence of sampled haplotypes. // represents a cut in hypothetical line of haplotype H16 mutation steps, indicating that 668 mutation from the central haplotype H1 took place in order to generate H16.
Demographic patterns
The demographic history of D. katsuwonicola was investigated using mismatch distribution. The goodness of fit test showed that no mismatch distributions for temporal datasets deviated significantly (P > 0·05) from predicted values under the sudden expansion model, thereby providing further evidence for population expansion in the past. Furthermore, the HRI values were low for both the designated didymozoid populations, indicating a significant fit between the observed and the expected distribution, thereby providing further evidence for population expansion in the past. SSD also supported a non-departure from the null hypothesis of population expansion. The θ-values for all cases confirmed the scenario of population expansions; the initial θ-values were always much smaller than the final θ-values (Supplementary Table S2). Some departure of the age expansion parameter (τ) between populations from different years suggest that the population expansion may date back to close historical period.
The results of Tajima's D test and Fu's Fs test are presented in Supplementary Table S2. For didymozoids Tajima's D values were negative for the Adriatic and Mediterranean, as well as overall populations, indicating an excess of rare nucleotide site variants compared with the expectation under a neutral model of evolution, with statistically significant differentiation from the model of neutral evolution. Fu's Fs test, which is based on the distribution of haplotypes, also showed negative values for designated and overall populations, indicating an excess of rare haplotypes over what would be expected under neutrality. Following these tests, the hypothesis of neutral evolution was rejected for both populations.
DISCUSSION
Analysing genetic structure of two D. katsuwonicola populations (2008 vs 2013) and further breaking down 2013 population in two randomly assigned groups (assigning each of two fellow occupants in the opposite group), we observed the existence of a single, genetically unstructured and panmictic population. Genetic structuring of parasite populations is usually attributed to parasite and host ecological characteristics as wide distributional range, fragmented nature of the habitat, low expected rate of long distance dispersal or hermaphrodite-type reproduction that fixes different alleles in different subpopulations, making them more structured (Blouin et al. Reference Blouin, Yowell, Courtney and Dame1995; Mes, Reference Mes2003). In case of didymozoid digeneans, our knowledge in terms of listed traits has been limited by their complex indirect life cycle and adulthood in large pelagic fish, resulting that even the development within the tissue cyst and dispersion of embryonated eggs remains speculative. On the temporal scale divided by 5-year period, our data suggest a constant gene flow that resulted in dominance of a single haplotype attributing to approximately 40% of sampled taxa, and the rest belonging to unique haplotypes, most likely originating from this central one. Comparing haplotype and nucleotide diversity between samples of two collecting years, 2008 seems to be less rich, indicating a relatively small number of didymozoid individuals that contributed to the offspring hatching. However, overall haplotype diversity is high compared with the overall nucleotide diversity, which can be a signature of a rapid demographic expansion from a small effective population size (Avise, Reference Avise2000), further supported by Tajima's D and Fu tests. Both 2008 and 2013 populations had a negative value indicative of a bias towards rare alleles, which is a landmark of recent population expansion and recent colonization event in the Adriatic. Observed haplotypes divergence (COXI 0·1–11·8%, average of 0·6%) is mostly in accordance with previous data for trematode species: Haliotrema aurigae (0·17–0·7%), Euryhaliotrematoides grandis (0·17–7·99%), Polylabroides (0·12–1·89%), Microcotyle (0·14–2·56%) (Li et al. Reference Li, Shi, Brown and Yang2011) and Clinostomum sp. metacercariae (Bonett et al. Reference Bonett, Steffen, Trujano-Alvarez, Martin, Bursey and McAllister2011).
When separating randomly hermaphroditic fellow occupants in two groups, we observed that only a negligible number of individuals (16·67%) within the same cysts belong to the same haplotype (H1, H25, H27 and H34). This fortifies the hypothesis that self-fertilization is less frequent type of reproductive strategy, but still does occur. In contrast, in 83·33% of cases different individuals contributed to the offspring production supporting the hypothesis of cross-fertilization between encysted didymozoids as more frequent reproduction type. Most paired individuals (83·33%) originated from different parent, belonging to a haplotype different from the haplotype of the second individual within the cyst. Such admixture of self- and cross-fertilized offsprings has been evidenced previously although in a cestodian model. In vitro studies evidenced that such strategy enables a higher infection success and faster development of cross-fertilized progeny when faced by a competitive situation where superior abilities of parasite to extract limiting resources from the host were crucial (Christen et al. Reference Christen, Kurtz and Milinski2002). Another in vitro study using Echinostoma caproni, a digenean that breeds in mammals or birds, showed that the percentage of self- and cross-fertilization was equalized (Trouvé et al. Reference Trouvé, Renaud, Durand and Jourdane1996), contrary to D. katsuwonicola. Authors suggest that selective pressure induced the type of reproduction observed, which indeed might explain why in vivo as in our study, self-fertilization took a minor part.
Reaching the definitive host after passing through the whole cycle is considered a relatively rare event, and further co-occurrence and physical encounter in the same host, at the same site and time of two didymozoids that will encapsulate together, is even more unusual (Poulin, Reference Poulin1998). To overpass such selective pressure, some didymozoids developed hermaphroditism as an ancestral reproductive trait, although few other didymozoid subfamilies evolved to more or less pronounced gonochorism (Pozdnyakov and Gibson, Reference Pozdnyakov, Gibson, Bray, Gibson and Jones2008). This indicates a tendency to switch towards dioecy, reinforced by the evidence that although a hermaphrodite, D. katsuwonicola still continues to encyst in pairs and cross-fertilize. Except inducing wider genetic heterogeneity in population that helps to overpass accumulation of deleterious recessive mutations (Maynard Smith, Reference Maynard Smith1989), dioecy has been suggested beneficial also because it implies a division of labour between sexes, consequently leading to higher fitness of both individuals (Poulin, Reference Poulin1998). Another example of dioecy in Digenea that derived from hermaphroditism is the case of Schistosoma spp. (Basch, Reference Basch1990). However, Schistosomatium douthitti apparently uses parthenogenesis regularly, even when paired with conspecific males (Jourdane et al. Reference Jourdane, Imbert-Establet and Tchuem Tchuenté1995). Such switch towards one or another reproductive strategy only reflects that at the certain evolutionary time, conditions favoured one towards the other strategy in the parasite life cycle.
Haplotype diversity in Adriatic didymozoid is indicative of a relatively large number of didymozoid individuals that contributed to the offspring hatching and suggests a large effective population size. Presence of multiple and diverse haplotypes within Adriatic didymozoid population and within a single cyst, suggests that cross-fertilization consequently leads to wider dispersion and fitness enabling panmictic distribution of the parasite. In most cases, we have observed two mature encysted individuals in respect to only two cases when an individual from the pair was immature, or the only one present. Lester (Reference Lester1980) observed in non-cyst forming didymozoids (Neometadidymozoon helicis, Nematobothrium spinneri and Didymozoon brevicolle) that most individuals were missing from their tissue cavity and instead a mass of eggs within worm remains was left over. Similarly, Timon-David (Reference Timon-David1937) noted haemorrhages in the gills surrounding D. wedli (syn. D. katsuwonicola), suggesting that the cyst reached a certain degree of maturity than spontaneously opened and evacuated its contents. Both authors conclude that the release of eggs through ulceration may be widespread among Didymozoidae. This was not confirmed in our study as no degraded adults or free egg masses were observed, as well as no cyst opening that would suggest the death of one of the paired didymozoid.
In conclusion, we evidenced that both self- and cross-fertilization mechanisms are employed by two mature didymozoids within the cyst, the later occurring more often, consequently stimulating higher gene flow and unstructured population genetics. Mostly, two adults occupy the single cyst, but on rare occasion immature juvenile may be alone in the cyst or coupled with a mature individual. Such traits result in a successful propagation of progeny that shows heterogeneous and panmictic population in the Adriatic, even though they develop through three or four hosts and undergo a very complex migration from the digestive system to periphery.
SUPPLEMENTARY MATERIAL
To view supplementary material for this article, please visit http://dx.doi.org/ 10.1017/S0031182015000839.
FINANCIAL SUPPORT
This work was supported by the Student Programme Grant, University of Split, Croatia to M. T.
