Phylogeography , species delimitation and population structure of a Western Australian short-range endemic mite harvestman ( Arachnida : Opiliones : Pettalidae : Karripurcellia )

The mite harvestmen of the genus Karripurcellia Giribet, 2003 are endemic to the tall, wet eucalypt forests of south-western Western Australia, a region known as a hotspot for biodiversity. Currently, there are two accepted species, K. peckorum Giribet, 2003 and K. sierwaldae Giribet, 2003, both with type localities within the Warren National Park. We obtained 65 COI mtDNA sequences from across the entire distributional range of the genus. These sequences, falling into two to three geographically separate groups, probably correspond to two species. Morphologically, all of the studied specimens correspond to K. peckorum, suggesting cryptic speciation within that species. A few common haplotypes occur in more than one population, but most haplotypes are confined to a single population. As a result, populations are genetically differentiated and gene flow after initial colonization appears to be very limited or completely lacking. Our study provides another example of short-range endemism in an invertebrate from the south-western mesic biome.


Introduction
Pettalidae is an old lineage of mite harvestmen (Arachnida, Opiliones, Cyphophthalmi) confined to the circum-Antarctic terranes, persisting on what are nowadays the former landmasses of southern Gondwana (i.e., southern South America, South Africa, Madagascar, Sri Lanka, Australia and New Zealand) (Shear 1980;Giribet 2003;Giribet et al. 2016).Karripurcellia Giribet, 2003 was erected for three species of Pettalidae from Western Australia due to their unique morphology and isolation from other Australian species, with which they may not form a clade (Giribet et al. 2016).The genus is endemic to the Southwest Australian Floristic Region (SWA-FR) or south-western Western Australian biodiversity hotspot (SWWA), a relatively small region on a temper-ate margin of the world's most arid and insular populated continent (Hopper and Gioia 2004;Rix et al. 2015).The region has been considered an island of biodiversity, a relatively wet continental refuge, bordered on two sides by ocean, and isolated by arid lands to the north, northeast, and east (Hopper 1979).Diversification in the region has been driven by a series of key geological and climatic events, including recent aridification and expansion and contraction of mesic environments (Rix et al. 2015).A number of relictual taxa still survive in mesic SWWA, including delicate animals sharing habitat with Karripurcellia, such as spirostreptid millipedes (Edward and Harvey 2010), spiders (Rix et al. 2009;Rix et al. 2010), pseudoscorpions (Harms 2018) and velvet worms (Reid 2002;Sato et al. in press), among others (see review by Rix et al. 2015).
Due to their low vagility, restriction to moist habitats, and usually isolated ranges, pettalids follow a typical distribution of short-range endemics (Jay et al. 2016), species which inhabit a relatively small area of less than 10,000 km 2 (Harvey 2002), as observed in many SWWA terrestrial arthropods (Rix et al. 2015;Harms 2018).The monophyly and validity of Karripurcellia and of one of its species have, however, been questioned by Karaman (2012), who erected the genus Milipurcellia Karaman, 2012 for K. sierwaldae Giribet, 2003, and synonymized K. harveyi Giribet, 2003with K. peckorum Giribet, 2003.Milipurcellia was later synonymized with Karripurcellia in a paper that included sequence data from multiple Karripurcellia specimens from different localities (Giribet et al. 2016).Unfortunately, little information is available about the broader population structure of the species of Karripurcellia, as the type localities of all three named species come from near the town of Pemberton.Here we examine the mitochondrial cytochrome c oxidase subunit I (COI) gene for 65 Karripurcellia specimens collected across a much expanded distribution range (Fig. 1), with the aim of elucidating species diversity and population structure within this relictual genus of Pettalidae.

