The banjo frogs or ‘pobblebonks’ are a distinctive and charismatic group of medium to large (34–94 mm body length) terrestrial and burrowing limnodynastid frogs endemic to Australia. Often referred to as the Limnodynastes dorsalis group (Martin, 1972; Roberts and Maxson, 1986; Schäuble et al., 2000), species in this group are recognizable by their conspicuous, single-note, resonant “bonk” or “tok” advertisement calls, which are often likened to the pluck of a banjo string. Eight taxa are recognized currently, including four species and five subspecies: Limnodynastes dorsalis, L. interioris, L. terraereginae, and L. dumerilii, which is further divided into the subspecies L. dumerilii dumerilii, L. d. fryi, L. d. insularis, L. d. variegatus, and L. d. grayi. Except for L. dorsalis, which is restricted to south-western Western Australia, all other taxa occur from southeastern to northeastern Australia where they share mostly parapatric distributions with some zones of sympatry (Martin, 1972). Hybrid zones are known to form at contact zones between multiple species, based on the presence of intermediate morphological and advertisement call phenotypes within these zones (Martin, 1972). The extent of hybridization between species at contact zones has been found to correspond particularly to divergence in advertisement calls: species with the most similar calls are thought to hybridize more frequently. Their distributions span a variety of mesic to semi-arid habitats, with two species, L. terraereginae and L. dumerilii grayi, even occurring in the highly acidic wallum wetlands of the eastern Australian coastline. Tadpoles of these species possess a remarkable tolerance for acidic waters and can complete their development in aquatic environments with pH levels as low as 3.0 (Hines and Meyer, 2011; Hird et al., 2022).

The phylogenetic relationships among species of the L. dorsalis group have been inferred previously using micro-complement fixation (Roberts and Maxson, 1986) and mitochondrial DNA sequences (Schäuble et al., 2000), with the findings of both studies suggesting the current taxonomy for the group, based primarily on morphology and advertisement call structure, may not best reflect the true evolutionary relationships among taxa. In particular, the genetic data indicate that L. dumerilii grayi is deeply divergent from topotypic L. dumerilii (Roberts and Maxson, 1986) and may be more closely related to L. terraereginae (Schäuble et al., 2000).

In this study, we aimed to clarify the phylogenetic relationships among taxa of the L. dorsalis group and re-evaluate the systematics and taxonomy of L. terraereginae sensu lato. We use comprehensive geographic sampling to explore divergence in mitochondrial and nuclear DNA, morphology, and advertisement calls between the taxa. Further, using targeted surveys of potential contact zones, we sought to document the interactions between taxa and quantify the extent of admixture occurring between them.

MATERIALS AND METHODS

Sampling.—We carried out targeted surveys for all eastern Australian members of the Limnodynastes dorsalis group in New South Wales (NSW), Queensland (QLD), Victoria (VIC), and Tasmania (TAS) between 2020–2022. Surveys were

conducted during spring and summer months and were timed to coincide with rainfall, when the species are most detectable. Our fieldwork focused on collecting specimens, genetic samples, and acoustic recordings of male advertisement calls to facilitate integrative taxonomic analyses. Newly collected samples were combined with a comprehensive coverage of tissues and specimens available in the collections of the Australian Museum (AMS), South Australian Museum (SAMA), Australian Biological Tissue Collection (ABTC), Queensland Museum (QM), and Museums Victoria (NMV). Collection locations of samples used in the genetic analyses are shown in Figure 1 and their details are listed in Table 1. Details of samples used for the morphological analyses are presented in Supplementary Table S1 (see Data Accessibility).

Fig. 1

Distribution of genotyped samples for the eastern Australian Limnodynastes dorsalis group examined in this study. Purple circles = CYP lineage; red circles = WS lineage; green circles = EC lineage; gray circles = L. dumerilii (including all subspecies); and black circles = L. interioris. Surface hydrology in blue.

Table 1

Limnodynastes dorsalis group specimens examined for molecular genetic analyses. Institution codes: ABTC—Australian Biological Tissue Collection, South Australian Museum; AMS—Australian Museum; NMV—Museums Victoria; SAMA—South Australian Museum; QM—Queensland Museum; TMAG—Tasmanian Museum and Art Gallery. HS—homestead, NP—National Park, SF—State Forest, NR—Nature Reserve. ND4 column refers to GenBank accession numbers.

Table 1

Continued.

Table 1

Continued.

Table 1

Continued.

Our fieldwork focused on two primary objectives: (1) genotyping unsampled populations to confirm the distribution of lineages, including collecting topotypic samples to enable interpretation of the nomenclature, and (2) identifying and sampling areas of actual or potential contact between lineages to evaluate admixture between them. In addition to the advertisement calls recorded during fieldwork, we used call recordings obtained from the Australian Museum's national citizen science project, FrogID (Rowley et al., 2019;  https://www.frogid.net.au).

