Phylogeographic Structure in Penguin Ticks across an Ocean Basin Indicates Allopatric Divergence and Rare Trans-Oceanic Dispersal

The association of ticks (Acarina) and seabirds provides an intriguing system for assessing the influence of long-distance dispersal on the evolution of parasitic species. Recent research has focused on host-parasite evolutionary relationships and dispersal capacity of ticks parasitising flighted seabirds. Evolutionary research on the ticks of non-flighted seabirds is, in contrast, scarce. We conducted the first phylogeographic investigation of a hard tick species (Ixodes eudyptidis) that parasitises the Little Blue Penguin (Eudyptula minor). Using one nuclear (28S) and two mitochondrial (COI and 16S) markers, we assessed genetic diversity among several populations in Australia and a single population on the South Island of New Zealand. Our results reveal two deeply divergent lineages, possibly representing different species: one comprising all New Zealand samples and some from Australia, and the other representing all other samples from Australian sites. No significant population differentiation was observed among any Australian sites from within each major clade, even those separated by hundreds of kilometres of coastline. In contrast, the New Zealand population was significantly different to all samples from Australia. Our phylogenetic results suggest that the New Zealand and Australian populations are effectively isolated from each other; although rare long-distance dispersal events must occur, these are insufficient to maintain trans-Tasman gene flow. Despite the evidence for limited dispersal of penguin ticks between Australia and New Zealand, we found no evidence to suggest that ticks are unable to disperse shorter distances at sea with their hosts, with no pattern of population differentiation found among Australian sites. Our results suggest that terrestrial seabird parasites may be quite capable of short-distance movements, but only sporadic longer-distance (trans-oceanic) dispersal.


