Long- term biogeographical processes dominate patterns of genetic diversity in a wingless grasshopper despite substantial recent habitat fragmentation

Low- vagility species may hold strong genetic signatures of past biogeographical pro - cesses but are also vulnerable to habitat loss. Flightless grasshoppers of the morabine group were once widespread in southeastern Australia, including Tasmania, but are becoming restricted to remnant patches of vegetation, with local ranges impacted by agriculture and development as well as management. Habitat fragmentation can generate genetically differentiated “island” populations with low genetic variation. However, following revegetation, populations could be re- established, and gene flow increased. Here we characterize single nucleotide polymorphism- based genetic varia - tion in a widespread chromosomal race of the morabine Vandiemenella viatica (race 19) to investigate the genetic health of remnant populations and to provide guidelines for restoration efforts. We update the distribution of this race to new sites in Victoria and Tasmania, and show that V . viatica populations from northern Tasmania and eastern Victoria have reduced genetic variation compared to other mainland populations. In contrast, there was no effect of habitat fragment size on genetic variation. Tasmanian V . viatica populations fell into two groups, one connected genetically to eastern Victoria and the other connected to southwestern Victoria. Mainland populations showed isolation by distance. These patterns are consistent with expectations from past biogeographical processes rather than local recent population fragmentation and emphasize the importance of small local reserves in preserving genetic variation. The study highlights how genomic analyses can combine information on genetic variability and population structure to identify biogeographical patterns within a species, which in turn can inform decisions on potential source populations for translocations.

increased. Here we characterize single nucleotide polymorphism-based genetic variation in a widespread chromosomal race of the morabine Vandiemenella viatica (race 19) to investigate the genetic health of remnant populations and to provide guidelines for restoration efforts. We update the distribution of this race to new sites in Victoria and Tasmania, and show that V. viatica populations from northern Tasmania and eastern Victoria have reduced genetic variation compared to other mainland populations. In contrast, there was no effect of habitat fragment size on genetic variation. Tasmanian V. viatica populations fell into two groups, one connected genetically to eastern Victoria and the other connected to southwestern Victoria. Mainland populations showed isolation by distance. These patterns are consistent with expectations from past biogeographical processes rather than local recent population fragmentation and emphasize the importance of small local reserves in preserving genetic variation. The study highlights how genomic analyses can combine information on genetic variability and population structure to identify biogeographical patterns within a species, which in turn can inform decisions on potential source populations for translocations.

K E Y W O R D S
biodiversity, biogeography, invertebrate, population structure, SNP markers

