Uncovering inbreeding, small populations, and strong genetic isolation in an Australian threatened frog, Litoria littlejohni

The status of many amphibian populations remains unclear due to undetected declines driven by disease and difficulties in obtaining accurate population estimates. Here, we used genome complexity reduction-based sequencing technology to study the poorly understood Littlejohn’s treefrog, Litoria littlejohni across its fragmented distribution in eastern Australia. We detected five identifiable genetic clusters, with moderate to strong genetic isolation. At a regional scale, population isolation was likely driven by population crashes, resulting in small populations impacted by founder effects. Moderate genetic isolation was detected among populations on the Woronora Plateau despite short distances between population clusters. Evidence of recent declines was apparent in three populations that had very small effective population size, reduced genetic diversity and high inbreeding values. The rates of inbreeding detected in these populations combined with their small size leave these populations at elevated risk of extinction. The Cordeaux Cluster was identified as the most robust population as it was the largest and most genetically diverse. This study exemplifies the value of employing genetic methods to study rare, cryptic species. Despite low recapture rates using traditional capture-recapture demographic methods, we were able to derive population estimates, describe patterns of gene flow, and demonstrate the need for urgent conservation management.


Introduction
Amphibians are the most threatened vertebrate group in the world, with many species declining or facing extinction (Stuart et al. 2004;Mendelson et al. 2006).Frogs are threatened by a wide range of processes including climate change, pollution, habitat alteration and introduced species (Laurance 2008;Alton & Franklin 2017;Hayes et al. 2010;Cushman 2006;Kats & Ferrer 2003).Infectious disease is the greatest threat to amphibians since chytridiomycosis has caused the decline of over 500 species globally (Bower et al. 2017;Scheele et al. 2019).This disease is caused by the chytrid fungus, Batrachochytrium dendrobatitidis, and at least 43 frog species in Australia have declined in population size due to this pathogen (Berger et al. 1998;Skerratt et al. 2007;Scheele et al. 2017;Scheele et al. 2019).Long term population trajectories are highly variable for Australian frogs threatened by chytrid.Some species continue to decline while others appear to be stabilizing or recovering.For other Australian species, the severity of initial declines and changes to distribution in temperate Australia remain unclear due to a lack of base-line information (Scheel et al. 2017).It is also not clear to what extent poor monitoring of species has caused some declining populations to appear stable (Scheele et al. 2017).With so much uncertainty, there is a pressing need to obtain accurate estimates of frog population sizes and estimate gene flow for poorly monitored threatened species.
Within conservation research, DNA sequencing technologies are becoming increasingly popular as these technologies enhance and complement traditional population monitoring techniques (Hohenlohe et al. 2021).Genetic methods typically require lower sampling effort and can provide insight into behaviours, such as dispersal and breeding, making them highly suitable for species that are difficult to observe 1 3 (Amato et al. 2009;Wang et al. 2009;McCartney-Melstad & Shaffer. 2015).In recent years, nuclear genome wide sequencing to identify single nucleotide polymorphisms (SNPs) have been used to greatly expand our knowledge of population connectivity, dispersal patterns, effective population sizes, mating systems and genetic diversity for several rare and cryptic species (Allendorf et al. 2010;Bradford et al. 2018;Harrison et al. 2019;Miller et al. 2018).By employing these methods alongside traditional fieldwork, researchers can develop more specific and informed conservation actions for species facing extinction (Frankham et al. 2019).
Littlejohn's treefrog (Litoria littlejohni) is one Australian amphibian of conservation concern for which population dynamics are uncertain.This moderate sized (adult body length 42 mm) brown tree frog occurs in the temperate climatic region of south-east Australia.Recently, L. littlejohni was redefined and taxonomically split into two species, and subsequently its conservation status was assessed as Endangered (EN) (International Union for the Conservation of Nature Red List criteria, IUCN 2012) by Mahony et al. (2020).What is now considered L. littlejohni is found in three isolated regions in the state of New South Wales (NSW) (Mahony et al. 2020).The geographic range is restricted to the Sydney Basin bioregion (Thackway & Cresswell 1995) with apparently isolated populations on the Woronora Plateau in the southeast, to the Blue Mountains in the west, and to the Central Coast Range (Watagan Mountains) in the northeast (Fig. 1).It is unknown whether the populations found in the Blue Mountains, Central Coast Range and Woronora Plateau are connected or isolated from one another.Furthermore, the degree of gene flow between occupied sites within these populations is unknown.
For species, such as L. littlejohni, that have a fragmented distribution across their range, it is important to consider patterns of gene flow and genetic processes such as Isolation by Distance and Isolation by Environment.Isolation by Distance refers to the process in which genetic differences increase with geographical distance as dispersal tends to decrease with distance (Wright 1943).Populations can also become more genetically distinct due to increases in environmental differences independent of geographical distance, which is termed Isolation by Environment (Wang & Summers 2010;Bradburd et al. 2013;Sexton et al. 2014;Wang & Bradburd 2014).When gene flow is restricted between populations, several genetic changes can occur in the resulting isolated populations, including reductions in effective population size, loss of genetic variation and increased local inbreeding within the isolated population (Frankham 2010).Geographically isolated populations are more likely to be vulnerable to stochastic environmental events including drought and fire, and stochastic demographic events including non-random mating, thus experiencing higher extinction risk (Frankham 2010).
The rarity of L. littlejohni has presented a scientific quandary.Determining the cause of rarity requires an understanding of population sizes and dynamics, which is difficult when little is understood about the species.Typically, L. littlejohni is found in dry sclerophyll or heath forests that are widespread community types across the east coast of NSW and the precise vegetation requirements remain unknown (White et al. 1994;Lemckert 2010;Klop-Toker et al. 2021).Additionally, exact breeding preferences remain poorly understood (Klop-Toker et al. 2021).It has only recently been confirmed that breeding peaks in late winter to early spring (daily mean minimum 6 °C and maximum 18 °C on Woronora plateau) (Gill et al. 2021).This poor understanding of habitat preferences and calling behaviours is thought to have led to mis-matched survey efforts and false absences in historical surveys conducted across the range (Gill et al 2021).Populations of L. littlejohni also do not appear to be large and densities at breeding ponds are relatively low compared to other Australian tree frogs, resulting in low capture rates at occupied sites (Mahony et al. 2020;Klop-Toker et al. 2021;Gillespie et al. 2016).Whether this low density is due to inherent low population density or due to declines is unknown (Lemckert 2004;Mahony et al. 2020).Litoria littlejohni is susceptible to chytrid with prevalence levels approximately 30% across all L. littlejohni populations (Klop-Toker et al. 2021), but there is no pre-chytrid population data for this species.There is evidence of gradual declines, in individual populations, particularly in the Blue Mountains where only three breeding locations on the King's Tablelands have been confirmed in the past decade despite intensive survey efforts (Mahony et al. 2020).While there is evidence of L. littlejohni declines, there has been no formal estimate of population size for any population.
In this study, we aimed to use genetic data (SNPs) to 1) elucidate patterns of gene flow across fragmented populations, 2) determine the size of L. littlejohni populations 3) estimate genetic diversity and inbreeding within populations, and 4) outline priority conservation and research actions to promote long-term persistence of L. littlejohni.