DNA extraction, PCR amplification and sequencing
All individuals were fixed in absolute ethanol upon collection and stored at -20°C in the collections of the Museum of Comparative Zoology (MCZ) or the Western Australian Museum (WAM) (Table 1).DNA extraction was performed with the Qiagen Blood and Tissue Kit (Qiagen, Valencia, USA) using a single leg per specimen (usually the third right leg).Total DNA was eluted in 50 µl and stored at -20°C.PCR reactions to amplify a fragment of COI had a total volume of 20 µl comprising 0.2 µl of each primer (100 µM each), 0.4 µl dNTPs (0.1 nM; Invitrogen, Carlsbad, USA), 0.1 µl GoTaq PCR Polymerase (Promega, Madison, USA) and 4 µl of associated buffer, 13.1 µl H 2 0 and 3 µl DNA extract.Primers used were LCO1490 (5'-GGT CAA CAA ATC ATA AAG ATA TTG G -3') and HCO2198 (5' -TAA ACT TCA GGG TGA CCA AAA AAT CA -3') (Folmer et al. 1994).The PCR program consisted of an initial denaturation step of 94 °C for 1 min, followed by 38 amplification cycles (94 °C for 1 min, 46 °C for 30 s, 70 °C for 1 min) and a final elongation step of 70 °C for 5 min, and PCRs were run on an Eppendorf Mastercycler Pro (Hamburg, Germany).PCR products were visualized on 1.5 % agarose/TAE gel and treated with ExoSAP-IT (Affymetrix, Santa Clara, USA) to remove residual primers.Sequencing reactions for cleaned PCR products were performed using the BigDye terminator 3.1 (Applied Biosystems, Carlsbad, USA) and one of the PCR primers, respectively, purified with Sephadex (Amersham Biosciences, Amersham, UK) and sequenced on an ABI Prism 3730xl Genetic Analyzer (applied Biosystems, Carlsbad, USA) at the Bauer Core Facility at Harvard University.The raw sequences were assembled and quality checked in GENEIOUS Pro 7.1.7(Biomatters Limited, Auckland, New Zealand).All COI sequences were submitted to GenBank (accession numbers MH133125-MH133188; Table 1).

Alignment, phylogenetic analyses and haplotype network
All sequences were aligned with GENEIOUS Pro 7.1.7.Seven Karripurcellia COI sequences were available on GenBank (accession numbers: KU207393-KU207397, DQ518106, DQ518107).These sequences originated from the same collections as studied herein, and we could not rule out that we resequenced some of these previously sequenced individuals.To avoid analyzing duplicates, we included only one (KU207397 from IZ-134726) of the published sequences in our analyses as this sequence differed by at least one mutation from all our newly generated sequences, suggesting its originality.The other six published sequences were all identical with sequences generated herein from the respective locality.The alignment was translated in MEGA7 (Kumar et al. 2016) to ascertain that no stop codons are present and that no sequence has an excess of amino acid substitutions.
Phylogenetic relationships were reconstructed with Maximum Parsimony in MEGA7 and with Bayesian Inference in MrBayes 3.2.6 (Ronquist et al. 2012) on the CIP-RES Science Gateway (Miller et al. 2010).The maximum parsimony analysis was run with subtree-pruning-regrafting with ten starting trees and 500 bootstrap replicates.For the Bayesian analysis, four runs with six chains were run for 10 6 generations each, sampling every 1200 th tree.Convergence was verified with Tracer 1.7 (Rambaut et al. 2018) and the first 10% of sampled trees were discarded as burnin.A Kimura 2-parameter (K2P) +I model of nucleotide evolution was employed.MEGA7 had determined the T92 +I model as the best-fit model under the BIC criterion, but this is not implemented in MrBayes and the K2P is the most similar model available.In both analyses Aoraki calcarobtusa westlandica Forster, 1952 (EU673667), Austropurcellia clousei Boyer, Baker & Popkin-Hall, 2015 (KJ796956), Rakaia stewartiensis Forster, 1948 (DQ518117) and Neopurcellia salmoni Forster, 1948 (DQ825638) were included as outgroups.To delimit potentially cryptic species, automated barcoding gap discovery (ABGD; Puillandre et al. 2012) was run (simple distance; steps=100; X = 0.1).
A median-joining haplotype network of all sequences was calculated with Network 5 (Fluxus Technology Ltd.) and re-drawn to scale with Adobe Illustrator CC 2017.Pairwise uncorrected p-distances were computed with MEGA7.Population genetic indices, pairwise φ ST (including only localities from which at least four sequences were available) and AMOVA were calculated in Arlequin 3.5.2.2 (Excoffier and Lischer 2010).P-values for pairwise φ ST were Bonferroni corrected to account for multiple comparisons.The AMOVA was run twice, once grouping all samples into two ("northern-central" and "southern") and once into three ("northern", "central" and "southern") groups (see results), respectively.All individuals were included in the AMOVA, with each locality corresponding to a single population, which were assigned to the respective groups.

