Introduction

The Australian Black Summer wildfires of 2019/2020 were unusual in both size and intensity, burning an estimated 30 million hectares across Australia (Dickman 2021). Most of these fires were concentrated in the eastern portions of New South Wales, Queensland, and Victoria. Due to the scale and severity of these fires—with an estimated 3 billion terrestrial vertebrates killed or displaced (van Eeden et al. 2020; Legge et al. 2021)—the Commonwealth Threatened Species Commission sought rapid conservation assessments for animal species that were potentially impacted by the wildfires. While most of these conservation efforts focused on vertebrates, a report released by the Wildlife and Threatened Species Bushfire Recovery Expert Panel (DAWE 2020) included 191 invertebrate species in need of urgent conservation, with an additional 147 species flagged for additional conservation assessment.

Traditional methods of conservation assessment, such as those used by the IUCN that calculate an Area of Occupancy (AOO) and Extent of Occurrence (EOO), rely on rich historical datasets and in-depth knowledge about species’ former and current ranges, information which is often severely lacking for most invertebrate species. Conservation assessments of the priority invertebrate species identified by DAWE (2020) using AOO and EOO were inconclusive (Cassis et al. 2022; Foon et al. 2022; Reid et al. 2022). A supplement to conservation of species is conservation of their genetic diversity. Population genetics analyses can provide concise information on whether any threatening process has adversely affected the genetic variation within and between localities. This is essential for short-term avoidance of inbreeding, and medium- to long-term adaptation to continual environmental change (Allendorf and Luikart 2007; Frankham et al. 2010).

We hypothesised that the wildfires would have a detrimental impact to the genetic structure of invertebrate species due to the direct degradation of habitat, habitat fragmentation, and the deaths of countless individual invertebrates. We selected the flightless dung beetle Amphistomus primonactus Matthews 1974 as our target; this species was chosen because of its relative abundance where present, its low vagility, its limited geographic range in the centre of the north-eastern forests of New South Wales—an area highly impacted by the 2019/2020 wildfires—and its supposed restriction to rainforest (Matthews 1974), a habitat type with historical resilience to fire, but which was burnt to the largest extent ever recorded in the 2019/2020 wildfires (Legge et al. 2021). This species was additionally well-represented in existing museum collections, and well sampled throughout its range. These factors mean that A. primonactus is relatively well studied compared to most other insect species, and additionally represent biological traits that are highly sensitive to ecological disturbance.

To test our hypothesis, we examined the genetic diversity of A. primonactus within populations that experienced differing burn intensities during the 2019/2020 wildfires. We also examined the genetic diversity between populations of A. primonactus that were separated by habitat that had been burnt at varying intensities to determine whether habitat fragmentation caused by burnt landscapes had an impact on the migration of the species. As well as these contemporary investigations, it would be ideal to have a temporal comparison, so we attempted to genotype A. primonactus collected prior to the Black Summer wildfires in 1993 (Ferrier et al. 1999); many 1993 sites were revisited in 2021 (Reid et al. 2022). However, unfortunately, the 1993 samples could not be genotyped (see Results).

Methods

Specimen sampling methods are described in-detail in Reid et al. (2022). Amphistomus primonactus was known to be restricted to high elevation in the area east of 152° longitude and between 29.70° and 31.85° latitude. We therefore focussed on sampling from as many burnt and unburnt forested sites in that area as possible, within constraints of time and access. Each site was photographed along its trapline, and notes made on the habitat and fire damage, if any. Coordinates and elevation were recorded at the point nearest to the site that had a GPS signal. Trapping was undertaken in February-March 2021, approximately 14 months after the fires.

The traps were cylindrical pots placed flush with the ground, with a cover to keep out rain. Traps were baited with either fresh kangaroo faeces (three traps per site) or crushed button mushrooms (one trap per site). Each trap contained 2–3 cm depth of ethylene glycol. They were placed approximately 4 m apart in a line at each site. Traps were removed approximately 24 h after placement and their contents sieved directly into 100% ethanol. Amphistomus primonactus was collected in 27 of 119 trap sites (Appendix 1). All specimens used in this study are accessioned and held at the Australian Museum (Appendix 2).

