Evidence for cryptic speciation in directly transmitted gyrodactylid parasites of Trinidadian guppies.

Cryptic species complexes are common among parasites, which tend to have large populations and are subject to rapid evolution. Such complexes may arise through host-parasite co-evolution and/or host switching. For parasites that reproduce directly on their host, there might be increased opportunities for sympatric speciation, either by exploiting different hosts or different micro-habitats within the same host. The genus Gyrodactylus is a specious group of viviparous monogeneans. These ectoparasites transfer between teleosts during social contact and cause significant host mortality. Their impact on the guppy (Poecilia reticulata), an iconic evolutionary and ecological model species, is well established and yet the population genetics and phylogenetics of these parasites remains understudied. Using mtDNA sequencing of the host and its parasites, we provide evidence of cryptic speciation in Gyrodactylus bullatarudis, G. poeciliae and G. turnbulli. For the COII gene, genetic divergence of lineages within each parasite species ranged between 5.7 and 17.2%, which is typical of the divergence observed between described species in this genus. Different lineages of G. turnbulli and G. poeciliae appear geographically isolated, which could imply allopatric speciation. In addition, for G. poeciliae, co-evolution with a different host species cannot be discarded due to its host range. This parasite was originally described on P. caucana, but for the first time here it is also recorded on the guppy. The two cryptic lineages of G. bullatarudis showed considerable geographic overlap. G. bullatarudis has a known wide host range and it can also utilize a killifish (Anablepsoides hartii) as a temporary host. This killifish is capable of migrating overland and it could act as a transmission vector between otherwise isolated populations. Additional genetic markers are needed to confirm the presence of these cryptic Gyrodactylus species complexes, potentially leading to more in-depth genetic, ecological and evolutionary analyses on this multi-host-parasite system.

Introduction

Comparative phylogeographic studies have demonstrated co-evolution between parasites and their hosts [1, 2], yet co-speciation does not appear to be a major factor in host-parasite evolution [3]. Many confounding factors such as multi-host systems, extinction and host-switching can cause incongruence between host and parasite phylogenies [4, 5]. Parasites typically evolve faster than their hosts [6-8] related to their shorter generation times that may result in higher substitution rates [9]. Although hosts and parasites may be involved in a co-evolutionary arms race in which adaptive evolution is driven by natural selection, parasite transmission often results in serial population bottlenecks and founder events, which can have a significant role in driving parasite evolution [6-8]. Such demographic dynamics will leave a distinct phylogenetic and population genetic signature, possibly leading to an increase in the genetic differentiation of parapatric/allopatric parasite populations, which can ultimately result in insipient speciation [10].

Monogeneans are particularly suitable parasites for revealing novel insights into host ecology and evolution [11, 12], due to their direct life-cycles and relatively high level of host-preference. Species from the genus Gyrodactylus belong to the most studied group of monogeneans, renowned for their impact on aquaculture and conservation planning. They parasitize many wild and farmed fish, including Atlantic salmon (Salmo salar) and brown trout (Salmo trutta) [13], and because they reproduce in situ on the host and are transmitted during host contact, epidemics can sweep quickly through fish populations. Although more than 400 species have been described for the genus, Gyrodactylus species are morphologically conserved and this is probably related to the adaptations associated with viviparity and progenesis [14, 15]. For this reason, molecular data is often essential for confirmation of species descriptions [16, 17] and discrimination [15-19]. Genetic studies have uncovered otherwise cryptic species [20] and demonstrated that Gyrodactylus diversity may be underestimated [17]. This is likely to be the case also for the well-studied guppy-gyrodactylid system in Trinidad [21].

Poecilia reticulata (guppy)—Gyrodactylus bullatarudis/G. turnbulli associations from Trinidadian rivers have been extensively examined in order to understand how parasitism can impact a host, namely affecting body size [22], reproduction [23-25], immune response and survival [26-28]. Trinidadian guppy-gyrodactylid dyads have also been studied to evaluate parasite transmission dynamics in social hosts [29, 30] and to test the impacts of infections on conservation strategies involving reintroduction from captive-breeding programs [31, 32]. Although the two parasite species can infect the same host species, they are very distinct phylogenetically [15, 33]. G. bullatarudis can also use a broader host range than G. turnbulli, both under laboratory [34, 35] and field [36] conditions. Relatively little is known about the biology of other poeciliid gyrodactylids, including G. poeciliae, which is closely related to G. bullatarudis [15]. Despite the obvious importance of this group of parasites to the evolutionary ecology of their hosts, very little is known about their molecular ecology. In fact the only record present in GenBank for G. bullatarudis, G. poeciliae and G. turnbulli (four sequences in total) were generated to confirm specific status and address phylogenetic questions [15, 21, 37].