Mitochondrial DNA extraction and analysis.—Nucleotide sequences of the mitochondrial NADH subunit 4 (ND4) gene were obtained from all members of the Limnodynastes dorsalis group, including L. terraereginae (n = 56), L. interioris (n = 5), L. dorsalis (n = 2), and representative topotypic samples of all subspecies of L. dumerilii (n = 21), with a focus on the eastern coastal NSW subspecies L. dumerilii grayi (n = 18). DNA was extracted from ethanol-preserved tissues (liver, muscle, or toe tip) using a DNeasyt Blood and Tissue Kit (QIAGEN GmbH, Hilden, Germany), following the manufacturer's protocols for purification of genomic DNA from animal tissues. A fragment of the ND4 gene was PCR amplified and directly sequenced using the primers: 5′–TGA CTA CCA AAA GCT CAT GTA GAA GC–3′ and 5′–GGT YAC GAG YAA TTA GCA GTT CT–3′. PCR was carried out in 25 µL reactions with a final concentration of 2,000 ng of template DNA, 1X MyTaq^TM Red Reaction Buffer, 2 pmol of each primer, 0.5 units of Bioline MyTaqTM Red DNA Polymerase, and 16.8 µL of autoclaved water. Thermocycling was performed on an Eppendorf Mastercycler EpS (Eppendorf, Hamburg, Germany) under the following cycling protocol: (1) initial denaturation at 94°C for 3 min, (2) 10 cycles involving a denaturation step of 94°C for 45 seconds, annealing at 60°C for 1 min, and extension at 72°C for 1 min, with the annealing temperature decreased by 1°C per cycle, (3) 34 cycles of 94°C for 45 sec, 50°C for 1 min, and 72°C for 1 min, and (4) a final extension step of 72°C for 6 min with samples kept at a holding temperature of 11°C. Amplification products were visualized on 1.5% agarose gels, purified using ExoSap-ITTM (USB Corporation, Cleveland, Ohio, USA), and sequenced in both directions at Macrogen (Seoul, South Korea). Sequence chromatograms were edited and checked for quality using Geneious Prime (v. 2023.1.2;  https://www.geneious.com). Sequences were deposited in GenBank (GenBank accession numbers listed in Table 1).