Introduction
Dispersal is a driving force in the isolation and subsequent evolution of species on both local and global scales [1]. However, long-distance dispersal events are often sporadic and difficult to predict, complicating the direct testing of hypotheses [1,2]. Molecular data can provide Only two major genetic studies have been conducted on hard ticks (Ixodes uriae) taken from penguin colonies [41,42]. These studies found that the I. uriae ticks in the northern and southern hemispheres are divergent [42], and that some intraspecific genetic isolation is discernable on much smaller geographic scales [41]. More generally, genetic studies of this tick species have supported a trend towards less genetic connectivity of hard ticks with greater geographic distance between penguin colonies. The driving forces behind evolution in Ixodid ticks have also been the subject of a number of genetic studies, again using I. uriae as a model [41][42][43][44][45][46]. These studies on I. uriae provide the closest analogue to the Little Blue Penguin-tick system, but I. uriae is a generalist seabird tick associated with over 50 host species [47][48][49][50][51], including many flighted seabirds, and is therefore likely to have more dispersal opportunities than the more specialised I. eudyptidis and I. kohlsi. Host specialisation can be affected by the phylogenetic and ecological similarities of potential host species, and the extent to which penguin ticks are capable of exploiting other seabirds breeding sympatrically with their primary host remains somewhat unclear [41]. Nonetheless, although I. eudyptidis and I. kohlsi have occasionally been found on a number of seabird species other than the Little Blue Penguin [33,[52][53][54][55], such observations have generally been associated with paralysis in the host, suggesting that these are not the tick species' usual hosts [54].
In order to infer the dispersal capacity of hard ticks (I. eudyptidis) associated with Little Blue Penguins in southeastern Australasia, we used multi-gene molecular approaches to assess phylogeographic structure in the parasites over a range of spatial scales. We hypothesised that, if the parasites are frequently able to disperse across long distances with their hosts (primary or otherwise), no clear phylogeographic structure would be observed.

Ethics statement
All fieldwork undertaken for this project was approved by the Australian National University Animal Experimentation Ethics Committee. Animal Ethics Protocol Number: A2012/57; approved 18 th October 2012.

Sampling and taxonomy
Four sites were selected for the study, and were chosen based on their distance from each other and the size of the host colony. The three Australian sites were Brush Island (35.5292°S, 150.4167°E), Montague Island (36.2500°S, 150.2167°E) in New South Wales and Phillip Island (38.4833°S, 145.2333°E) in Victoria (Fig 1). The samples for the New Zealand site were obtained from Oamaru (45.0842°S, 170.9806°E) on the South Island (Fig 1). Sampling was undertaken during the breeding season of the Little Blue Penguin (November-January) in 2012-2013, and involved briefly taking Little Blue Penguins (usually late-stage chicks, but occasionally younger chicks or adults) from the nest environment, and removing attached ticks with forceps. In some cases, where human-made nest boxes were inhabited by penguins, ticks were also collected from under the lid of the box. Ticks were obtained from 6 nests on Brush Island, 15 on Montague Island, 14 on Phillip Island, and 18 from Oamaru. In some cases, multiple ticks were taken from single burrows (see S1-S4 Tables for details). Ticks were preserved immediately after collection in 96% ethanol. Little Blue Penguins are reported to be parasitised by three hard tick species within their distribution (Ixodes eudyptidis, I. kohlsi and I. uriae), and these can be morphologically separated by the shape of the anal groove [33]. Ticks were therefore identified under a dissecting microscope, primarily using the morphology of the anal groove as described in Roberts (1970) to assess whether they were most likely I. eudyptidis, I. kohlsi or I. uriae (see S1-S4 Tables, Taxonomic Methods and S1 Fig) [33].

DNA sequencing
Tick tissue was removed from the posterior half of each tick, with the anterior retained as a voucher, except in the case of extremely small ticks, where the entire animal was used for extraction. DNA extractions were performed using a standard Chelex procedure [56]. Three genetic markers were amplified: mitochondrial cytochrome c oxidase subunit I (COI), 16S rRNA, and nuclear 28S rRNA, all of which have been shown to be informative for phylogeographic research [57][58][59], particularly COI. Furthermore, all primers specifically amplified the invertebrate DNA, thus avoiding any contamination by the penguin DNA remaining in the ticks. PCR amplifications were carried out in 25 μl volumes, each containing 2.5 μl of DNA (concentration not assessed prior to amplification), 0.5 μM of each primer, 1 x buffer, 0.8 mM of dNTPs, 1.5 mM MgCl 2 and 1 U of EconoTaq DNA Polymerase (Lucigen Corporation, Middleton, Wisconsin, United States of America). See Table 1 for a list of PCR primers used.
Amplification was performed in an Eppendorf Mastercycler (epgradient S, Eppendorf, Hamberg Germany) using the following profile: 94°C for 2 minutes; 40 cycles of 15 s at 94°C, 30 s at the specific annealing temperature (see Table 1), 1 min at 72°C, followed by a final 4 min extension at 72°C. PCR products were then purified using the GE Healthcare illustra GFX PCR DNA and Gel Band Purification Kit, Protocol 5.3 with Elution buffer type 6 (GE Healthcare UK Limited, Buckinghamshire, UK), and sequenced by Macrogen Inc. Standard Sequencing Service (Guman-sugan, Korea). Geneious version 6.1.6 (created by Biomatters, available from http://www.geneious.com/) was used to process, align and check the sequence data.

Analyses
Both Maximum Likelihood (ML) and Bayesian phylogeographic analyses were performed, as consistency between the topologies of two differing approaches enhances confidence in the interpretations of patterns. Both methods are considered appropriate for assessing evolutionary relationships [61]. Maximum Likelihood phylogenetic analyses were carried out using PhyML 3.0 [62] with evolutionary model parameters as estimated by the Akaike Information Criterion (AIC) of jModeltest2 [63,64]. Phylogenetic analyses incorporated published sequences from several congeneric species (see S7 Table for a list of NCBI accession numbers) [65] and Amblyomma americanum as outgroups. ML analyses were then performed with a GTR + I + G model for COI (as selected by jModeltest2. Base frequencies A = 0.27805, C = 0.18587, G = 0.15270, T = 0.38339 and gamma shape parameter: 1.032; proportion of invariant sites: 0.414), and a GTR + G model for 16S (A = 0.4350, C = 0.0620, G = 0.1280, T = 0.3750 and gamma shape parameter = 0.3410). Although a TPM1uf+G model was selected by jModeltest2 for 16S, PhyML does not support this model, so GTR was implemented instead. Analysis of a concatenated dataset of all three markers was carried out with a GTR + I + G model. Support for each node was assessed by bootstrapping, with heuristic analysis of 1000 replicate datasets.
Bayesian phylogenetic analyses were carried out for COI and 16S individually and for a mixed model Bayesian concatenated analysis of 16S, COI and 28S using MrBayes [66] These analyses used the same model parameters as for ML analyses, and Markov Chain Monte Carlo (MCMC) searches were executed with a total of four chains of 5,000,000 generations, with trees sampled every 100 generations. The first 10,000 trees were discarded as burn-in. PartitionFinder v1.1.1 [67,68] was used to test the best-fit models of molecular evolution for all markers. The concatenated analysis was carried out using the model parameters proposed by PartitionFinder and the parameters outlined above for the individual analyses (COI = TrNef+G, nst = 6, G = 1.032; 16S = TPM1uf+G, nst = 6, G = 0.3410; 28S = HKY model, nst = 2).
Unrooted statistical parsimony networks were created with TCS 1.21 [69] at the 95% confidence limit. Haplotype diversity was calculated in GenAlEx 6.5 [70,71] using raw sequence data. Haplotype accumulation rarefaction curves were estimated for all genes, using the vegan package for R [72]. Pairwise distances and diversity indices were calculated as proportion of differing nucleotides (p-distance) using MEGA 5 [73]. Analyses of molecular variance (AMOVA) were conducted in Arlequin v3.5.1.3 [74] to evaluate the relative importance of within and among-population genetic variance in the 16S and COI sequence datasets, using both haplotype frequencies (F ST ) and p-distances between sequences (F ST ), with significance (P < 0.01) based on 1023 random permutations of the data. Hierarchical AMOVAs were conducted for each marker to partition genetic variation between the two major clades (AUST and OAMA) and among the populations within each of these clades. Separate AMOVAs were then used for the AUST and OAMA clades to identify differences in structure among populations within clades. Pairwise F ST and F ST values were also estimated among all populations for both markers with significance (P < 0.01) based on 1023 random permutations of the data in Arlequin v3.5.1.3 [74].

Sequence data
A total of 96 individual ticks (from 53 separate burrows) were sequenced for a 612 base pair (bp) fragment of COI, 95 for a 351 bp fragment of 16S, and 32 for a 768 bp fragment of 28S (see Table 2 and S1-S4 Tables for details of the occurrence of multiple samples taken from one nest box). All unique sequences have been deposited in GenBank (accessions KM488485-KM488532). The most phylogenetically informative marker was COI, with 177 variable sites and 30 haplotypes, followed by 16S with 129 variable sites and 12 haplotypes, and 28S with 62 variable sites and 5 unique sequences. The transition / transversion ratio was 1.28 for COI, 0.52 for 16S and 4.34 for 28S. Initial 28S analysis showed the marker to have extremely low variability, so rather than sequence all samples, analysis was completed for a subset of samples chosen to represent all unique lineages in the most variable marker (COI).

Phylogeographic structure
Phylogenetic trees constructed using ML and Bayesian approaches were largely topologically consistent, with only outgroup relationships differingonly the ML trees are shown in Figs 2 and 3, with ML bootstrap and Bayesian posterior probability (PP) values shown on any branch receiving greater than 50% / 0.5 support, respectively. Due to the lack of variation within 28S (only 5 unique sequences detected), building a phylogenetic tree would have been uninformative, so only network analysis was carried out for this marker. Phylogeographic analyses for COI ( Fig 2) and 16S (Fig 3), and the concatenated analysis (see S2 Fig) revealed two distinct clades. The largest clade, which we refer to from here on as 'AUST', contained most individuals from  Tables for individuals). For both markers the individuals found in each of the two main clades were consistent. As this divergence could indicate the presence of multiple species, each clade was analysed separately for network and AMOVA analyses. All phylogeographic analyses confirmed the monophyly of the samples collected during this study, and that they were distinct from Ixodes uriae sequences available in GenBank (Figs 2 and 3). In addition, both markers identified an individual that was most closely related to the Amblyomma americanum 'outgroup' sequence (see Figs 2 and 3). The OAMA clade showed some phylogeographic structure across all three markers, with Australian and New Zealand grouping separately (see Figs 2 and 3). However, while COI and 16S analyses both indicated that Oamaru and Australian sample groups within the OAMA clade were paraphyletic, the direction of paraphyly differed for each marker.
Hierarchical AMOVA analyses indicated that the majority of genetic variation was partitioned among clades (F CT ) rather than among populations within the clades (F SC ) (see Table 3). However, patterns were more complex in the clade-specific AMOVA analyses. F ST , F ST and clade-specific AMOVA analyses for both COI and 16S (see Table 3) indicated that genetic variation between geographically disparate sites within the OAMA clade was significantly higher than the variation within the sites. Differences were driven by population pairwise F ST  For the network analyses, across all markers (COI, 16S and 28S), the AUST clade samples could not be connected to the OAMA samples at the 95% confidence limit, and for COI, the Australian and New Zealand samples within the OAMA clade could also not be connected (see Fig 4). Diversity was generally higher for COI than for 16S within sites or clades. For COI, the OAMA clade was only found to have marginally lower haplotype diversity than AUST  S5 Table). COI diversity at each population was compared by rarefaction of haplotype diversity for a standard sample size (n = 10) (see Table 4). Estimates differed between sites but were generally higher in Australian populations than the New Zealand site (Oamaru = 5.1, Brush Island = 4.3, Montague Island = 6.7 and Phillip Island = 6.4) (see S3 Fig).
For 16S, the samples from Oamaru were found to exhibit very low haplotype diversity (h = 0.236, n = 22); especially when compared to the less numerous six Australian individuals from the OAMA group (h = 0.278, n = 6) and the OAMA and AUST groups as a whole (h = 0.495, n = 28; h = 0.473, n = 59). Indeed, the small sample size of the Australian OAMA individuals (n = 6) may have affected the precision of these figures, but the high diversity of the Australian samples as a whole is clear. While Brush Island was found to have low haplotype diversity (h = 0.240, n = 15), the site had a small sample size (n = 15). Furthermore, the haplotype diversity calculated for Brush Island (n = 15) is still higher than that calculated for Oamaru despite the higher sample size at the New Zealand site (n = 22). Again, F ST values for 16S indicate population structure within the OAMA clade (F ST = 0.711, P = 0.000) whereas non-significant F ST values for 16S within the AUST clade populations suggest ongoing gene flow (F ST = 0.020, P = 0.351 and F ST = -0.017, P = 0.432) ( Table 3). In addition, pairwise population F ST values for 16S within the OAMA clade indicate almost no gene flow between Oamaru and Montague populations (F ST = 0.937, P = 0.000) (see S6 Table). Rarefaction estimates of 16S allelic diversity (for n = 10) differed between sites but generally remained higher in Australian than New Zealand sites (Oamaru = 1.9, Brush Island = 2.3, Montague Island = 3.6 and Phillip Island = 4.5) (see Table 4, S4 Fig).
We identified only five unique 28S sequences, precluding detailed phylogeographic analysis, but supporting the divisions between the two major clades. For 28S, no division was found between Australian and New Zealand samples in the OAMA clade (only a single unique sequence was recovered from samples in this clade). 28S diversity at each population was compared by rarefaction of haplotype diversity for a standard, albeit small, sample size (n = 3) (see Table 4). Estimates followed the same trend of COI and 16S, where Australian diversity was generally greater than diversity in the New Zealand site (Oamaru = 1, Brush Island = 1, Montague Island = 2.2 and Phillip Island = 2.0) (see S5 Fig).

