Wildfires in forests globally have become more frequent and intense due to changes in climate and human management. Shrub layer fuels allow fire to spread vertically to forest canopy,  creating high-intensity fires. Our research provides a deep-time perspective on shrub fuel loads in fire-prone southeastern Australia. Comparing 2,833 records for vegetation cover, past climate, biomass burning, and human population size across different phases of human occupation, we demonstrate that Indigenous population expansion and cultural fire use resulted in a 50% reduction in shrub cover, from approximately 30% from the early-mid Holocene (12-6 ka) to 15% during the late-mid Holocene (6-1 ka). Following British colonization, shrub cover has increased to the highest ever recorded (mean of 35% land cover), increasing the risk of high-intensity fires.

Vegetation reconstruction

We applied REVEALS (32) to quantify past vegetation using 31 pollen records in southeast Australia (Table S1), covering the Holocene (n=29) and Last Interglacial periods (n=2). The REVEALS model was run to convert raw pollen data (counts) into estimates of land cover (%) by correcting for biases in pollen production (i.e. different plant species produce different amounts of pollen) and pollen dispersal (i.e. dispersal patterns differ in response to pollen grain properties) (32). Pollen productivity estimates (PPEs) required for REVEALS for the 19 most abundant pollen taxa were derived from Mariani et al. These 19 taxa cover a large proportion of the vegetation and pollen counts, for example, across >275 vegetation quadrats for the state of Victoria (5), a median of 81% were target taxa (± 24%). Based on the modern (moss) pollen counts, a median of 96% of moss samples were made up by target taxa (±22%). For fossil samples across the region, these 19 target taxa constitute approx. 75% of pollen counts (median value, ±14%). The missing % are usually within the Proteaceae and Fabaceae families (shrub layer), which means our reconstruction of shrub cover is on the conservative side and reconstructed values might be slightly higher. REVEALS was executed using the R package disqover version 0.9.09 accessible through (https://github.com/MartinTheuerkauf/disqover.), using the Lagrangian stochastic model (LSM) for pollen dispersal parameterized as previous studies in Europe and Australia (5, 58–61). The median across all REVEALS results combining each site per 200 year time bins were compiled using MS Excel and the R code provided. The results of the Holocene estimates were compared with post-colonial vegetation estimates previously published (n= 51; Mariani et al.).

Statistical analyses were undertaken to assess significance of shrub cover % change amongst the throughout the reconstruction period. A pairwise t-test was conducted on square-root transformed percentage data for shrub cover (Table S1). Square-root transformation was required prior to the t-test, as the dataset did not have a normal distribution and parametric tests assume normality. The square-root transformation provided a normal distribution (Figure S11a,b) and the autocorrelation of the shrubs time series is presented in Figure S11c. To further support the results from the t-test, we undertook a Kruskal-Wallis test (62) with Wilcoxon pairwise comparisons (63). This test is non-parametric, hence the raw data (non-normal) were supplied and results are presented in Figure S11d.

Human population

A total of 2,358 radiocarbon ages of archaeological evidence of past human occupation across southeast Australia were used in this study. Initially 6,522 radiocarbon ages were extracted from the SahulArch database, accessible through OCTOPUS v.2 (https://octopusdata.org/); Codilean et al. (936 and subject to screening for region of interest, Holocene age range, and appropriate associations with archaeological deposit (if context was indicated as ‘sterile’ or ‘non-occupation’ by original study, these were excluded). The resulting 2,368 radiocarbon ages were then calibrated using SHCal20, and the summed probability density (SPD) of calibrated ages was calculated using the thinning approach in the rcarbon package using version 1.5.1 to infer past human population changes. We acknowledge that the number of radiocarbon dates and associated uncertainties can influence summed probability estimates, especially for estimates aimed at detecting short-term variations and rates of change. Our study is focused on long-term multi-millennia-scale changes and does not consider rates of change. Bayesian bounded population growth models were further used to assess the fit of the SPD of radiocarbon ages using the nimbleCarbon package version 0.2.5 in R. The models were fitted through Markov Chain Monte-Carlo and ranked using the Watanabe–Akaike information criterion (WAIC). The fitted population growth models include exponential, double exponential, and exponential logistic models (Fig. S4; see Crema and Shoda, for method details). Among the fitted models for the SPD, the exponential logistic model ranked as the top model with the lowest WAIC, while the double exponential model ranked the lowest with the highest WAIC (Fig. S4).

Biomass burning

Sedimentary charcoal records from 108 lakes and wetlands across southeastern Australia were collated from the Global Charcoal Database and Neotoma (Table S2). Charcoal concentrations were converted to charcoal accumulation rates using the available chronological information. As elsewhere in the world, researchers in Australia have used various methods to quantify charcoal, necessitating data transformations to extract regional-scale palaeofire trends. Because these transformations tend to mask inter-site variability, we grouped records by analysis method prior to min-max rescaling of all records within each group and square-root transformation. As woody charcoal is preferentially preserved and accumulated in the sedimentary records, as opposed to grass charcoal (38), we can interpret our charcoal influx trends as woody biomass burned. The method is fully described in Mariani et al. and Rowe et al.

Palaeoclimate

Five Holocene terrestrial palaeomoisture records (see Table S3) were compiled in this investigation comprising three precipitation/evaporation (P/E) proxy reconstruction types. Lake level reconstructions where lower lake levels represent diminished lake recharge into closed lakes through precipitation. Palaeosalinity records where higher salinities indicate periods of increasing lake desiccation and reduced regional recharge (drier conditions). A mean annual rainfall reconstruction for Swallow Lagoon, Stradbroke Island, directly reconstructing precipitation. Whilst the palaeoclimate records span a large climatic gradient (Fig. S2), the individual trends were found to be coherent, justifying their compositing (Fig. S6). All data were individually z-scored and then merged using MS Excel, before GAM smoothing (k=100), to create a composite record for the southeastern Australian region.

Lake level reconstructions from sediment grain size analysis at the maar crater Lakes Keilambete and Gnotuk were used as indicators of regional precipitation patterns from westerly wind circulation. Both lakes have relatively simple hydrological inputs and are considered good indicators of evaporation-precipitation oscillations. Palaeosalinity reconstructions as Total Dissolved Solids (TDS g/L) for Lake Keilambete were originally procured through ostracod Modern Analogue Technique (MAT) reconstructions using an analogue database of 491 samples. A depauperate fossil record (n=3 species) of ostracods at Lake Gnotuk meant only grain size lake levels were included in the index. Palaeosalinity reconstructions (Log g/L-1) for Lake Jacka and NW Jacka (72) were calculated from ostracod assemblages using a weighted-averaging transfer function with 119 modern analogue samples. The main control on modern assemblage composition was the total salinity of the lakes. At Blue Lake, palaeosalinity (Log10TDS g/L) were originally derived from the Weighted Averaging of Modern Analogues Technique (WMAT) utilizing 534 analogue samples. The rainfall reconstruction at Swallow Lagoon (56) was included to consider ENSO-derived palaeomoisture signals. The rainfall reconstruction of Swallow Lagoon used δ13C ratios from ancient Melaleuca quinquenervia leaf fragments as a proxy of historical rainfall (mm), calibrated against a 12-year monthly record of Melaleuca quinquenervia litter δ13C and recorded rainfall. The record reflects mean annual rainfall for the total ENSO system rather than El Niño/La Niña events, where the 1cm samples represent an average of 24.4 years of data.

Generalized linear modelling

Generalized linear model (GLM) was used to identify the main driver(s) of ladder fuels (shrub) cover with palaeoclimate synthesis, biomass burned, and SPD of archaeological ages set as predictors. Variables were randomly sampled 100 times without replacement at a lower resolution to remove the effect of autocorrelation. The lower resolution includes 30, 60 and 90% of datasets to check the consistency of the results (Table S3). A separate GLM was also fitted to predict mid-late Holocene shrub cover changes using early-mid Holocene GLM (predictors: human population –SPD, palaeoclimate index and charcoal influx) as a training set. The mid-late Holocene model was intended to reflect shrub cover changes under the scenario of no human influence and so, only palaeoclimate index and charcoal influx were included as predictors. GLMs were fitted using the MASS package version 7.3 in R (74).