In contrast, the genetic diversity and differentiation of guppy populations from Northern Trinidad has been the focus of many studies over the last three decades; based on this knowledge we can make specific predictions about the genetic diversity of their parasite species. Originally, based on allozyme and mitochondrial data, Shaw et al. [38] and Fajen and Breden [39] argued that guppies from Northern Trinidad form at least two genetically distinct groups, possibly originating from two different colonisation events of the Northern-Caroni and Oropouche drainages, and from different ancestral stocks. Recent studies, however, using microsatellites and SNPs indicate that colonisation of Trinidadian rivers was much more complex [40, 41]. In general, all studies are in agreement in detecting high levels of genetic differentiation between populations from different drainages due to genetic drift and adaptive selection [41]. Population differentiation has also been detected between upland and lowland sites in the Marianne River (Northern drainage) and rivers from the Caroni drainage, related to genetic drift caused by geographic distance and physical barriers such as waterfalls [42, 43]. Guppy migration within drainages mainly occurs in a downstream direction with headwater populations generally being more isolated and less genetically diverse [42, 43]. Migration does occur between rivers, but typically in lowland stretches during wet season flooding, as demonstrated by high levels of genetic diversity of lowland populations [38, 43].

We hypothesized that, due to the high level of phylogeographic structure previously reported amongst host populations and the high levels of host-specificity reported for Gyrodactylus species, parasite speciation might have led to multiple cryptic species. Furthermore, as G. bullatarudis uses alternative hosts and has a higher potential for host switching and migration, it is likely to show a lower level of genetic differentiation between rivers when compared to more host-specific Gyrodactylus species. Thus, the present study aimed to examine the genetic diversity and phylogeographic patterns of Gyrodactylus species infecting guppies (P. reticulata) from five rivers of North Trinidad using partial sequences of the mitochondrial cytochrome oxidase II (COII) gene. In addition, the phylogeographic structure of guppies was evaluated through sequences of a portion of the first domain (HVR) of mitochondrial control region (D-loop) to compare with that of the parasites and to test the levels of host-specificity. Our results seem to confirm both our hypotheses and provide the first evidence that cryptic speciation might have occurred within the three Gyrodactylus spp. studied.

Material and Methods

Ethics Statement

This work was conducted using the guppy (Poecilia reticulata)–Gyrodactylus turnbulli/ G. bullatarudis/ G. poeciliae model systems. All fish were collected, handled and killed [44] according to UK Home Office Project license (PPL 30/2357) regulations and approved by the Cardiff University Ethics Committee. In Trinidad there is no legislation restricting collection of fishes from public areas (and none of our sites were located or accessed via privately owned land) and P. reticulata is not an endangered or protected species. Map reference coordinates for sampled sites are as follows: Lower Aripo 694410E, 1177783N; Upper Aripo 694030E, 1182128N; Lower Lopinot 683553E, 1175663N; Upper Lopinot 683520E, 1182443N; Lower Marianne 685890E, 1193642N; Upper Marianne 685891E, 1192747N; Lower Oropouche 0704394E, 1178967N; Upper Oropouche 0700467E, 1183194N; and Upper Yarra 683427E, 1189518N.

Fish collection and morphological identification of Gyrodactylus spp.

