Lineage diversity within a widespread endemic Australian skink to better inform conservation in response to regional‐scale disturbance

Abstract Much attention is paid in conservation planning to the concept of a species, to ensure comparability across studies and regions when classifying taxa against criteria of endangerment and setting priorities for action. However, various jurisdictions now allow taxonomic ranks below the level of species and nontaxonomic intraspecific divisions to be factored into conservation planning—subspecies, key populations, evolutionarily significant units, or designatable units. Understanding patterns of genetic diversity and its distribution across the landscape is a key component in the identification of species boundaries and determination of substantial geographic structure within species. A total of 12,532 reliable polymorphic SNP loci were generated from 63 populations (286 individuals) covering the distribution of the Australian eastern three‐lined skink, Bassiana duperreyi, to assess genetic population structure in the form of diagnosable lineages and their distribution across the landscape, with particular reference to the recent catastrophic bushfires of eastern Australia. Five well‐supported diagnosable operational taxonomic units (OTUs) existed within B. duperreyi. Low levels of divergence of B. duperreyi between mainland Australia and Tasmania (no fixed allelic differences) support the notion of episodic exchange of alleles across Bass Strait (ca 60 m, 25 Kya) during periods of low sea level during the Upper Pleistocene rather than the much longer period of isolation (1.7 My) indicated by earlier studies using mitochondrial sequence variation. Our study provides foundational work for the detailed taxonomic re‐evaluation of this species complex and the need for biodiversity assessment to include an examination of cryptic species and/or cryptic diversity below the level of species. Such information on lineage diversity within species and its distribution in the context of disturbance at a regional scale can be factored into conservation planning regardless of whether a decision is made to formally diagnose new species taxonomically and nomenclaturally.

taxonomic ranks below the level of species and nontaxonomic intraspecific divisions to be factored into conservation planning-subspecies, key populations, evolutionarily significant units, or designatable units. Understanding patterns of genetic diversity and its distribution across the landscape is a key component in the identification of species boundaries and determination of substantial geographic structure within species. A total of 12,532 reliable polymorphic SNP loci were generated from 63 populations (286 individuals) covering the distribution of the Australian eastern three-lined skink, Bassiana duperreyi, to assess genetic population structure in the form of diagnosable lineages and their distribution across the landscape, with particular reference to the recent catastrophic bushfires of eastern Australia. Five well-supported diagnosable operational taxonomic units (OTUs) existed within B. duperreyi. Low levels of divergence of B. duperreyi between mainland Australia and Tasmania (no fixed allelic differences) support the notion of episodic exchange of alleles across Bass Strait (ca 60 m, 25 Kya) during periods of low sea level during the Upper Pleistocene rather than the much longer period of isolation (1.7 My) indicated by earlier studies using mitochondrial sequence variation. Our study provides foundational work for the detailed taxonomic re-evaluation of this species complex and the need for biodiversity assessment to include an examination of cryptic species and/or cryptic diversity below the level of species. Such information on lineage diversity within species and its distribution in the context of disturbance at a regional scale can be factored into conservation planning regardless of whether a decision is made to formally diagnose new species taxonomically and nomenclaturally.

