Royal Society Publishing

Climate, not Aboriginal landscape burning, controlled the historical demography and distribution of fire-sensitive conifer populations across Australia

Shota Sakaguchi, David M. J. S. Bowman, Lynda D. Prior, Michael D. Crisp, Celeste C. Linde, Yoshihiko Tsumura, Yuji Isagi


Climate and fire are the key environmental factors that shape the distribution and demography of plant populations in Australia. Because of limited palaeoecological records in this arid continent, however, it is unclear as to which factor impacted vegetation more strongly, and what were the roles of fire regime changes owing to human activity and megafaunal extinction (since ca 50 kya). To address these questions, we analysed historical genetic, demographic and distributional changes in a widespread conifer species complex that paradoxically grows in fire-prone regions, yet is very sensitive to fire. Genetic demographic analysis showed that the arid populations experienced strong bottlenecks, consistent with range contractions during the Last Glacial Maximum (ca 20 kya) predicted by species distribution models. In southern temperate regions, the population sizes were estimated to have been mostly stable, followed by some expansion coinciding with climate amelioration at the end of the last glacial period. By contrast, in the flammable tropical savannahs, where fire risk is the highest, demographic analysis failed to detect significant population bottlenecks. Collectively, these results suggest that the impact of climate change overwhelmed any modifications to fire regimes by Aboriginal landscape burning and megafaunal extinction, a finding that probably also applies to other fire-prone vegetation across Australia.

1. Introduction

During the last glacial period, the Australian continent experienced a series of dramatic environmental changes. Global climate oscillation increased aridity, causing inland lakes and rivers to dry, sand dunes to become active and interior deserts to expand [1]. Fire is also a critical factor affecting vegetation and ecosystems across Australia. Although lowered CO2 and rainfall reduced the biomass of grass fuel, leading to a lower fire intensity than today [2], fire activity would have been high in areas that were seasonally wet and dry, especially in the northern monsoonal tropics [3,4]. Moreover, since human arrival in Australia (ca 50 kya) [5], landscape burning by Aboriginal people has been linked to significant changes in the geographical range and demographic structure of many vegetation types [6,7]. Indeed, burning is postulated as the key driver of the mass extinction of megafauna in Australia around 45–50 kya through altering vegetation to fire-adapted plant communities [8]. On the other hand, there is a lack of congruence between human activity and fire records during 20–40 kya, a period of consistently low Aboriginal populations [9] and large changes in fire regimes [10], and recent studies argue for significant roles of climate change on vegetation shift and megafauna extinction [1114]. Furthermore, studies of sedimentary records from the humid tropics and southeastern Australia suggest an ecological feedback, where relaxed herbivore pressure following extinction of megaherbivores resulted in a switch to flammable sclerophyllous vegetation, which increased fire activity [15,16]. Thus, there is still a controversy concerning whether climate or fire impacted vegetation more strongly, and the role of fire regime changes owing to human activity and megafaunal extinction during the last glacial period.

One of the most critical problems in resolving this uncertainty is the limited availability of palaeoecological records with sufficient coverage of space and time. Most palaeoproxies representing vegetation and fire regime changes in the last glacial period have come from the temperate or tropical fringes of the continent [10,17], where sedimentary cores are relatively well preserved. In the humid tropics, during the Last Glacial Maximum (LGM; ca 20 kya) there was a reduction in rainforest angiosperms and a transition to sclerophyllous woodland, possibly caused by reduced rainfall and accelerated by Aboriginal burning [18]. Southern temperate regions saw an expansion of drought-tolerant communities dominated by Myrtaceae, Asteraceae and Chenopodiaceae in response to cooler and drier conditions [1921]. In the vast area of the arid zone, however, the palaeovegetation data come mainly from shallow time-depth records at salt lakes and springs, which do not provide a continuous record of vegetation change throughout the LGM [22,23]. This is especially the case in Central Australia and the Pilbara region, where reconstruction of the LGM vegetation relies on limited phytolith and charcoal deposits at a few research sites [24,25]. Therefore, to understand the impacts of climate change and fire, and to tease apart the effects of Aboriginal landscape burning and megafaunal extinctions, integrative modelling is required to fill the gaps between the spatially and temporally sparse palaeontological records.

