Unrevealing the leaf frogs Cerrado diversity: A new species of Pithecopus (Anura, Arboranae, Phyllomedusidae) from the Mato Grosso state, Brazil

The Neotropical frog genus Pithecopus comprises currently 10 species. A recent molecular phylogeny suggested the existence of two subclades within it, one of them including P. palliatus, P. azureus, P. hypochondrialis, and P. nordestinus (lowland species). Herein we describe a new species of this subclade from Pontal do Araguaia, in the Brazilian Cerrado in the Mato Grosso state. Recognition of the new species is supported by adult morphology, advertisement call and molecular data. The new species differs from Pithecopus highland species by its smaller head width and lack of the reticulate pattern on flanks. From lowland species, the new form differs by being significantly smaller in snout vent-length, advertisement call with the greatest number of pulses, and high genetic distance. Interestingly, we also report on occurrence of P. hypochondrialis (its sister species) at an adjacent site (about 3km). Also, we report on the occurrence of the new species in the Chapada dos Guimarães and Santa Terezinha, both also in the Mato Grosso state.


Introduction
The increase of the fieldwork sampling efforts and the use of multiple databases in taxonomical studies open the doors to the fascinating biological diversity present in the central region of Brazil [1]. The Cerrado domain occupies the largest portion of central Brazil, and is composed by seasonal savannas with corridors of mesic gallery forests [2]. The increasing in anuran species richness recognized to Cerrado have been grown fast in the last decades and pointed to that species composition of this area need to better documentation [1]. In contrast, the a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 Types and/or vouchers are deposited in the Museu de Zoologia (ZUEC), at Universidade Estadual de Campinas (UNICAMP); Célio F. B. Haddad amphibian collection (CFBH) at the Universidade Estadual Paulista (UNESP), Rio Claro, both in São Paulo state; Museu Nacional do Rio de Janeiro (MNRJ), at Universidade Federal do Rio de Janeiro (UFRJ), Rio de Janeiro, Rio de Janeiro state; and in collection of frogs of the Museu de Biodiversidade do Cerrado (AAG-UFU), Universidade Federal de Uberlândia (UFU), Uberlândia, Minas Gerais state, all in Brazil. Sixteen adult males from Chapada dos Guimarães and three adult males from Santa Terezinha (Fig 1) of the new species were also measured as additional specimens. The complete list of vouchers examined is available on S1 File (Appendix A).
Because the new species most resembles P. hypochondrialis and this taxa have a large geographical distribution in Brazil with morphological variation already reported [5,7], we considered a large number of specimens from several localities in our analyses. We treated separately specimens of P. hypochondrialis from Pará and Amapá as "P. hypochondrialis North" because their proximity to Suriname (type-locality of P. hypochondrialis) and the specimens from Minas Gerais, Goiás, Mato Grosso, and Tocantins as "P. hypochondrialis South", according to genetic differences pointed out by Bruschi et al. [7] and intrageneric relations recovered in our mitochondrial topology (see results section). Our sample of "P. hypochondrialis South" contains specimens, calls and mtDNA sequences from Barra do Garças, Mato Grosso state (Fig 1), a neighboring locality to the type-locality of the new form. Details for examined specimens in S1 File (Appendix A).

Morphometry
Morphometric traits of adult males were measured using a Mitutoyo Absolute digital caliper CD-6" CSX to the nearest 0.1 mm, except for finger and toe discs which were measured under a stereomicroscope (Zeiss Stemi 2000) coupled to an ocular micrometer.
The calls of some species of Pithecopus often have a few weak isolated pulses at the end (more common) or at the beginning (unusual) (see results section). To compare call traits equivalently among specimens/species, we used the main (stronger) group of pulses of each call (hereafter called "core portion") of the new species, P. hypochondrialis, and P. nordestinus. Traits such as pulse duration and pulses per second were measured within the core portion (further details in S1 Table-Appendix B). The presence or absence of isolated pulses were recorded for each individual; when they were absent, we state it as "0", when at least one was present, as "1", and two pulses as "2". We calculated means and standard deviations (SD) for each individual and then the overall mean and SD was calculated based on those values, whereas the range encompassed the minimum and maximum values for the whole sample. We analyzed calls using Raven Pro 1.5, 64 bit version [16] with the following settings: window type = Hann, window size = 256 samples, 3 dB filter bandwidth = 270 Hz, brightness = 50%, contrast = 50%, overlap = 85% (locked), color map = "Cool", DFT size = 1024 samples (locked), grid spacing (spectral resolution) = 46.9 Hz. We analyzed temporal traits in oscillograms and spectral traits in spectrograms [15]; we used the 'Peak Frequency' function to determine the peak of dominant frequency. We generate sound figures through package Seewave v.1.6 package [17], R platform (version 3.3.1) [18]. Seewave settings for the spectrograms were: Hanning window, 85% overlap, and 256 points resolution (FFT). Analyzed sound files are listed in S2 Table (Appendix C).