Sample Collection and Selection
We conducted surveys for L. littlejohni between May 2018 and June 2021 in three regions of NSW that span the known distribution of the species-the Blue Mountains, Central Coast Range, and Woronora Plateau (Fig. 1).Within the Woronora Plateau, 18 × 100 m transects were distributed across four water catchments-O'Hare's Creek, Cordeaux River, Avon River, and the Cataract River (Fig. 1).Within the Blue Mountains and Central Coast Range, we established two and nine sites respectively at fire dams or large ponds, rather than 100 m transects.Ponds at these two locations were selected instead of stream transects because L. littlejohni have been detected more commonly at dams and ponds in these two regions and were therefore considered more reliable sites for catching adult frogs or tadpoles.
Surveys were conducted using standardized methods across all seasons with the total number at a site ranging from one to 13 occasions.Sites with low survey occasions were typically established late into the project upon discovery of new occupied sites.Surveys for adults were conducted at night and involved a timed active search of stream/pond bank habitat aided by spotlights.This initial search was followed by a repeated active search in conjunction with call-playback to help locate male frogs.All adults caught were microchipped to reduce the likelihood of double sampling in future surveys and to allow the collection of recapture data.A hand-held GPS was used to collect the precise location (latitude and longitude) of frog collection to the nearest 5 m or the latitude and longitude of the site.We collected genetic samples by taking a 4 mm biopsy from the rear foot webbing of frogs with a snout to vent length (SVL) over 40 mm and stored samples in 70% ethanol.
We also collected genetic samples and GPS locations from L. littlejohni found opportunistically outside of surveyed transects or ponds to boost our sample size and geographical reach.During the day we visually inspected ponds for tadpoles at each site.When tadpoles were detected, we collected a subset of tadpoles by dip-net and tail tips were collected for genetic analysis.Tail tips were only collected from tadpoles longer than 18 mm in body length.Tadpole samples were included in genetic analysis if they were from sites where the adult sample size was low (< 6).We selected tadpole samples from different ponds along 100 m transects to reduce the likelihood of sequencing related individuals.Additionally, we included 26 samples that were collected and sequenced by Mahony et al. (2020) to extend the geographic reach and increase sample size.The main difference in methodology between the sample collection of Mahony et al. (2020) and that of the present study was the addition of microchipping adults in the present study (Table S2).While most of the samples from Mahony et al. (2020) came from eight additional sites, we included samples from three sites also surveyed during the present study to boost sample size at the locations.Overall, samples were collected from 37 locations-15 in Cordeaux River, 9 in the Central Coast range, 6 in the O'Hare's Creek, 3 in the Blue Mountains, 2 site in Avon River, 1 in the Nepean River and 1 in the Cataract River (Fig. 1).Additionally, we included two double blinds to help detect whether we double sampled individuals due to microchip loss.
Litoria watsoni is the most closely related species to L. littlejohni (Mahony et al. 2020) and we used 12 Litoria watsoni samples from the NSW range as an outgroup.These samples were collected across 5 national parks between 2018 and 2020 using the same methods as described for L. littlejohni or were collected as per Mahony et al. 2020 (Fig. 1; Table S1).All genetic samples were stored in 70% ethanol.