Between October 2003 and November 2004, 230 specimens of P. reticulata were collected from 9 different locations of the Caroni (Lower and Upper Aripo River and Lower and Upper Lopinot River), Northern (Upper Yarra, Upper and Lower Marianne River) and the Oropouche drainages (Upper and Lower Oropouche), Trinidad, and stored individually in eppendorfs pre-filled with 90% ethanol for subsequent molecular analysis (see Table 1 for details, and Fig. 1 in Willing et al. [41] for drainages and Fig. 1 in Barson et al. [43] for details of Caroni rivers). The samples were transported back to Cardiff and specimens of Gyrodactylus were collected from the fins and body of guppies, individually placed in a drop of water on a glass slide and bisected in half using fine needles under a dissecting-microscope. The posterior half of the worm was partially digested using the method detailed by Paladini et al. [45], and the specimens initially identified according to Harris & Cable [21] and Cable et al. [37]. A total of 131 Gyrodactylus were identified to species level in this way and the anterior half of each parasite was kept for molecular analysis in 90% ethanol (Table 1). Our initial plan was to obtain rDNA Internal Transcribed Spacers (ITS) and microsatellites data from each individual parasite, but due to difficulties encountered with the mtDNA sequencing, which required multiple PCRs per worm due to high PCR failure rates from individual worms, this was not possible given the limited amount of DNA remaining. For some specimens, DNA was extracted only from the anterior part of the parasite, measuring less than 0.2 mm and containing fewer than 200 cells [15, 46], sufficient in theory for three PCRs. If this failed, DNA was extracted from the remainder of the parasite.

Table 1

Number of individuals sampled (N), number of haplotypes (H), Haplotype diversity (Hd) and nucleotide diversity (π) per locality for each lineage within Gyrodactylus morphospecies G. turnbulli (GtM, GtC and GtO), G. poeciliae (GpY and GpCM), and G. bullatarudis (Gb1 and Gb2), and for their fish host, Poecilia reticulata.

Morphospecies	Lineage		ARIPO	Upper Aripo	Lower Aripo	LOPINOT	Upper Lopinot	Lower Lopinot	OROPOUCHE	Upper Oroupoche	Lower Oroupoche	MARIANNE	Upper Marianne	Lower Marianne	UPPER YARRA	Total
G. turnbulli	GtC	N	7	6	1	5	2	3	-	-	-	-	-	-	-	12
		H	2	1	1	3	2	1	-	-	-	-	-	-	-	4
		Hd	0.29	0	-	0.7	0.5	0	-	-	-	-	-	-	-	0.731
		π	0.006	0	-	0.011	0.004	0	-	-	-	-	-	-	-	0.024
	GtO	N	-	-	-	-	-	-	1	-	1	-	-	-	-	1
		H	-	-	-	-	-	-	1	-	1	-	-	-	-	1
		Hd	-	-	-	-	-	-	-	-	-	-	-	-	-	-
		π	-	-	-	-	-	-	-	-	-	-	-	-	-	-
	GtM	N	-	-	-	-	-	-	-	-	-	12	-	12	-	12
		H	-	-	-	-	-	-	-	-	-	2	-	2	-	2
		Hd	-	-	-	-	-	-	-	-		0.409		0.409	0.409	0.409
		π	-	-	-	-	-	-	-	-	-	0.000	-	0.000	-	0.000
G. poeciliae	GpCM	N	12	3	9	6	5	1	-	-	-	40	17	25	-	60
		H	2	1	1	3	2	1	-	-	-	3	1	3	-	6
		Hd	0.41	0	0	0.6	0.4	-	-	-	-	0.5	0	0.24	-	0.485
		π	0.006	0	0	0.008	0.002	-	-	-	-	0.002	0	0.003	-	0.007
	GpY	N	-	-	-	-	-	-	-	-	-	-	-	-	2	2
		H	-	-	-	-	-	-	-	-	-	-	-	-	1	1
		Hd	-	-	-	-	-	-	-	-	-	-	-	-	0	0
		π	-	-	-	-	-	-	-	-	-	-	-	-	0	0
G. bullatarudis	Gb1	N	16	5	11	5	-	5	11	-	11	-	-			32
		H	2	1	2	1	-	1	1	-	1	-	-	-	-	2
		Hd	0.125	0	0.182	0	-	0	0	-	0	-	-	-	-	0.315
		π	0.001	0	0.007	0	-	0	0	-	0	-	-	-	-	0.001
	Gb2	N	5	5	-	1	-	1	6	6	-	-	-	-	-	12
		H	1	1	-	1	-	1	1	1	-	-	-	-	-	2
		Hd	0	0	-	-	-	-	0	0	-	-	-	-	-	0.167
		π	0	0	-	-	-	n.a	0	0	-	-	-	-	-	0.001