Divergence between clades
16S uncorrected p distances between OAMA and AUST ranged from 8.7-10%, and the sample that was grouped with the Amblyomma species (Phillip Island Sample 3) was found to be 21-23.1% divergent from the two groups. COI distances between OAMA and AUST individuals ranged from 13.7-17.3%. The divergence between the Oamaru population samples, and the OAMA-clade Australian samples (from Montague and Phillip Islands), was between 4.5-5%. The sample that was grouped with the Amblyomma species (A. americanum) for both markers (Phillip Island Sample 3) was found to be 22.9-25.4% divergent from the two groups for COI.

Morphology
Minor morphological differentiation in the shape of the anal groove was not found to correlate with the major genetic divisions (see S1-S4 Tables). Almost all samples within the highly divergent OAMA and AUST clades (see S1-S4 Tables and S1 Fig for details) were identified as I. eudyptidis based on anal groove morphology. Although some variation of anal groove shape was noted among the samples, these did not appear to correspond to described differences among Ixodes species (see S1 Fig) [33]. No samples collected during this study were found to correspond genetically with available references sequences of I. uriae, and the only other described Ixodes species parasitising E. minor are I. eudyptidis and I. kohlsi. If both of these species were present in our samples, the differences in anal groove shape are clearly not adequate to distinguish among the species, and taxonomic revision may be warranted.

