We follow the unified species concept, considering a species as an independent evolutionary lineage with consistent diagnostic traits (e.g., morphological, behavioral, genetic; de Queiroz 2005, 2007). Species delimitation based solely on genetic data is likely inadequate (Carstens et al. 2013). To address this, taxonomic decisions were made based on a combination of information sources, including mitochondrial, morphological, and acoustic data, which represent independent lines of evidence.

To assess the taxonomic status of Chaco populations of L. fuscus that could be assignable to L. gualambensis and its phylogenetic position within Leptodactylus, we included representatives of all species groups recognized for the genus by de Sá et al. (2014), with special emphasis on the L. fuscus species group. Considering the broad distribution of L. fuscus and its recognition as a species complex, we included sequences of specimens from several localities covering the geographic distribution of the species, including the type locality. Given our aim of testing the taxonomic validity of the Chaco populations of L. fuscus, we focused efforts on sampling from the South American Gran Chaco and its surroundings (including samples from the L. gualambensis type locality), following the geographic distribution of the species suggested by Gallardo (1964) (Fig. 1; File S1). Representation of other taxa from the L. fuscus group included samples from Leptodactylus apepyta Schneider et al., 2019, Leptodactylus avivoca Carvalho et al., 2021, Leptodactylus barrioi Silva et al., 2020, Leptodactylus bufonius Boulenger, 1894, Leptodactylus camaquara Sazima & Bokermann, 1978, Leptodactylus cunicularius Sazima & Bokermann, 1978, Leptodactylus cupreus Caramaschi, Feio & São Pedro, 2008, Leptodactylus didymus Heyer, García-Lopez & Cardoso, 1996, Leptodactylus elenae Heyer, 1978, Leptodactylus fragilis (Brocchi, 1877), Leptodactylus furnarius Sazima & Bokermann, 1978, L. gracilis, Leptodactylus jolyi Sazima & Bokermann, 1978, Leptodactylus kilombo Silva et al., 2020, Leptodactylus labrosus Jiménez de la Espada, 1875, Leptodactylus laticeps Boulenger, 1918, Leptodactylus latinasus Jiménez de la Espada, 1875, Leptodactylus longirostris Boulenger, 1882, Leptodactylus marambaiae Izecksohn, 1976, Leptodactylus mystaceus (Spix, 1824), L. mystacinus, Leptodactylus notoaktites Heyer, 1978, Leptodactylus oreomantis Carvalho, Leite & Pezzuti, 2013, Leptodactylus plaumanni Ahl, 1936, Leptodactylus poecilochilus (Cope, 1862), Leptodactylus syphax Bokermann, 1969, Leptodactylus tapiti Sazima & Bokermann, 1978, Leptodactylus troglodytes Lutz, 1926, Leptodactylus ventrimaculatus Boulenger, 1902, and Leptodactylus watu Silva et al., 2020. From the other groups of Leptodactylus we included: Leptodactylus insularum Barbour, 1906 and Leptodactylus latrans (Steffen, 1815) from the Leptodactylus latrans group; Leptodactylus colombiensis Heyer, 1994 and Leptodactylus podicipinus (Cope, 1862) from the Leptodactylus melanonotus group; Leptodactylus labyrinthicus (Spix, 1824) and Leptodactylus pentadactylus (Laurenti, 1768) from the Leptodactylus pentadactylus group. External to Leptodactylus we included samples of Adenomera thomei (Almeida & Angulo, 2006) and Lithodytes lineatus (Schneider, 1799) (Leptodactylinae); Physalaemus cuvieri Fitzinger, 1826 and Engystomops freibergi (Donoso-Barros, 1969) (Leiuperinae), Paratelmatobius poecilogaster Giaretta & Castanho, 1990, Crossodactylodes septentrionalis Teixeira et al., 2013, and Scythrophrys sawayae (Cochran, 1953) (Paratelmatobiinae), and to root the trees, Hyalinobatrachium fleischmanni (Boettger, 1893).

We extracted total genomic DNA from samples conserved in 100% ethanol (muscle or liver) using the DNeasy extraction kit (Qiagen, Valencia, CA, USA) following the manufacturer’s protocol. We used gene fragments of the 16S ribosomal RNA gene (16S; up to 580 base pairs [bp] per terminal) and of the cytochrome c oxidase subunit I (COI; up to 597 bp per terminal), two markers commonly used in taxonomic approaches (Vences et al. 2005) and proposed as barcoding fragments (see Lyra et al. 2017). For 16S we used the primers 16S-AR (CGCCTGTTTATCAAAAACAT; Palumbi et al. 1991) and 16S-Wilk2 (GACCTGGATTACTCCGGTCTGA; Wilkinson et al. 1996), and for COI, AnF1 (ACHAAYCAYAAAGAYATYGG; Lyra et al. 2017) and AnR1 (CCGGTCTGAACTCAGATCACGT, Lyra et al. 2017). For fragment amplifications through polymerase chain reactions (PCR) we used a commercial kit (Master Mix, Thermo Fisher; Walthan, MA, USA). For the 16S fragment, the protocol included an initial denaturation step of 3 min at 94°C; followed by 35 cycles consisting in 30 s at 94°C for denaturation, 30 s at 50°C for annealing, and 90 s at 72°C for extension; the final step consisted in 7 min at 72°C for extension. For the COI fragment, the protocol consisted in an initial denaturation step of 3 min at 95°C; followed by ten cycles consisting in 20 s at 95°C for denaturation, 20 s at 50°C (increasing +0.5°C in each subsequent cycle) for annealing, 50 s at 60°C for extension; the next step consisted in 25 cycles including 15 s at 95°C for denaturation, 20 s at 50°C for annealing, 50 s at 60°C for extension; the final step consisted in 5 min at 60°C for extension (Lyra et al. 2017). We sent PCR products for purification and sequencing to Macrogen Inc. (Seoul, South Korea). We checked chromatograms and edited the sequences in CodonCode Aligner 3.5.4 (Codon Code Corporation). See File S1 for GenBank accession numbers.

