Sea-level changes and palaeo-ranges: reconstruction of ancient shorelines and river drainages and the phylogeography of the Australian land crayfish Engaeus sericatus Clark (Decapoda: Parastacidae)

Molecular Ecology (2008) 17, 5291–5314

doi: 10.1111/j.1365-294X.2008.03996.x

Sea-level changes and palaeo-ranges: reconstruction of ancient shorelines and river drainages and the phylogeography of the Australian land crayfish Engaeus sericatus Clark (Decapoda: Parastacidae) Blackwell Publishing Ltd

M A R K B . S C H U LT Z ,* D A N I E L A . I E R O D I A C O N O U ,† S A R A H A . S M I T H ,* P I E R R E H O RW I T Z ,‡ A L A S TA I R M . M . R I C H A R D S O N ,§ K E I T H A . C R A N D A L L ¶ and C H R I S T O P H E R M . A U S T I N ** *Arafura Timor Research Facility, School of Environmental and Life Sciences, Charles Darwin University, PO Box 41775, Casuarina, Northern Territory 0811, Australia, †Deakin University, School of Life and Environmental Sciences, PO Box 423, Warrnambool, Victoria 3280, Australia, ‡School of Natural Sciences, Edith Cowan University, 100 Joondalup Drive, Joondalup, Perth, Western Australia 6027, Australia, §School of Zoology, University of Tasmania, Private Bag 5, Hobart, Tasmania 7001, Australia, ¶Department of Biology, Brigham Young University, 675 Widstoe Building, Provo, UT 84602-5181, USA, **School of Environmental and Life Sciences, Charles Darwin University, Darwin, Northern Territory 0909, Australia

Abstract Historical sea levels have been influential in shaping the phylogeography of freshwaterlimited taxa via palaeodrainage and palaeoshoreline connections. In this study, we demonstrate an approach to phylogeographic analysis incorporating historical sea-level information in a nested clade phylogeographic analysis (NCPA) framework, using burrowing freshwater crayfish as the model organism. Our study area focuses on the Bass Strait region of southeastern Australia, which is marine region encompassing a shallow seabed that has emerged as a land bridge during glacial cycles connecting mainland Australia and Tasmania. Bathymetric data were analysed using Geographical Information Systems (GIS) to delineate a palaeodrainage model when the palaeocoastline was 150 m below present-day sea level. Such sea levels occurred at least twice in the past 500 000 years, perhaps more often or of larger magnitude within the last 10 million years, linking Victoria and Tasmania. Inter-locality distance measures confined to the palaeodrainage network were incorporated into an NCPA of crayfish (Engaeus sericatus Clark 1936) mitochondrial 16S rDNA haplotypes. The results were then compared to NCPAs using present-day river drainages and traditional great-circle distance measures. NCPA inferences were cross-examined using frequentist and Bayesian procedures in the context of geomorphological and historical sea-level data. We found distribution of present-day genetic variation in E. sericatus to be partly explained not only by connectivity through palaeodrainages but also via present-day drainages or overland (great circle) routes. We recommend that future studies consider all three of these distance measures, especially for studies of coastally distributed species. Keywords: freshwater crayfish, Geographical Information Systems (GIS), historical biogeography, mitochondrial 16S rDNA, nested clade phylogeographic analysis, palaeodrainages, parastacid Received 1 February 2008; revision received 19 August 2008; accepted 14 October 2008

Introduction Phylogeographic studies using molecular data have provided useful insights into historical biogeographic processes and are relevant to both conservation and evolutionary Correspondence: Christopher M. Austin, Fax: +61-(0)8-8946-6700; E-mail: [email protected] © 2008 The Authors Journal compilation © 2008 Blackwell Publishing Ltd

biology (Avise 2000). For species that display significant phylogeographic structure, a challenge often arises in seeking plausible explanations for the observed geographical patterns of genetic variation. Freshwater-dependent species generally show strong phylogeographic structure; therefore, existing river system and river basin connectivity provide a logical framework for interpreting phylogeographic

5292 M . B . S C H U LT Z E T A L . relationships. There are often strong correlations between contemporary riverine structure and phylogeography of freshwater species (Murphy & Austin 2004); but frequent exceptions imply ancient riverine connections or drainage re-arrangements (Jerry & Woodland 1997; McGlashan & Hughes 2001; Nguyen et al. 2004; Waters et al. 2007). Understanding such discrepancies between phylogeographic structure and contemporary drainage systems (e.g. Burridge et al. 2007) is often limited by an absence of supporting palaeogeographical or palaeoclimatological evidence (see Waters et al. 2007). Australia has an old and ecologically diverse freshwater crayfish fauna, which includes both the world’s largest and one of the world’s smallest species (Crandall et al. 1999, 2000). Of special note are a group of ecologically distinctive land crayfish placed in the genus Engaeus Erichson, 1846, which are characterised by a diminished life-history requirement for surface water compared to other groups of freshwater crayfish (Horwitz & Richardson 1986; Crandall et al. 1999). Engaeus are strongly burrowing freshwater crayfish endemic to the southeast of Australia, principally southern Victoria and northern and western Tasmania (Horwitz 1990, 1994). The genus represents a substantial portion of Australia’s crayfish diversity and contains 35 species, 13 of which are listed as vulnerable or endangered (IUCN 2007). Short-range endemism is a feature of species in the genus (Horwitz 1990). The range of the genus is subject to largescale anthropogenic change and includes zones identified as conservation priorities for freshwater crayfish (Whiting et al. 2000). Southeastern Australia is geographically and geologically complex, characterised by volcanism, uplift, subsidence and sedimentation (e.g. Douglas & Ferguson1976; Bird 1993). The area is a mosaic of lowlands, swamps, plains and highlands, with a vegetation cover ranging from grasslands to temperate rainforest. Since European settlement, just over 200 years ago, the landscape has undergone anthropogenic change resulting in extensive drainage of wetland systems and deforestation for agricultural enterprises (Ierodiaconou et al. 2005). Adding to the complexity of this area are historical sea-level changes driven by glacial cycling, leading to the uniting and disjoining of the Victorian and Tasmanian coastlines across Bass Strait (Horwitz 1988; Bird 1993; Lambeck & Chappell 2001; Unmack 2001). Bass Strait is now a shallow seaway approximately 350 km wide and 500 km long with an average depth of 60 m (Harris et al. 2005). During marine lowstand, Victoria and Tasmania were linked through a terrestrial landscape containing freshwater habitats and during marine highstand, this land bridge was flooded, dividing contemporary mainland Australia and Tasmania (Jennings 1959; Horwitz 1988; Bird 1993; Lambeck & Chappell 2001). Based on present-day bathymetric data, lowstands large enough to link Victoria and Tasmania were likely restricted to within the last 10 million

years. Levels as low as 150 m below the present-day sea level occurred at least twice in the past 500 000 years, perhaps more often and of even larger magnitude within the last 10 million years (Haq et al. 1987; Rohling et al. 1998; Lambeck & Chappell 2001; Rabineau et al. 2006). Highstands up to 20 m above the present-day sea level occurred at least three times in the last 500 000 years (Lambeck & Chappell 2001 and references therein) and up to 80 m between 500 000 and 5.5 million years (Haq et al. 1987). Developments in Geographic Information Systems (GIS) and marine benthic mapping data (e.g. Harris et al. 2005; Ierodiaconou et al. 2007; Heap & Harris 2008) allow the construction of palaeocoastline and corresponding drainage models, enabling tests of hypotheses of faunal contractions and migrations across present-day shallow seas (Lambeck & Chappell 2001). Engaeus sericatus Clark, 1936 and the closely related E. merosetosus Horwitz 1990 are distributed on the northwestern landmass bordering Bass Strait with a largely ‘coastal’ (i.e. on a continental scale) distribution. The taxonomic validity of E. merosetosus is uncertain and will be treated elsewhere (M.B. Schultz et al., unpublished data); thus, for the purposes of this study, we refer to E. sericatus and E. merosetosus as the E. sericatus species complex, or simply as E. sericatus. Given its distribution, a phylogeographic study of this species complex provides an opportunity to investigate the extent to which ancient sea levels have shaped its presentday distribution and, potentially, the distribution of other similarly distributed aquatic species. Engaeus species represent a useful model for investigating the complex historical biogeography of southeastern Australia due to their freshwater-limited coastal distribution on the landmasses surrounding Bass Strait and their expected low dispersal capacity (Horwitz 1988, 1990). If environmental tolerances of coastally distributed Engaeus have remained similar over time, it is likely that historical coastlines and freshwater palaeodrainage networks have significantly influenced phylogeographic patterns within the genus (and within other similarly distributed freshwater-dependent species). These patterns are likely to have included contractions of distributions into higher elevations or protected areas during phases of marine transgression, as well as expansions during marine regressions (see Horwitz 1988). A commonly used analytical tool for the study of phylogeography is nested clade phylogeographic analysis (NCPA; Templeton et al. 1995; Templeton 1998, 2004; Posada et al. 2006; Petit 2008). NCPA seeks to explain population history by examining the geography of population structure and genealogical divergence. Frequently, the NCPA procedure is carried out using the program GeoDis (Posada et al. 2000), which automatically calculates interlocality distances in the form of great-circle distances; however, for freshwaterdependent, riparian or coastal species with constrained dispersal routes (e.g. Engaeus spp.), interlocality distances © 2008 The Authors Journal compilation © 2008 Blackwell Publishing Ltd

C R AY F I S H P H Y L O G E O G R A P H Y A N D PA L A E O D R A I N A G E S 5293 are better described as a matrix of riverine or coastal distances (Posada et al. 2000). Fetzner & Crandall (2003) demonstrated the advantages of using river distances over standard great-circle distances in NCPA but researchers continue to predominantly use great-circle distances for freshwater-restricted species (e.g. Daniels et al. 2006; Gouws et al. 2006; Ponniah & Hughes 2006; Alexander et al. 2007). This probably stems from the fact that river distances are laborious to calculate, or sometimes indeterminable when drainages are connected subsurface or submerged due to sea-level rise (e.g. Buhay & Crandall 2005; Finlay et al. 2006; Buhay et al. 2007). Hitherto, no study (to our knowledge) has used NCPA to examine phylogeographic patterns in the context of GIS-calculated palaeorivers in a palaeolandscape now submerged by marine waters. The novelty of this approach is that bathymetric sea-floor data are employed in addition to the more widely used terrestrial topographical data. The bathymetric data allow reconstruction of an ancient terrestrial landscape during times of lowered sea levels, facilitating the calculation of pairwise interlocality distances through palaeodrainage systems which are now components of the marine seascape. Fetzner & Crandall (2003) and Turner et al. (2000) used river distances within an NCPA framework but did not use marine-inundated palaeodrainage distances. Burridge et al. (2007) investigated the phylogeographic significance of palaeodrainages but did not use river network distances in an NCPA framework. Funk et al. (2005), Hankison & Ptacek (2008) and Vignieri (2005) incorporated interlocality distances as riverine distances, and Spear et al. (2005) used wetland likelihood distances, but none of these studies incorporated the geographical distances into NCPA (see Storfer et al. 2007). And while Kozak et al. (2006) used NCPA to investigate palaeodrainage connections, they used greatcircle distances rather than river distance measures. In this study, the first of a series, we present an approach to estimating river distances for incorporation into phylogeographic analysis of freshwater organisms. Using GIS procedures, we simulate former coastlines and delineate palaeodrainages in the northwest Bass Strait region. We then present a phylogeographic analysis of the burrowing freshwater crayfish E. sericatus species complex as an example to show how different modes of population connectivity can be investigated using a spatial framework. Specifically, we attempt to discriminate between crosscatchment migration, dispersal via present-day drainage connections, and dispersal via reconstructed palaeodrainage systems when coastlines were 150 m below the presentday sea level. The mitochondrial 16S rDNA gene region has been successfully utilised to examine phylogeography, evolution and taxonomy of crayfish and invertebrates (Fetzner & Crandall 2003; Munasinghe et al. 2003; Nguyen et al. 2004; Ponniah & Hughes 2004, 2006; Buhay & Crandall © 2008 The Authors Journal compilation © 2008 Blackwell Publishing Ltd

2005; Buhay et al. 2007; Schultz et al. 2007), and we utilise sequences of a fragment of this gene to: (i) resolve the history of phylogeographic diversification of the Engaeus sericatus species complex in the context of the geomorphology and geological history of the Victorian western district and western Bass Strait using NCPA, and (2) compare and contrast the utility of palaeodrainage, present-day riverine and great-circle overland geographical distance measures as models to explain phylogeographic patterns.

Materials and methods Sampling, laboratory procedures and data collection Specimens of Engaeus sericatus were collected throughout the geographical range of the species complex in southwestern Victoria (see Horwitz 1990), mostly by excavating burrows, but sometimes with a dip-net or bait-trap. Live specimens or pereopod tissue samples were placed in labelled plastic vials or zip-lock bags, chilled on ice during transport and then frozen at –80 °C on return to the laboratory. Tissue samples from private, university and museum collections, preserved in either a mixture of 75% ethanol and 5% glycerol (Horwitz 1990), or 70 to 75% ethanol, were used to supplement field collections. One GenBank sequence was included in the data set. Positional data were recorded at all field sites using hand-held GPS, or locality descriptions were made GIS-compatible through retrospective georeferencing [using Geoscience Australia (www.ga.gov.au), Google Maps (www.maps.google.com) and Google Earth (version 4.2.0198.2451 beta)]. Appendix and Fig. 1 give details of all samples used in this study. Using standard methods, total genomic DNA was isolated from specimens for use as template in polymerase chain reaction (PCR) amplification of the mitochondrial 16S rDNA. PCR products were sequenced directly. Detailed laboratory procedures are as described in Schultz et al. (2007).

