Sample collection and molecular methods
Tissue samples were collected from 1982 to 2004, from 31 locations throughout the species' range in southeastern Australia (n = 279; Table 1, Fig. 1), by live-trapping (e.g., Browning et al. 2001; Eldridge et al. 2004; Hazlitt et al. 2006) and opportunistically (e.g., road kill). Total cellular DNA was extracted from frozen and alcohol-preserved tissue using standard techniques (Sambrook et al. 1989). Samples were genotyped using 11 polymorphic microsatellite loci. Six loci were derived from the allied rock-wallaby (P. assimilis: Pa55, Pa297, Pa385, Pa593, Pa595, Pa597; Spencer et al. 1995) and five from the tammar wallaby (Macropus eugenii: Me2, Me14, Me15, Me16, Me17; Taylor and Cooper 1998). Individual genotypes were detected using either α33-P labeling visualized by autoradiography (described in Spencer et al. 1995; Taylor and Cooper 1998) or using fluorescent labeling resolved using an automated Amersham Biosciences MegaBACE 500 capillary sequencer (described in Hazlitt et al. 2006). Fifteen to 35 individuals per locus were scored using both methods to ensure consistency, with an overall genotyping scoring error of 1.3% estimated between methods (further details described in Hazlitt et al. 2006).
Table 1. Brush-tailed rock-wallaby sample sites, sample sizes, numbers of individuals sequenced for the control region and numbers of mitochondrial DNA haplotypes detected. Site numbers correspond to those in Figure 1.
|1||Yarraman Creek (YC-Q)||1||1+||1|| AY040890 |
|2||Cooyar Creek (CC-Q)||4||4||1|| AY040890 |
|3||Nukinenda Falls (NF-Q)||1||1+||1|| AY040889 |
|4||Sommerset Dam (SD-Q)||1||1+||1|| AY040890 |
|5||Crows Nest NP (CN-Q)||12||12||2||KJ396276, AY040890|
|6||Perseverance Dam (PD-Q)||16||11||2||KJ396276, KJ396277|
|7||Emu Creek (EC-Q)||10||7||2||EU887006, EU887009|
|8||Farm Creek (FC-Q)||11||9||1|| EU887005 |
|9||Farm Creek East (FCE-Q)||10||4||3||EU887005, EU887010, EU887011|
|10||Hurdle Creek (HC-Q)||54||47||5||EU887004-EU887008|
| ||New South Wales|
|11||Bonalbo (Bon-N)||1||1+||1|| AF357277 |
|12||Armidale (Arm-N)||12||12||2||AF357279, AY040887|
|13||Warrumbungles (War-N)||1||1+||1|| AY040884 |
|14||Woko National Park-2 (Wo2-N)||2||2b||1|| KJ396285 |
|15||Woko National Park-1 (Wo1-N)||2||2b||1|| KJ396284 |
|16||Martindale (Mar-N)||7||7||2||KJ396280, KJ396281|
|17||Yellow Rock (YR-N)||1||1+||1|| KJ396283 |
|18||Drews Creek (DC-N)||20||20||1|| AF357281 |
|19||Ingles Road (IR-N)||29||28||1|| AF357282 |
|20||Bowmans Road (BR-N)||20||20||1|| AF357282 |
|21||St Albans (StA-N)||8||8||3||KJ396278, KJ396279, KJ396282|
|22||Winmalee (Win-N)||1||1+||1|| AY040886 |
|23||Jenolan Caves (Jen-N)||30||30||1|| AF348699 |
|24||Taralga (Tar-N)||2||2*||1|| AF357280 |
|25||Kangaroo Valley (KV-N)||4||4||1|| AF357278 |
|26||Rocky Plains Creek (RPC-V)||8||8||1|| AF357272 |
|27||Little River Gorge (LRG-V)||4||3a||2||AF3572723, AF357276|
|28||Farm Creek (FC-V)||2||2a||1|| AF3572723 |
|29||Gelantipy Creek (GC-V)||2||2a||1|| AF357275 |
|30||Currie Creek (CC-V)||1||1+||1|| AF357274 |
|31||Grampians (Gra-V)||2||2*||1|| AF357271 |
Mitochondrial DNA control region (CR) was amplified using conserved marsupial primers (L15999M and H16498M Fumagalli et al. 1997), and individuals assigned to haplotypes using SSCP (as previously described Sunnucks et al. 2000; Eldridge et al. 2001a) (Table 1). Sequence data were obtained for each unique haplotype using BigDye termination chemistry and resolved using automated capillary sequencers. Over 500 base pairs (bp) of CR sequence was obtained from 3 to 11 individuals of each SSCP haplotype except where a unique haplotype was identified in only one or two individuals. Some mtDNA CR data were also available from previous studies (see Table 1 for GenBank Accession Numbers). In addition, sequences from four other rock-wallaby species were included: two from P. herberti the putative sister species (Potter et al. 2012a)(AF357284 and AY040892), one from P. assimilis (a northerly species from the same species complex), and one sequence from each of P. lateralis (AF348694) and P. purpureicollis (AY057377) for use as outgroups.
Population genetic structure inferred from mtDNA
The program MEGA v5 (Tamura et al. 2007) was used to check sequences and create an initial alignment (using default parameters in ClustalW), which was then adjusted by eye. Phylogenetic relationships among unique CR haplotypes were reconstructed using Bayesian methods implemented in the program MrBayes v3.1.2 (Ronquist and Huelsenbeck 2003). For this analysis, a GTR+I+G model was selected using Modeltest 3.06 (Posada and Crandall 1998). Indels were included as a second data partition within the same analysis, coded as binary data with a variable ascertainment bias. Rate variation within this second partition was initially modeled using a gamma distribution, but an examination of the posterior distribution indicated the data were not informative with respect to this parameter, so in the final model, rates for this data partition were set to equal. In all MrBayes analyses, four chains per run and two independent runs were used. A temperature setting of 0.2 and run length of 6,000,000 generations allowed adequate mixing among chains and convergence between runs. Parameters and trees were sampled every 1000 generations with tree topology and node support assessed over the final 500,000 generations. Convergence between runs, convergence of parameters, and appropriate levels of chain swapping were assessed using Tracer v1.4 (Rambaut and Drummond 2007).
We used Bayes factors, estimated using twice the difference in the natural log of the harmonic mean of model likelihoods of each model (2∆lnHML), assessed following Kass and Raftery (1995), to assess the applicability of a molecular clock. Models run under the three strict clock models implemented in MrBayes (uniform, birth-death and coalescent) were compared with a nonclock model. In all cases, the data strongly supported the use of a molecular clock (all 2∆lnHML > 47.2). The software BEAST v1.4.8 (Drummond and Rambaut 2007) was used to assess node ages, using only the substitution data matrix. The topology of the haplotype tree, with multiple clusters of closely related sequences at the ends of long branches, suggested that these data are at the interface between data best-modeled using coalescent demographic processes and those using speciation processes. Applying a single demographic model across the whole tree would have been inappropriate, as each cluster represented an independently evolving unit, and similarly, models which only incorporate lineage speciation/extinction rates would have been inappropriate at shallower levels in the tree. Hence, we analyzed these data using a Bayesian skyline coalescent model, allowing population size to fluctuate over the tree, therefore eliminating the constraints any particular demographic/speciation model would impose. Bayes factor comparison supported a relaxed lognormal over a strict clock (2∆lnHML = 26.96). All individuals, rather than just a single representative of each haplotype, were entered into the BEAST analyses to allow a more realistic approximation of the coalescent model of sequence divergence and hence lineage divergence times. Two independent runs of each model, each of 12,000,000 generations sampled every 2000 generations, were performed with parameter estimates based on the final 2,000,000 generations. Because estimates of a clock rate for mammalian CR vary widely (6–38% pairwise divergence per million years: Troy et al. 2001; Savolainen et al. 2004; Saarma et al. 2007), we used a mid-range rate of 15% pairwise divergence per million years (Birungi and Arctander 2000); however, given this variation, the resulting time estimates should be treated with caution. The median age and 95% highest posterior probability bounds for major nodes were calculated in Tracer v1.4 and then mapped onto the MrBayes CR haplotype consensus tree.
The hypothesis that three differentiated geographic groups exist within brush-tailed rock-wallabies was tested using the hierarchical analysis of genetic differentiation in ARLEQUIN 3.1 (AMOVA Excoffier et al. 2005). Population pairwise ΦST was calculated using the HKY+G distance model, described below, with significance evaluated with 10,100 permutations. Few haplotypes were shared among populations; therefore, standard FST values (based only on haplotype frequency) would have been unlikely to provide reliable estimates of genetic differentiation. Some of the 31 populations were excluded or pooled with neighboring sites for the AMOVA, as the ΦST distance estimates derived from “populations” with very small sample sizes of only one or two samples were likely to be poor estimates of the true values. Populations with sample sizes less than four were excluded (Yarraman Creek, Warrumbungles, Winmalee, Yellow Rock, Bonalbo, Sommerset Dam, Nukinenda Falls, Currie Creek) unless the small sample size was reflective of the census size (Grampians, Taralga). Populations with n < 4 but within 5 km or less of a neighboring site were pooled (Woko National Park sites 1 and 2, Farm Creek – VIC, Little River Gorge and Gelantipy Creek), resulting in 20 populations (Table 1).
Isolation by distance (IBD) was tested across the entire range and within two of the identified phylogeographic groups (Northern and Central lineage). Sampling for the Southern lineage was not adequate for an IBD analysis. We compared geographic distances (km) and corrected average pairwise population differences, calculated in ARLEQUIN, using Mantel tests with 10,000 permutations in GenAlEx v6.1 (Peakall and Smouse 2006). Pairwise differences were based on a HKY (Tamura in ARLEQUIN) model of haplotype distance with a gamma distribution value of 0.136 (accounting for +G, determined in MODELTEST). Similar results were obtained if alternative distance models were used. Evidence for historic demographic expansion events was tested by mismatch analysis (sum of squared deviations (SSD) and Harpending's raggedness index (R)), as well as Tajimas's D and Fu's Fs tests of neutrality in ARLEQUIN, examining all haplotypes from across the species' range and within each of the three major lineages (Southern, Central, Northern).
Population genetic structure inferred from microsatellites
Microsatellite analyses were conducted using data from 14 sites (n = 247; Table 2). These data included 13 colonies (individuals inhabiting a discreet habitat patch) with sample sizes of n ≥ 7 (Table 2), and three colonies with smaller sample sizes (n ≤ 4) that were merged together to create a single population because they were within a 5 km distance of each other (Little River Gorge, Gelantipy Creek, and Farm Creek, Victoria, n = 8; Table 2), as dispersal has been detected over this distance (Eldridge et al. 2001ab). Exact tests for deviations from Hardy–Weinberg equilibrium for each locus and linkage disequilibrium between loci were carried out for each population in GENEPOP 3.1 (Raymond and Rousset 1995) using the Markov chain method with 1000 iterations. When performing multiple comparisons, we adjusted the statistical significance level using the sequential Bonferroni procedure at α = 0.05 (Rice 1989). Observed and expected heterozygosity (HO and HE) for all loci were estimated using the program POPGENE 1.3.2 and allelic richness (AR, the average number of alleles per locus standardized for unequal sample sizes between sites) was calculated for each sampled colony using FSTAT version 2.9 (Goudet 1995).
Table 2. Genetic diversity (mean ± SE) at 11 microsatellite loci in 14 brush-tailed rock-wallaby populations from southeastern Australia. See Supplementary Table S1 for population allele frequencies.
|5||Crows Nest NP (CN-Q)||12||3.6 (0.53)||2.2 (0.17)||0.52 (0.06)||0.52 (0.07)|
|6||Perseverance Dam (PD-Q)||16||3.7 (0.41)||2.3 (0.13)||0.52 (0.07)||0.57 (0.05)|
|7||Emu Creek (EC-Q)||10||4.4 (0.56)||2.5 (0.14)||0.67 (0.06)||0.64 (0.05)|
|8||Farm Creek (FC-Q)||11||4.5 (0.43)||2.4 (0.17)||0.66 (0.07)||0.61 (0.04)|
|9||Farm Creek east (FCE-Q)||10||3.9 (0.56)||2.5 (0.12)||0.73 (0.09)||0.58 (0.07)|
|10||Hurdle Creek (HC-Q)||54||6.1 (0.61)||2.5 (0.11)||0.68 (0.05)||0.66 (0.04)|
| ||New South Wales|
|12||Armidale (Arm-N)||12||3.9 (0.37)||2.5 (0.14)||0.71 (0.07)||0.62 (0.05)|
|16||Martindale (Mar-N)||7||4.4 (0.41)||2.7 (0.10)||0.71 (0.05)||0.68 (0.03)|
|18||Drews Creek (DC-N)||20||4.3 (0.51)||2.4 (0.13)||0.64 (0.06)||0.60 (0.05)|
|19||Ingles Road (IR-N)||29||5.2 (0.35)||2.6 (0.08)||0.75 (0.03)||0.69 (0.02)|
|20||Bowmans Road (BR-N)||20||5.1 (0.42)||2.6 (0.12)||0.74 (0.06)||0.67 (0.04)|
|23||Jenolan Caves (Jen-N)||30||3.3 (0.33)||2.0 (0.14)||0.56 (0.08)||0.49 (0.06)|
|26||Rocky Plains Creek (RPC-V)||8||3.3 (0.25)||1.6 (0.11)||0.43 (0.09)||0.33 (0.06)|
|27–29||aLittle River (4), Farm Creek (2), Gelantipy Creek (2) (Mer-V)||8||2.1 (0.36)||2.3 (0.15)||0.45 (0.07)||0.54 (0.06)|
| ||All Populations||247||13.4 (1.01)||3.2 (0.07)||0.63 (0.02)||0.59 (0.02)|
We calculated pairwise values of FST for all colonies and tested for significance with FSTAT version 2.9 (Goudet 1995). Population genetic structure was also inferred using a Bayesian model-based clustering analysis in the program STRUCTURE 2.3.1 (Pritchard et al. 2000). STRUCTURE was run under the admixture model with alpha inferred from the data, allele frequencies uncorrelated and lambda set to 1.0. After a burn-in of 100,000 and 200,000, iterations were performed. For the whole data set, we tested the number of genetic clusters (populations, K) present using values of K between 1 and 14, with 10 replicates of each. The inferred number of populations within the sample was deduced using both maximum posterior probability (L(K) Pritchard et al. 2000), and maximum delta log likelihood (ΔK Evanno et al. 2005) implemented in STRUCTURE HARVESTER 0.6.93 (Earl and vonHoldt 2012). Each identified cluster was subsequently rerun to test for additional substructuring within clusters. Finally, we created an unrooted neighbor-joining tree to visualize genetic similarity among the sample sites. The unrooted neighbor-joining tree was based on the average allele sharing genetic distance (Dps) among the 14 populations. Dps values and bootstrap iterations (1000) were calculated in MICROSAT and constructed using the NEIGHBOR and CONSENSUS subroutines in PHYLIP version 3.5 (Felsenstein 1995), with the tree created in TREEVIEW version 1.5.