| INTRODUC TI ON
Species are often regarded as fundamental units of conservation concern and correct species delimitation as essential for an unbiased evaluation of biodiversity in a region or a country (Bickford et al., 2007). As a consequence, much attention is paid to species concepts, to ensure comparability across studies and regions when classifying taxa against criteria of endangerment and setting priorities for action. Lack of comparability lies in part in differing opinions on what should be regarded as species and what should be regarded as lineages within species (Sukumaran & Knowles, 2017). Some authors consider all substantive lineages to be species and that speciation occurs at the point lineages first diverge (Fujita et al., 2012); others follow the biological species concept (Butlin & Stankowski, 2020) and its criterion of reproductive isolation (Mayr, 1942). There are many opinion between these extremes (de Queiroz, 1998), points of view that are often held for operational reasons. In particular, reproductive isolation in allopatry can rarely be definitively demonstrated in context and so presents a fundamental operational challenge to the application of the biological species concept. Regarding lineages as species brings a welcome level of objectivity, and so is favored by some citing operational imperatives, but this approach admits potentially numerous ephemeral entities as species, that is, those lineages that will rapidly admix should the allopatric populations again come into contact Kautt et al., 2016;Matute et al., 2020;Momigliano et al., 2017). In this paper, we take the view most recently espoused by Sukumaran and Knowles (2017) that there is a distinction between species in the biological sense and lineage diversity within species.
Were it not for the focus on species-level taxonomy by legislators, biodiversity assessment could entirely sidestep the issue of species delimitation (Miller et al., 2018). Lineage diversity, whether among or within species, and taking into account phylogenetic weighting (Faith & Baker, 2006), is sufficient as a framework for assessing biodiversity and comparing it through time and across regions. Biodiversity assessment and comparability across geographic and temporal scales has become more prominent as human-induced impacts extend from the local to regional, continental, and global scales (Pecl et al., 2017). Various jurisdictions now allow entities below the level of species to be factored into conservation planning, as subspecies, key populations, evolutionarily significant units, or designatable units (Green, 2005), and policy development continues in this area (Hoban et al., 2020), but the definitions of these entities are often vague, confusing for managers, and often ignored when it comes to implementation of priorities (Coates et al., 2018;Mouquet et al., 2012;Vane-Wright et al., 1991).
An excellent case in point of the importance of complete and up-to-date taxonomic assessments in conservation planning is that of the regional catastrophic fires in Australia. In the Austral summer of 2019-2020, approximately 97,000 km 2 of vegetation across southern and eastern Australia was scorched by fires of intensity unprecedented in modern times (Godfree et al., 2021;Ward et al., 2020). Over a billion animals perished in the fires (WWF, 2020), and many species, already endangered, were brought closer to the brink of extinction. But what of cryptic diversity, the diversity represented by lineages within currently accepted species? Species-level taxonomy lags behind the demands of data used to set priorities for conservation (the "taxonomic impediment"); diversity of lineages within species is arguably even more poorly documented for many species.
Yet many of these lineages within species are deeply divergent and some can be regarded as incipient or undescribed species. They are an essential component of biodiversity and their distribution in the context of regional scale disturbance is an important component of conservation planning. Because their ranges are typically narrower than the species to which they are currently assigned, such substantive lineages are particularly vulnerable to regional-scale catastrophic events, such as widespread bushfires.
Geographical barriers, past paleoclimatic incidents, climatic/ environmental factors, and populations surviving glacial maxima in disconnected refugia have played a vital role in the diversification of Australian reptile fauna Chapple, Hoskin, et al., 2011;Dubey & Shine, 2010;Pepper et al., 2018;Rosauer et al., 2015). The family Scincidae represents unsurpassed diversity among vertebrates within the Australian continent (Cogger, 2014;Mitchell et al., 2019). In particular, the forests of eastern Australia are a global biodiversity hotspot and in the top 2.5% of global species richness for lizards (Roll et al., 2017;Williams et al., 2011). Advances in genetic techniques and evolutionary models have accelerated biological studies into isolation and divergence of populations. Current progress in the fields of genomics has opened new opportunities to better understand relationships among the populations within species and their phylogeography.
The eastern three-lined skink, Bassiana duperreyi (Gray, 1838), is a medium-sized, oviparous scincid lizard, distributed broadly across southern and south-eastern Australia (Cogger, 2014), including the areas most impacted by the recent fires (Godfree et al., 2021). The widespread nature of this species across a range of bioregions, including cool Alpine, woodland, and heaths to coastal habitat,

T A X O N O M Y C L A S S I F I C A T I O N
Population genetics provides an excellent model to examine lineage diversity within a single species and relate this to habitat alteration and loss. We generated genome-wide data using the complexity reduction method DArTseq™ (DArT Pty Ltd, Canberra, Australia) which combines nextgeneration sequencing to generate a genome-wide sample of singlenucleotide polymorphisms (SNPs) (Kilian et al., 2012). This technique has recently become a popular tool for understanding genetic diversity, gene flow, lineage phylogenies, species delimitation, and evolutionary history of a range of organisms for which there is little or no prior genomic information Melville et al., 2017).
Our specific goal was to assess lineage diversity within B. duperreyi in south-eastern Australia and make an informed decision on which lineages should be regarded as species and which should be regarded as representing substantial diversity within species. We couple our nuclear DNA data with previously published mtDNA data (Dubey & Shine, 2010) to identify substantial lineage diversity within this species of relevance to management, which will provide a foundation for a comprehensive taxonomic evaluation in the future.
We also consider the implications for management in the context of widespread fires of lineage diversification within this species.