Discussion
The strong phylogenetic divisions identified between ticks collected from Little Blue Penguins in New Zealand and those from Australia suggest that dispersal events between the two countries (a distance of at least 2000 km) have occurred, but that these have been too rare to maintain gene flow and so have led to allopatric diversification. Indeed, the large divergences among clades (consistent across both the presumably maternally-inherited mtDNA markers, and nuclear 28S) suggest the presence of multiple species, despite the apparent morphological concordance of almost all samples with the description of I. eudyptidis (see S1-S4 Tables, Taxo-nomic Methods and S1 Fig). In contrast to the evidence for limited dispersal of penguin ticks across the Tasman Sea, our results show no significant genetic differentiation (within the major clades) among Australian populations, even on considerable geographic scales (e.g. Brush Island versus Phillip Island, separated by about 600 km of coastline). These results suggest that penguin tick populations within Australia have remained well connected, at least on evolutionary time scales.
Rare trans-oceanic dispersal as a driver of evolutionary diversification The genetic differentiation between the New Zealand and Australian samples from within the OAMA clade was considerable (4.5-5% for COI), suggesting long-term separation of the populations and allopatric divergence. Successful trans-Tasman (Australia -New Zealand) dispersal opportunities for the Little Blue Penguin ticks have therefore probably been extremely rare. That we found low diversity among Oamaru samples in the OAMA clade suggests that the population may be relatively young, perhaps reflecting recent colonisation of Oamaru from an Australian source, for example by dispersal of a common ancestor of OAMA and AUST clades to New Zealand, followed by speciation and dispersal of some of these ticks back to Australia. Alternatively, divergence between the clades may have happened elsewhere, for example in Australia, with recent one-way dispersal of some members of the OAMA clade to New Zealand. In either case, we infer that trans-Tasman dispersal of penguin ticks has taken place. There are two main Little Blue Penguin clades in New Zealand [32], both of which are present at the Oamaru site [21]. The more northern clade (north of the Canterbury Bight) is thought to have been in New Zealand for a long time, whereas the southern clade is thought to represent a recent recolonisation of southern New Zealand [24,31,32] possibly even after human arrival [77]. These southern birds are genetically similar to Australian populations, supporting a hypothesis of trans-Tasman movement of Little Blue Penguins and possibly some occasional gene flow among populations [31]. Rare trans-Tasman dispersal of both host and parasite is thus supported by this and previous studies.
Although we interpret the low genetic diversity in Oamaru tick samples to indicate that they may represent a relatively recent dispersal event, these diversity values might also result from sampling bias (with only 18 nests sampled and host genetics not analysed in this study) or from the relatively cool New Zealand climate limiting the effective population size of the parasites (ticks are generally more numerous in warmer weather: [78]). Future research should assess penguin tick diversity on broader spatial scales and with larger sample numbers per site to attempt to resolve these issues.