For population genetics analysis, we selected 13 sites (Fig. 1) from the 2021 survey that had approximately 30 individuals each to ensure adequate sample sizes; seven sites were burnt in 2019/20, two were partially burnt (i.e., one of the pitfall traps in the array was in a burnt area, with the others in unburnt areas), and four were unburnt. We selected five sites from the 1993 survey which were in close proximity to our 2021 sites; four of these sites had 30 individuals each, and one site had 29 individuals.

Fig. 1
figure 1

A) Map of Amphistomus primonactus sample localities in NE NSW used for the population genetics analyses, overlayed with FESM burn intensity data. B) Increased-scale map of A. primonactus localities used in the reduced dataset, overlayed with FESM burn intensity data. Pairwise comparisons between localities used in IBD, FST, G"ST, I, and pairwise burn intensities are represented by the black lines. Map data sources: 1 second SRTM Derived DEM-S (Geosciences Australia 2011); Fire Extent and Severity Mapping (FESM) 2019/20 (State Government of NSW and Department of Planning, Industry and Environment, 2020)

Full size image

Amphistomus primonactus was described by Matthews (1974), based on material from just two sites, and Matthews provided the only published key to species of this genus. All species in Amphistomus follow a morphological species concept based on male secondary and primary sexual characters, especially the shape of the parameres. Material collected during fieldwork was compared to these descriptions, and additionally was compared to the type material of the species (held at the Australian National Insect Collection, CSIRO, Canberra). Examination of male genitalia was achieved by removal of the abdomen from the specimen, followed by the extraction of the genital capsule from within the abdomen. All identifications were made by Chris Reid.

We used one leg from each specimen as the source of genetic material for DNA sequencing. These were each placed in an individual well of 96-well PCR plates with 100% ethanol and sent to Diversity Arrays Technology, Inc (Canberra, DArT) for single-nucleotide polymorphism (SNP) sequencing via DArTseq 1.0. DArT 2-row format data were imported into RStudio 2022.02.1 + 461 running R 4.1.2 (Core Team 2021) using dartR 1.9.9.1 (Gruber et al. 2018) and converted to an adegenet 2.1.5 (Jombart and Ahmed 2011) genlight object. The following filters were applied to the dataset using dartR, and are listed in the order that they were used: gl.filter.secondaries, gl.filter.monomorphs, gl.filter.reproducibility, gl.filter.callrate. The secondaries filter examines amplicons (strings of amplified DNA up to 64 base pairs in length) for multiple SNPs based on the DArT Clone ID (a unique value assigned by DArT to each amplicon for identification); within each amplicon, the SNP with the highest reproducibility (proportion of replicate assay pairs for which the SNP call is consistent) and highest average polymorphism information content (avgPIC) was retained while the other SNPs were discarded to reduce the effect of linked genes. The monomorph filter removed any loci that were monomorphic (i.e., only a single nucleotide) across all loci. The threshold for the reproducibility filter was set to 100%. The call rate (the proportion of samples for which the genotype call was a 1 or 0 and not missing data) filter was set to 90%, except for one instance where the threshold was set to 50% (see Results for explanation). The FIS of each sample locality was calculated using gl.basic.stats in dartR to determine whether the filtered dataset exhibited any significant departures from HWE. The finalised dataset was then exported to GenAlEx format for further analysis.

