Shared phylogeographic patterns between the ectocommensal flatworm Temnosewellia albata and its host, the endangered freshwater crayfish Euastacus robertsi

Study area

The study area was located in the Daintree rainforest, which is the largest continuous rainforest on the continent. Situated within the wet tropics in northern Queensland, Australia, the tropical climate has hot wet summers and cool dry winters. The three mountaintops inhabited by T. albata and E. robersti are Mount Pieter Botte (elevation 1,009 m), Thornton Peak (1,375 m) and Mount Finnigan (1,083 m) (Fig. 1, Table S1). The area which is <100 km², is believed to contain the entire range of both species (Morgan, 1997; Ponniah & Hughes, 2006; Coughran & Furse, 2010). Samples were collected, as per Ponniah & Hughes (2004), in stream reaches at elevation >750 m. On each mountaintop one site was sampled except for Mount Finnigan where two stream reaches were sampled.

DNA extraction, amplification, and sequencing

Total genomic DNA was extracted from the leg tissue of E. robertsi and from whole samples of T. albata as per methods outlined in Carini & Hughes (2006). A final edited 610 base pair fragment (E. robertsi) and 603 base pair fragment (T. albata) of the mtDNA COI gene was produced after polymerase chain reaction (PCR) using the COI primer set LCO-1490 and HCO-2198 of Folmer et al. (1994). PCR conditions were: denaturation of DNA occurred at 95 °C for 5 min, followed by 30 cycles of 94 °C denaturing for 1 min, 55 °C annealing for 30 s, and 72 °C extension for 1 min, followed by a final 68 °C extension step for 5 min. Dye terminator cycle sequencing reactions were used for sequencing (Perkin Elmer, Foster City, CA) as per manufacturer’s instructions. Sequencing was carried out on an Applied Biosystems (Foster City, CA) 3130xl automated sequencing machine.

We also wanted to include nuclear gene data in the analysis to compare phylogenetic patterns for the two species. Therefore, we included a 734 base pair edited fragment (E. robertsi) and a 692 base pair edited fragment (T. albata) of 28S ribosomal DNA. We used the primers Rd1a and Rd4b (Crandall, Harris & Fetzner, 2000). PCR conditions followed those of Crandall, Harris & Fetzner (2000).

The nucleotide sequences for COI and 28S were aligned and edited with SEQUENCHER v4.9 (Gene Codes Corporation). The mtDNA sequences were visually assessed for the occurrence of nuclear mitochondrial pseudogenes (numts) using techniques described in Bensasson et al. (2001) and none were found.

Networks, phylogenetic trees and divergence estimates

For the 28S data we constructed haplotype networks using TCS v1.21 (Clement, Posada & Crandall, 2000) for E. robertsi and T. albata. Phylogenetic trees were constructed for the COI data using unique haplotypes. For the E. robertsi tree, E. fleckeri was selected as an out-group, as it has previously been shown to be the sister species of E. robertsi (Ponniah & Hughes, 2004). Temnosewellia aphyodes was chosen as an out-group for the T. albata tree, as it is the resident flatworm of E. fleckeri (Sewell, Cannon & Blair, 2006). We used jModeltest v0.1 (Posada, 2008) to choose the best-fit substitution model for each COI dataset. Using the Akaike information criterion the model selected for T. albata was TPM2uf + I + G and for E. robertsi was TIM3 + I. Tree construction for each data set was run using Bayesian analyses. MrBayes v3.1.2 (Huelsenbeck & Ronquist, 2001) was used for tree topology comparison, and BEAST v1.7.5 (Drummond & Rambaut, 2007; Drummond et al., 2012) was used to construct rooted ultrametric trees for comparison of node divergence times.