Here, we present a continent-wide phylogeography and population demographic analysis, coupled with palaeodistribution reconstruction, of a widespread conifer clade, the Callitris columellaris species complex (Cupressaceae) [26]. The species complex is an ideal model system given its geographical range spanning the monsoon tropical, temperate and arid zones of Australia (figure 1a,b). Although the species complex includes one of the world's most drought-tolerant plants [27], successive wet years are necessary for seedling recruitment and population maintenance [28]. As the populations are also sensitive to fire [29,30] and mammalian herbivory [31], their population demography records environmental changes [31]. Reconstruction of the Callitris populations’ histories could, therefore, serve as a proxy for regional vegetation responses to the varying environmental changes through time and across the continent.

Figure 1.

(a) The climate zones of Australia with the biogeographic regions defined for this study (see more details in the caption of electronic supplementary material, table S1). (b) Distribution records of the C. columellaris species complex, superimposed on the predicted distribution at present. (c) The predicted distribution of each species at the LGM. The map for the LGM was created by extracting the consensus pixels where the each species’ presence was supported under all of the three climate simulations, with thresholds of the maximum training sensitivity plus specificity used. The predictions at the Nullarbor Plain are masked as not being filled by any Callitris species, reflecting the edaphic unsuitability of this arid limestone plain.

We examine the alternative hypotheses concerning the effect of either climate or human-caused fire disturbance on Callitris distributions during the last glacial period. Under the hypothesis that climate had a major role in Callitris distribution change, it can be expected that the distribution ranges in the arid zone would have contracted or even retreated to coastal regions because of desertification under decreased rainfall during the LGM [32,33]. On the other hand, the human-caused burning hypothesis predicts that, since human arrival, the fire-prone Callitris populations would have collapsed most severely in the monsoon tropics, where grassy fuels and human activity were the highest [9]. In this study, the alternative hypotheses are tested based on three approaches: reconstruction of palaeodistributions using species distribution modelling; phylogeographic analysis of relationships among regional populations and explicit modelling of historical demography of populations based on genetic data. Geographical patterns of estimated demographic changes, temporal changes in demography, coupled with the reconstructed palaeodistribution patterns, are then compared with the expectations under the alternative hypotheses to assess the relative importance of environmental factors.

2. Material and methods

(a) Palaeodistribution reconstruction