The number of clusters in the data with minimum departure from Hardy-Weinberg Equilibrium (HWE; K) was determined using STRUCTURE 2.3.4 (Pritchard et al. 2000; Falush et al. 2003, 2007; Hubisz et al. 2009). We ran two sets of STRUCTURE analyses: one for the full set of data with less stringent callrate filtering, and one for a reduced dataset of seven sample localities and stringent filtering (see Results for explanation). Each K between 1 and 20 was calculated for the full dataset, while K was calculated from 1 to 12 for the reduced dataset. These calculations were made to assist in determining the optimal K represented in each dataset. Results were analysed with CLUMPAK (Kopelman et al. 2015) using the main pipeline. This pipeline uses two methods to calculate the best K from the STRUCTURE runs: those of Pritchard et al. (2000) and Evanno et al. (ΔK, 2005). Graphical representations of optimal K clusters were drawn by CLUMPAK using distruct (Rosenberg 2004).

We obtained fire severity data from the Fire Extent and Severity Mapping (FESM) 2019/20 dataset (State Government of NSW and Department of Planning, Industry and Environment, 2020) (Fig. 1). Within-locality fire intensities were determined by placing 50 m, 100 m, and 200 m buffers around each sample locality in ArcMap 10.8, then calculating the mean fire severity within each buffer polygon using the Zonal Statistics tool in QGIS 3.26. We determined fire severity between locality pairs by drawing a line between each pair of localities (Fig. 1B), then placing 50 m, 100 m, and 200 m buffers around these lines and calculating the mean fire severity within each buffer polygon using Zonal Statistics.

We calculated the following measures of within-locality genetic diversity using GenAlEx 6.51b2 (Peakall and Smouse 2006, 2012): observed heterozygosity (Ho, Halliburton 2004), unbiased expected heterozygosity (uHe, Halliburton 2004), Shannon’s Information (1 H, Sherwin et al. 2017), and the inbreeding coefficient (FIS, Halliburton 2004). We ran regressions in Microsoft Excel to determine whether within-locality fire severity means were correlated to the within-locality genetic variation measures Ho, uHe, 1 H, and FIS.

Genetic differentiation measures calculated between localities were FST (Halliburton 2004), G"ST (Meirmans and Hedrick 2011) and Shannon’s Mutual Information (I, Sherwin et al. 2017) using GenAlEx. We ran a Mantel test (Mantel 1967) using vegan 2.5-7 to determine whether the results of G"ST and I significantly correspond and therefore provide the same information.

To determine the effect of fire on between-locality genetic differentiation, we needed to account for the effects of other factors on genetic differentiation, namely geographic distance. We ran an isolation-by-distance (IBD) test (Slatkin 1993) using dartR to determine whether geographic and linearised genetic distances correspond significantly. To correct for any effect of geographic distance between localities, we used partial Mantel tests (Legendre and Legendre 2012) for each combination of genetic diversity measure, burn intensity, and geographic distance. The partial Mantel test examines the effect of fire on genetic diversity while removing the effect of geographic distance. We aim to use the linear geographic distance as measured along the surface of the earth. The geographic distance between points was extracted from the IBD calculation in R, which transforms spherical coordinates in decimal degrees to a log-linearised distance. We reversed this conversion by exponentiating the values, then multiplying by a scale factor of 0.875 to account for the geographic transformation used by adegenet in the conversion from spherical coordinates to a linear distance. We used Bonferroni corrections (p = 0.05, n = number of localities, α = p/(n – 1) for all tests with multiple comparisons.

Results

We submitted 512 specimens of A. primonactus to DArT for sequencing; unfortunately, none of the specimens from the 1993 survey had viable DNA for the sequencing process and were therefore omitted from the remainder of this study. Our initial dataset as received from DArT comprised 363 individuals and 27,524 SNPs with 44.44% missing data; filtering the dataset reduced the number of SNPs to 14 with 8.13% missing data. We deemed this an inadequate number of loci for analysis, and we adjusted the callrate filter from 0.95 to 0.5 to produce a larger dataset with 5,472 SNPs and 35.89% missing data.

An excessive amount of homozygosity was present in the dataset (mean FIS of all sample localities = 0.5018). The results from STRUCTURE indicated values of K = 2 (Evanno, ΔK) and K = 3 (Pritchard) (Fig. 2A,B). For both Ks, intermixing between populations was virtually non-existent (Fig. 3A,B). We subdivided the raw dataset to reduce the effect of the bias against heterozygous loci. This is because the inclusion of multiple allopatric taxa within a single dataset during the sequencing process may generate a bias towards markers that discriminate those taxa, leading to an artificial decrease in heterozygosity (Andrezj Kilian, pers. comm.). We did this by excluding the populations that possessed the less frequent genotype (AMPR_06, AMPR_09 – AMPR_13) and retaining the populations with the most frequent genotype (AMPR_01 – AMPR_05, AMPR_07, AMPR_08). We then reanalysed the resultant dataset after applying stringent filtering (i.e., callrate = 0.9); Selechnik et al. (2020) demonstrated that stringent filtering of reduced representation genomic datasets (such as DArT) was important in detection of populations that were slightly differentiated (as in our new, smaller dataset), but stringent filtering was not required for detection of populations that were highly differentiated (as in our original dataset).

Fig. 2
figure 2

Determination of the optimum number of populations (K) represented within each A. primonactus dataset. (A) Estimate of K following Pritchard et al. (2000) for all 13 sample localities. (B) Estimate of K following Evanno et al. (2005) for all 13 sample localities. (C) Estimate of K following Pritchard et al. (2000) for the reduced dataset. (D) Estimate of K following Evanno et al. (2005) for the reduced dataset

Full size image
Fig. 3
figure 3

DISTRUCT graphs for each K from Fig. 2; these graphs report the population membership probability for each individual (represented by vertical bars within each graph). (A) K = 2 for all 13 populations. (B) K = 3 for all 13 populations. (C) K = 3 for the reduced dataset. (D) K = 7 for the reduced dataset. (E) K = 8 for the reduced dataset. (A) and (B) show strong genetic isolation between populations, while C), D), and E) show strong admixture between populations