Sequence alignment and phylogenetic analysis Raw nucleotide sequence data were edited and assembled in Codoncode Aligner version 2.0.4 (Codoncode Corporation). Taxa were ordered randomly with MacClade version 4.08 (Maddison & Maddison 2005) and aligned with Clustal_X 1.83.1 (Thompson et al. 1997) using default settings (gap opening, 15; gap extension, 6.66; delay divergent sequences, 30%; and DNA transition weight, 0.50). Ambiguously aligned sites were fine-tuned through manual alignment and all nucleotide sites were included in the analysis. Phylogenetic analyses of E. sericatus 16S rDNA haplotype sequences were performed using Bayesian phylogenetic inference, implemented in MrBayes version 3.1.2 (Huelsenbeck

5294 M . B . S C H U LT Z E T A L .

Fig. 1 Collection localities for samples used in this study. Sample codes are described in the Appendix.

& Ronquist 2001; Ronquist & Huelsenbeck 2003). Haplotypes were determined using tcs (see below), treating gaps as missing data. Outgroup taxa were determined from a phylogenetic analysis of mitochondrial 16S rDNA taken from a near-complete taxon sample of Engaeus (M.B. Schultz et al., unpublished data). The outgroup taxa were E. quadrimanus Clark (GenBank Accession no. EU977376), E. fultoni Smith & Schuster (GenBank Accession no. EU977356), E. karnanga Horwitz (GenBank Accession no. EF493105) and E. tayatea Horwitz (GenBank Accession no. EU977388). Model selection and Bayesian analyses were employed as outlined in Schultz et al. (2007). The best-fit model of evolution GTR + I + G was used for the Bayesian analysis. Using four runs and sampling every 1000 cycles, the analysis ran for 2 × 106 generations with a temperature setting of 0.1.

Nested clade phylogeographic analysis Excluding the outgroup taxa, NCPA (Templeton 1993, 1998, 2001, 2004; Templeton et al. 1995) was performed on all 16S rDNA sequences, using the programs tcs version 1.21 (Templeton et al. 1992; Clement et al. 2000), GeoDis

version 2.5 (Posada et al. 2000) and Templeton’s inference key dated 11 November 2005. First, tcs was used to create the haplotype network, treating gaps in the nucleotide alignment as missing data (see Hillis et al. 1996) and using a 95% connection limit. Loops in the network were then resolved, according to the predictions of coalescent theory outlined in Posada & Crandall (2001). Haplotypes in the network were then hierarchically nested into clades (Templeton et al. 1987; Templeton & Sing 1993; Crandall 1996; Templeton 1998, 2002; Panchal 2007). The root of the haplotype network was determined, so that the interior-tip (I-T) status of nested clades could be defined in the GeoDis input file, by recalculating the tcs network in the presence of the outgroup taxa while allowing connections outside of the 95% connection limit. The point of coalescence between the outgroup and in-group was deemed the root (i.e. the most interior, or oldest, node) of the in-group. After determining the root, nested clades showing both genetic and geographical variation were examined by GeoDis (Posada et al. 2006). GeoDis NCPA was run for 1 × 106 random permutations to test the null hypothesis of no association between haplotype distributions and geography (i.e. panmixia). As part of NCPA, GeoDis performs chi-squared (χ2) exact permutational © 2008 The Authors Journal compilation © 2008 Blackwell Publishing Ltd

C R AY F I S H P H Y L O G E O G R A P H Y A N D PA L A E O D R A I N A G E S 5295 contingency tests in which rows are nesting clades and columns are geographical locations. These are the simplest assessments of geographical association of clades as they do not include information on geographical distances between localities (Posada et al. 2000, 2006). A more sophisticated analysis that does incorporate geographical distances between localities is also performed by GeoDis. For that part of the analysis, geographical distances were input as user-defined, GIS-calculated pairwise distances between sample localities (see below for method of distance calculation). The geographical distance NCPA is used to calculate four main statistics for each clade in the nesting hierarchy, with values being either significantly small or significantly large: Dcl, Dnl, I-Tcl and I-Tnl. Significant Dnl values indicate significant interclade pairwise geographical distances, significant Dcl values indicate significant pairwise within-clade distances, and significant I-Tcl or I-Tnl indicate significant values for average interior distance minus average tip distance for both types of distances (Posada et al. 2000, 2006). For any clade with at least one statistically significant Dcl, Dnl, I-Tcl or I-Tnl value, H0 is rejected and Templeton’s inference key is applied (i.e. inferences of demographic history are made).

Statistical correction of false-positive NCPA inferences In the recent literature, there have been a number of criticisms made of NCPA (Knowles & Maddison 2002; Petit & Grivet 2002; Panchal & Beaumont 2007; Petit 2008). Many of these criticisms have been neutralised by followup responses (Garrick et al. 2008; Templeton 2004, 2008); however, debate continues over the elevated rate, causes and means of avoidance of false-positive inferences (Beaumont & Panchal 2008; Templeton 2008). Panchal & Beaumont (2007) stated that the chance of wrongly rejecting H0 and making a false-positive inference (i.e. making a type I error) increases with increasing number of summary statistics within a clade, since various correlations exist among the summary statistics that make them nonindependent. Templeton (2008) responded to this statement by pointing out that each nesting clade yields only a single inference in NCPA, and thus, no multiple test correction is actually needed for tests within clades. Instead, Templeton (2008) suggested that for tests within clades, a Bonferroni-corrected P value should be applied (i.e. α/k) rather than the 5% criterion. This is calculated using α/k, where α is 0.05 and k is the total number of nesting clades analysed by GeoDis. Beaumont & Panchal (2008) responded by stating that correction for multiple tests is not straightforward because each clade yields many statistics, only one of which needs to be significant to make an inference. They show that the relationship between the probability of obtaining at least one significant test statistic in any clade tested and the number of summary statistics © 2008 The Authors Journal compilation © 2008 Blackwell Publishing Ltd

within that clade is not as would be predicted by simple mathematical expectations: Templeton’s (2008) Bonferroni method is conservative when there are less than 14 test statistics within a clade but inadequate when the number of statistics within a clade is greater than or equal to approximately 15 (Panchal & Beaumont 2007; Beaumont & Panchal 2008). In this study, the maximum number of test statistics within any of the clades analysed by GeoDis was 14 (see Results, clade 2-2; Tables 1, 2 and 3), which is within the ‘conservative’ region presented by Beaumont & Panchal (2008). For this particular clade, a Bonferroni-corrected significance cut-off value would be P = 0.05/14 = 0.0036. If we were to apply Templeton’s (2008) Bonferroni correction, the significance cut-off would be even more conservative: a total of 18 clades were analysed by GeoDis, and thus, the cut-off would be P = 0.05/18 = 0.0028. Hence, we applied Templeton’s (2008) Bonferroni method so there was no danger of making false-positive inferences in our analyses. Despite our over-conservative approach to determining significant clade statistics, we chose also to heed the recommendations of Garrick et al. (2008) by employing a range of supplementary procedures. We cross-examined NCPA inferences (after Buhay et al. 2007; Crandall et al. 2008; Garrick et al. 2008) for the total cladogram and the highestlevel clades before the total cladogram as outlined below in the section titled ‘Cross-examinations of NCPA inferences’.

Calculation of geographical distance measures The E. sericatus species complex is distributed in freshwater systems. The species constructs ‘type 2–3’ burrows (connected to the groundwater-table), and is expected to have a low dispersal capacity (Horwitz & Richardson 1986; Horwitz 1988, 1990). Therefore, a present-day wetland–river– coastline network may more realistically represent interlocality distances than pairwise great-circle geographical distances (Posada et al. 2000; Fetzner & Crandall 2003). E. sericatus are distributed within an approximately 150-kmwide band following the southwestern Victorian coastline, on the northwest of the marine divide Bass Strait. This coastline, confining the southern extent of the present-day distribution (Horwitz 1990), has been enormously variable on a geological timescale (Bird 1993) and historical migration Table 1 Mantel test results, testing for correlations between geographical distance matrices

Great circle Present day Palaeodrainage

Great Circle

Present Day

— Rxy = 0.519, P = 0.010* Rxy = 0.566, P = 0.010*

— — Rxy = 0.926, P = 0.010*

*indicates significance at Bonferroni-corrected P < 0.016.

5296 M . B . S C H U LT Z E T A L . Table 2 Results of the NCPA of the Engaeus sericatus-complex 16S haplotypes based on 1 × 106 permutations and using the palaeodrainage network as the pairwise distance matrix between populations Haplotypes Clade

One-step clades Dnl

Dcl

Clade

Dcl

Two-step clades Dnl

Clade

14 (tip) 0 (int) I-T

4.70 203.24 198.53

285.99 215.43 –70.56

1-1 (int)

223.41

218.78

17 (tip) 4 (int) I-T

0.00 196.59 196.59

117.38 184.53 67.16

1-2 (tip)

175.66

198.56

47.74

20.22

I-T 15

1-6

2 (tip) 24 (int) I-T

1-16 (int)

9.10

204.98

38.81 0.00 –38.81

44.49 53.29 8.80

1-3 (tip)

s*46.65

L*323.53

0.00 7.90 7.90

15.87 11.88 –3.98

1-8 (tip)

13.21

225.49

1-9 (tip)

s*2.84

229.35

23 (tip) 22 (int) I-T

0.00 0.00 0.00

23.02 23.02 0.00

1-12 (tip)

23.02

206.41

12 (tip) 21 (int) I-T

0.00 0.00 0.00

68.49 68.49 0.00

1-11 (tip)

68.49

211.91

5 (tip) 9 (int) I-T 1

Dnl

Clade

2-1 (tip)

214.11

210.51

2-4 (int)

0.00

148.36

–214.11

–62.15

I-T 13

Dcl

Three-step clades

I-T –22.01 –53.67 1,19,20,2,3,5,6,7,RGF/dispersal but with some LDD

2-2 (tip)

s*251.84

s*274.20

10

1-7

2-5 (int)

243.55

300.82

3

1-18(int)

2-3 (tip)

s*28.68

L*328.48

25 26

1-4(tip) I-T

25.72

28.65

0.00

28.73

25.72

–0.08

I-T 36.34 1,2,3,5,6,7,RGF/dispersal, with some LDD 8 (tip) 16 (tip) 20 (int) I-T

13.58 0.00 0.00 –10.87

34.96 76.41 58.71 15.46

1-5(int)

45.82

64.55

18 (tip) 11 (int) I-T

0.00 0.00 0.00

142.74 142.74 0.00

1-10(tip)

142.74

115.84

I-T

–96.91

–51.29

1-13(int)

0.00

550.46

1-14(tip)

4.51

352.34

–4.51

198.12

19 6 (tip) 7 (int) I-T

0.00 5.64 5.64

2.82 4.23 1.41

I-T

15.76

2-6 (int)

s*77.37

s*382.63

2-7 (tip)

429.39

L*691.88

I-T –352.01 1,19,20,2,11,12,CRE

s*–309.24

Dcl

Dnl

3-1 (tip)

s*207.13

s*345.86

3-3 (int)

s*287.29

377.54

3-2 (tip)

499.10

L*485.21

4.25

–4.55

I-T 1,2,11,12,CRE

Significantly large (L*) or small (s*) Dcl, Dnl and I-Tcl or I-Tnl values are indicated using a Bonferroni-corrected significance value (P < 0.0028, see text for details). Tip, (tip); interior, (int). Bold-font test statistics indicate significantly small values before Bonferroni correction and bold-font italic statistics © indicate 2008 The significantly large values before correction. Chains of inference are given below I-T values. See caption in Table 4 for abbreviations used in inferences.

Authors Journal compilation © 2008 Blackwell Publishing Ltd

C R AY F I S H P H Y L O G E O G R A P H Y A N D PA L A E O D R A I N A G E S 5297 Table 3 Results of the NCPA of the Engaeus sericatus complex 16S haplotypes based on 1 × 106 permutations and using the present-day river network as the pairwise distance matrix between populations Haplotypes Clade

Dcl

One-step clades Dnl

Clade

Dcl

Two-step clades Dnl

Clade

14 (tip) 0 (int) I-T

7.26 128.28 121.02

151.99 132.30 –19.69

1-1 (int)

134.53

125.86

17 (tip) 4 (int) I-T

0.00 53.14 53.14

37.45 50.75 13.30

1-2 (tip)

48.99

88.76

I-T

85.53

37.10

15

1-6

13

5 (tip) 9 (int) I-T

38.81 0.00 –38.81

44.49 53.29 8.80

0.00 7.90 7.90

15.87 11.88 –3.98

1

1-16 (int)

9.10

122.43

1-3 (tip)

46.65

L*190.51

1-8 (tip)

13.21

125.64

1-9 (tip)

s*2.84

112.46

23 (tip) 22 (int) I-T

0.00 0.00 0.00

10.42 10.42 0.00

1-12 (tip)

s*10.42

120.19

12 (tip) 21 (int) I-T

0.00 0.00 0.00

37.32 37.32 0.00

1-11 (tip)

37.32

122.29

Dnl

Clade

2-1 (tip)

117.30

115.86

2-4 (int)

0.00

90.67

–117.30

–25.19

I-T

2 (tip) 24 (int) I-T