| INTRODUC TI ON
Many insects are thought to be declining around the world (Black & Vaughan, 2009;Hafernik, 1992;Leidner & Neel, 2011;Sands, 2018), with particularly well-documented cases of reductions in flying insects (Hallmann et al., 2017). The causes of this decline are multifactorial and include potential factors such as pesticide applications and climate change, but an overriding factor is habitat destruction as agricultural production intensifies, and other factors such as residential developments and infrastructure projects impinge on natural habitats (Harvey et al., 2022;Sánchez-Bayo & Wyckhuys, 2019).
This decline in turn threatens the provision of ecosystem services, including pollination, food provision and soil nutrient turnover, provided by insects.
Although much of the attention on insect declines has so far been on flying species that rely on natural habitat for components of their life cycle, dramatic effects are also likely in low-vagility insects that are resident in natural habitats, such as flightless herbivores and ground-dwelling insects that feed on detritus. Declines in these species are expected to be related directly to the loss of habitat. However, much less information is available on such species, given that they often tend to be cryptic (Baur et al., 2020) and poorly known taxonomically (New & Samways, 2014). Low-vagility invertebrates also have potential to reveal details of past biogeographical events due to the strong genetic signature they provide on population fragmentation and expansion events (Hugall et al., 2002;Sunnucks et al., 2006). However, genetic effects of habitat fragmentation may obscure such signals (Myburgh & Daniels, 2015).
Flightless matchstick grasshoppers (family Morabidae) represent one unique group of low-vagility insects which have an unusually long history of study. The group is unique to Australia and has diversified across a wide range of habitats (Key, 1977). In southeastern Australia, the decline in members of this species and particularly the species Keyacris scurra was already noted in the 1950s (White, 1956).
As a consequence of agricultural land clearing, K. scurra rapidly became restricted to small "islands" of habitat such as maintained in cemeteries where the native vegetation associated with this species (native daisies, native grasses and particularly Themeda triandra) persisted (Hoffmann, White, et al., 2021;Rowell & Crawford, 1995;White, 1956). Sporadic surveys highlighted the continued decline of this species (Hoffmann, White, et al., 2021;Rowell & Crawford, 1995) which now has a national listing as being "endangered" and of conservation concern.
A few other related matchstick grasshoppers occur in southeastern Australia and are likely to be equally threatened by native habitat destruction. These include other undescribed species in the genus Keyacris as well as the closely related genus Vandiemenella (morabine genera are described in Key, 1977). While most of these have been poorly studied, Vandiemenella viatica has received attention in evolutionary studies because it is composed of a series of geographically separate chromosomal races that inspired initial ideas around modes of speciation (Kawakami, 2008;Kawakami et al., 2011;White et al., 1964). White et al. (1967) noted that "… they only survive, somewhat precariously in many cases, as small, isolated colonies separated by large areas of wheat and pasturelands." This species has disappeared from many of its former habitats because of land clearing, weed invasion and other developments, although it remains quite abundant in local pockets of woodland, grassland and heath (Kawakami, 2008).
Although low-vagility insect species are easily damaged by habitat destruction, they also have the potential to persist in small remnant habitat patches that remain suitable. For instance, surveys of K. scurra have uncovered their presence in small patches of tens of square metres that have persisted within modified landscapes (Hoffmann, White, et al., 2021). When habitat requirements are known, this also makes such species suitable for translocation efforts for revegetated areas even when these are quite small (as is typically the case).
To inform management of low-vagility species and their future translocations, it is important to undertake genetic studies that can infer patterns of past connectivity among populations, evidence of inbreeding in remaining populations and future adaptive potential based on the persistence of genetic variation generally (Weeks et al., 2016;Willi et al., 2006). So far, population genetic studies have been undertaken on a few morabines, including recent population genomic studies on some Warramaba species (Kearney et al., 2022) and K. scurra (Hoffmann, White, et al., 2021). The latter showed that many K. scurra populations retain a surprisingly high level of genetic variation and show a lack of inbreeding even when found in small patches of habitat.
In this study, we have undertaken a detailed study of one of the chromosomal races of V. vandiemenella, namely the "viatica 19" race found in eastern parts of South Australia, as well as in Victoria and Tasmania (Kawakami et al., 2011;White et al., 1964). Kawakami et al. (2009) found that there was quite deep structure across populations of this chromosomal race when assessed through allozyme markers, particularly between Kangaroo Island and other populations from South Australia and Victoria/Tasmania ( Figure 1). In contrast, a different picture emerged from COI variation, where the viatica 19 race included two genetically distinct groups despite having a similar cytotype.
Here we focus on the eastern viatica 19 group of populations. This group covers half the range of the species in an area that has been particularly negatively affected by landscape modifications in the past 200 years (Lunt, 1991). We provide a detailed population genetic assessment of this group with a focus on the impacts of recent glacial climatic cycle effects in the region, including the formation and loss of the Bassian land bridge across Bass Strait connecting Victoria and Tasmania most recently at 25-18 thousand years ago when woodland first developed and was then replaced by grassland (Adeleye et al., 2021). This study also helps inform population translocations taking place with this group around Victoria to reintroduce these insects into managed areas (Hiromi Yagui, unpublished).
We use single nucleotide polymorphism (SNP) markers and mitochondrial DNA (mtDNA) to examine patterns of genetic variation within and between the viatica 19 populations, and we include samples from newly discovered fragments that extend the distribution of this race substantially. We consider the following issues. (1) Is there isolation by distance in this chromosomal race, or is isolation connected to biogeographical factors, such as the formation of the Bassian land bridge? (2) Is there evidence for inbreeding which might associate with the size of remaining remnants? (3) Is there an association between current habitat size and genetic variation as measured by SNP-based heterozygosity? (4) Is there an association between climate and genetic diversity? Answers to these questions provide a picture of landscape processes affecting this species in southeastern Australia and inform ongoing translocations involving revegetated areas.