Full size image

The new filtered dataset with reduced number of sample localities comprised 193 individuals with 90 SNPs and 6.11% missing data. Values for FIS indicated less departure from HWE than the prior dataset, with only AMPR_02 (FIS = 0.182, SE = 0.038) and AMPR_05 (FIS = 0.101, SE = 0.036) indicating a large degree of excess homozygosity. Results from STRUCTURE were mixed, with ΔK returning a bimodal result of K = 3 and K = 8 (Fig. 2D) and Pritchard returning K = 7 (Fig. 2C). In each instance of K the populations exhibit a large degree of intermixing (Fig. 3 C,D,E).

The results for the within-locality genetic diversity measures show there was no significant effect of fire intensity on Ho, uHe, 1 H, or FIS at 50 m, 100 m, or 200 m buffer distances (Tables 1 and 2).

Table 1 Within-locality genetic diversity measures. 1 H = Shannon’s Information; Ho = observed heterozygosity; uHe = unbiased expected heterozygosity = (2 N / (2 N-1)) * He; FIS = Inbreeding Coefficient = (He - Ho) / He = 1 - (Ho / He)
Full size table
Table 2 Results of regression analyses between within-locality genetic diversity measures and mean burn intensity within 50 m, 100 m, and 200 m zones around each locality
Full size table

The highest FST and G"ST values were between localities AMPR_03 and AMPR_08 (Table 3), while the lowest FST and G"ST values were between AMPR_04 and AMPR_05 (Table 3). This differs from I, with the highest value between AMPR_02 and AMPR_03 and the lowest between AMPR_04 and AMPR_05 (Table 3).

Table 3 Pairwise genetic differentiation comparisons of localities using FST, G”ST, and Mutual Information (I), plus the geographic distance between each pair of localities in metres. Top triangles are the p-values for significant difference from zero. Differentiation values that are significantly different from zero after Bonferroni correction (p = 0.008) are bolded. See Fig. 1 for locality codes
Full size table

Mantel tests comparing FST to I (r = 0.4756, p = 0.069, repetitions = 999) and G"ST to I (r = 0.5514, p = 0.053, repetitions = 999) were not significant. A Mantel test comparing FST to G"ST was highly significant (r = 0.9362, p < 0.001, repetitions = 999).

