Population genomics of a predatory mammal reveals patterns of decline and impacts of exposure to toxic toads

Abstract Mammal declines across northern Australia are one of the major biodiversity loss events occurring globally. There has been no regional assessment of the implications of these species declines for genomic diversity. To address this, we conducted a species‐wide assessment of genomic diversity in the northern quoll (Dasyurus hallucatus), an Endangered marsupial carnivore. We used next generation sequencing methods to genotype 10,191 single nucleotide polymorphisms (SNPs) in 352 individuals from across a 3220‐km length of the continent, investigating patterns of population genomic structure and diversity, and identifying loci showing signals of putative selection. We found strong heterogeneity in the distribution of genomic diversity across the continent, characterized by (i) biogeographical barriers driving hierarchical population structure through long‐term isolation, and (ii) severe reductions in diversity resulting from population declines, exacerbated by the spread of introduced toxic cane toads (Rhinella marina). These results warn of a large ongoing loss of genomic diversity and associated adaptive capacity as mammals decline across northern Australia. Encouragingly, populations of the northern quoll established on toad‐free islands by translocations appear to have maintained most of the initial genomic diversity after 16 years. By mapping patterns of genomic diversity within and among populations, and investigating these patterns in the context of population declines, we can provide conservation managers with data critical to informed decision‐making. This includes the identification of populations that are candidates for genetic management, the importance of remnant island and insurance/translocated populations for the conservation of genetic diversity, and the characterization of putative evolutionarily significant units.


| INTRODUC TI ON
Declines of vertebrate species are occurring globally due to habitat degradation, invasive species, land clearing and harvesting (Díaz et al., 2019). The Australian continent and surrounding islands harbour highly unique faunas and have been exceptionally vulnerable to such changes, with 33 mammal species becoming extinct since European colonization in 1788 (Roycroft et al., 2021;Woinarski et al., 2019). As species continue to decline, population genomic analysis provides an opportunity to identify historical patterns of gene flow and population structure, as well as investigate and quantify the impacts of decline on genetic diversity. With genetic diversity a critically important component of global biodiversity measures (Hoban et al., 2021), such analyses can benefit threatened species management by providing data on populations that are candidates for genetic intervention, quantifying the importance of remnant island and insurance/translocated populations for conserving genetic diversity, or characterizing putative evolutionarily significant units (Frankham et al., 2017;Ralls et al., 2018;von Takach, Penton, et al., 2021;Wright et al., 2020). Broadly, consideration of these factors can assist in the conservation of adaptive potential of threatened populations, by promoting genetic diversity and adaptive variation (Kelly & Phillips, 2016;Madsen et al., 1999;von Takach Dukai et al., 2019).
The invasion of exotic species into new environments, and the associated novel selection pressures that invasive species impose, can induce a range of adaptive evolutionary responses in native species (Berthon, 2015;Strauss et al., 2006). Importantly, when substantial declines of native populations occur, the resulting loss of standing genomic variation within populations is likely to hinder adaptation to subsequent environmental change (i.e., loss of adaptive potential). Few studies have quantified the impacts of invasive species on population genomic diversity on declining native species, with most relevant studies investigating the genetic consequences of disease outbreaks (Schoville et al., 2011;Serieys et al., 2015).
One species that is well suited for such analysis is the northern quoll (Dasyurus hallucatus). The northern quoll is a small (<1500 g) generalist marsupial predator (Oakwood, 1997) that was, until relatively recently, abundant across a broad distribution and range of ecosystems in northern Australia (Oakwood, 2008), and spans several biogeographical breaks that drive the population genetic structure of many taxa (Bowman et al., 2010;Catullo et al., 2014;Edwards et al., 2017;Eldridge et al., 2012;Rollins et al., 2012). Since European colonization of Australia, from 1788 onwards, the northern quoll has declined across much of its former distribution (Moore et al., 2019), partly due to pernicious and compounding threats such as land clearing, inappropriate fire and grazing regimes, and the impact of feral cats (Braithwaite & Griffiths, 1994). From 1935 onwards, the decline of the northern quoll has been greatly accelerated by the introduction and continuing spread of cane toads (Rhinella marina) in Phillips et al., 2003;Sabath et al., 1981). Like many Australian taxa, the northern quoll has no evolutionary history of exposure to toad (Bufonidae) toxins (cf. rodents; Cabrera-Guzmán et al., 2015), and are therefore extremely sensitive to the toxic defences of cane toads (Phillips et al., 2003;Smith & Phillips, 2006). As with some large Australian predators (e.g., freshwater crocodiles Crocodylus johnstoni; Letnic et al., 2008), northern quolls mistake cane toads for suitable prey and suffer precipitous population-level impacts upon their arrival (Shine, 2010). Consequently, northern quolls have been extirpated from many localities (Burnett, 1997;Moore et al., 2019;Woinarski et al., 2010;Woinarski, Legge, et al., 2011).
The northern quoll is now listed as Endangered under state, national and international listings (IUCN Red List), and populations are being actively managed by various jurisdictions. Because the spread of cane toads has so far proved impossible to halt, current management actions are aimed at mitigating the impact of toads.
For example, the Northern Territory government established insurance populations on two offshore islands in 2003 (Rankmore et al., 2008), and there is ongoing research into the feasibility of increasing the frequency of toad-avoidance traits in populations via training (Indigo et al., 2021;Jolly et al., 2020) or translocation of heritable "toad-smart" genotypes into toad-naïve populations ("targeted gene flow") Kelly & Phillips, 2019).
Critically, these management strategies require understanding and managing various genetic aspects for them to be successfully implemented.
Here, we provide a comprehensive assessment of the population genomic structure of the northern quoll, illustrating how patterns of genomic diversity may be impacted by cane toad invasion and various management strategies, to improve the efficacy of decisionmaking. By assembling one of the most comprehensive sets of samples of any Australian mammal species to date, including toad-naïve, toad-exposed, island and translocated populations, we investigate patterns of genomic diversity and population genomic structure within and among these treatments, as well as demonstrate impacts of the ongoing spread of cane toads. We also identify putatively adaptive loci that have higher-than-expected levels of differentiation (suggestive of diversifying selection) as well as loci that are significantly associated with the presence of cane toads. We predict that our genomic analyses will identify (i) population genomic structure resulting from vicariance across biogeographical barriers; (ii) patterns of drift associated with isolated island populations and (iii) large impacts of cane toads, or other drivers of population decline, on population genomic diversity, with population bottlenecks in toad-exposed populations potentially contributing to a reduction in genetic diversity. Our temporal replication will also enable us to directly observe the effect of establishing insurance populations on offshore islands in the short term (16 years). We suggest that this study will provide a quantitative example of the impacts that invasive species can have on the patterns of genomic diversity in native species, the conclusions of which will benefit attempts to mitigate the impacts of invasive species.