In MrBayes, a MCMC chain of 2,000,000 iterations was used with a sample frequency of 100. The first 25% of iterations were discarded as burn-in. MEGA v5.10 (Tamura et al., 2011) was used to calculate uncorrected percentage divergence between clades. In BEAST a lognormal relaxed clock model was first used to estimate divergence times of clades. However, the data could not reject a strict clock (ucld.stdev included zero); therefore, a strict clock model was used along with a coalescent constant size tree prior. Owing to lack of fossil calibration points and uncertainties in transferring molecular clock rates across taxa, we chose to incorporate a range of rates from the literature to place an approximate time-frame on COI divergences within the E. robertsi and T. albata datasets. Clock rates were used to describe a lognormal prior for the estimated clock rate, where 95% of the probability density was contained within highest and lowest values taken from the literature. For the temnocephalans these values were 0.0027 and 0.015 substitutions/site/lineage/million years (Schmidtea mediterranea, Platyhelminthes: Lazaro et al., 2011; Dugesia, Platyhelminthes: Sola et al., 2013). For the crayfish the values were 0.0083 and 0.012 substitutions/site/lineage/million years (Chirocephalus, Crustacea: Ketmaier, Argano & Caccone, 2003). Convergence, mixing and effective sample size of model parameters (>200) was assessed using the program Tracer v1.5 (Drummond & Rambaut, 2007) after running the analysis for 10⁸ generations.

To investigate the magnitude of genetic divergence within each taxon, without reliance on transformation using molecular clock rates, we fitted the COI and 28S datasets to a two population isolation-with-migration model (IM), implemented in the software IM (v.12/17/2009; Hey & Nielsen, 2004). Three pair-wise population comparisons (among the three mountaintop populations) were made for E. robertsi and T. albata. To ensure that results were consistent, each pair-wise comparison was run a minimum of three times (18 h/run) with different random number seeds. Model parameters of interest were taken from the peaks of the estimated distributions. These were population splitting time (t) scaled by the (unknown) geometric mean of the mutation rates for COI and 28S, and between-population migration rates (m1, m2).

Temnosewellia albata

A total of 63 T. albata individuals were taken from 20 crayfish hosts sampled across five locations on three mountaintops (>700 m above sea level) (see Supplemental Information). Sixty one T. albata were sequenced for 603 bp of the COI mtDNA region. Seventeen unique haplotypes were identified (GenBank Accession: Table 1, see Supplemental Information). The target fragment contained no gaps and was variable at 91 sites (15%) of which 87 were parsimony informative (14%, Table 2). Between one and ten T. albata were sequenced for COI per crayfish (mean = 3.2). Average heterozygosity of T. albata sampled from one individual crayfish was not consistently lower compared to the average heterozygosity of T. albata sampled from a number of different crayfish. Due to sequencing issues, a small subset of the T. albata was used in the sequencing of the 28S ribosomal DNA region. For the 28S region, eight samples were sequenced for 692 base pairs and two unique haplotypes were identified (Table 1). The target fragment contained one gap at site 20 and had nine (1.3%) variable sites, all of which were parsimony informative (Table 2).

Euastacus robertsi

For E. robertsi, 610 base pairs of the COI region were available for 64 individuals (including 16 from Genbank) (see Supplemental Information). Six unique haplotypes were discovered (GenBank accession: Table 3, see Supplemental Information). The target fragment contained no gaps and was variable at 31 sites (5%) of which 30 sites (97%) were parsimony informative (Table 2). The out-group consisted of three sequences from E. fleckeri (including two from GenBank). Twenty five of the E. robertsi individuals were sequenced for 733 base pairs of 28S ribosomal DNA and, four unique haplotypes were identified (Table 3). The target fragment contained no gaps and had 12 variable sites (1.6%; Table 2) 11 sites were parsimony informative (92%).

Genetic variation

Comparison of COI tree topology between species

Tree topologies for both taxa featured two or three deeply divided in-group clades (Fig. 2A). The main difference between taxa was, for T. albata there was division into three well-supported clades (clades A, B, C), whereas for their crayfish hosts E. robertsi, two well-supported clades (D, E) were identified. The three T. albata clades corresponded strongly with the three separate mountaintop populations, although one clade (B) included a mixture of samples from Thornton Peak and Mt Finnigan. The spatial pattern of the two clades in the E. robertsi tree grouped together samples from Mt Pieter Botte and Thornton Peak together (clade D). Clade E was comprised mostly of samples from Mt Finnigan.