| Surveys
Grasshoppers were collected across 2019-2021 with aspirators from an area of >20 m 2 per site and preserved individually in 100% ethanol in Eppendorf tubes before being stored at −20°C. We targeted areas where the chromosomal race had previously been found in woodland, coastal heath and grassland habitats, based on records in the literature (Kawakami, 2008;White et al., 1964) Figure 1 with details of samples provided in Table 1. Because we extended the known distribution of Vandiemenella viatica to the west in Tasmania to Narawntapu National Park, we checked the karyotype of this population following White et al. (1967) to confirm that the Narawntapu population has the standard "viatica 19" karyotype of two large pairs of autocentric autosomes, one large pair of metacentric autosomes and one metacentric (or submetacentric) × chromosome, and a considerable number of smaller autosomes to a total of 19 chromosomes.

| DArT-Seqprocessing
A high-throughput genotyping method using the DArT-Seq technology at Diversity Arrays Pty Ltd was used. Here, complexity reduction is used to enrich nuclear genome representations with active genes and low-copy sequences through combinations of restriction enzymes and reduction methods (https://www.diver sitya rrays.com/ techn ology -and-resou rces/darts eq/ and Kilian et al., 2012). Implicit fragment size selection and next-generation sequencing of representatives were subsequently performed with HiSeq2000 (Illumina) (Georges et al., 2018;Kilian et al., 2012). This technology was con-   Grasshopper hind limb tissue (upper half) was supplied to Diversity Arrays Pty Ltd for high-density (~2.5 million sequences per sample used in marker calling) DArTSeq assays. Eight samples were first tested with multiple restriction enzyme combinations, and an "optimal" set was determined based on the fraction of the genome represented, while controlling average read depth and the number of polymorphic loci (https://www.diver sitya rrays.com/ techn ology -and-resou rces/darts eq/). DArTSeq DNA extraction, sequencing and SNP genotyping methods are explained in detail elsewhere (Georges et al., 2018;Kilian et al., 2012). DArT uses a hybridization approach to focus on low-copy sequences rather than repetitive sequences (Kilian et al., 2012, pp. 69-71) which is an important consideration in grasshoppers with typically large genomes.

| Bioinformatics
Raw reads from N = 171 samples were trimmed using the program cutadapt version 2.8 (Martin, 2011) to remove adaptors and barcodes from the beginning of the sequence (16-18 bp) (−b command), and then bases from the end were trimmed to create a final maximum length of 63 bp (−l 63). This ensures that the same sequences were compared downstream despite variable-length barcodes. Barcodes were those given in the metadata supplied by DArTSeq in the "barcode" column, prepended to the connector sequence "TGCAG." Trimmed reads were processed with the denovo_map.pl pipeline of stacks version 2.55 (Catchen et al., 2013). This pipeline as- The 2,300,939 variant sites were initially filtered to retain one SNP per locus with a minor allele count (MAC) = 1 at a missingness threshold of no more than 15% missing individuals per site via vcftools (-mac 1, -max-missing 0.85, -thin 5000) (this thinning works here as loci are never longer than 500 bp in this data set). Using the resulting data set of 833 SNPs, missingness data for each individual were generated (vcftools -missing-indv) with vcftools version 0.1.13 (Danecek et al., 2011). Following this step, lower quality duplicates were removed from the data set and remaining unique individuals with more than 40% of sites missing were also removed. Most of these corresponded to individuals which had failed DArTSeq quality checks.
The final data set was constructed by varying individuals per population, MAC (in the case of heterozygosity), whether all SNPs or one SNP per sequence were included in the analysis, and whether all nucleotides including monomorphic nucleotides were included in estimates of heterozygosity. MAC rather than minor allele frequency was used, as recently advocated for population structure analysis (Linck & Battey, 2019). We used a different filtering approach when considering estimates of population variation (Schmidt et al., 2021) vs. population structure. An MAC of 1 was used when characterizing variability within populations. V. viatica individuals (n = 150) were filtered to include biallelic SNPs present in >95% of individuals, one per RAD tag, and an MAC of 1. A total of 1752 SNPs were retained across the 150 individuals. The mean sequence depth per site across individuals and the overall site depth spectrum are provided in Figure S1. The median site coverage per individual is around 30; only 11 have a coverage of more than 60 per individual, and only one site has a value of more than 90. Coverage is quite widely spread with no strong right bias.
For multivariate analyses, missing data were interpolated by mean, first within sites, then across the entire data set.

| Populationgenomicanalysis
Observed heterozygosities (H O ) and inbreeding coefficients (F IS ) were calculated for each site with the stacks version 2.55 program "populations" (Catchen et al., 2013;Rochette et al., 2019). For each population, all nonmissing polymorphic and monomorphic sites were used for these calculations rather than only polymorphic sites.
This use of nonmissing sites for H O calculation follows recent recommendations (Schmidt et al., 2021).
Discriminant analyses of principal components (DAPCs) were performed for each filtered not available-interpolated SNP data set with the "adegenet" R package (Jombart, 2008). Clusters were assigned using the "find.clusters" function. The function "xvalDapc" was then used to identify the lowest number of PCs required to generate a relatively low root mean square error (RMSE) for assignment.
Finally, this number with the relevant clustering was used to generate the final DAPC with the "dapc" function. An analysis of molecular variance (amova) was carried out by the R package "ad4" to identify the proportion of variance associated with regions identified in the PCA and DAPC (see below) in comparison to the variation among collection sites within regions.

| Isolationbydistance
For all mainland samples, isolation by distance (IBD) was tested with a Mantel test with 10,000 permutations using the function "mantel.randtest" from the "ade4" R package (Dray & Dufour, 2007). For the Mantel tests, Euclidean genetic distances between individuals were used. Geographical distances were calculated using the "pointDist" function from the "enmSdm" R package (Smith, 2022), which was set to use the "distGeo" function from the "geosphere" R package (Hijmans et al., 2019). IBD analysis was repeated on population-level rather than individual-level data.

| mtDNAandanalysis
We screened 53 individuals from 16 populations for a 700-bp section of the mitochondrial cytochrome oxidase subunit 1 (COI) gene using primers and conditions described in Hoffmann, White, et al. (2021). Mitochondrial haplotype analyses were conducted with the R package "pegas" (Paradis, 2010) following initial trimming and alignments with the program mega11 (Tamura et al., 2021). We included sequences from Victorian and Tasmanian collections from Kawakami et al. (2007) in the mtDNA trees for comparison.

| Vegetationmapping,andassociations
We looked at the possible relationship between genetic diversity and habitat size across many geographical scales (Hoffmann, White, et al., 2021). This link was expected because We also compared H O against the latitude, longitude and annual mean temperature at each site. Temperature was extracted from the Worldclim data set (Fick & Hijmans, 2017), and we expected that that cooler areas would have lower H O if the species' southern distribution limit is related to temperature.

| Geneticstructureofpopulations
The PCA identified six or seven main population groups when considered along the two major axes (Figure 2), with South Australia  Figure S4). This was also evident in the structure analysis with k = 10 ( Figure S5; from the k = 15 data).
An amova indicated that 69.4% of the molecular variance was accounted for by the regions identified above. This compares with a lower value of 10.8% for variation between the sample sites within a region and 19.8% for variation among individuals within the sample sites. Variation among regions and sites within regions was highly significant (p < .001).
A Mantel test indicated a significant association between geographical and genetic distance (r = .5487, p < 1e-5, 10,000 permutations). IBD was evident in comparisons between individuals from the different collection sites but clear groupings formed in these comparisons ( Figure 3). These are presented separately with and without comparisons for the Tasmanian individuals, which formed separate groups closely related to the different Victorian sites, as shown in Figure 2. The top graph in Figure 3 shows  showing low levels of genetic variation, although some confidence intervals for these populations was quite wide.

