Genetic diversity of a flightless dung beetle appears unaffected by wildfire

The wildfires of Australia’s Black Summer in 2019/2020 caused a massive loss of wildlife and habitats, but the effects of the fire on invertebrate species post-burn are unknown. We hypothesised that the fires would negatively affect the genetic diversity of invertebrate species by impeding movement between populations due to habitat loss. We studied the genetic diversity of a flightless dung beetle, Amphistomus primonactus Matthews 1974, to determine the impact of the wildfires on this species. We examined 90 SNPs from 193 individuals across seven localities impacted by the wildfires in north-eastern New South Wales. We used STRUCTURE to determine the overall population structure of the seven localities. We calculated four within-locality genetic diversity measures (observed heterozygosity (Ho), unbiased expected heterozygosity (uHe), Shannon’s Information (1H), and the inbreeding coefficient (FIS). We calculated three between-locality genetic diversity measures (Fixation Index (FST), Hedrick’s G”ST, and Shannon’s Mutual Information (I). We used partial Mantel tests to compare the between-locality genetic diversity measures with the mean fire intensity along each pairwise linear transect, while accounting for genetic variation due to geographic distance. We compared the within-locality genetic diversity measures to the mean fire intensity at each site. STRUCTURE showed a large degree of intermixing between localities. We found no significant effect of fire on any within-locality genetic diversity measure, or on any between-locality genetic diversity measure. We suggest that the genetic diversity of A. primonactus was not significantly affected by the Black Summer wildfires. Implications for insect conservation: Our results show that the 2019/2020 wildfires had a negligible impact on the genetic structure of A. primonactus. This offers a promising outlook for the species in its recovery from the fires.


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 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.
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 F IS 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. 2003Falush et al. , 2007Hubisz 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. 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 F IS 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 nonexistent (Fig. 3A,B). We subdivided the raw dataset to We calculated the following measures of within-locality genetic diversity using GenAlEx 6.51b2 Smouse 2006, 2012) Sherwin et al. 2017), and the inbreeding coefficient (F IS , 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 H o , uH e , 1 H, and F IS .
Genetic differentiation measures calculated between localities were F ST (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 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).
The new filtered dataset with reduced number of sample localities comprised 193 individuals with 90 SNPs and 6.11% missing data. Values for F IS indicated less departure 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 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 F ST , G" ST , I, or F ST ' at 50 m, 100 m, or 200 m buffer distances (Table 4).

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 HWE than the prior dataset, with only AMPR_02 (F IS = 0.182, SE = 0.038) and AMPR_05 (F IS = 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 H o , uH e , 1 H, or F IS at 50 m, 100 m, or 200 m buffer distances (Tables 1 and 2).
The highest F ST and G" ST values were between localities AMPR_03 and AMPR_08 (Table 3), while the lowest F ST 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).
Mantel tests comparing F ST 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 F ST to G" ST was highly significant (r = 0.9362, p < 0.001, repetitions = 999). 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 F ST and G" ST was an expected result, given that the latter measure is a 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 Table 3 Pairwise genetic differentiation comparisons of localities using F ST , 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   derivation of the former. The lack of significant correspondence between F ST and I, and G" ST and I, indicates that geographic pattern of common alleles (which are highly influential for F ST ) 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 longterm 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.