Divergence estimates and isolation-with-migration model

Percent divergence of the T. albata COI clades (shown in Fig. 2) was 12% for both clades A–C and B–C, and 10% for clades A–B. The median divergence time calculated by BEAST for clades A–B was ∼11 mya (Table 4). Percent divergence for E. robertsi for clades D–E was 5%. The median divergence time calculated by BEAST for clades D–E was ∼2.6 mya. Evaluation of the population divergence time parameter (t), incorporating both COI and 28S data in an IM model revealed no difference between T. albata and E. robertsi for the three among-mountaintop comparisons (Fig. 3; Table 5). Mean point estimates of t were 0.303 for crayfish and 0.338 for temnocephalans with broadly overlapping 95% credibility intervals. Note that these values are unscaled parameter estimates. Conversion into units of real time would require scaling these values by the (unknown) substitution rates for each locus and by the (unknown) generation time of each species.

Table 3:

GenBank accession numbers for Euastacus robertsi. All the new sequences that were generated as part of this study are highlighted in boldface type.

Molecular marker	Haplotype ID & GenBank accession number	Mountain top location
COI	FI-A; DQ006368, DQ006372, DQ006377,	Mt Finnigan/Thornton Peak
	DQ006378, AY800362, DQ006369, AY324346
COI	FI-P; KJ939254	Mt Finnigan
COI	TP1; DQ006370, DQ006376	Thornton Peak
COI	PB1; DQ006373, AY800364, AY324347	Mt Pieter Botte
COI	TP3; AY800363	Thornton Peak
COI	TP2; KJ939253	Thornton Peak
28S	F1; EU920988	Mt Finnigan/Thornton Peak
28S	P1; KJ941016	Mt Pieter Botte
28S	T1; KJ941015	Thornton Peak
28S	T2; KJ941017	Thornton Peak

DOI: 10.7717/peerj.552/table-3

Table 4:

Divergence estimates between clades (see Fig. 2A).

Comparison is between sequenced COI haplotypes, for Temnosewellia albata and Euastacus robertsi in Queensland, Australia.

	Pairwise comparison	Substitutions/site/lineage/million years		Diverged (mya)	±95% HPD (mya)
		Lower	Upper
Temnosewellia albata	Clade A–B	0.0027	0.015	11	3.4–25
	Clade A–B–C	0.0027	0.015	15	4.4–32
Euastacus robertsi	Clade D–E	0.0083	0.012	2.6	1.5–4.3

DOI: 10.7717/peerj.552/table-4

Table 5:

Parameter estimates from an isolation with migration model for three pairwise mountaintop population comparison.

	Temnosewellia albata				Euastacus robertsi
Comparison	Parameters	HiPt	HPD90Lo	HPD90Hi	Parameters	HiPt	HPD90Lo	HPD90Hi
	qF	1.702	0.361	6.448	qF	0.1700	0.0154	1.126
1. Mt Finnigan	qP	0.147	0.049	1.522	qP	0.0163	0.0163	0.992
–	qA	44.82	19.75	103.1	qA	30.85	6.423	32.50
2. Mt Pieter Botte	t	0.308	0.053?	14.99?	t	0.173	0.022?	12.69?
	mF	0.005	0.005	0.515	mF	0.005	0.005	0.995
	mP	0.005	0.005	1.285	mP	0.005	0.005	1.355
	qF	0.965	0.170	4.596	qF	0.1	0.02	0.86
1. Mt Finnigan	qT	3.440	0.313	41.802	qT	1.661	0.2518	7.100
–	qA	38.46	13.03?	103.1	qA	27.545	10.52?	91.49?
2. Thornton Peak	t	0.878	0.068?	7.238?	t	0.6225	0.0375?	14.99?
	mF	0.005	0.005?	8.325?	mF	0.005	0.0050?	5.755?
	mT	1.225	0.005	7.27	mT	0.275	0.015	3.855
	qT	3.563	0.6600	16.75	qT	3.020	0.3583	16.63