DNA sequencing and initial filtering
Samples were sent to Diversity Arrays Technology Pty Ltd (Canberra, Australia) for DNA extraction following methods outlined in Kilian et al. (2012).SNP genotyping and discovery was carried out by Diversity Arrays Technology using the DArTseq protocol (Georges et al. 2018).Some of the major advantages of the DartSeq method is that it can handle lower DNA input and DNA quality than other sequencing methods.This sequencing method does not require a reference genome allowing for poorly understood species (non-model animal and plant species), such as our target L. littlejohni, to be studied and is therefore also cost effective.Briefly, this method employs PstL and Sphl restriction enzymes, custom proprietary bar-coded adapters, PCR amplification, and sequencing using Illumina HiSeq2500 following the methods detailed in Georges et al. (2018).Sequences generated were processed using a proprietary DArT analytical pipeline described by Georges et al. (2018).Samples from this study were co-analysed with Mahony et al. (2020) to ensure that the datasets could be merged.
The DArTSeq sequencing, and filtering pipelines provided a total of 62,368 SNPs.SNPs were further filtered using the DartR V2.3 package in Rstudio V2022.7.0.548 to obtain the highest quality data (Mijangos et al. 2022;RStudio Team 2022).We used gl.filter.locmectric to filter the read depth for both the reference and alternative allele to > 7. We then filtered SNPs to fit the following criteria: reproducibility of > = 99%, no monomorphic loci, SNP call rate > 97% (SNPs not found in 3% of the individuals), and no secondary SNPs.Secondary SNPs are those found within the same sequencing fragment and are likely to be linked (Gruber et al. 2018).Secondaries were removed by filtering out all but the first SNP with the same CloneID.We also excluded individuals with an individual call rate of < 95%, (i.e., individuals with > 5% of SNPs missing).We filtered for minor a frequency using a threshold of 0.00169 to remove overall rare alleles from the dataset.
After initial filtering for SNP quality, we assessed whether there were closely related individuals and double sampled individuals (due to microchip loss) that needed to be removed from the dataset.Related individuals and double sampling of individuals can heavily bias population structure as models typically assume individuals to be unrelated (Patterson et al. 2006;Anderson and Dunham 2008;O'Connell et al. 2019).To conduct this analysis, we created 5 small datasets based on sample region to help easily identify related individuals (Blue Mountains, Central Coast Range, Cordeaux/Avon, O'Hare's/Cataract and L. watsoni).For each of these small datasets, we ran the initial filtering steps and then employed the bitwise.distfunction from the poppr V2.9.3 R package suing default settings.This function calculates pairwise genetic distances between samples (Kamvar et al. 2014(Kamvar et al. , 2015)).We used the double-blind samples to calculate a threshold value for genetic distances.For pairs that differed by less than 0.01%, we removed one of the individuals from our overall data set.Once double sampled individuals were removed, we refiltered the data and we used a genetic relatedness network to assess relationships.We employed the gl.grm.networkfunction from DartR, which uses a similar approach to that used by Goudet et al. (2018).This method is suited to populations for which little is known as it uses the kinship of all pair uses the average inbreeding coefficient as a reference value and is then subtracted from the inbreeding coefficient of each pair of distinct individuals (Goudet et al. 2018).We removed individuals so that no pair of individuals had a relatedness factor ≥ 0.25.This threshold was selected as sibling or parent-offspring relationships have a kinship value around 0.25 (Speed & Balding 2015).When all regions had been assessed, we made two datasets for remaining analysis: Dataset 1, which has the outgroup L. watsoni samples included, and Dataset 2, which has only L. littlejohni samples.

Population structure
To investigate population structure, we implemented the Bayesian clustering approach in STRU CTU RE on Dataset 1 (Outgroup included) using the admixture model and correlated allele frequency due to short sampling distance between the different river catchments on the Woronora Plateau.We assessed values of K from 4 to 10 with 3 independent runs with 20, 000 burnin and 50 000 MCM iterations for each value of K. Output files from Structure were processed with StructureSelector (Li & Liu 2018) and we used the Puechmaille method to select the best K (Puechmaille 2016).The Puechmaille method employs four estimators (MedMeaK, MaxMeaK, MedMedK and MaxMedK) to determine the best number of clusters (Puechmaille 2016).This method helps to alleviate the issue of uneven sampling (Puechmaille 2016), which is present in the current dataset due to small sample size in the Blue Mountains region.An individual is assigned to a cluster if its arithmetic mean (MedMeak and MaxMeak) or its median (MedMeDK and MaxMedK) membership coefficient to that cluster is greater than the threshold used.To assess the performance of the estimators and prevent over estimation of the true number of clusters we ran replicates of the Puechmaille method using increasing thresholds (0.5, 0.6, 0.7, 0.8).Once best K was selected, we ran STRU CTR E again with the preferred K for 20, 000 burnin and 100, 000 MCMC iterations.A plot identifying the population clusters was created using CLUMPAK through the StructureSelector web interface (Kopelman et al. 2015).
We conducted a genetic test of neutrality by calculating Tajima's D for all samples and for each individual population in Dataset 2 (Tajima 1989).This test was conducted using the hierfstat V0.5-10 package in R studio and the TajimaD.dosagefunction after transforming the data into dosage data that reports the number of copies of each allele using the fstat2dos function (Goudet 2005).
Following the test of neutrality, we applied further filtering using the DartR package to Dataset 1 and Dataset 2 to remove loci that were significantly not in Hardy-Weinberg Equilibrium and increased minor allele frequency threshold to 0.02 (RStudio Team 2022; Mijangos et al. 2022).To test for Hardy-Weinberg Equilibrium, we used gl.filter.hwefrom the DartR package with individuals grouped by population cluster.This function uses observed frequencies of reference homozygotes, heterozygotes and alternate homozygotes to filter out loci that show departure from HWE. Regarding settings for this function, we employed the exact test as it is recommended in most cases (Wigginton et al. 2005) and used the selome method to estimate p-values with an alpha value of 0.05.The selome method computes p = values as the sum or probabilities of all samples less or equally likely as the current sample (Graffelman 2015).P-values were corrected for multiple tests using the BY method based on Benjamini & Yekutieli (2001), which controls the false discovery rate.The function gl.report.ld.map, from the Dart R package, was employed to test for linkage disequilibrium within populations.This function reports pairwise linkage disequilibrium (LD) for SNPs for which the LD measure is > 0 in all populations and uses a threshold of R 2 = 0.20 (Delourme et al. 2013;Li et al. 2014).

Population genetic differentiation
We computed pairwise F ST (Weir and Cockerham 1984) using the clusters from the STRU CTU RE analysis to further assess the degree of genetic differentiation between population clusters.This was computed using gl.fst.popfrom the DartR package and with 1000 bootstrap replications to evaluate whether F ST values were significantly different from zero (Gruber et al. 2018;Mijangos et al. 2022).

Isolation by distance and isolation by environment
We conducted a Mantel test of genetic differentiation vs geographic distance between all pairs of L. littlejohni populations, using gl.ibd with 999 permutations from the DartR package on Dataset 2 in R studio (species ingroup only) (Mijangos et al. 2022;RStudio Team 2022).For this test we used Euclidean distance, and we transformed the geographic distances using log(Dgeo + 0.01) as some individuals have identical GPS coordinates that will otherwise result in a geographical distance of 0.
We conducted a partial redundancy analysis in R studio to investigate the influence of environment (Isolation by Environment) and spatial features (Isolation by Distance) on the distribution of genetic variation in Dataset 2 (species ingroup only) (Orsini et al. 2013;Capblanc & Forester 2021: RStudio Team 2022).Partial redundancy analysis cannot handle missing values, therefore, the 0.58% of loci with missing data were substituted with values taken from the nearest neighbouring individual without a missing value using a population-by-population basis.This substitution was conducted using the gl.impute function from the Dart R package.The "neighbour" method as it is recommended for small populations and is nearest neighbour is the individual with the smallest Euclidean distance from the focal individual (Mijangos et al. 2022).The response variable was individual genotypes and we used 3 sets of variables (1) environmental, (2) geography, and (3) genetic structure as explanatory or conditioning variables.To select environmental variables, we downloaded 22 Bioclim variables and elevation data from Atlas of Living Australia Spatial Portal for each individual GPS point (Belbin 2011;Xu & Hutchinson 2011).To determine which environmental variables to include in our partial redundancy analysis we ran a forward selection model using the rda and ordiR2step functions from the R package vegan V2.5.7 to identify variables significantly associated with genetic variation (Oksanen et al. 2013) (Table S2).We then removed environmental variables that were strongly correlated with each other (Pearson's correlation > 0.60) with preference for those that had higher R2 values and AIC in the forward selected model.The final environmental variables were radiation seasonality, precipitation of driest month, minimum temperature of coldest month, and isothermally.Geographic variables used were the GPS coordinates of individual frogs (Latitude and Longitude).Genetic structure variables were the scores along the first five axes of a genetic Principal Coordinate Analysis conducted on Dataset2 using gl.pcoa from the DartR package (Gruber et al. 2018).We selected the first five axis as they all explained > 1% of variance, totalling to 21.8%.Once variables were finalised, we ran RDAs for a full model (Climate, Geography, Genetic Structure) and one for each set of explanatory variables conditioned by the remaining sets to variables (e.g., Climate conditioned by geography and genetic structure).This allows us to partition the percentage of genetic variance explained by each of the explanatory variables.R squared values were adjusted and an ANOVA with 9999 permutations was run for each model to obtain a p-value.

Genetic diversity
We employed gl.report.heterozygosityfrom the DartR package using Dataset 2 (species ingroup only) to calculate the observed heterozygosity, expected heterozygosity, unbiased expected heterozygosity, and F IS .Standard deviation for F IS was obtained using boot.ppfisfrom the hierfstat V0.5-10 package with 100 bootstraps (Goudet 2005).

Mean kinship and inbreeding
Mean kinship is considered one of the best metrics for managing genetic diversity and can help identify the best source populations for augmentation (Frankham et al. 2019).We employed the related V0.8 R package to estimate within population kinship (Pew et al. 2015).There are several relatedness estimators that use the genetic information differently and can yield different estimates (Wang 2011), thus, to determine the best relatedness estimator for our data, we used the function compareestimator from the related package.We implemented this function using a subset of 200 randomly selected loci to compare the seven different relatedness estimators.This function simulates individuals of known relatedness using the loci provided, which then allow users to evaluate the correlation between observed and expected relatedness (Pew et al. 2015).The quellergt estimator, described in Queller and Goodnight (1989) performed the best (Pearson's r = 0.894), thus used to estimate relatedness.We also employed the lynchrd estimator as it calculates individual inbreeding where as the quellergt estimator does not.This measure of inbreeding (F) is based on the probability of alleles at a random locus being identical by descent (Lynch and Ritland 1999), while F IS indicates the proportion of the variance in a subpopulation contained in the individual.Genetic relatedness among all pairs of L. littlejohni in Dataset 2 were estimated using the quellergt and lynchrd estimators using the function coancestry with 1000 bootstraps and inbreeding allowed.We averaged relatedness to get within population mean kinship and we averaged F for each population.

Effective population size
Effective population sizes were estimated using Dataset 2 and the program NE estimator V (Do et al. 2014).We ran NE Estimator using the Linkage Disequilibrium method with random mating and reported the estimates at a critical allele frequency of 0.05 for each sampling region.

Survey results
Over 185 survey occasions, we caught 265 unique individual frogs with 31 frogs (11.7%) recaptured at various times (Table S1).Recapture rates were the highest in the Central Coast Range with two sites accounting for 48.3% of all recaptures.Twenty sites had no recaptures over the study period.The highest number of individuals caught at a site was 30.No adult frogs were ever found at 2 sites, but tadpoles were detected there.

Population structure
Dataset 1 (outgroup included) had 1,049 loci and 312 individuals post initial filtering.All four of the Puechmaille estimators indicated that the optimal K was six clusters (Fig. 2; Fig. S2).This did not change when using higher threshold values.The six clusters are 1) L. watsoni, 2) Blue Mountains,3) Central Coast Range, 4) O'Hare's Catchment, 5) Cataract Catchment, and the final cluster is formed by 6) Cordeaux Catchment, Avon Catchment and the Nepean Catchment (referred to as Cordeaux Cluster).Tajima's D was between zero and one for all populations.The Cordeaux Catchment and Central Coast Range populations had values below 0.50 (Table 3).Post additional filtering, Dataset 2 had 981 loci and 299 individuals.