| mtDNAvariation
The mtDNA variation did not always segregate well with nuclear variation; while western district locations (Kiata, Little Desert) clustered together, other western populations (Hindmarsh, Junction Dam) were separated from these and connected to them through Truganina which was from a different SNP group (Figure 6a,b). Sandringham from the Melbourne group and Kangaroo Ground from the northern Melbourne group were also not separated based on mtDNA data. The Melbourne and western Victoria groups were particularly diverse with respect to mtDNA variation, and the two Tasmanian population groups identified via SNPs were also separated from each other. There was similarity between the mtDNA variants from Kawakami et al. (2007) and our sequences collected from nearby locations (e.g., Grampians and Portland with Hindmarsh) but also separation as in Tasmanian

| Vegetationandclimateversus heterozygosity/inbreeding
We found no significant association between H O and the size of the patches from which grasshoppers were sourced at different buffer distances (Figure 7; Figure S7). There was a tendency for heterozygosity to decrease in larger areas but this was driven by the Tasmanian samples which tended to have low heterozygosity but come from the more vegetated sites that encompassed larger coastal heath reserves. Inbreeding levels (F IS ) were also not connected significantly to areas of vegetation. There was a strong positive correlation (r = .81) between H O and annual mean temperature (Figure 8), which was similar in strength to the correlation with latitude (r = .80) although the correlation with longitude was even stronger (r = .94).