| Sampling
We  Table   S1). Up to 10 samples per locality were collected when available; for samples obtained from museum collections and field collections, we assigned individuals captured within a 20-km radius to a single locality. We conducted fieldwork in areas where the lizards were most  Table S1.

| DNA extraction and SNP genotyping
Tissue samples were sent to Diversity Array Technology Pty Ltd. and SphI (5'-GCATG|C-3') was selected for the complexity reduction by double digestion.
DNA was digested then processed following Kilian et al. (2012), with two different adaptors annealed to the two restriction enzyme overhangs. The library was then subjected to competitive PCR and sequenced using an Illumina Hiseq2500. The sequencing (singleend) was run for 77 cycles. A full account of the DArTseq™ process used to generate sequences for these samples is given by Georges et al. (2018). The data comprised a matrix of SNP loci by individuals, with the contents stored as 0, homozygote, reference state; 1, heterozygote; and 2, homozygote, alternate state.

| Additional SNP filtering
The SNP data and associated locus metadata were read into a genlight object (R Package adegenet; Jombart, 2008) to facilitate processing with package dartR v.1.5.533 . Only loci with >99% repeatability (repAvg) were chosen for subsequent analysis. Further filtering was undertaken based on the call rate (>95%), and secondary SNPs where multiple SNPs occurred within a single sequence tag, only one was retained at random. We did not filter for minor allele frequency because this would potentially inflate the counts of fixed allelic differences, particularly given the sample sizes per locality of only n ≤ 11. The number of samples per locality was small (n ≤ 11), so we could not filter loci for departures from Hardy-Weinberg equilibrium (HWE) or linkage disequilibrium; the sparse sampling of loci across the genome allows the reasonable assumption of little or no linkage between loci. We regard the data remaining after filtering as highly reliable.

| Visualization and qualitative analysis
Genetic similarities for individuals and populations were visualized using principal coordinates analysis (PCoA) as implemented in the gl.pcoa and gl.pcoa.plot functions of dartR. A scree plot of eigenvalues (Cattell, 1966), taken in the context of the average percentage variation explained by the original variables (using the diagnostics provided by gl.pcoa function in dartR), guided the number of informative axes to examine.

| Genetic diversity
Observed heterozygosity was used as a measure of relative genetic diversity. Heterozygosity was obtained for each population from allele frequencies using the gl.report.heterozygosity function of dartR.

| Fixed difference analysis
To examine the possibility that more than one taxon (Operational Taxonomic Unit, OTU) might exist within the geographic distribution of B. duperreyi sensu lato, a fixed-difference analysis was done using the scripts gl.fixed.diff and gl.collapse in dartR. An OTU is defined here as an aggregation of populations that can be differentiated by other such aggregations by one or more diagnostic characters. A fixed difference between two populations at a locus occurs when the populations share no alleles at that locus. Accumulation of fixed differences between the two populations is a strong indication of a lack of gene flow. Fixed differences were summed over populations taken pairwise, and when two populations had no fixed differences, they were combined, and the process repeated until there was no further reduction. The resultant OTUs are, by definition, putatively diagnosable at one or more SNP loci. The set of putative F I G U R E 1 Location of B. duperreyi populations SNP genotyped in this study from across the range of the species in south-eastern Australia and including the location of recognized biogeographic barriers. Color scheme is consistent with other figures and OTUs as described in Figure 2. Underlying map generated using ArcGIS 10.5.1 (http://www.esri.com) and data from the Digital Elevation Model (Geoscience Australia) made available under Creative Commons Attribution 3.0 Australia (https://creat iveco mmons.org/licen ses/by/3.0/au/ legal code, last accessed 9-Jul-20) OTUs arising from the above fixed-difference analysis was then tested for significance (using the test = TRUE option in gl.collapse. recursive.r in dartR, , and pairs for which the number of fixed differences was not statistically significant (number less than expected false positives given the sample sizes) were further amalgamated.

| Phylogeny
To evaluate the relationships among individuals of B. duperreyi, we performed SVDquartets analysis on the individuals grouped by locality (including the individuals from locality: P60 and P59) after further filtration. We selected SVDquartets analysis (Chifman & Kubatko, 2014) because our dataset was composed of short reads with single variable sites per locus (Chou et al., 2015). Heterozygous SNP positions were represented in the dataset by standard ambiguity codes. We used the implementation of SVDquartets in PAUP* v.4.0a165 (Swofford, 2002) with parameters evalQuartets = random, bootstrap = standard, nreps = 10,000 and ambigs = distribute, and designated the B. platynota and B. trilineatus as the outgroup.
We then assessed the relationships concerning their geographic origin and compared the clades identified with their geographic distributions.    Table 1).

