Integrative taxonomic reassessment of Odontophrynus populations in Argentina and phylogenetic relationships within Odontophrynidae (Anura)

Amphibians are the most vulnerable vertebrates to biodiversity loss mediated by habitat destruction, climate change and diseases. Informed conservation management requires improving the taxonomy of anurans to assess reliably the species’ geographic range. The genus Odontophrynus that is geographically refined to Argentina, Bolivia, Brazil, Uruguay and Paraguay includes currently 12 nominal species with many populations of uncertain taxonomic assignment and subsequently unclear geographic ranges. In this study, we applied integrative taxonomic methods combining molecular (mitochondrial 16S gene), allozyme, morphological and bioacoustic data to delimit species of the genus Odontophrynus sampled from throughout Argentina where most species occur. The combined evidence demonstrates one case of cryptic diversity and another of overestimation of species richness. The populations referred to as O. americanus comprise at least three species. In contrast, O. achalensis and O. barrioi represent junior synonyms of the phenotypically plastic species O. occidentalis. We conclude that each of the four species occurring in Argentina inhabits medium to large areas. The Red List classification is currently “Least Concern”. We also propose a phylogenetic hypothesis for the genus and associated genera Macrogenioglottus and Proceratophrys (Odontophrynidae).


INTRODUCTION
Distribution pattern of Neotropical amphibian diversity is not well understood because of incomplete information on taxonomy and distribution (Vieites et al., 2009;Winter et al., 2016). Yet amphibians are of high conservation concern, with almost 43% the currently known species being globally threatened and another 25% data deficient (Stuart et al., 2004). Taxonomic uncertainty stems partially from the prevalence of the morphospecies concept in most original descriptions of amphibian species (Frost, 2018).
It remains controversial, if diploids of the O. americanus group derived from tetraploids or tetraploids several times independently from diploids (Beçak & Beçak, 1974;Beçak, 2014). With respect to these issues and the validity of the phenetic groups within Odontophrynus, the most recent molecular phylogeny of Odontophrynus is inconclusive being based on only three of the 12 nominal taxa (Dias et al., 2013). Earlier phylogenies proposed by Amaro, Pavan & Rodrigues (2009), Pyron & Wiens (2011) as well as the recent one by Feng et al. (2017) agree in that Odontophrynidae is monophyletic and that Macrogenioglottus and Odontophrynus are sister taxa.
Consequently, a reliable taxonomic and geographic delimitation of Odontophrynus species requires an integrative approach critically evaluating information available on genetic differentiation and phenotypic plasticity. In this study we adopt the consensus protocol for integrative taxonomy to delimit species . Geographically, we focus on Argentina where most currently recognized Odontophrynus species and several populations of still undetermined taxonomic status occur. The character complexes included in the re-assessment of taxa are quantitative morphometrics, advertisement call features, allozyme differentiation and partial mitochondrial 16S rRNA sequences, all providing meaningful taxonomic information. Data refer to 34 populations, among them those at the type localities for reference. Sites of sympatry (O. americanus/ O. occidentalis, O. cordobae/O. occidentalis) are contrasted with those in narrow contact zones (O. americanus/O. cordobae, O. achalensis/O. occidentalis) and sites of allopatry. Additional data on the molecular Odontophrynus diversity in Brazil are used for a broader phylogenetic view on Odontophrynidae (Amaro, Pavan & Rodrigues, 2009;Dias et al., 2013;Lyra, Haddad & De Azeredo-Espin, 2017). Specifically, we test the following hypotheses: (1) Phenotypic plasticity within and among Odontophrynus taxa is associated with corresponding genetic differentiation; (2) The phenetic groups within Odontophrynus represent distinct phylogenetic lineages; (3) The current Red List classification does not reflect the risks upon each species.