New sequences were combined with sequences downloaded from GenBank (originally published by Schäuble et al., 2000). The phylogeny was rooted using two sequences each for Limnodynastes depressus, L. fletcheri, L. peronii, L. salmini, and L. tasmaniensis. Sequences were aligned using MAFFT (Katoh et al., 2002), implemented in Geneious Prime v. 2023.1.2. A best-fit partitioning model was inferred using ModelFinder (Kalyaanamoorthy et al., 2017) using the IQ-TREE webserver (Trifinopoulos et al., 2016;  http://iqtree.cibiv.univie.ac.at), using the Bayes Information Criterion (BIC). Maximum-likelihood phylogenetic analyses were also performed with IQ-TREE with branch support assessed using 1,000 ultrafast bootstrap replicates. We considered branches receiving ≥70% bootstrap support to be well supported following Hillis and Bull (1993).

Net average sequence divergence between lineages (dA) was calculated in MEGA 11 version 11.0.13 (Tamura et al., 2021) as: dA = dXY – (dX + dY)/2, where dXY is the average distance between groups X and Y, and dX and dY are the within-group means.

SNP data filtering.—Samples were submitted to Diversity Arrays Technology (DArT Pty Ltd, Canberra, ACT, Australia) for commercial DNA extraction and DArTseq^TM 1.0 genotyping (Kilian et al., 2012). Diversity Arrays Technology uses a combination of genome complexity reduction methods and next generation sequencing platforms. DNA samples were processed in restriction enzyme digestion/ligation reactions using a combination of the PstI/SphI restriction enzymes, and ligated fragments were PCR amplified and sequenced as described by Kilian et al. (2012) and Mahony et al. (2021a, 2021b).

The data were converted to a matrix of SNP loci by individuals, with the contents stored as integers 0, homozygote, reference state; 1, heterozygote; and 2, homozygote for the alternate state. DNA sequences and statistics such as call rate, polymorphic information, heterozygosity, read depth, and reproducibility for all loci and individuals were also reported. Diversity Arrays Technology reports associated with this study include DFr21-6065, DLimno21-6282, and DLimno22-6847.

Due to slight differences in the methods used by Diversity Arrays Technology in the sequencing of each submission, we observed initial biases in the data correlated with each sequence submission. To account for this, we reassembled the short read sequences obtained from Diversity Arrays Technology to reduce any potential adapter contamination using Trimmomatic v0.39 (Bolger et al., 2014) and Stacks v2.62 (Catchen et al., 2013; Rochette et al., 2019). Initially, we used Stacks::process_radtags to remove bar-code sequences from each file. We then used Trimmomatic to filter adapter sequences, using the following parameters: ILLUMINACLIP:TruSeq3-SE.fz:2:20:10; LEADING:5, SLIDING-WINDOW:4:5; MINLEN:68; CROP:68. To standardize read lengths among sequencing submissions, we truncated each read to 68 base pairs using the crop parameter as read lengths were far longer in the more recent sequencing submission. Following this, we completed the core Stacks pipeline de novo, using mostly default parameters, apart from at the ‘ustacks’ stage setting ‘m’ (minimum read depth required to make a stack) to 4, and in the ‘cstacks’ stage setting ‘n’ (the number of mismatches allowed between loci of different samples when assembling the loci catalog) to 2, based upon recommendations in Paris et al. (2017). We output a VCF file during the ‘populations’ process, which was then converted to be compatible with the R package dartR (Gruber et al., 2018) for further analyses.

The SNP data and associated metadata were read into a genlight object (Jombart et al., 2010) to facilitate processing with dartR. Only loci with 100% repeatability (reproducibility) were chosen for subsequent analysis. Further filtering was undertaken based on having a call rate <95% and the locus being present in at least 70% of individuals. We retained only one SNP from each locus at random. Any monomorphic loci arising because of the removal of individuals were also deleted. Given the low within-population sample sizes (n ≤ 15), we did not filter loci for departures from Hardy-Weinberg equilibrium or linkage disequilibrium.

SNP analyses.—We used several approaches to visualize genetic similarity among individuals and detect potential admixture in the SNP dataset. First, we visualized genetic clusters using the principal coordinates analysis (PCoA) ordination method implemented in the gl.pcoa and gl.pcoa.plot functions of dartR. We used a scree plot of eigenvalues to determine the number of informative PC axes to examine using the gl.pcoa.scree function. We then constructed a phylogenetic representation of the SNP data using the neighbor-joining method (Saitou and Nei, 1987) via the gl.tree.nj function of dartR.

In instances where the SNP data appeared to cluster geographically within a lineage, we tested for isolation by distance (IBD). Most localities were represented by a single individual in our dataset, rendering interpopulation analyses using F_ST/1–F_ST as the genetic distance measure unreliable. Instead, we calculated an individual-based genetic distance matrix based on the proportion of shared alleles (D_SA) between individuals, which has been shown to be a powerful statistic for detecting IBD (Sere et al., 2017). D_SA was calculated for individuals with the gl.propShared function of dartR, with log-transformed Euclidean geographic distances calculated between collection localities for individuals in the Mercator projection with the dist function of the stats package in R. We then performed a Mantel test between pairwise geographic and genetic distance matrices via the mantel.randtest function of the adegenet package (Jombart et al., 2010), assessing significance through 999 permutations.

Secondly, we used the Bayesian model-based clustering algorithm implemented in STRUCTURE (Pritchard et al., 2000) to identify population structure and detect admixture. We tested hierarchical support for the lineages identified in the mtDNA and PCoA analyses using a balanced subset of individuals from a spread of the range of each taxon. This was done to minimize bias introduced when using an uneven sample size (i.e., Puechmaille, 2016). STRUCTURE runs were implemented via the gl.run.structure function in dartR, using the uncorrelated allele frequency and admixture ancestry models to assess values of K from 1 to 5, performing three independent runs with 20,000 burn-in and then 50,000 MCMC iterations for each value of K. The preferred K value was determined using L(K) and the change in the second order of likelihood, ΔK (Evanno et al., 2005), obtained using the gl.evanno function. We then ran ten independent runs with the preferred K for 20,000 burnin and 100,000 MCMC iterations and summarized the individual ancestries across all ten runs in CLUMPAK (Kopelman et al., 2015) implemented via the gl.plot.structure function.

We further assessed divergence between clusters identified in the PCoA and STRUCTURE analyses by determining the number of loci showing fixed allelic differences between them. Fixed difference at a locus occurs when two populations share no alleles. When many loci are examined, and sample sizes are finite, fixed differences will occur through sampling error. We used simulations implemented in dartR (Georges et al., 2018) to estimate the expected false positive rate in pairwise comparisons using the gl.fixed.diff. We used tloc = 0.05, meaning that SNP allele frequencies of 95,5 and 5,95 percent were regarded as fixed when comparing two populations at a locus.

Finally, we used the program NewHybrids (Anderson and Thompson, 2002) to identify F₁, F₂, or backcrossed hybrid individuals in the dataset. Two parental lineages were compared per run with parental reference states identified through the PCoA and STRUCTURE clustering. We selected a subset of 200 loci that were most informative in assessing hybridization, namely loci that showed fixed differences between the parental populations, using the “AvgPic” method in the gl.nhybrids function of dartR.

Adult morphology.—We examined preserved specimens held in the AMS collection, including the type series for Limnodynastes dorsalis var. terraereginae, and all genotyped specimens in the SAMA collection. In addition, we examined the type of Platyplectrum superciliare at Zoologisches Forschungsmuseum Alexander Koenig, Bonn, Germany (ZFMK) and high-resolution images of the type of Heliorana grayi held at the Natural History Museum Vienna, Austria (NHMW). We measured 23 morphometric characters (adapted from Watters et al., 2016) with digital calipers to the nearest 0.1 mm (Table 2) and documented variation in a range of external morphological features deemed by preliminary investigation to be potentially diagnostic (i.e., extent of webbing on the hind foot, ventral pattern, and presence/absence of a vertebral stripe), thereby extending our analysis of morphological traits to include poorly preserved specimens for which we were unable to obtain reliable morphometric measurements. Abbreviations for all morphometric traits are listed in Table 2. Sex of adults was determined by directly observing testes and darkened nuptial pads or vocal sacs in males, or ovarian follicles and distended finger spatulae in females.

Table 2

Definition of morphometric traits measured for the Limnodynastes dorsalis group. Asterisks indicate characters not defined by Watters et al. (2016).

To compare differences in geometric shape among taxa, we used a multivariate linear discriminant function analysis (DFA). Due to sexual size dimorphism, male and female samples were analyzed independently. Potentially confounding variation associated with differing body sizes and allometric growth was minimized by adjusting measurements to the values they would assume if they were of a mean body size for that sex using the allometric growth equation of Thorpe (1976): Y_i* = log₁₀Y_i – b(log₁₀SVL_i – log₁₀SVL_mean), where Y_i* is the adjusted value for character Y of the ith specimen; Y_i is the raw/unadjusted value for character Y; b is the mean of the regression coefficients for Y_i against SVL_i estimated independently for each taxon from logarithmically transformed values of Y_i and SVL_i; SVL_i is the measured snout–vent length (SVL) of the ith specimen; and SVL_mean is the pooled mean SVL.

DFAs were conducted after the measurements had been adjusted for size/growth and log-transformed as described above using the lda function from v7.3-40 of the R package MASS (Venables and Ripley, 2002) in RStudio v 4.2.1. We allocated specimens a priori into three “taxa” for the DFA based on individual genotypes or, if genetic data were not available, then based on whether the collection location fell within the geographic distribution of a genetic group. The raw mensural data and DFA results are presented in Supplementary Table S1 (see Data Accessibility). Because the geographical provenance of the holotype of Platyplectrum superciliare is uncertain (=Australia), we assigned it to a taxon using the predict function for unknowns in lda.

Advertisement calls.—In-field recordings were made using a Zoom H5 Handy Recorder with a RODE NTG2 shotgun microphone. Acoustic recordings obtained from the FrogID dataset were 20–60 second recordings (MPEG AAC audio file) made using smartphones, with a sampling rate of 44.1 kHz.

Acoustic analyses were conducted using Raven Pro 1.6 (Center for Conservation Bioacoustics, 2019), with a fast-Fourier transformation of 512 points, and 50% overlap. We measured dominant frequency (kHz), call duration (s), and call rate (calls/min) for up to five calls per individual frog, where call rate was calculated using the formula: (e.g., Mitchell et al., 2020). To Sample duration (s) visualize the harmonics of each call more accurately, we then used a fast-Fourier transformation of 1,024 points, and 50% overlap to measure fundamental frequency (kHz). Following the definitions of Köhler et al. (2017), dominant frequency refers to the frequency of a call which contains the highest energy, while fundamental frequency refers to the base frequency of a call. We obtained mean values for each call parameter, where the unit of replication was the individual.

Ambient temperature was recorded in the field for the four acoustic recordings taken during fieldwork. As ambient temperature is not recorded in the FrogID app (Rowley et al., 2019), we estimated temperature for all FrogID recordings using historical weather data from the Australian Bureau of Meteorology (BOM). Using the package “chillR” (Luedeling, 2021), we obtained an hourly estimate of temperature based on the minimum and maximum temperatures of the 15 days before and 15 days after each recording was made (e.g., Mitchell et al., 2020). These temperature data were estimated from the BOM weather station closest to each recording. To determine whether the interaction between call parameters and taxonomic lineage was influenced by ambient temperature, we conducted linear regression models.

We assigned taxonomic group based on location and analyzed recordings from genotyped specimens and throughout the core ranges of each lineage. To determine whether call parameters differ among the lineages, we conducted one-way analysis of variance. Given the significant relationship between temperature and fundamental frequency, we conducted a one-way analysis of covariance for fundamental frequency, with ambient temperature as a factor. To test whether differences existed between groups, we conducted pairwise comparisons using Tukey post hoc tests of honest significant difference. All statistical analyses were conducted in R (R Core Team, 2020).

Conservation assessments.—Due to an absence of population data, we assessed the conservation status of the lineages redescribed herein using Criterion B of the IUCN Red List guidelines (IUCN Standards and Petitions Committee, 2022) by estimating the geographic range of each taxon through calculation of their area of occupancy (AOO) and extent of occurrence (EOO). Under the IUCN Red List guidelines, to qualify as threatened under Criterion B, taxa must not only meet the minimum distribution threshold (AOO: <2,000 km², EOO: <20,000 km²) but also at least two of three other conditions, specifically: (a) severely fragmented or number of locations ≤10; (b) continuing decline observed, estimated, inferred, or projected in any of: (i) EOO, (ii) AOO, (iii) area extent and/or quality of habitat, (iv) number of locations or subpopulations, (v) number of mature individuals; or (c) extreme fluctuations in any of: (i) EOO, (ii) AOO, (iii) number of locations or subpopulations, (iv) number of mature individuals. Calculations of AOO and EOO were made using the Atlas of Living Australia (ALA) spatial portal ( https://www.ala.org.au; accessed online August 2022). For the assessments, we estimated the distribution of each taxon by tracing polygons in QGIS (v. 3.10.9; QGIS.org, 2023) around occurrence records obtained from genotyped specimens, morphologically examined museum specimens, validated acoustic recordings from FrogID, and occurrence records obtained from ALA (accessed August 2022). These distribution polygons were then uploaded into the ALA spatial portal and used as the basis for assessment of AOO (2 × 2 km grid resolution) and EOO (minimum convex hull).

RESULTS

Mitochondrial DNA.—The sequences in the ND4 alignment varied from 408 to 715 bp in length. The best-fit nucleotide-substitution model identified was TIM3+F+I+G4. The maximum-likelihood phylogenetic analyses (Fig. 2) provided strong support for two major clades within Limnodynastes terraereginae sensu lato, including a lineage restricted to Cape York Peninsula in Far North QLD (hereafter referred to as the CYP lineage) and a widespread lineage (hereafter the WS lineage) covering the remainder of the species range. The WS lineage was further subdivided into northern (mid to north QLD) and southern sub-clades; however, the southern sub-clade did not receive strong bootstrap support. Samples of L. dumerilii grayi from eastern coastal NSW (hereafter the EC lineage) formed a highly distinct and well-supported sister lineage to the CYP and WS L. terraereginae clades. Limnodynastes interioris and L. dumerilii (including the subspecies dumerilii, fryi, insularis, and variegatus) form a divergent ancestral clade which is closest to the EC lineage, wherein L. interioris is the sister sub-clade to the L. dumerilii subspecies dumerilii and variegatus, rendering L. dumerilii paraphyletic.

Fig. 2

Maximum likelihood phylogenetic analyses of the Limnodynastes dorsalis group inferred using mitochondrial ND4 sequences. Tree is rooted with Limnodynastes fletcheri, L. peronii, L. depressus, L. tasmaniensis, and L. salmini. Dots at nodes indicate bootstrap support values: >95% black dots; 70–95% gray dots. Scale bar represents substitutions/nucleotide site. See Data Accessibility for tree file.

Net average sequence divergence between the CYP, WS, and EC lineages ranged from 4–9% (Table 3). The EC lineage was most divergent, with net sequence divergence of 7% from the WS lineage and 9% from the CYP lineage. Net sequence divergence of the CYP lineage from the northern and southern sub-clades of the WS lineage ranged from 4–5%, with the northern and southern sub-clades only 2% divergent from each other. The CYP lineage is diagnosed by 16 apomorphic nucleotide states, WS lineage by 5 diagnostic states, and EC by 17 diagnostic states (Table 4).

Table 3

Net average mitochondrial ND4 sequence divergence between species of Limnodynastes.

Table 4

Comparison of 36 apomorphic diagnostic nucleotide states (in bold) for the CYP, WS, and EC lineages compared with other members of the Limnodynastes dorsalis group. Site numbers (vertical) reflect position of the site within the ND4 sequence alignment analyzed.

SNP analyses.—A total of 12,179 binary SNP loci were scored for 125 individual samples. After filtering by a call rate of 95% and locus presence in at least 70% of individuals, 1,368 SNP loci for 121 individuals remained, with a total of 2.2% missing data across the dataset. Four samples (A005543, EBU115629, EBU115671, and R171924) were removed due to not meeting the call-rate threshold.

In the initial PCoA clustering analysis involving the total SNP dataset, the proportion of variation explained by the PC axes were: 1^st axis—40%, 2^nd axis—20%, 3^rd axis—10%, and 4^th axis—4%. Three distinct genetic clusters were present in the PCoA, including: (1) EC lineage; (2) WS/CYP lineages; and (3) L. interioris and the L. dumerilii subspecies cluster (Fig. 3A). One sample (ABTC142322) appeared to be intermediate between the EC lineage and L. dumerilii and was later confirmed by STRUCTURE and NewHybrids to be a backcrossed hybrid. To further examine relationships within the lineages of L. terraereginae, we performed a secondary PCoA clustering analysis comparing only the WS (n = 51) and CYP (n = 3) samples. In this analysis, northern and southern samples of the WS lineage clustered across a continuum which appeared to correspond with geographic distribution (possibly indicating a pattern of IBD), whereas the CYP samples were clearly separated from the WS lineage on both the 1^st and 2^nd axes of the PCoA (Fig. 3B). To further examine the influence of geographic distance on divergence within the WS lineage, we performed an independent IBD analysis for this group which confirmed a significant (P value = 0.001) linear correlation between geographic and genetic distances matrices for individuals. This suggests that IBD is an important driver of the divergence observed within the WS lineage.

Fig. 3

Results of clustering analyses of SNP data for the eastern Limnodynastes dorsalis group: (A) initial PCoA ordination with representative samples of all eastern Limnodynastes dorsalis group members, (B) PCoA of CYP and WS lineages only, (C) STRUCTURE barplot for a balanced subset of individuals showing preference for K = 3 ancestry model.

The STRUCTURE analyses of the balanced subset including the CYP, WS, and EC lineages identified K = 3 as the optimal ancestry model. The plot of individual ancestry coefficients, based on K = 3, clustered northern and southern individuals of the WS lineage as a single clade, with the CYP and EC lineages clustering as separate distinct clades (Fig. 3C). The number of loci showing fixed allelic differences between the WS, CYP, and EC lineages ranged from 23–131 (out of 1,574 loci compared; Table 5) and comparisons between each taxa were significant after simulation (P < 0.001).

Table 5

Numbers of loci with fixed allelic differences between lineages. Above diagonal: number of loci showing a fixed difference; below diagonal: expected count of false positives for each comparison, by simulation. Pairwise comparisons were based on 1,368 loci for WS and EC lineages and 1,327 for the CYP lineage. Values in bold were significant after simulation. Values in parentheses denote sample size per clade included in the pairwise comparisons.

In the neighbor-joining (NJ) analysis of relationships among individuals (Fig. 4), WS, CYP, and EC samples each formed distinct lineages consistent with the SNP clustering findings (Fig. 3) and the relationships among their mtDNA sequences (Fig. 2). In both the mtDNA phylogenetic analysis and the SNP NJ tree, WS and CYP are sister lineages. In the SNP NJ analysis, L. interioris is the sister lineage to all the subspecies of L. dumerilii which contrasts with the relationships among their mitochondrial DNA sequences wherein L. interioris is the sister lineage to only the dumerilii and variegatus subspecies of L. dumerilii.

Fig. 4

Neighbor-joining tree based on analyses of the SNP dataset for the eastern Limnodynastes dorsalis group. Scale bar represents substitutions/site. See Data Accessibility for tree file.

Admixture analyses.—Our surveys confirmed that the eastern species of the Limnodynastes dorsalis group variously occupy allopatric, parapatric, and sympatric distributions throughout their range, with boundaries between species corresponding with shifts in habitat and/or major biogeographic boundaries. We found the WS and CYP lineages to be allopatric over a distance of at least 110 km between the Einasleigh Uplands/Wet Tropics and Cape York Peninsula bioregions in northern QLD with no evidence of admixture between them.

Similarly, the WS and EC lineages appear to be narrowly allopatric in northern NSW. Analyses of the WS–EC subset did not detect admixture, despite inclusion of samples collected within ∼70 km between South West Rocks (EC lineage) and Coffs Harbour (WS lineage). Morphologically verified museum specimens indicate that the EC lineage may occur as far north as the Nambucca River, which would reduce the gap between lineages to ∼49 km. However, numerous targeted surveys conducted in 2020–2022 failed to locate either species in closer proximity than the samples available to us and so we consider the taxa to be allopatric at present.

In contrast, the WS lineage and L. dumerilii were found to be parapatric/sympatric in a variety of ecotonal habitats in northern and western NSW. We did not detect admixture between the subsets of WS–L. dumerilii, despite sampling of the lineages from close geographic proximity in north-eastern NSW (∼23 km), Northern Tablelands (∼15 km), and Central Tablelands region (∼40 km).

Admixture was detected between the EC lineage and L. dumerilii in samples from an ecotone between the species in the Central Coast region of NSW (Fig. 5). The NewHybrids analyses classified ABTC142322 with 95% posterior probability for being a backcrossed hybrid (Supplementary Table S2; see Data Accessibility). This was the same individual identified in the PCoA and neighbor-joining tree (Fig. 4) as intermediate between the EC lineage and L. dumerilii. In addition, the STRUCTURE analyses classified the same individual as admixed and indicated potential admixture in a further six samples from the same region (ABTC140573, ABTC140438, ABTC1175, ABTC25842, ABTC140441, and ABTC140571; Fig. 5). The overall level of hybridization occurring in this area appears to be relatively low given no F₁ or F₂ hybrids were detected, and an additional seven samples collected at the same site or within 32 km of the admixed individuals were classified as pure EC lineage. We also found no further evidence of admixture between the EC lineage and L. dumerilii in other regions of NSW where the taxa have been suspected previously to hybridize (Martin, 1972). This included samples of each lineage collected within a straight-line distance of 40 km in southern Sydney, and 17 km in the Blue Mountains across an elevational transition of 404 m.

Fig. 5

Map of sampled contact zones between the EC lineage (green) and L. dumerilii (gray), with STRUCTURE ancestry proportions for individuals indicated in pie charts. Labeled polygons are IBRA7 major biogeographic subregions.

Adult morphology.—There was substantial variation in overall size among the taxa (Fig. 6; Table 6), with the CYP lineage tending to be largest, EC lineage the smallest, and the WS lineage generally falling into a size class between the two (Fig. 6). We therefore relied on analyses of the SVL-corrected morphological measurements to detect variation in overall geometric shape among the taxa. The DFAs for males and females each returned two linear discriminant functions (LD), with the majority of variation explained by LD1 (88% for males and 74% for females; Fig. 7). For males (n = 94), the overall predictive accuracy was 0.88, while the overall predictive accuracy for females (n = 51) was 0.92. For males, the traits with the highest coefficients for each of the two linear discriminants were for LD1: PL, SL, and Toe4W and for LD2: SVL and TEY. For females, the traits with the highest coefficients for each of the two linear discriminants were for LD1: UAL, AL, and TEY, and for LD2: HDD and FOL. Finally, the holotype for Platyplectrum superciliare (ZFMK 28331) was assigned by the DFA to the EC lineage morphogroup with a posterior probability of 0.99.

Fig. 6

Comparison of morphological measurements between adult specimens of the Limnodynastes dorsalis group lineages redescribed herein.

Table 6

Summary of metric variation (mean ± SD, and range in mm) for 23 morphometric traits between members of the Limnodynastes dorsalis group lineages redescribed herein. Sample size in parentheses beside sex (m or f).

Fig. 7

Linear discriminant function analysis (DFA) scatterplot of adult morphometric characters for specimens of the Limnodynastes dorsalis group lineages redescribed herein.

Advertisement calls.—We analyzed the advertisement calls of 65 individual male banjo frogs (Table 7): WS lineage (n = 41), EC lineage (n = 19), and CYP lineage (n = 5). Ambient temperature was not significantly correlated with dominant frequency (R² = 0.30, P = 0.31), call duration (R² = 0.10, P = 0.72), or call rate (R² = 0.01, P = 0.33), although a positive relationship existed with fundamental frequency (R² = 0.22, P = 0.02). The advertisement calls of each species were similar, with considerable overlap in the measured call parameters (Figs. 8, 9; Table 7). While dominant frequency (F_2,62 = 14.44, P < 0.001), call duration (F_2,62 = 5.78, P = 0.005), and fundamental frequency (F_2,62 = 7.00, P = 0.02) differed significantly among taxa, call rate did not differ significantly (F_2,62 = 1.17, P = 0.331).

Table 7

Summary of advertisement call data for members of the Limnodynastes dorsalis group lineages redescribed herein, expressed as means ± standard deviation and the range of values for each taxon.

Fig. 8

Comparison of the male advertisement calls of (A) WS lineage—Rockhampton, QLD (FrogID 312613), (B) EC lineage—Harrington, NSW (AMS R.188164), (C) CYP lineage—Finch Bay, Cooktown, QLD (QM J97862), (D) Limnodynastes dumerilii dumerilii—Tenterfield Creek, NSW (AMS R.188096), and (E) L. interioris —Narrandera, NSW (FrogID 276467). Audiospectrograms and oscillograms were produced in Raven Pro 1.6 with a fast-Fourier transformation of 512 points, and 50% overlap.

Fig. 9

Boxplot comparison of variation in the advertisement call traits between the Limnodynastes dorsalis group lineages redescribed herein, depicting variation in (A) dominant frequency (kHz), (B) call duration (s), (C) fundamental frequency (kHz), and (D) call rate (calls/min). Small dots represent outliers.

The advertisement call of the EC lineage was most distinct from the other taxa (Figs. 8, 9; Table 7), with both a significantly higher dominant frequency and fundamental frequency than the WS lineage (dominant frequency—P < 0.001; fundamental frequency—P = 0.003) and CYP lineage (dominant frequency—P = 0.003; fundamental frequency—P = 0.030). While both the dominant and fundamental frequencies of CYP lineage calls were marginally lower than those of the WS lineage (dominant frequency—P = 0.784; fundamental frequency—P = 0.741), these differences were not significant. Calls of the EC lineage were shorter in average duration than both the WS lineage (P = 0.025) and CYP lineage (P = 0.012). The call durations of the CYP lineage were only marginally longer than the WS lineage (P = 0.273). We observed no differences between the advertisement calls of the northern and southern mitochondrial sub-clades of the WS lineage, with respect to dominant frequency (P = 0.909), call duration (P = 1.000), and fundamental frequency (P = 0.556), although the northern clade did call at a faster rate than the southern clade (P = 0.007).

Systematic implications.—Based on broad congruence between our datasets, we conclude that the WS, CYP, and EC lineages warrant recognition as distinct species under the evolutionary species concept (sensu de Queiroz, 1998, 2007). Our evidence for lineage separation is based on the following operational criteria:

Monophyly: Our analyses of sequences of the mitochondrial ND4 gene reveal the presence of three, well-supported, reciprocally monophyletic mitochondrial groups (Fig. 2), with a level of sequence divergence of 4–11% between lineages (Table 3), comparable to other limnodynastid and myobatrachid species groups such as Heleioporus (4–22%, Mahony et al., 2021b), Philoria (5–15%, Mahony et al., 2022), and Assa (6%, Mahony et al., 2021a). Our analyses of the SNP dataset revealed consistent clustering of samples in both the PCoA ordination, neighbor-joining tree, and Bayesian methods (Figs. 3, 4), with the SNP clusters corresponding fully to the lineages observed in the ND4 phylogeny.

Reproductive isolation: Representative sampling across the range of each proposed taxon identified deep genetic breaks among the lineages across relatively short geographic distances. Interactions among the taxa can be variously characterized as allopatric, parapatric, and sympatric. Low levels of admixture were detected at an ecotone between the parapatric EC lineage and Limnodynastes dumerilii (Fig. 5), suggesting occasional hybridization between these non-sister species occurs. We suspect hybridization occurs due to accidental mismating following the sporadic migration of individuals between otherwise discrete habitats.

Diagnosability: The lineages are diagnosable based on a combination of: a) a range of external morphological features (i.e., adult body size, extent of foot webbing [Fig. 10] and aspects of color/pattern); b) differences in geometric shape as identified in the linear discriminant function analyses (Fig. 7); c) differences in the male advertisement call (i.e., dominant frequency, call duration), particularly between the EC lineage and other taxa; and d) 5–17 apomorphic nucleotide character states (Table 4) and fixed allelic differences at significant numbers of SNP loci (Table 5).