1. Thornton Peak	qP	0.087	0.087	1.478	qP	0.0217	0.0217	0.5852
–	qA	49.21	13.03?	171.6?	qA	23.18	8.446	89.94
2. Mt Pieter Botte	t	0.292	0.068?	7.328?	t	0.5175	0.0225?	14.99?
	mT	0.055	0.005	1.685	mT	0.265	0.005	3.565
	mP	0.005	0.0050?	5.845	mP	0.005	0.0050?	8.795?

DOI: 10.7717/peerj.552/table-5

Notes:

HiPt

the value of the bin with the highest count

HPD90Lo

lower bound of the estimated 90% highest posterior density (HPD) interval. A question mark ‘?’ indicates unreliable or limit due to flat or incomplete posterior probability distribution sampled

HPD90Hi

upper bound of the estimated 90% highest posterior density (HPD) interval

q

the effective population size, population indicated by the letter (F, T, P)

qA

ancestral population

m

the migration rate per gene copy per generation, letters indicate the population (F, T, P)

t

a divergence estimate (not transformed to years)

Assessment of among-mountaintop migration using the IM model indicated no migration was compatible with the data for most pairwise comparisons (Table 5). However, the data did support non-zero migration from the Mt Finnigan population into the Thornton Peak population for both taxa: crayfish (m_FtoT = 0.2; i.e., parameter estimate converted to demographic units representing effective number of migrants per generation), and temnocephalans (m_FtoT = 1.9). Non-zero migration was also detected for crayfish from the Mt Pieter Botte population moving into the Thornton Peak population (Table 5). However, the temnocephalan data did not mirror this pattern.

Shared phylogeographic patterning

At two independent loci we identified haplotype sharing between mountains for the flatworm and its host. We cannot say with absolute certainty if these shared haplotypes are the result of gene flow or retention of ancestral haplotypes, but analysis using the IM model was compatible with low levels of unidirectional gene flow after population isolation. The pattern shared by crayfish and flatworms was for migration from Mt Finnigan into Thornton Peak. We established that, in both datasets, Mt Finnigan and Mt Pieter Botte were isolated from each other due to a lack of haplotype sharing with no evidence of migration. It has long been postulated that intervening lowlands are effective barriers to dispersal for Queensland Euastacus (Ponniah & Hughes, 2004; Ponniah & Hughes, 2006). It is also likely that the lack of migration for this species could be attributed to these two mountains being on completely separate ridges. Furthermore, we consider it possible that Mt Finnigan and Thornton Peak may have once shared a ridge making historical connections between them more likely. These connections may have been present either overland or through historical stream connections. Historical connections between these two mountains have been found for the beetle Philipis (Baehr, 1995). As Euastacus are known to be able to survive for long periods out of water they have the ability to traverse over land (Furse & Wild, 2002), although overland dispersal may be rarer in some species. Intervening high points along ridge lines may have allowed for historical migration pathways. Current elevations between these mountains are no lower than 350 m at some places. Therefore, although it has been shown to be a rare occurrence, it is possible that migration between sites is possible, at least between two of the mountain ridges.

The congruence that we observed in the phylogeographic pattern of T. albata and its crayfish host suggests that their evolutionary histories are spatially linked; therefore, if hosts are capable of overland dispersal, so are the flatworms. The exact mechanism of dispersal for Temnosewellia is unknown. As the genus is generally considered to be host-specific it is expected that, like other temnocephalans, they undergo their entire life cycle and subsequent generations on a single host crayfish (Sewell, Cannon & Blair, 2006). The mechanisms that allow them to survive the moult phase of their hosts are not known; however, observations by Haswell (1983) and Nichols (1975), on closely related temnocephalans, noted that they may be able to survive for some time in the absence of a host. Even though the flatworm’s ectoderm is somewhat prone to desiccation (Haswell, 1909; C Hurry, pers. obs., 2014), flatworms may still be able to disperse overland with their host due to the durability of their unhatched eggs. Temnosewellia will lay tens to hundreds of eggs which stick firmly to the exoskeleton of the crayfish (Wild & Furse, 2004; C Hurry, pers. obs., 2014). Eggs are enclosed in a tough outer coating and have a large fluid filled cavity (Haswell, 1909), which may prevent desiccation of the unhatched young allowing long distance movement in the absence of water. As so little is known of Temnosewellia this hypothesis has not been tested and further work is needed to determine dispersal mechanisms in these flatworms.