Results
The 528 bp alignment comprised 65 sequences of Karripurcellia, had 25 variable sites and the translated amino acid alignment featured two amino acid substitutions in total (not considering the outgroup): a change of Thr to Ile in position 136 in three out of nine individuals from IZ-134726; and a change from Ser to Asn (or vice versa) in position 130 between the individuals from the north and south populations.The latter amino acid change is thus diagnostic for these clades.No stop codons were present.
The examined Karripurcellia specimens fall into two or three geographically separated groups (Figs 1, 2, Table 2).One includes all localities sampled south of Northcliffe (herein referred to as "southern"), the second all localities to the north ("northern-central").At least seven mutational steps (1.3 -2.5 % uncorrected p-distance) separate the COI sequences of these two groups (Fig. 1).The localities north of Northcliffe are geographically and genetically further subdivided: the localities between Northcliffe and Pemberton ("central") and the ones north of Pemberton (IZ-134719, IZ-134722, IZ-134723; "northern").However, here the genetic differentiation is less pronounced with as few as two mutational steps (0.2 -1.3 % uncorrected p-distance).It is noteworthy, that the type localities of K. peckorum and K. sierwaldae are located in the area where the "central" group occurs.The differentiation between "southern" and "northern-central" was also recovered in the ABGD analysis for genetic distances between 0.20 -0.39 %, smaller distances grouped each haplotype separate and larger distances collapsed all into a single group.Given the low number of mutations and the overlap of genetic distances (Table 3), it is not surprising that the groups "central" and "northern" were not differentiated in the ABGD analysis.In the phylogenetic analyses, the groups "southern" (bootstrap support [BS] = 96%; posterior probabilities [pp] = 1) and "northern" (BS = 75; pp = 1) were each monophyletic.The groups "central" and "northern-central" were non-monophyletic in both phylogenetic analyses, with members of the "central" group forming a paraphyletic grade with respect to the "northern" specimens.In the   AMOVA analyses, a larger portion of the observed genetic variation was assigned to differences among groups in the analyses with all three groups (80 %; F CT = 0.80) than with two groups (76.5 %; F CT = 0.766) (Table 4).Three haplotypes were shared among some localities, one haplotype each for the "northern", "central" and "southern" populations (Fig. 1).Populations with shared haplotypes were separated at most by 17 km (i.e., IZ-134724 & IZ-134720).However, other populations with similar geographic distances did not share haplotypes (e.g., IZ-134722 & IZ-134725).These three haplotypes were also the dominant haplotypes in their respective groups.However, the majority of haplotypes were recorded at a single locality only.Even between geographically closely associated localities not more than one haplotype is shared (Fig. 1).Pairwise φ ST as well as F SC derived from the AMOVA suggested relatively high and significant levels of genetic differentiation among and between the individuals inhabiting the various localities within each of the groups (Tables 4, 5).Tajima's D as well as Fu's Fs are negative for Karripurcellia for each of the groups.This might indicate population expansion; however, none of these values were statistically significant (Table 2).