Fig. 10

Illustration of hind foot webbing extent, a useful diagnostic character for identification of the Limnodynastes dorsalis group. (A) Well-developed webbing (extends to beyond half-way between the 1^st subarticular tubercle and tip of Toe 1), (B) moderate webbing (extends one-quarter to half-way between the 1^st subarticular tubercle on Toe 1), (C) vestigial webbing (does not or only slightly extends beyond the 1^st subarticular tubercle on the 1^st Toe). Illustration copyright Alana de Laive.

Nomenclatural implications.—To resolve long-standing confusion surrounding the provenance and identity of the holotypes of Heliorana grayi and Platyplectrum superciliare we examined the type specimens and GMS investigated their provenance. We briefly discuss the nomenclatural implications of our findings below and provide a more detailed summary of the provenance of the type specimens in Appendix 1.

Provenance of the Heliorana grayi type: The collection locality for the Heliorana grayi type (NHMW 4695) was initially vaguely specified by Steindachner (1867) as Neu-Süd-Wales (i.e., NSW). Inspection of the NHMW type catalogue and an investigation into the provenance of the type specimen (see Appendix 1) has since revealed the true collection locality for the type of Heliorana grayi to be Rockhampton, QLD (Häupl and Tiedemann, 1978; Gemel et al., 2019).