Population genetic differentiation
Pairwise F ST values indicate that the Central Coast Range is the most distinctive L. littlejohni population (Table 1).The Central Coast Range population had the highest pairwise F ST with all other L. littlejohni populations.Interestingly, the F ST value comparing the Central Coast Range and the O'Hare's catchment indicated that these populations are almost as genetically differentiated as L. watsoni and O'Hare's.The same pattern was seen for Central coast range and the Cataract cluster.The Blue Mountains is the second most genetically distinctive L. littlejohni population with pairwise F ST values all above 0.10.The three populations on the Woronora Plateau had pairwise F ST values over 0.05.This value indicates that these are moderately genetically different from each other, however, we detected some movement between population clusters in the structure plot.Three individuals sampled from the O'Hare's catchment had approximately 50% membership to the Cataract clusters.

Isolation by distance and isolation by adaptation
The Mantel test indicated a significant association between geographical and genetic distances (Mantel's r = 0.658,

Genetic diversity, effective population size and kinship
Effective population size was relatively low across all populations with the highest value below 200 in the Cordeaux Cluster (Table 3).The Blue Mountains, Central Coast Range and Cataract populations all had effective population sizes below 50, with the Blue Mountains having the lowest.These smaller populations also had lower observed heterozygosity than those with effective population size over 50 (Cordeaux and O'Hare's).For all populations, the expected heterozygosity was higher than observed, but this pattern was more pronounced for populations with effective population size below 50.Mean kinship, F, and F IS were highest in the Central Coast Range, Blue Mountains, and Cataract catchment.

Discussion
Obtaining estimates of population size and connectivity is a critical step within conservation research, however, for many threatened amphibians this is a difficult task due to the cryptic nature or rarity of such species.In this study of L. littlejohni, where we were unable to obtain sufficient recaptures to estimate population sizes through traditional mark-recapture methods, however, we were able to estimates effective population size in five genetically distinct populations using genetic methods.Additionally, we revealed the restriction of gene flow at the regional and local scale.This additional knowledge would be unlikely obtained through standard capture-recapture methods.Two population clusters are clearly geographically and genetically isolated, but the other three are located on the same plateau.Climate and geography were not strong drivers of genetic isolation.High levels of inbreeding, small population size and reduced genetic diversity were detected in the two geographically isolated populations and one of the populations on the Woronora plateau.These results provide strong evidence that the rarity of this frog is likely due to small population sizes rather than low detection rates and urgent conservation action is needed to ensure the persistence of this threatened frog.

Population structure
When comparing the major sampling regions, strong genetic isolation was observed in L. littlejohni.Based on the results of the pRDAs, we suggest that genetic drift due to population contractions is the main driver of genetic variation across populations, rather than large climatic difference or geographic distance between populations (Capblancq & Forester 2021;Coleman et al. 2013).While Isolation by Distance was significant, the partial RDA indicated that only a small percentage of genetic variation could be attributed to geography.In the partial RDA models, population genetic structure was the strongest predictor of genetic variation indicating that demographic history rather isolation by adaption or isolation by distance have driven patterns of genetic variation in L. littlejohni (Capblancq & Forester 2021).However, these variables were strongly confounded and over 63% of variation in the full model could not be exclusively attributed to climate, geography or genetic structure alone.This suggests any future work to investigate loci under selection or conduct further landscape genetic analysis in L. littlejohni must consider ways to account for this collinearity (Frichot et al. 2015;Hoban et al 2016;Capblancq & Forester 2021).The test of neutrality further supports the notion of genetic drift following population crash as the best explanation of modern population differentiation in L. littlejohni.All populations had positive Tajima's D values, which may indicate very balancing selection or sudden population contraction (Tajima 1989).However, all values had quite low absolute value indicating that any selection would be weak.
Increased genetic structuring of populations has been linked to disease-related population contractions in numerous vertebrate species.In these species, population contractions led to increased population fragmentation and exacerbated genetic drift causing populations to become more genetically differentiated (Phillips et al. 2020;Serieys et al. 2015;Lachish et al. 2011).Both the Central Coast Range and Blue Mountains have effective population sizes below 50, high inbreeding values and lower genetic diversity compared to other L. littlejohni populations, indicating that these populations have indeed suffered from population bottlenecks.Between these two regions are largely continuous natural areas consisting of seemingly suitable habitat where there were historical sightings of L. littlejohni, thus these regions were once connected (Mahony et al. 2020).Due to a lack of pre-chytrid genetic samples, we cannot confirm chytrid as the driver of L. littlejohni declines in this study, although this disease is the most likely driver of sudden decline in this winter active, stream and pond frog (Hero 2005).However, it seems likely that these two populations have become further isolated both geographically and genetically due to disease driven population declines with the resulting small population sizes exacerbating the consequences genetic drift.
Large scale damming and population crashes due to disease may explain the strong genetic isolation between the Blue Mountains and Woronora Plateau populations.Damming of riverine systems alters population connectivity by fragmenting suitable habitat (Heggenes & Roed 2006;Emel & Storfer 2012;Hansen et al. 2014;Werth et al. 2014).In other frog species, river regulation has led to changes in frog species composition, increased genetic isolation, and loss of genetic diversity (Naniwadekar & Vasudevan 2014;Peek et al. 2021, Grummer & Leache 2017).We hypothesize that the construction of the Warragamba Dam in the 1960s may have disrupted gene flow between the Blue Mountains and Woronora populations.The Warragamba dam creates Lake Burragorang, however, L. littlejohni does not occupy large lakes.Rather, this frog is typically associated with first order streams and small fire dams (Mahony et al. 2020;Lemckert 2010).Thus, the lake may act as barrier due to unsuitable habitat.
The degree of genetic structuring observed on the Woronora Plateau was unexpected as the landform is contiguous and there is little urban development beyond major roads.This region has the largest population of L. littlejohni, thus it was anticipated that gene flow would occur between the sampled river catchments.Instead, we detected three genetic clusters with low to moderate genetic differentiation (O'Hares, Cataract and Cordeaux; Fig. 2).Distinct genetic structure at a local scale is typically seen in species that have high site fidelity and inherently low dispersal abilities (Zickovich & Bohonak 2007;Richardson et al. 2021), however, we do not think this is the case for L. littlejohni.The geographic size of the Cordeaux cluster (> 84km 2 ) indicates that L. littlejohni is not limited by its ability to disperse across substantial distances.We propose three factors that may explain the fine-scale structuring of L. littlejohni populations: 1) roads that prevent movement between catchments; 2) dispersal is restricted to creek lines; and 3) habitat disturbance that has isolated populations locally.As discussed with the other populations, the genetic structuring may also be exacerbated by general distribution-wide declines due to chytrid.
Roads (1) restrict movement of a variety of vertebrate species, leads to genetic isolation of populations (Hale et al. 2013;Holderegger & Di Giulio 2010;Serieys et al. 2015;Youngquist et al. 2017).For example, major urban infrastructure (housing and dual carriage roads) has been linked to the genetic sub-division of southern bell frog (L.raniformis) (Hale et al. 2013), and the presence of highways is predictive of genetic distance in Blanchard's cricket frogs (Acris blanchardi) (Youngquist et al. 2017).The Woronora Plateau has several large dual carriageway roads which traverse it (Picton Road and Appin Road) (Fig. 1), which are possible barriers to dispersal in L. littlejohni.
Dispersal primarily along creek lines (2) rather than through terrestrial landscapes may contribute to isolation of populations.For some aquatic and semi aquatic species, patterns of gene flow reflect the hierarchical nature of river systems and/or the isolation of populations to independent river basins or drainages (Brauer et al. 2016;Coa et al. 2020;Fonseca et al. 2021).Whilst the river catchments on the Woronora Plateau are geographically close to each other, they are parts of different major river drainages (Fig. 1).The Cordeaux Cluster and Cataract Catchment flow into the Hawkesbury-Nepean River system, while the O'Hare's Catchment flows into the Georges River/Tucoerah River.O'Hare's and the Cataract site is essentially the border of these two major river catchments.Although we did detect movement between Cataract and O'Hare's sites, the moderate genetic differentiation of these two clusters indicates movement is not common despite the short distance between them (less than 5 km between Cataract and closest O'Hare's site).As L. littlejohni tadpoles are aquatic, larval dispersal along creeks (particularly during high rain periods) is highly probable.Additionally, during our surveys we often encountered adult frogs utilizing the sandstone bedrock of creek lines to move between water bodies.We only detected a frog more than thirty metres from water once, despite accessing sites by walking.There have been reports of L. littlejohni found under rocks away from streams during other fauna surveys (G.Daly unpubl.data), thus some terrestrial movement does occur, but it is suspected to be limited.
In addition to roads, the main habitat disturbance (3) occurring in the Woronora is longwall coal mining.Longwall mining can lead to the diversion of water underground due to subsidence causing cracks in the surface rock (Booth 2006).The Cataract River Catchment has been undermined extensively, however, to what degree this mining has impacted frog species is unclear due to a lack of research (WaterNSW 2016).In the last ten years, sightings of L. littlejohni within Cataract have been limited in number and are restricted to small pockets within the catchment.It is plausible that mining activity has fragmented populations locally as drying of creeks may result in limited successful dispersal events.However, longwall mining also occurs within the Cordeaux Cluster, which has no population sub-structuring, thus there is no clear trend as to how creek drying due to mining impacts frog movement.

Small populations and inbreeding
The small effective population size across all L. littlejohni populations indicate that the rarity of this species is indeed due to historical declines and not just strict habitat preference or cryptic behaviour.These estimates are also of conservation concern as the fate of fragmented populations is determined by the effective population size of the isolated population, rather than the species as a whole (Frankham et al. 2019).While the effective population sizes only reflect that of the sampled areas, as we sampled a large proportion of known currently occupied sites, thus it is unlikely that these estimates would change dramatically with additional samples.The three smallest populations (Central Coast range, Blue Mountains, and Cataract) also have higher inbreeding coefficients for both F and F IS , lower observed heterozygosity and higher mean kinship than the O'Hare's and Cordeaux clusters.Together these factors put populations at risk of reduced adaptive potential, inbreeding depression, and increased sensitivity to stochastic events (Frankham et al. 2019).Additionally, populations may be impacted by mutational meltdown, a type of extinction vortex (Keller & Waller 2002).Mutational meltdown is a process in which genetic and demographics processes mutually reinforce one another as deleterious mutations become fixed and accumulate within the population.These mutations lead to further reductions in population size and in turn increases the rate that mutations are accumulated and the speed at which populations decline in size (Lynch & Gabriel 1990).Alternatively, deleterious alleles may be purged from populations with high inbreeding and small population size.Inbreeding is expected to increase the expression of deleterious alleles as more individuals will be homozygotes for recessive alleles.Natural selection should select against such individuals in the population, leading to less individuals overall with the recessive deleterious alleles.However, the extent of genetic purging depends on many genetic factors and the environment in which it takes place (Keller & Waller 2002).
The coefficient of inbreeding recommended by Frankham et al. (2019), F, uses identity by descent methods to estimate cumulative inbreeding over generations (Frankham et al. 2019).On the other hand, F IS is used more widely in population genetic research, but only detects deviations from random mating in recent generations.In the present study, we have populations in three different states (1) High F and high F IS , (2) Low F and low F IS , and (3) High F but low F IS .The Cordeaux cluster is in State 1 (Low F and F IS ) and has the largest population size, indicating that this population may not have faced large population declines or the population maintained a large enough size to prevent changes in random mating patterns.The Central Coast Range, Blue Mountains and Cataract are all in State 2 (High F and F IS ) and have small population sizes, indicating that these populations have declined and not recovered.The final population, O'Hare's, is in State 3 (High F, but low F IS ) and has a small but larger population than those in State 2. The difference in inbreeding values may indicate that the current population is not deviating from random mating, but that in the past the population was small enough for inbreeding to be present.Thus, O'Hare's likely experienced decline and has recovered to some extent.There are several examples of chytrid susceptible species that have declined in some parts of their range but persist or have recovered in others (Osborne et al. 1996;Retallick et al. 2004;Wassens 2008;Hamer et al 2010;Mahony et al 2013;Newell et al. 2013;Lips 2016).Habitat quality, changes to disease dynamics, recruitment rates and population connectivity have all been linked to stabilizations and recoveries of frog species (Scheele et al. 2014;McDonald et al. 2005;Scheele et al. 2017;Mc Knight et al. 2019).As L. littlejohni populations appear to differ in recovery history they may provide a useful model for further investigation into drivers of recovery.

Management recommendations
Through the application of genetic methods, conservation managers will be able to make decisions for L. littlejohni with more clarity going forward.The small effective populations, high inbreeding, reduced genetic diversity and high mean kinship detected in some L. littlejohni population indicates that urgent on-ground action is needed to increase population sizes.Determining habitat preferences and protecting vital breeding habitat will be important for supporting successful recruitment within populations.Recommended research includes assessing whether inbreeding depression is present in populations with high F IS and F. Inbreeding depression in other frogs has been linked to reduced sperm quality and reduced offspring survival (Hinkson & Poo 2020;Anderson et al. 2004), which obviously has implications for population viability of threatened species.Additionally, we recommend assessing whether assisted gene flow will aid in genetic rescue or lead to outbreeding depression (Frankham et al. 2019).The Cordeaux and O'Hare's clusters are recommended as potential source populations as they contain the highest genetic diversity, lowest inbreeding, and largest populations.Additionally, the metrics for coefficient of inbreeding and mean kinship have both been used in other species to assess or plan translocations (Cowen et al. 2021;Farquharson et al. 2021), thus the values reported here provide a baseline for future conservation management.

Global significance/conclusions
This study demonstrates the power of pairing genetic methods with traditional surveys to study rare and difficult to track species.Without employing genetic methods in the present study, we could have not resolved the issue of low sampling rates and low recapture rates that prevented successful mark-recapture methods.Additionally, the SNP data provided insight into potential inbreeding depression, which highlights an urgent need for conservation actions and further research, that would not have been detected through visual encounters surveys alone.Incorporating genetics into conservation research can have a profound impact on both our understanding and our ability to provide adequate management of threatened species.
need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http:// creat iveco mmons.org/ licen ses/ by/4.0/.

Fig. 1 a
Fig. 1 a GPS points for all L. littlejohni and L. watsoni samples included in population genetic analysis represented by circles and triangles, respectively.b Zoomed in map of Woronora plateau showing river catchments in black lines, creeks in grey and sampled L. littlejohni individuals as triangles c Map of Woronora Plateau showing major roads in dark grey and black, sampled L. littlejohni indi-

Table 1
Pairwise F ST values between population clusters identified in structure analysis Variance is partitioned into pure climatic, pure geography and pure structure.The Proportion of Explainable Variance corresponds to the partitioned variance relative to the constrained variance of the full RDA model (Spatial + Climate + Structure) ***p ≤ 0.001

Table 3
Genetic Diversity metrics for Dataset 2 L. littlejohni with population clusters determined by fastStructure The number of individuals followed by number of loci, observed heterozygosity (H O ), expected heterozygosity (H E ), Unbiased expected heterozygosity (H E -U), inbreeding coefficient (F IS ), Co-efficient of inbreeding identity by descent method (F), Mean Kinship, and Effective population size (N E ) .Values in brackets indicate standard deviation (SD) or 95% confidence interval (CIs)