Dispersal of penguin ticks along the Australian coastline
The negligible phylogeographic structure between disparate AUST tick populations in Australia suggests the parasites are dispersing, possibly by travelling at sea with their hosts, and that this occurs often enough to maintain some level of gene flow. However, tick survival may be generally limited by length of time at sea, with shorter distances and stepping-stone dispersal along the coast more likely to result in successful dispersal events than long trans-oceanic voyages. The capacity of hard ticks to survive long sea journeys remains disputed [53, 55,79]. Some ticks are capable of surviving several months in fresh water under experimental conditions [79]. Whether penguin ticks could remain attached to birds over-wintering at sea for [6][7] months has yet to be tested, though it seems unlikely given the short attachment durations of the parasites [40].
Other hosts may also play a major role in the dispersal of the penguin parasites. Although I. eudyptidis is primarily associated with Little Blue Penguins [80], the species has occasionally been found on some flighted seabirds [54] and thus may be able to disperse aerially (see further discussion below). Future studies should assess host-species specificity by investigating the ticks of both penguins and sympatric flighted seabirds. In addition, future studies should ideally use more rapidly evolving markers (such as microsatellites or SNPs) to shed light on the frequency and extent of I. eudyptidis dispersal.

Dispersal of ticks with penguins or other hosts
Although Little Blue Penguin foraging trips during breeding seasons do not usually exceed 30 km from the colony [25,34,38], failed breeders or young birds prospecting for new nest sites or colonies make longer trips [41,81] and are a likely mechanism for the dispersal processes underlying the patterns of genetic structure in south-eastern Australian populations. During sampling for this study, it was noted that a common location for ticks on adult Little Blue Penguin individuals was inside the ear. Settling in the ear cavity may allow ticks to survive sea journeys, as ears would be protected from seawater, even during dives.
Another possible explanation for dispersal of penguin ticks across the Tasman Sea is that they may be able to travel with other hosts, such as flying seabirds. Although seabird-associated Ixodes species have shown strong host-specificity in some studies [42][43][44][45]82] (but see [41,52]), Ixodes species have been recorded from shearwaters [5,83], which nest sympatrically with Little Blue Penguins [24]. While the primary host of I. eudyptidis is the Little Blue Penguin, the tick species has been recorded in association with 17 species of seabird in New Zealand, including gulls, gannets, and shags [54]. Future research should use genetic techniques to assess whether ticks on flighted seabirds are from the same species and genetic lineages as those on sympatric penguins.

