DOI: 10.5061/dryad.3xsj3txp7
title: Shrub cover declined as indigenous populations expanded across southeast Australia

description: 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.
description: <strong>Vegetation reconstruction</strong>
 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 &gt;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.
 <strong>Human population</strong>
 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).
 <strong>Biomass burning</strong>
 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.
 <strong>Palaeoclimate</strong>
 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. <strong>Generalized linear
 modelling</strong> 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).
description: # Shrub cover declined as Indigenous populations expanded across southeast
 Australia
 [https://doi.org/10.5061/dryad.3xsj3txp7](https://doi.org/10.5061/dryad.3xsj3txp7) ## Description of the data and file structure ### 'REVEALS data.zip' The archive 'REVEALS data.zip' contains all the vegetation cover data from all sites (by taxa) used in the regional compilation. This is a .zip folder with 33 Excel spreadsheets with Holocene land cover for 19 plant taxa recorded every 200 years (column = 'Age*'*). The land cover % values were obtained by applying the REVEALS model (see Methods section) on raw pollen data from the Indo-Pacific database (see Herbert et al., 55). Cells with no data correspond to time intervals (200 year bins) with no pollen counts available. 'Age' =  time in cal yr BP 'TR'  = median tree cover % 'SH' =  median shrub cover % 'HG' = median herbs and grasses cover % ### 'Compilations\_results.xlsx' The file 'Compilations_results.xlsx' contains the output compilation data published in this work (1,000-11,700 cal yr BP). This Excel spreadsheet is organised in four worksheets: 'Population', 'Charcoal', 'Palaeoclimate', 'Vegetation'. In all worksheets: 'Age' =  time in cal yr BP **The 'Population' worksheet contains:** 'SPD(raw) ' = raw summed probability density of archaeological data 'SPD(exponential logistic model)' = modelled summed probability density using an exponential-logistic model (see Methods section) **The 'Charcoal' worksheet contains:** 'Charcoal influx (z-scores) ' = composite charcoal record in z-scores (see Methods section for compositing methodology) **The ' Palaeoclimate' worksheet contains:** 'Composite (z-scores)' = composite palaeomoisture record in z-scores (see Methods section for compositing methodology) 'Upper' = upper limit of the palaeomoisture composite GAM 'Lower' = lower limit of the palaeomoisture GAM **The 'Vegetation' worksheet contains:** 'TR'  = median tree cover % 'SE_TR' = standard error for tree cover % 'SH' =  median shrub cover % 'SE_SH' = standard error for shrub cover % 'HG' = median herbs and grasses cover % 'SE_HG' = standard error for herbs and grasses cover % ### 'FEE2022REVEALS\_Data to share.xlsx' The file 'FEE2022REVEALS_Data to share.xlsx' contains the **post-colonial vegetation cover %** data previously published in Mariani et al. (2022). 'Age' =  time in cal yr BP 'TR'  = median tree cover % 'SH' =  median shrub cover % 'HG' = median herbs and grasses cover % ## Sharing/Access information All raw data used in this analysis are freely available through the following online sources: * Archaeological data from SahulARCH see Saktura et al. (54) - ([https://ro.uow.edu.au/data/71/](https://ro.uow.edu.au/data/71/)) * Pollen data from the Indo-Pacific Pollen database see Herbert et al. (*55*) - ([https://www.neotomadb.org/data/constituent-databases](https://www.neotomadb.org/data/constituent-databases)) * Charcoal data from the Global Palaeofire Database - ([https://www.paleofire.org](https://www.paleofire.org)) * Palaeoclimate data from the NOAA Paleoclimatology Database - ([https://www.ncei.noaa.gov/products/paleoclimatology](https://www.ncei.noaa.gov/products/paleoclimatology)) and for Barr et al. (*56*) - ([https://figshare.com/s/b4b5431fd9577afd95ef](https://figshare.com/s/b4b5431fd9577afd95ef)) see Table S3. ## Code/Software The following R codes used for the compilations were uploaded: 'Human population_SPD models code for pub.R' = contains the R code used to generate the summed probability densities to track human population expansion. 'GLMs code for pub.R' = contains the R code used to generate the generalised linear model for shrub cover % in the absence of human population. Palaeoclimate compilation.R' =  contains the R code used to obtain the palaeomoisture composite. 'Vegetation cover compilation.R' = contains the R code used to obtain a regional compilation of the land cover % data. Charcoal compilation_R code extension for Mariani Science ms.R =  contains the R code used to obtain the charcoal influx compilation. No new packages were generated.