Study area and field sampling
We identified and sampled 34 local populations of toads pertaining to the genus Odontophrynus in Argentina (Table S1). The type localities of the nominal taxa O. achalensis Di Tada et al., 1984 (Pampa de Achala, Cordoba province), O. barrioi  (Aguadita springs, Sierra de Famatina, La Rioja province), O. cordobae Martino & Sinsch, 2002 (Villa General Belgrano, Cordoba province) and O. lavillai Cei, 1985 (Villa de la Punta, Santiago del Estero province) were sampled to obtain topotypical individuals for taxonomic comparison. Unfortunately, the exact type localities of the most wide-spread species O. americanus (Duméril & Bibron, 1841) and O. occidentalis (Berg, 1896) are unknown because the original descriptions only state that the holotype of O. americanus was "sent from Buenos Aires" and that the holotype of O. occidentalis was collected in "Arroyo Agrio" in the Neuquén province (Frost, 2018). According to Lescure et al. (2002), the type locality of O. americanus is probably "rives du Rio Negro, Patagonie" in the Río Negro province (Argentina; D'Orbigny, 1847) at the southernmost part of the current geographical range. Still, populations of tetraploid O. americanus were readily distinguished from those of the diploid taxa by erythrocyte size (Rosset et al., 2006;Grenat, Salas & Martino, 2009). Populations of uncertain taxonomic assignment were tentatively referred to as O. cf. achalensis (Locality: La Carolina, San Luis province) or O. cf. barrioi (Localities: Aguada de Molle, Huerta de Guachi, San Juan province; Table S1). Material and data collected at the study sites were: (1) blood smears for ploidy assessment; (2) adult specimens for morphometric measurements, (3) records of advertisement calls, (4) muscle and liver homogenates for allozyme analyses and (5) alcohol preserved tissue for phylogenetic analyses (partial sequences of the mitochondrial 16S rRNA gene). The carcasses of specimens studied were deposited in museum collections; Table S1. The Córdoba Environment Agency (A.C.A.S.E.), Environmental Secretary of Córdoba Government (A01-2013), authorized our study and issued research and collecting permits.
Each variable was standardized by subtracting the sample mean and dividing by the sample standard deviation. The matrix of standardized variables was subjected to a principal component analysis with a fixed number of three PCs extracted. By this means, we explored the morphometric variability independent of taxonomic assignment and reduced the information to statistically unrelated factors. PC1 represents size-related features, PC2 and PC3 shape-related ones. Separate PCAs were run on the taxa of the phenetic groups. Assignment of populations to a phenetic group was based on the advertisement call structure (O. americanus-group: simple pulsed calls; O. occidentalis-group: complex calls consisting of several pulse groups; Salas & Di Tada, 1994;Martino & Sinsch, 2002). The morphospace built by three PC-axes was used to evaluate partitioning among taxa. A discriminant analysis (Procedure: backward selection) with a priori taxon assignment was applied to quantify the partitioning of morphospace for each phenetic group. Removal of one or more morphometric variables is based on an F-to-remove test. If the least significant variable has an F-value less than 4.0, it will be removed from the model. We used the discriminant analyses to assess the magnitude of correct taxon classification for the individuals of each predefined group. Due to the partial sympatry among species of the O. occidentalis group we tested for clinal variation of PCs along latitudinal and altitudinal gradients by a multiple regression analysis (Procedure: backward selection at F = 4.0). Significance level was set to a = 0.05. All calculations were performed using the statistical package statgraphics centurion, version XVIII (Statpoint Inc., Warrenton, Virginia, USA).

Bioacoustic data
Advertisement calls given by 302 individuals (series of 11-116 calls per individual) were recorded in the field using a DAT recorder Sony TCD-100 © with stereo microphone ECM-MS907 Sony © and tapes TDK DA-RGX 60 © (Table 1). Ambient temperature (to the nearest 0.5 C) was registered at the individual calling sites (usually shallow water near shore) immediately after recording. Short advertisement call series are available as Audio S1-S8. Whenever possible, specimens were collected to obtain tissue samples and for morphometric measurements. Oscillograms, sonograms and power spectra of the call series were prepared with the Medav Mosip 3000 Signal Processing System or the PC program Adobe Audition 1.0. Each call series was characterized by nine parameters which were measured in three calls per series (terminology and procedure according to Martino & Sinsch (2002) and Schneider & Sinsch (2007)): (1) call duration (ms); (2) number of pulse groups per call (N); (3) duration of pulse group (ms); (4) interval between pulse groups; (5) pulses per pulse group (N); (6) pulse duration (ms); (7) interpulse interval (ms); (8) pulse rate (pulses/s); (9) dominant frequency (Hz). Bioacoustic raw data are available in Dataset S2.
The arithmetic means of these call parameters were calculated for each series (=individual) and used for further analyses. Thus, the basic data set describing the advertisement calls of the populations studied consisted of 10 variables (nine call parameters and the corresponding ambient temperature) with N = 302 observations. As several call variables co-vary with ambient temperature, we calculated linear regression models of call parameter vs. temperature and used the standardized residuals to obtain a temperature-adjusted data set for further analysis. Analogous to the treatment of morphometric data, a principal component analysis was run on call data subsets of populations with homologous call structure (simple calls with six variables vs. complex calls with nine variables) to explore the bioacoustic differentiation among the taxa of each phenetic group. The three PCs explaining most of the variance were extracted to describe the sound space utilized by Odontophrynus and its partitioning among taxa. Moreover, a discriminant analysis (procedure: backward selection at F = 4.0) was applied on standardized call variables to quantify the partitioning of among-taxon sound space. We used the discriminant analyses to assess the magnitude of correct taxon classification for

Molecular phylogenetic analysis
We compared the partial sequence of the mitochondrial 16S rRNA gene of the samples from the different localities in Argentina to assess the number of species present in the country and their phylogenetic relationships (Table S3). The 16S barcoding gene is widely used, it amplifies with great success for most amphibians, there is a lot of available information already, and it is generally informative for tree-tips relationships, that is, among closely related species (Vences et al., 2005;Zhang et al., 2013). DNA was extracted using Qiagen DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany) following the manufacturer's protocol. Polymerase chain reaction (PCR) was used to amplify fragments of approximately 560 bp of 16S mitochondrial rRNA using the standard primers 16SAL (5′-CGCCTGTTTACTAAAAACAT-3′), and 16SBH (5′-CCGGTCTGAACTCAGAT CACGT-3′). Amplification followed the standard PCR conditions (Palumbi, 1996) with the following thermal cycle profile: 120 s at 94 C, followed by 33 cycles of 94 C for 30 s, 49 C (12S)/53 C (16S) for 30 s, and extension at 65 C for 60 s. All amplified PCR products were verified using electrophoresis on a 1.4% agarose gel stained with ethidium bromide. PCR products were purified using Highpure PCR Product Purification Kit (Roche Diagnostics, Risch-Rotkreuz, Switzerland). Sequencing of both strands was performed with the DYEnamic ET Terminator Cycle Sequencing Premixkit (GE Healthcare, Munich, Germany) for sequencing reactions run on a MegaBACE 1000 automated sequencer (GE Healthcare, Munich, Germany). Chromas lite 2.1.1 software (Technelysium Pty Ltd, South Brisbane, Australia; http://www.technelysium.com.au) was used to check and read the chromatograms of the sequences. The obtained sequences were compared with those in GenBank using a standard nucleotide-nucleotide BLAST search. Homologous sequences of Odontophynus as well as from species of the closely related genera Macrogenioglottus and Proceratophrys were downloaded from GenBank and incorporated in an alignment. Sequences of Pleurodema somuncurensis (Leptodactylidae), Rhinella marina (Bufonidae) and Ceratophrys cornuta (Ceratophrydae) were used as outgroups (Table S3) to represent related anuran families (Feng et al., 2017). The sequences were aligned using the MUSCLE algorithm (Edgar, 2004) implemented in MEGA 7 (Kumar, Stecher & Tamura, 2016). The final alignment consisted of 552 bp. Pairwise distances were calculated in MEGA7. Distances >1% were considered indicative for differentiation at species level (Sáez et al., 2014).
The general time-reversible model with proportion of invariable sites and gamma-distributed rate variation among sites (GTR + I + G) was chosen as the best-fitting model of sequence evolution on the basis of the Akaike information criterion as implemented in jModelTest 2 (Darriba et al., 2012) and was applied in maximum likelihood (ML) and Bayesian inference (BI) analyses. ML was performed in MEGA 7 with heuristic searches with stepwise addition and TBR branch-swapping algorithm, generating 1,000 bootstrap replicates. BI was performed using MrBayes 3.2.5 (Ronquist et al., 2012). Two independent Metropolis-coupled Monte Carlo Markov Chain (Larget & Simon, 1999) analyses were run for 10 Million generations, each with one hot and three cold chains and the temperature set at 0.2. Trees were sampled every 5,000 generations. The first 500 samples of each run were discarded as burn-in, and the remaining trees from both runs were used to calculate a consensus tree and Bayesian posterior probabilities.

