Introduction
An understanding of the evolutionary processes that underpin host–parasite associations is a central aim in parasitology (Šimková et al., Reference Šimková2004). An early attempt at understanding these processes led to the development of the host–parasite co-speciation hypothesis, at the beginning of the 20th century (Fahrenholz, Reference Fahrenholz1913). According to this theory, each parasite species evolves with its host species and traces its evolutionary lineage, undergoing speciation whenever its host does (fig. 1a). Co-speciation is largely constrained by phylogenetic factors (Fahrenholz, Reference Fahrenholz1913). However, congruence between parasite and host phylogenies is rarely observed in practice (Hafner and Nadler, Reference Hafner and Nadler1990; Paterson and Banks, Reference Paterson and Banks2001; Mendlová et al., Reference Mendlová2012; de Vienne et al., Reference de Vienne2013) and is not static (Nylin et al., Reference Nylin2017; Jorge et al., Reference Jorge2018). Hence, co-speciation is not the sole process shaping host–parasite evolutionary histories. For instance, congruence between host and parasite phylogenies can be eroded by a variety of ecological or environmental features (Page, Reference Page and RDM2003; fig. 1) leading to duplication (intrahost speciation, fig. 1b); or lineage sorting (parasite lineage absent from a host lineage) as a consequence of the extinction of the parasite (failure to colonize any of the descendants of a population, fig. 1c) or “missing the boat” (failure to colonize all descendants of a host, fig. 1d); and/or by spreading (failure to speciate altogether, fig. 1e) (e.g. Banks and Paterson, Reference Banks and Paterson2005). Incongruence can also result from host-switching events (fig. 1f). Host-switching is a sudden colonization of a previously unused host species by a small founder population of a parasite species, and may lead to speciation, whereas the source population of the parasite remains unaffected by this event and continues exploiting its original host species (Rózsa et al., Reference Rózsa, Tryjanowski, Vas, Morand, Krasnov and Littlewood2015).
Host–parasite co-evolutionary scenarios: (a) co-speciation, (b) independent speciation, (c) extinction, (d) “missing the boat”, (e) failure to speciate, and (f) host switching. The solid green region shows the lineage of the host, and the black line within shows the lineage of the parasite. Adapted from Page (Reference Page and RDM2003).
In the natural world there are many examples of where these co-evolutionary scenarios have shaped host–parasite interactions. The evidence of co-speciation is apparent in the classic study of pocket gophers (Rodentia: Geomyidae) and their chewing lice (Phthiraptera, Trichodectidae) parasites, where the gopher phylogenetic tree appears to be a mirror image of that of its host (Reed and Hafner, Reference Reed and Hafner1997). On the other hand, speciation following host-switching is thought to be the cause of diversification of the modern species of nematodes from the genus Trichinella, in the past 20 years (Cirtwill et al., Reference Cirtwill2016). Given that co-speciation and host-switching represent two extremes of a continuum (Page, Reference Page and RDM2003), it is not surprising to come across examples where speciation in parasite lineages is independent of the host as a result of duplication, independent assortment, and/or spreading. For instance, as many as eight such speciation events are reported to have occurred in the co-evolution of parasitic lice (Phthiraptera) and their petrel and penguin bird hosts (Paterson et al., Reference Paterson, Gray and Wallis1993). Lineage extinction as well as “missing the boat” can occur when a few founding members colonize a new region, leaving some of their original parasite species behind, or when a population undergoes a massive reduction in size or “bottleneck”. In this specific example, either of the two co-evolutionary scenarios is predicted to explain the absence of the two parasitic lice species Austromenopon himantopi and Saedmundssonia platygaster from the New Zealand black stilts Himantopus novaezealandiae despite these parasites being present in Himantopus himantopus, the closest relative of H. novaezealandiae (Paterson and Banks, Reference Paterson and Banks2001).
Marine skates of the order Rajiformes (referred to as skates hereafter) and their tapeworm parasites are a promising system in which to test the importance of both co-speciation and all other co-evolutionary scenarios in structuring host–parasite relationships. Skates are ecologically and economically important species, but they are now under threat as a result of overfishing coupled with other anthropogenic stressors, including habitat degradation and climate change (Ferretti et al., Reference Ferretti2013; Sguotti et al., Reference Sguotti2016; White et al., Reference White, Last, Séret and Last2016). As skates are definitive hosts for many tapeworm species, understanding these relationships can help us assess the risk of co-extinction faced by tapeworm parasites if vulnerable skate species such as Dipturus batis and Dipturus laevis become extinct (Strona, Reference Strona2015; White et al., Reference White, Last, Séret and Last2016). Furthermore, understanding the ecological parameters that facilitate co-evolutionary events in this system will provide some insight into the potential for persistence of a parasite species in the absence of its original host species, provided these ecological parameters remain undisturbed (Page, Reference Page and RDM2003; Dunn et al., Reference Dunn2009). It is crucial to evaluate the risks faced by marine parasites as they are indispensable for the health and stability of marine ecosystems. Parasites constitute a large biomass in marine ecosystems and a tremendous amount of energy flows through them via trophic interactions (Wood et al., Reference Wood2013).
The predominance of co-speciation in the skate–tapeworm system is ambiguous (Olson et al., Reference Olson1999; Caira and Jensen, Reference Caira and Jensen2001; Palm et al., Reference Palm2009; Bernot et al., Reference Bernot, Caira and Pickering2016). Co-speciation is considered important because several species of elasmobranch tapeworms are host-specific and are unable to exploit a wide range of parasite species – a characteristic linked to co-speciation (Caira and Jensen, Reference Caira and Jensen2014). On the other hand, tapeworm parasites are known to closely follow the ecology of their hosts. For instance, the life cycles of tapeworms of the orders Tetraphyllidea (sensu lato), Diphyllidea and Tetrabothriidea are concentrated in the shallower neritic zones where most of their potential definitive and intermediate hosts are present (Caira and Reyda, Reference Caira, Reyda and K2005); a constraint likely to erode the congruence between skate and tapeworm phylogenies. Accordingly, the aim of this study was to investigate the role of both phylogenetic relatedness and ecological similarity in structuring the tapeworm assemblages of marine skates.
Materials and methods
Skates sampled from the Falkland Island Shelf, Bay of Fundy, the North Sea and the south-western Pacific Ocean were included in the study (fig. 2a). A detailed map of the sampling locations in the Falkland Islands, where the majority of the host species were sampled, is provided in fig. 2b. In addition, samples from the Bay of Fundy were collected from the West Isles (approximately 47.9500°N, 67.0167°W), samples from the North Sea at 59.2486°N, 1.5033°E, and those in the south-western Pacific were retrieved from Nuggets, off the shore of the Catlins, New Zealand (approximately 46.4481°S, 169.8147°E). In total, 20 species of skates were sampled, and 59 species of tapeworm parasites were recovered from these specimens. DNA fragments consisting of 1044 bp of the NADH2 gene marker and 650 bp of the COI gene marker were sequenced from the skate tissue samples (table 1). A fragment of 1800 bp of the parasite 28S (large subunit ribosomal DNA) gene marker, encompassing the D1–D3 domains was sequenced for the tapeworm species (GenBank accession numbers MH686154, MH686155 and MH688697–MH688752). Additionally, sequences of some of the sampled host species already available in GenBank (Clark et al., Reference Clark2016) were included (table 1). Voucher specimens of parasites infecting Zearaja nasuta are deposited at Otago Museum. Accession numbers IV5692 to IV5695 were used for parasite taxa (all were hologenophores).
(a) Locations of skate collections for this study: the Falkland Islands Shelf, near the Falkland Islands; the North Sea, from the United Kingdom to Norway; the Bay of Fundy, near Canada; and the Pacific Ocean, near New Zealand. The sampling location in New Zealand was near the Catlins, in the South Island. (b) Sampling locations on the Falkland Islands Shelf.
Source of sequences of the skate species used in this study. N = number of specimens sampled in the study; (H) and (P) in the voucher number column indicate whether the voucher deposited with Otago Museum is a hologenophore or a paragenophore. Accession numbers in bold correspond to those sequences generated as part of this study.
DNA extraction and sequencing
DNA was extracted from skate and cestode tissues using the Chelex extraction protocol described in Casquet et al. (Reference Casquet, Thebaud and Gillespie2012), with an additional centrifugation stage at 36,000 rpm for ten minutes prior to an overnight incubation. For all three genes, polymerase chain reactions (PCRs) were performed using a total volume of 20 μl and the 5x MyTaq Red Reaction Buffer Kit (Bioline [Australia] Pty. Ltd, 2016). PCRs were performed at different concentrations for different gene markers. The PCR protocol for 0.5 μl of DNA template included 5 μl of 5x buffer, 0.1 μl of Taq, and 0.35 μl each of forward and reverse primers at 10 μm concentration. The primers used in DNA amplification and sequencing are described in table 2. Between 0.5 and 3 μl of NADH2 gene template was amplified in a Mastercycler pro S. The amplification protocol was as follows: initial denaturation for 4 minutes at 94°C followed by 38 cycles of denaturation for 30 s at 94°C, annealing for 30 s at 50°C, extension for 2 minutes at 72°C, and the final extension for 7 minutes at 72°C (Naylor et al., Reference Naylor2012). A volume of 3 μl of COI gene template was amplified as follows: initial denaturation for 2 minutes at 94°C followed by 35 cycles of denaturation for 30 s at 94°C, annealing for 4 minutes at 52°C, extension for 1 minute at 72°C, and the final extension for 10 minutes at 72°C (Spies et al., Reference Spies2011). A volume of 0.5 μl of 28S gene template was amplified using the following protocol (Randhawa et al., Reference Randhawa2008): initial denaturation for 4 minutes at 94°C, 34 cycles of denaturation for 30 seconds at 94°C, annealing for 30 s at 50°C, extension for 90 s at 72°C, and the final extension for 7 minutes at 72°C. The PureLink Quick Gel Extraction Kit was used to purify the DNA amplified from all the three gene markers by following the manufacturer's protocol (Life Technologies, 2011). Purified DNA was sequenced at Genetic Analysis Services at the Department of Anatomy, University of Otago, Dunedin, New Zealand.
List of PCR and sequencing primers and their sources for all the three gene markers used in this study.
Phylogenetic analysis
Sequences were aligned and edited in Geneious (v. 8.0.5) (Kearse et al., Reference Kearse2012). For the construction of the skate phylogenetic tree, the sequences generated from the NADH2 gene marker, the COI gene marker, and those obtained from GenBank for the two gene markers were multiple aligned using the MUSCLE algorithm. All parasite gene sequences were multiple aligned using the MUSCLE algorithm. For the skate phylogenetic tree, Hydrolagus colliei (spotted rat fish: Chimaeriformes) and Heptranchias perlo (sharpnose sevengill shark: Hexanchiformes) were used as outgroups (Aschliman et al., Reference Aschliman2012). For the tapeworm phylogenetic tree, Gyrocotyle rugosa (a tapeworm of chimaeras from the cestodarian order Gyrocotylidea) and Monostephanostomum nolani (a digenetic trematode) were used as outgroups. Given the diversity of tapeworms of elasmobranchs, including skates, members of the phylum Platyhelminthes other than true tapeworms or Eucestoda, such as a cestodarian and a digenean parasite, are considered by parasitologists to be suitable outgroups for phylogenetic tree construction (Caira et al., Reference Caira2014).
The aligned host and parasite sequences were analysed in PartitionFinder (v. 2.1.1), using the AICc selection criteria to identify appropriate substitution models for tree construction (Lanfear et al., Reference Lanfear2012). Following the outcome of PartitionFinder, the trees were constructed using GTR substitution model with gamma-distributed rate variation (GTR + G). The data sets were partitioned by gene and codon positions 1, 2 and 3. Phylogenetic trees were constructed using Bayesian inference in MrBayes (Ronquist et al., Reference Ronquist2012). For both skate and tapeworm phylogenetic trees, MrBayes was run for 50,000,000 generations, sampling every 500 generations following a burn-in of 10,000,000 generations. 50% majority rule consensus trees were then calculated from the posterior distribution of trees. These trees were compared to trees constructed using the Maximum Likelihood approach in RAxML (Stamatakis, Reference Stamatakis2006). As all the species used in the study are extant, the resultant phylogenetic trees were made ultrametric using the chronopl function in the R package “ape” with the lambda parameter set to 1 (Paradis et al., Reference Paradis, Claude and Strimmer2004; R Core Team, 2016).
Ecological data
Three ecological variables – host diet, depth range and average size – for each skate species were obtained from the literature to allow tests for ecologically mediated co-evolutionary events. Geographical location was recorded as the region in which each skate species was sampled. Diet composition, depth and size were extracted from published literature (Bizikov et al., Reference Bizikov2004; Ebert and Bizzarro, Reference Ebert and Bizzarro2007; Belleggia et al., Reference Belleggia2014) and FishBase (Froese and Pauly, Reference Froese and Pauly2016). These variables were then converted to ecological distance matrices between skate species, both individually and combined across all four variables. Pair-wise ecological distances were calculated based on the approach described in Geange et al. (Reference Geange2011):
(1) Diet: The proportion by weight of each of the important prey items such as fish, polychaetes, molluscs, amphipods and echinoderms were obtained for each skate species. Dietary niche overlap was calculated as the sum of minima for each dietary value, subtracted from 1 to obtain a distance measure:
where i and j represent the indices of two hosts, M is the number of dietary components, and 0 ≤ pi,m ≤ 1 is the proportion of component m in the diet of species i.
$$D_{i,j}^{diet} = 1 - \sum\nolimits_{m = 0}^M {\min \{ p_{i,m},p_{\,j,m}\}}, $$(2) Depth range: The intersection between the maximum and the minimum foraging depth of a species was calculated, and then divided by their union. A distance measure was obtained by subtracting from 1:
where each pair a i ≤ b i specifies the depth range of species i.$$D_{i,j}^{depth} = 1{\rm -} \max \left\{ {\displaystyle{{\min \lcub {b_i,b_j} \rcub {\rm -} \max \lcub {a_i,a_j} \rcub } \over {\max \lcub {b_i,b_j} \rcub {\rm -} \min \lcub {a_i,a_j} \rcub }},0} \right\},$$(3) Size: The average total length of males of each species was log-transformed to accommodate for the large variation in sizes. The pairwise differences between species were normalized by dividing by the largest difference to obtain values between 0 and 1:
$$C_{i,j} = \left\vert {\log s_i{\rm - }\log s_j} \right\vert$$$$D_{i,j}^{size} = \displaystyle{{C_{i,j}} \over {\max C}},$$where s i is the average length of species i.
(4) Geographical location: Geographical distance between species was calculated as a binary value: 1 if the species were sampled from different regions and 0 if they were found in the same regions:
where l i is a discrete variable taking on different values for different geographical clusters. Use of great circle distance was not suitable as the sampling regions were too scattered to form a representative cluster. Moreover, great circle distance can be misleading, as swimming distance may not reflect the physical distance between two locations.$$D_{i,j}^{geography} = \left\{ {\matrix{ 0 & {if\,l_i = l_j} \cr 1 & {otherwise} \cr}} \right.,$$
The four ecological variables were also combined into a single distance matrix using a Euclidean metric. This was done by calculating the square root of the sum of the squares of the four component distance values for a pair of species:
$$D_{i,j} = \sqrt {D_{i,j}^{diet^2} + D_{i,j}^{depth^2} + D_{i,j}^{size^2} + D_{i,j}^{geography^2}} $$Co-speciation and host-switching
A host–parasite association matrix was created based on the specimens found in this study, and a tanglegram produced with phytools (Revell, Reference Revell2012) was used to visually assess the degree of congruence between host and parasite phylogenetic trees. A similar tanglegram was produced to compare the parasite phylogenetic tree to a dendrogram produced using a hierarchical clustering analysis of hosts based on their ecological distance.
The program Jane 4 was then used to estimate the frequency of co-speciation and host-switching events (Conow et al., Reference Conow2010). As Jane uses an event-based approach, assigning costs to various events, a range of different cost values was tried. The cost of co-speciation was never higher than the cost of duplication, host-switch or loss, and was always set to 0. The costs of failure to diverge and loss varied between 0 and 2. The cost of duplication varied between 0 and 1, and the cost of host-switching and duplication also varied between 0 and 1, with the constraint that the cost of host-switching and duplication was never less than the cost of duplication (Míguez-Lozano et al., Reference Míguez-Lozano, Rodriguez-Gonzalez and Balbuena2017). The cost scheme with all costs set to 0 was also excluded.
The PACo method was employed to statistically test for significant congruence between the host and parasite phylogenies. This method carries out a principal coordinate analysis on host and parasite phylogenetic distance matrices separately, uses Procrustes superimposition to maximize alignment, then permutes the host–parasite association matrix to test for significant congruence (Balbuena et al., Reference Balbuena, Míguez-Lozano and Blasco-Costa2013, Reference Balbuena2017).
A slightly modified PACo analysis was used to test the importance of ecological factors in skate–tapeworm associations. This was done by using the host ecological distance matrix described above in place of the host phylogenetic distance matrix. If all co-evolutionary events except co-speciation are dependent on the ecological similarity between hosts, it is expected that closely related parasite species will occur in ecologically similar hosts, not necessarily in closely related hosts. The ecological distance matrix, combining diet, distribution depth, size and geographical location, was used along with the parasite phylogenetic distance matrix and the host–parasite association matrix to test for an association with the parasite phylogeny. Additionally, host diet and depth were used separately to investigate the individual importance of these ecological parameters. Geographical location could not be tested separately due to the use of binary data with very few data points and coarse measurement. PACo could not be run using host size alone due to insufficient dimensionality.
A Mantel test was used to test for phylogenetic signal in the ecological variables by comparing the host phylogenetic and the ecological distance matrices combining the host diet, distribution depth, average size and geographical location. A Mantel test was a more suitable choice than other measures of phylogenetic signal such as Pagel's lambda and Blomberg's K, as it can be used for composite distance matrices or for variables such as diet composition that cannot be expressed as a single continuous trait (Mantel, Reference Mantel1967; Pagel, Reference Pagel1999; Blomberg et al., Reference Blomberg, Garland and Ives2003).
Results
Both RAxML and MrBayes (supplementary figs S1 and S2) generated trees with the same structure for skates as well as for the tapeworms, with all known taxa recovered as monophyletic. In contrast with previous suggestions by Healy et al. (Reference Healy2009) and Caira et al. (Reference Caira2014) that Zyxibothrium might share affinities with the Rhinebothriidea, the lone Tetraphyllidea (sensu stricto) (Caira et al., Reference Caira, Jensen, Ruhnke, Caira and Jensen2017) in our analyses, Zyxibothrium kamienae, resolved as the sister group to Diphyllidea + Trypanorhyncha + Onchoproteocephalidea + Phyllobothriidea with ≥ 95% posterior probability and not as a member of the Rhinebothriidea. This species diverged significantly from the other Phyllobothriidea genera, Guidus, Phyllobothrium, and Phyllobothriidea New Genus. The phyllobothriidean P. piriei from the north-east Atlantic clustered within the New Genus recovered from skates in the south-west Atlantic and resolved within the Phyllobothriidea with a posterior probability value of ≥ 95%.
Visual examination of the phylogenetic trees using a tanglegram (fig. 3a, b), as well as the history reconstructed by Jane (fig. 4a, b, c), indicated a substantial number of co-speciation events. However, other co-evolutionary scenarios also exist in this host–parasite system. Different co-evolutionary scenarios based on 26 different cost combinations were re-constructed by Jane. The number of occurrences of each of the five co-evolutionary events within three groups of the significant cost combinations is shown in fig. 4a, b, c. The number of losses was always high. There was always a reasonable number of co-speciation events, except in scenarios where the cost for duplication and duplication with host-switching were both set to zero. In these scenarios there was no incentive for the Jane algorithm to find co-speciation where it was present, as host-switching and duplication were always able to explain the outcome in these cost schemes without any penalties. In the remaining scenarios there was always a considerable number of host-switching events.
(a) Host–parasite associations between the skate hosts and tapeworm parasites used in this study. Both host and parasite trees were generated using MrBayes using a GTR + G model. Analysis was run for 50,000,000 generations, sampling every 500 generations following a burn-in of 10,000,000 generations. (b) Host–parasite associations based on the host ecological characteristics, where the host dendrogram is built by clustering pairwise ecological distance measures using hierarchal clustering with the hclust function in phytools. The parasite phylogenetic tree was built in MrBayes using a GTR + G model. Analysis was run for 50,000,000 generations, sampling every 500 generations following a burn-in of 10,000,000 generations. See Appendix for key to parasite and host species abbreviations.
Box-plots showing the number of occurrences of different co-evolutionary scenarios, namely co-speciation (C), host-switching (HS), duplication (D), failure to diverge (FD), and losses (L), produced by Jane. The cost of co-speciation was always set to 1. Cost for other evolutionary scenarios varied in each block represented in the image. (a) Results from the following cost combinations: D = 0, HS = 0, L = 1–2 (values ranging between 1 and 2), FD = 0–2; (b) results from D = 0–1, HS = 1, L = 1–2, FD = 0–2; (c) results from D = 1, HS = 1, L = 0, FD = 0–2.
Similarly, the test for co-speciation using PACo yielded a residual sum of squares value (m 2) of 20.44 with a permutational value of p = 0.0002 for the overall global fit of tapeworms with their skate hosts, indicating more congruence between the phylogenies than expected by chance alone. The Procrustean plot (fig. 5a) indicated that the hosts and the tapeworm parasites were roughly segregated into three distinct clusters. The first cluster was formed by all Psammobatis spp. and their parasites, the second cluster was formed by Leucoraja, Malacoraja, Zearaja, Amblyraja and Raja species with their parasites, and all Bathyraja species formed the third cluster. The Procrustean plot (fig. 5a) also indicated that although co-speciation was important, it was not the sole influence on host–parasite assemblages. The jack-knife plot (fig. 5b), obtained using the phylogenetic distance in PACo indicated that some members of Phyllobothriidea New Gen. sp. colonizing Bathyraja multispinus, B. albomaculata, B. macloviana, B. scaphiops and B. cousseauae, as well as Rhinebothriidea New Gen. spp. infecting the Psammobatis spp., were among the parasites that closely fit the co-speciation hypothesis.
(a) Procrustean superimposition plot produced by phylogenetic distance between hosts and parasite species in PACo. Blue squares represent host species and red dots indicate parasite species. The black lines show the links between hosts and parasites. The length of lines represents the residuals on the first and second axes; shorter lines and more spread out host–parasite clusters mean more co-speciation. (b) Jack-knife plot using the phylogenetic distance between host and parasite species, generated using PACo analysis. The plot represents the contribution of each host–parasite specimen to the sum of squares error. Those under the red line fit the co-speciation hypothesis more closely. See Appendix for key to parasite and host species abbreviations.
The parasite phylogeny fit the host ecological distance matrix combining the four variables, with an m 2 value of 33.81 and a p value of 0.016. The Procrustean plot (fig. 6a) did not illustrate any clear clustering of species. The jack-knife plot (fig. 6b) showed that Z. kamienae and Pseudanthobothrium hanseni from Malacoraja senta, Grillotia sp. from B. brachyurops, P. piriei from Leucoraja naevus, P. hanseni and Echeneibothrium canadensis from Amblyraja radiata, and Echeneibothrium sp. from Bathyraja griseocauda and B. scaphiops provided a close fit to the ecological data. Rhinebothriidea New Gen. spp. from Psammobatis spp. also fit the host-switching hypothesis. On fitting the parasite phylogenetic tree with the ecological variables separately (excluding geographical location, and average length for the reason described above), only diet was found to be significant (p = 0.032), whereas depth did not yield a significant result (p = 0.451). The Mantel test indicated a moderate positive correlation between phylogenetic and ecological distances (r = 0.17, p = 0.023).
(a) Procrustean superimposition plot generated using ecological distance measures, combining all four ecological variables, between host and parasite species in PACo. Blue squares represent host species and red dots indicate parasites. The lines between hosts and parasites indicate host–parasite associations. Shorter lines indicate that the host–parasite relationship is more dependent on ecological parameters. (b) Jack-knife plot using ecological distance measure between hosts and parasites in PACo, showing the contribution of each data point to the result. The host–parasite specimens under the red line fit the ecological data more closely. See Appendix for key to parasite and host species abbreviations.
Discussion
To the best of our knowledge, this was the first study to examine incongruence in a skate–tapeworm system, and to pinpoint some of the key ecological parameters leading to this incongruence. As hypothesized, co-speciation was found to be important in the evolution of skates and their tapeworm parasites. The importance of co-speciation in tapeworm–elasmobranch relationships was highlighted by Olson et al. (Reference Olson1999), who studied the evolutionary history of tetraphyllidean and lecanicephalidean tapeworms. Recently, Bernot et al. (Reference Bernot, Caira and Pickering2016) reported congruence between tapeworms of the genus Symcallio and their shark hosts of the genus Mustelus. However, while co-speciation clearly occurs in the skate–tapeworm system, it does not explain all host–parasite associations. In another study, Olson et al. (Reference Olson2010) found that several species of tapeworms of the order Trypanorhyncha did not appear to co-speciate with their elasmobranch hosts. They reported that a clade containing species of Grillotia and Pterobothrium contained species parasitizing rays and species parasitizing sharks, indicating shifts between distantly related host species. Accordingly, the present study confirmed that co-speciation in the skate–tapeworm system is not ubiquitous, and that other co-evolutionary events are relatively common. There is some evidence in the literature that strict host-specificity in the elasmobranch–tapeworm host–parasite system may not be indicative of strict co-speciation (Caira and Jensen, Reference Caira and Jensen2001). It is also possible that the lack of knowledge regarding the diversity of some tapeworm taxa could be exaggerating the evidence of their host specificity (Bernot et al., Reference Bernot, Caira and Pickering2016).
Our results suggest that the incongruence between skate and tapeworm phylogenies is influenced by a combination of the four ecological factors tested in this study. The influence of host diet in determining infection by tapeworms is evident in the eggs of some tapeworm species that imitate the important prey items of their hosts, increasing the prospect of their transmission (Caira and Reyda, Reference Caira, Reyda and K2005). Given that tapeworms are transmitted trophically, host diet intuitively plays a key role in determining their chances of transmission (Cirtwill et al., Reference Cirtwill2016). Likewise, the influence of the depth distribution of host species is evident in life-history characteristics of several marine tapeworm species of the order Tetraphyllidea (sensu lato) (Caira and Reyda, Reference Caira, Reyda and K2005). These parasites have egg pockets that bear either filaments for attachment to vegetation at the bottom of the sea or floating devices for floating to the top, depending on whether their host species is benthic or pelagic (Caira and Reyda, Reference Caira, Reyda and K2005). Similarly, an influence due to geographical characteristics of a host species was reported in a study by Randhawa and Poulin (Reference Randhawa and Poulin2010), in which tapeworm assemblages in elasmobranch hosts were influenced by the host's latitudinal range. Finally, host body size can determine the parasitic assemblages of the host species, as in the case of the European bitterling Rhodeus amarus (Dávidová et al., Reference Dávidová2008). Although the four ecological parameters were significant in combination, the relatively small number of host species in this study may have precluded the detection of significant effects for most individual traits. First-hand measurement of ecological variables at the time and site of sampling of host species could also improve the ecological data by incorporating variation among individuals and changes in diet through ontogeny.
Exploring the role of ecological parameters in shaping skate–tapeworm relationships has provided valuable insight into the vulnerability of these parasites to co-extinction with their host species. Skates are fished worldwide for food, which has subjected them to immense fishing pressure, driving several species, including Dipturus batis, to the brink of extinction (White et al., Reference White, Last, Séret and Last2016). This study predicts that some tapeworms, such as Echeneibothrium vernetae, that co-speciate with their Leucoraja hosts are likely to face secondary extinction in the event of their hosts going extinct. On the other hand, those more vulnerable to perturbations in the biology or ecology of their respective host species include Z. kamienae and P. hanseni from M. senta, Grillotia sp. from B. griseocauda, E. canadensis from A. radiata, and Echeneibothrium sp. from B. griseocauda and Bathyraja scaphiops. Host diet, distribution depth, size and geographical location together are important in influencing patterns of co-evolution. Environmental perturbations, such as loss of important prey species, could be detrimental to the fate of those tapeworms that are more likely to be influenced by other co-evolutionary scenarios, such as host-switching (Brooks and Hoberg, Reference Brooks and Hoberg2007). Parasites constitute a copious biomass in marine ecosystems, facilitating a tremendous flow of energy (Wood et al., Reference Wood2013). A marine ecosystem bestowed with a large diversity, displaying complex trophic interactions, is healthier and more resilient (Hudson et al., Reference Hudson, Dobson and Lafferty2006; Worm et al., Reference Worm2006). Therefore, parasite diversity is indispensable to the health and stability of marine ecosystems. Highlighting the vulnerability of skate-infecting tapeworms to co-extinction and environmental perturbations, this study has implications for marine conservation the world over (Wood et al., Reference Wood2013).
Although several gaps regarding the role of co-evolution have been filled by this study, further questions have also been raised. This study represents the first application of PACo to investigate the role of ecological as well as phylogenetic distances between hosts. However, further extensions of this or other methods will be required to measure the relative importance of co-speciation and other co-evolutionary events, and to robustly account for the potentially confounding effect of phylogenetic signal in ecological variables (this study found a moderate positive correlation between phylogenetic and ecological distances). The development of methods that estimate the importance of phylogenetic and ecological factors simultaneously would be a valuable advance that would allow tests for a more diverse array of host–parasite co-evolutionary scenarios (Clark and Clegg, Reference Clark and Clegg2017). This study could have been improved by access to a larger number of species occurring across a greater breadth of geographical locations. However, PACo is reasonably robust to incomplete phylogenetic trees (relative to Jane and other techniques), so the missing data should not have resulted in bias (Balbuena et al., Reference Balbuena, Míguez-Lozano and Blasco-Costa2013).
In conclusion, both phylogenetic and ecological parameters are important in shaping the evolutionary history of skates and their tapeworm parasites. The ecology of skate hosts has emerged to be important in shaping skate–tapeworm relationships, with the host diet, depth, size and location being the key players in facilitating host-switching. Given the importance of parasites in ensuring a healthy, resilient and stable marine ecosystem, it is essential to understand these ecological characteristics. Conservation of skate species that harbour co-speciating parasites is needed to prevent host–parasite co-extinction. Marine ecosystems globally are threatened by multiple stressors, including climate change, habitat degradation and fishing pressures. This study may provide valuable information for scientists and managers that could facilitate conservation of entire global marine ecosystems, including conservation of parasite species.
Supplementary material
To view supplementary material for this article, please visit https://doi.org/10.1017/S0022149X18000809
Acknowledgements
We would like to acknowledge Professor Robert Poulin, Dr Bronwen Presswell, Dr Tania King and Dr Fátima Jorge for advice that contributed to the successful completion of this project. Samples were collected on board the following vessels: W. B. Scott R/V and the CCGS Pandalus III R/V (Canada), Scotia R/V (North Sea), Echo F/V (New Zealand), and Dorada R/V (Falkland Islands). Assistance with collection of samples by crew and officers from these vessels is gratefully acknowledged, as are the contributions of scientists and fisheries observers of the Falkland Islands Fisheries Department (Falkland Islands) and Fisheries Research Service (Scotland, UK) who assisted HSR with sample collection.
Financial support
We acknowledge the contributions of Otago Museum in granting the Otago Museum Zoological Research Scholarship to AB to partially fund this research project, and housing the molecular vouchers in their collection. Additional financial support was provided by the Departments of Botany and Zoology of the University of Otago to AB in the form of a student budget for this project. The Falkland Islands Fisheries Department funded sequencing performed by HSR. Collection trips by HSR to the North Sea and Falkland Islands were funded by a John S. Little International Study Fellowship (University of New Brunswick) and the Falkland Islands Fisheries Department, respectively.
Conflict of interest
None.
Appendix