Without examining the type, authors of subsequent taxonomic revisions (i.e., Parker, 1940; Martin, 1972) were led to assume that the type of Heliorana grayi was collected in NSW and thus incorrectly referred the name Limnodynastes (dorsalis) dumerilii grayi to the distinctive population of frogs occurring in eastern coastal NSW (here referred to as the EC lineage). Our examination of high-resolution images of the type of Heliorana grayi confirmed that the specimen corresponds in morphology (medium-large size, robust build, moderate foot webbing, and aspects of dorsal color/pattern) with genotyped specimens of the WS lineage which occurs in Rockhampton, and we therefore correctly apply the name Limnodynastes grayi to this taxon.

The original collection locality for the type of Platyplectrum superciliare (ZFMK 28331) was even more vaguely stated by Keferstein (1867) as Australien (i.e., Australia). Based on the investigation into the provenance of the type specimen (see Appendix 1), the results of the group assignment by the DFA, and the presence of several consistent diagnostic morphological features (i.e., small size, vestigial foot webbing [Fig. 10] and aspects of dorsal and ventral color/pattern [Fig. 11]), we conclude that the type specimen of Platyplectrum superciliare represents the distinctive EC lineage currently incorrectly referred to as Limnodynastes dumerilii grayi. We hereby apply the name Limnodynastes superciliaris to this taxon.

