Evolutionary radiation of earless frogs in the Andes: molecular phylogenetics and habitat shifts in high-elevation terrestrial breeding frogs

Study area and field surveys

We surveyed multiple montane sites between April 2012 and September 2014, and collected specimens and ecological data on 11 species of Phrynopus. The surveyed areas included the type locality of Phrynopus badius Lehr, Moravec & Cusi, 2012, Phrynopus curator Lehr, Moravec & Cusi, 2012, Phrynopus daemon Chávez et al., 2015, Phrynopus interstinctus Lehr & Oróz, 2012, Phrynopus inti Lehr et al., 2017, Phrynopus juninensis (Shreve, 1938), Phrynopus montium (Shreve, 1938), Phrynopus peruanus Peters, 1873, Phrynopus unchog Lehr & Rodríguez, 2015, Phrynopus vestigiatus Lehr & Oróz, 2012, Phrynopus spI (R von May, 2018, unpublished data). In addition to collecting primary data on these species, we reviewed published studies to assemble a dataset containing minimum and maximum elevation, elevational midpoint, and habitat (Table S1) for those species included in phylogenetic comparative analyses (see below). Research and collecting permits were issued by the Dirección General Forestal y de Fauna Silvestre (DGFFS) and the Servicio Nacional Forestal y de Fauna Silvestre (120-2012-AG-DGFFS-DGEFFS, 064-2013-AG-DGFFS-DGEFFS, 292-2014-AG-DGFFS-DGEFFS, R.D.G. No 029-2016-SERFOR-DGGSPFFS, R.D.G. 405-2016-SERFOR-DGGSPFFS, and Contrato de Acceso Marco a Recursos Genéticos, No 359-2013-MINAGRI-DGFFS-DGEFFS) and the Servicio Nacional de Areas Naturales Protegidas (No 001-2012-SERNANP-JEF). Export permits were issued by the Ministerio del Ambiente, Lima, Peru. Use of vertebrate animals was approved by the Animal Care and Use Committee of the University of California (ACUC #R278-0412, R278-0413, and R278-0314). Specimens collected during this study were deposited in the herpetological collection at the Museo de Historia Natural Universidad Nacional Mayor de San Marcos (MUSM), Lima, Peru.

Morphological data

We compiled available morphological data for 21 species of Phrynopus included in our phylogenetic analysis (data were not available for Phrynopus sp, an undescribed species from Pasco region, central Peru; Lehr et al., 2017). Most of these data were extracted from the original publications describing these species (Hedges, 1990; Lehr, Köhler & Ponce, 2000; Lehr, 2001; Lehr & Aguilar, 2002; Lehr & Aguilar, 2003; Lehr, Lundberg & Aguilar, 2005; Lehr, Fritzsch & Müller, 2005; Chaparro, Padial & De la Riva, 2008; Duellman & Hedges, 2008; Lehr & Oróz, 2012; Lehr, Moravec & Cusi, 2012; Chávez et al., 2015; Lehr & Rodriguez, 2017; Lehr & von May, 2017; Rodríguez & Catenazzi, 2017). Additionally, for three species (P. juninensis, P. peruanus, and Phrynopus spI), we took measurements of specimens collected in recent expeditions. Our dataset included only adult individuals (N = 76); for all of these, sex and maturity of specimens were determined by observing gonads through dissections under a stereomicroscope (Table S1; see original species descriptions for further detail). Our dataset included the following measurements: snout–vent length (SVL), tibia length (TL, distance from the knee to the distal end of the tibia), foot length (FL, distance from proximal margin of inner metatarsal tubercle to tip of Toe IV), head length (HL, from angle of jaw to tip of snout), head width (HW, at level of angle of jaw), horizontal eye diameter (ED), interorbital distance (IOD), upper eyelid width (EW), internarial distance (IND), eye–nostril distance (E-N, straight line distance between anterior corner of orbit and posterior margin of external nares). Information on species tympanum condition was obtained from Duellman & Lehr (2009) and from subsequent species descriptions (Lehr & Oróz, 2012; Lehr, Moravec & Cusi, 2012; Chávez et al., 2015; Lehr & Rodriguez, 2017; Lehr et al., 2017; Rodríguez & Catenazzi, 2017).

Micro-computed tomography

To determine the condition of the columella, a key element of the tympanic middle ear, we obtained X-ray µCT images for one eared and one earless species, Phrynopus peruanus and Phrynopus montium, respectively. Two voucher specimens stored in ethanol were scanned in the Micro-CT Core facility at the University of Michigan. Specimens were placed in a 34 mm diameter specimen holder and scanned using a microCT system (µCT100; Scanco Medical, Bassersdorf, Switzerland). Scan settings were: voxel size 11.4 µm, 55 kVp, 145 µA, 0.5 mm AL filter, 1,000 projections around 180°, integration time of 1,000 ms and 3 average data. Data were exported to DICOM files using Scanco’s proprietary software. We obtained images after segmentation using isosurface representations and three-dimensional renderings using the Amira-Avizo software. According to the predictions of Pereyra et al. (2016) the presence of tympanic membrane and tympanic annulus in P. peruanus implies that the columella is present. Thus, we expected to observe the columella in a CT scan of P. peruanus. In contrast, the absence of both tympanic membrane and tympanic annulus in P. montium did not imply that the columella would be absent. Thus, CT scanning this specimen was particularly critical to determine if the columella was present (or not) in this species.

Molecular phylogenetic analysis

Our analysis included DNA sequence data from Phrynopus species that were available in GenBank (as of 1 August 2017; Table S2) as well as sequences from other closely related genera (Lynchius, Oreobates) and more distantly related ones (Ischnocnema guentheri, Hypodactylus brunneus, and Hypodactylus dolops) as outgroups following the results of Padial, Grant & Frost (2014). Additionally, newly produced sequences include those obtained from 11 species: P. badius, P. curator, P. daemon, P. interstinctus, P. juninensis, P. montium, P. peruanus, P. unchog, P. vestigiatus, Phrynopus inti, and Phrynopus spI. A notable addition to this phylogeny is P. peruanus, the type species of the genus.

The mitochondrial genes were a section of the 16S rRNA gene, a section of the 12S rRNA gene, and the protein-coding gene cytochrome c oxidase subunit I (COI). The nuclear genes were the recombination-activating protein 1 (RAG1) and Tyrosinase precursor (Tyr). Extraction, amplification, and sequencing of DNA followed protocols previously used for Neotropical terrestrial breeding frogs (Lehr, Fritzsch & Müller, 2005; Hedges, Duellman & Heinicke, 2008). Primers used are listed in Table S3. We employed the following thermocycling conditions to amplify DNA from each gene using the polymerase chain reaction (PCR). For 16S, we used: 1 cycle of 96 °C/3 min; 35 cycles of 95 °C/30 s, 55 °C/45 s, 72 °C/1.5 min; 1 cycle 72 °C/7 min. For 12S, we used: 1 cycle of 94 °C/1.5 min; 35 cycles of 94 °C/45 s, 50 °C/1 min., 74 °C/2 min; 1 cycle 72 °C/10 min. For COI, we used: 1 cycle of 96 °C/3 min; 35 cycles of 95 °C/30 s, 55 °C/45 s, 72 °C/1.5 min; 1 cycle 72 °C/7 min. For RAG1, we used: 1 cycle of 96 °C/2 min; 40 cycles of 94 °C/30 s, 52 °C/30 s, 72 °C/1.5 min; 1 cycle 72 °C/7 min. For Tyr, we used: 1 cycle of 94 °C/5 min; 40 cycles of 94 °C/30 s, 54 °C/30 s, 72 °C/1 min; 1 cycle 72 °C/7 min. We completed the cycle sequencing reactions by using the corresponding PCR primers and the BigDye Terminator 3.1 (Applied Biosystems, Foster City, CA, USA), and obtained sequence data by running the purified reaction products in an ABI 3730 Sequence Analyzer (Applied Biosystems, Foster City, CA, USA). We deposited the newly obtained sequences in GenBank (Table S2).

We used Geneious R6, version 6.1.8 (Biomatters, 2013; http://www.geneious.com/) to align the sequences with the built-in multiple alignment program. Prior to conducting phylogenetic analysis, we used PartitionFinder, version 1.1.1 (Lanfear et al. 2012) to select the appropriate models of nucleotide evolution and used the Bayesian information criterion (BIC) to determine the best partitioning scheme and substitution model for each gene. We divided the partitions a priori; each mitochondrial gene was considered a data block whereas each nuclear gene was divided by codon position (i.e., three data blocks per gene). The best partitioning schemes and models of evolution are included in the Results section.

We employed a Bayesian approach using MrBayes, version 3.2.0 (Ronquist & Huelsenbeck, 2003) to infer a molecular phylogeny. Our analysis included 72 terminals and a 2,646-bp concatenated partitioned dataset. We performed an MCMC Bayesian analysis that consisted of two simultaneous runs of eight million generations, and we set the sampling rate to be once every 1,000 generations. Each run had three heated chains and one “cold” chain, and the burn-in was set to discard the first 25% samples from the cold chain. Following the completion of the analysis, we used Tracer 1.6 (Rambaut & Drummond, 2007) to verify convergence. Subsequently, we used FigTree (http://tree.bio.ed.ac.uk/software/figtree/) to visualize the majority-rule consensus tree and the posterior probability values to assess node support.

We also used a multispecies coalescent approach implemented in *BEAST v1.6.2 (Drummond & Rambaut, 2007) to infer a species tree to be used for phylogenetic comparative analyses (see below). Our analyses only depend on the relative branch lengths of the tree, but we preferred to illustrate our tree in approximate units of time. Therefore, we used an uncorrelated relaxed molecular clock with the rate of nucleotide substitution for 16S set at 1% per million years as done in a recent study (von May et al., 2017). However, as in von May et al. (2017), we note that the dates associated with the tree should only be viewed as very approximate and that there can be multiple sources of error when calibrating phylogenies (Arbogast et al., 2002). The analysis in *BEAST included two independent runs, each with 100 million generations and sampled every 10,000 generations. Following the completion of the analysis, we used Tracer v1.5 (Rambaut & Drummond, 2007) to examine effective sample sizes, verify convergence of the runs, and to ensure the runs had reached stationarity. Observed effective sample sizes were sufficient for most parameters (ESS > 200) except for substitution rates for a few partitions. We discarded the first 10% of samples from each run as burn-in. Subsequently, we used LogCombiner v1.6.2 to merge all remaining trees from both runs and used TreeAnnotator v1.6.2 (Drummond & Rambaut, 2007) to summarize trees and obtain a Maximum Clade Credibility tree. We visualized the Maximum Clade Credibility (MCC) tree and the associated node support values in FigTree (http://tree.bio.ed.ac.uk/software/figtree/).

Morphological and ecological divergence

We mapped the tympanum condition on to the species tree to infer whether the loss of hearing structures took place once or multiple times. The hypothesized ancestral tympanic condition in Phrynopus along the internal branches of a phylogenetic tree including members of two other closely related genera (Lynchius and Oreobates) was determined using Bayesian stochastic character mapping (Huelsenbeck, Nielsen & Bollback, 2003) implemented in the phytools package (Revell, 2012). A mapping of tympanum presence/absence in the phylogeny was run on 500 stochastically mapped trees to take into account the uncertainty in the ancestral tympanum condition.

We also mapped the elevational range on to the species tree to visually assess the patterns of elevational distribution, phylogenetic relatedness, and presence/absence of tympanic membrane and annulus. For visualization purposes, we modified the range of species known from a single elevation or those known from a very narrow range (<50 m) to cover a minimum elevational band of 80 m. For species known from a single elevation, this elevation was assumed to be the midpoint and this number was used to expand the range on both directions (−40 and +40) to be 80 m. This allowed us to plot a color-coded bar representing the main habitat of each species. For species known to cover a narrow elevational band (<50 m), we first calculated the elevational midpoint (i.e., maximum elevation—minimum elevation) and used this value to expand the range on both directions (−40 and +40) to be 80 m. These adjusted ranges were used in only one figure (i.e., original data were used in all other figures and analyses).

Additionally, we examined if body size and body shape vary across habitats or elevations. Given that terrestrial breeding frogs living in high-elevation habitats appear to have larger body size and different body shape than species living at lower elevations (Hedges, 1999; Gonzalez-Voyer et al., 2011; Lehr & von May, 2017), we hypothesized that body size and body shape in Phrynopus are associated with the type of habitat use and elevation. First, we used phylogenetic principal components analysis (PCA) to examine morphological diversity between eared and earless species and among habitats. To remove the possible confounding effect of body size difference among species while controlling for non-independence resulting from shared phylogenetic history, we performed a phylogenetic size-correction using generalized least-squares (GLS) regression analysis (Revell, 2009). In this analysis, we considered SVL as the independent variable and all other morphological data were considered dependent variables. To perform the phylogenetic PCA, we pruned the species tree to only include 21 species of Phrynopus with corrected morphological data. We projected the first two PC axes on a phylomorphospace using the function phylomorphospace in phytools (Revell, 2012). We also used a phylogenetic ANOVA (Garland Jr et al., 1993) to quantify shape differences between eared and earless species and to quantify shape differences among habitats. For the latter comparison, we defined three groups (Andean grassland, montane forest, both habitats) and selected a posthoc test to compare means among groups; we used the “holm” method to adjust the P-values to account for multiple comparisons. These analyses were also implemented using phytools.

We used a scatterplot pairwise matrix to examine the relationships among body size, other body size-corrected morphometric variables, and elevation. We used the R package phylolm (Ho & Anné, 2014a; Ho & Anné, 2014b) to fit phylogenetic generalized linear regression models (PGLS), implemented under three models of evolution: (a) pure ordinary least squares (OLS) regression (i.e., no phylogenetic signal), (b) pure Brownian motion (BM) model, and (c) lambda transform model, which allows intermediate levels of phylogenetic signal between pure BM and pure OLS (we reported Pagel’s lambda coefficient, λ, for the third model). We used the AIC value to identify the model that best explains the variation in observed data (Ho & Anné, 2014b).

Micro-computed tomography

We visualized three-dimensional reconstructions of X-ray µCT images of the skull of the eared taxon Phrynopus peruanus, and confirmed the presence of the columella in this species (Fig. 2). In contrast, the columella has been completely lost in the earless P. montium (Fig. 2). This observation, in combination with previous findings (Trueb & Lehr, 2008), confirms that P. montium has lost the three main elements of the middle ear (tympanic membrane, tympanic annulus, and columella).

Molecular phylogenetic analysis

The best partitioning scheme included six subsets (BIC value: 34189.93). The first partition subset included both the 12S and 16S sequences and the best fitting substitution model was GTR + I + G. The remaining five subsets were partitioned according to codon positions as follows (substitution model in parenthesis): one set including the 1st codon position of COI (K80 + I), one set including the 2nd codon position of COI (HKY + I), one set including the 3rd codon position of COI (GTR + G), one set including the 1st and 2nd codon position of RAG1 and the 1st and 2nd codon position of Tyr (HKY + I), and one set including the 3rd codon position of both RAG1 and Tyr (HKY + G). At the end of the MCMC run, the average standard deviation of split frequencies was 0.002538. We found that Phrynopus is a well-supported monophyletic group that is closely related to Oreobates and Lynchius (Fig. 3). The earliest divergence within Phrynopus occurred between a clade formed by two species, P. peruanus and an undescribed species (Phrynopus spI), and all other species. Our analysis suggested that two nominal species, Phrynopus curator and Phrynopus nicoleae, are junior synonyms of Phrynopus tribulosus; a proposed formal taxonomic change is included below (New synonymy).

Phylogenetic relatedness and trait-and-habitat data

We used the species tree obtained with *BEAST to visually assess tympanum presence/absence (Fig. 4). Our analysis identified a single evolutionary transition that involved the loss of both the tympanic membrane and tympanic annulus (Fig. 4).

We used the species tree to assess the patterns of elevational distribution and phylogenetic relatedness (Fig. 5). We observed substantial differences in elevational distribution among closely-related species, with some species pairs occurring in distinct habitats (Fig. 5). In four cases, close relatives had non-overlapping elevational distributions (P. peruanus–Phrynopus spI; P. bracki–P. badius; P. miroslawae–P. tautzorum; P. unchog–P. daemon). Additionally, in three cases, close relatives were associated with either Andean grassland or montane forest (P. peruanus—Phrynopus spI; P. barthlenae–P. horstpauli; P. miroslawae–P. tautzorum).

Our analysis of body-size corrected morphological data showed that Phrynopus species exhibit a considerable variation in body shape. This variation appeared not to be associated with the presence or absence of tympanic middle ear structures (Figs. 6A and 6B); instead, species differences in body shape were associated with the habitat they occupy (Figs. 6C and 6D). Specifically, species that live in Andean grassland occupied a smaller area in the phylomorphospace projection that species that live in montane forest or those that occupy the two habitats.

Considering data from all adult specimens (males and females), the first seven principal components were found to explain 97.5% of the variance within the size-corrected morphological dataset following a phylogenetic size-correction using GLS regression (Table 1). Considering data from female specimens only, the first seven principal components were found to explain 96.7% of the variance within the size-corrected morphological dataset following a phylogenetic size-correction using GLS regression (Table 2). We found no shape differences between eared and earless species (phylogenetic ANOVA, all tests P > 0.10). We found that habitat had a significant effect on shape, notably in those variables (FL, HL, HW, IOD, E-N) that had greater contribution on PC2 (phylogenetic ANOVA, F = 3.943, P = 0.036).

We used a pairwise scatterplot matrix to visualize the cross-correlations among body-size corrected morphological data, SVL, and elevation (Fig. S1). Our PGLS analyses showed that there is a significant pattern of increasing body size at higher elevations, and that species living at higher elevations tend to develop shorter limbs, shorter head, and shorter snout than those species living at lower elevations (Table 3, Fig. 7). All tested relationships were significant both when controlling for phylogeny under the lambda transform model and also under pure OLS regression, and OLS was the best-supported model in five cases. One relationship, total length vs. elevational midpoint, was best supported under the Brownian motion (BM) model (and it was also significant under both the OLS and Lambda models).

New synonymy

Here we synonymize both Phrynopus curator Lehr, Moravec & Cusi, 2012 and Phrynopus nicoleae Chaparro, Padial & De la Riva, 2008 with Phrynopus tribulosus Duellman & Hedges, 2008.

Phrynopus tribulosus Duellman & Hedges, 2008

Phrynopus nicoleae, Chaparro, Padial & De la Riva, 2008; p. 53, Figs. 5 and 6 (MHNC 6441(holotype) collected by JC Chaparro, A Quiroz and D Salcedo at Santa Bárbara (10°20′36.3″S, 75°38′17.9″W; 3,589 m a.s.l.), Distrito de Huancabamba, Provincia de Oxapampa, Departmento Pasco, Peru, 2007)—Lehr & Oróz, 2012, p. 55, 59; Lehr, Moravec & Cusi, 2012, p. 51, 52, 55, 59, 61, 65, 66, Table 1; Padial, Grant & Frost, 2014, p. 124; Chávez et al., 2015, p. 21; De la Riva et al., 2017, p. 5, 31, 44; Lehr & Rodriguez, 2017, p. 333, 334; Rodríguez & Catenazzi, 2017, p. 383, 386, 390, 399, 401, 403, Table 2.

Phrynopus curator Lehr, Moravec & Cusi, 2012; p. 63, Figs. 5 and 7, 8, 10 (MUSM 31106 (holotype) collected by E Lehr, J Moravec, and JC Cusi at Quebrada Yanachaga (10°22.772′S, 75°27.717′W; 3,000 m a.s.l.), Yanachaga-Chemillén National Park (Sector San Daniel), Distrito de Huancabamba, Provincia de Oxapampa, Departamento de Pasco, Peru, 2012)—Padial, Grant & Frost, 2014, p. 124; Chávez et al., 2015, p. 18, 21, 22; De la Riva et al., 2017, p. 31; Lehr & Rodriguez, 2017, p. 333, 338; Rodríguez & Catenazzi, 2017, p. 383, 390, 399, 402, 403, Table 2.

Duellman & Hedges (2008) described Phrynopus tribulosus based on a single male collected in 1987, Chaparro, Padial & De la Riva (2008) described P. nicoleae based on a single female collected in 2007, and Lehr, Moravec & Cusi (2012) described P. curator based on a single female collected in 2012. Phrynopus tribulosus is known from three localities and an elevation from 2,600 to 3,589 m. Even though Chaparro, Padial & De la Riva (2008) and Lehr, Moravec & Cusi (2012) compared in their diagnoses the new species with P. tribulosus, sufficient morphological differences were noted to justify the descriptions. Some of the characters mentioned for P. tribulosus by Duellman & Hedges (2008) such as skin on venter smooth (areolate according to Chaparro, Padial & De la Riva, 2008 and Lehr, Moravec & Cusi, 2012) or dorsolateral folds absent (elongate dorsolateral warts forming a long discontinuous row that does not fuse to form a fold according to Chaparro, Padial & De la Riva, 2008, ridges forming discontinuous dorsolateral fold according to Lehr, Moravec & Cusi, 2012) are likely the result of long storage of the specimen in preservative. Whereas other differences such as absence of dentigerous processes of vomers as noted by Duellman & Hedges (2008) and Lehr, Moravec & Cusi (2012) but noted as present by Chaparro, Padial & De la Riva (2008) likely reflect intraspecific variation. Singleton species descriptions can be problematic when additional discoveries of specimens are difficult to assign because variation for the species is not covered or the species was not accurately described (see Lim, Balke & Meier, 2012; Köhler & Padial, 2016; Lehr & Rodriguez, 2017). Therefore, we have updated the species diagnosis of P. tribulosus based on Duellman & Hedges (2008), Chaparro, Padial & De la Riva (2008), and Lehr, Moravec & Cusi (2012) as follows.

Diagnosis

A species of Phrynopus having the following combination of characters: (1) Skin on dorsum shagreen with small scattered tubercles, middorsum with prominent ridges, often Y- or X-shaped, flanks shagreen with small scattered tubercles, skin on venter weakly areolate; ridges forming discontinuous dorsolateral fold; discoidal fold absent, weak thoracic fold present; (2) tympanic membrane and tympanic annulus absent; (3) snout rounded in dorsal and lateral views; (4) upper eyelid with or without enlarged tubercles; width of upper eyelid narrower than IOD; cranial crests absent; (5) dentigerous processes of vomers present or absent; (6) vocal slits and nuptial pads absent; (7) Finger I shorter than or equal to Finger II; tips of digits rounded; (8) fingers without lateral fringes; (9) ulnar tubercles absent; (10) heel with a distinct conical tubercle; inner tarsal fold absent; outer edge of tarsus with row of subconical tubercles; (11) outer metatarsal tubercle rounded, larger or equal in size of ovoid inner metatarsal tubercle; supernumerary plantar tubercles weakly defined; (12) toes without lateral fringes; basal webbing absent; Toe V slightly shorter or equal to Toe III; toe tips rounded, about as large as those on fingers; (13) in life, dorsum green with brown markings, tan with black irregular stripes and bluish-gray tubercles or reddish brown with dark gray and yellowish-brown mottling; venter gray, marmorated with small, brown, tan and metallic blue blotches or venter gray with pale gray mottling and brownish-cream flecks around posterior half of belly; groin with dark brown blotch, tan with abundant bluish-white spots and an orange spot or groin brown and gray mottled; (14) SVL in single male 15.2 mm, in females 20.7–21.2 (n = 2) mm.