Dcl

Three-step clades

I-T –16.64 –24.05 1,19,20,2,3,5,6,7,RGF/Dispersal but with some LDD

2-2 (tip)

s*143.48

s*163.80

10

1-7

2-5 (int)

168.67

193.84

3

1-18(int)

25.72

28.65

25 26

1-4(tip)

0.00

28.73

25.72

–0.76

2-3 (tip)

s*28.68

L*216.29

I-T

I-T 48.16 1,2,3,5,6,7, RGF/Dispersal but with some LDD 8 (tip) 16 (tip) 20 (int) I-T

13.58 0.00 0.00 –10.87

34.96 76.40 58.71 15.46

1-5(int)

45.83

64.55

18 (tip) 11 (int) I-T

0.00 0.00 0.00

142.74 142.74 0.00

1-10(tip)

142.74

115.84

I-T

–96.91

–51.29

1-13(int)

0.00

444.07

1-14(tip)

4.51

284.64

–4.51

159.43

19 6 (tip) 7 (int) I-T

0.00 5.64 5.64

2.82 4.23 1.41

I-T

19.54

2-6 (int)

s*77.37

s*307.41

2-7 (tip)

346.64

L*541.70

I-T –269.26 1,19,20,2,11,12,CRE

s*-234.28

Dcl

Dnl

3-1 (tip)

s*114.49

s*200.01

3-3 (int)

177.08

220.98

3-2 (tip)

L*395.65

L*346.63

I-T –10.51 –17.15 1,2,3,5,15,PF and/or LDC. Since tcs-network branch lengths are short, colonisation is inferred, perhaps with a recent expansion. Significantly large (L*) or small (s*) D , D and I-T or I-T values are indicated using a Bonferroni-corrected significance value (P < 0.0028, see text for details).

cl nl cl nl © Authors Tip,2008 (tip);The interior, (int). Bold-font test statistics indicate significantly small values before Bonferroni correction and bold-font italic statistics indicate significantly large values before correction. Chains of inference are given below Journal compilation © 2008 Blackwell Publishing Ltd I-T values. See caption in Table 4 for abbreviations used in inferences.