The IBD analysis indicated that there was no significant relationship between geographic distance and genetic information (Mantel, r = 0.1675, p = 0.283). There was no significant effect of fire intensity on between-locality genetic differentiation measures FST, G"ST, I, or FST’ at 50 m, 100 m, or 200 m buffer distances (Table 4).

Table 4 Results of partial Mantel tests for effect on between-locality genetic diversity measures due to mean burn intensity within 50 m, 100 m, and 200 m buffer zones along straight-line transects between each locality pair, removing the effect of geographic distance between each locality pair
Full size table

Discussion

We examined the genetic diversity between populations of A. primonactus that were fragmented by the 2019/2020 wildfires to determine the impact of the fires on dispersal and migration of the species. We investigated the impact of differing burn intensities within sites and between sites to determine whether fire intensity or connectivity had an effect on genetic diversity. We also attempted to genotype A. primonactus collected in 1993 to assess temporal population genetic structure.

The STRUCTURE analysis of the first dataset (inclusive of all sample localities) for both K = 2 and K = 3 showed a distinct pattern of allopatry with little or no gene flow between the populations located north of a line running from Coffs Harbour, through Dorrigo, to Ebor (AMPR_01 – AMPR_09) and those to the south of this line (AMPR_10 – AMPR_13).

Close examination of the male genitalia of all populations showed a slight morphological difference between the northern and southern populations; this morphological difference combined with the distinct genetic separation may warrant the description of the southern population as a new species, because the holotype of A. primonactus was collected from Dorrigo National Park and matches the morphology of our northern specimens. Reid et al. (2022) treated the southern population as a separate morphospecies and found that the habitat preferences between the two taxa was different, with the southern populations restricted to rainforest, and A. primonactus sensu stricto exhibiting more generalist behaviour.

The significant correspondence between FST and G"ST was an expected result, given that the latter measure is a derivation of the former. The lack of significant correspondence between FST and I, and G"ST and I, indicates that geographic pattern of common alleles (which are highly influential for FST) was different to the pattern for I (which weights each allele by its proportion in the population).

Despite our expectation that there would be an effect of fire intensity on genetic diversity of the flightless dung beetle A. primonactus, we did not find any effect on either within-locality diversity or between-locality differentiation. We expected to observe an immediate effect on the genetic diversity of the species due to direct mortality from the fire, incompatible matrix habitat caused by the burnt landscape, or a combination of these factors and others as a result of the environmental disturbance caused by the fires. Because we did not detect an immediate effect of the fires on the genetic diversity of A. primonactus, it is our opinion that any effect exerted on the genetic diversity of the species will likely have a long time lag and would necessitate further sampling and analysis.

Reid et al. (2022) noted the broad range of habitats occupied by this species, including both burnt and unburnt sites, suggesting that it is fire tolerant. Amphistomus primonactus buries faeces in a tunnel in the soil, where the larvae feed (Matthews 1974). The depth of such tunnels is undescribed, but soil temperature between 5 and 10 mm soil depth during and following wildfires of all intensities decreased precipitously in Western Australian jarrah forest (Burrows 1999); burrows greater than 15 mm depth would likely contribute to the survival of individual beetles during a fire event. The ability of the species to take shelter and avoid direct mortality from the fire front, combined with an apparent ability to persist in the fire-altered landscape appear to have minimised the impact of the fires on the species’ reproductive and dispersal ability.

We have determined that there are no implications for the conservation status of A. primonactus, barring any long-term effects that we could not detect. This agrees with the conservation assessment by Reid et al. (2022) of 10 Scarabaeinae species (inclusing two Amphistomus) which found little evidence of any negative impact by the 2019/2020 bushfires on Scarabaeinae populations. Further assessments analysing other invertebrate species in similar scenarios (i.e., pre- and post-wildfire) would help illuminate any patterns of impact of wildfires on the population genetics of invertebrates.