F I G U R E 4
Autosomal heterozygosity in samples arranged by different regions from west to east but with Tasmania separate. The Tasmanian populations were all from areas with an annual temperature below 13.5°C and showed no correlation with temperature on their own (p = .251), while the correlation for mainland populations considered by themselves was r = .59 (p = .001).

| DISCUSS ION
Through additional sampling, we have expanded the known range of Vandiemenella viatica in an easterly and northerly direction in Victoria and a westerly direction in Tasmania. We provide the first record of this species in box-ironbark forest at Chiltern, a vegetation type that extends from central Victoria just across the border into south-central New South Wales (http://www.virid ans.com/ECOVE G/box-ironb ark.htm). We also extended the range of the species eastwards along the Victorian coastline into Gippsland, particularly around coastal heath. In Tasmania, we identified several new areas of coastal heath where the species occurs (Table 1), which includes heath around Musselroe Bay and to the west at Narawntapu.
Attempts to locate the species from other sites between these areas separated by more than 120 km were unsuccessful, with most of the coastal areas now cleared for agriculture and development. We also F I G U R E 6 Haplotype network for mtDNA variation in individuals tested from different sites identified to samples (a) or groups defined in Figure 2 as defined by SNP variation (b). See Table 1 for sites, not all of which were tested for mtDNA variation. We have also included samples from Kawakami et al. (2007) from Victoria and Tasmania.
failed to detect the species in heath west of Narawntapu such as at Rocky Cape, so Narawntapu may represent the northwestern limit of the species in northern Tasmania.
Despite the persistence of V. viatica in many sites, in contrast to Keyacris scurra which is now only found at very few sites in Victoria, V. viatica is now restricted to habitat fragments of suitable where an extensive area of heath had been burnt a year previously and failed to detect V. viatica even though the species persisted in heathland in an adjacent unburnt area. Inappropriate prescribed burns are a major threat to insect conservation in Australia but have been poorly studied (Sands, 2018).
Despite fragmentation of the species, we found that there was little impact at this stage of fragment size on levels of local genetic variation or inbreeding. In this sense, our findings with V. viatica are akin to those we obtained with K. scurra where a similar analysis also failed to find any association with fragment size. In the case of V. viatica, we note the persistence of genetic variation and lack of inbreeding in small areas such as Truganina and Diamond Creek (both  (Willi et al., 2006). The population size of both K. scurra (Hoffmann, White, et al., 2021) and V. viatica may also be quite large even in relatively small isolated fragments, with up to two or three individuals per m 2 recorded in some of our surveyed areas. We explored these issues further through One interpretation of this pattern is that the proto-viatica race originated in South Australia where there is a diversity of chromosomal races (Kawakami et al., 2009;White et al., 1964) and then moved eastwards, with bottlenecks occurring during this colonization process that has resulted in lower levels of variation towards the east. Given the blunted nuclear differentiation across the chromosomal races reflecting ongoing gene flow among the races (Kawakami, 2008;Kawakami et al., 2007Kawakami et al., , 2009 The low level of heterozygosity in Tasmanian populations but distinctness of the Narawntapu population from the other east coast populations points to two colonization events from the mainland into Tasmania. Based on the SNP genetic patterns it appears Gippsland was the source in one event and a region southwest of Melbourne in the other event. Separate colonization events into Tasmania are well established in plants (Adeleye et al., 2021) and have been suggested in a few invertebrate species, but data are more limited. In particular, phylogenetic data involving mtDNA and one nuclear gene suggest that northwestern populations of the soil millipede Pogonosternum nigrovirgatum may have been colonized from central Victoria, whereas northeastern Tasmania was colonized by an undefined source from eastern Victoria and probably from Gippsland (Decker, 2016). Several local invertebrate faunal breaks have also been recognized in Tasmania that include separation of the western and eastern zones in northern Tasmania (Mesibov, 1994) and could contribute to differences among populations. It also appears that the Yarra River in the Melbourne area ( Figure 1) has acted as a barrier between lineages that may have recolonized the region from different refugia after the Last Glacial Maximum. We found differences in mtDNA variation that linked partly to the geographical location of the samples, such as the western Victorian samples which had diverged from other samples and to some extent from each other, perhaps because of past gene flow with populations from other chromosomal races in South Australia.
In other cases, there was a poor link between the mtDNA and SNP information. This included the high diversity present in the samples from the Melbourne region linked to several other groupings.
Obviously, there is much more power in using many SNP markers across the genome when establishing population structure, although mtDNA sequences remain important in distinguishing species and assessing patterns deeper in the phylogeny, and mtDNA data have been used successfully in other morabine grasshopper studies (Kawakami et al., 2011;Kearney & Blacket, 2008). The mtDNA data do highlight the substantial separation between V. viatica 19 samples from different areas; population separation was also evident in K. scurra particularly with respect to isolation of one region from a montane region, although variability among COI haplotypes was lower in this species (compare Figure 6 with Figure 4 in Hoffmann, White, et al., 2021).
The information we have collected here is crucial to ongoing translocations aimed at re-establishing the species in revegetated and managed sites across the wider Melbourne region and beyond. With increasing interest from councils, attempts to restore natural ecosystems are starting to turn from canopy vegetation and ground cover to the fauna that can be supported by the vegetation (Lindenmayer et al., 2012). V. viatica and other grasshoppers can provide important ecosystem services in such areas and are of inherent interest in their own right given the uniqueness of the morabine group. When translocating grasshoppers, it is important to consider both the level of genetic variation and inbreeding in populations acting as sources and the genetic distance between them and populations that might have persisted in the revegetated target sites. The groupings we have identified can assist in this process, by using individuals in translocations that have a high level of heterozygosity but that are nevertheless related to these populations genetically (Hoffmann, Miller, & Weeks, 2021). For instance, if the isolated population at Narawntapu became extinct due to bushfires or other events, it would make sense to source translocation material from genetically heterogeneous populations from the central Victorian group rather than from eastern Tasmania.
In conclusion, the SNP data we have produced provide a clear

ACK N O WLE D G E M ENTS
We thank Katie Robinson for assistance with the molecular work, Amy Liu and Sarah Jaboor with assistance with sampling, and Ann Stocker for checking the cytology of the Narawntapu population.

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors declare that there are no conflicts of interest associated with this work.