Morphological variation
All nominal taxa of Odontophrynus resemble each other considerably in coloration and external morphology reflecting their semi fossorial mode of living (Fig. 1). Quantitative morphometric analyses demonstrated a significant morphological variation among some taxa. The three principal component representing the axes of morphospace explained 77.2% of total variance in the O. americanus group and 80.1% in the O. occidentalis group, respectively (Table 1). The morphospace of the O. americanus group was partitioned between O. lavillai on one side and the indistinguishable pair O. americanus/O. cordobae on the other side (Fig. 2). The discriminant analysis confirmed a significant separation of O. lavillai with 85.7% of individuals correctly assigned to this  Table 2). Classification success of discriminant function varied between 70% and 80% (Table 2). A significant proportion of morphometric variability among individuals assigned to the O. occidentalis group was caused by a clinal variation along altitudinal and latitudinal gradients. Size-related variation (PC1) was significantly correlated with altitude and latitude (Multiple regression model, R 2 = 32.1%, F 2,102 = 24.03, P < 0.00001), that is, size of individuals increased with elevation and from south to north. PC2 (position of nares and eyes) was significantly correlated with latitude (Multiple regression model, R 2 = 16.6%, F 1,103 = 20.52, P < 0.00001), PC3 (HL) with altitude (Multiple regression model, R 2 = 10.0%, F 1,103 = 11.42, P = 0.001).