Fig. 11

Examples of typical ventral color pattern, a useful diagnostic trait between members of the eastern Limnodynastes dorsalis group. (A) Limnodynastes terraereginae (CYP lineage), (B) L. grayi (WS lineage), (C) L. superciliaris (EC lineage), (D) L. dumerilii dumerilii, (E) L. dumerilii insularis, and (F) L. interioris. Specimen registration labels represent 25 mm.

Notes on the Limnodynastes dorsalis var. terrae-reginae type series: There is no conjecture over the original collection locality of the Limnodynastes dorsalis var. terrae-reginae type (AMS R.4525, Somerset, Cape York Peninsula, QLD). Our examination of the holotype confirmed it corresponds in morphology (large size, moderate foot webbing [Fig. 10]; aspects of color/pattern [Fig. 11]) with genotyped individuals of the CYP lineage and so we apply the name Limnodynastes terraereginae to the distinctive lineage of banjo frogs restricted to the Cape York Peninsula bioregion. Of an additional 57 paratypes held in the Australian Museum, we consider only one (AMS R.4526, Somerset, Cape York Peninsula, QLD) represents true L. terraereginae. The remaining paratypes, collected around Eidsvold in the Upper Burnett River region of mid-eastern QLD, can be assigned to L. grayi on the basis of size, aspects of color/ pattern, and proximity to genotyped samples.