Time of divergence

As low levels of migration were detected it should be easier to detect founding events. Our results show that, for both these species, isolation and divergence among refugial mountaintop populations was old enough to have resulted in accumulation of mutational differences. Comparison of COI divergence times for a node marking the split between Mt Finnigan and Thornton Peak for both taxa suggested the temnocephalan divergence may be older than the crayfish (∼11 mya compared to ∼2.6 mya). However this result is contingent on calibration using molecular clock rates from the literature. Numerous studies have highlighted variability in substitution rates between taxa (Wilke, Schultheiß & Albrecht, 2009; Lanfear, Welch & Bromham, 2010), so caution is required in interpreting the divergence times presented here. By taking a different population-based approach—using the multilocus, isolation-with-migration model—we showed that the population splitting parameter t was indistinguishable between the two species. This comparison incorporates data from another locus in addition to COI, and does not depend on application of molecular clock rates (which may not be appropriate for our study species). However it does have the drawback of not being expressed in units of absolute time. Weighing-up both of these results leads us to conclude that either (1) the mountaintop divergences of both species did occur contemporaneously, but that a greater number of substitutions have become fixed in the mitochondrial genome of the flatworm compared to the crayfish (i.e., their divergence rates are different) or (2) T. albata did diverge earlier than their hosts. The second assumption is entirely possible if in the past T. albata had a different host which has now become extinct. A different host could either be an ancestral Euastacus which was not sampled as part of this study or another crustacean host altogether. If we were able to confirm the authenticity of the single T. albata sampled upon the crayfish Cherax depressus ∼515 km south of Mount Pieter Botte (Sewell, Cannon & Blair, 2006); we may find that in the past the distribution of this species was much wider. However, as this hypothesis is based upon just one sample we are cautious in offering this interpretation. Our calculations support a separation sometime during the Pliocene. As previously stated, conditions during the Pliocene are believed to have had a significant impact on population distribution and structure of fauna and flora in the north and the coastal east of Australia due to vicariance events.

Due to the close association that Temnosewellia share with their host they may be, in future studies, considered to be a suitable proxy in resolving phylogeographic patterning in their hosts. Nieberding & Olivieri (2007) tell us that ‘parasites’ that act as suitable proxy species are without intermediate hosts and have no phase of living independently of their host. Equally they are individuals which display smaller N_e at the population level and exhibit lower gene flow than their hosts among populations. These factors combined allow them to display a stronger population structure than their host, making them especially useful in cases where hosts are rare or hard to sample in large numbers. As Temnosewellia satisfies many of these criteria their role in future studies may be to help resolve host phylogenies or phylogeographic history. The existence in this study of three deeply divided COI lineages in the temnocephalan compared with two in the crayfish may indeed indicate that the ectocommensal genealogy records part of the crayfish history that is lost due to stochastic sorting of lineages and extinction/recolonisation events. One possibility is that a third crayfish mtDNA lineage did exist on Mt Pieter Botte, but was replaced by the Thornton Peak mtDNA lineage following a colonisation event. In this scenario the ectocommensal history acts as proxy for the crayfish history. However the extra temnocephalan mtDNA lineage might also be explained by lineage retention or by failure to detect a corresponding third lineage in our crayfish sample.

We were able to demonstrate that the association between E. robertsi and T. albata has likely persisted over several million years. The results from our study are applicable to host-commensal relationships worldwide, as they show that shared histories between such close commensal species may span millions of years. A growing number of examples in the literature are demonstrating that symbionts can be used to infer host history for conservation gains (Colwell, Dunn & Harris, 2012), which highlight the importance of studying symbiotic species alongside their hosts. We suggest that future phylogeographic studies exploit host-commensal interactions to provide objective measures of biodiversity, population subdivisions and phylogeographic information in the host. These interactions should be considered in management plans for crayfish species, especially as this technique may prove useful when host numbers are small, due to rarity or low catch rates.