Molecular Ecology (2010) 19, 2531–2544

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

Patterns of persistence and isolation indicate resilience to climate change in montane rainforest lizards RAYNA C. BELL,* JUAN L. PARRA,† MARIA TONIONE,* CONRAD J. HOSKIN,‡ J A S O N B . M A C K E N Z I E , * S T E P H E N E . W I L L I A M S § and C R A I G M O R I T Z * *Museum of Vertebrate Zoology, University of California, Berkeley, CA 94720, USA, †Department of Ecology and Evolution, State University New York, Stony Brook, NY 11794, USA, ‡Research School of Biology, The Australian National University, Canberra, Australian Capital Territory 0200, Australia, §Centre for Tropical Biodiversity and Climate Change, School of Marine & Tropical Biology, James Cook University, Townsville, Queensland 4810, Australia

Abstract Globally, montane tropical diversity is characterized by extraordinary local endemism that is not readily explained by current environmental variables indicating a strong imprint of history. Montane species often exist as isolated populations under current climatic conditions and may have remained isolated throughout recent climatic cycles, leading to substantial genetic and phenotypic divergence. Alternatively, populations may have become contiguous during colder climates resulting in less divergence. Here we compare responses to historical climate fluctuation in a montane specialist skink, Lampropholis robertsi, and its more broadly distributed congener, L. coggeri, both endemic to rainforests of northeast Australia. To do so, we combine spatial modelling of potential distributions under representative palaeoclimates, multi-locus phylogeography and analyses of phenotypic variation. Spatial modelling of L. robertsi predicts strong isolation among disjunct montane refugia during warm climates, but with potential for localized exchange during the most recent glacial period. In contrast, predicted stable areas are more widespread and connected in L. coggeri. Both species exhibit pronounced phylogeographic structuring for mitochondrial and nuclear genes, attesting to low dispersal and high persistence across multiple isolated regions. This is most prominent in L. robertsi, for which coalescent analyses indicate that most populations persisted in isolation throughout the climate cycles of the Pleistocene. Morphological divergence, principally in body size, is more evident among isolated populations of L. robertsi than L. coggeri. These results highlight the biodiversity value of isolated montane populations and support the general hypothesis that tropical montane regions harbour high levels of narrow-range taxa because of their resilience to past climate change. Keywords: Lampropholis, lineage divergence, montane tropical diversity, Pleistocene refugia Received 22 January 2010; revision received 24 March 2010; accepted 30 March 2010

Introduction Diversity in tropical rainforests is strongly influenced by past geological and climatic effects on distributions of habitats and species that depend on them (Whitmore & Prance 1987; Williams & Pearson 1997; Moritz et al. 2000). Montane tropical forests are characterized by especially high endemism (Blackburn & Measey 2009) that is Correspondence: Rayna C. Bell, Fax: +1 607 255 8088; E-mail: [email protected]  2010 Blackwell Publishing Ltd

not readily explained by current environmental variables such as productivity, topographic complexity and seasonality (Jetz et al. 2004), but rather may reflect responses to past climatic (Graham et al. 2006; Voelker et al. 2010) or geological dynamics (Guarnizo et al. 2009). Tropical species are especially sensitive to climatic fluctuations because their narrow thermal tolerances and elevational ranges can restrict their ability to persist in, or disperse across, alternate habitats (Janzen 1967; Ghalambor et al. 2006; Deutsch et al. 2008; McCain 2009). Climatically stable regions that persisted throughout the Quaternary

2532 R . C . B E L L E T A L . act as refugia for such species, from which populations expand during permissive climates (Hugall et al. 2002; Hewitt 2004; Carnaval et al. 2009). Thus, these regions can represent centres of diversification and accumulation of low dispersal species (Fjeldsa & Lovett 1997). Montane tropical species often exist as multiple isolated populations under the current, relatively warm conditions. It is possible, however, that these populations extended to lower elevations during colder climates (such as existed during glacial periods of the Pleistocene) allowing for more continuous ranges and gene flow (Wiens 2004). Alternatively, populations of tropical montane specialists may have persisted in separate refugia, resulting in high levels of genetic, and perhaps phenotypic, divergence (Garcia-Paris et al. 2000; Bowie et al. 2006; Cadena et al. 2007). Phylogeographic analysis can distinguish between these hypotheses as the former predicts shallow divergence among current geographically isolated populations, whereas the latter will result in pronounced phylogeographic structure with divergence dates predating the last glacial maximum (LGM). Phylogeographic analysis of montane species is most powerful when phylogeography and estimation of demographic parameters extend to multiple, independent loci (Brito & Edwards 2008; Knowles 2009; but see Zink & Barrowclough 2008), and are interpreted in the context of spatially explicit distribution models under representative palaeoclimates (Knowles et al. 2007) with comparative evidence from related species with broader elevation ranges (Bermingham & Moritz 1998). Here we adopt these methods to analyse responses of montane vs. more broadly distributed species of Lampropholis litter skinks, both endemic to the rainforests of northeast Australia—the ‘Australian Wet Tropics’ bioregion. The rainforests of the Australian Wet Tropics have emerged as a model system for understanding effects of climate-driven habitat fluctuation on species and genetic diversity (Williams & Pearson 1997; Schneider et al. 1998; Hugall et al. 2002; Graham et al. 2006). Comparative mitochondrial phylogeographies of mid-montane (> 800 m) vs. elevationally widespread congeners reveal stronger phylogeographic structure in the former (Bell et al. 2007; Moussalli et al. 2009). Palaeodistribution models predict that higher elevation taxa expanded their ranges under colder (cold-dry, glacial; cold-wet, early Holocene) climates, but contracted under the warm-wet phase of the mid-Holocene, before expanding to their current distributions. In contrast, rainforest specialists with broader elevational ranges contracted most strongly under the glacial conditions that dominated much of the last million years (Kershaw et al. 2005; Hocknull et al. 2007).

As yet, there is no phylogeographic analysis of geographically dispersed, truly montane-restricted (> 1000 m) species in the Australian Wet Tropics. Montane habitats are limited and fragmented; thus, extirpation of local populations during climatically restrictive periods may have occurred followed by re-colonization from adjacent, larger refugia during more favorable climates (e.g. Joseph et al. 1995). If populations persist in small montane refugia, then drift or selection-mediated divergence in phenotypes may occur, in contrast to the commonly observed pattern of morphological conservatism among long-isolated lineages in the Wet Tropics (Schneider & Moritz 1999; Hoskin et al. 2005; Hoskin 2007). The high local endemism observed in these montane rainforests is likely due to long-term persistence and divergence between isolated mountaintop populations (Yeates et al. 2002; Hoskin 2004). In this study we compare responses to past climate change of a montane rainforest skink, Lampropholis robertsi (1000–1600 m) with that of its more broadly distributed congener, L. coggeri. Lampropholis robertsi occurs as disjunct populations on the highest mountains of the region, whereas L. coggeri has one of the broadest elevational (50–1300 m) and geographic ranges of any rainforest-restricted species from this region (Fig. 1; Williams et al. 2010). For each species, we first model potential distributions under current and representative past climates. Second, we use multi-locus (mitochondrial and seven to eight nuclear loci) phylogeographic data to test predictions from the palaeodistribution models. Finally, we analyse variation in morphological traits to test for divergence among major intraspecific lineages. To the extent that populations persisted in disjunct refugia, we expect to see marked divergence among populations and this should be most pronounced in the montane specialist.