Taxonomy.—An updated diagnostic key to the species of the eastern Australian Limnodynastes dorsalis group is provided below.

Fig. 12

Holotype of Heliorana grayi, NHMW 4695. (A) Left-side profile, (B) hind foot, (C) dorsal profile, (D) ventral profile. Images copyright Alice Schumacher, NHM Vienna.

Fig. 13

Images in life of Limnodynastes grayi. (A) AMS R.188151, female, Dirty Creek, north coast, NSW. (B) AMS R.188350, female, Tyndale, north coast, NSW. (C) AMS R.185843, male, Yetman, Northern Tablelands, NSW. (D) QM J97851, male, Watsonville, Atherton Tablelands, Qld. (E) QM J97848, male, Hervey Range, Qld. (F) QM J97857, female, Mount Carbine, Qld.

Fig. 14

Holotype of Platyplectrum superciliare, ZFMK 28331. (A) Left-side profile, (B) hind foot, (C) dorsal profile, (D) ventral profile. Images copyright Morris Flecks, ZFMK, Germany.

Fig. 15

Images in life of Limnodynastes superciliaris. (A) AMS R.188161, male, South West Rocks, NSW. (B) AMS R.188207, male, lower Blue Mountains, NSW. (C) AMS R.188165, male, Agnes Banks, Cumberland Plains, NSW. (D) AMS R.188164, Harrington, mid-northern coast, NSW. (E) AMS R.188163, female, Laurieton, mid-northern coast, NSW. (F) AMS R.188135, juvenile, Hat Head National Park, mid-northern coast, NSW.