P. reticulata		N	61	36	25	57	25	32	50	24	26	62	27	35	-	230
		H	6	5	3	7	3	5	4	3	4	7	5	5	-	22
		Hd	0.60	0.39	0.35	0.79	0.64	0.69	0.76	0.68	0.75	0.60	0.64	0.55	-	0.905
		π	0.003	0.001	0.004	0.023	0.002	0,019	0.003	0.003	0.003	0.004	0.004	0.003	-	0.023

Figure 1

Median joining haplotype networks estimated for (A) Poecilia reticulata, and (B-D) for each Gyrodactylus species from nine localities in the Northern Mountain Range of Trinidad.

DNA extraction, PCR amplification and sequencing

For P. reticulata, genomic DNA was extracted from the caudal fin using the HotSHOT protocol [47]. A total of 414 bp of the first domain (HVR) of the mitochondrial (mt) D-loop was amplified using primer pairs L15926 and H16498 [48, 49] according to the following PCR conditions: 95°C for 3 minutes and 35 cycles at 95°C for 30 seconds, 54°C for 30 seconds and 72°C for 1 minute. For Gyrodactylus spp. genomic DNA was extracted using DNeasy Blood & Tissue QIAGEN kit. A 262 bp fragment of the mt protein coding COII gene was amplified using the QIAGEN Multiplex PCR kit and primer pairs COX2 F1 (TACATAYCGCCCGTCAATYT) and COX2 R1 (TCARTAYCACTGDCGDCCYA), developed specifically for this study, and following PCR protocol of initial denaturation temperature of 95°C for 15 minutes and 40 cycles at 95°C for 30 seconds, 50°C for 90 seconds and 72°C for 90 seconds. All PCR products were cleaned using Exo/AP and sequenced by a commercial company. All sequences were deposited in GenBank (Accession Numbers KP168263-KP168415).

Phylogenetic analysis

The COII alignment was translated into protein to examine the presence of stop codons using DnaSP software v5.10 [50]. Median-Joining haplotype networks were reconstructed for the guppy host and separately for each Gyrodactylus spp. using Network 4.611 software [51]. Maximum Likelihood and Bayesian Inference phylogenetic tree reconstructions were estimated for parasites using the software Garli 0.96 [52] and MrBayes v.2.1 [53], respectively. G. salaris and G. thymalli (GenBank acc. nos. NC008815 and NC009682, respectively) were used to root the Gyrodactylus spp. tree. Appropriate model of sequence evolution was chosen using jModeltest2 and the AIC criteria [54], and the GTR+G (-lnL = 1451) model was selected for phylogenetic tree reconstruction. Finally, to assess intraspecific and interspecific divergences, uncorrected p-distances were calculated within Gyrodactylus spp. and between each major clade using the software MEGA version 5 [55]. To assess typical levels of intraspecific diversity in Gyrodactylus spp., uncorrected p-distance was additionally calculated for the same 262 bp fragment for all published COII from other Gyrodactylus spp. available on GenBank (acc. nos. GU131204; GU131200; GU131198; GU131220, GU131214; GU131210; EU293891; NC009682 and NC008815). Substitution saturation was tested and rejected for the three current Gyrodactylus spp. datasets (results not shown) using a test by Xia et al. [56, 57] implemented in software DAMBE 5.2.78 [58]. For guppies, uncorrected p-distances were calculated using the option to treat gaps as pairwise deletions.

Population structure analysis

Haplotype and nucleotide diversities were calculated for the guppy and main phylogenetic lineages of Gyrodactylus spp. for each sampling locality whenever n ≥ 3 in DnaSP software v5.10 [50]. F measures were calculated for guppy and widespread lineages of Gyrodactylus spp. (see Results section) for each sampling locality whenever n ≥ 3 using Arlequin v3.5.1.3 [59].

Assessment of host specificity levels

To test for any correlation between genetic differentiation of host and parasite populations, F values for each of the selected lineages of Gyrodactylus spp. were plotted against those of the guppy host and Pearson correlation coefficient was calculated for each parasite-host dyad in Minitab ver. 12.1.