5298 M . B . S C H U LT Z E T A L . routes were probably lost, fragmented or altered, by sealevel change (Horwitz 1988). The closely related Engaeus cunicularius (Erichson), E. laevis (Clark) and Geocharax gracilis Clark all have trans-Bass Strait distributions (Horwitz 1988, 1990; Schultz et al. 2007), and thus, it is very likely that palaeodrainages now submerged by Bass Strait have been influential in shaping freshwater crayfish distributions. Thus, marine-inundated palaeodrainage networks could be a significant predictor of present-day phylogeographic structure. To investigate these possibilities, three NCPAs were performed, using palaeodrainage distances, present-day riverine distances and great-circle distances. This approach extends the methodology of Posada et al. (2000), Turner et al. (2000) and Fetzner & Crandall (2003). Pairwise geographical distances through the reconstructed palaeodrainage network when sea levels were 150 m lower than today were determined for all sample locations (56 localities, resulting in 1596 combinations). The Australian Bathymetry and Topography 9 Arc Second Grid (0.0025 or 250-m grid cell size at the Equator) was used to model the palaeodrainage network (Webster & Petkovic 2005). This provides the highest resolution data set available for the submerged Bass Strait region. The negative 150-m (relative sea level) contour was identified and used to define the palaeocoastline. Palaeodrainages were calculated using the Stream Order extension in ArcGIS 9.2 (ESRI Inc.; http://arcscripts.esri.com). This script uses the elevation data set to determine flow direction and accumulation and identifies the most likely path for palaeodrainages. For comparison, the stream network was also modelled for current land drainages. A high degree of uniformity was observed between the spatial alignment of the predicted and actual land stream networks. Further, the predicted palaeodrainage stream networks were highly concordant with those presented by Harris et al. (2005). This gives confidence in the continuity and connectivity of the modelled palaeodrainage network within the limits of the resolution of the ocean bathymetry data available for model input. An assumption was made that migration between drainage systems would be along the coastline if individuals were to traverse between populations in independent basins (after Horwitz 1988; Fetzner & Crandall 2003). Therefore, the negative 150-m contour, representing the palaeocoastline, was incorporated into the network to allow coastal distance calculations between independent basins. The assumption of migration along the coastline is justified for two reasons. First, coastal regions are generally lower and flatter than other areas of the catchment and, although even more so during times of pre-European settlement, dominated by coastal wetlands. Second, although freshwater crayfish are not found in estuaries, they are likely to have some capacity to tolerate elevated saline conditions over short periods (Mills & Geddes 1980; Pinder et al. 2005). A study of Engaeus fossor (Erichson) and E. cisternarius Suter showed

salinity tolerances of up to 29‰ without mortality, and mortalities of only 13% and 25%, respectively, at the salinity of seawater (~35‰) within 48 h (Suter 1975). A species closely related to Engaeus, Geocharax gracilis (Crandall et al. 1999; Schultz et al. 2007), shows a fair degree of salinity tolerance, with 50% mortality at approximately 35‰ (Kefford et al. 2003). Present-day river distances were calculated between sample localities using a 1:25 000 digital drainage network. Before analysis, the drainage layer used as the network was scrutinised with the addition of 1:100 000 digital elevation data to ensure connectivity of the network data set and to ensure no false holes or sinks were included, which could result in an alternative longer route. The assumption was again made that migration between drainage systems would be along the current coastline. Therefore, a 1:25 000 coastline layer was incorporated into the network to allow coastal distance calculations between sample locations in independent basins. Network calculation for distance between sample locations using the current and palaeodrainage networks were automated using ArcView GIS Network analyst (ESRI Inc.) and Shortest Network Path V.1 extension available at ESRI Arcscripts (http://arcscripts. esri.com). Great-circle distance between sample localities (i.e. the shortest path on the surface of a sphere) was calculated using the Distance and Azimuth Matrix V.2 extension (Jenness 2005) in ArcView GIS 3.3 (ESRI Inc.). Mantel tests for matrix correspondence (Mantel 1967) following the methods of Smouse et al. (1986) and Smouse & Long (1992) were performed with the program GenAlEx version 6 (Peakall & Smouse 2006), using 99 random permutations, to test for correlations between the three geographical distance matrices. We also tested for significant differences between the means of the matrices; however, as the locality set is identical across the three matrices, there was a problem of nonindependence of observations. Therefore, we randomly subsampled the distance matrices by assigning each pairwise comparison a random number between 0 and 1. The matrices were filtered to retain values assigned a random number between 0.5 and 1. Sample sizes (n) of the filtered matrices were: palaeodrainages = 766, present-day riverine = 770, and great-circle distances = 799. anova was used to test for significant differences between the means of the three subsampled matrices. Two-sample t-tests were used with a Bonferroni-corrected cut-off value (P = 0.0167) to perform post hoc comparisons of the means of the three subsampled matrices.

Geographical mapping of nested clade boundaries Based on the results of the nesting design of the tcs haplotype network, the locations of the clades were plotted © 2008 The Authors Journal compilation © 2008 Blackwell Publishing Ltd

C R AY F I S H P H Y L O G E O G R A P H Y A N D PA L A E O D R A I N A G E S 5299 geographically. A tab-delimited text layer was compiled, containing locality names and corresponding clade allocations, and locality coordinates in decimal degrees. Using ESRI ArcGIS 9.2, the three-layer types were assembled and clade boundaries digitised.

Cross-examinations of NCPA inferences (genetic diversity and demography) Using 16S haplotypes, Tajimas’s D (Tajima 1989) statistic was used to examine deviations from neutrality. Significantly negative values of D are consistent with bottlenecks. Mismatch analysis (Rogers & Harpending 1992; Rogers et al. 1996; Harpending et al. 1998) of 16S haplotypes was also performed to look for signatures of range expansion events. Coalescent simulations using 1000 replicates were used to determine a raggedness index of mismatch distributions (Harpending et al. 1993). Coalescent simulations using 1000 replicates were also used to calculate Fu and Li’s F* and D* (Fu & Li 1993), and Fu’s Fs (Fu 1997). Nonsignificant values of F* and D* in combination with significant Fs is indicative of recent population growth (Fu 1997). Nucleotide diversity (π; Nei & Li 1979) was also estimated from 16S rDNA sequences, which is a measure of the average number of nucleotide differences between two sequences that does not depend on sample size (Nei & Kumar 2000). The above tests were performed using the program DnaSP version 4.50.3 (Rozas et al. 2003). Tests for isolation by distance were made using Mantel tests for matrix correspondence, which assessed the correlation between the within-clade pairwise genetic distance matrices and the palaeodrainage, present-day riverine and great-circle geographical distance matrices. Genetic distances were calculated using the maximum-likelihood criterion in paup*, using the HKY + I best-fit model of evolution, with settings: lset base = 0.3418, 0.0877, 0.2168; nst = 2; t ratio = 9.8393; rates = equal; pinvar = 0.7451. Model parameters were selected using the Akaike information criterion (AIC; Akaike 1974) in MrModeltest version 2.0 (Nylander 2004), implemented in paup* 4.0b10 (Swofford 2003). Population dynamics over time were explored using the Bayesian skyline plot model (Drummond et al. 2005), implemented in beast (Drummond & Rambaut 2006, 2007) using the raw 16S rDNA sequence data. The Bayesian skyline plot uses sequence data and a Markov chain Monte Carlo (MCMC) sampling procedure to generate a posterior distribution of effective population size through time (Drummond et al. 2005). The outgroup taxa were excluded from this analysis. The sample sequences were deemed suitable for coalescent analysis because they were from a geographically diverse selection (nearly the entire range of the species), samples were approximately random, population subdivision was not extreme (see Results) and an independently estimated rate of nucleotide substitution is © 2008 The Authors Journal compilation © 2008 Blackwell Publishing Ltd

available (Drummond et al. 2005). Performing eight independent runs, MCMC samplings were run for 2 × 107 cycles, sampling every 1000 cycles. Outputs were analysed using Tracer version 1.4 (Rambaut & Drummond 2008) to ensure that runs were performing correctly and that each run was sampling from the same distribution. Burn-in was 2 × 106 cycles of each run. Analyses employed the best-fit nucleotide substitution model HKY + I (as selected with AIC). Using a strict molecular clock model, the clock (substitution) rate was calibrated using a normally distributed prior, with the mean equal to 0.009 (i.e. 0.9% per million years; Sturmbauer et al. 1996) and the standard deviation of the mean (0.0045) chosen to reflect the 0.5 million years uncertainty in this estimate (after Ho 2007). The underlying constant skyline model used a group size set of 21, which represents one-quarter of the sequences in the data set. upgma was used to construct a starting tree for the analyses and operators were auto-optimised.

Results We amplified an approximately 541 base-pair (bp) fragment of the mitochondrial 3’ 16S rDNA from 85 samples of the Engaeus sericatus species complex. After sequence editing, an alignment of 464 characters (including gaps) was analysed. Sequences are deposited in NCBI GenBank under Accession nos EU313341–EU313404, EU313406–EU313420 and EU313422–EU313424 (Appendix). Of the 464 aligned characters, 428 were constant, and 30 of 36 variable characters were parsimony-informative. Twenty-seven haplotypes were identified (labelled 0 to 26) from the 56 sampled localities. The Bayesian phylogenetic analysis showed strong support for monophyly of all samples but the resolution among haplotypes was low (Fig. 2). One individual was sampled from each of 33 localities. More than two individuals were sampled from each of the remaining 23 localities: two individuals were sampled from each of 19 localities, three individuals were sampled from each of three localities, and five individuals were sampled from one locality. Genetic variation between individuals within localities was minimal: haplotypes were identical within 21 of 23 multisample localities and haplotypes differed by only 1 bp within two of 23 multisample localities (Appendix; Fig. 3). Given the repeated sampling of zero within-locality variability (i.e. in 21/23, or 91.3%, of cases), increasing sample sizes within localities is unlikely to dramatically alter our estimate of within-locality diversity; thus, we expect that our sampling design has captured a good representation of within-locality variability. In pairwise comparisons of 16S sequences from the 85 sampled individuals, the mean number of nucleotide differences was 5.47 ± SE 1.15 and the mean proportion of nucleotide differences (p-distance) was 0.01 ± SE 0.00.

5300 M . B . S C H U LT Z E T A L .

Fig. 2 Phylogeny (majority-rule consensus tree) of Engaeus samples, recovered from Bayesian analysis of 464 nucleotide sites of the mitochondrial 16S rDNA. Numbers above branches are Bayesian posterior probabilities (Pp). Haplotype numbers are described in Appendix I.

The root of the haplotype network (Fig. 3) is an inferred haplotype that is positioned one mutational step from haplotype 10 and three mutational steps from haplotypes 3 and 13. Geographically, this implies that the oldest haplotype should be found in the Gellibrand River with closely related, and perhaps slightly younger, haplotypes to be found in the Curdies River and upper Barwon River. These three rivers are in the central part of the range for the species complex and all rivers have headwater tributaries in the Otways Ranges. A single haplotype network connected all samples with one ambiguous loop that was resolved by applying predictions of coalescent theory (Fig. 3; see Posada & Crandall 2001). Haplotypes separated by up to nine mutational steps had a probability ≥ 0.95 of being connected under parsimony. Nesting of clades resulted in 16 one-step clades, 7 two-step clades, 3 three-step clades and 1 four-step clade (total cladogram). At the highest nesting level before the total cladogram, the network shows three geographically distinct clades (Fig. 4). These correspond to a southwestern clade (clade 3-1), a central clade (3-2) and a southeastern clade (3-3). The southwestern clade (3-1) is geographically widespread and situated primarily over low-lying, swampy, exorheic terrain. It comprises five haplotypes (Fig. 3) that differ by

up to 8 bp. The central clade (3-2) is situated in the higher elevation, endorheic ‘Corangamite lakes district’ and comprises eight haplotypes that differ by up to 7 bp. The southeastern clade (3-3) is situated around the foothills of the mountainous, exorheic Otways district, containing 14 haplotypes that differ by up to 9 bp. If the root of the haplotype network is the oldest haplotype, then the southeastern clade is the oldest clade. Predominantly, the three higher-level clades (southwestern 3-1, central 3-2 and southeastern 3-3) are allopatrically distributed (Fig. 4) but contact zones between the clades do occur at regions of both high and low levels of genetic divergence. Low-level genetic divergences occur between the southwestern and southeastern clades across the Hopkins River/Curdies River catchment-divide, specifically between haplotypes 0 or 15 and haplotype 13 (4 bp difference). Low-level divergences between the southeastern and central clades occur across the Lake Corangamite/upper Barwon River catchment-divide, between haplotypes 20 and 3 (1 bp difference). High-level genetic divergences (8–10 bp divergence) occur between the southwestern clade and the central clade along the northeastern boundary of the southwestern clade 3-1 (Fig. 4). Between the southeastern and central clades, high-level divergence (5 bp) occurs at the southern margin of the central clade 3-2. The nesting algorithm produced 18 clades with both genetic and geographical variation that could be examined by GeoDis. Chi-squared (χ2) exact permutational contingency tests detected significant geographical associations for clades 1-1, 2-1, 2-2, 3-3 and 4-1 (Table 4). Of the three high-level clades, the chi-squared statistic for the southwestern clade (3-1) was nonsignificant.

Comparison of palaeodrainage, present-day and great-circle distance treatments Pairwise interlocality distances ranged from 0.26 to 1025.55 km (mean 387.33 km) for palaeodrainages, 0.26 km to 855.24 km (mean 237.24 km) for present-day rivers, and 0.26 to 240.31 km (mean 68.95 km) for great-circle measures. Mantel tests showed the three distance matrices to be significantly correlated, especially the palaeodrainage and present-day riverine distance matrices (Table 1). anova showed highly significant differences between the means of the three subsampled matrices (Fdf1=2, df2=2332 = 787.35, P = 0.0000). Two-sample t-tests (Bonferroni-corrected cut-off P = 0.0167) showed the that the mean subsampled palaeodrainage distance was 156.60 km greater than the mean subsampled present-day riverine distance (Tdf=1424 = 15.49, P = 0.0000), which in turn, was 169.11 km greater than the mean subsampled great-circle distance (Tdf=860 = 26.98, P = 0.0000). The mean subsampled palaeodrainage distance was 325.71 km greater than the © 2008 The Authors Journal compilation © 2008 Blackwell Publishing Ltd

C R AY F I S H P H Y L O G E O G R A P H Y A N D PA L A E O D R A I N A G E S 5301

Fig. 3 tcs network of 16S rDNA haplotypes showing nesting design as used in the NCPA. Haplotypes (numbers in circles) and numbers of individuals per haplotype are described in Appendix I. Circles indicate haplotypes and circle sizes are proportional to the number of individuals contained within. Larger circles indicate more individuals. Empty circles represent inferred intermediate haplotypes. Rectangles are nesting clades and numbered to reflect hierarchy. For example, ‘1–4’ indicates the fourth level-1 nesting clade, with level-1 being the first nesting category above haplotype. Nesting of clades continues until the total cladogram is reached (level 4). Clades with significant Dcl, Dnl, I-Tcl or I-Tnl values for great-circle, present-day and palaeodrainage distance measures are denoted with *, whereas clades with significant values for great-circle and present-day but not palaeodrainage distance measures are denoted with†; clades 1-1, 2-1, 2-2, 3-3 and 4-1 were significant for the chi-square (χ2) permutational contingency test. The dashed line represents the ambiguous loop, which was removed, and the root of the network is marked with arrows.

mean subsampled great-circle distance (Tdf=816 = 39.70, P = 0.0000). The geographical distance analyses identified similar but not identical geographical associations to those described for the chi-squared tests. The measure of geographical distance and the application of a Bonferronicorrected P value affected the frequency and outcome of inferences. Using the widely accepted significance cut-off of P < 0.05, inferences could be made for six clades for the palaeodrainage distance analysis (Table 2), seven clades for the present-day riverine analysis (Table 3) and seven clades for the great-circle distance analysis (Table 5). But after applying a Bonferroni-corrected P value (Templeton 2008) of 0.0028, the number of inferences was reduced to four for the palaeodrainage analysis, four for the present-day riverine analysis and five for the great-circle distance analysis. There were 15 significant Bonferroni-corrected Dcl, Dnl, I-Tcl or I-Tnl values for palaeodrainage distances, 15 for present-day distances and 19 for great-circle distances, © 2008 The Authors Journal compilation © 2008 Blackwell Publishing Ltd

with more significantly small values for great-circle distances than for present-day riverine and palaeodrainage distances (Table 6). There was one more significantly large Dnl statistic for palaeodrainage and present-day riverine matrices than there was for great-circle distances. With the increase in distance from great-circle to present-day riverine distances, the number of significant values decreased, but with the increase from present-day riverine to palaeodrainage distances, the number of significant values remained constant. Of the three highest-level clades before the total cladogram, the geographical distance analysis showed that, regardless of the distance treatment, there was no significant association between haplotype distributions and their geographical localities for the southwestern clade (3-1); however, significant associations were found for the central (3-2) and the southeastern (3-3) clades. Contiguous range expansion was inferred for the central (3-2) clade, regardless of distance measure. Restricted gene flow or dispersal (with some long-distance dispersal) was inferred

5302 M . B . S C H U LT Z E T A L . Table 4 Summary of clades with significant geographical associations grouped by geographical distance measure Distance measure Palaeodrainage

Present-day

Great circle

Clade

χ2

Probability

Chain of inference

Inference

Haplotypes in 1-1 One-step clades in 2-1 One-step clades in 2-2 Two-step clades in 3-1 Two-step clades 3-2 Two-step clades 3-3 Three-step clades in 4-1 Haplotypes in 1-1 One-step clades in 2-1 One-step clades in 2-2 Two-step clades in 3-1 Two-step clades 3-2 Two-step clades 3-3 Three-step clades in 4-1

27.0000 35.0000 120.0000 37.0000 13.0000 70.0000 170.0000 as above as above as above as above as above as above as above

0.0178* 0.0002* 0.0000* 0.0560 0.2355 0.0003* 0.0000* as above as above as above as above as above as above as above

H0 cannot be rejected H0 cannot be rejected 1, 19, 20, 2, 3, 5, 6, 7 H0 cannot be rejected 1, 19, 20, 2, 11, 12 1, 2, 3, 5, 6, 7 1, 2, 11, 12 H0 cannot be rejected H0 cannot be rejected 1, 19, 20, 2, 3, 5, 6, 7 H0 cannot be rejected 1, 19, 20, 2, 11, 12 1, 2, 3, 5, 6, 7 1, 2, 3, 5, 15

Haplotypes in 1-1 One-step clades in 2-1 One-step clades in 2-2 Two-step clades in 3-1 Two-step clades 3-2 Two-step clades 3-3 Three-step clades in 4-1

as above as above as above as above as above as above as above

as above as above as above as above as above as above as above

1, 19, 20 H0 cannot be rejected 1, 19, 20, 2, 3, 5, 6, 7 H0 cannot be rejected 1, 19, 20, 2, 11, 12 1, 2, 11, 12 1, 2, 3, 4, 9

— — RGF/dispersal but with some LDD — CRE RGF/dispersal, with some LDD CRE — — RGF/dispersal but with some LDD — CRE RGF/dispersal but with some LDD PF and/or LDC. Since tcs -network branch lengths are short, colonisation is inferred, perhaps with a recent expansion IGS — RGF/dispersal but with some LDD — CRE CRE AF

An asterisk (*) indicates significance of the chi-square (χ2) permutational contingency test at the P < 0.05 level. Inferences were made by manually applying the inference key dated 11 November 2005. AF, allopatric fragmentation, CRE, contiguous range expansion; IBD, isolation by distance; IGS, inadequate geographical sampling; I-T, interior-tip; LDC, long-distance colonisation; LDD, long-distance dispersal; RE, range expansion; RGF, restricted gene flow; PF, past fragmentation; SF, subsequent fragmentation.

for the southeastern clade (3-3) using palaeodrainages and present-day riverine distances, but contiguous range expansion was inferred for this clade using great-circle distances. At the total cladogram (4-1), contiguous range expansion was inferred from palaeodrainage distances, long-distance colonisation perhaps with a recent range expansion was inferred from present-day riverine distances, and allopatric fragmentation was inferred from great-circle distances.

Cross-examinations of NCPA inferences Cross-examinations of inferences for the southwestern (3-1), the central (3-2), the southeastern (3-3) and the total cladogram (4-1) helped to clarify demographic histories of the 16S rDNA fragment. For the southwestern clade (3-1), Tajima’s D (–0.8554) was nonsignificant (P > 0.10), which suggests that the clade did not recently experience a bottleneck. The mismatch analysis for the southwestern clade (initial θ = 0, final θ = 1000, expansion parameter τ = 2μt = 3.800) showed a ragged bimodal distribution (Fig. 5). The raggedness

index (r) was 0.1600. Higher values of r are consistent with stationary populations and lower values of r are consistent with expanding populations (Harpending et al. 1993; Harpending 1994). The probability that r simulated from a stationary population (expected) was less than the r observed from the data was high (Pr-expected < r-observed = 0.2043); therefore, r for clade 3-1 is consistent with a stationary population. Fu and Li’s D* (–0.8554) and F* (–0.8993), and Fu’s Fs (–1.8053), were nonsignificant (P > 0.10), which, in combination, do not indicate population growth. Nucleotide diversity (π) was low (0.0026 ± SD 0.0009). Mantel tests did not support isolation by distance in clade 3–1 (R = –0.0670 to 0.0530; P = 0.2300 to 0.4200). These tests were in support of the NCPA, which did not reject H0. For the central clade (3-2), Tajima’s D (–0.3204) was nonsignificant (P > 0.10). The mismatch analysis (initial θ = 0, final θ = 1000, expansion parameter τ = 2μt = 3.607) showed a slightly ragged bi-modal distribution (Fig. 5). However, the raggedness index (r) was significantly low (r = 0.0485; Pr-expected < r-observed = 0.0493), indicating a distribution that is largely consistent with an expanding population. Fu and Li’s D* (–0.2019) and F* (–0.2543) were © 2008 The Authors Journal compilation © 2008 Blackwell Publishing Ltd

C R AY F I S H P H Y L O G E O G R A P H Y A N D PA L A E O D R A I N A G E S 5303 Table 5 Results of the NCPA of the Engaeus sericatus complex 16S haplotypes based on 1 × 106 permutations and using great-circle distances as the pairwise distance matrix between populations Haplotypes Clade 14 (tip) 0 (int) I-T 1,19,20,IGS 17 (tip) 4 (int) I-T

Dcl

One-step clades Dnl

Clade

6.14 s*33.90 27.77

L*85.64 s*40.73 s*−44.91

1-1 (int)

0.00 16.20 16.20

19.83 16.75 –3.08

15

Dcl

Two-step clades Dnl

Clade

45.81

44.45

1-2 (tip)

17.16

35.63

I-T

28.65

8.81

1-6

2 (tip) 24 (int) I-T 5 (tip) 9 (int) I-T

1-16 (int) 20.98 0.00 –20.98

23.36 26.67 3.13

0.00 4.71 4.71

5.62 5.16 –0.45

1

1-3 (tip)

4.55

59.02

s*24.17

L*83.54

1-8 (tip)

5.31

s*50.59

1-9 (tip)

s*2.77

59.22

23 (tip) 22 (int) I-T

0.00 0.00 0.00

2.06 2.06 0.00

1-12 (tip)

s*2.06

58.21

12 (tip) 21 (int) I-T

0.00 0.00 0.00

15.18 15.18 0.00

1-11 (tip)

15.18

55.81

Dnl

Clade

2-1 (tip)

42.41

41.96

2-4 (int)

0.00

33.89

–42.41

–8.07

I-T 13

Dcl

Three-step clades

I-T –7.79 –7.71 1,19,20,2,3,5,6,7,RGF/ Dispersal but with some LDD

2-2 (tip)

L*65.74

L*62.56

10

1-7

2-5 (int)

s*19.40

s*44.72

3

1-18(int)

12.88

14.55

25 26

1-4(tip)

0.00

14.80

12.88

–0.24

2-3 (tip)

s*14.63

49.84

I-T s*–36.12 1,2,11,12,CRE

s*–15.29

I-T

8 (tip) 16 (tip) 20 (int) I-T

9.11 0.00 0.00 –7.28

19.16 39.15 34.73 11.57

1-5(int)

25.09

27.31

18 (tip) 11 (int) I-T

0.00 0.00 0.00

48.35 48.35 0.00

1-10(tip)

48.35

35.08

–23.26

–7.77

1-13(int)

0.00

20.22

1-14(tip)

1.54

13.57

–1.54

6.65

I-T 19 6 (tip) 7 (int) I-T

0.00 1.92 1.92

0.96 1.44 0.48

I-T

2-6 (int)

s*29.25

51.73

2-7 (tip)

s*16.15

62.40

I-T 13.10 1,19,20,2,11,12,CRE

–10.67

Dcl

Dnl

3-1 (tip)

s*41.52

66.61

3-3 (int)

57.82

76.99

3-2 (tip)

55.75

70.48

I-T 1,2,3,4,9,AF

12.60

9.38

Significantly large (L*) or small (s*) Dcl, Dnl and I-Tcl or I-Tnl values are indicated using a Bonferroni-corrected significance value (P < 0.0028, see text for details). Tip, (tip); interior, (int). Bold-font test statistics indicate significantly small values before Bonferroni correction and bold-font italic statistics indicate significantly

© 2008 The Authors large values before correction. Chains of inference are given below I-T values. See caption in Table 4 for abbreviations used in inferences. Journal compilation © 2008 Blackwell Publishing Ltd

5304 M . B . S C H U LT Z E T A L .

Fig. 4 Level 3 clades from the haplotype network (Fig. 3) and stream networks (bold) used to calculate distance matrices. Marine-inundated palaeodrainages are shown in the blue region, which shows bathymetry down to 150 m below present-day sea level. Land with elevation less than +80 m is stippled red to highlight likely areas of historical marine transgression. Red indicates clade 3-1 (the southwestern clade), yellow indicates clade 3-2 (the central clade) and green indicates clade 3-3 (the southeastern clade).

Table 6 Summary of Dcl, Dnl and I-T values for clades showing a significant Bonferroni-corrected P-value (Bonf. P = α/k, where α = 0.05 and k = number of statistics within a nesting clade), grouped by distance measure. (S/L), the number of values that were significantly small and significantly large, respectively. Summaries calculated from Table 2, Table 3 and Table 5

Statistic Dcl Dnl I-T Total

Great circle (S/L)

Present day (S/L)

10 (9/1) 6 (3/3) 3 (3/0) 19 (15/4)

7 (6/1) 7 (3/4) 1 (1/0) 15 (10/5)

Palaeodrainage (S/L) 7 (7/0) 7 (3/4) 1 (1/0) 15 (11/4)

nonsignificant (P > 0.10) and Fu’s Fs (– 4.9581) was highly significant (P = 0.0020), which is a pattern of significance that is consistent with population growth. Nucleotide diversity (π) was 0.0071 ± SD 0.0009, which is larger than

that for clade 3-1. Mantel tests showed that some isolation by distance has occurred in clade 3-2 (R = 0.6870 to 0.7550; P = 0.0100). For this clade, the Mantel test value for R increased with the increase from great-circle to present-day to palaeodrainage distances, meaning larger geographical distance matrices gave stronger correlations with genetic distance. In summary, these statistics support the NCPA inferences of ‘contiguous range expansion’ for the central clade (3-2), perhaps expanding from a source region. For the southeastern clade (3-3), Tajima’s D (–0.9625) was nonsignificant (P > 0.10). The mismatch analysis (initial θ = 0, final θ = 1000, expansion parameter τ = 2μt = 4.593) showed a smooth unimodal distribution (Fig. 5). The raggedness index (r) was significantly low (r = 0.0198; Pr-expected < r-observed = 0.0140), indicating expansion. Fu and Li’s D* –0.3765) and F* (–0.6135) were nonsignificant (P > 0.10) and Fu’s Fs (–11.2842) was highly significant (P = 0.0000), which, in combination, indicate population growth. Nucleotide diversity (π) was 0.0095 ± SD 0.0006, © 2008 The Authors Journal compilation © 2008 Blackwell Publishing Ltd

C R AY F I S H P H Y L O G E O G R A P H Y A N D PA L A E O D R A I N A G E S 5305

Fig. 5 Mismatch distribution for 16S haplotypes in clades 3-1 (southwestern clade), 3-2 (central clade), 3-3 (southeastern clade) and 4-1 (total cladogram). The frequency of pairwise nucleotide differences between samples in each clade is represented by diamonds (blue), whereas the frequency of expected pairwise differences under the expansion model is represented by squares (red).

which was the higher than the nucleotide diversity of clades 3-1 and 3-2. Mantel tests indicated the presence of isolation by distance within clade 3-3 (R = 0.2740 to 0.3310; P = 0.0100). The Mantel test value for R increased with the increase from great-circle to present-day to palaeodrainage distances, meaning larger geographical distance matrices gave stronger correlations with genetic distance. In summary, the combination of these statistics with the NCPA inferences of ‘restricted gene flow or dispersal, with some long-distance dispersal’ (palaeodrainage and present-day riverine), and ‘contiguous range expansion’ (great circle), suggest expansion from a source region, perhaps followed by isolation. At the scale of the total cladogram (4-1), Tajima’s D (–1.24373) was nonsignificant (P > 0.10). The mismatch analysis (initial θ = 0, final θ = 1000, expansion parameter τ = 2μt = 5.869) showed a smooth unimodal distribution (Fig. 5). The raggedness index (r) was significantly low (r = 0.0154; Pr-expected < r-observed = 0.0350), indicating expansion. Fu and Li’s D* (–0.7270) and F* (–1.0481) were nonsignificant (P > 0.10) and Fu’s Fs (–28.7558) was highly significant (P = 0.0000), which, together, are consistent with population growth. Nucleotide diversity (π) was 0.0118 ± SD 0.0005. Mantel tests indicated that some isolation by distance has occurred within the total cladogram © 2008 The Authors Journal compilation © 2008 Blackwell Publishing Ltd

(R = 0.4120 to 0.4440; P = 0.0100). The Mantel test value for R remained relatively constant with the increase from great-circle to present-day to palaeodrainage distances. The mean of the 95% HPD of Bayesian skyline plot (Fig. 6) suggests gradual population expansion since 25 000 years before present; however, the error associated with this estimate is considerable, especially within the time frame of the detected expansion. In summary, the cross-examinations combined with the NCPA inferences of ‘contiguous range expansion’ (palaeodrainage), ‘long-distance colonisation perhaps with recent expansion’ (present-day riverine) and ‘allopatric fragmentation’ (great circle) suggest a complex demographic history at the finer scale with a recent, gradual population expansion at the broad scale.

Discussion Geographical distance measures and phylogeography of Engaeus sericatus species complex Advancements in data extraction techniques using GIS, in particular semi-automated methods to calculate pairwise river network distances, are now feasible for much larger data sets than has previously been possible using less-exact

5306 M . B . S C H U LT Z E T A L . Fig. 6 Bayesian Skyline Plot showing population size as a function of time before the present. The solid black line depicts the mean population size, the grey line depicts the median population size and the light blue lines indicate the upper and lower bounds of the 95% highest posterior density (HPD) interval. The 95% HPD is the shortest interval that contains 95% of all values sampled from the posterior. Population size is in units of effective population size (Ne) multiplied by generation time (tau, τ).

and time-consuming manual approaches. GIS methods can be used to simulate ancient coastlines and infer submerged palaeodrainages from bathymetric data, providing a powerful framework for interpreting contemporary phylogeographic patterns. Fetzner & Crandall (2003) demonstrated the utility of river distances in the phylogeographic study of freshwater-dependent species and recommended that river distances be used in addition to great-circle distances in NCPA. Our study heeds the recommendations of Fetzner & Crandall (2003) but takes the analysis one step further, being the first to incorporate palaeodrainage geographical distance data in an NCPA phylogeographic framework. Harris et al. (2005) provided a broad-scale overview of palaeodrainages on the Australian continental shelf. Our GIS-calculated palaeodrainages are consistent with theirs but our implementation means that we can extract pairwise interlocality distance matrices from the drainage models at any simulated sea level. If palaeodrainages provided favourable freshwater habitat during times of lowered sea level, then our inferred drainages were likely to have influenced the present-day distribution of freshwater crayfish in southwestern Victoria, Australia. Indeed, if catchment boundaries are considered an obstacle to dispersal, which is often the case for freshwaterdependent species (Unmack 2001; Fetzner & Crandall 2003; Murphy & Austin 2004), then coastal samples in the southeastern clade (3-3; occupying the Curdies/Gellibrand and Barwon River systems) show a more obvious connection in the context of palaeodrainages than they do in the context of present-day drainages (Fig. 4). Under a palaeodrainage model, samples either side of Cape Otway are within the same catchment boundary (with only 3-bp difference), but under the present-day drainage model, the easternmost and westernmost coastal samples are separated by 31 coastal catchments. Sea levels as low as

140 m below the present are required to connect the palaeodrainages on either side of Cape Otway. Sea-level oscillations of this magnitude are known to have occurred numerous times in the last 10 million years (Haq et al. 1987; Lambeck & Chappell 2001). Expansion of effective population size of the Engaeus sericatus complex over the last 25 000 years coincides with the most recent postglacial sea level rise in which sea levels rose from a level that was low enough to allow linkage of palaeodrainages around Cape Otway to the present-day level that completely isolates rivers either side of Cape Otway (Lambeck & Chappell 2001). Data from additional loci are needed to further explore this finding as multiple loci are better at recovering rapid successions of fluctuations in population size (A.J. Drummond et al., unpublished data). The palaeodrainage model offers much to explaining the present-day distribution of haplotypes in the southeastern clade (3-3). Schultz et al. (2007) found a similar distribution of the closely related freshwater crayfish species Geocharax gracilis around Cape Otway to that which we found for the southeastern clade (3-3) of the E. sericatus complex in this study. Of particular note is Schultz et al.’s (2007) finding that the Bryant Creek (Gellibrand River, Victoria) G. gracilis and King Island (Tasmania) G. gracilis (Figs 3 and 4, Schultz et al. 2007) were very similar genetically, despite being separated by the Bass Strait. In the context of our palaeodrainage model, the Gellibrand River and river systems on King Island are connected within (i.e. headwaters of) the Cape Otway palaeocatchment. To further explore this observation, we calculated a tcs haplotype network from the Schultz et al. (2007) G. gracilis 16S rDNA sequences, and we rooted the network using the closely related ancestral lineage ‘G. sp. nov. 2’ (Schultz et al. 2007). Interestingly, we found the root, or oldest haplotype, of the G. gracilis lineage was represented by a sample from Kennedy Creek. This creek © 2008 The Authors Journal compilation © 2008 Blackwell Publishing Ltd

C R AY F I S H P H Y L O G E O G R A P H Y A N D PA L A E O D R A I N A G E S 5307 is a tributary of the Gellibrand River. The second oldest haplotype in the G. gracilis network was represented by a sample from the Curdies River. The Gellibrand and Curdies rivers drain from or are geographically proximal to the Otways Ranges, and these rivers are also the inferred root localities of the E. sericatus haplotype network. Spatial concordance of the oldest haplotypes across two genera and three species is consistent with the hypothesis of Horwitz (1988) that refugia in the mountainous Cape Otway regions have played an important role during the evolution of freshwater crayfish in southeastern Australia. The palaeodrainage model does not, however, always provide the more parsimonious explanation of genetic relationships. For example, it does not always reduce the number of catchment boundaries between geographical extremes of the southwestern clade (i.e. coastally in clade 3-1) or of the central clade (3-2), and there are instances where riverine connections between localities are the same for both the present-day and palaeodrainage models. Overland (great-circle) dispersal can also effectively explain some aspects of present-day distributions. For example, within the central clade, two subclades occur: one surrounding Lake Corangamite and another in the northwest (Penshurst) region. Palaeodrainage connections place samples from the northwest over 1000 km from the samples west of Lake Corangamite, and present-day river connections show a distance of more than 850 km. A great-circle connection between these areas is approximately 90 km across flat and historically moist and swampy terrain (Ollier & Joyce 1964). The relatively low genetic divergence and the substantially shorter great-circle distance between crayfish from these two areas indicate that this connection most likely arose through overland expansion or dispersal. Lake Corangamite is the largest permanent lake on the Australian mainland and is the largest permanent saline lake in the country (Williams 1995). In the pre-basaltic period [lower Pleistocene, ~1.6 million years ago (Ma); Ollier & Joyce 1964; Douglas et al. 1976], Lake Corangamite was far more extensive than it is today covering an area of up to 1800 km2. The lake extended westwards to meet the upper reaches of Mount Emu Creek and eastwards to the upper reaches of the Barwon River (Currey 1964). The extent of this large palaeolake provides an explanation for the secondary contact shown between the central clade and the southwestern clade at Mount Emu Creek. Between Early (lower) Pleistocene times to the recent (~1.6 Ma to ~0.01 Ma), volcanism poured out lava over much of the western plains, altering drainage patterns of the lake catchment (Currey 1964; Ollier & Joyce 1964; Jenkin et al. 1976). In the post-basaltic period (< ~0.01 Ma), lake levels receded and connections with Mount Emu Creek ceased; however, flows continued intermittently between Lake Corangamite and the upper reaches of the Barwon River (Currey 1964; Williams 1995). The former lake extent and the more recent © 2008 The Authors Journal compilation © 2008 Blackwell Publishing Ltd

intermittent flows provide an explanation for genetic similarities across the Lake Corangamite/Barwon River catchment-divide. The Otways ranges are much older than, external to, and higher than the neighbouring volcanic province (Ollier & Joyce 1964; Douglas et al. 1976). These factors offer further support to the hypothesis that the Otways Ranges have been an important refuge/locus during the evolution of southwestern Victorian freshwater crayfish. The southwestern clade (3-1) generally follows low-lying terrain from the Portland region in the west to the Curdies River in the east. The submerged palaeodrainages in the Portland region are also across flat terrain. Low-level sister clades nested within the southwestern clade are either widely spread throughout the present-day coastal region or concentrated in the west. Flat terrain and numerous lowlying wetlands in the region have likely facilitated dispersal (see Horwitz 1988). The phylogeographic signal of 16S haplotypes in clade 3-1 is weak, making it difficult to recover the demographic history of this clade; thus, more data are required in this regard. It is noteworthy, however, that the genetic boundary at the eastern edge of the southwestern clade (near the Curdies River) is supported by an almost identical genetic boundary between two members of the freshwater crayfish genus Geocharax (Schultz et al. 2007). Spatial concordances across two genera and three species indicate the presence of a significant barrier to gene flow rather than an apparent barrier detected purely as an artefact of coalescent stochasticity. The barrier has likely been in effect for considerable time; still acting in the present, because within Geocharax, it acts at the interspecific level and in the E. sericatus complex, it acts nearer the intraspecific level.

Modes of dispersal in E. sericatus Although difficult to model within a statistical framework, there remains the possibility that the strong burrowing habit of Engaeus species affords the genus special opportunities in terms of cross-catchment dispersal and habitat exploitation. Engaeus species often construct burrows independent of the water table (Horwitz & Richardson 1986) in a wide range of soil types and drainage environments (Horwitz 1990). Specifically, E. sericatus constructs burrows connected to the water table (type 2) — although sometimes the species is found above the water table — deriving water from both groundwater and surface runoff (Horwitz & Richardson 1986; Horwitz 1990). During times of high rainfall and major episodic flooding events, especially in low-lying swampy environments subject to inundation, the species will emerge from burrows possibly to traverse catchment boundaries of limited relief (Horwitz 1985, 1988). Our data support greater dispersal capacities for this species in flatter, swampier terrain (e.g. around Lake Corangamite in the central clade) and, although not

5308 M . B . S C H U LT Z E T A L . statistically supported by the NCPA or cross-validations, widespread haplotypes in the southwestern clade seem to adhere to this hypothesis of a recent rapid dispersal (see Excoffier & Ray 2008). This contrasts with the more likely mode of dispersal via ancient riverine connections in highly dissected terrain, which is apparent from the distribution of E. sericatus across the western and eastern side of the mountainous Otway Ranges. Dispersal across catchment boundaries has been documented in other freshwater crayfishes (Lodge et al. 2000), such as in species of the genus Cherax Erichson (Nguyen et al. 2004; Gouws et al. 2006) and Paranephrops White (Apte et al. 2007). Euastacus Clark, on the other hand, shows a high degree of drainage endemicity (Ponniah & Hughes 2006), as does Orconectes luteus (Creaser) (Fetzner & Crandall 2003). It is possible that for coastally distributed freshwater crayfish (and other freshwater dependent species), the incidence of cross-catchment dispersal has previously been overemphasised due to the absence of information on marine-inundated palaeodrainage connections (e.g. Gouws et al. 2006; Apte et al. 2007), which may provide equally plausible or more parsimonious explanations for migration pathways. Dispersal via river drainages typically occurs over a much larger geographical distance than via more direct, swampy routes. And changes in sea level occur on a timescale of tens of thousands of years with large-magnitude fluctuations occurring less frequently than small-magnitude fluctuations (Haq et al. 1987; Yokoyama et al. 2000; Lambeck & Chappell 2001; Rabineau et al. 2006). Therefore, it is expected that the time required for dispersal around the mountainous Cape Otway should be far greater than the time required for dispersal across the flatter southwestern and central regions, and the former event is expected to have occurred less often than the latter. Future studies should investigate migration rates of these and similarly distributed taxa, using multiple loci, to explore the period required for migrations via palaeodrainages around Cape Otway. Falling sea levels are expected to leave at least partly saline habitats in their wake, which might be expected to block dispersal in freshwater-distributed taxa. However, the salinity tolerance detected in Engaeus and other closely related freshwater crayfishes (Suter 1975; Mills & Geddes 1980; Kefford et al. 2003; Pinder et al. 2005) implies that short-term, mildly saline environments are not likely to block dispersal in E. sericatus.

Conclusions By quantitatively examining phylogenetic relationships in the context of palaeo- and present-day geography, we are helping to answer the call in the recent literature asking researchers to put the geography into phylogeography (Kidd & Ritchie 2006; Storfer et al. 2007; Kozak et al. 2008). Our phylogeographic study demonstrates a novel application

of GIS for modelling palaeo-ranges and historical dispersal pathways, using as a model a freshwater-limited organism (Engaeus sericatus) that has limited gene flow and a largely coastal distribution in southeastern Australia. Digital elevation and bathymetric data were highly successful for quantitatively reconstructing probable palaeohabitat networks, providing a framework to evaluate genealogical relationships of present-day crayfish populations residing in these networks. The combination of GIS and phylogenetic analysis in this study allowed the testing of a qualitative biogeographic hypothesis proposed by Horwitz (1988), that posits that the disjunct present-day distribution of freshwater crayfish circum-Bass Strait is the result of widespread ancestral lineages being recently fragmented by rising sea levels. This hypothesis was strongly supported by identification of the large southwesterly flowing palaeodrainage system that linked a large part of central west-coast Victoria with northwestern Tasmania and nearby islands, implicating the Otway Ranges as a likely historical refuge area. Our results support the conclusions of Fetzner & Crandall (2003) that the NCPA framework for freshwater species does best to include both great-circle and presentday riverine distances as neither model of population connectivity is likely to solely explain extant genetic variation. However, in addition, we recommend that for studies of coastally distributed freshwater species, marine-inundated palaeodrainage connections need to be considered too as these are likely to have contributed important historical dispersal pathways. In this regard, we found that a single network model did not explain ancestral linkages between Engaeus sericatus samples overall. Instead, we found that a combination of the three models was required. Judging by concordance between the results of our cross-examinations, the Bonferroni-corrected NCPAs of 16S haplotypes performed well in terms of recovering the demographic history of E. sericatus (despite recent criticisms of NCPA). This finding is in agreement with Garrick et al. (2008), who state that NCPA is a useful tool in the phylogeographers toolbox. Our reconstructed coastline and palaeodrainage model can readily be expanded to the study of palaeo-ranges and phylogeography of other species in the Bass Strait region. Linking multiple organism studies within a single spatial framework will likely strengthen biogeographic hypotheses and the power to test biogeographic scenarios (Kidd & Ritchie 2006). The integration of GIS and phylogenetics may also provide important insights into predicting species responses to climate change, helping to identify priority refuge areas for protection. The method in this study is applicable to submerged landscapes in general, such as the Arafura Sea between northern Australia and New Guinea, through southeastern Indonesia, between the North Island and South Island of New Zealand, and southeastern USA. Future studies incorporating multiple species and multiple © 2008 The Authors Journal compilation © 2008 Blackwell Publishing Ltd

C R AY F I S H P H Y L O G E O G R A P H Y A N D PA L A E O D R A I N A G E S 5309 genes (e.g. Burridge et al. 2007) using GIS-based palaeoriver and landscape models offer the prospect of new and exciting approaches to the study of evolutionary and earth history of freshwater species bordering shallow seas.

Acknowledgements This work was supported by ARC Discovery grant number DP0557840 to A.M.M.R., C.M.A., P.H. and K.A.C., and partly funded by Charles Darwin University. For access to museum specimens, we thank Dr Gary Poore and Dr Joanne Taylor (Museum of Victoria). For provision of samples and/or assistance with fieldwork, we thank Sabine L. M. Pircher, Clinton T. Hill, Céleste Jacq, Dr Stuart Linton, Dr Claire F. McClusky, Leon B. Meggs, Dr Adam D. Miller, Dr Hemali Munasinghe, Dr Binh Than Thai and Darren Welsh. Field sampling was carried out under the Department of Sustainability and Environment Victoria Permit Number 10003070 (Dr Sue Hadden). Special thanks to Dr Simon Y. W. Ho for assistance with the Bayesian Skyline Plot analysis. Part of this work was carried out by using the resources of the Computational Biology Service Unit from Cornell University, which is partially funded by Microsoft Corporation. We would also like to thank the teams at the Institute of Medical and Veterinary Studies (Adelaide, South Australia) and Bioscience North Australia (Darwin, Northern Territory Australia) for supporting sequencing. GIS analyses were undertaken at Deakin University, Warrnambool, Victoria, GIS Laboratory facility. We are grateful for the comments offered by three anonymous reviewers, which have helped us to improve this publication.

References Akaike H (1974) A new look at the statistical model identification. IEEE Transactions on Automatic Control, 19, 716–723. Alexander HJ, Taylor JS, Sze-Tsun Wu S, Bredena F (2007) Parallel evolution and vicariance in the guppy (Poecilia reticulata) over multiple spatial and temporal scales. Evolution, 60, 2352– 2369. Apte S, Smith PJ, Wallis GP (2007) Mitochondrial phylogeography of New Zealand freshwater crayfishes, Paranephrops spp. Molecular Ecology, 16, 1897–1908. Avise JC (2000) Phylogeography, the History and Formation of Species. Harvard University Press, Cambridge, Massachusetts. Beaumont MA, Panchal M (2008) On the validity of nested clade phylogeographical analysis. Molecular Ecology, 17, 2563– 2565. Bird ECF (1993) The Coast of Victoria: the Shaping of Scenery. Melbourne University Press, Melbourne, Victoria. Buhay JE, Crandall KA (2005) Subterranean phylogeography of freshwater crayfishes shows extensive gene flow and surprisingly large population sizes. Molecular Ecology, 14, 4259–4273. Buhay JE, Moni G, Mann N, Crandall KA (2007) Molecular taxonomy in the dark: Evolutionary history, phylogeography, and diversity of cave crayfish in the subgenus Aviticambarus, genus Cambarus. Molecular Phylogenetics and Evolution, 42, 435–448. Burridge CP, Craw D, Waters JM (2007) An empirical test of freshwater vicariance via river capture. Molecular Ecology, 16, 1883–1895. Clement M, Posada D, Crandall KA (2000) tcs: a computer program to estimate gene genealogies. Molecular Ecology, 9, 1657–1659. © 2008 The Authors Journal compilation © 2008 Blackwell Publishing Ltd

Crandall KA (1996) Multiple interspecies transmissions of human and simian T-cell leukemia/lymphoma virus type I sequences. Molecular Biology and Evolution, 13, 115–131. Crandall KA, Fetzner JW, Lawler SH, Kinnersley M, Austin CM (1999) Phylogenetic relationships among the Australian and New Zealand genera of freshwater crayfishes (Decapoda: Parastacidae). Australian Journal of Zoology, 47, 199–214. Crandall KA, Harris DJ, Fetzner JW (2000) The monophyletic origin of freshwater crayfish estimated from nuclear and mitochondrial DNA sequences. Proceedings of the Royal Society B: Biological Sciences, 267, 1679–1686. Crandall KA, Robison HW, Buhay JE (2008) Avoidance of extinction through nonexistence: the use of museum specimens and molecular genetics to determine the taxonomic status of an endangered freshwater crayfish. Conservation Genetics, doi: 10.1007/S10592008-9546-9 Currey DT (1964) The former extent of Lake Corangamite. Proceedings of the Royal Society of Victoria, 77, 377–386. Daniels SR, Gouws G, Crandall KA (2006) Phylogeographic patterning in a freshwater crab species (Decapoda: Potamonautidae: Potamonautes) reveals the signature of historical climatic oscillations. Journal of Biogeography, 33, 1538–1549. Douglas J, Ferguson JA (1976) Geology of Victoria. In: Geological Society of Australia Special Publication No. 5. Geological Society of Australia, Melbourne, Victoria. Douglas JG, Spencer-Jones D, Bennett WJ, Wail S, McInnes DW (1976) Victoria: Geological Map. Geological Survey of Victoria, Melbourne, Australia. Drummond AJ, Rambaut A (2006) BEAST Version 1.4. Available online at http://beast.bio.edsac.uk/ [verified October 2007]. Drummond AJ, Rambaut A (2007) beast Bayesian evolutionary analysis by sampling trees. BMC Evolutionary Biology 7. Drummond AJ, Rambaut A, Shapiro B, Pybus OG (2005) Bayesian coalescent inference of past population dynamics from molecular sequences. Molecular Biology and Evolution, 22, 1185–1192. Excoffier L, Ray N (2008) Surfing during population expansions promotes genetic revolutions and structuration. Trends in Ecology & Evolution, 23, 347–351. Fetzner JW, Crandall KA (2003) Linear habitats and the nested clade analysis: an empirical evaluation of geographic versus river distances using an Ozark crayfish (Decapoda: Cambaridae). Evolution, 57, 2101–2118. Finlay JB, Buhay JE, Crandall KA (2006) Surface to subsurface freshwater connections: phylogeographic and habitat analyses of Cambarus tenebrosus, a facultative cave-dwelling crayfish. Animal Conservation, 9, 375–387. Fu YX (1997) Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics, 147, 915–925. Fu YX, Li WH (1993) Statistical tests of neutrality of mutations. Genetics, 133, 693–709. Funk WC, Blouin MS, Corn PS et al. (2005) Population structure of Columbia spotted frogs (Rana luteiventris) is strongly affected by the landscape. Molecular Ecology, 14, 483–496. Garrick RC, Dyer RJ, Behergaray LB, Sunnucks P (2008) Babies and bathwater: a comment on the premature obituary for nested clade phylogeographical analysis. Molecular Ecology, 17, 1401– 1403. Gouws G, Stewart BA, Daniels SR (2006) Phylogeographic structure of a freshwater crayfish (Decapoda: Parastacidae: Cherax preissii) in south-western Australia. Marine and Freshwater Research, 57, 837–848.

5310 M . B . S C H U LT Z E T A L . Hankison SJ, Ptacek MB (2008) Geographical variation of genetic and phenotypic traits in the Mexican sailfin mollies, Poecilia velifera and P. petenensis. Molecular Ecology, 17, 2219–2233. Haq BU, Hardenbol J, Vail PR (1987) Chronology of fluctuating sea levels since the Triassic. Science, 235, 1156–1167. Harpending HC (1994) Signature of ancient population growth in a low-resolution mitochondrial DNA mismatch distribution. Human Biology, 66, 591. Harpending HC, Sherry ST, Rogers AR, Stoneking M (1993) The genetic structure of ancient human populations. Current Anthropology, 34, 483–496. Harpending HC, Batzer MA, Gurven M et al. (1998) Genetic traces of ancient demography. Proceedings of the National Academy of Sciences, USA, 95, 1961–1967. Harris PT, Heap AD, Passlow V et al. (2005) Geomorphic Features of the Continental Margin of Australia. Geoscience Australia Record 2003/30. Geoscience Australia, Canberra, Australia. Heap AD, Harris PT (2008) Geomorphology of the Australian margin and adjacent seafloor. Australian Journal of Earth Sciences, 55, 555–585. Hillis DM, Moritz C, Mable BK (1996) Molecular Systematics, 2nd edn. Sinauer Associates, Sunderland, Massachusetts. Ho SYW (2007) Calibrating molecular estimates of substitution rates and divergence times in birds. Journal of Avian Biology, 38, 409–414. Horwitz P (1985) Aspects of the life history of the burrowing freshwater crayfish Engaeus leptorhyncus at Rattrays Marsh, northeast Tasmania. The Tasmanian Naturalist, 82, 1–5. Horwitz P (1988) Sea-level fluctuations and the distributions of some freshwater crayfishes of the genus Engaeus (Decapoda; Parastacidae) in the Bass Strait area. Australian Journal of Marine and Freshwater Research, 39, 497–502. Horwitz P (1990) A taxonomic revision of species in the freshwater crayfish genus Engaeus Erichson (Decapoda: Parastacidae). Invertebrate Taxonomy, 4, 427–614. Horwitz P (1994) A new species of the freshwater crayfish genus Engaeus Erichson (Decapoda: Parastacidae) from north-western Tasmania. Memoirs of the National Museum of Victoria, 54, 439– 445. Horwitz PHJ, Richardson AMM (1986) An ecological classification of the burrows of Australian freshwater crayfish. Australian Journal of Marine and Freshwater Research, 37, 237–242. Huelsenbeck JP, Ronquist F (2001) MrBayes: Bayesian inference of phylogeny. Bioinformatics, 17, 754–755. Ierodiaconou D, Laurenson L, Leblanc M et al. (2005) The consequences of land use change on nutrient exports: a regional scale assessment in south-west Victoria, Australia. Journal of Environmental Management, 74, 305–316. Ierodiaconou D, Burq S, Reston M, Laurenson L (2007) Marine benthic habitat mapping using multibeam data, georeferenced video and image classification techniques in Victoria, Australia. Journal of Spatial Science, 52, 93–104. IUCN (2007) IUCN Red List. Available online at: http://www.iucnredlist.org [verified October 2007]. Jenkin JJ, Ollier CD, Joyce EB (1976) Chapter 10: Geomorphology. In: Geology of Victoria (eds Douglas JG, Ferguson JA). Geological Society of Australia, Melbourne, Victoria, Australia. Jenness J (2005) Distance Matrix (dist_mat_jen.avx) Extension for ArcView 3.x, Version 2. Available at: http://www.jennessent.com/ arcview/dist_matrix.htm. [verified July 2007]. Jennings JN (1959) The submarine topography of Bass Strait. Proceedings of the Royal Society of Victoria, 71, 49–72.

Jerry DR, Woodland DJ (1997) Electrophoretic evidence for the presence of the undescribed ‘Bellinger’ catfish (Tandanus sp.) (Teleostei: Plotosidae) in four New South Wales mid-northern coastal rivers. Marine and Freshwater Research, 48, 235–240. Kefford BJ, Papas PJ, Nugegoda D (2003) Relative salinity tolerance of macroinvertebrates from the Barwon River, Victoria, Australia. Marine and Freshwater Research, 54, 755–765. Kidd DM, Ritchie MG (2006) Phylogeographic information systems: putting the geography into phylogeography. Journal of Biogeography, 33, 1851–1865. Knowles LL, Maddison WP (2002) Statistical phylogeography. Molecular Ecology, 11, 2623–2625. Kozak KH, Blaine RA, Larson A (2006) Gene lineages and eastern North American palaeodrainage basins: phylogeography and speciation in salamanders of the Eurycea bislineata species complex. Molecular Ecology, 15, 191–207. Kozak KH, Graham CH, Wiens JJ (2008) Integrating GIS-based environmental data into evolutionary biology. Trends in Ecology & Evolution, 23, 141–148. Lambeck K, Chappell J (2001) Sea level change through the last glacial cycle. Science, 292, 679–686. Lodge DM, Taylor CA, Holdich DM, Skurdal J (2000) Nonindigenous crayfishes threaten North American freshwater biodiversity: Lessons from Europe. Fisheries, 25, 7–20. Maddison DR, Maddison WP (2005) MacClade Version 4.08: Analysis of Phylogeny and Character Evolution. Sinauer Associates, Sunderland, Massachusetts. Mantel N (1967) The detection of disease clustering and a generalized regression approach. Cancer Research, 27, 209–220. McGlashan DJ, Hughes JM (2001) Genetic evidence for historical continuity between populations of the Australian freshwater fish Craterocephalus stercusmuscarum (Atherinidae) east and west of the Great Dividing Range. Journal of Fish Biology, 59, 55–67. Mills BJ, Geddes MC (1980) Salinity tolerance and osmoregulation of the Australian freshwater crayfish Cherax destructor Clark (Decapoda: Parastacidae). Marine and Freshwater Research, 31, 667–676. Munasinghe DHN, Murphy NP, Austin CM (2003) Utility of mitochondrial DNA sequences from four gene regions for systematic studies of Australian freshwater crayfish of the genus Cherax (Decapoda: Parastacidae). Journal of Crustacean Biology, 23, 402– 417. Murphy NP, Austin CM (2004) Phylogeography of the widespread Australian freshwater prawn, Macrobrachium australiense (Decapoda, Palaemonidae). Journal of Biogeography, 31, 1065– 1072. Nei M, Kumar S (2000) Molecular Evolution and Phylogenetics. Oxford University Press, New York. Nei M, Li WH (1979) Mathematical model for studying genetic variation in terms of restriction endonucleases. Proceedings of the National Academy of Sciences, USA, 76, 5629–5273. Nguyen TTT, Austin CM, Meewan MM, Schultz MB, Jerry DR (2004) Phylogeography of the freshwater crayfish Cherax destructor Clark (Parastacidae) in inland Australia: historical fragmentation and recent range expansion. Biological Journal of the Linnean Society, 83, 539–550. Nylander JAA (2004) MrModeltest Version 2.2. Program distributed by the author, Evolutionary Biology Centre, Uppsala University, Uppsala, Sweden. Ollier CD, Joyce EB (1964) Volcanic physiography of the western plains of Victoria. Proceedings of the Royal Society of Victoria, 77, 357–376. © 2008 The Authors Journal compilation © 2008 Blackwell Publishing Ltd

C R AY F I S H P H Y L O G E O G R A P H Y A N D PA L A E O D R A I N A G E S 5311 Panchal M (2007) The automation of nested clade phylogeographic analysis. Bioinformatics, 23, 509–510. Panchal M, Beaumont MA (2007) The automation and evaluation of nested clade phylogeographic analysis. Evolution, 61, 1466– 1480. Peakall R, Smouse PE (2006) GenAlEx 6: genetic analysis in Excel. Population genetic software for teaching and research. Molecular Ecology Notes, 6, 288–295. Petit RJ (2008) The coup de grâce for the nested clade phylogeographic analysis? Molecular Ecology, 17, 516–518. Petit RJ, Grivet D (2002) Optimal randomization strategies when testing the existence of a phylogeographic structure. Genetics, 161, 469–471. Pinder AM, Halse SA, McRae JM, Shiel RJ (2005) Occurrence of aquatic invertebrates of the wheatbelt region of Western Australia in relation to salinity. Hydrobiologia, 543, 1–24. Ponniah M, Hughes JM (2004) The evolution of Queensland spiny mountain crayfish of the genus Euastacus. I. Testing vicariance and dispersal with interspecific mitochondrial DNA. Evolution, 58, 1073–1085. Ponniah M, Hughes JM (2006) The evolution of Queensland spiny mountain crayfish of the genus Euastacus. II. Investigating simultaneous vicariance with intraspecific genetic data. Marine and Freshwater Research, 57, 349–362. Posada D, Crandall KA (2001) Intraspecific gene genealogies: trees grafting into networks. Trends in Ecology and Evolution, 16, 37–45. Posada D, Crandall KA, Templeton AR (2000) GeoDis: a program for the cladistic nested analysis of the geographical distribution of genetic haplotypes. Molecular Ecology, 9, 487–488. Posada D, Crandall KA, Templeton AR (2006) Nested clade analysis statistics. Molecular Ecology Notes, 6, 590–593. Rabineau M, Berne S, Olivet J-L et al. (2006) Paleo sea levels reconsidered from direct observation of paleoshoreline position during Glacial Maxima (for the last 500 000 yr). Earth and Planetary Science Letters, 252, 119–137. Rambaut A, Drummond AJ (2008) Tracer Version 1.4: MCMC Trace Analysis Tool. Available online from http://tree.bio.ed.ac.uk/ software/tracer/ [verified 2 March 2008]. Rogers AR, Fraley AE, Bamshad MJ, Watkins WS, Jorde LB (1996) Mitochondrial mismatch analysis is insensitive to the mutational process. Molecular Biology and Evolution, 13, 895– 902. Rogers AR, Harpending H (1992) Population-growth makes waves in the distribution of pairwise genetic-differences. Molecular Biology and Evolution, 9, 552–569. Rohling EJ, Fenton M, Jorissen FJ et al. (1998) Magnitudes of sea-level lowstands of the past 500 000 years. Nature, 394, 162–165. Ronquist F, Huelsenbeck JP (2003) MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics, 19, 1572– 1574. Rozas J, Sanchez-De I, Barrio JC, Messeguer X, Rozas R (2003) DnaSP, DNA polymorphism analyses by the coalescent and other methods. Bioinformatics, 19, 2496–2497. Schultz MB, Smith SA, Richardson AMM et al. (2007) Cryptic diversity in Engaeus Erichson, Geocharax Clark and Gramastacus Riek (Decapoda: Parastacidae) revealed by mitochondrial 16S rDNA sequences. Invertebrate Systematics, 21, 569–587. Smouse PE, Long JC (1992) Matrix correlation-analysis in anthropology and genetics. Yearbook of Physical Anthropology, 35, 187–213. © 2008 The Authors Journal compilation © 2008 Blackwell Publishing Ltd

Smouse PE, Long JC, Sokal RR (1986) Multiple regression and correlation extensions of the Mantel test of matrix correspondence. Systematic Zoology, 35, 627–632. Spear SF, Peterson CR, Matocq MD, Storfer A (2005) Landscape genetics of the blotched tiger salamander (Ambystoma tigrinum melanostictum). Molecular Ecology, 14, 2553–2564. Storfer A, Murphy MA, Evans JS et al. (2007) Putting the ‘landscape’ in landscape genetics. Heredity, 98, 128–142. Sturmbauer C, Levinton JS, Christy J (1996) Molecular phylogeny analysis of fiddler crabs: test of the hypothesis of increasing behavioral complexity in evolution. Proceedings of the National Academy of Sciences, USA, 93, 10855–10857. Suter PJ (1975) The taxonomy and ecology of two Engaeus species from Tasmania: Engaeus fossor Erichson, and Engaeus cisternarius sp. nov. B. Sc. (Hons) Thesis. University of Tasmania, Tasmania, Australia. Swofford DL (2003) PAUP*. Phylogenetic Analysis Using Parsimony (*and Other Methods), Version 4. Sinauer Associates, Sunderland, Massachusetts. Tajima F (1989) Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics, 123, 585–595. Templeton AR (1993) The ‘Eve’ hypotheses: a genetic critique and reanalysis. American Anthropologist, 95, 51–72. Templeton AR (1998) Nested clade analyses of phylogeographic data: testing hypotheses about gene flow and population history. Molecular Ecology, 7, 381–397. Templeton AR (2001) Using phylogeographic analyses of gene trees to test species status and processes. Molecular Ecology, 10, 779–791. Templeton AR (2002) Out of Africa again and again. Nature, 416, 45–51. Templeton AR (2004) Statistical phylogeography: methods of evaluating and minimizing inference errors. Molecular Ecology, 13, 789–809. Templeton AR (2008) Nested clade analysis: an extensively validated method for strong phylogeographic inference. Molecular Ecology, 17, 1877–1880. Templeton AR, Sing CF (1993) A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping. IV. Nested analyses with cladogram uncertainty and recombination. Genetics, 134, 659–669. Templeton AR, Boerwinkle E, Sing CF (1987) A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping. I. Basic theory and an analysis of alcohol dehydrogenase activity in Drosophila. Genetics, 117, 343– 351. Templeton AR, Crandall KA, Sing CF (1992) A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping with DNA sequence data. III. Cladogram estimation. Genetics, 132, 619–633. Templeton AR, Routman E, Phillips CA (1995) Separating population structure from population history: a cladistic analysis of the geographical distribution of mitochondrial DNA haplotypes in the tiger salamander, Ambystoma tigrinum. Genetics, 140, 767– 782. Thompson JD, Gibson TJ, Plewniak F, Jeanmougin F, Higgins DG (1997) The Clustal_X Windows interface: flexible strategies for multiple sequence alignment aided by quality analysis tools. Nucleic Acids Research, 25, 4876–4882. Turner TF, Trexler JC, Harris JL, Haynes JL (2000) Nested cladistic analysis indicates population fragmentation shapes genetic diversity in a freshwater mussel. Genetics, 154, 777–785.

5312 M . B . S C H U LT Z E T A L . Unmack PJ (2001) Biogeography of Australian freshwater fishes. Journal of Biogeography, 28, 1053–1089. Vignieri SN (2005) Streams over mountains: influence of riparian connectivity on gene flow in the Pacific jumping mouse (Zapus trinotatus). Molecular Ecology, 14, 1925–1937. Waters JM, Rowe DL, Apte S et al. (2007) Geological dates and molecular rates: rapid divergence of rivers and their biotas. Systematic Biology, 56, 271–282. Webster MA, Petkovic P (2005) Australian Bathymetry and Topography Grid. Geoscience Australia, Canberra, Australia. Whiting AS, Lawler SH, Horwitz P, Crandall KA (2000) Biogeographic regionalization of Australia: assigning conservation priorities based on endemic freshwater crayfish phylogenetics. Animal Conservation, 3, 155–163. Williams WD (1995) Lake Corangamite, Australia, a permanent saline lake: conservation and management issues. Lakes and Reservoirs: Research and Management, 1, 55–64. Yokoyama Y, Lambeck K, De Deckker P, Johnston P, Fifield LK (2000) Timing of the Last Glacial Maximum from observed sea-level minima. Nature, 406, 713–716.

Mark Schultz has recently completed his PhD research on molecular evolution and biogeography of land yabbies. Daniel Ierodiaconou’s research interests include the application of GIS and remote sensing technologies to of terrestrial and marine ecosystems. He is particularly interested in the development of cross-disciplinary geospatial approaches for ecological applications. Sarah Smith is interested in understanding patterns of evolution and diversity, her work is focused on scincid lizards and parastacid crayfish in Australia, New Caledonia and New Guinea. Pierre Horwitz’ research interests extend to the areas of systems ecology and human health, environmental history, the role of science in decision making, and community-based participatory approaches to natural resource and catchment management issues. Alastair Richardson is interested in the biology of terrestrial and freshwater crustaceans, and the role they play as ecosystem engineers. Keith Crandall’s research interests include phylogenetics (theoretical and empirical), bioinformatics, computational biology, population genetics, crustacean systematics, visual pigment evolution, molecular evolution, conservation genetics, biodiversity and bacterial genetics. Christopher Austin is interested in inverstigating, discovering and understanding of evolutionary relationships and biogeographic patterns in freshwater invertebrates. With a major emphasis on the use of molecular genetic tools.

© 2008 The Authors Journal compilation © 2008 Blackwell Publishing Ltd

C R AY F I S H P H Y L O G E O G R A P H Y A N D PA L A E O D R A I N A G E S 5313

Appendix Specimen collection localities and NCBI GenBank Accession numbers

ID no.

Locality

1 2 3* 4 5 6 7 8

Merri Cutting, Dennington Merri Cutting, Dennington Merri Cutting, Dennington Drysdale Creek, Gordons Lane Drysdale Creek, Gordons Lane Mount Emu Creek, McKinnons Bridge Glenfyne, southwest of Cobden Gnotuk, Lake Colongulac southernmost inflow creek Gnotuk, Lake Colongulac southernmost inflow creek Hopkins River, Deakin Uni Lubra Creek, N of Caramut Lubra Creek, N of Caramut Lubra Creek, N of Caramut Lake Pertobe, Warrnambool Lake Pertobe, Warrnambool Lake Pertobe, Warrnambool Lake Pertobe, Warrnambool Lake Pertobe, Warrnambool Panmure, Mount Emu Creek Panmure, Mount Emu Creek Union Creek, E of Woolsthorpe Moyne R. trib., 6 km E Willatook Moyne R. trib., 6 km E Willatook Warrnambool, stormwater drain at Wollaston Bridge Aire River, Great Ocean Road Aire River, Great Ocean Road Aire River, Great Ocean Road Aire River, Hordern Vale Barwon River, S of Stonehaven Barwon River, S of Stonehaven Cargerie Creek, 12 km E of Meredith Waurn Ponds Creek, approx 2 km N Devon Waurn Ponds Creek, Deakin Uni Waurn Ponds Creek, Deakin Uni Birregurra Creek, Princes Hwy Birregurra Creek, Princes Hwy Colac, Barongarook Creek, N of Princes Hwy Colac, Barongarook Creek, N of Princes Hwy Brucknell Creek, Halls Bridge, S of Laang Brucknell Creek, Halls Bridge, S of Laang Brucknell Creek, S of Panmure Curdies River tributary, E of Brucknell Curdies River tributary, E of Brucknell Hopkins River, under Hopkins Falls Hopkins River, under Hopkins Falls Bryant Creek, W of Devondale Burchetts Creek, 2 km E of Caramut Burchetts Creek, 2 km E of Caramut Caramut wetland Camperdown, N of town

9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

© 2008 The Authors Journal compilation © 2008 Blackwell Publishing Ltd

Sample identification

Haplotype no.

GenBank Accession no.

Latitude, S

Longitude, E

MRC1_1 MRC1_H5 — DRY1_1 DRY1_2 EMU1_1 GLE1_1 GNO1_1

0 0 0 0 0 0 0 0

EU313341 EU313342 AY223713 EU313370 EU313371 EU313375 EU313380 EU313381

–38.374201 –38.374201 –38.374201 –38.148927 –38.148927 –38.215904 –38.404636 –38.210382

142.439743 142.439743 142.439743 142.655159 142.655159 142.985903 143.000567 143.144117

GNO1_2

0

EU313382

–38.210382

143.144117

HOP2_1 LBR1_1 LBR1_3 LBR1_4 LPE1_1 LPE1_H94 LPE1_H96 LPE1_H98 LPE1_H99 PAN1_3 PAN1_4 UNN1_1 WHM1_2 WHM1_3 WNB1_2

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

EU313387 EU313389 EU313390 EU313391 EU313395 EU313396 EU313397 EU313398 EU313399 EU313402 EU313403 EU313410 EU313415 EU313416 EU313419

–38.397727 –37.946080 –37.946080 –37.946080 –38.394481 –38.394481 –38.394481 –38.394481 –38.394481 –38.336897 –38.336897 –38.173601 –38.146768 –38.146768 –38.363497

142.539006 142.505975 142.505975 142.505975 142.476213 142.476213 142.476213 142.476213 142.476213 142.727216 142.727216 142.517928 142.324684 142.324684 142.492342

AIR1_1 AIR1_2 AIR1_3 AIR3_1 BAR1_1 BAR1_2 CRG1_1 WPC1_1 WPC2_1 WPC2_3 BIR1_1 BIR1_2 COL2_1

1 1 1 1 2 2 2 2 2 2 3 3 3

EU313343 EU313344 EU313345 EU313346 EU313347 EU313348 EU313365 EU313420 EF493153 EU313422 EU313349 EU313350 EU313362

–38.764167 –38.764167 –38.764167 –38.801258 –38.135927 –38.135927 –37.845199 –38.194389 –38.203327 –38.203327 –38.301443 –38.301443 –38.339935

143.473056 143.473056 143.473056 143.480266 144.264229 144.264229 143.948766 144.221133 144.302672 144.302672 143.770808 143.770808 143.593355

COL2_2

3

EU313363

–38.339935

143.593355

BRC1_1 BRC1_2 BRC2_1 CUR1_1 CUR1_2 HOP1_1 HOP1_2 BRY1_1 BUR1_1 BUR1_2 CAR1_1 CAM1_H21

4 4 4 4 4 4 4 5 6 7 7 8

EU313351 EU313352 EU313353 EU313366 EU313367 EU313385 EU313386 EU313354 EU313355 EU313356 EU313358 EU313357

–38.392735 –38.392735 –38.388333 –38.469211 –38.469211 –38.334601 –38.334601 –38.652948 –37.968327 –37.968327 –37.963876 –38.215796

142.810371 142.810371 142.777778 142.925623 142.925623 142.617315 142.617315 143.233037 142.538275 142.538275 142.517106 143.147959

5314 M . B . S C H U LT Z E T A L . Appendix Continued

ID no.

Locality

51

Gnotuk, Lake Colongulac middle latitude inflow creek Gnotuk, Lake Colongulac northernmost inflow creek Unnamed creek, 16.5 km S of Lismore Chapple Creek, 5 km W of Chapple Vale Kennedy Creek, just N of Kennedy Creek town, Otways Carlisle River, E of bridge Chapple Vale (wetland 4 km N) Floating Islands Reserve, on S of Princes Hwy Gellibrand River, 3 km E of Gellibrand Loves Creek, Kawarren Picnic Ground Barongarook Creek, Wallace St, Colac Curdies River trib, 3 km W of Timboon Curdies River trib, 3 km E of Curdie Vale Fenton Creek, W of Timboon West Mosquito Creek, 4 km W of Lower Heytesbury Duttons Creek, Portland Wattle Hill Creek, Portland Wattle Hill Creek, Portland Ellerslie, Stony Creek Ellerslie, Stony Creek Gnarkeet Chain of Ponds, 8.5 km SW of Wallinduc Lake Gillear, creek under Buckleys Road Lismore, Browns Waterholes Penshurst, wetland park Penshurst, wetland park Pirron Yallock, 5 km S of Nalingil Port Campbell, 9 km SE of Lower Heytesbury Spring Creek, 2 km W of Lower Heytesbury Spring Creek, 2 km W of Lower Heytesbury Wallaby Creek, 2 km S of Lower Heytesbury Wallaby Creek, 2 km S of Lower Heytesbury Winchelsea, Barwon River, Princes Hwy Winchelsea, Barwon River, Princes Hwy Yan Yan Gurr, 9.5 km W of Birregurra Yan Yan Gurr, 9.5 km W of Birregurra

52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85

Sample identification

Haplotype no.

GenBank Accession no.

Latitude, S

Longitude, E

GNO2_1

8

EU313383

–38.193525

143.160845

GNO3_1

8

EU313384

–38.188359

143.177792

NNA1_1 CHC1_1 KCR1_1

8 9 9

EU313401 EU313359 EU313388

–38.073852 –38.626598 –38.595539

143.238146 143.277179 143.240291

CLS2_1 CPV1_1 FIR1_1 GEL2_1 LOV1_1 COL1_1 CUR2_1 CUR3_1 FEN1_1 MSQ1_2

10 10 10 10 10 11 12 13 13 13

EU313360 EU313364 EU313377 EU313379 EU313394 EU313361 EU313368 EU313369 EU313376 EU313400

–38.557592 –38.611906 –38.349140 –38.532005 –38.480198 –38.341494 –38.479549 –38.523771 –38.506660 –38.555470

143.398160 143.328234 143.419890 143.563363 143.582360 143.595617 142.950586 142.858886 142.899338 142.878710

DUT1_H1 WHC1_1 WHC1_2 ELL1_1 ELL1_2 GCP1_2

14 14 14 15 15 16

EU313372 EU313413 EU313414 EU313373 EU313374 EU313378

–38.378319 –38.341817 –38.341817 –38.160035 –38.160035 –37.909813

141.636475 141.562006 141.562006 142.689850 142.689850 143.447161

LGC1_1 LIS1_1 PEN1_3 PEN1_4 PIY1_1 PTC1_1 SPC1_5 SPC1_6 WAL1_1 WAL1_2 WIN1_1 WIN1_2 YYG1_2 YYG1_4

17 18 19 19 20 21 22 22 23 23 24 24 25 26

EU313392 EU313393 EU313404 EF493125 EU313406 EU313407 EU313408 EU313409 EU313411 EU313412 EU313417 EU313418 EU313423 EU313424

–38.424167 –37.952362 –37.873010 –37.873010 –38.339936 –38.612200 –38.563929 –38.563929 –38.583870 –38.583870 –38.242588 –38.242588 –38.328669 –38.328669

142.583889 143.347721 142.290512 142.290512 143.457309 142.993094 142.901673 142.901673 142.916886 142.916886 143.992018 143.992018 143.869343 143.869343

All specimens are from Victoria, Australia. †Sample obtained from Museum Victoria, Melbourne, Australia; ‡sample obtained from Australian Museum, Sydney, Australia. *Sequences obtained from NCBI GenBank, not generated by this study. **Sequence from M.B. Schultz et al. (unpublished). Specimens with ID number ‘—’ were not included in the phylogeographic study, but were used in the Bayesian analysis.

© 2008 The Authors Journal compilation © 2008 Blackwell Publishing Ltd