We aligned both 16S and COI sequences with MAFFT (Katoh and Toh 2008). The 16S sequences were aligned with the Q-INS-i strategy (which considers RNA secondary structure) and the COI sequences with the G-INS-i strategy (which considers global homology). All other parameters were set as default. Sequence files were merged with SequenceMatrix (Vaidya et al. 2011).

To reconstruct phylogenetic trees, we used Bayesian inference (BI) and maximum parsimony (MP) as criteria. While parsimony allows the usage of gaps as a fifth state of character (i.e., it considers the value of indels in the evolutionary process; Giribet and Wheeler 1999; Goloboff and Pol 2005), the Bayesian approach permits the implementation of distinct evolutionary models for distinct DNA fragments, acknowledging that these fragments may be subject to different evolutionary pressures (Ronquist and Deans 2010). To construct trees under BI, we first used PartitionFinder v2.1.1 (Lanfear et al. 2016) to infer the best partition scheme and the evolutionary model that best fits each partition. Branch lengths were treated as linked, and the analysis considered only the models employed by MrBayes 3.2 (Ronquist et al. 2012). We analyzed 16S and codon positions of COI as distinct partitions. The best models were selected under the corrected Akaike information criterion (AICc) using the “greedy” heuristic search algorithm and PhyML software (Guindon et al. 2010). We performed the BI analysis in MrBayes v3.2, implementing the inferred models of nucleotide substitution on two independent runs, each one with four chains, sampling every 3,000 generations for 30 million generations and discarding the first 2,500 trees as burn-in. We verified convergence with TRACER v1.6 (Rambaut et al. 2014), and by examining the standard deviation of split frequencies between independent runs (<0.01), Estimated Sample Size (values >100), and Potential Scale Reduction Factor (values ~ 1). Node support was estimated as posterior probability (PP).

For MP analyses, we used T.N.T. Willi Hennig Society Edition (Goloboff et al. 2003, 2008; Goloboff and Catalano 2016) in the command line version. Searches were conducted with new technology search level 5, finding the best score 100 times. Gaps were treated as fifth state. To infer branch support, we calculated bootstrap (BT) as absolute frequencies (Farris et al. 1996) using new technology search under level 2 and hitting the best score 5 times, totaling 1,000 replicates. These parameters (level of search and best score hit number) were chosen based on achieving the best length in preliminary analyses. BI and MP trees were edited with FigTree (Rambaut 2014). We considered MP bootstrap (BT) values ≥ 70% and BI posterior probabilities (PP) above 0.95 to indicate strongly supported clades (Hillis and Bull 1993; Alfaro et al. 2003; Erixon et al. 2003).

We calculated uncorrected p distance (proportion of nucleotide sites at which two compared sequences are different) for both gene fragments, within and between clades containing L. fuscus specimens using MEGA 7 (Kumar et al. 2016), with the “pair-wise deletion” option of gaps/missing data treatment.