| Sample collection
We collated tissue samples from across the entire species distribution from multiple personnel who had been involved in mammal trapping between 1993 and 2019, covering the regional jurisdictions of Queensland (QLD), the Northern Territory (NT) and Western Australia (WA) (Figure 1). While some localities had all samples collected in a single year, others were sampled during multiple visits across one or two years, with low capture rates making large sample sizes in any one session difficult to obtain. Trapping efforts within a locality were variable among jurisdictions and field teams, but typically involved sampling a set of several small grids (e.g., 1 ha) of cage traps spread over 10-20 km. All samples taken from the set of grids within a locality were given an identifying name for population genomic analysis. In total, 459 samples were collected from 41 localities, with five localities in QLD, 17 localities in the NT and 19 localities in WA ( Figure S1). Four of the localities were trapped on separate occasions, which we have treated independently (Table S1). Of these 46 independent sampling periods/localities, 21 were represented by 10 or more individuals and 30 were represented by six or more individuals. The mean pairwise geographical distance between all sampled localities was 883 km, with a maximum distance between sampled individuals of 3220 km.
For each sampling locality, we determined whether cane toads were present at the date of collection by cross-referencing with published timelines of cane toad spread (Brown et al., 2015;Doody et al., 2019;Phillips et al., 2007;Urban et al., 2008). Localities were identified as belonging to one of four treatments, including (i) those that had never been exposed to cane toads ("toad-naïve"); (ii) those that have been exposed to cane toads ("toad-exposed"); (iii) those naturally occurring on toad-free islands ("island"); and (iv) those established via translocations onto toad-free islands ("translocated").
Islands were considered separately from the toad-naïve treatment because long periods of isolation and small effective population sizes at these locations have generally resulted in low levels of genetic diversity (Cardoso et al., 2009;Spencer et al., 2017). Translocated populations of northern quolls on Astell Island and Pobassoo Island have been present since 2003, when the species was introduced from the NT mainland to the islands as insurance populations (Cardoso et al., 2009). In total, 24 of the 46 independent sampling periods/locations were classed as toad-naïve, 10 were classed as toad-exposed, and eight were from naturally toad-free islands (NT and WA only). One location in the NT, "East Alligator," had >10 samples collected both prior to, and after, the arrival of cane toads to the region, allowing for preliminary before-and-after comparisons at this locality to be made.