Fig. 16

Holotype of Limnodynastes dorsalis var. terrae-reginae, AMS R4525. (A) Right-side profile, (B) hind foot, (C) dorsal profile, (D) ventral profile.

Fig. 17

Images in life of Limnodynastes terraereginae. (A) QM J97862, female, Cooktown, Qld. (B) QM J97860, male, Cooktown, Qld. (C) QM J97861, male, Cooktown, Qld. (D) Unvouchered individual, Kutini-Payamu (Iron Range) National Park, Qld (photo copyright: Stephen Zozaya).

DISCUSSION

Our study represents the first range-wide assessment of variation in genetic diversity, morphology, and calls of Limnodynastes terraereginae sensu lato. Our results provide strong evidence for the recognition of three evolutionarily distinct taxa within this species group. These include the massive and superbly patterned L. terraereginae that is endemic to Cape York Peninsula, the medium-sized and widely distributed L. grayi, and a diminutive sister taxon, L. superciliaris (formerly L. dumerilii grayi), which is endemic to sandy habitats of eastern coastal NSW.

We found that the species described herein occupy narrowly allopatric ranges which replace each other in a north–south series over a distance of more than 2,800 km, extending from the Sydney Basin to the northernmost tip of Cape York Peninsula. The distributional gap between L. grayi and L. terraereginae occurs across a distance of 110 km between the Carbine Uplands and Cooktown in southern Cape York Peninsula. The intervening region is characterized as a lowland expanse of hot and humid open woodland and savannah which appears to present a barrier to dispersal between the species today. In tropical north QLD, L. grayi is restricted to upland sclerophyll woodlands, open forests, and heathland above 700 m a.s.l., and the data suggest that it is absent from hot and humid coastal lowlands, suggesting a physiological intolerance to elevated temperature and humidity. In contrast, L. terraereginae appears to be well adapted to these conditions given its distribution is entirely restricted to eastern lowland coastal habitats within Cape York Peninsula.

In northern NSW, L. grayi and L. superciliaris appear to be separated by a gap of approximately 70 km between Coffs Harbour (L. grayi) and South West Rocks (L. superciliaris). We conducted several targeted surveys within this gap and failed to locate either species within closer proximity. Of note, this region also represents a southerly distributional limit for the Wallum Sedge Frog (Litoria olongburensis) and gap in the distribution for the Wallum Froglet (Crinia tinnula), suggesting that dispersal across this zone may be restricted by an absence of suitable wallum habitat.

In contrast, L. grayi and L. superciliaris were found to occur in sympatry and/or parapatry with L. dumerilii in various locations throughout the species' range. The distributions of L. dumerilii and L. grayi overlap broadly in north-eastern NSW and south-eastern QLD where they have been found breeding side-by-side without evidence of hybridization (i.e., no morphological intermediates; Martin, 1972). The advertisement calls of each taxon are markedly distinct; however, whether this character displacement has reinforced species boundaries in sympatry remains unclear. Martin (1972) found the dominant frequency of the call of L. dumerilii in sympatry with L. grayi is significantly lower (i.e., more different from L. grayi) than those of allopatric populations. However, the opposite trend was seen for L. grayi, which has calls more similar to L. dumerilii in sympatry compared to allopatric populations. We found no evidence of admixture between these taxa in our SNP dataset despite sampling from close geographic proximity in north-eastern NSW (∼23 km), Northern Tablelands (∼15 km), and Central Tablelands region (∼40 km), affirming the species maintain reproductive isolation in parapatry/sympatry.

In NSW, L. superciliaris is predominantly distributed along the east coast, with L. dumerilii occurring mainly within and west of the GDR. The presence of possible morphological intermediates between the species in southern Sydney (i.e., around Wattamolla and Stanwell Tops) previously has raised the possibility of a hybrid zone among the taxa, although insufficient data were available to conclude on the extent of hybridization (Martin, 1972). We detected admixture between L. dumerilii and L. superciliaris at an ecotonal contact zone on the Central Coast of NSW, at the junction of the Pittwater, Yengo, and Wyong sub-bioregions (Fig. 5). The overall level of admixture among the taxa was found to be relatively low given no F₁ or F₂ hybrids were detected, and an additional seven samples collected at the same site, or within 32 km of the admixed individuals, were classified as pure L. superciliaris. We also found no evidence of admixture between L. dumerilii and L. superciliaris in other regions of potential contact. This included samples collected within a straight-line distance of 40 km in southern Sydney, and 17 km in the Blue Mountains across an altitudinal transition of 404 m. We found the taxa to occupy distinct habitats: L. superciliaris occurs within lowland sandy heaths and L. dumerilii is associated with wet and dry sclerophyll forest, agricultural, and suburban areas on heavier soils. There is marked differentiation in morphology and advertisement calls between the taxa which likely acts as a pre-mating reproductive isolation barrier. Based on these findings, we suspect that occasional hybridization occurs as a consequence of mismating following sporadic migration of individuals between otherwise fairly discrete habitats. Further studies of the Central Coast contact zone would shed light on the reproductive compatibility of the taxa and the role of pre- and post-mating barriers in maintaining contemporary species boundaries.

DATA ACCESSIBILITY

Supplemental material is available at  https://www.ichthyologyandherpetology.org/h2023025. Unless an alternative copyright or statement noting that a figure is reprinted from a previous source is noted in a figure caption, the published images and illustrations in this article are licensed by the American Society of Ichthyologists and Herpetologists for use if the use includes a citation to the original source (American Society of Ichthyologists and Herpetologists, the DOI of the Ichthyology & Herpetology article, and any individual image credits listed in the figure caption) in accordance with the Creative Commons Attribution CC BY License.