| Fixed difference analysis
These diagnosable OTUs are in broad agreement with the structure evident in the PCoA plot ( Figure 2). The fewest fixed differences were observed between Kangaroo Island OTU and Mt Lofty Ranges-Fleurieu Peninsula OTU, whereas the most were observed between South-eastern highland and Australian alps OTU and the Kangaroo Island OTU.

| Phylogenetic inference
The fixed difference analysis directs consideration of diagnosability on to the phylogeny to identify those lineages that are diagnostic

| DISCUSS ION
This study provides quantification of population structure for the eastern three-lined skink B. duperreyi and delineates the distribution of structural elements across the landscape to provide valuable TA B L E 1 Genetic diversity of Bassiana duperreyi. Matrix of Euclidean genetic distances (above diagonal) and fixed genetic differences (below diagonal) between the final set of operational taxonomic units to arise from a fixed difference analysis applied to the ingroup data set. Comparisons were based on an average of 12,532 loci after filtering for call rate >95%. All fixed differences were significant at p < .0001 guidance to conservation efforts to protect linage level species diversity in the face of regional scale divergence. Accurate species delimitation of species and structural elements within species are critical for evaluating biodiversity on regional or global scales and for informing effective conservation strategies (Mace, 2004). Failure to appropriately identify species boundaries, structure within species, and the distribution of these across the landscape can have significant consequences for conservation management of imperilled species by diverting finite conservation resources (Agapow et al., 2004).

South eastern highland and Australian alps
Here, we studied B. duperreyi as a case study for how within-species diversity and its distribution across the landscape in the context of regional-scale disturbance can substantially alter perspectives and priorities set for conservation.
The eastern three-lined skink B. duperreyi is currently regarded as a single species (Hutchinson et al., 1990). Mitochondrial variation within this species (Dubey & Shine, 2010)  is, sister to the remaining lineages to the west (Figure 3). Our largest OTU from the fixed difference analysis combines the 16 populations from the south-eastern Alpine region (see Table 1) (Schäuble & Moritz, 2001;Symula et al., 2008), marsupial dunnarts (Cooper et al., 2000), and grasshoppers (Kawakami et al., 2009). It is a concordant pattern that is believed to be the result of Pleistocene. This agrees with the sea-level-driven dispersal opportunities identified by Dubey and Shine (2010).
The level of structure we observed across the range of B. duperreyi is comparable to other widely distributed species within a geographical range of similar size (Chapple, Hoskin, et al., 2011;Smissen et al., 2013;Symula et al., 2008). This presumably arises as a response to common ecological or geomorphological barriers-such as the Great Dividing Range, the Murray River or habitat structure across the landscape and, as we have discussed, episodic isolation from the mainland of Tasmania, Flinders Island, and Kangaroo Island by rising sea levels. The Murray River appears to be a barrier to dispersal for a range of species (reviewed by Ansari et al., 2019), including lizards, since it extensively dries and breaks into discrete pools nearly once every century (Close, 1990). Bassiana duperreyi fits this pattern, and dating studies for this (Dubey & Shine, 2010) and other species suggest a Plio-Pleistocene diversification in response to barriers afforded by Lake Bungunnia which formed ca 2.5 Mya and persisted to ca 700 Kya when the modern Murray River was established (McLaren et al., 2011;Stephenson, 1986). It appears that attributes of the modern Murray River and the habitats supported in its basin have been sufficient to maintain signatures of the more ancient divergences in both nuclear and mitochondrial genes of B. duperreyi (present study) and Tiliqua rugosa (Ansari et al., 2019), though the haplotype distribution of the 11 nuclear genes in T. rugosa were less definitive than the SNP markers were for B. duperreyi.
Does B. duperreyi comprise more than one species? We have shown the species to comprise five lineages that have diverged to the point of being diagnosable by one or more corroborated allelic fixed differences. Under some species concepts, diagnosability of a lineage is sufficient to warrant recognition as a species. We take the view, however, that genetic diagnosability, while necessary for a lineage to be considered as a named taxon, is not sufficient . That is, we admit the possibility of substantial is not supported and will hopefully prompt revision supported by a combined genetic and morphological analysis. A formal description of the species is beyond the scope of our current study.