Statistical analysis
Considering the morphological and acoustic (multivariate) datasets separately, we sought for discrimination between populations/species by applying two functions: (1) randomForest (RF) (radomForest v. 4.6-12 package) [19] and (2) [20,21]. RandomForest algorithm constructs many (e.g. 500) classification trees using bootstrap samples of the data (each split using the best predictors among those randomly chosen at each node) then generating classifiers and aggregating results by voting to classes [19]. The classic Discriminant Analysis (DA) depends on multivariate normality [22] and on a larger number of objects than variables. The multivariate normality of the original data was evaluated through the function mardiaTest (MVN package) [23]. The DAPC performs analyses on the Principal Component scores [20,21]. The application of a DA on a few axes (preserving about 95% of the variance) of a Principal Component Analysis, as performed by DAPC, improves the imbalance between objects and traits [21]. Despite the lack of normality in both of our datasets (details not shown), the results of DAPC are presented within an exploratory context to assess the congruence between it and random-Forest discriminations. The direct or indirect packages related the application of both discriminant functions were run in R.
For the multivariate analysis and statistical tests we used all the morphometric features detailed above. In addition to the types, additional specimens of the new form from Chapada dos Guimarães and Santa Terezinha (MT) were also considered. For the acoustic analyses we used: call duration (entire call), number of pulses, pulse duration, interpulse interval within core, core duration (main stronger group of pulses), pulses per core, pulses per second, isolated pulse and peak of dominant frequency. Considering that both multivariate analyses, to both datasets, were concordant in species discrimination (see Results section) we present the RF classification results in tables and DAPC in scatter plots.
The acoustic and morphometric traits were tested for statistical significance of the differences among population/species through the Exact Wilcoxon Mann Whitney Rank Sum Test, function wilcox_test of the package Coin (Resampling Statistics model) [24] in R. As these tests were done between species/populations pairs, the significance levels ("p") were adjusted considering the number of pairings through the method of Holm (p.adust function in R).