| DNA extraction and sequencing
Ear tissue samples (~3 mm 2 ) were prepared and extracted in plate format (Qiagen DNeasy 96 Blood & Tissue Kit) following the standard protocol with an extended lysis (incubation at 56°C for 2 h then reduced to 37°C overnight). Samples consisting of dried blood spots on filter paper were extracted with a slight modification to the lysis step (Choi et al., 2014), where ~25 mm 2 of filter paper containing dried blood was incubated at 37°C overnight in 400 μl of 1× PBS buffer, prior to completion of the Qiagen protocol from the addition of ethanol at step three.
Double-stranded DNA concentrations were quantified using a Qubit 3.0 Fluorometer and normalized to 200 ng DNA in 25 μl, and then arranged on 96-well plates, where each plate contained two within-plate technical replicates, one technical replicate from each of the other plates and one negative control (blank) well.
Plates were sent for double-digest restriction-associated DNA (ddRAD) sequencing at the Australian Genome Research Facility in Melbourne, Victoria. An optimal combination of two restriction enzymes was determined using three establishment samples (broadly representative of the species distribution), with PstI and NlaIII considered most suitable for achieving the best level of amplification and minimizing repetitive sequences. Briefly, the library preparation protocol consisted of (i) digestion using PstI and NlaIII,

| Bioinformatic pipeline and SNP filtering
Raw sequence data were demultiplexed using the process_radtags function of stacks (Catchen et al., 2013), with checking for intact RAD sites and reads truncated to 135 bp to meet a quality threshold of 30. Of the 1.647 billion total reads, 95.8% were retained after demultiplexing. These were then mapped to the draft assembly of the northern quoll reference genome (scaffold N50 89 Mb) (S. The program angsd version 0.93 (Korneliussen et al., 2014) was used to filter reads and call single nucleotide polymorphisms (SNPs) to output an SNP-by-sample matrix. Reads were only used to call SNPs if the mapping quality score was ≥20 (thus excluding reads that mapped poorly or mapped to repeat regions of the genome).
Loci were only retained if they (i) were polymorphic, as inferred by a likelihood ratio test p value ≤1 × 10 −5 (Kim et al., 2011);(ii)  To ensure that relationships between individuals could be accurately inferred from the data, we used these SNPs and samples to construct a hierarchical clustering dendrogram based on genetic distance, with visual examination of the dendrogram confirming that all technical replicates paired closely together on long branches ( Figure S2). We also made a higher-level bootstrapped dendrogram by using genetic distances among sampling localities instead of individuals ( Figure S3). The percentage difference between called genotypes of technical replicates was also used to confirm that genotyping error rates were low (mean 99.99% ± 0.003% SE similarity between replicates). One of each pair was then removed from all further analyses. The data set was then further pruned to remove highly related individuals and SNPs showing signals of putative selection, both of which can bias analyses that were developed using the assumption of unrelated individuals and neutrally evolving loci.
To estimate relatedness between individuals, we used the recently developed method-of-moments estimate for the kinship coefficient (Goudet et al., 2018;Weir & Goudet, 2017), implemented in the "hierfstat" package (Goudet, 2005). This approach estimates kinship values between pairs of individuals relative to the average kinship values of all pairs of individuals in the sampled population (i.e., within each locality). This avoids the problematic requirement for estimation of reference allele frequencies. We found 23 pairs of highly related individuals (kinship coefficients >0.25, i.e., full sib or higher) in a total of six sampling localities, and randomly removed one of each pair, retaining a total of 352 samples for use in the remaining analyses.

| Putatively adaptive loci
We investigated the SNP panel for genetic signatures of local adaptation using two methods: identification of SNPs with higher-thanexpected levels of differentiation (F ST outliers), and identification of SNPs that were significantly associated with the presence of cane toads (via a latent factor mixed model). To identify outlier SNPs that showed signatures of diversifying selection, we first computed an ancestry matrix for each individual using the methods implemented in the LEA package (Frichot & François, 2015) and calculated an F ST statistic based on this ancestry matrix Martins et al., 2016). The F ST statistic was then transformed into squared z-scores, with p-values computed using a chi-squared distribution (Frichot & François, 2015;Weir, 1996). The Benjamini-Hochberg algorithm (Benjamini & Hochberg, 1995), with a false discovery rate of 1 in 1000, was used to correct for issues associated with multiple testing. To identify SNPs associated with the presence of cane toads, we first imputed all missing genotype data, via the "impute" function of the LEA package (method = "mode"), based on the ancestry matrix computed above. As with many genotype-environment association methods, latent factor mixed models require no missing data. The estimated length of time (number of years) that each population had been exposed to cane toads was used as the environmental predictor variable with which genotypes were associated. We generated latent factor mixed models (via the "lfmm2" exact least-squares function of the LEA package) to identify allele frequencies that were correlated with cane toad occurrence (Caye et al., 2019;Frichot & François, 2015). This method controls for population structure via a number of latent factors equal to the number of ancestral populations (von Takach, Ahrens, et al., 2021). The p-values for each SNP were adjusted using the in-built robust estimate of the genomic inflation factor to improve the shape of the histogram  and subjected to a Benjamini-Hochberg algorithmic correction (Benjamini & Hochberg, 1995) to ensure a low rate of false discovery (1 in 1000).
We recorded the genomic locations of all SNPs under putative selection, investigated the extent of overlap in results between the two methods, and then removed all SNPs (485 in total) from the data set that were identified using either method. This produced a selectively neutral data set for population genomic analyses, many of which rely on assumptions of neutrality. While future work may attempt to further elucidate and understand the patterns and genomic architecture of selection in response to cane toad invasion, the set of candidate loci identified here (Table S4) provides a starting point for understanding adaptation across the northern quoll genome.
After all filtering steps had been conducted, we retained a total of 10,191 SNPs and 352 individuals for the remaining analyses. The overall level of missing data for the final filtered SNP-by-individual matrix was 1.2%.