Methods General sampling We sampled a total of 668 L. coggeri and 78 L. robertsi representing the full geographic and elevational range of each species. The exact location of the point of capture was noted using a global positioning system (GPS) unit and all lizards were measured, sampled, and then released. A small terminal portion of the tail was removed and placed in 95% ethanol for genetic analysis.

Niche modelling We modelled the current geographic distribution of both skinks using contemporary presence records and a  2010 Blackwell Publishing Ltd

P H Y L O G E O G R A P H Y O F T W O R A I N F O R E S T S K I N K S 2533 A

B

C

Fig. 1 (A) Distribution of wet forest in the Australia Wet Tropics. Dashed line represents the buffer used in palaeomodelling analyses. Sampling localities used for modelling analyses relative to the full-parameter modelled distribution of suitable habitat in the current environment for L. robertsi (B) and L. coggeri (C). Arrows point to regions of over-prediction.

maximum entropy algorithm (MAXENT v. 3.2.19, Phillips et al. 2006; Phillips & Dudik 2008). We used a database comprising 335 unique localities for L. coggeri and 39 for L. robertsi at an 80 m2 resolution (Fig. 1; Williams et al. 2010). Contemporary distribution models were developed using seven environmental variables with correlation coefficients of less than 0.6: annual mean temperature (C), mean diurnal range (C), temperature seasonality (unitless), temperature annual range (C), annual precipitation (mm), precipitation seasonality (unitless), and radiation of the driest quarter (MJ ⁄ m2 ⁄ day). These variables were modelled using Anuclim 5.1 (McMahon et al. 1995) coupled with an 80 m resolution digital elevation model. As a background area we included all moist forests and a buffer including areas of sclerophyll forests (VanDerWal et al. 2009; Fig. 1). To evaluate model performance, we set aside a random subset of 25% of the total unique records and measured the area under the curve of the receiver operating characteristic, a threshold-independent measure of performance. We also estimated the potential distribution of skinks at three historic time periods that represent extreme conditions experienced during the late Quaternary: the  2010 Blackwell Publishing Ltd

LGM characterized by cooler and drier climates, the Pleistocene ⁄ Holocene transition (PHT) representing cooler and wetter climates, and the Holocene climate optimum (HCO), which was warmer and wetter (Kershaw & Nix 1988). We used three environmental layers available for all three time periods: annual mean temperature (C), annual precipitation (mm), and precipitation of the driest quarter (mm). These environmental layers were generated following VanderWal et al. (2009) in which layers representing current climatic conditions were adjusted to reflect overall increases or decreases in temperature and precipitation characteristic of the historic time periods (e.g. annual mean temperature differences of )3.5, )2.0, +2.0C for the LGM, PHT and HCO respectively). As in previous studies (Hugall et al. 2002; Moussalli et al. 2009), we developed a current model based on the same contemporary record localities and projected it onto the past climate layers.

Stability and connectivity To evaluate whether measures of historical stability and connectivity predict the distribution of genetic diversity in montane and elevationally-widespread species we