F I G U R E 3
Phylogenetic analyses of Dartseq SNPs with SVDquartets (left) compared to a published phylogeny of two partial mitochondrial genes (ND 2 and ND 4 ) (not to scale) (see Dubey & Shine, 2010 Without adequate morphological data, we cannot resolve the taxonomic issues here. However, we propose that there are two putative species within B. duperreyi. The first is distributed in southeastern highland and Australian alpine region (including Wilsons Promontory, Flinders Island, and Tasmanian) and the second is an aggregation of diagnosable lineages (ESUs) occupying the lower elevation regions and coastal regions (including Kangaroo Island, the type locality) (Figure 3). They have broadly parapatric distributions ( Figure 1). The south-eastern highland and alpine taxon is of particular interest, because it has a system of sex determination that involves both differentiated sex chromosomes and sex reversal of the XX genotype to a male phenotype at a frequency that varies predictably with elevation and nest temperatures (Dissanayake et al., 2020. Our study has clear management implications in the context of regional catastrophic events and progressive habitat fragmentation and modification at regional levels. The recent fires in Australia were restricted to the eastern mainland portion of the range of B. duperreyi as currently defined (and Kangaroo Island) (Figure 4).
Concern for the impact of the fires is ameliorated somewhat by the existence of substantial populations of B. duperreyi to the west that were unimpacted by fire. When we consider biodiversity within the species, and in particular the existence of five diagnosable lineages, two of which could be considered putative species, then concerns are reignited. The mainland habitat of the distinctive eastern highlands and Australian Alps OTU (putative species 1 and potentially unnamed) has been severely impacted by fire, which raises concerns for its conservation. Tasmania and Flinders Island provide insurance, but the widespread impact of fire on the habitat of this eastern highlands and Australian Alps lineage is now of serious management concern. Similarly, our demonstration of the distinctive lineage of B. duperreyi on Kangaroo Island, also seriously impacted by the recent fires, raises its priority for its assessment and management. Our genome-wide SNP analysis of structure within a widespread species highlights the need for biodiversity assessment at the regional scale to include an examination of cryptic diversity below the level F I G U R E 4 The distribution of Bassiana duperreyi in relation to the intensity and extent of the Australian megafire event, which occurred from 1 st July 2019 to 11 th February 2020. Refer to Supporting Information Table S1 and Figures 1 and 3 for the corresponding population details. The fire intensity and distribution data were obtained from Godfree et al. (2021). The name of the each megafire bracketed value after each megafire is the fire area in millions of hectares. Underlying map generated using ArcGIS 10.5.1 (http://www.esri.com) and data from the Digital Elevation Model (Geoscience Australia) made available under Creative Commons Attribution 3.0 Australia (https://creat iveco mmons.org/licen ses/by/3.0/au/legal code, last accessed 9-Jul-20)

Mid-North
Wollemi ( of species and perhaps also a re-evaluation of species delimitation.
Such information on lineage diversity within species and the distribution of those lineages across the landscape can be factored into conservation planning regardless of whether a decision is made to define new species with formal taxonomic decisions.

ACK N OWLED G M ENTS
For their assistance with museum collection access and tissue loans,

CONFLICT OF INTEREST
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
All data not available in the paper or supplementary materials are deposited in the Dryad Digital Repository (https://doi.org/10.5061/ dryad.tx95x 69zc) together with the scripts necessary to replicate the analysis.