Molecular divergences
To evaluate the reciprocal monophyly of new species respect to the sister species/taxon, we performed phylogenetic analyses of mitochondrial 16S ribosomal sequences (1490 characters) of paratopotypes (ZUEC 21657, 21659 and 21660), and additional specimens (CFBHt 4638, 4684, 4687) of the new Pithecopus species. The data matrix was completed with sequences of at least one individual of each of the Phyllomedusidae species available in GenBank database. We also included in our dataset sequences of the individuals from Pithecopus from Pará and Amapá states ("P. hypochondrialis North") and the specimens from Minas Gerais, Maranhão, Goiás and Tocantins states ("P. hypochondrialis South"), according the topology of Bruschi et al. [7]. We chose Agalychnis granulosa as the outgroup based on Faivovich et al. [5]. A list of species, sequences and GenBank accession numbers is provided in S3 Table. Genomic DNA was extracted from liver or muscle tissue using the TNES method, as applied by Bruschi et al. [25]. The tissues are stored at -70˚C in the tissue bank of the Departamento de Biologia Estrutural e Funcional at UNICAMP, São Paulo state, Brazil. The mitochondrial 16S ribosomal genes were amplified using the primers MVZ 59(L), MVZ 50(H), 12L13, Titus I (H), Hedges16L2a, Hedges16H10, 16Sar-L, and 16Sbr-H (for primer sequences, see [26]. The amplified PCR products were purified using Exonuclease I (10 units) and SAP (1 unit), with a 45-min incubation at 37˚C and a 10-min denaturation at 85˚C, then used directly as templates for sequencing in an automatic ABI/Prism DNA sequencer (Applied Biosystems, Foster City, CA, USA) with the BigDye Terminator kit (Applied Biosystems, Foster City, CA, USA), as recommended by the manufacturer. The DNA samples were sequenced bidirectionally and edited in the software CodonCode Aligner 3.7.1 (Codon Code Corporation, Dedham, MA, USA).
We aligned sequences using CLUSTALW [27] implemented in CodonCode Aligner 4.0. We used three different gap penalties (5, 10, and 15) for each gene to identify potential sites of ambiguous homology [28]. Gap length was kept constant (0.20) and all other parameters were set at default settings. The dataset was used to phylogenetic reconstruction by Bayesian inferences (BI) and Maximum Parsimony criterion.
Bayesian inference was based on a Markov chain Monte Carlo (MCMC) analysis conducted in MrBayes 3.1.2 [29] with two independent runs, each with four chains and sampling every 1000 generations for six million generations. An adequate burn-in (the first 25% trees were excluded) was determined by examining a plot of the likelihood scores of the heated chain for convergence and stationarity. The evolutionary model most appropriate (GTR+R+I) was selected by MrMODELTEST [30] using the Akaike Information Criterion (AIC). The trees were sampled every 100 generations, excluding the first 25% of the trees as burn-in, determined by examining a plot of the likelihood scores the heated chain for convergence and stationarity. Tracer software version 1.5 [31] was used to confirm the quality of the parameters of the Bayesian inferences. To the Maximum Parsimony criterion was implemented in TNT v1.1 software [32] using a heuristic search method with tree bisection-reconnection (TBR) swapping and 100 random additional replicates. The bootstrap values of the branches inferred in this analysis were calculated with 1000 non-parametric pseudoreplicates.
Uncorrected genetic distances (p-distances) among the sequences of the Pithecopus species were calculated using the maximum composite likelihood model [33] implemented in MEGA5 [34]. Gaps and missing data were eliminated in this analysis and all parameters were left as in the default settings.
To evaluate if geographic distances could explain the genetic differences among populations from Pontal do Araguaia + Chapada dos Guimarães and Barra do Garças, we performed a Mantel test in Alleles In Spaces (AIS) software. For identifying biogeographic boundaries or areas with largest genetic differences among samples collected in the Mato Grosso state and allocated in two different clades in our phylogenetic inferences, we compute the Monmonier's algorithm of maximum-differences [35] in Alleles In Spaces (AIS) software.

Nomenclatural acts
The electronic edition of this article conforms to the requirements of the amended International Code of Zoological Nomenclature, and hence the new names contained herein are available under that Code from the electronic edition of this article. This published work and the nomenclatural acts it contains have been registered in ZooBank, the online registration system for the ICZN. The ZooBank LSIDs (Life Science Identifiers) can be resolved and the associated information viewed through any standard web browser by appending the LSID to the prefix "http://zoobank.org/". The LSID for this publication is: urn:lsid:zoobank.org:pub:00E56B5F-1EBD-47C8-8A54-2EC50EAD70A9. The electronic edition of this work was published in a journal with an ISSN, and has been archived and is available from the following digital repositories: PubMed Central, LOCKSS.    Head length 6.8 ± 0.5 (6.1-
From its closer relatives (the lower land species), the new species is significantly smaller than P. nordestinus, P. azureus and P. hypochondrialis "North" and "South" in snout-vent length, head length, head width, axilla-groin length, eye-nostril distance, internarial distance, hand length, thigh length, tibia length and foot length (Exact Wilcoxon-Mann-Whitney Test: p < 0.01). The randomForest model on morphometric variables classified 35 individuals of P. araguaius sp. n. correctly (Table 3), only one male (the largest paratopotype AAG-UFU 3449; SVL = 33.8 mm) was misclassified. However, this same male is a call voucher specimen and was classified within new species group in acoustic multivariate analyzes (see results below). Accordingly, the DAPC discriminated the new species, with a greater separation along axis 1 (LD1 = 55%), but the axis 2 (LD2 = 28%) also contribute to separation. Tarsus (32%), tibia   (17%), and thigh (14%) lengths mainly accounted for species separation along LD1 ( Fig 5A); while eye-nostril (29%) and internarial (19%) distances, and head (10%) and hand (10%) lengths along LD2 (Fig 5A). The new species had the lowest values for all these variables along the axes 1 and 2, all Exact Wilcoxon-Mann-Whitney Tests with p < 0.01. The specimen AAG-UFU 3101 (SVL = 32.2 mm) of P. hypochondrialis South from Araguari (MG) is the smallest male for this group, due to this it was classified within new species cluster. Raw morphological measurement data of all examined specimens in S4 Table. In comparison with the lowland species, the new species differs from Pithecopus palliatus by having a higher dominant frequency (2239-3316 vs. 1580 Hz) and by having a single note call, and not release as double notes as in P. palliatus [41]. The new species can be differentiated from P. hypochondrialis North and South, P. azureus and P. nordestinus by having a higher peak of dominant frequency (p <0.01, < 0.01, 0.01, and 0.01; respectively), a higher number of pulses per call (p <0.01, <0.01, 0.01, and 0.01), a longer core duration (p <0.01, <0.01, <0.01, and 0.01), and a higher number of pulses per core (p <0.01, <0.01, 0.01, and 0.01). From P. hypochondrialis North, the new species can also be differentiated by having a longer call duration (p < 0.01). In addition, P. araguaius sp. n. is differentiated from P. hypochondrialis South by its higher pulse rate (p < 0.01).
Description of the holotype. General aspect slender (Fig 2A and 2B); snout truncate in lateral ( Fig 3A) and dorsal view (Fig 3B). Head wider than long; loreal region slightly concave; canthus rostralis rounded; nostrils small, subcanthal, placed latero-frontally, closer to snout tip than to eyes; internarial distance longer than eye-nostril distance and tympanum diameter, but smaller than eye diameter; eyes lateral-frontally positioned; tympanum nearly circular, with annuli undefined at the superior border; tympanum diameter less than half eye diameter; supratympanic dermal fold present, beginning on top of tympanum and ending nearly at the mouth corner; dorsolateral macrogland undistinguished; no vocal slits or external vocal sac; tongue nearly ovoid, free posteriorly, longer than wide; rudimentary vomerine teeth present, bordering choanae internally; choanae small, located laterally, slight rounded.
Upper arm thin and forearm robust; no finger webbing; comparative finger length when adpressed I<II<IV<III; finger discs poorly developed, smaller than the tympanum diameter; finger I enlarged at the base; nuptial asperity covering most of the dorsal surface of finger I, except the tip; palmar tubercles poorly developed (Fig 3C), subarticular tubercles developed, rounded and distinct from the supernumerary ones; inner and outer metacarpal tubercles undifferentiated; comparative toe length when adpressed II<III<I<V<IV; no toe webbing; plantar callosities poorly developed (Fig 3D), inner and outer metatarsal tubercles undifferentiated; subarticular tubercle developed, single and rounded; supernumerary tubercles rounded and poorly developed; legs slender, thigh slightly longer than tibia; heel reaches the posterior border of tympanum when a leg is adpressed to body; dorsal skin smooth; ventral skin granulated on belly and throat and smooth on thigh and foot. Left thigh slightly damaged ventrally due to tissue sampling. Color in life (n = 17 male paratopotypes). Dorsal surfaces of forearm, tibia, tarsus and portions of foot and hand green (Fig 4A and 4B). A green stripe on dorsal surface of the thigh, which never reaches the final proximal thigh portion (Fig 4). Vertical dark stripe on flanks, dorsal surface of arm, medial surface of forearm and hidden part of the leg orange (Fig 4). Border of upper lip and upper eyelid with a white line that, on lip, never reaches the border of lower eyelid. Thin dark lines surrounding the upper surface of the fingers and toes (Fig 4C and 4D). In preservative the green areas become pale bluish gray and the orange parts become light cream or whitish and the dark stripes, upper lip and upper eyelid lines remain unchanged (Figs  2 and 3). All type specimens have a thin white line at the edge of the upper lip which reaches the border of lower eyelid and the incomplete green strip in the dorsal face of the thigh (Fig 4).

Measurements of holotype (mm, abbreviations as in M&M
Variation. Size variation in Table 1 Etymology. The epithet araguaius it is masculine latinized form of the indigenous Tupi word "araguaia", a reference to the Araguaia River, which cross the type-locality of the new species. Advertisement call. See also S1 Table (Appendix B) for further call traits definitions. Seven males of the new species were recorded, 54 calls analyzed. Quantitative traits are summarized in Table 2. The advertisement call of Pithecopus araguaius sp. n. consists of a single pulsed note emitted sporadically (Fig 6). Generally, the notes are composed of a long and strong group of pulses (core), and rarely by isolated weak final pulses. The core portion lasted 28-48 ms (mean = 39.0; SD = 5.0; n = 54) and have 5-8 pulses per core (mean = 6.0; SD = 0.5; n = 54). Also considering the isolated pulses, notes lasted 28-63 ms (mean = 41.0, SD = 4.9; n = 54) and have 5-8 pulses (mean = 6.0, SD = 0.6, n = 319). When present (only one male), isolated pulses were limited to one and lasted 2-5 ms (mean = 3.5; n = 5). The interval between the core and isolated pulses ( Table 2) varied from 5-8 ms (mean = 6.0; n = 5). The pulse duration varied from 2-17 ms (mean = 5.1, SD = 1.3, n = 398), emitted at rates of 114-206 pulses/ second (mean = 155.0, SD = 20.0, n = 54). The arrangement of pulses had the following patterns: (1) core with six pulses and no isolated pulse (41%; n = 22 calls, noticed in all 7 recorded males; see Fig 6A); (2) core with five pulses and no isolated pulse (33%; n = 18 calls, noticed in 4 males); (3) core with seven pulses and no isolated pulse (13%; n = 7 calls, in 4 males); (4) core with six pulses and one isolated pulse (6%; n = 3 calls, in 1 male); (5) core with seven pulses followed by one isolated pulse (see Fig 6B) and (6) core with eight pulses and no isolated pulse also was found and had the same frequency as other calls (4%; n = 2 calls, in 1 male). The peak of dominant frequency varied from 2240-3316 Hz (mean = 2540 Hz, SD = 308; n = 54).
In our sample, the advertisement call of the sister species of the new taxon, Pithecopus hypochondrialis South (33 males, 180 calls), consists of a short note, with sporadic emission (quantitative call traits are summarized in Table 2; Fig 7). The core duration varied from 19-61 ms (mean = 33.2, SD = 3.2; n = 180) and the intervals between core and isolated pulses varied from 5-21 ms (mean = 12.5, SD = 3.4; n = 64). The number of pulses per core varied from 3-5 pulses (mean = 4.0, SD = 0.2; n = 180) ( Fig 6); with intervals (or no interval) within core pulses from 1-7 ms (mean = 2.1, SD = 0.9; n = 539). Pulses were arranged in (1) core with four pulses and no isolated pulse (53%; n = 95 calls) (see Fig 7A and Fig 7B); (2) core with four pulses followed by an isolated pulse (31%; n = 55 calls); (3) core with five pulses and no isolated pulse (7%; n = 13 calls); (4) core with three pulses and no isolated pulse (4%; n = 8 calls); (5) core with three pulses followed by one isolated pulse (4%; n = 7 calls) and (6)  Phylogenetic relationships. The topologies inferred by parsimony criterion and Bayesian inference were congruent one each other and similar to the phylogenetic relations recovered by Faivovich et al. [5], Bruschi et al. [7,39] and Duellman et al. [4]. In all our analyses, the paratopotypes of the P. araguaius sp. n. were nested with additional specimens from Chapada dos Guimarães and Santa Terezinha (both MT), forming a highly supported monophyletic group in both phylogenetics methods (Fig 8). This clade is the sister group of a clade composed by specimens of P. hypochondrialis North and South.
Specimens from the adjacent municipality of Barra do Garças (Fig 1) were recovered within P. hypochondrialis clade, nested with P. hypochondrialis South populations (specimens AAG-UFU 1083 and 1084; see topology in Fig 8). The genetic distances between P. araguaius sp. n. and other Pithecopus species ranged from 4 to 9% (Table 5); within the lowland species the genetic divergence levels varied from 4 to 7% (Table 5). Besides the high genetic distance (Table 5) and distinct phylogenetic position (Fig 8), morphometric and bioacoustics data of P. hypochondrialis from Barra do Garças are in our morphometric and acoustic multivariate analyses, and in both cases, they were not recovered within the clusters of the new species (Fig 5).
Distribution. Based on the morphological and genetic similarities between populations from Pontal do Araguaia, Chapada dos Guimarães and Santa Terezinha (all MT) [7], the specific identity between all three is well supported. Chapada dos Guimarães and Santa Terezinha are about 380 km west and 630 km north from the type-locality of P. araguaius sp. n., respectively. Based on the distribution provided by Bruschi et al. [7], the range of the new species overlaps with that of P. hypochondrialis South, but not with those of P. azureus and P. nordestinus (see their Fig 2).
The genetic distance between the topotypes and specimens from Chapada dos Guimarães which are around 400 km apart shown low levels of divergence (P-distance 1%) in comparison with the P-distance of the 4% between individuals from Barra dos Garças (Table 5). The Mantel test scores (r = 0,083) (S2 File) formally reveled that geographic distances did not plays an important role in the distribution of genetic variation detected among samples and represent a strong evidence of the complete genetic isolation between P. hypochondrialis South and P. araguaius sp. n.. Monmonier algorithm detected biogeographic boundaries between populations from Pontal do Araguaia + Chapada dos Guimães versus Barra dos Garças that could be limited the gene flow among populations and represent putative origin of high genetic variation detected, however what are these biogeographical boundaries should be carefully evaluated in the future (S2 File).

Discussion
A large amount of Brazilian Anuran hidden diversity has been revealed by the application of molecular genetics and its integration into taxonomic researches [39,42,43]. Phenotypic and genotypic evidences allowed the recognition of the new species here described. Pithecopus araguaius sp. n. have an evolutionary history related to the remaining lowland species of the Brazilian Cerrado and possibly has been historically confused to P. hypochondrialis due the taxonomic problems related to its definition and sympatry. Particularly, the inclusion of P. hypochondrialis samples from Barra dos Garças in our analyses was a crucial evidence to postulate P. araguaius sp. n. as a valid new species, as evidenced by the high level of genetic distance between them (noteworthy is that P. hypochondrialis AAG-UFU 1083 is a molecular, acoustic and morphometric voucher).
Indeed, the routine morphological identification of the Pithecopus species is not easy, mostly due to the high level of conservative morphology observed among P. hypochondrialis, P. azureus and P. nordestinus. Even to experienced herpetologists it has been difficult to separate these taxa when compared to the "nomimal species", what has led to a repeated history of nomenclature changes in this group [8]. For example, the diagnostic characters proposed by Caramaschi [8] to distinguish P. hypochondrialis from P. azureus are too variable within and among populations [7] and in some cases has failed to allow its discrimination. Ours findings also re-emphase the notion that the species richness within lowland Pithecopus species was underestimated.
Haga et al. [10] pointed that the similarities observed among the advertisement calls of the P. hypochondrialis, P. azureus and P. nordestinus, and the resulting difficulty in discriminating them, indicates that acoustic characters are uninformative for these closely related species. However, these same authors emphasized that molecular and cytogenetic datasets supported their independent specific identities, as suggested by Bruschi et al. [7]. Our results corroborated those highlighted by Haga et al. [10], once Pithecopus araguaius sp. n. was the only species with complete discrimination in the randomForest model and DAPC on acoustic dataset. It is noteworthy that the P. azureus and P. nordestinus (sister species) are morphometrically indistinguishable from one another according to these same analyzes. Therefore, even though the acoustics and morphometric differences found in our statistical analyses may be deemed unable to fully diagnose the new species from the other three species abovementioned, due to the slight overlaps in their traits, they provide sufficient evidence to support our hypothesis that this lineage is evolving independently.