2534 R . C . B E L L E T A L . focused on topographically defined regional populations; 11 in L. coggeri and five in L. robertsi. We estimated habitat stability throughout the late Quaternary by calculating the geometric mean of presence probabilities across the four time periods so that pixels with probabilities of zero suitability at any time period remained as zero in the calculated stability surface. To estimate connectivity, we used the stability surface as a conductance grid (each pixel value representing an index of the amount of movement into the pixel), and measured connectivity among the centroids of each potential refuge using Circuitscape v.3.3 (McRae 2006; McRae & Shah 2009). This measure takes into account all possible paths between two centroids and is superior to the least cost path method (McRae & Beier 2007; McRae et al. 2008). To test the hypothesis that population persistence and connectivity was affected by climate-induced habitat change, we tested for correlation of nucleotide diversity and population size with indicators of habitat stability, including the size of individual refugia (number of pixels with stability > zero) and refuge quality (average non-zero pixel value through time). We then tested whether historical measures of connectivity among stable areas coincide with measures of historical gene flow across populations by performing partial mantel tests of FST and geographic distance and historical connectivity (arithmetic mean of connectivity across the four time periods). We performed 999 iterations to assess significance of correlations. All statistic analyses were conducted in R v. 2.9.2 (http://www.r-project.org).

Genetic data and analysis Sampling details. For mitochondrial (mtDNA) analyses, samples were selected to cover the range of each species. This resulted in 53 L. robertsi samples from six subregions and 240 L. coggeri samples from 26 subregions (Fig. 3B and 4B, Table S1, Supporting Information). For nuclear analyses we selected a subset of the mitochondrial samples that comprised the geographic range of each species as well as the mtDNA diversity within subregions. This resulted in 34 samples of L. robertsi from six subregions and 72 samples of L. coggeri from 26 subregions (Table S1, Supporting Information). Genetic data collection. We extracted total genomic DNA from tail tips preserved in 95% ethanol (Aljanabi & Martinez 1997). We amplified and sequenced one mitochondrial fragment (ND4) in each species as well as eight nuclear loci in L. robertsi and seven in L. coggeri. To amplify ND4, we used published primers (ND4 ⁄ leu; Arevalo et al. 1994) to amplify approximately 950 bases in L. robertsi and 850 bases in L. coggeri. L. robertsi sam-

ples from Mt. Bellenden Kerr (BK) did not amplify with published primers so we designed species-specific primers (Table S2, Supporting Information). PCR reactions included 20 ng of total DNA in a final volume of 10 lL containing: 1X Buffer, 2.5 lM MgCl2, 0.5 lM of each primer, 0.5 lM dNTP mix, and 0.1 units of Taq Polymerase. Amplification was carried out as follows: initial denaturation for 3 min at 94C, followed by 35 cycles (60 s denaturation at 94C, 60 s annealing at 53C, 90 s extension at 72C) and a final extension for 7 min at 72C. For nuclear genes, we used six published loci including a-enolase (Friesen et al. 1997), b-globin (Dolman & Phillips 2004), b-fibrinogen (Prychitko & Moore 1997), Cmos (Saint et al. 1998), Gapd (Friesen et al. 1997) and R35 (Leache 2009). To obtain additional nuclear loci we sequenced clones from a cDNA library for L. coggeri and blasted resulting sequences against Gallus and Xenopus genomes to determine the positions of exon and, by inference, intron regions. Primers, designed using Primer 3 v.0.4.0 (Rozen & Skaletsky 2000), were anchored in conserved regions flanking intron sequence. Details for new loci, which also amplify samples from related genera (Saproscincus and Carlia), are available in Table S2 (Supporting Information). PCR reactions included 20 ng of total DNA in a final volume of 12.5 lL containing: 1.25 lM 10· Buffer, 2 lM Betaine, 0.75 lM of each primer, 0.1 lM dNTP mix, and 0.1 units of Taq Polymerase. Amplification was carried out with initial denaturation for 3 min at 96C, followed by 35 cycles (30 s denaturation at 95C, 30 s annealing at 48– 55C (Table S2, Supporting Information) and 90 s extension at 72C) and a final extension for 7 min at 72C. We amplified a total of approximately 3163 base pairs across eight loci in L. robertsi and 2348 base pairs across seven loci in L. coggeri. PCR products were visualized on an agarose gel, purified using ExoSAP-IT (USB), and sequenced using BigDye v. 3.1 (Applied Biosystems) on an ABI 3730 automated DNA sequencer. DNA sequences were edited using Sequencher 4.7 (Gene Codes Corporation) and are accessioned in GenBank (HM029372–HM030393). Mitochondrial data analysis. Sequences were aligned using ClustalX v. 2.0.10 (Larkin et al. 2007). We implemented MrModelTest v. 2.3 (Nylander 2004) to establish that the HKY+G and TIM+1 + G models best represented the substitution processes for the L. robertsi and L. coggeri datasets respectively. Bayesian phylogenetic analyses were conducted using MRBAYES v. 3.1.2 (Huelsenbeck & Ronquist 2001) with datasets partitioned by codon position. We allowed four incrementally heated Markov chains to proceed for 20 million generations, sampling every 1000 generations. Bayesian posterior  2010 Blackwell Publishing Ltd

P H Y L O G E O G R A P H Y O F T W O R A I N F O R E S T S K I N K S 2535 probability values were estimated from the sampled trees remaining after 20000 burn-in samples were discarded. We used Arlequin v. 3.1 (Excoffier et al. 2005) to calculate Fu’s Fs (using the full mitochondrial dataset) to test for population expansion within distinct mitochondrial lineages. Using the subset of mitochondrial samples also sequenced for nuclear loci, we calculated FST, nucleotide diversity (hs, hp), sequence divergence among populations (Dxy and Da using the Tamura-Nei model) and hierarchical F statistics for the northern and southern lineages, as well the topographically defined subregions. Nuclear data analysis. To resolve haplotypes from heterozygous individuals, we used PHASE v. 2.1 (Stephens et al. 2001). We used multi-locus genotypes, with phased haplotypes coded as allelic states, to infer population structure and assign individuals to genetic populations with the Bayesian assignment program STRUCTURE v. 2.2.3 (Pritchard et al. 2000). Based on preliminary analyses, we ran ten iterations of each value of K from one through eight in L. coggeri and K from one to six in L. robertsi. We used a burn-in of 1 000 000 steps, MCMC length of 3 000 000 steps, correlated allele frequencies and the admixture ancestry model. To determine the optimal number of genetic clusters in our dataset, we used the method described by Evanno et al. (2005). STRUCTURE outputs were plotted with DISTRUCT v1.1 (Rosenberg 2004). For the northern and southern lineages, as well as the geographic subregions, we calculated FST, nucleotide diversity (hs, hp), hierarchical F statistics, and average and net sequence divergence among populations (Dxy and Da using the Tamura–Nei model). All calculations were conducted in Arlequin v. 3.1 (Excoffier et al. 2005). To visually represent overall divergence patterns, we clustered geographic subregions according to pairwise Da values using the neighbor-joining algorithm in PAUP v. 4 (Swofford 2003). We also produced individual gene trees for all nuclear loci using the neighbor-joining method in PAUP v. 4 (Figs S1 and S2, Supporting Information). We used LAMARC v. 2.1.3 (Kuhner 2006) to estimate the growth rates of the genetic populations identified in the STRUCTURE analyses. For each lineage we used the phased nuclear (nuDNA) datasets to estimate the growth parameter (g). We used the Bayesian method with 100 000 final recorded parameter sets, a burn-in of 40 000 steps, sampling every 100, and adjusted the growth rate range from –5000 to 15 000. The F84 model of nucleotide substitution was used for all loci. Combined data analysis. We used IMa (Hey & Nielsen 2007) to estimate divergence time, ancestral population  2010 Blackwell Publishing Ltd

and current population sizes. We found no evidence of recombination using the difference of sums of squares method (TOPALi v. 2.5, Milne et al. 2004) and implemented IMa in M mode with two independent runs for each analysis. Each run consisted of seven geometrically heated chains with 20 000 000 steps and a burn-in of 100 000 such that effective sample size estimates were greater than 200 and autocorrelations of parameter values were close to zero. In L. robertsi, we performed the analyses in pairwise comparisons of the four major evolutionary lineages. Estimates of population size and divergence time were consistent whether migration was included or excluded. For L. coggeri, analyses were performed between the two main lineages. We do not report results of this analysis because estimates of the divergence time parameter did not converge. Compound parameter estimates for L. robertsi were transformed assuming a mtDNA (ND4) divergence rate of 0.65% per Myr per lineage (Macey et al. 1998), mutation scalars for nuclear genes as obtained from posterior distributions in IMa, and a generation time of 1 year.

Phenotypic divergence We collected the following measurements from live animals in the field: snout to vent length (SV); head width (HW); head length (HL, anterior edge of the tympanum to tip of snout); and hind limb length (HLL, length of tibia with leg bent). Analyses included only mature individuals, defined as the minimum size of confirmed males: SV > 39 mm for L. robertsi and SV > 31.5 mm for L. coggeri. We tested for morphological divergence between major lineages revealed by genetic analyses. Sexual dimorphism in body size and shape was evident for both species, so each sex was analysed separately. The methods for each analysis follow that outlined in Hoskin et al. (2005). We conducted multivariate (PCA and MANOVA) analyses of size and shape divergence among lineages, removing effects of elevation as appropriate (detailed methods in Data S1).

Results Palaeomodelling Predicted distributions in current environment. The full parameter Maxent model of L. robertsi correctly predicted observed presence (AUCtest = 0.999 ± 0.001), but also over-predicted in some areas. In particular, the isolated high elevation areas of the Lamb Uplands (LU) and Windsor Uplands (WU) were predicted to be suitable though repeated surveys detected L. coggeri and not L. robertsi in these regions (Fig. 1B). The model predicted a current distribution of L. coggeri that closely

2536 R . C . B E L L E T A L . matches the observed distribution (AUCtest = 0.891 ± 0.016) with slight over-prediction to the west of WU (Fig. 1C). Total annual precipitation and annual mean temperature contributed similarly to the L. coggeri model (38% and 24% respectively) while the L. robertsi model was primarily driven by annual mean temperature (89%; Table S3, Supporting Information). The performances of the three variable models of L. robertsi and L. coggeri for the current climate were similar to the full models (AUCtest = 0.997 ± 0.001, AUCtest = 0.882 ± 0.016, respectively; Table S3, Supporting Information), again with precipitation variables dominant for L. coggeri and temperature for L. robertsi. Predicted distributions under past environments and regions of habitat stability. The predicted distribution of L. robertsi was severely restricted during the warm-wet period of the mid-Holocene for which only the highest peaks were predicted to be suitable (Fig. 2A). The cold climates of the LGM (cold-dry) and early Holocene (cold-wet) were more permissive with predicted distributions that extended beyond the observed current distribution, especially across the Atherton Uplands (AU). During the most permissive climate (LGM) the model predicted potential connectivity between Mt. Bartle Frere (BF) and the AU, which now are isolated. In contrast, the model predicted that other current populations Mt. Bellenden Ker (BK), the Thornton Uplands (TU), and the Carbine Uplands (CU), remained isolated throughout. Only the highest peaks in the central (western AU, BK, BF) and the northern (TU, CU, WU) Wet Tropics showed strong predictions of habitat stability (Fig. 2A). In contrast to L. robertsi, the predicted distribution of L. coggeri was most restricted during the LGM, during which the suitable range was most extensive in the central Wet Tropics with scattered predicted distributions to the north and south (Fig. 2B). The model predicted high connectivity across the latitudinal range during the early Holocene (cold-wet) period. Connectivity was predicted to decrease during the mid-Holocene (warmwet) period, with slight expansion to the present. Predicted refugia were more extensive and continuous than for L. robertsi. Relatively large refugia were predicted in Finnegan Uplands (FU), TU, CU, WU and the Atherton Tablelands (Fig. 2B). Smaller stable areas were predicted in LU in the north and Hinchinbrook Island (HI), Kirrama Uplands (KU), and Spec Uplands (SU) in the south.

Phylogeography Mitochondrial and nuclear phylogeography. The mitochondrial phylogeny of L. robertsi (Fig. 3A) recovered

exceptionally strong geographic structuring of diversity, with a primary split between northern and southern lineages (Da = 8.8%; Table 1) and further deep divisions within the northern and the southern populations (Fig. 3B). The BF population clustered closely with the western Atherton populations, rather than the geographically proximal BK population. For the nuclear genes, the STRUCTURE analysis consistently recovered four demes in L. robertsi (TU, CU, BK and BF + AU) with no evidence of admixture (Fig. 3C). The population phylogram, a visual representation of pairwise divergence (Da) between geographic regions, clustered the northern populations (TU and CU), and BF with AU populations; the BK lineage was the most distinct (Fig. 3D). Remarkably, populations of L. robertsi frequently had monophyletic alleles at individual nuclear loci (Fig. S1, Supporting Information). Pairwise estimates of population divergence (Da) based on nuclear and mitochondrial analyses were correlated (r = 0.88, P < 0.05), though estimates of molecular diversity within populations (hp) were not (r = –0.02, P > 0.05). Local subdivision in L. robertsi was pronounced with an overall FST of 0.97 (P < 0.001) for mtDNA, 0.89 (p < 0.001) for nuDNA, and mean net divergence between populations of 11.7% for mtDNA and 1.3% for nuDNA (Tables 1, S4A&B, Supporting Information). Correspondingly, hierarchical F statistics were high and significant, except among grouped northern and southern populations (Table 1). The mitochondrial phylogeny of L. coggeri (Fig. 4A) recovered two major lineages with substantial genetic divergence (Da = 9.4%; Table 1) that correspond to the northern and southern Wet Tropics, and are currently parapatric at the northern AU (Fig. 4B; AUCE). Both lineages were highly structured with four sub-lineages in the north and five in the south. Two peripheral populations, Bakers Blue (BB) west of the CU rainforest buffer and Hervey Range (TV) at southern limit of the range, had phylogenetically distinct mitochondrial lineages. The STRUCTURE analysis recovered six demes that correspond to three geographic regions in the north and three in the south (Fig. 4C). These demes were largely concordant with the mitochondrial lineages though some peripheral populations clustered with the adjacent geographic regions. The two northernmost demes exhibited limited admixture (the WU and TU populations) while none was observed between the three southern demes (Fig. 4C). Compared to L. robertsi, individual nuclear gene trees for L. coggeri appeared less geographically structured (Fig. S2, Supporting Information). Nonetheless, the L. coggeri phylogram of pairwise nuDNA divergences (Fig. 4D and Table S5C, Supporting Information) was generally consistent with the STRUCTURE results and the  2010 Blackwell Publishing Ltd

P H Y L O G E O G R A P H Y O F T W O R A I N F O R E S T S K I N K S 2537 A

B

Fig. 2 (A) (L. robertsi) and (B) (L. coggeri) three parameter models of potential distributions under current and past climates. Shown are the LGM (cold-dry) the PHT (cold-wet), the HCO (warm-wet), and current environments. The stability surface was obtained by taking the geometric mean of the suitable pixels through all four time periods.

mtDNA phylogeny. Heterogeneity was more pronounced in the south where the three population groupings in the phylogram corresponded precisely to the STRUCTURE demes. However, several peripheral localities (i.e. MT, BB and TV) were more distinct in the nuDNA phylogram (which considers sequence divergence among haplotypes) than the STRUCTURE analysis (which does not). As in L. robertsi, pairwise estimates of population divergence (Da) in L. coggeri, based on nuclear and mitochondrial analyses, were highly correlated (r = 0.83, P < 0.001), whereas estimates of molecular diversity within populations (hp) were not (r = –0.43, P > 0.05). For both mtDNA and nuDNA, population subdivision and regional structure were substantial (Table 1) with an average pairwise FST between geographic populations of 0.70 for mtDNA and 0.80 for nuDNA (Table S5D, Supporting Information). In contrast to L. robertsi, the proportion of molecular diversity among populations within northern and southern lineages (FSC) was less than that between these areas (FCT; Table 1). Diversity & historical demography. Overall, stable refuge size (number of non-zero pixels) was much smaller in L. robertsi than in L. coggeri, though refuge quality (geo 2010 Blackwell Publishing Ltd

metric mean of pixel value through time) was slightly greater (Fig. 2). There was no relationship between refuge quality and corresponding nuclear diversity (hp) in either L. robertsi (r = –0.08, P > 0.10) or L. coggeri (r = 0.21, P > 0.10). Refuge size and nuclear diversity were marginally correlated in L. robertsi (adjusted r2 = 0.82, P < 0.10) but not in L. coggeri (adjusted r2 = 0.36, P > 0.10). As predicted, both species exhibited a negative correlation between pairwise population divergence (nuDNA FST) and geographic connectivity though the relationship was stronger in L. robertsi (adjusted r2 = 0.819, P < 0.10 than L. coggeri (adjusted r2 = 0.365, P < 0.05). We detected population growth in several populations of L. robertsi and L. coggeri based on analyses of both mtDNA (Fu’s Fs) and nuDNA (LAMARC). In L. robertsi, both northern demes (TU and CU) and one southern deme (BK) showed signatures of expansion with estimates of g = 14,083 (4177–14389), 1375 (438–4653), and 2595 (510–13197), respectively. Consistent with these results, the TU and CU demes had significantly negative values of Fu’s Fs for mtDNA (TU, Fs = –13.46, P < 0.05; CU, Fs = –6.85, P < 0.05). The Atherton tableland deme (AU ⁄ BF) also had significant values of Fu’s

2538 R . C . B E L L E T A L . A

B

C

D

Fig. 3 (A) L. robertsi mitochondrial tree including all samples sequenced. Posterior probabilities greater than 95% are denoted by asterisks. (B) Distribution of sampling localities in L. robertsi relative to the distribution of current suitable habitat. (C) Individual assignment probabilities of L. robertsi samples for K = 4 in STRUCTURE analysis. (D) Population phylogram based on net divergence across nuclear loci for L. robertsi. In all cases, sampling localities are colored according to the demes inferred from STRUCTURE analysis (C) of nuclear loci. Table 1 Summary of nucleotide divergence estimates (Dxy, Da, FST) and analysis of molecular variance (AMOVA) for the mitochondrial and nuclear datasets of L. coggeri and L. robertsi. Results are shown at the level of northern vs. southern lineages as well as geographic groups, which refer to the topographic subregions (five in L. robertsi and eleven in L. coggeri) found within the primary northern and southern lineages North vs. South

L. robertsi L. coggeri

mtDNA nuDNA mtDNA nuDNA

Geographic groups

Dxy (%)

Da (%)

FST

Mean Dxy (%)

Mean Da (%)

Mean FST

FSC

FCT

FST

14.2 1.7 13.0 1.7

8.8 1.2 9.4 1.1

0.62*** 0.66*** 0.71* 0.67*

11.9 1.4 9.2 1.3

11.7 1.3 7.4 1.2

0.96* 0.89* 0.70 0.80

0.97*** 0.94*** 0.48*** 0.78***

0.32 0.40 0.68** 0.61***

0.98*** 0.96*** 0.84*** 0.91***

Significance of * P < 0.05, ** P < 0.01 and *** P < 0.001.

Fs (–20.95, P < 0.05) whereas the BK population did not. For L. coggeri, analyses of the nuDNA loci revealed significant population growth in the northern-most deme (FU ⁄ TU; g = 563, 231–1300) and in the central deme (AU ⁄ CC ⁄ MT; g = 969, 525–2024). Significantly negative values of Fs were obtained for the northern CU ⁄ WU deme (Fs = –3.53, p < 0.05), while values for the northern FU ⁄ TU and central Atherton demes were negative, but not significant (Fs = –1.64 and Fs = -1.82 respectively). LAMARC estimates of ‘g’ for the remaining demes were either zero or the confidence intervals included negative values. Coalescent (IMa) estimates of pairwise divergence time between the four L. robertsi demes indicate that these populations persisted in isolation since the Plio-

cene (2.3–4.6 Myr; Table 2), with minimum divergence estimates mostly predating the Pleistocene (lower 95% limits > 1.6 Myr). Estimates of within-lineage population size are generally consistent across pairwise estimates and indicate that the populations at CU and AU are approximately equivalent, and are larger than those at TU and BK. Estimates for L. coggeri are not reported because parameter estimates did not converge.

Morphological divergence among lineages The major lineages in both species exhibited morphological divergence, though the analysis for L. robertsi should be regarded as preliminary due to small sample sizes and a lack of data for males of the TU lineage.  2010 Blackwell Publishing Ltd

P H Y L O G E O G R A P H Y O F T W O R A I N F O R E S T S K I N K S 2539 A

B

C

D

Fig. 4 (A) L. coggeri mitochondrial tree including only unique sampling localities or unique haplotypes. Posterior probabilities greater than 95% are denoted by asterisks. (B) Distribution of sampling localities in L. coggeri relative to the distribution of current suitable habitat. (C) Individual assignment probabilities of L. coggeri samples for K = 6 in STRUCTURE analysis. (D) Population phylogram based on net divergence across nuclear loci for L. coggeri. In all cases, sampling localities are colored according to the demes inferred from STRUCTURE analysis (C) of nuclear loci.

Table 2 Converted parameter estimates from L. robertsi IMa analyses including current population sizes (q1, q2), ancestral population sizes (qA), all shown in thousands, and divergence time (T), shown in millions of years. Conversions apply the estimated mutational scalar for ND4 and assume a mtDNA substitution rate of 0.65% per lineage per million years. Values shown are peak probabilities with 95% confidence intervals q1

q2

qA

T

Populations

High Pt

95% L

95% H

High Pt

95% L

95% H

High Pt

95% L

95% H

High Pt

95% L

95% H

TU(1)–AU(2) CU(1)–AU(2) CU(1)–TU(2) AU(1)–BK(2) CU(1)–BK(2) TU(1)–BK(2)

383 762 1368 631 830 243

128 544 870 348 362 243

638 1423 2878 1050 2021 553

894 762 371 213 117 243

638 544 371 71 117 78

1404 1198 870 348 362 398

6510 4480 627 2029 117 1669

2680 980 128 213 117 243

21307 14739 5116 8044 6063 10720

5.74 5.13 5.37 3.55 7.02 5.24

3.19 2.96 2.81 1.56 2.77 1.75

11.27 10.73 9.46 7.52 11.70 10.67

In both species, morphological differentiation occurred primarily in body size (PC1) with little detectable divergence in body shape (see Supplemental Materials for details). For L. robertsi, size divergence was evident among lineages but the pattern of results differed between sexes. For females, the Atherton Tableland lineage was significantly smaller, whereas for males the CU lineage was significantly larger (Fig. 5A). Additionally, some shape differences were detected among L. robertsi lineages for females, but not for males. For L. coggeri, females from the southern lineage were significantly larger than those of the northern lineage, but  2010 Blackwell Publishing Ltd

size of males did not differ (Fig. 5B). The two lineages did not differ in body shape for either sex.

Discussion This study extends previous phylogeographic analyses by comparing spatio-temporal patterns of persistence and isolation between montane and elevationally-widespread species that are otherwise ecologically similar. By combining assessments of multi-locus diversity with spatial modelling of current and palaeodistributions we aim to shed light on the historical processes that con-

2540 R . C . B E L L E T A L . A

Fig. 5 Boxplots of PC1 (body size) divergence between major genetic lineages in (A) L. robertsi females and males, and (B) L. coggeri females and males. In all four analyses PC1 accounts for the majority of variation in the morphological data and represents body size. The box plots show the median, 25th and 75th quartiles, and minimum and maximum data. Populations that differ significantly in size are denoted by asterisks.

B

tribute to high local endemism of montane rainforests. The results reveal strong geographic structuring of diversity in both species, consistent with persistence in small refugia and low dispersal overall, but also extreme, genealogically structured divergence among most populations of the montane specialist. For these species, mtDNA alone identifies major historical lineages (Zink & Barrowclough 2008), and these are corroborated by the analysis of multiple nuclear loci. Yet, the discordance between mitochondrial and nuclear estimates of diversity within lineages confirms that the latter are critical to obtain reliable estimates of historical demographic parameters (population size and divergence time) that provide insight into underlying processes (Brito & Edwards 2008). Palaeomodels predict that the higher elevation species, L. robertsi, was most restricted during warmer climates. This result parallels comparative modelling for another set of rainforest skinks (Saproscincus; Moussalli et al. 2009) and beetles (Bell et al. 2007). The palaeomodels indicate that L. robertsi populations were severely isolated during warmer climates, yet were potentially more widespread under colder climates, including glacial conditions, as occurred for much of the past million years. In contrast, like most taxa from the Australian Wet Tropics (Hugall et al. 2002; Bell

et al. 2007; Moussalli et al. 2009) and the rainforest itself (Graham et al. 2006; VanDerWal et al. 2009), models of the more broadly distributed congener L. coggeri predict restriction during the cold-dry LGM with increasing suitability through the Holocene to current climate. Despite pronounced current fragmentation, several analyses indicate that there was genetic exchange between populations during permissive climates in both species. In L. robertsi, spatial models predict that the currently disjunct BF population was contiguous with populations on the adjacent AU during colder climates. This prediction is supported by the close genetic similarity between these populations relative to the continuously isolated northern (CU, TU) and BK populations. Similarly, L. coggeri is predicted to have been more contiguous in the northern regions during the cold-wet early Holocene than at present, which is consistent with the overall genetic similarity among these areas and indications of admixture from the analysis of nuDNA and the mtDNA gene tree. Overall, the spatial models adequately predict current distributions and predictions under past climates are broadly consistent with inferred demographic responses. Molecular diversity within isolated L. robertsi populations is low and correlated with predicted refuge size. Smaller refugia in the north (TU) and south (BK)  2010 Blackwell Publishing Ltd

P H Y L O G E O G R A P H Y O F T W O R A I N F O R E S T S K I N K S 2541 likely experienced population reductions during the mid Holocene and these populations, along with CU, show the clearest signatures of population expansion. In contrast, predicted L. coggeri refugia are much larger and less isolated, thus the effects of habitat contraction on molecular diversity are less substantial, though detectable for several populations. Pairwise population divergence (FST) and predicted connectivity between populations are negatively correlated in both L. coggeri and L. robertsi, which confirms that modelled low suitability areas limit migration between populations. The relationship is maintained across all scales of geographic connectivity and indicates that this method of estimating connectivity between historic refugia can serve as a proxy for expected molecular divergence in moderately to severely isolated populations. The most obvious discrepancy in the spatial models is that they strongly predict L. robertsi on WU throughout the current and historic climatic conditions, yet repeated surveys do not detect the species in this region. Though L. robertsi are not present on WU, L. coggeri are commonly found and their elevational range here, in the absence of L. robertsi, is higher than throughout the rest of their distribution (to 1300 m). Where the two species are locally sympatric, along the western AU, the upper limit of L. coggeri is approximately 1000 m and where they co-occur, the species partition the microhabitat such that L. coggeri is found along warmer ridgelines and L. robertsi is found in cooler, shady gullies. This pattern raises the possibility of competitive interaction between the species and competitive release as an explanation of the increased elevational range of L. coggeri in WU (cf. Anderson et al. 2002). Alternatively, microhabitats above 1000 m on WU may be more favourable for L. coggeri than L. robertsi such that the species’ distributions simply result from divergent microhabitat (e.g. thermal or moisture) preferences. For both species, the presence of strong phylogeographic structure across both mtDNA and multiple nuclear loci attests to high persistence in historically isolated rainforest patches. Genetic diversity in Lampropholis is partitioned at a finer scale and with deeper divergences than most other vertebrates studied in this system (reviewed in Moritz et al. 2009), and is on par with that in snails (Hugall et al. 2002). In particular, phylogeographic structure in L. coggeri is much stronger than in Carlia rubrigularis, whose geographic range overlaps extensively with L. coggeri (Dolman & Moritz 2006). Lampropholis coggeri persists in small rainforest patches in the south (TV) and west (BB) where C. rubrigularis (and most other rainforest herpetofauna) are absent, which reflects higher persistence of L. coggeri in small, peripheral habitat fragments (Kanowski et al.  2010 Blackwell Publishing Ltd

2006). Despite this ability to persist in eco-geographically marginal rainforest areas, L. coggeri, like most rainforest species, was strongly affected by long-term separation of rainforests in the northern and central wet tropics—the ‘Black Mountain Barrier’ (Schneider et al. 1998; Joseph et al. 2005). Most remarkable is the extreme molecular divergence among current and historical populations of the montane species, L. robertsi. Coalescent analyses place divergence times of the four major L. robertsi lineages in the Pliocene, with minimum estimates predating the early Pleistocene (> 1.6 Myr). Thus, these montane rainforest skinks persisted as multiple isolated populations throughout the high amplitude and long-duration climate cycles of the Pleistocene. Though we are unable to obtain direct estimates of divergence time between the northern and southern lineages of L. coggeri, the magnitude of genetic divergence is on par with that observed among L. robertsi lineages, suggesting that divergence between northern and southern lineages also predates the Pleistocene. Though divergence times predate the available palaeoclimatic surfaces (LGM to the present), models of historic climate cycles provide spatially explicit predictions of exchange between, and persistence of, geographically isolated populations. Here, the palaeomodels predict overall patterns of persistence and isolation through the late Pleistocene as well as also over much longer time periods (to the Pliocene), indicating that previous climatic oscillations resulted in spatially congruent processes. This consistency is expected given the overall geological antiquity and stability of the system and implies that late Pleistocene models provide a reasonable surrogate for spatial dynamics over much longer timescales. Long-term population isolation, as observed between the northern and southern lineages of L. coggeri and most populations of L. robertsi, allows for morphological divergence due to divergent selection and ⁄ or genetic drift. Divergence in body size between northern and southern females of L. coggeri, though modest, is noteworthy as most studies of broadly distributed Wet Tropics species detect little morphological variation among genetically divergent lineages (Schneider & Moritz 1999; Schneider et al. 1999; Hoskin et al. 2005; Hoskin 2007). Though limited by small sample sizes, morphological differentiation between the four genealogically distinct populations of L. robertsi appears more pronounced. This preliminary result indicates that morphological divergence may be most evident in populations that persist in small refugia and remain isolated for extended periods of time. The evidence from multi-locus phylogeography and spatial modelling of potential palaeodistributions demonstrates persistence and isolation of most populations

2542 R . C . B E L L E T A L . of the montane species, L. robertsi, throughout the strong climate oscillations of the late Pleistocene, and likely extending back to the Pliocene. Other than transient connections between the adjacent BF and AU refugia, the results contrast with the expectation that montane populations would be reconnected during colder climates (Wiens 2004). These isolated refugia, especially BK, TU and CU, are particularly rich in narrowly endemic species, which is attributed to their relative stability during recent episodes of climate change (Williams & Pearson 1997; Yeates et al. 2002; Graham et al. 2006; VanDerWal et al. 2009). The presence of genetically divergent populations of the more broadly-distributed L. coggeri in geographically marginal rainforest patches, as previously observed for some other species in this system (Moussalli et al. 2009), also points to the importance of such areas as nuclei of unique diversity. Together, the above observations support the general hypothesis that isolated tropical montane regions harbour high levels of narrow-range taxa because of their resilience to past climate change (Fjeldsa & Lovett 1997; Jetz et al. 2004). Montane tropical areas are thought to be highly susceptible to effects of future humaninduced global warming (Williams et al. 2007; Wright et al. 2009), and the same is true for the Australian Wet Tropics (Williams et al. 2003; Williams & Hilbert 2006). At first sight, species such as L. robertsi would seem especially prone to local extinction and loss of considerable genetic diversity with any further warming; yet, these populations and those of other high-montane endemic species (Cophixalus frogs; Hoskin 2004) have evidently persisted through past warming events. Understanding what features of these organisms, or their habitats (Shoo et al. 2010), enabled persistence will improve predictions of future climate change impacts in this and other highly endemic tropical systems (Williams et al. 2008).

Acknowledgments We thank the National Science Foundation (DEB 0416250 and 0817035) and Queensland Smart State scheme for support, the Moritz lab, Zamudio lab, and two anonymous reviewers for valuable comments on previous versions of the manuscript, and Sean Rovito for advice on spatial modelling. STRUCTURE, IMa, and MRBAYES were performed at Cornell Computational Biology Service Unit, a facility partially funded by Microsoft Corporation.

References Aljanabi SM, Martinez I (1997) Universal and rapid saltextraction of high quality genomic DNA for PCR-based techniques. Nucleic Acids Research, 25, 4692–4693.

Anderson RP, Peterson AT, Gomez-Laverde M (2002) Using niche-based GIS modeling to test geographic predictions of competitive exclusion and competitive release in South American pocket mice. Oikos, 98, 3–16. Arevalo E, Davis SK, Sites JW (1994) Mitochondrial DNA sequence divergence and phylogenetic relationships among eight chromosome races of the Sceloporus grammicus complex (Phrynosomatidae) in Central Mexico. Systematic Biology, 43, 387–418. Bell KL, Moritz C, Moussalli A, Yeates DK (2007) Comparative phylogeography and speciation of dung beetles from the Australian Wet Tropics rainforest. Molecular Ecology, 16, 4984–4998. Bermingham E, Moritz C (1998) Comparative phylogeography: concepts and applications. Molecular Ecology, 7, 367–369. Blackburn DC, Measey GJ (2009) Dispersal to or from an African biodiversity hotspot? Molecular Ecology, 18, 1904–1915. Bowie RCK, Fjeldsa J, Hackett SJ, Bates JM, Crowe TM (2006) Coalescent models reveal the relative roles of ancestral polymorphism, vicariance, and dispersal in shaping phylogeographical structure of an African montane forest robin. Molecular Phylogenetics and Evolution, 38, 171–188. Brito PH, Edwards SV (2008) Multilocus phylogeography and phylogenetics using sequence-based markers. Genetica, 135, 439–455. Cadena CD, Klicka J, Ricklefs RE (2007) Evolutionary differentiation in the Neotropical montane region: molecular phylogenetics and phylogeography of Buarremon brushfinches (Aves, Emberizidae). Molecular Phylogenetics and Evolution, 44, 993–1016. Carnaval AC, Hickerson MJ, Haddad CFB, Rodrigues MT, Moritz C (2009) Stability predicts genetic diversity in the Brazilian Atlantic Forest Hotspot. Science, 323, 785–789. Deutsch CA, Tewksbury JJ, Huey RB et al. (2008) Impacts of climate warming on terrestrial ectotherms across latitude. Proceedings of the National Academy of Sciences, USA, 105, 6668–6672. Dolman G, Moritz C (2006) A multilocus perspective on refugial isolation and divergence in rainforest skinks (Carlia). Evolution, 60, 573–582. Dolman G, Phillips B (2004) Single copy nuclear DNA markers characterized for comparative phylogeography in Australian wet tropics rainforest skinks. Molecular Ecology Notes, 4, 185–187. Evanno G, Regnaut S, Goudet J (2005) Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Molecular Ecology, 14, 2611–2620. Excoffier L, Laval G, Schneider S (2005) Arlequin (version 3.0): an integrated software package for population genetics data analysis. Evolutionary Bioinformatics, 1, 47–50. Fjeldsa J, Lovett JC (1997) Biodiversity and environmental stability. Biodiversity and Conservation, 6, 315–323. Friesen VL, Congdon BC, Walsh HE, Birt TP (1997) Intron variation in marbled murrelets detected using analyses of single-stranded conformational polymorphisms. Molecular Ecology, 6, 1047–1058. Garcia-Paris M, Good DA, Parra-Olea G, Wake DB (2000) Biodiversity of Costa Rican salamanders: implications of high levels of genetic differentiation and phylogeographic structure for species formation. Proceedings of the National Academy of Sciences, USA, 97, 1640–1647.

 2010 Blackwell Publishing Ltd

P H Y L O G E O G R A P H Y O F T W O R A I N F O R E S T S K I N K S 2543 Ghalambor CK, Huey RB, Martin PR, Tewksbury JJ, Wang G (2006) Are mountain passes higher in the tropics? Janzen’s hypothesis revisited Integrative and Comparative Biology, 46, 5– 17. Graham CH, Moritz C, Williams SE (2006) Habitat history improves prediction of biodiversity in rainforest fauna. Proceedings of the National Academy of Sciences, USA, 103, 632– 636. Guarnizo CE, Amezquita A, Bermingham E (2009) The relative roles of vicariance versus elevational gradients in the genetic differentiation of the high Andean tree frog, Dendropsophus labialis. Molecular Phylogenetics and Evolution, 50, 84–92. Hewitt GM (2004) Genetic consequences of climatic oscillations in the Quaternary. Philosophical Transactions of the Royal Society of London Series B-Biological Sciences, 359, 183–195. Hey J, Nielsen R (2007) Integration within the Felsenstein equation for improved Markov chain Monte Carlo methods in population genetics. Proceedings of the National Academy of Sciences, USA, 104, 2785–2790. Hocknull SA, Zhao JX, Feng YX, Webb GE (2007) Responses of Quaternary rainforest vertebrates to climate change in Australia. Earth and Planetary Science Letters, 264, 317–331. Hoskin CJ (2004) Australian microhylid frogs (Cophixalus and Austrochaperina): phylogeny, taxonomy, calls, distributions and breeding biology. Australian Journal of Zoology, 52, 237–269. Hoskin CJ (2007) Description, biology and conservation of a new species of Australian tree frog (Amphibia: Anura: Hylidae: Litoria) and an assessment of the remaining populations of Litoria genimaculata Horst, 1883: systematic and conservation implications of an unusual speciation event. Biological Journal of the Linnean Society, 91, 549–563. Hoskin CJ, Higgie M, McDonald KR, Moritz C (2005) Reinforcement drives rapid allopatric speciation. Nature, 437, 1353–1356. Huelsenbeck JP, Ronquist F (2001) MRBAYES: Bayesian inference of phylogenetic trees. Bioinformatics, 17, 754–755. Hugall A, Moritz C, Moussalli A, Stanisic J (2002) Reconciling paleodistribution models and comparative phylogeography in the Wet Tropics rainforest land snail Gnarosophia bellendenkerensis (Brazier 1875). Proceedings of the National Academy of Sciences, USA, 99, 6112–6117. Janzen DH (1967) Why mountain passes are higher in the tropics. American Naturalist, 101, 233–249. Jetz W, Rahbek C, Colwell RK (2004) The coincidence of rarity and richness and the potential signature of history in centres of endemism. Ecology Letters, 7, 1180–1191. Joseph L, Moritz C, Hugall A (1995) Molecular support for vicariance as a source of diversity in rainforest. Proceedings of the Royal Society of London Series B-Biological Sciences, 260, 177–182. Kanowski JJ, Reis TM, Catterall CP, Piper SD (2006) Factors affecting the use of reforested sites by reptiles in cleared rainforest landscapes in tropical and subtropical Australia. Restoration Ecology, 14, 67–76. Kershaw AP, Nix HA (1988) Quantitative palaeoclimatic estimates from pollen data using bioclimatic profiles of extant taxa. Journal of Biogeography, 15, 589–602. Kershaw AP, Moss PT, Wild R (2005) Patterns and causes of vegetation change in the Australian Wet Tropics region over the last 10 million years. In: Tropical Rainforests: Past, Present,

 2010 Blackwell Publishing Ltd

and Future (eds. Bermingham E, Dick CW, Moritz C), pp. 374–400. University of Chicago Press, Chicago, IL. Knowles LL (2009) Statistical phylogeography. Annual Review of Ecology, Evolution and Systematics, 40, 593–612. Knowles LL, Carstens BC, Keat ML (2007) Coupling genetic and ecological-niche models to examine how past population distributions contribute to divergence. Current Biology, 17, 940–946. Kuhner MK (2006) LAMARC 2.0: maximum likelihood and Bayesian estimation of population parameters. Bioinformatics, 22, 768–770. Larkin MA, Blackshields G, Brown NP et al. (2007) Clustal W and clustal X version 2.0. Bioinformatics, 23, 2947–2948. Leache AD (2009) Species tree discordance traces to phylogeographic clade boundaries in North American Fence Lizards (Sceloporus). Systematic Biology, 58, 547–559. Macey JR, Schulte JA, Ananjeva NB et al. (1998) Phylogenetic relationships among agamid lizards of the Laudakia caucasia species group: testing hypotheses of biogeographic fragmentation and an area cladogram for the Iranian Plateau. Molecular Phylogenetics and Evolution, 10, 118–131. McCain CM (2009) Vertebrate range sizes indicate that mountains may be ‘higher’ in the tropics. Ecology Letters, 10, 55–560. McMahon JP, Hutchinson MF, Nix HA, Ord KD (1995) ANUCLIM User’s Guide v. 1.0. Centre for Resource and Environmental Studies, Australian National University, Canberra. McRae BH (2006) Isolation by resistance. Evolution, 60, 1551– 1561. McRae BH, Beier P (2007) Circuit theory predicts gene flow in plant and animal populations. Proceedings of the National Academy of Sciences, USA, 104, 19885–19890. McRae BH, Shah VB (2009) Circuitscape User’s Guide (http:// www.circuitscape.org). The University of California, Santa Barbara, CA. McRae BH, Dickson BG, Keitt TH, Shah VB (2008) Using circuit theory to model connectivity in ecology, evolution, and conservation. Ecology, 89, 2712–2724. Milne I, Wright F, Rowe G et al. (2004) TOPALi: software for automatic identification of recombinant sequences within DNA multiple alignments. Bioinformatics, 20, 1806–1807. Moritz C, Patton JL, Schneider CJ, Smith TB (2000) Diversification of rainforest faunas: an integrated molecular approach. Annual Review of Ecology and Systematics, 31, 533–563. Moritz C, Hoskin CJ, MacKenzie JB et al. (2009) Identification and dynamics of a cryptic suture zone in tropical rainforest. Proceedings of the Royal Society B-Biological Sciences, 276, 1235– 1244. Moussalli A, Moritz C, Williams SE, Carnaval AC (2009) Variable responses of skinks to a common history of rainforest fluctuation: concordance between phylogeography and palaeo-distribution models. Molecular Ecology, 18, 483–499. Nylander JA (2004) MrModeltest v. 2.0. Program distributed by the author. Evolutionary Biology Centre, Uppsala University. Phillips SJ, Dudik M (2008) Modeling of species distributions with Maxent: new extensions and a comprehensive evaluation. Ecography, 31, 161–175. Phillips SJ, Anderson RP, Schapire RE (2006) Maximum entropy modeling of species geographic distributions. Ecological Modelling, 190, 231–259.

2544 R . C . B E L L E T A L . Pritchard JK, Stephens M, Donnelly P (2000) Inference of population structure using multilocus genotype data. Genetics, 155, 945–959. Prychitko TM, Moore WS (1997) The utility of DNA sequences of an intron from the beta-fibrinogen gene in phylogenetic analysis of woodpeckers (Aves: Picidae). Molecular Phylogenetics and Evolution, 8, 193–204. Rosenberg NA (2004) DISTRUCT: a program for the graphical display of population structure. Molecular Ecology Notes, 4, 137–138. Rozen S, Skaletsky HJ (2000) Primer3 on the WWW for general users and for biologist programmers. In: Bioinformatics Methods and Protocols: Methods in Molecular Biology (eds Krawetz S, Misener S), pp. 365–386. Humana Press, Totowa, NJ. Saint KM, Austin CC, Donnellan SC, Hutchinson MN (1998) C-mos, a nuclear marker useful for squamate phylogenetic analysis. Molecular Phylogenetics and Evolution, 10, 259– 263. Schneider C, Moritz C (1999) Rainforest refugia and evolution in Australia’s Wet Tropics. Proceedings of the Royal Society of London Series B-Biological Sciences, 266, 191–196. Schneider CJ, Cunningham M, Moritz C (1998) Comparative phylogeography and the history of endemic vertebrates in the wet tropics rainforests of Australia. Molecular Ecology, 7, 487–498. Schneider CJ, Smith TB, Larison B, Moritz C (1999) A test of alternative models of diversification in tropical rainforests: ecological gradients vs. rainforest refugia. Proceedings of the National Academy of Sciences, USA, 96, 13869–13873. Shoo LP, Storlie C, Williams YM, Williams SE (In press) Potential for mountaintop boulder fields to buffer species against extreme heat stress under climate change. International Journal of Biometeorology, doi: 10.1007 ⁄ Soo484009-0286-4. Stephens M, Smith NJ, Donnelly P (2001) A new statistical method for haplotype reconstruction from population data. American Journal of Human Genetics, 68, 978–989. Swofford DL (2003) PAUP*: Phylogenetic Analysis Using Parsimony (* and Other Methods), v. 4.0. Sinauer Associates, Sunderland, MA. VanDerWal J, Shoo LP, Williams SE (2009) New approaches to understanding late Quaternary climate fluctuations and refugial dynamics in Australian wet tropical rain forests. Journal of Biogeography, 36, 291–301. Voelker G, Outlaw RK, Bowie RCK (2010) Pliocene forest dynamics as a primary driver of African bird speciation. Global Ecology and Biogeography, 19, 111–121. Whitmore TC, Prance GT, eds. (1987) Biogeography and Quaternary History in Tropical America. Oxford University Press, Oxford UK. Wiens JJ (2004) Speciation and ecology revisited: phylogenetic niche conservatism and the origin of species. Evolution, 58, 193–197. Williams SE, Hilbert D (2006) Climate change threats to the biodiversity of tropical rainforests in Australia. In: Emerging

Threats to Tropical Forest. (eds Laurance WF, Peres C.), Chicago University Press, Chicago, IL. Williams SE, Pearson RG (1997) Historical rainforest contractions, localized extinctions and patterns of vertebrate endemism in the rainforests of Australia’s wet tropics. Proceedings of the Royal Society of London Series B-Biological Sciences, 264, 709–716. Williams SE, Bolitho EE, Fox S (2003) Climate change in Australian tropical rainforests: an impending environmental catastrophe. Proceedings of the Royal Society of London Series BBiological Sciences, 270, 1887–1892. Williams JW, Jackson ST, Kutzbacht JE (2007) Projected distributions of novel and disappearing climates by 2100 AD. Proceedings of the National Academy of Sciences, USA, 104, 5738–5742. Williams SE, Shoo LP, Isaac JL, Hoffmann AA, Langham G (2008) Towards an integrated framework for assessing the vulnerability of species to climate change. Plos Biology, 6, 2621–2626. Williams SE, VanDerWal J, Isaac J et al. (2010) Distributions, life history specialization and phylogeny of the rainforest vertebrates in the Australian Wet Tropics. Ecology, in press. Wright SJ, Muller-Landau HC, Schipper J (2009) The future of tropical species on a Warmer Planet. Conservation Biology, 23, 1418–1426. Yeates DK, Bouchard P, Monteith GB (2002) Patterns and levels of endemism in the Australian wet tropics rainforesty: evidence from flightless insects. Invertebrate Systematics, 16, 605–661. Zink RM, Barrowclough GF (2008) Mitochondrial DNA under siege in avian phylogeography. Molecular Ecology, 17, 2107– 2121.

Supporting Information Additional supporting information may be found in the online version of this article. Data S1 Detailed statistical methods and results for morphological analyses. Table S1 Sampling localities Table S2 Primer sequences and PCR conditions Table S3 Variable weights for palaeomodels Table S4 Mitochondrial and nuclear diversity for topographically defined regions Tables S5 Pairwise estimates of divergence and FST for topographically defined regions Figs. S1 and S2 Individual gene trees of phased haplotypes. Please note: Wiley-Blackwell are not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

 2010 Blackwell Publishing Ltd