The study group comprises the morphospecies C. columellaris F. Muell., C. glaucophylla J. Thompson & L. A. S. Johnson, C. intratropica R. T. Baker & H. G. Sm., C. gracilis R. T. Baker and C. verrucosa (Endl.) F. Muell. s.l. We followed the taxonomy of Hill [34], except that we included C. tuberculata R. T. Baker & H. G. Sm. and C. preissii Miq. in C. verrucosa because these three taxa appeared to intergrade. To construct the species distribution model for the five morphological species, the distribution records (figure 1b and electronic supplementary material, S1) were collated with bioclimatic and topographic variables using the software Maxent v. 3.3.3e [35]. Six temperature and precipitation variables (Bio1, Bio12, Bio16, Bio17, Bio18 and Bio19) were chosen as likely to influence the distribution range of the species complex [36], and the climate layers in 2.5 arc-min resolution were downloaded from WorldClim ( The variable slope, which is potentially correlated with Callitris distribution through its association with topographic fire protection [30], was calculated based on 2.5 arc-min digital elevation model using ArcGIS (Esri, USA). Modelling was performed for each species using the same set of variables, with the background designated as Australia's overall land area during the time period modelled and the other settings held as default. The area under the receiver operating characteristic (ROC) curve (area under the curve, AUC) was then calculated for each run, which is an indicator of the accuracy of model prediction [37]. The established models were projected onto the climatic condition at present and the LGM. The LGM climatic layers were created following the method in Sakaguchi et al. [38] based on the simulation summaries of three global climate models (CCSM3.0 [39], FGOALS-g1.0 [40] and MIROC3.2 [41]), which are available at the Paleoclimate Modelling Intercomparison Project site (

(b) Plant material, microsatellite genotyping and DNA sequencing

We collected foliage samples of 1434 trees from 62 populations (see electronic supplementary material, table S1), well representing the known distributions. Genomic DNA was extracted from ca 1.0 cm of dried foliage tissue using a slightly modified hexadecyltrimethylammonium bromide (CTAB) method [42]. Thirty nuclear-expressed sequence tags—simple sequence repeat (EST-SSR) markers applicable to the whole species complex [43]—were used to genotype all the samples (see electronic supplementary material, table S2). The PCR was carried out following the standard protocol of the Qiagen multiplex PCR kit (Qiagen, Hilden, Germany) in a final volume of 10 µl. PCR products were loaded onto an auto sequencer (3100 genetic analyser; Applied Biosystems, USA) to assess fragment lengths using Genotyper software (Applied Biosystems). In addition, a chloroplast fragment of trnL (UAA) 5′ exon-trnF (GAA) [44], which are confirmed to be inherited paternally for these species [45], was PCR amplified and sequenced for 314 samples (5.1 samples per population) following the protocol of Sakaguchi et al. [45]. Sequencing of PCR products was conducted using both forward and reverse primers with an ABI Prism BigDye terminator cycle sequencing ready reaction kit v. 3.1 (Applied Biosystems), and electrophoresed on an ABI Prism 3100 genetic analyser. The chloroplast sequence data were edited and aligned using bioEdit v. [46].

(c) Population genetic diversity and structure analysis

Genetic differentiation among all populations at the chloroplast locus was estimated with GST [47] using the software Contrib v. 1.02 [48]. A median-joining network with Callitris endlicheri (GenBank: AB723699) as an outgroup was inferred using Network v. [49], with the indels coded as the fifth state.

Allelic richness [50] and private allelic richness at 30 EST-SSR loci were calculated after rarefaction of each population to 20 gene copies using HP-rare1.1 [51] and expected heterozygosity [52] using FSTAT v. 2.9.3 [53]. Two summary statistics, mean allele size variance across loci and mean M index across loci [54], which can capture past population demographic changes, were calculated using DIYABC v. beta [55,56]. Genetic differentiation at 30 EST-SSRs was evaluated by GST [57] and FST [58] for each locus and by GST for overall samples using MSA v. 4.05 [59]. Population structure based on the nuclear markers was examined by both distance- and model-based methods. Pairwise genetic distances of Nei's DA [60] among all the populations were calculated using Populations v. 1.2.30 beta [61]. This matrix was then used to reconstruct a split network [62] using SplitsTree4 v. 4.10 [63,64]. The model-based Bayesian algorithm, Structure v. 2.3 [65] was used to detect population structure, with the correlated allele frequencies model [66] and the sampling location specified as prior information [67]. Twenty-five independent simulations with 100 000 burn-in steps, followed by 100 000 Markov chain Monte Carlo steps, were performed for the whole dataset (K = 1–20) and the MCMC chain convergences were confirmed for each parameter. To determine the number of genetic clusters (K), the log-likelihood of the data L(K) and the second-order rate change of L(K) (ΔK) [68] are plotted against K assumed in the Structure analysis [65].

(d) Demographic analysis

Approximate Bayesian computation (ABC) implemented in DIYABC v. beta [55,56] was used to model the demographic history of Callitris populations through the last glacial period. ABC is a recent flexible class of computation algorithms designed to perform complex model-based inferences and can bypass the exact likelihood calculations by using summary statistics and massive computer simulations [56,69]. In this study, we assumed a simple demographic scenario, where the ancestral population size (NA) changes to current size (NC) at a time (t) during the last glacial period. Consideration of this scenario allows estimation of both population size expansion and any bottleneck, and the rate of size change for each population. The alternative is a ‘stable scenario’, where no population size change occurred since the last glacial period. One million simulations (0.5 million for each scenario) were performed by drawing the parameter values from prior distribution combinations, which were pre-evaluated by principal component analysis (PCA) plotting to confirm that the observed datasets were well surrounded by simulated datapoints. The prior distributions used here were as follows: population size (N) log-uniform [10, 10 000] with no constraints on the relationship between NC and NA; the time parameter (t) uniform [100, 10 000] in generations and sufficiently spans the last glacial periods by assuming a generation time for the species complex of 20–70 years (see discussion in the electronic supplementary material, S2). The prior distribution of mutation rates of the autosomal SSRs was assumed to be broad because EST-SSRs evolve more slowly than anonymous SSRs in general [70]; the mean mutation rate was uniformly distributed [1×10–5, 1×10–3] and individual mutation rate was also uniformly distributed [1×10–6, 1×10–3] (per generation). The Gamma distribution [1×10–1, 3×10–1] was specified for the geometric coefficient P, to generate the datasets following the generalized stepwise mutation model [71]. A series of summary statistics available for a single-population model in DIYABC v. beta (mean number of alleles, mean genetic diversity [52], mean size variance and mean M index across loci [54]) were used in the model checking (or assessment of the goodness-of-fit of a model to posteriors) [55] and estimation of posterior distribution of parameters. Model checking was performed by comparing the deviation of summary statistics of the posterior predictive distribution with the observed data. A posterior distribution of demographic parameters was estimated for each demographic scenario, based on 0.1% of the simulated data closest to the observed data via logit transformations. Likewise, logistic regression of 1% of simulated data closest to the observed was used to calculate the posterior probabilities of the scenarios to be compared. When the 95% confidence interval of posterior probabilities of the two scenarios did not overlap, we considered that the population scenario with the higher probability was statistically selected. Alternatively, if the posterior distributions did not significantly differ between the scenarios, or if significant deviations of simulated summary statistics from the observed were detected in both scenarios in model-checking procedures, then we conservatively considered them to be ‘indecisive’. The demographic analyses were performed for each population, and also for the populations grouped by the biogeographic regions and morphospecies. The latter analyses were motivated by a simulation study that showed that factors such as population substructure, genetic diversity and the sampling scheme play a role in generating false bottleneck signals, even for demographically stationary populations, and that this confounding effect can be countered by analysing the whole metapopulation to test whether a bottleneck is still detected [72].

While the ABC is a generally valid and practical approximation method for Bayesian inference, information loss owing to the approximation based on a limited number of suitable summary statistics could lead to poor parameter estimation and model choice with low confidence [73,74]. To confirm the validity of using ABC to analyse our dataset, we compared the posterior distributions of demographic parameters for the grouped population data, estimated under the population size change scenario by ABC computation in DIYABC v. beta and a full-likelihood Bayesian method implemented in msvar1.3 [75]. The msvar1.3 is shown to be efficient to detect past population decline and expansion, outperforming the other moment-based methods [76]. The demographic model with exponential growth/decline in msvar1.3 includes the demographic parameters of population size (NC, NA), the time parameter (t) and mutation rate (μ), which can vary among loci under the hyperpriors in the hierarchical model. The priors (mean, s.d.) in log-scale were set as follows: NC (3, 2), NA (3, 2), t (3, 1) and μ (−4.5, 1), and each Markov chain was run for 1×109 steps with parameter values recorded every 2×104 steps.

Although msvar1.3 is a powerful method for detecting past population size change, it is known to bias estimators of demographic parameters for the expanded populations [76]. In this study, we simulated 10 replicates of genotype data using DIYABC v. beta, in which 10-fold population expansion was assumed to have occurred under the following settings (Nc = 5000, NA = 500, t = 1000, mean mutation rate = 1×10–4 over 30 microsatellite loci), mostly compatible with the real dataset in Callitris populations. The simulated dataset was then analysed by msvar1.3 and DIYABC v. beta under the same prior combinations as above, and the posterior estimates of demographic parameters were checked to see whether the different algorithms could successfully recover the assumed population demography.

3. Results

(a) Palaeodistribution reconstruction

The species distribution models for the five morphological species showed high predictability, with test AUC values (25% of the records used for model testing) from 0.829 (C. glaucophylla) to 0.996 (C. columellaris s.s.). The predicted present distribution patterns were mostly concordant with the observed records (figure 1b). The LGM distribution of C. glaucophylla was predicted to have been severely fragmented and mostly confined to the major mountain ranges in arid Central Australia, the Pilbara and temperate parts of the continent (figure 1c). On the other hand, the potential LGM distribution was estimated to have been stable for C. columellaris s.s., C. gracilis and C. intratropica and to have shifted somewhat northward for C. verrucosa.

(b) Phylogeographic structure and population genetic diversity

Sequencing of the chloroplast fragment [44] for 314 trees from 62 populations recovered 12 haplotypes within the species complex (GST = 0.934; figure 2a and electronic supplementary material, table S3). The haplotype network showed three genetic groups, each separated by one hypothetical haplotype node, corresponding with the species delineation and geography; C. intratropica (monsoon tropics), C. columellaris s.s. (Central Eastern Coast) and the other three species occurring in the southern regions (see electronic supplementary material, figure S1a). Phylogeographic barriers were detected within the northern tropical region, demarcating the disjunct populations in the ranges of the Kimberley, the Top End and Cape York Peninsula (figure 2a). In the southern group, C. verrucosa’s haplotypes were distinguished from those belonging to C. glaucophylla and C. gracilis. A single haplotype (yellow) was widespread in the central to southwestern populations across climate zones and shared by C. glaucophylla and C. gracilis.

Figure 2.

Geographical distribution of (a) chloroplast haplotypes and (b) nuclear SSR genetic clusters estimated by Structure [65]. The C. gracilis populations are distinguished by the pie charts having dotted outlines and C. verrucosa populations are indicated by ‘V’ characters next to their pie charts. (c) A population network based on Nei's DA distance [60], with the proportions of ancestry of each population superimposed at the node tips. The size of pie charts is proportional to rarefied population genetic diversity. The characters indicate the sampling regions (figure 1a). Colour blobs indicate the species delineation based on morphological characters.

The phylogeographic structure estimated by Structure analysis [65] of genotype data from nuclear EST-SSR loci collected for 1434 trees from 62 populations (286 alleles in total; GST = 0.422; electronic supplementary material, table S2) was broadly consistent with that of the chloroplast but the greater variability of these markers resolved finer structuring of C. glaucophylla and C. gracilis in central to southern arid regions (figure 2b; see discussion on determining the number of clusters in the electronic supplementary material, figure S2). Notably, the disjunct arid populations in Central Australia, Pilbara and the southern regions were shown to form distinct genetic clusters, which were also separated from each other in the distance-based population network (figure 2c).

Population genetic diversity varied among regional populations and morphological species (figure 2c and the electronic supplementary material, figure S3). The populations in the southern half of the continent showed higher genetic diversity than the tropical regions occupied by C. intratropica, whereas the population group of C. glaucophylla in the isolated Pilbara region, the narrow-ranging C. columellaris s.s. and C. gracilis harboured the lowest levels of diversity. Populations in arid Central Australia, the Pilbara and southern semi-arid regions showed a large variance in allele size (see electronic supplementary material, figure S3c) and remarkably also low M index values (see electronic supplementary material, figure S3d), which are below the common threshold (M = 0.68) for a genetic bottleneck suggested by Garza & Williamson [54].

(c) Demographic analysis

The log-transformed posterior distribution of the ratio of NC and NA represents the intensity of population size change under the population size change scenario (figure 3a). Clearly, variation in population demography can be explained by the climate zones in Australia, with contracted populations in arid regions, stable in the monsoonal tropics and stable to expanded in temperate regions. While we acknowledge that scenario selection based on the ABC method should be considered with caution, owing to its approximation error in the posterior probabilities of the models [73], we detected 24 populations where the population size change scenario received significantly more statistical support than the stable scenario: the deviation of summary statistics of the posterior predictive distribution from the observed data tended to be smaller than those of stable scenario (see electronic supplementary material, table S4). The significantly bottlenecked populations (N = 9) were located mostly in the arid zone and currently temperate zone where LGM annual precipitation is estimated to be less than 250 mm (see electronic supplementary material, figure S4 and table S4). The other two bottlenecked populations were located in the climatically temperate zone in the southwest, on the geographical periphery of the species complex's current distribution. The posterior distribution of the time parameter (t) for these populations showed clear peaks around 30 kya (figure 3b), with broad tails spanning the timings of LGM and human arrival in Australia. The other 15 populations whose population sizes were shown to have increased were located exclusively within the southern temperate zones (figure 3 and electronic supplementary material, figure S3). By contrast, the monsoon tropical populations were inferred to have been stable throughout the last glacial period, except for one population that appeared to have expanded, though not strongly (figure 3).

Figure 3.

(a) Posterior probability density estimates for the intensity of changes in population sizes, represented as the log-transformed ratio of the current to the ancestral population size (NC/NA). Colours of density lines indicate the climate zone. Note that ‘arid’ includes the populations in the arid zone and in the currently temperate zone that received less than 250 mm annual precipitation during the LGM. (b) Posterior probability density estimates of the time parameter for the populations that are inferred to have experienced population size changes. The parameter was converted to absolute time by assuming that Callitris' generation time is 30 years. The two dotted lines represent the time envelopes when generation time was assumed as its minimum (20 years) and maximum (70 years), respectively. The broken line shows the prior distribution for the time parameter.

The population declines in the arid Pilbara and/or Central Australia regions were supported by the ABC and full-Bayesian analyses of the grouped data (see electronic supplementary material, figures S5 and S6), indicating that the signature of population bottlenecks in the single-population analyses is unlikely to be an artefact owing to population substructure within the regions. While the ABC analyses on the grouped data also showed compatible parameter estimates for regional populations and the single-population analyses, the temperate and monsoon tropics populations were estimated to have experienced slight population declines by msvar1.3 (see electronic supplementary material, figure S6), rather than stability or expansion as estimated by ABC computation. This discrepancy was shown by simulation-based analyses to have been resulted from biased parameter estimation in msvar1.3 for the expanded populations (see electronic supplementary material, figure S7), similar to a finding by Girod et al. [76]. In the analyses, the population size change ratios estimated by msvar1.3 showed sharp peaks around zero on the log-scale, whereas DIYABC v. beta recovered the 10-fold size expansions assumed in the data simulations. This led us to prefer the ABC analyses for not only detecting past population size change but also estimating demographic parameters.

4. Discussion

Given the slow mutation rate estimated for Cupressaceae species [77], the strong phylogeographic structure detected in the species complex should reflect long histories of regional lineages during the Pleistocene, dating back prior to the last glacial period. Our finding of distinct lineages in regional mountain ranges suggests that the Callitris species complex could have persisted in long-term refugia probably for multiple glacial–interglacial cycles during the late Pleistocene, a period of dramatic environmental change. Furthermore, multiple distinct lineages detected from nuclear loci in the arid interior indicate long-term resilience and persistence of Callitris populations in arid habitats. This finding supports the hypothesis that the drought-tolerant species survived in LGM refugia under an extremely arid climate, not only around the Puritjarra rock shelter in Central Australia [25] but also in the Pilbara and southern regions, contradicting the view that these arid ranges were treeless desert at that time [33].

Results from the ABC analyses clearly showed that the demographic change in Callitris populations was influenced by regionally differential climate change during the last glacial period. The majority of populations identified as being significantly bottlenecked during the period including the LGM using ABC analyses were located in the interior arid ranges. Consistently, our palaeodistribution reconstruction showed that the LGM distribution of C. glaucophylla in the arid zone contracted into regional mountain ranges. Although fire is one of the strongest influences on the current patchy distribution of C. glaucophylla in Central Australia [30], the estimated range/population contraction here around the LGM was probably caused by reduced rainfall rather than fire. A phytolith spectrum showed that grasses sharply declined in abundance at ca 40 kya and were mostly eliminated at the LGM [24], owing to a weakened influence of the summer monsoon and subsequent decrease in the frequency and intensity of fires in this region [78]. Therefore, Callitris populations in the arid Pilbara and continental interior were probably confined to microrefugia in canyons, on mountain ranges and inselbergs, where moisture stress is the least and orographic rainfall may occur. Confinement to such microrefugia through the LGM would have reduced the genetic diversity of arid populations by stochastic processes, which may explain the very low genetic diversity and M index of current Pilbara Callitris populations. Consistent with climate being the predominant control on Callitris distribution are pollen records from arid sites near Lake Frome and Lake Eyre, which suggest that Callitris woodlands rapidly increased soon after the LGM, and became most prominent around 10 kya as the climate became less arid [32].

In contrast to the arid populations, genetic signals of population expansion were detected in southern temperate populations. In the southeastern highlands of Australia, Callitris pollen has been recorded, even during the LGM, and it attained its maximum representation during the Holocene [21]. These palynological data on Callitris population expansion during the Holocene harmonize with our results from the demographic analysis as well as predictions from species distribution models that the populations persisted in the temperate zone during the glacial period. The population size increases in the temperate regions may have been accelerated by relaxed herbivory following megafaunal extinctions. Although fossil evidence that extinct mammals ate Callitris is inconclusive [79], some introduced mammalian herbivores (sheep, goats and rabbits) can cause recruitment bottlenecks [31]. Our study finds no evidence of adverse effects on temperate Callitris populations, for example from putative increased burning resulting from fuel build-up following megafaunal extinction [15]. If this had been the case, then the increased fire frequency should have genetically impacted Callitris Australia-wide, whereas we detected regionally differential demographic changes and these changes matched the modelling of species' distributions.

There remains controversy over the effects of Aboriginal landscape burning on the Australian biota [6]. This subject has been researched intensively for the fire-sensitive C. intratropica in the northern monsoon tropics, the region with the highest fire risk and longest interactions between humans and ecosystems in Australia. From charcoal and pollen records, the decline in fire-sensitive vegetation, including C. intratropica, and its replacement with more fire-tolerant vegetation in the tropical regions during the late Pleistocene, has been attributed to both Aboriginal landscape burning since ca 50 kya and long-term climatic fluctuation [18]. The demographic analysis of C. intratropica populations in this study showed that almost all populations have been stable since the beginning of the last glacial period, suggesting that the impact of fire during the period of Aboriginal colonization was not sufficiently severe to leave a genetic signature of bottlenecking in these populations. However, changed fire regimes in the present day are exerting differential effects on the demography of C. intratropica in northern Australia [29]. High intensity, frequent fires since European settlement have destroyed many populations in the monsoon tropics and arid zone, restricting surviving populations to sites with topographic fire protection [6]. By contrast, under frequent, low-intensity burning characteristic of Aboriginal management, healthy populations of C. intratropica are maintained [80]. The mechanism for the latter is a vegetation–fire feedback where C. intratropica groves with closed-canopies suppress grassy fuels, thereby enabling the persistence of patches of fire-sensitive Callitris woodland in highly flammable savannahs with heavy grass fuel loads, provided that the fires are of low intensity [81]. We therefore propose that the traditional Aboriginal landscape burning in the northern monsoon tropics, and possibly southern Australia, helped to maintain Callitris populations that may otherwise have contracted into fire-protected refugia during the LGM. Our results from Callitris provide evidence in direct opposition to the view that Australia was subjected to an ‘ecosystem collapse’ owing to sustained Aboriginal burning during the last glacial period [8].

Data accessibility

The EST-SSR genotype data are deposited at Dryad: doi:10.5061/dryad.6j777.

Funding statement

This research was partly supported by a Grant in Aid for the Japan Society for the Promotion of Science Fellows (no. 256059), a Grant in Aid for Scientific Research (no. 24248028) from the Ministry of Education, Culture, Sports, Science and Technology, Japan, and Commonwealth Environment Facilities Research Fund grant no. B0016193.


The authors are grateful to L. G. Cook, J. R. P. Worth and M. Yamasaki for their helpful discussions on this study. We also thank many people, especially P. Moser, L. McCaw, I. Radford, M. King and S. Nichols, for sample collection. We are grateful to the Northern Territory Parks and Wildlife Service, Kakadu National Park, NSW State Forests, NSW Department of Environment and Climate Change, Queensland Department of Environment and Resource Management, SA Department of Environment and Heritage, WA Department of Environment and Conservation, Victorian Department of Sustainability and Environment, Anindilyakwa Land Council, Northern Land Council, Tiwi College, the Arid Recovery Reserve, the Australian Wildlife Conservancy and many private landholders for help with site selection and permission to sample on their land.

  • Received August 21, 2013.
  • Accepted October 7, 2013.


View Abstract