Advertisement call variation
The taxa of the O. americanus group emit simple and short pulsed advertisement calls, whereas those of the O. occidentalis group produce long and complex advertisement calls consisting of a variable number of short pulse groups (Fig. 4). Quantitative analyses of the advertisement calls based on six temperature-adjusted variables in the O. americanus group showed a significant variation among the three taxa. Three PCs explained 85.3% of total variance represented the axes of sound space (Table 1). The sound space was partitioned into three discrete groups representing O. americanus, O. cordobae and O. lavillai individuals, respectively (Fig. 2). The discriminant analysis assigned all but five calls correctly to the corresponding taxon (Table 3).  Analogous to morphometric variation, sound space partitioning was low among the taxa of the O. occidentalis group, with O. occidentalis, O. achalensis and O. cf. achalensis being indistinguishable among each other ( Fig. 3; Table 3). The sound space occupied by O. barrioi and O. cf. barrioi differed from the continuum formed by the other taxa, but showed a slight overlap between each other. Still, the correct classification rates for O. barrioi and O. cf. barrioi were 100% and 95.4%, respectively (Table 3). Temperature-adjusted advertisement call variation was also influenced by geographical clines. PC1 (size-related dominant frequency) and PC2 (call duration) were significantly correlated with latitude (Multiple regression models: R 2 = 9.8%, F 1,74 = 7.94, P = 0.0062, and R 2 = 14.1%, F 1,74 = 11.97, P = 0.0009, respectively), and PC3 (pulse group duration) with altitude (R 2 = 9.1%, F 1,74 = 7.31, P = 0.0085).