Morphometric analysis of Gyrodactylus specimens

Following results from the phylogenetic analyses, 25 individual point-to-point morphometric characteristics of the hamuli, ventral bar and marginal hooks (as outlined by [60]) were measured from those specimens clearly identified within a particular genetic lineage, and for which the posterior half of the body was of suitable quality. The aim of the morphometric analysis was to test whether these taxa were truly morphologically cryptic. To visualise multivariate patterns of morphological variation, a Principal Component Analysis (PCA) was conducted. To test whether genetic lineages where morphologically differentiated, a multivariate ANOVA based on permutations of the Euclidean Distance matrix between lineages mean values was performed. Post-hoc pairwise comparisons were also used to test whether there were significant morphological differences between cryptic lineages of G. bullatarudis (Gb1 and Gb2), as samples sizes were too small to test the same for the different lineages of G. turnbulli and G. poeciliae. All analyses were performed using available statistical R packages.

Results

Morphological identification of Gyrodactylus specimens

Morphological analyses of 131 Gyrodactylus specimens allowed the identification of three species: G. poeciliae, G. turnbulli and G. bullatarudis. While the latter two species had been already reported to parasitize P. reticulata, this is the first record of G. poeciliae infecting this host species.

Phylogenetic analysis and macroevolutionary patterns of host and parasite species

Phylogenetic analysis of the guppy shows high phylogeographic structure of populations (see Fig. 1A). Three haplogroups are evident in haplotype network, which include: 1) most individuals from Aripo and 23 individuals from Lower Lopinot River (Caroni drainage); 2) individuals from Marianne River; and 3) individuals from the Oropouche, the remaining individuals from Lopinot and one individual from the Aripo River. Divergence amongst guppy individuals varied between 0 and 4.4%.

Although our phylogenetic tree supports the monophyly of the three Gyrodactylus taxa (Fig. 2), divergence within G. turnbulli, G. poeciliae and G. bullatarudis varied between 0.4–17%, 0.4–6.9% and 0–13%, respectively. Since the values found for typical interspecific divergence of Gyrodactylus spp. range between 4% (G. salaris vs G. thymalli) and 39% (e.g. G. thymalli vs G. corydori), the intra-specific divergence found in the present study is equivalent to the typical interspecific divergence in all three cases, suggesting the presence of cryptic lineages.

Figure 2

Best Maximum-likelihood tree obtained for the COII of the three Gyrodactylus study species.

Values of nodes correspond to bootstrap support and Bayesian posterior probability, respectively. The three major clades corresponding to Gyrodactylus morphospecies are surrounded by a dotted line. Cryptic lineages within these clades are signaled with black bars.

For morphospecies G. bullatarudis, two phylogenetic lineages and haplogroups (Gb1 and Gb2) were detected (Figs. 1B and 2), although in this case, lineages/haplogroups overlapped in geographical range, both occurring in Aripo, Lopinot and Oropouche Rivers. Nonetheless, the two cryptic phylogenetic lineages from G. bullatarudis were retrieved as two very well supported sister clades (Fig. 2). Divergence within these lineages was low, ranging from 0–0.4%, and between lineages was reasonably high (12.6–13.0%).

For morphospecies G. turnbulli, three geographically segregated phylogenetic lineages and haplogroups were identified by both phylogenetic analyses (Figs. 1C and 2) distinguishing individuals from the Marianne River (GtM); Caroni drainage (GtC); and, an individual from the Oropouche River (GtO). Divergence within lineages varied between 0.4–3.4%, and between lineages was 5.7–7.3% between GtC vs GtO, 15.7–16.9% between GtC vs GtM and 16.9–17.2% between GtM vs GtO.

For morphospecies G. poeciliae, two apparently allopatric phylogenetic lineages and haplogroups were recovered by both phylogenetic analyses (Figs. 1D and 2), which distinguished two identical individuals collected from upper Yarra River (GpY) and the remaining individuals from the Caroni drainage and Marianne (GpCM). Divergence within GpCM varied between 0.4–2.3%, and divergence between the two lineages was 6.1–6.9%.