To define independent evolving lineages, we used the multispecies coalescent model (Rannala and Yang 2003; Edwards 2009; Yang and Rannala 2014) implemented in BPP v.4.6.2 (Yang and Rannala 2010). We applied the analysis A11 that estimates the number of populations and topology of trees simultaneously. For this analysis we used a matrix that contained mitochondrial sequences of all L. fuscus specimens, and a guide tree based on the lineages found in phylogenetic analyses. We estimate priors for population size (θ) and divergence time (τ) with Minimalist BPP (https://brannala.github.io/bpps). We ran two replicates with the two rjMCMC algorithms (0 and 1) 100,000 generations, with 8,000 generations discarded as burn-in, with automatic finetune.

Based on the BPP delimitation, we estimated genetic differentiation (F_ST) and gene flow (Nm, derived from F_ST) between lineages, analyzing each gene fragment separately using DnaSP 5.10.01 (Librado and Rozas 2009). Following Wright (1965) we considered values of F_ST between 0 and 0.05 as indicative of no genetic differentiation, between 0.05 and 0.15 as indicative of limited genetic differentiation, between 0.15 and 0.25 as indicative of moderate genetic differentiation, and larger than 0.25 as indicative of substantial genetic differentiation. Regarding Nm, values <1 indicate that the differentiation between populations is due to genetic drift (i.e., no gene flow), >1 indicates high gene flow and low genetic differentiation, and >4 indicates nearly complete gene flow and very low differentiation (Grant and Bowen 1998).

We also estimated gene flow using the Bayesian algorithm implemented in Migrate-n 5.0.4 (Beerli and Felsenstein 1999, 2001). Both mtDNA fragments were treated as a single dataset, and migration was estimated in both directions among all lineages (i.e., a full migration model). Priors for the mutation-scaled population size (θ) were uniform, ranging from 0 to 0.1, and for the mutation-scaled migration rate (M), from 0 to 3,000. Following a burn-in of 50,000 steps, 50,000 steps were recorded every 500 steps, resulting in a total of 25,000,000 genealogies sampled. Three heated chains were run with temperatures of 1.0, 1.5, and 3.0, and 1,000,000. The number of migrants per generation (Nm) was calculated using the formula: Nm_(ij) = θ_j*M_(ij).

The datasets used in phylogenetic analyses, lineage delimitation, genetic differentiation, and gene flow are available at https://github.com/franbrusquetti/chacoan_leptodactylus_taxonomy.

As in our molecular phylogenetic analyses, we placed special emphasis on the South American Gran Chaco and its surrounding areas, including topotypes of L. gualambensis (213 specimens). We also included specimens of L. fuscus from various localities, including the type locality of the species (55 specimens; see File S2). We compared specimens assignable to L. gualambensis with all the other Leptodactylus species. For species for which it was not possible to study specimens, we based comparisons on the pertinent literature, including original species descriptions and taxonomic revisions. Comparisons with L. fuscus sensu stricto were restricted to specimens from Suriname (type locality) and surrounding areas, in concordance with our phylogenetic results. Specimens initially identified as L. fuscus, which, based on the phylogenetic results, were not considered as L. fuscus sensu stricto, were also compared, but as L. aff. fuscus (407 specimens).

We based the comparisons between species in structures frequently used in taxonomic studies in terms of presence/absence, shape, position, size relative to other structures, and color pattern, following the terminology of Duellman (1970), Cei (1980), Heyer (1978), Heyer et al. (1990), and de Sá et al. (2014). We took the following measurements based on Duellman (1970): Snout−vent length (SVL), from the tip of the snout to the vent; head length (HL), from the tip of the snout to the posterior edge of the jaw articulation; head width (HW), measured at the level of the jaw articulation; tympanum diameter (TD), measured horizontally; eye diameter (ED), measured horizontally; interorbital distance (IOD), distance between the interior margins of the eye bulges; eye−nostril distance (EN), from distal corner of the orbit to proximal margin of external nostril; internostril distance (IND), between the interior margins of the nostrils; eye−snout distance (ESD), from the tip of the snout to the distal corner of the eye; forearm length (FAL), from elbow to the proximal margin of the thenar tubercle; thigh length (THL), from vent to the knee; shank length (SHL), from knee to heel; and tarsus length (TAL), from the tibiotarsal articulation to the proximal margin of the internal tarsal tubercle. We measured specimens with digital calipers, to the nearest 0.1 mm under a dissecting microscope. To avoid ontogenetic morphological variation, we measured only adult specimens. Adult males were determined by visual examination of secondary sexual characters (i.e., presence of vocal sac and vocal slits). Adult females were determined by the presence of ovarian follicles (when possible, i.e., abdominal cavity opened) coupled with the absence of vocal sac and vocal slits. Alternatively, we considered specimens with over 40 mm SVL as adults, since all confirmed adults were larger than this. Final datasets for morphometric analyses consisted in 11 males and 10 females of L. fuscus sensu stricto, 49 males and 33 females of L. gualambensis, and 33 males and 33 females of L. aff. fuscus.

To assess differentiation between the clades, we performed a Principal Component Analysis (PCA) and Canonical Discriminant Analysis (CDA). PCA was applied as an unsupervised method to explore patterns of multivariate variation, while CDA aimed to maximize discrimination between a priori defined groups based on mitochondrial clades. CDA plots included individual scores, group centroids, and spider diagrams connecting individuals to their respective group centroids. To reduce dimensional redundancy, pairwise Spearman correlation coefficients were calculated, and variables showing high collinearity (ρ > 0.75) were excluded. Additionally, variables with very low contributions to the first canonical axes in preliminary CDA analyses were also removed, as they did not contribute meaningfully to group discrimination. We also determine differences of each measurement between the three clades. We first tested each variable for normality with the Shapiro-Wilk test (Shapiro and Wilk 1965). As not all variables had a normal distribution, we used the non-parametric Kruskal-Wallis test (Kruskal and Wallis 1952). For each variable with significant difference, we applied the Dunn test (Dunn 1961) with Bonferroni correction, as a post-hoc test to detect between-clade differences. All morphometric analyses were conducted with the MorphoTools2 package (Šlenker et al. 2022) in R (R Core Team 2023) separately for males and females.

We analyzed advertisement call records from specimens identified as L. fuscus from its type locality (Suriname), from the type locality of L. gualambensis (Urundel, Argentina) and from other localities through the known distribution of the species (File S3). The records were digitized with 16-bit resolution and 44,100 Hz sampling rate using the software Raven Pro 1.6 (Cornell Ornithology Lab). We produced audio-spectrograms and oscillograms, and measured the selected variables with Fast Fourier Transformation, using a window of 256, and contrast and brightness of 50%. We estimated call duration (in seconds), dominant frequency (as the peak frequency function of the software, in Hz), 5% frequency (in Hz), 95% frequency (in Hz), 90% bandwidth (difference between 95% and 5% frequency, in Hz), low frequency (in Hz), high frequency (in Hz), and the difference between high and low frequencies (delta frequency, in Hz). The terminology and the “call centered” approach were applied following Köhler et al. (2017). Figures were generated with the R package seewave 2.2.3 (Sueur et al. 2008).

To assess call differentiation between the clades, we first tested each variable for normality with the Shapiro-Wilk test. As no variable had a normal distribution, we used the non-parametric Kruskal-Wallis test to determine differences of each bioacoustics variable between the three clades. We applied the Dunn test as a post-hoc test for the variables with significant difference to detect between-clade differences. The mentioned tests were conducted with the MorphoTools2 package in R.

The best fit model for 16S and for COI codon 1 was GTR (General time reversible) + I (proportion of invariable sites) + G (gamma distribution), for COI codon 2 was F81 + I, and for COI codon 3 was GTR + G. Overall, both BI (Figs 2, S1) and MP (Fig. S2) analyses resulted in very similar topologies. The “Leptodactylus fuscus species complex” is recovered as monophyletic and well supported in both analyses (PP 1.0, BT 96%) and structured in three main clades, of which two were well supported and one poorly supported. The well supported clades correspond to two main geographic areas: Central America + northern South America (L. fuscus clade 1; PP 0.97, BT 79%) and southwestern South America (southeastern Bolivia, southwestern Brazil, western Paraguay, northern and northwestern Argentina; L. fuscus clade 2; PP 1.0, BT 84%). The unsupported clade contains the rest of the samples, which are from eastern Paraguay, northeastern Argentina, and northern, central and southeastern Brazil (L. fuscus clade 3; PP 0.86, BT < 50%). The relationships among the three clades were not well supported.

Regarding the L. fuscus clade 1, two samples, geographically distant from all the other samples of the clade, were clustered together, nested within in MP analysis, and as sisters of the rest of the clade in BI. One of the samples is from Beni department, northern Bolivia (MJ 1357) and the other is from Pernambuco state, northeastern Brazil (USNM 284551). Besides the above-mentioned samples, the L. fuscus clade 1 grouped all samples from Panama, Trinidad, Tobago, Venezuela, Guyana, French Guiana, Colombia, Roraima State in Brazil, and Suriname, type locality of L. fuscus. The L. fuscus clade 2 grouped all samples from Dry Chaco, western Humid Chaco, northern Pantanal, southern Cerrado, and Chiquitano dry forest, including samples from the L. gualambensis type locality (Urundel, Salta Province, Argentina). The L. fuscus clade 3 includes all samples from southern State of Bahia, Brazil (southern Caatinga and Bahia coastal forest), Cerrado and transition with southern Caatinga and Upper Parana Atlantic Forest.

The uncorrected pairwise genetic distances (p distances) between clades in both gene fragments were similar, with means around 3% in 16S and 10% in COI. Within the clades, we found larger internal genetic distances in the L. fuscus clade 3, with values up to 3.3% in 16S and 7.5% in COI (Table 1).

Lineage delimitation using BPP recovered the same three lineages identified in the phylogenetic analyses. All runs, two replicates using algorithms 0 and 1, recovered the L. fuscus lineages with posterior probabilities of 1.0. The same tree topology as in the phylogenetic analyses was also consistently recovered (Table 2).

Genetic differentiation (F_ST) among all populations showed similar values for both gene fragments, ranging from 0.55 to 0.68, indicating substantial genetic differentiation. Nm values based on F_ST were <1 for all lineage comparisons, suggesting that differentiation is mainly driven by genetic drift (i.e., no gene flow). Migrate-n results support low immigration rates, with mean Nm values among populations ranging from 0.47 to 0.59 individuals per generation (Table 3).

We identified some character states that distinguished specimens attributable to the different clades, corresponding to shape-related characters, such as snout shape in dorsal view, which is rounded to nearly rounded in individuals of the L. fuscus clade 2, but sub-elliptical to subovoid in the L. fuscus clade 1 (Fig. 3).

Regarding morphometric analyses, the CDA revealed clear morphological separation among genetic groups in both sexes. In males, the first canonical axis (Can1) explained 69.2% of the total variation and was primarily associated with HL (0.81), SVL (–0.68), and IND (0.44), while the second axis (Can2, 30.8%) was most influenced by SHL (–0.94), IOD (0.56), and ESD (–0.72) (Fig. 4A, B; Table 4). In females, Can1 explained 72.8% of the variation and was driven by HL (–0.75) and IOD (0.37), while Can2 (27.2%) was influenced by SHL (–0.78), ED (–0.58), and TD (–0.50) (Fig. 4C, D; Table 4). The PCA supported these general patterns, although group overlap was greater. In males, PC1 and PC2 together explained 54.7% of the variance, with HL (–0.39), ESD (–0.39), SHL (–0.37), and EN (–0.32) loading most heavily on PC1 (Fig. S3A, B; Table S1A). In females, PC1 and PC2 explained 51.0% of the variance, with HL (–0.44), EN (–0.46), and SHL (–0.37) as the dominant contributors to PC1 (Fig. S3C, D; Table S3B). In Kruskal-Wallis and Dunn tests we found statistically significant differences in some morphometric proportions among the members of the different clades. For example, members of the L. fuscus clade 2 tend to have a proportionally shorter and wider head than members of the L. fuscus clade 1 (p<0.01), and proportionally shorter shanks than members of the L. fuscus clade 3 (p<0.01; Fig. 5; Tables 5, 6).

The advertisement call corresponding to the three clades is structurally the same. It consists of a short whistle-like note with an increasing frequency throughout most of the call, except at the end, where it shows a slight drop. The call also has harmonic structure with up to three additional harmonic bands, being the first one the most energetic (Fig. 6). In the Kruskal-Wallis test we found statistically significant differences between the three clades in all the bioacoustics variables (Tables 7, 8). Based on the post-hoc test (Dunn-test) the L. fuscus clade 2 had advertisement calls with shorter duration than the others, as well as higher delta frequency. Leptodactylus fuscus clade 3 had lower frequency and longer duration than the others (Table 8).

From the morphometric and bioacoustics variables that significantly differentiated the clades, we proposed as diagnostic characters only those with null or low overlap values. These include head length, leg length, advertisement call duration, and advertisement call frequency (Figs 5, 7).

Based on the nesting position of the topotypes of L. fuscus (NZCS-A 2237–8 from Surinam) within the L. fuscus clade 1, we consider all samples of that clade as L. fuscus sensu stricto, restricting the distribution of the species to the northern South America and Central America, from northern Bolivia to Panama (Fig. 1). Similarly, all topotypes of L. gualambensis (LGE 23547–9 from Urundel, Salta, Argentina) were nested within the L. fuscus clade 2. Based on these results, together with morphological and bioacoustic evidence, we propose the revalidation of Leptodactylus gualambensis Gallardo, 1964, restricting its distribution to the Dry Chaco, its transition with the Yungas, the Humid Chaco, the Pantanal, and the southern Cerrado (Fig. 1). For the L. fuscus clade 3, which includes all the remaining samples, from eastern Paraguay and northeastern Argentina to northern Brazil (south of the Amazonas River), we maintain it as Leptodactylus aff. fuscus due to its weakly supported relationships and high internal genetic distances. We expect that a deeper sampling could lead to the detection of more than one independent lineage in this clade, which may correspond to some of the other names available under L. fuscus synonyms (see Discussion).

Gallardo (1964) compared L. gualambensis with L. fuscus (as L. sibilatrix) based on specimens from localities of the Argentinian provinces of Corrientes and Misiones, the Brazilian states of São Paulo and Rio Grande do Sul, and northern Uruguay, which, judging by the geographic distribution, could be attributed to L. aff. fuscus (L. fuscus clade 3). As diagnostic characters Gallardo mentioned that L. gualambensis has shorter legs (especially thighs and shanks); a more concave loreal region, with a major inclination to the outside; interocular spots situated further back; and with all tibial transversal stripes of equal width, and the middle one notably wider in what he considered L. fuscus. Although with some overlap, we also found that L. gualambensis differs from specimens of the clade L. aff. fuscus by leg length, but in this case mainly in the relative length of the shank. We found that the position of the interocular spots and the tibial color pattern were highly variable in our sample and, therefore, not useful as diagnostic characters.

As mentioned previously, the ratios related to shank length (SH/SVL) partially overlap among L. gualambensis and L. aff. fuscus specimens. Notably, most L. gualambensis specimens with overlapping values (SH/SVL >0.54, 15% of 49 male specimens; >0.51, 33% of 33 female specimens) are from the eastern distribution of the species, in the Humid Chaco and Cerrado, particularly from localities near transition areas between these ecoregions and the Atlantic Forest. In contrast, virtually all specimens (except for one female, IIBP-H 3634) from the western Humid Chaco and Dry Chaco have relatively short legs (SH/SVL <0.54, 85% of 49 male specimens; <0.51, 67% of 33 female specimens).

We suggest that the variation in leg length, which seems to reflect the geographical distribution of the specimens, could be associated with the exposure to different climatic conditions. Different environments may exert local selective pressure, resulting in phenotypic differentiation for species with a wide geographical range and populations facing diverse climatic conditions (Nali et al. 2023).

One possibility is that the variation in the hind limb length is related to burrowing behavior. In semiarid environments, a strategy used by terrestrial frogs to reduce the risk of desiccation during prolonged periods of drought is to bury themselves in self-made underground burrows (Carvalho et al. 2010). Heyer (1978) suggested that in species of the L. fuscus group the burrowing behavior is exclusive to males and used only for breeding purposes; burrows are used as incubation chambers. However, in the Caatinga, a semiarid ecoregion in northeastern Brazil, an individual of L. fuscus was found at a depth over 100 cm depth, aestivating in a skin cocoon (Varjão and Ribeiro 2018). A similar strategy may be employed by individuals of L. gualambensis to survive in drier parts of their distribution and this behavior could be reflected in the hind limb length, as in other species of head-first burrowers (e.g., Myobatrachidae, Vidal-García et al. 2014). Furthermore, although no data exists regarding digging behavior in the species, based on descriptions of burrow construction in other species of the L. fuscus group, an active function of hind limbs could be expected. Males of L. fuscus sensu stricto construct burrows by pressing their snouts against the ground, powered by their hind limbs (Martins 1988). Similar burrowing behavior was reported for L. bufonius, although a specific function of the hind limbs was not mentioned (Philibosian et al. 1974). We suggest that to reduce the risks of desiccation, individuals from populations inhabiting areas with more pronounced seasonal climate, such as the Dry Chaco, are forced to rely on fossorial behavior more frequently than individuals inhabiting areas with more homogeneous climatic conditions, such as those in the Humid Chaco and Cerrado. Based on the idea that hind limbs play an active role in digging behavior in L. gualambensis, and that shorter legs are more efficient in assisting digging in head-first burrowers, shorter leg lengths could be related to an increased fossoriality in populations inhabiting drier areas.

Another possible explanation for the intraspecific variation in L. gualambensis hind limb length is also related to the different climatic conditions faced by populations of the species, specifically with the developmental effects influenced by the duration of larval periods. Some species from arid or semi-arid environments exhibit a positive correlation between the larval period and hind limb length, potentially resulting in intraspecific morphological variation (Gómez-Mestre and Buchholz 2006). A suggested explanation of this phenomenon is that during pre-metamorphosis, a short larval period determines a short elongation phase of the hind limbs (Gómez-Mestre and Buchholz 2006). The larval period is influenced by diverse factors (see Székely et al. 2017), for pond-breeding species inhabiting semi-arid environments one of the main factors is the duration of temporary ponds. Some level of plasticity in larval development time in response to pond duration is an important adaptation to survive in an environment where breeding sites availability is a highly unpredictable resource (Newman 1992).

Leptodactylus gualambensis deposits eggs in a foam nest inside the underground nuptial chamber, constructed at the water´s edge (e.g., temporary ponds). During rains, the tadpoles are washed out of the nest and develop as free larvae in the ponds (Fabrezi 2011, as L. fuscus). Differences in the availability and duration of superficial water throughout its geographic distribution are marked. In the Humid Chaco and Cerrado, where higher rain regimes and clayish soils prevail, permanent or semi-permanent superficial waters are very common. In contrast, a large portion of the Dry Chaco, characterized by less rain and sandy soils, offers mainly temporary ponds that may last only a few days. This could accelerate the metamorphosis and reduce the duration of the elongation phase of the hind limbs. We suggest that in individuals from populations in the drier portions of the species’ distribution, where the reproduction occurs mainly in temporary ponds of short duration, comparatively shorter legs could result from a shorter elongation phase.

Other variation in morphological characters could also be related to inhabiting different environments in this species. Rhinella granulosa, with a similar geographic distribution pattern that includes contrasting environments, shows differences in ventral skin width and vascularization. Individuals from more xeric environments show thinner and more vascular ventral skin, a combination of traits more effective for water uptake (Navas et al. 2004). Leptodactylus gualambensis, with a geographic distribution including dry and highly seasonal environments as well as humid and relatively homogeneous environments, is a promising model for studying morphological adaptations associated with surviving in semiarid environments at intraspecific level and for better understanding the history of the Chaco anurofauna.

Although our main goal was to assess the taxonomic validity of L. gualambensis, the available names for different populations within L. fuscus deserve some discussion. As mentioned above, several synonyms have been attributed to L. fuscus, and throughout these years, with the goal of clarifying the taxonomic situation of the L. fuscus species complex, samples from type localities or localities attributable to some of these names have been included in different studies (e.g., Wynn and Heyer 2001; Heyer and Reid 2003; Camargo et al. 2006) but never discussed.

Before addressing the other names, we need to clarify an issue detected in the type series of L. gualambensis. Gallardo (1964) designated the specimen MACN 9752 as the holotype and included a drawing of a head in dorsal view in page 47, and a photograph of the complete body of an unlabeled specimen in the Plate II. In both figures the legend refers to the specimen represented as the holotype, with the following sentence: Leptodactylus gualambensisnov. sp., MACN 9752 (typus). By the position of the legs is possible to note that the photographed specimen is the one currently labeled MACN 9753 (Fig. 9). To determine whether this was an error in the photograph or caption of Plate II in Gallardo (1964), or an error made when labeling the holotype, we examined the species description. However, it appears that only a few parts of the description are based solely on the holotype (e.g., measurements), while others, such as the reference to the medio-dorsal stripe being “frequently present”, suggest an observation based on a group of specimens. Based on new measurements of the specimens MACN 9752 (SVL 43.76 mm) and MACN 9753 (SVL 45.53 mm), along with the total body length reported by Gallardo (1964) as 47 mm, it is evident that the specimen designated as holotype in the original description of L. gualambensis corresponds to the one currently labeled as MACN 9753, reinforcing the idea of a mistake when labeling the specimens. Given that both the photograph and the measurements indicate that the current MACN 9753 is the actual holotype described, despite the original description attributing the number MACN 9752 to it, we recommend switching the current labels to assign MACN 9752 to the correct specimen.

Regarding the oldest names, such as R. fusca and Rana typhonia Daudin, 1802, we do not have new information, and the synonymy with L. fuscus is currently well accepted (see Heyer 1968, 1978; de Sá et al. 2014). We relate those names to the L. fuscus clade 1 because both share Suriname as their type locality (Heyer 1978). On the other hand, the synonymy of Rana virginica Laurenti, 1768 with L. fuscus is controversial. Both Heyer (1968) and Dubois and Ohler (2009) considered that the figure in Seba (1734: pl. 75, fig. 4), type of the species, represents a specimen of the genus Lithobates Fitzinger, 1843 (as Rana pipiens Schreber, 1782 in Heyer 1968) because of its webbed feet. However, in de Sá et al. (2014), this information was ignored, and R. virginica was included in the list of synonyms of L. fuscus without any comment. Furthermore, besides the doubts about the correspondence between Seba’s (1734) figure and L. fuscus, the type locality of R. virginica was not designated (de Sá et al. 2014). Considering the mentioned evidence, we suggest the removal of the name Rana virginica Laurenti, 1768 from the synonymy of Leptodactylus fuscus and designate it as a nomen dubium until new evidence clarifies its application.

Another name that could be considered in the L. fuscus clade 1 is Cystignathus schomburgkii Troschel, 1848 (type locality Guyana). Boulenger (1882) synonymized this species with Leptodactylus typhonius (Daudin, 1802). Heyer (1978), following some traits from the original description, such as the species being closest to C. gracilis and having a uniform brown dorsal color, suggested that L. fuscus, L. longirostris (based on the C. gracilis similarity), and northern populations of L. mystaceus (based on dorsal color pattern), are the only known species from Guyana that could correspond to Troschel’s description.

In the absence of type material (considered as lost), Heyer (1978) tentatively treated the name C. schomburgkii Troschel, 1848 as a L. fuscus synonym. Recently, the syntypes were found in the Museum für Naturkunde, Berlin, Germany. They were both originally labelled ZMB 3337, but one of them was renumbered as ZMB 75292 in November 2010 (Frank Tillack pers. comm.). The locality, according to the original catalog entry and jar label, is “Guiana”, both collected by Robert Herman Schomburgk. We had access to pictures of both specimens (Fig. 10). One of them shows serious skin and muscle damage on the left side of the dorsal surface (ZMB 3337).

Both specimens lack toe fringes, restricting their taxonomic identity to the species of the L. fuscus and L. pentadactylus groups (all species of the L. latrans group and L. melanonotus group have toe fringes; de Sá et al. 2014). Furthermore, both are adult males (evident vocal sac; Fig. 10A, B) and lack chest and thumb spines, which restricts them to the L. fuscus group (all species in the L. pentadactylus group have chest and/or thumb spines; de Sá et al. 2014). According to Cole et al. (2013), distribution maps provided by de Sá et al. (2014), and previous reports by Heyer (1978), only three species of the L. fuscus group are present in Guyana; L. fuscus, L. longirostris, and L. mystaceus.

Both syntypes of C. schomburgkii, in general, appear faded; however, the dorsal dermal folds are still evident (Fig. 10C, D). These folds are present only in L. fuscus and in individuals of L. longirostris with median dorsal stripe (Heyer 1978; de Sá et al. 2014; Silva et al. 2020). Another distinguishing character between L. fuscus and L. longirostris compared to L. mystaceus is the shape of the fingertips: L. fuscus and L. longirostris have narrow and thin fingertips, whereas L. mystaceus has globular fingertips (Heyer 1978; Toledo et al. 2005; de Sá et al. 2014; Silva et al. 2020; PDPP pers. observation).

Furthermore, in the original description of C. schomburgkii, Troschel (1848) mentioned that the vocal sac is very notable when inflated, forming a bubble-like extension on each side of the head. In L. fuscus and L. longirostris the vocal sac is double and lateral, while in L. mystaceus it is subgular and not expanding laterally (Heyer 1978; Toledo et al. 2005; de Sá et al. 2014). Both syntypes exhibit a double, lateral vocal sac (Fig. 10A, B), dorsal folds, and narrow, thin fingertips (Fig. 10C, D). Those traits align with L. fuscus and L. longirostris, distinguishing them from L. mystaceus.

By discarding L. mystaceus from the list of potential synonyms with C. schomburgkii, two possibilities remain: L. fuscus and L. longirostris. According to Heyer (1978) and de Sá et al. (2014) there are no clear morphological characteristics that differentiate those species, as the intraspecific variation of both overlaps. Our examination of specimens of both species collected in Guyana also revealed no morphological differences (see Material Examined). This presents a challenge in resolving the taxonomic identity of C. schomburgkii, rendering its synonymization with L. fuscus or L. longirostris inconclusive.

In light of this uncertainty, we follow the ICZN (1999) recommendation on the stability of nomenclature, maintaining the name Cystignathus schomburgkii Troschel, 1848, as a junior synonym of Leptodactylus fuscus (Schneider, 1799), as designated by Heyer (1978). However, it would be important to study the population genetics of L. fuscus and L. longirostris in Guyana. Given their sympatry and morphological similarities, the possibility of a complex interaction between them, including potential hybridization, cannot be ruled out.

Besides L. gualambensis, another L. fuscus synonym with a well-defined type locality is Leptodactylus raniformis Werner, 1899. In the phylogenetic analyses, we included a sample from the region of the type locality (previously included by Guarnizo et al. 2015; from Meta, Colombia, Andes-A 1763). Heyer (1978) mentioned that the holotype of L. raniformis (ZFMK 28484, adult male, Fig. S4), although with indistinct dorsolateral folds, presents a combination of traits characteristic of L. fuscus, such as the dorsal spots pattern, and tarsal and foot surfaces smooth with light pigmented spots. Based on the conclusion of Heyer (1978) and on our resulting topology (Fig. 2), with all samples from Colombia grouped together with the samples from Suriname (type locality of L. fuscus) in the L. fuscus clade 1 (i.e., L. fuscus sensu stricto), we corroborate that L. raniformis is a junior synonym of L. fuscus.

Two names are available for populations of L. fuscus from Brazil: Rana pachypus var. 2 Spix, 1824 from Pará State and Rana sibilatrix Wied, 1824 from the coast of southern Bahia state. The only character of R. pachypus var. 2 mentioned by Spix (1824) that could be attributable to L. fuscus is the presence of six longitudinal folds. However, the synonymy was suggested by Peters (1872), who according to Hoogmoed and Gruber (1983) was the only one who examined the type-specimen, currently lost. We were not able to include any sample that could be assigned to R. pachypus var. 2 in our analyses. Vanzolini (1981), based on Spix and Martius’ itinerary, restricted the possible type locality to the surroundings of the city of Belém (State of Pará). Our nearest sampled locality is Barreirinhas, Maranhão, which is 650 km E of Belém (Fig. 1). In view of the type specimen being lost and the lack of samples attributable to R. pachypus var. 2, we suggest maintaining the synonymy as proposed by Peters (1872).

Regarding R. sibilatrix, the short morphological description, comments about its vocalization, and the figure (see Wied 1824), correspond with a specimen that could be part of the L. fuscus species complex. Additionally, Heyer (1978) pointed out that L. fuscus is the only species with the characteristics described on the coast of southern Bahia. Heyer (1978) also mentioned some uncertainties about the correspondence between the only available syntype (AMNH A-485) and the figure published in Wied (1824). Vanzolini and Myers (2015), based on a direct comparison of the syntype and Wied’s figure, agree with Heyer (1978) that it is not the specimen drawn. This idea is reinforced by size differences between the syntype AMNH A-485 (SVL of about 36 mm, estimated from the photograph in Vanzolini and Myers 2015) and the size mentioned by Wied (~ 42.8 mm, “an inch plus seven lines”), although without a reference to a specific specimen in Wied (1824). Vanzolini and Myers (2015) recommended that, if necessary, either the only existing syntype (AMNH A-485) or the lost specimen drawn in Wied (1824) could be designated as lectotype.

We refrain from suggesting any taxonomic decisions regarding R. sibilatrix given the lack of support for the L. fuscus clade 3, and our sampling missing any locality at least near its type locality. Müller (1927) restricted the type locality of R. sibilatrix to the first locality mentioned by Wied (1824), Peruipe River, Nova Viçosa (Peruhype bei Villa Viçoz). The other localities mentioned by Wied are all proximal to Nova Viçosa; Caravelas to the north and Mucurí to the south, both also on the coast of southern Bahia. Our sequences from the closest localities to these mentioned are from Linhares, Espírito Santo (CFBH 32037), 195 km S of the type locality, and from Itacaré (CFBH 21108) and Uruçuca (CFBH 32435), Bahia, 400 and 350 km N from Nova Viçosa, respectively.

The geographic distances between our samples and the type locality are not the only challenges in associating any of the samples with this name. The samples mentioned were grouped into different sub-clades, which, although not well-supported, show relatively high genetic distance among them (16S p distance from 1.4 to 3.3%). Within sub-clades, genetic distances among some localities, such as Buritizeiro, Minas Gerais state, Brazil (LHUFCG 491) and those from Bahia state, Brazil (CFBH 21078, CFBH 21086, and CFBH 21108), were also relatively high (16S p distance from 1.6 to 2.3%).

As mentioned above, our main goal was to test the taxonomic validity of L. gualambensis, focusing our efforts on sampling the Chaco and surroundings areas. With our current sampling, we are not able to clarify the taxonomic status of all populations of L. fuscus and the still available names; however, we believe this study significantly reduces the problem. With the recognition of L. gualambensis as a valid species for Chaco populations and the restriction of L. fuscus sensu stricto (and its synonyms) to the populations of northern South America and Central America, only two names remain unresolved: Rana sibilatrix and Rana pachypus var. 2.