Discussion
Only a handful of studies have looked at population structure among Cyphophthalmi, but the few available have focused on pettalids.The New Zealand species Aoraki denticulata is known as an example of an old species with extreme population-level variation, without a single haplotype shared between localities as distant as the ones sampled here and with p-distances of up to 22% for the same marker (Boyer et al. 2007;Fernández and Giribet 2014).While population-level studies are not yet available, Karripurcellia seems to follow a pattern more similar to that of the eastern Australian Austropurcellia (see Jay et al. 2016), whose species are also short-range endemics with distributions rarely spanning more than 50 km and which have been suggested as a model for biogeographic studies of forest corridors and contraction.
The sampled Karripurcellia show low genetic variation but clear structure and suggest the presence of a species boundary between the "northern-central" and "southern" samples in the ABGD analysis, a pattern congruent with both the haplotype network and the phylogenetic trees.The observed variation between the "northern" and "central" clades probably reflects intraspecific genetic variation between populations.Of seven of the herein studied populations the conserved nuclear ribosomal marker 18S rRNA was available on GenBank from a previous study on pettalid phylogeny (Giribet et al. 2016) IZ-134721 [KU207241] from the "southern" group).This marker has been used to differentiate species in other arthropods and onychophorans (e.g., Edgecombe and Giribet 2008;Vélez et al. 2012;Sato et al. in press).Interestingly, the two sequences from the "southern" group and three sequences from the "central" group are identical for this marker, but the only available "northern" sequence (IZ-134719) differed from these in one nucleotide position.This fragment of 18S rRNA is apparently too conserved to allow differentiation between these two species, although the sequence available for IZ-134726 ("central") differed from the other sequences in four nucleotides.
The type localities of K. peckorum and of K. sierwaldae were both reported as from the Warren National Park (which corresponds with the distribution of the "central" group), but no more precise data were available from these specimens collected in 1980.These two species are clearly distinguishable morphologically, among other salient features by the presence of a strong ventral process in the palp trochanter of K. sierwaldae (see Giribet 2003).We included samples from multiple additional localities nearby, including the Dave Evans Bicentennial Tree, in Warren NP.The sample from Bicentennial Tree (GenBank accession DQ518106) is thus identified as K. peckorum, both based on its morphology and location.The third described species, K. harveyi, was collected by pitfall traps in 1976 and the locality was specified as "Crowea, Forest Coupe, c. 12 km S.E. of Pemberton, Table 3. Uncorrected p-distances within and between groups in percent.
"northern" "central" "southern" "northern" 0.0 -0.4 "central" 0.2 -1.3 0.0 -1.1 "southern" 1.7 -2.7 1.3 -2.5 0.0 -0.1 Table 4. AMOVA.In the AMOVA the populations correspond to the individuals collected at each locality, respectively, and the groups correspond to the geographically separated groups A) "northern-central" and "southern" (F SC = 0.554 (p = 0.000), F ST = 0.896 (p = 0.000) and F CT = 0.766 (p = 0.01)) and B) "northern", "central" and "southern" (F SC = 0.367 (p = 0.000), F ST = 0.874 (p = 0.000) and F CT = 0.800 (p = 0.000)).ridge site (34°28′S, 116°10′E)" but the exact locality has not been clearly identified, and seems to be close enough to the K. peckorum localities, confirming the synonymy of K. harveyi proposed by Karaman (2012).Original differences proposed for K. peckorum and K. harveyi are thus probably due to the status of the pitfall-collected specimens, and not species-specific.It is clear that morphologically none of the identified samples correspond to the species K. sierwaldae, known only from Warren NP and from a small adjacent area, Brockman National Park.We expected that the genetically different specimens from locality IZ-134726 would be the sympatric K. sierwaldae, but morphologically these specimens are indistinguishable from K. peckorumlacking clear autapomorphies of K. sierwaldae, such as the ventral process in the palp trochanter.Likewise, the southern specimens resemble K. peckorum morphologically, but these, in agreement with all the sources of genetic information examined here, probably correspond to a different cryptic species of Karripurcellia.