Population structure and microevolutionary analysis of host-parasite coevolution

Haplotype and nucleotide diversity for each Gyrodactylus lineage and guppy are depicted in Table 1. Haplotype diversity of morphospecies G. bullatarudis was generally low compared with the other two morphospecies. For the guppy, F values showed high genetic differentiation between populations from different rivers (S1 Table). Comparing the pairwise F values of parasite populations (S1 Table) to the pairwise F of host populations showed that the level of genetic differentiation was approximately similar for host and parasites (all pairwise t-tests: t<1.75 p≥0.141). Interestingly, there was a significant positive correlation between F values of G. poeciliae GpCM and guppies (Pearson correlation coefficient r = 0.479, p = 0.027), and a marginally significant correlation (Pearson correlation coefficient r = 0.775, p = 0.070) between F values of GtC and guppies (see Fig. 3), but not for Gb1 (Pearson correlation coefficient r = -0.230, p = 0.522).

Figure 3

F plots for guppy host vs parasites of (A) lineage Gb1 of Gyrodactylus bullatarudis, (B) lineage GpCM of G. poeciliae, and (C) lineage GtC of G. turnbulli.

Pearson’s correlation depicted if apparent. Additionally, mean F and SEM (standard error of the mean) are depicted in (D) for the guppy host and all 3 parasites.

Based on 25 measurements from the attachment hooks of each specimen (see S2 Table), PCA results highlighted high levels of morphological variation between Gyrodactylus lineages (GpCM, GtC, GtO, Gb1 and Gb2) examined, and PC1 and PC2 explained 42% and 18%, respectively, of the morphological variation found (Fig. 4) with PC3 and PC4 combined only contributing a further 16%. Differences between G. turnbulli (GtO and GtC) and the remaining Gyrodactylus species were clearly captured by PC1, whereas variation between GpCM, Gb1 and Gb2 were more evident along the PC2 axis, supporting the previous findings of Harris & Cable [21]. There was significant morphological differentiation between genetic lineages (ANOVA, F4,15 = 6.26, p<0.001).

Figure 4

PCA plot showing morphological variation along PC1 and PC2 axes for Gyrodactylus lineages GtC, GtO, GpCM and Gb1 and Gb2.

Discussion

Despite several studies having previously examined the genetic variability and structure of the guppy (Poecilia reticulata), the present dataset is the most comprehensive mitochondrial dataset of guppies in Trinidad. Our results confirm previous mtDNA studies and show that guppy population connectivity between different drainages is limited [38, 39]. There is evidence of migration from the Oropouche into the Caroni drainage, which has also been reported in recent SNP and microsatellite studies [40, 41]. Within rivers, gene flow appears to be in the downstream direction, with the lowlands populations of the Caroni drainage acting as a sink for immigrants from the upland habitats, which is consistent with previous studies [41, 43].

Regarding the guppy parasites, the three Gyrodactylus species examined each possess two to three genetic lineages, diverging by 5.7 to 17.2%. This is in the range of the divergence observed between published sequences for described species in this genus (4% divergence between G. salaris vs G. thymalli, and 39% between G. thymalli vs G. corydori [61-63]), which confirms that cryptic speciation is common within this genus [17]. In the case of G. bullatarudis, the two cryptic lineages had overlapping geographical ranges in both the Caroni and Oropouche drainages. Furthermore, for these two genetic lineages of G. bullatarudis there were no significant morphometric differences, confirming their cryptic nature. For G. turnbulli, the three cryptic lineages appeared to be geographically isolated in the Caroni, Oropouche and Marianne drainages. Within G. poeciliae, one lineage was found in the Yarra River (Northern drainage), and the other in both the Marianne River (Northern drainage) and Caroni drainage. However, while the divergence between lineages is evident in the haplotype network with both lineages being separated by a large number of mutations, phylogenetic support was low. Furthermore, as the G. poeciliae from the Yarra River were not included in the morphological analysis, the support for a cryptic species is much more tentative in G. poeciliae than for both other parasites. Remarkably, both for G. poeciliae and G. turnbulli, the level of genetic differentiation (expressed as pairwise population F) within lineages was comparable to that of the guppy host. Cryptic speciation and genetic diversification of the parasites may have occurred due to a variety of mechanisms including historical geographic isolation, co-evolution with multiple hosts and host switching, as will be discussed in the following section.