| Population genomic diversity and structure
Using the filtered SNP data set, we calculated mean values of stand- (P E ). The East Alligator site was separated into toad-naïve (2003) and toad-exposed (2014 and 2017) subsets prior to calculations being made. Similarly, the historical (2005) samples from Astell Island were separated from the more recent (2016) samples. Calculations were made using the "gstudio" package (Dyer, 2016) with the standard small sample size correction used to account for bias in estimates of expected heterozygosity resulting from differences in sample size between localities. As heterozygosity calculations based on filtered SNP data can potentially show biases due to sample sizes and filtering parameters, we also calculated observed and expected values of autosomal heterozygosity for each locality using the methods of Schmidt et al. (2021), which considers both monomorphic and polymorphic nucleotides. This included building aligned sequences into a stacks catalogue via the "ref_map" pipeline, with the filtered and sorted BAM files as inputs (using only previously filtered individuals), and analysing the data set using the core program "Populations" with all missing sites removed. The heterozygosity estimates and standard errors in the subsection of the summary output titled "# All positions (variant and fixed)" were recorded, and Pearson's product moment correlation coefficient was used to determine the strength of the relationship between values of autosomal heterozygosity and SNP heterozygosity.
We visualized geographical patterns of population genomic structure using a principal coordinate plot of pairwise genetic distances between individuals. The first two principal coordinate dimensions of the genetic distance matrix were calculated using the cmdscale function, and individuals plotted on the coordinates. To investigate patterns of within-region population structure, we also performed independent principal coordinate analyses on subsets of individuals that formed the major clusters in the comprehensive plot.
To assign individuals to ancestral population genomic clusters, investigate patterns of admixture between populations and assess hierarchical population structure, we used the alternating projected least squares algorithm implemented in the "tess3r" package (Caye et al., , 2018. This method applies a model of genetic structure featuring a discrete number (k) of ancestral populations, allowing for independent investigation of values for k that have low cross-entropy metrics (Frichot et al., 2014;Frichot & François, 2015). It also incorporates the spatial location of sampling, to remove bias associated with patterns of isolation-by-distance. Cross-entropy criteria were package (Pembleton et al., 2013), which calculates significance values by bootstrapping across loci ("nboot = 1000"). The G ′′ ST values were calculated using the pairwise_Gst_Hedrick function of the "mmod" package (Winter, 2012). The G ′′ ST metric is an F ST analogue that uses bias-corrected estimates of within-population gene diversity (H S ) and total gene diversity (H T ) to improve inference of demographic history among populations (Hedrick, 2005;Meirmans & Hedrick, 2011). To test the level of correlation between the F ST and G ′′ ST metrics, we conducted a Mantel test on the pairwise matrices, using 10,000 permutations.
To visualize patterns of isolation-by-distance we constructed a scatterplot of standardized F ST values (Rousset, 1997) against geographical distance for all pairwise comparisons of mainland localities (where n ≥ 6). A Mantel test was used to test the significance of geographical distance on standardized genomic differentiation, using 10,000 permutations. Isolation-by-distance among all mainland individuals was also investigated using spatial autocorrelation of multilocus genotypes, plotted on a correlogram. Pairwise geographical distances between individuals were calculated from latitude and longitude observations for each individual via the earth.dist function of the "fossil" package (Vavrek, 2011). Correlation and significance values were made using the genetic_autocorrelation function of the "gstudio" package (Dyer, 2016) with 999 permutations, and the results plotted in R.

| Demographic and temporal trends
To test the relative influence of vicariance and cane toads on genomic diversity among our natural treatments (toad-naïve, toadexposed and island localities), we constructed four generalized least square models using the "nlme" package (Pinheiro et al., 2020). The response variables were A E , F IS , H E and P E , with treatment as the predictor variable, and a Gaussian spatial correlation structure included to account for spatial autocorrelation of predictor values within geographical regions. We note that toad-exposed populations cluster towards the eastern side of the species distribution, due to the invasion front spreading from the far east. Tukey multiple comparison tests were then used to identify whether mean values for each response were significantly different among treatments. To test the effect that isolation on islands (i.e., Astell Island) or toad invasion and establishment (i.e., East Alligator) had on genomic diversity metrics, we constructed individual generalized linear models for each metric, where loci were samples and the predictor variable was a factor corresponding to the year of the survey.
To identify patterns of historical demography at a locality, including the potential influence of cane toads on population size, we investigated past changes in effective population size (N e ) using three methods, of which one is appropriate for estimating recent (i.e., tends to hundreds of generations) trends and two that are most useful for deeper time (i.e., thousands to hundreds of thousands of generations) trends.
To estimate recent historical trends in effective population size (N e ), we used the version 1.1 software tool (Barbato et al., 2015), which analyses SNP data based on the relationship between linkage disequilibrium (LD) and N e , with corrections for sample size and recombination rate. A variant call format (VCF) file output from angsd was split into individual VCF files for each sampling locality, with plink version 1.9 (Purcell et al., 2007) used to generate input files for snep. We used Sved and Feldman's (1973) mutation rate modifier for correcting the recombination rate, a sample size correction for unphased genotypes, and a recombination rate of 9.1 × 10 −9 (Feigin et al., 2018). Effective population sizes for the past 200 generations (or to the limit calculated by snep) were plotted and the x-axis scaled to time assuming a mean generation time of 2 years (Pacifici et al., 2013). To increase the available data for the toad-exposed QLD localities, we combined the localities of Hope Vale and Black Mountain, as these were geographically close and clustered together in the principal coordinates plot. For historical context, we added the timing of the first major European colonization event on the Australian continent (1788). For populations that had been exposed to cane toads at the time of sampling, we also added the timing of cane toad invasion to the locality.
To investigate deeper-time demographic history, we used two independent techniques: one (the "stairway plot") that estimates demographic histories based on the composite likelihood of a population's folded site frequency spectrum (Liu & Fu, 2015, 2020, and one (the "extended Bayesian skyline plot" [EBSP]) that uses coalescent theory to trace multilocus lineages backwards through time (Heled & Drummond, 2008).
To build stairway plots for each locality (where n ≥ 8), we constructed a folded site frequency spectrum using genotype likelihoods estimated by the angsd software package. All filtering criteria based on mapping and read qualities were the same as previously described, with the filters for polymorphic loci removed to capture information on all monomorphic and polymorphic nucleic sites. By using likelihoods, the site frequency spectrum could be built from between 121 million and 565 million nucleic sites per locality. The site frequency spectrum was loaded into the folded blueprint file provided with the software, with an appropriate mutation rate per site per generation (Feigin et al., 2018), and 200 input files (cf. bootstraps) used to produce the median trend and 95% confidence intervals. The results for each locality were read into R and all plotted on identical scales, with the approximate timing of the last glacial maximum (19,000 years before present) added to each panel for historical context.
We constructed extended Bayesian skyline plots for a subset of localities, to compare results with the stairway plot method and assess their feasibility for use with ddRAD data. We followed the methods of Trucchi et al. (2014), selecting a small number of variable RAD loci (via a custom R script) each 135 bp long, and converting these to NEXUS format for analysis with beast (Suchard et al., 2018).
We were able to identify 10 loci for all localities except Weipa, for which only nine loci were found. To minimize bias due to linkage between RAD loci, the selected loci were typically located on separate scaffolds of the reference genome (Table S2), with the minimum distance between any pair of selected loci >11 million base pairs. As in Trucchi et al. (2014), the number of parameters in the model was limited to achieve convergence of the chains. One site model, HKY with four Gamma categories, was defined for each locus class, containing 2-7 SNPs (Table S3). As no information was available to calibrate the substitution rate, we plotted unscaled reconstructions.
The Markov chain Monte Carlo (MCMC) analyses, with lengths of 15-40 billion iterations, were run in duplicate for each population to check for adequate convergence. Convergence was checked using tracer  and consistency in estimates of the posterior distribution of the effective sample size. The posterior distributions of the number of population size changes ("sum(indicators. alltrees)" posterior mean) and 95% highest posterior density interval were also recorded to identify whether the constant population size hypothesis could be rejected for each locality.

| Putatively adaptive loci
Together, our outlier and latent factor mixed model analyses identi-  (Table S4), with 19 SNPs physically located <100 bp from another candidate SNP.

| Population genomic diversity and structure
The principal coordinates plot identified three predominant clusters (with some substructure, Figure 2 Figure  S4). The largest drop in cross-entropy occurred between k = 1 and k = 2, followed by the drop between k = 2 and k = 3, suggesting that a value of 2 or 3 is most useful for describing the high-level genetic structure across the species. At k = 2, the WA populations separated from the combined NT and QLD populations (Figure 3).
Temporally replicated localities showed some significant changes through time, corresponding with isolation on islands or toad invasion and establishment (Table 1) (Table S4). While we were unable to scale the temporal axis for the EBSPs, the WA populations showed weaker declines than suggested by snep, suggesting that, as for the stairway plots, the deeper-time analysis using ESBPs may not be observing very recent declines (e.g., last 50-100 years). This was most pronounced for populations in the highly remote and rugged areas of the western Kimberley, where threats may have developed more recently.

| DISCUSS ION
While biodiversity declines are well recognized as a global problem, the consequences of these declines for patterns of genomic diversity are often poorly known. Using genomic data to infer recent and long-term historical demographic trends can help identify the extent to which loss of genomic diversity itself is threatening species' persistence. Northern Australia's vast savannas are the location of some of the world's most dramatic recent declines of mammals  Woinarski, Legge, et al., 2011), and here we illustrate patterns of genomic diversity and long-term trajectories for a keystone declining species, the northern quoll. Two major influences appear to be driving the observed patterns of northern quoll population genomic structure and diversity across the Australian continent. First, major known biogeographical barriers are responsible for strong hierarchical population structure resulting from patterns of vicariance and long-term isolation. Second, it seems that the effective population size of many populations has been declining over much of the past century, in some cases exacerbated by the rapid range expansion of the introduced and lethally toxic cane toad.
While ongoing longitudinal analysis of northern quoll populations will help disentangle the impacts of historical biogeography from the impacts of cane toad invasion on population genomic diversity, it appears as though multidecadal exposure to the cane toad is associated with a considerable reduction in heterozygosity. This is approximately equivalent to the loss of diversity on islands that have been isolated from the mainland for >10,000 years. Together, these results warn that ongoing and future declines of mammal populations in northern Australia resulting from anthropogenic pressures and introduced species may have long-term impacts on species-wide genomic diversity, despite remnant populations persisting in some localities.
This data set represents the most comprehensive tissue sampling of a mammal across northern Australia to date and is thus an important quantification of the genomic impacts of regional mammal declines. Our data suggest that effective population sizes have been declining even in the absence of cane toads, suggesting that cane toads represent a compounding threat for northern quoll populations. While many areas are lacking quantitative historical data on population sizes, our results agree with a small number of studies and anecdotal data suggesting that northern quoll populations were indeed declining prior to the introduction and spread of cane toads (Braithwaite & Griffiths, 1994;Woinarski et al., 2008). Such declines have also been occurring in a range of other native mammal species across northern Australia Penton et al., 2021;von Takach et al., 2020;Woinarski et al., 2010). This massive conservation problem appears to be driven by a suite of interacting factors, including habitat degradation due to feral herbivores and pastoralism (Legge et al., 2011), increased predation due to feral cats (Felis catus) , frequent or high-intensity fires (von Takach et al., 2020, and possibly climate change (Kutt et al., 2009;Traill et al., 2011).
These broad drivers of decline contribute to isolation of remnant populations and reduction of gene flow, both of which can lead to reductions in effective population size. Unfortunately, the introduction and spread of cane toads appears to have exacerbated these threats and is causing a major loss of diversity as cane toads spread across the Australian mainland. This loss of diversity probably results from the combination of strong selection on genotypes expressing avoidance behaviours (Kelly & Phillips, 2017) as well as the severe reductions in northern quoll population sizes when cane toads arrive in an area (Burnett, 1997;Woinarski et al., 2010).
In this context, future monitoring and tissue sampling of north- While we did not observe significant increases in the inbreeding coefficient (F IS ) in toad-exposed populations, it is worth noting that this metric does not provide information on the cumulative level of F I G U R E 4 Effect of three treatments on five northern quoll population genomic parameters. Parameters include the effective number of alleles (A E ), inbreeding coefficient (F IS ), expected heterozygosity (H E ), observed heterozygosity (H O ) and the locus polymorphic index (P E ). Treatments include populations that naturally occur on toad-free islands (Island), occur on the mainland and have not been exposed to cane toads (Naïve), or occur on the mainland but have been exposed to cane toads (Exposed). inbreeding in a population, but rather detects nonrandom mating in the most recent generation (Waples, 2015). With further longitudinal sampling of impacted populations, investigation of alternative metrics (e.g., mean kinship coefficients within a population) may be warranted. Our Kimberley population genomic data provide a comprehensive baseline understanding of toad-naïve populations in the region, and can be used in a before-and-after-impact analysis framework. Ongoing sampling in the region will clarify the immediate impacts of cane toad invasion, but sampling over the next two to three decades will be vital for understanding the temporal pattern of the loss of genomic diversity.
The implication that mammal declines may have been occurring for long periods of time (prior to the start of widespread wildlife monitoring efforts) could help quantify the relative importance of recent (<50 years) changes in landscape management as drivers of decline. For example, while feral cats have probably been present for ~150 years , the impacts of climate change have probably been occurring for much longer periods (Brüniche-Olsen et al., 2014;White et al., 2018). Using genomic data to obtain information on historical effective population sizes is one valuable way to disentangle these drivers of demographic change. With genome-wide sequencing data becoming cheaper and more readily accessible, we suggest that an important avenue of future work will be to conduct similar analyses on additional species with similar distributions, allowing us to compare and contrast patterns of decline between species.
As anthropogenic stressors impact ecological communities, declining species tend to show contractions in their realized niche space (Scheele et al., 2017). Populations within a subset of environmental conditions can thus exhibit stronger patterns of persistence . For many Australian mammal species, the drivers of decline appear to be either mitigated or better tolerated on offshore islands (Legge et al., 2018). For example, while cane toads have shown uncontrolled spread across the mainland, they have not yet managed to establish viable populations on numerous islands.
Some of these islands serve as useful conservation refuges for populations of northern quolls (How et al., 2009), including two islands onto which northern quolls have been translocated as a means of conservation management (Rankmore et al., 2008).
Our data show that islands certainly have a strong role to play in the conservation of genomic diversity, with WA and the NT home to remnant island populations of northern quolls that contain a portion of the species-wide diversity. For many island lineages, this probably includes a small unique contribution of alleles to overall allelic richness (von Takach, Penton, et al., 2021). However, as is the case for many species and islands around the world, these populations tend to show reduced genetic diversity (Cardoso et al., 2009;Frankham, 1997;Spencer et al., 2017  is somewhat limited, and there may be a case for translocations between islands, and from the mainland, if these islands are to act as refugia for substantial portions of the available genetic diversity.
Islands are, however, also not immune to invasion by cane toads.
A sobering example is the extirpation of northern quolls on the Sir Edward Pellew Island group, after cane toads established populations by rafting on freshwater plumes and debris associated with a large flooding event on the McArthur River (Woinarski, Ward, et al., 2011).
While remnant island populations exhibit lower genetic diversity, our data also show that the process of translocating northern quolls to toad-free islands as insurance populations has not led to a substantial reduction in genome-wide diversity metrics. While changes in some metrics, such as effective population size, can take over 20 generations to be detectable in RADseq data (Nunziata & Weisrock, 2018), islands such as Astell Island have carrying capacities in the thousands of individuals (Griffiths et al., 2017), making drift a weaker force and allowing for the retention of high levels of genomic diversity. Such insurance populations require ongoing management and genetic monitoring, however, and even short-term isolation of island populations from predators can result in the rapid loss of critical behavioural traits (Jolly et al., 2018;. This may prevent successful reintroduction of island individuals back to the mainland into their former realized niche. Further, successful establishment into historically occupied niche space may be made more difficult due to the ecological release of interspecific competitors or the invasion of exotic species into vacated niche space (Doherty & Driscoll, 2018), factors that require consideration for future management actions.
Broadly, our data concur with previously observed phylogeographical studies, showing patterns of isolation resulting from biogeographical barriers such as the Great Sandy Desert, the Bonaparte Gulf/Ord Arid Region and the Carpentarian Gap. These barriers have variously been observed to limit dispersal in plants (Edwards et al., 2017), frogs (Catullo et al., 2014), reptiles (Melville et al., 2011), birds (Lee & Edwards, 2008) and mammals (Eldridge et al., 2014;Potter et al., 2012). Early collections of northern quolls, led Thomas (1926) to recognize four distinct forms, including Dasyurus hallucatus hallucatus, D. h. nesaeus, D. h. exilis and D. h. predator, with type localities from the Top End (NT), Groote Eylandt (NT), eastern Kimberley (WA) and Cape York (QLD), respectively (Mahoney & Ride, 1988). These subspecies have not been recognized by subsequent authors due to the apparent lack of distinct morphological differentiation (Jackson & Groves, 2015;Oakwood, 2008;Viacava et al., 2020), despite molecular evidence of two to four divergent lineages (Firestone, 2000;Woolley et al., 2015). Here, we present data showing strong hierarchical population structuring of these divergent lineages, which conforms to the disjunct nature of the northern quoll distribution across the continent. Ongoing morphological assessment will help to clarify differences in phenotype between F I G U R E 6 Stairway plots of estimated historical effective population size (N e ) at 14 sampling localities of the northern quoll. The solid black line in each panel represents the median estimate of N e , with shaded regions representing 95% confidence intervals. Colours match those used in Figure 3 (k = 5)  these lineages (Umbrello, 2018). Another direction currently being investigated is to use exon capture data to investigate in greater detail the level and timing of divergence between the major lineages presented here (Eldridge et al., 2020). Further investigation of the functional and adaptive significance of putatively adaptive loci will complement such findings, with future genome annotation probably able to provide compelling evidence for agents of selection, physiological processes, phenotypes and the occurrence of selective sweeps (Pardo-Diaz et al., 2015).
Together, our data suggest that several critical aspects need consideration for the management of northern quoll populations across the Australian continent. These include (i) that biogeographical barriers separate subregions each containing substantial amounts of unique genomic diversity; (ii) that remnant island populations harbour only a subset of the diversity present within each subregion; (iii) that large insurance populations on islands or in intensively managed/fenced safe havens can maintain high levels of genomic diversity (although these populations may lose adaptive traits and require ongoing genetic management) and (iv) that declines induced by cane toads are dramatic, but play out against broader declines resulting from other anthropogenic stressors and threats. With the exception of cane toad impacts, these considerations are likely to be similar for many of the widespread faunal species within the Australian monsoonal tropics (Bowman et al., 2010).
The dramatic impact of toads on population sizes of northern quolls points to the potential value of using offshore islands as insurance populations, and also to the value of preventing cane toads from colonizing the Pilbara region of Western Australia (Southwell et al., 2017). Our work also raises questions about whether there may be value in mixing populations across islands or between mainland and islands to secure the genetic diversity of quolls. Future work on quolls and other taxa will help to clarify the extent to which intraspecies conservation management needs to operate within and among lineages or evolutionarily significant units across northern Australia.

AUTH O R CO NTR I B UTI O N S
SB, BvT and BLP designed the research with input from KO and DF.
CPB, SFC, RH, EK, BLP, IJR, PBSS and GJT contributed to data collection and/or sample preparation. BvT and LR carried out the data analysis. BvT and SB interpreted the results with input from all authors. BvT wrote the paper with input, advice and contributions from all authors with respect to manuscript structure, framing and intellectual content.

ACK N OWLED G EM ENTS
We acknowledge the Traditional Custodians of the lands on which the research presented here was undertaken, and pay our respects to their Elders past, present and emerging. We also acknowledge the contribution of the Oz Mammals Genomics Initiative consortium (https://ozmam malsg enomi cs.com/conso rtium/) in the generation of data used in this publication. Librarians.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest. construction.

B EN EFIT-S H A R I N G S TATEM ENT
A research collaboration was developed with researchers around Australia, who provided tissue samples for genetic analysis. All collaborators were offered the opportunity to contribute to the manuscript as co-authors, and the results of research are available for the broader scientific community through online databases and open access publishing. The research addresses a priority concern, in this case the conservation genetics and management of the threatened northern quoll, and findings will be incorporated into management actions for the species. Bioplatforms Australia is a nonprofit organization that supports Australian Life science research by investing in state-of-the-art infrastructure and expertise in genomics, proteomics, metabolomics and bioinformat- Stephen Frankenberg https://orcid.org/0000-0002-2497-7356