Sum of squares Percentage of variation
The genetic structure within each group and species suggest that dispersal and gene flow are very low or non-existent between populations, even between those separated by only a few kilometers.The only haplotypes shared among populations are the most common haplotypes.All other haplotypes are probably derived from these common haplotypes.Possibly expansion commenced from one "northern -central" and one "southern" population, which featured the respective common haplotype, into the surrounding localities.At each locality new haplotypes arose and there is no indication of any genetic exchange following the initial expansion and colonization.As a result, populations are highly structured and with high levels of endemicity at the species as well as genetic level.Similar distribution patterns of species or intraspecific genetic lineages in this part of the SWWA have been observed for several other taxa from similar localities and with similar life styles, like pseudoscorpions (Harms 2018), trapdoor spiders (Cooper et al. 2011) and onychophorans (Sato et al. in press).In particular, the results of the pseudoscorpion species of Pseudotyrannochthonius (see Harms 2018) are intriguing as these were in part collected at the same sites as the Karripurcellia specimens studied herein.The observed structure could thus be a result of common historic events (e.g., Plio-Pleistocene climatic fluctuations) or geographic barriers limiting dispersal of these taxa with relatively similar lifestyles.For example, the "northern" population is separated from all others by the Donnelley River and the Warren River (though the Warren River partly runs among the "central" populations) and the "northern-central" population is separated from the "southern" population by the Dombakup Brook (which is a tributary to the Warren River).However, the level of divergence was much higher among pseudoscorpions, trapdoor spiders and onychophorans, suggesting that the distribution of Karripurcellia species and their intraspecific genetic lineages is shaped by more recent colonization events.One of the pseudoscorpion spe-cies, Pseudotyrannochthonius giganteus lineage PE/SH (see Harms 2018), and one of the onychophorans (Sato et al. in press) have comparable levels of genetic differentiation across the Dombakup Brook as Karripurcellia.This suggests that physical barriers like rivers do not pose unsurmountable obstacles to dispersal.Highly structured populations and high levels of endemicity at the genetic level have been observed for other Cyphophthalmi (Boyer et al. 2007;Fernández and Giribet 2014) as well; though the level of genetic diversification was much lower in Karripurcellia.Their overall low vagility in combination with priority effects and competitive exclusion (Waters 2011), which limit gene flow between established populations via numerical or fitness advantages of the resident populations, might explain why these species are able to expand into previously uncolonized habitats but exhibit little to no gene flow among established populations.

Conclusion
Our study provides genetic evidence for geographic structuring among the populations of the SWWA endemic genus Karripurcellia, identifying two putative species, K. peckorum and a new "southern" species.Genetic data also suggest the presence of a cryptic species sympatric with K. peckorum that is different from K. sierwaldae, a species that remains elusive and that is only known from the original collections by S. Peck in 1980.Our results also corroborate the synonymy of K. harveyi with K. peckorum previously proposed by Karaman (2012).Future research and sampling in the many unsampled (for Cyphophthalmi) small patches of forest of this biodiversity hotspot should add further clarification to the taxonomic status of these minute and cryptic harvestmen.

Figure 1 .
Figure 1.Distribution of haplotypes.A) Map of sampling localities.B) Median-joining haplotype network.Each mutation is indicated by a dash on the lines connecting haplotypes, mutations resulting in an amino acid changes are marked in red.The size of haplotypes correspond to their observed frequency as indicated by the scale and the colors correspond to the respective localities.

Figure 2 .
Figure 2. Phylogenetic relationships.The depicted phylogenetic tree is based on the maximum likelihood analysis.Bootstrap support from maximum likelihood as well as posterior probabilities from the Bayesian analysis are depicted on the respective branches.Colors correspond to the respective localities (compare Fig. 1).

Table 1 .
Details on individuals and localities.Collection number for the Museum of Comparative Zoology (MCZ) and the Western Australian Museum (WAM) are provided.WA, Western Australia.

Table 2 .
Genetic indices.Standard genetic indices are reported for the whole species as well as for each of the groups.

Table 5 .
Pairwise genetic differentiation.Genetic differentiation was measured with pairwise φ ST .Only populations with at least four sequenced individuals are shown.P-values Bonferroni corrected to account for multiple comparisons.