Unrecognised species
A general global mitochondrial DNA mutation rate for arthropods has been calculated at approximately 1.17-2.3% per million years (myr) [84][85][86]. These figures are, however, based on arthropods that are only distantly related to the tick, and the mutation rates of parasitic organisms are believed to be higher than that of their hosts [87], which may mean higher rates for ticks in comparison to other arthropod species. Nonetheless, even when considering the most conservative figure of just under 9% divergence between the two main tick clades in our study, the separation of the OAMA and AUST clades (up to 17% uncorrected p distance) is likely to be an ancient one (several millions of years). Furthermore, that some sites had both OAMA and AUST clade individuals occurring sympatrically suggests that these clades are reproductively isolated, and most probably represent distinct species. Future taxonomic work should assess morphological and other differences between specimens from these clades to identify unrecognised species.

Conclusion
Our results indicate that long-distance oceanic dispersal of the penguin tick I. eudyptidis has occurred, and may be ongoing among Little Blue Penguin breeding colonies on the east coast of Australia. These findings are particularly relevant in light of other molecular biogeographic studies that have indicated that long distance dispersal events have been more important in driving the current distribution patterns of southern biota than traditionally assumed [1,3,88,89]. Our research contributes to the growing body of literature relating to the importance of long-distance dispersal mechanisms as drivers of evolution, and has important implications for the conservation of penguin populations in terms of understanding disease transmission vector dynamics.