Genetic distances: allozymes and barcoding
A total of 14 putative loci were scored in the nominal taxa (Table S2). Two loci (mAAT, mMDH) were monomorphic in all taxa. The overall degree of allele fixation was high and varied between five loci in O. americanus and 11 in O. lavillai. Five taxa possessed   Argentina were at species level (1.  (Fig. 5).  Note: Analyses were run separately on the two phenetic Odontophrynus groups. The subset include the smallest combination of measured variables that maximizes classification success. For details see text.
Morphological diagnosis of tadpoles. In the exotrophic, lentic and benthic Odontophrynus spp. tadpoles, the body is ovoid in dorsal and ventral views and slightly depressed in lateral view, the snout is circular in dorsal and ventral views with an anteroventral oral disc, the eyes and nostrils placed dorsally, and the spiracle sinistral, short, and with an opening posterodorsally directed (Do Nascimento et al., 2013). Odontophrynus occidentalis differs from the other species by being considerably larger at stages 36-38 (average: 71.7 mm; maximum: 75 mm (Cei, 1987). Do Nascimento et al. (2013) and González et al. (2014) state that external tadpole morphology is very similar in the achalensis, barrioi and occidentalis ecotypes, unlike Cei & Crespo (1982). Shared characters apart from the general Odontophrynus spp. features are a tail, that is, larger than half of the total length with robust musculature and 2(2)/3(1) as labial tooth row formula. Among the distinguishing characters mentioned by , Cei & Crespo (1982) and Cei (1987) for the barrioi ecotype, only the diverging form of the tip of tail (rounded contrary to acuminate in occidentalis and achalensis) has been confirmed by González et al. (2014). The significance of this difference remains unclear because phenotypic plasticity of tail morphology in response to syntopic predators is well-documented (Van Buskirk, 2009).

Bioacoustic diagnosis. Odontophrynus occidentalis is unique among all currently known
Odontophrynus species by giving series of complex advertisement calls consisting of a variable number of pulse groups (Figs. 4B and 4C). Rosset et al. (2007) state that advertisement call features in the ecotypes achalensis, barrioi and occidentalis are very similar. The bioacoustic distinction of achalensis from occidentalis (Salas & Di Tada, 1994) and of barrioi from the other ecotypes (Rosset et al., 2007) is based mainly on differences in the dominant frequency (related to body size), in call duration (sensitive to environmental temperature), and in intercall interval (sensitive to behavioral interactions such as chorusing; Schneider & Sinsch, 2007). In fact, temperature-adjusted dominant frequency Standardized Dominant Frequency (SDF) was significantly correlated with SVL of the caller demonstrating that ecotype distinction reflects mainly size variation among populations (Linear regression model: SDF = 2.405-0.065 Ã SVL, R 2 = 0.623; F 1,74 = 120.6, P < 0.0001; Fig. S1). SVL-and temperature-adjusted call duration did not differ significantly among the ecotypes (ANOVA, F 4,74 = 2.14, P = 0.084).
Geographic and habitat range. Populations are present in a latitudinal range from 41 S in the Rio Negro Province to 27.5 S in the north of the Catamarca Province Rosset et al., 2007). The corresponding altitudinal range is from 690 m above sea level to 1,854 m in the Sierra Pie de Palo, 1,965 m in the Sierra de San Luis, 2,149 m in the Sierra de Achala, and 2,200 m in the Sierra de Famatina (Table S1). Terrestrial habitats include pastureland, montane grassland, shrubs, forests and rock outcrops. Aquatic habitats used for reproduction and larval development are exclusively creeks and slowly running rock pool sections of permanent streams.

DISCUSSION
Lines of evidence obtained from phenotypic and genotypic character complexes in Odontophrynus toads exemplify the common dilemma of taxonomy-which level of character differentiation requires taxonomic consequences? Our case study demonstrates that phenotypic plasticity may result in an overestimation of species number (O. occidentalis group), whereas molecular data may reveal unexpected cryptic diversity in morphologically uniform populations (referred to as O. americanus). The following discussion of the three hypotheses basic to our investigation will present a revised view on the actual Odontophrynus species richness and propose a model of the phylogenetic relationships within the genus Odontophrynus.
Hypothesis 1: Phenotypic plasticity within and among Odontophrynus taxa is associated with corresponding genetic differentiation Applying the consensus protocol for integrative taxonomy  we evaluate the support for the genetically delimited species O. americanus, O. cordobae, O. lavillai and O. occidentalis by the among-taxon variation of allozyme, morphology and bioacoustic character complexes. Congruence of genetic and phenotypic differentiation clearly delimits O. lavillai from the congeneric taxa in Argentina. Independent support for its distinction from the only sympatric Odontophrynus species stems from karyological features, that is, O. lavillai is diploid whereas topotypic O. americanus is tetraploid (Barrio & Pistol de Rubel, 1972;Cei, 1985;Rosset & Baldo, 2014). Natural hybrids between these species have not been detected.
The significant differentiation of 16S RNA gene sequence distinguishes between the tetraploid O. americanus and the diploid O. cordobae as does the advertisement call structure (Martino & Sinsch, 2002;this study). Congruent allozyme and morphometric differentiation seem to be absent (Barrio & Pistol de Rubel, 1972; this study). However, adjusting size by age reveals that the two species differ significantly in their ontogenetic growth trajectories (Martino & Sinsch, 2002) supporting the molecular distinction. Premating isolation mechanisms allow the two species to maintain reproductive isolation in the narrow zone of sympatry, as indicated by a very low incidence of triploid hybrids at only two localities (Grenat et al., 2018). Thus, cumulating available evidence on character differentiation (except for the allozyme pattern) supports the species status of these two taxa.
To our surprise, the 16S RNA gene sequence of a topotypical tetraploid O. americanus from the Buenos Aires province, Argentina and that of five specimens of unknown ploidy referred to as O. americanus from different localities in Brazil (Amaro, Pavan & Rodrigues, 2009;Lyra, Haddad & De Azeredo-Espin, 2017) differed so profoundly that we consider them as a complex of three distinct species. One sample from Rio Grande, Brazil proved to be very similar to the topotypical tetraploid O. americanus indicating that this species extends over Uruguay to southern Brazil. A second group of two specimens originating from southern Minas Gerais is termed here as O. aff. americanus 1 and a distinct species from topotypical tetraploid O. americanus. Its distinction from or similarity with the diploid O. juquinha from northern Minas Gerais remains to be studied in the future. Finally, a third group of three specimens from Sao Paulo and Santa Catarina provinces, termed here as O. aff. americanus 2, is a distinct species from Minas Gerais taxon and true O. americanus. One of these specimens originated from Irai which is about 100 km distant from the "O. americanus"-populations reported from the province Misiones, Argentina (Rosset et al., 2006). It seems reasonable to assume that the geographical range of O. aff. americanus 2 extends to Argentina in the west, but further studies are needed to resolve the relationship to the Misiones populations.
In contrast, we did not detect any significant genetic differentiation among the specimens belonging to different nominal taxa of the O. occidentalis group despite a considerable phenotypic plasticity among and within populations. The case of O. barrioi is more complicated because the absence of significant genetic differentiation from O. occidentalis contrasts with considerable morphometric, bioacoustic and allozyme differentiation. The O. barrioi populations vary morphometrically only by their greater size whereas shape variation is the same as in the O. occidentalis/ O. achalensis continuum. Within-species altitudinal and latitudinal size variation is well known in anurans (Sinsch, Pelster & Ludwig, 2015), and SVL differences alone appear to be poor indicators of species distinction (Rakotoarison et al., 2012;Rojas et al., 2016). Advertisement call variation is mainly based on differences in dominant frequency, again an indicator of size of the calling individual (Fig. S1) and thus, of low taxonomic significance. The features of qualitative external morphology proposed by Rosset et al. (2007) to diagnose O. achalensis, O. barrioi and O. occidentalis represent the extremes of a continuum between lowland and highland ecotypes, and between eastern and northern variants of the same species. For the same reason González et al. (2014) failed to detect significant morphological differences among the tadpoles within the O. occidentalis group refuting the claim of Cei & Crespo (1982). Moreover, defensive behavior of adults is also indistinguishable (Borteiro et al., 2018). Thus, the phenotypic peculiarities of O. barrioi seem responses to local environmental conditions rather than indicating taxonomically relevant information. Applying the Principle of Priority of the International Code of Zoological Nomenclature (International Commission on   (Fig. 5). The basal splitting of lineages separates O. occidentalis, the only species with complex advertisement call consisting of several pulse groups, from the two lineages with simple pulsed calls. The ancestral character state of advertisement call structure in Odontophrynidae is most likely a simple pulsed call, as also present in the sister group Macrogenioglottus alipioi (Abravaya & Jackson, 1978) and in most Proceratophrys. The species occurring in Argentina and Bolivia, the diploids O. cordobae and O. lavillai, and the topotypical tetraploids O. americanus are closely related, but the sister species relationship between O. americanus and O. cordobae is well resolved possibly indicating an autopolyploid origin of these tetraploids.
In conclusion, we verify hypothesis 2 with respect to the genetic base of the phenetic groups. Our reconstruction of phylogenetic relationships among these groups suggests that O. occidentalis evolved from the ancestral stock before the diversification of the O. americanus and O. cultripes group occurred.
Hypothesis 3: The current Red List classification does not reflect the risks upon each species The geographical distribution of O. occidentalis is larger than previously appreciated extending to north (barrioi ecotype) and to east (achalensis ecotype) (Fig. 6). Odontophrynus occidentalis is endemic to Argentina inhabiting many localities in eight provinces covering about 16% of the territory. This species is highly adaptable to a wide altitudinal range, and tolerant to local sympatry with O. americanus and O. cordobae. Odontophrynus aff. americanus (pale blue label). Misiones province; probably identical with an undescribed americanus-like species from Brazil. Details on localities are given in Table S1.
Full-size  DOI: 10.7717/peerj.6480/ fig-6 Thus, the Red List classification "Least Concern" seems justified, whereas the associated ecotypes "achalensis" and "barrioi" ("Vulnerable" and "Data Deficient") do not deserve classifications apart. With respect to the species referred to as O. americanus our study suggests strongly that there is more than one species involved. The nominal species O. americanus is certainly widespread in Argentina (16 provinces and ca. 67% of the territory) and extends to Bolivia and Paraguay in the north, and to Uruguay and southern Brazil in the east. The status "Least Concern" seems appropriate. The exact range of this taxon in Brazil remains to be assessed using barcoding for species identification. Most probably, the easternmost locality in Misiones pertains rather to O. aff. americanus 2 of Brazil than to the nominal taxon of Argentina. If this assumption proves true, the still undescribed Odontophrynus taxon would constitute the fifth species of this genus in Argentina. O. cordobae has the smallest area of distribution of the four species, occurs exclusively in the central part of the Cordoba province, and thus, is endemic to Argentina (Fig. 6). Recent assessment of localities inhabited demonstrates that there is a viable network of probably connected populations (Grenat et al., 2018). Therefore, we propose the classification "Least Concern" as long as there is no further shrinkage of its geographical range. Finally, O. lavillai inhabits eight provinces of Argentina as does O. occidentalis, but its range extends further north to Bolivia and Paraguay (Rosset & Baldo, 2014). The classification "Least Concern" seems reasonable for this species as well. The Red List status of newly described species from Brazil and Uruguay and those of the O. cultripes group are outside the scope of this study.
In conclusion, we falsify hypothesis 3 with respect to genetically delimited Odontophrynus species of Argentina.

CONCLUSIONS
Integrative taxonomy has proved to be the appropriate tool to cope with distinct levels of character differentiation in the morphologically highly conserved genus of Odontophrynus toads. Genotypic variation among the nominal taxa of the O. occidentalis group did not correspond to the phenotypic plasticity in response to altitudinal and latitudinal gradients found in the ecotypes "achalensis", "barrioi" and "occidentalis". Consequently, molecular evidence melts down the O. occidentalis group to a single, polymorphic species O. occidentalis. Whereas the number of species was grossly overestimated in this case, considerable genetic divergence between O. americanus originating from a topotypical population (Argentina) and specimens of unknown ploidy from Brazil indicates that specimens referred to as O. americanus in Argentina and in Brazil are distinct species except for the southernmost populations in Brazil. Phylogenetic relationships among Odontophrynus species suggest that O. occidentalis evolved from the ancestral stock before the diversification of the O. americanus and O. cultripes group occurred. Reliable taxonomic delimitation of Odontophrynus taxa allows for a precise assessment of the corresponding geographical ranges and for an informed basis of the Red List classification. The four species occurring in Argentina do not seem endangered currently, but the small geographic range of O. cordobae may require a future reassessment of the species' status.