Macroevolution in Gyrodactylus spp. resulted in cryptic variation

Poulin and Morand [64] postulated that the main reason for the high level of diversification among monogenean parasites is their small body size, which confers the ability for adaptation to different host micro-habitats, leading to sympatric speciation on a host. Such sympatric speciation has indeed been confirmed in Dactylogyrus gill parasites [65]. In the current study, however we assessed the genetic variation of viviparous rather than oviparous monogeneans [18]. Gyrodactylus spp. are renowned for their “Russian-Doll” reproductive strategy with one embryo developing inside another within the mother’s uterus [16]. Unlike Dactylogyrus and other oviparous monogeneans, they lack a free-living dispersal stage and the parasite infrapopulation can grow exponentially on a host with transmission occurring during host-host contact [16]. Cryptic lineages were not found to be co-infecting the same fish, and in the case of G. turnbulli and G. poeciliae, lineages did not even overlap geographically, although both G. bullatarudis lineages were found in the same host populations (but never on the same individual). The current results therefore do not support the hypothesis that intra host niche partitioning has resulted in sympatric speciation of Gyrodactylus spp. on the host.

The data from G. turnbulli and G. poeciliae seem to be most consistent with allopatric speciation, given that their lineages appeared to be geographically isolated. However, as G. poeciliae was originally described from Poecilia caucana [21], and the current study is the first to describe this parasite species on the guppy (P. reticulata), it is possible that Trinidadian populations of G. poeciliae may have also coevolved (and hence differentiated) on multiple host species.

The two cryptic lineages of G. bullatarudis had overlapping geographical ranges in both the Caroni and Oropouche drainages. This morphospecies is known to infect multiple hosts under laboratory conditions [34], and in Trinidad at least one lineage of G. bullararudis is able to infect killifish (Anablepsoides hartii [syn. = Rivulus hartii]) [36]. Individuals of G. bullatarudis can survive on killifish out of water for over an hour [36], which may confer a dispersal advantage as these fish can migrate overland. Unlike the guppy, gene flow between G. bullatarudis populations is not restricted by natural barriers in rivers, such as waterfalls and weirs [66], because the parasite could (theoretically) migrate overland when infecting a killifish. Hence, although the two cryptic lineages could have arisen due to allopatric isolation on a single host, it is also possible that divergence arose from co-evolution with different hosts. Secondary contact of the two lineages may have been easily facilitated by infection of other fish, such as the killifish, which would explain their current geographic overlap. Whichever was the case, both hypotheses are consistent with our observation that unlike the other parasite species, there was no positive correlation between all pairwise population F values of guppies with those of the G. bullatarudis populations.

Conclusions

The current study represents the first evidence of cryptic genetic differentiation within Trinidadian Gyrodactylus species. We hypothesize that allopatric isolation and/or co-evolution with different hosts accounts for the extant species complexes. However, we fully acknowledge that such interpretations are only tentative at this stage because of the small sample sizes. Additionally, reported levels of differentiation should be confirmed with nuclear data to exclude the possibility of introgressive hybridisation [67]. The presence of multiple cryptic species could explain the difficulties previously encountered during development of microsatellite libraries for Gyrodactylus species of guppies [68]. The current findings will have important implications for future research using guppy-Gyrodactylus spp. dyads as models to test the impacts of parasites on host evolution [22, 69] and conservation planning [22, 31].

Pairwise F for (A) guppies and (B) lineages of Gyrodactylus poeciliae (GpCM), G. bullatarudis (Gb1 and Gb2), and G. turnbulli (GtC), respectively.

The symbol † corresponds to the only pairwise F estimated for lineage Gb2. Non-significant F values are highlighted in bold.

(DOCX)

Click here for additional data file.

Morphological measurements of Gyrodactylus bullatarudis, G. poeciliae and G. turnbulli (mean±1 standard deviation followed by the range in parentheses).

Measurements are provided in micrometres. Abbreviations: H = hamulus; VB = ventral bar; MH = marginal hook.

(DOCX)

Click here for additional data file.