AN ADAPTIVE RADIATION OF FROGS IN A SOUTHEAST ASIAN ISLAND ARCHIPELAGO

Living amphibians exhibit a diversity of ecologies, life histories, and species-rich lineages that offers opportunities for studies of adaptive radiation. We characterize a diverse clade of frogs (Kaloula, Microhylidae) in the Philippine island archipelago as an example of an adaptive radiation into three primary habitat specialists or ecotypes. We use a novel phylogenetic estimate for this clade to evaluate the tempo of lineage accumulation and morphological diversification. Because species-level phylogenetic estimates for Philippine Kaloula are lacking, we employ dense population sampling to determine the appropriate evolutionary lineages for diversification analyses. We explicitly take phylogenetic uncertainty into account when calculating diversification and disparification statistics and fitting models of diversification. Following dispersal to the Philippines from Southeast Asia, Kaloula radiated rapidly into several well-supported clades. Morphological variation within Kaloula is partly explained by ecotype and accumulated at high levels during this radiation, including within ecotypes. We pinpoint an axis of morphospace related directly to climbing and digging behaviors and find patterns of phenotypic evolution suggestive of ecological opportunity with partitioning into distinct habitat specialists. We conclude by discussing the components of phenotypic diversity that are likely important in amphibian adaptive radiations.

Our DNA dataset is based on sequences from 172 specimens of frogs of the family Microhylidae, including members of 17 outgroup genera. Because our goals were to estimate the phylogeny of the Southeast Asian genus Kaloula and previous phylogenetic analyses have not definitively established the sister taxon to this genus (e.g., Bossuyt et al. 2006;Frost et al. 2006;van Bocxlaer et al. 2006;van der Meijden et al. nearly all known allopatric populations of the widespread Kaloula baleata, including samples from Java (the type locality), Bali, Sulawesi, the Togian Islands, Peninsular Malaysia, Vietnam, Borneo, and Palawan Island (Philippines). Finally, we include samples of the widespread species Kaloula pulchra from throughout the entirety of its native range (China, Peninsular Malaysia, Laos, Thailand, Vietnam, Sumatra, and Sulawesi). However, our sampling also includes a number of undescribed species.
Further details on specimens and localities can be found in the Supplementary Table 1.

DNA DATA COLLECTION
We used either Qiagen DNeasy Tissue Kit (Cat. No. 69506) or a non-commercial guanidine thiocyanate method (Esselstyn et al. 2008) to extract genomic DNA from liver (or occasionally muscle) samples. Four primer pairs (MVZ59-Val, 12L1-16Sh, 12Sm-16Sa, and 16Sc-16Sd; Moriarty and Cannatella 2004;Darst and Cannatella 2004) were used to amplify via polymerase chain reaction (PCR) an approximately 2400 bp region of the mitochondrial genome spanning most of the 12S and 16S ribosomal RNA (rRNA) genes and the intervening transfer RNA (tRNA) for Valine. We used the now standard reaction mixtures and thermal cycle profiles for both PCR and sequencing reactions, using the same primers and following the methods of Moriarty and Cannatella (2004) and Evans et al. (2003). Samples were purified using Qiagen gel preps or Exosap purification protocols (USB Corp., Cleveland, OH, USA; Esselstyn et al. 2008). Sequencing reactions were conducted with identical undiluted PCR primers, using ABI Big Dye terminator chemistry (Perkin-Elmer, Boston, MA, USA) and Sephadex clean-ups (GE Healthcare, Uppsala, Sweden). Sequencing was performed on an ABI 3130xl automated PRISM sequencer (Applied Biosystems, Foster, CA, USA). All sequences are deposited in GenBank (Supplementary Table 1); a multiple alignment file is deposited in Dryad (doi:10.5061/dryad.dj342).
Because PCR products were amplified with varying success, the full length of the target sequence (approx. 2400 bp) could not be obtained for some specimens. In most cases this resulted from failure to amplify the 5' end of our target region, namely the MCZ59-Val region; the resulting sequences are approx. 1900 bp in length. Although missing data can bias or mislead phylogenetic inference, especially when concentrated in particular taxa or certain clades, taxa with missing data still add information to phylogenetic inference and missing data will not necessarily mislead phylogenetic analyses (Wiens 2006;Wiens and Moen 2008). Our goal in including specimens with incomplete data is to place particular samples with statistical support into the phylogeny to facilitate species delimitation and identification; in this way, we increased the number of specimens represented in both our molecular and morphological datasets. We included a given sample if we obtained sequence data for at least one of the four overlapping targeted gene regions; in most cases, these incomplete sequences derived from degraded tissue samples for which only 600-700 bp of the 16Sc-16Sd region could be amplified and sequenced. Taxa with incomplete gene representation were included in our analyses to aid in the confident identification of all specimens; analyses excluding taxa with missing data provided largely similar results (data not shown).
All fragments were sequenced in both directions. Sequences were assembled and manually vetted in Sequencher 4.5 (Genecodes, Ann Arbor, MI, USA), and then initially aligned using default parameter values in Clustal X 1.83.1 (Thompson et al. 1997). Non-overlapping DNA fragments from a single specimen were aligned separately and then merged in Mesquite 2.5 (Maddison and Maddison, 2008) to form a single terminal taxon for phylogenetic analyses. Alignments were examined exhaustively by eye and manually adjusted to minimize the number of parsimony informative change across sites (Moriarty and Cannatella 2003). Because excluding ambiguously aligned regions and autapomorphic insertion-deletions ("indels") from phylogenetic analyses resulted in minimal differences in topology and bootstrap values in exploratory parsimony searches (data not shown), we chose to maintain these sites in our analysis. The resulting alignment of 2567 bp corresponds to positions 2062-4620 of the Xenopus laevis mitochondrial genome (GenBank NC-001573). Uncorrected pairwise distances (pdistances) were calculated using MEGA 4.0.1 (Tamura et al. 2007).

PHYLOGENETIC ANALYSES
We used maximum-likelihood and Bayesian approaches to estimate phylogenetic relationships. A maximum-likelihood tree was obtained using RAxML v.7.0.4 (Stamatakis 2006) using a random starting tree, the faster rapid hill-climbing algorithm of (Stamatakis et al. 2007), and a GTR + I + Γ model of evolution (with four discrete categories for Γ); the data were not partitioned into subsets. Three hundred replicate searches were conducted to avoid basing our phylogenetic inference on a local search optimum; the tree with the best -ln likelihood from these replicate searches was selected as our maximum likelihood estimate of the phylogeny. Support for the preferred ML topology was estimated using 1000 independent nonparametric bootstrap replicates (ML BS), using the same model of sequence evolution, a random starting tree, and one search replicate per bootstrap replicate. Split-support for clades was summarized using SumTrees in DendroPy (Sukumaran and Holder 2010) and clades present in ≥ 70% of the bootstrap trees were considered well-supported (Hillis and Bull 1993).
A Bayesian estimate of phylogeny was obtained with MrBayes v.3.1.2 (Ronquist and Huelsenbeck 2003) using the GTR + I + Γ model and unpartitioned dataset. Four replicates of four MCMC chains were run for 20 million generations, sampled every 2000 generations, using a temperature of 0.2 and default priors. We examined trends and distributions of log-likelihoods and parameter values using Tracer 1.5 (Rambaut and Drummond 2009). We assessed convergence by examining correlations of split frequencies and chain variability among independent runs using Are We There Yet? (AWTY; Nylander et al. 2008). Based on examination of trends in Tracer and correlations of split frequencies, stationarity and convergence was achieved after one million generations. However, because of patterns of chain variability during the first few million generations, we took a conservative approach and discarded trees sampled during the first ten million generations as burn-in. Based on the cumulative set of post-burn-in trees from each run, all effective samples sizes (ESS) were greater than 400. The topology and posterior probabilities (PP) were then summarized based on the post-burnin trees using SumTrees and Tracer. We considered topologies with posterior probabilities ≥ 0.95 to be well supported (Huelsenbeck et al. 2001).
We generated a time-calibrated estimate of phylogenetic relationships to test patterns of lineage diversification and morphological evolution. Analyses were conducted using the uncorrelated relaxed clock method (Drummond et al. 2006) implemented in BEAST v.1.6.2 (Drummond et al. 2010). Phylogenetic relationships were estimated using the Yule pure-birth speciation prior for the tree shape and a GTR + I + Γ model of sequence evolution. For comparative analyses, we generated chronograms with relative divergence times by setting a normal prior of 100 (mean: 100; standard deviation: 0.001) on the most recent common ancestor (MRCA) of all tips in the tree; the distant outgroup Callulina was excluded from these analyses. MCMC analyses were run for 200 million generations with MCMC steps and divergence times recorded every 5000 generations; to avoid autocorrelation between MCMC steps, the sampling was further thinned to every 10,000 generations. We then assessed stationarity of parameter estimates using Tracer and thus discarded the first 5 million generations as burn-in; all ESS values were > 190.

ECOTYPES AND MORPHOLOGICAL VARIATION
Our characterization of ecotype summarizes data on ecology and natural history, including reproductive and larval biology, and microhabitat preferences. These assessments utilize a half-century of accumulated information on the life history and biology of the frogs of interest (Inger 1954(Inger , 1966Inger andStuebing 1989, 1997;Zhao and Adler 1993;Dutta and Manamendra-Arachichi 1996;Manthey and Grossman 1997;Fei 1999;Inger et al. 1999;Fei and Ye 2000;Diesmos et al. 2002;Schleich and Kastle 2002;Malkmus et al. 2002), as well as our collective personal experience with many of these species.
We selected measurements of morphological features likely related to habitat utilization in anurans. For instance, we expect the relative size of limbs, toe pads, and metatarsal tubercles to be directly related to behaviors of habitat specialists such as climbing or burrowing. These measurements are standard in Asian anuran taxonomy (Inger 1954(Inger , 1966Taylor 1962;Matsui 1984Matsui , 1994Inger andStuebing 1989, 1997;Zhao and Adler 1993;Dutta et al. 1996;Fei and Ye 2000;Brown and Guttman 2002;Malkmus et al. 2002) with the exception of our measure of webbing (see below). The measurements used in our analyses are as follows: snout-vent length; head width (measured at the widest point of the head, near the jaw articulation); snout length (measured from the most anterior margin of the eye to the tip of the rostrum); forearm length (measured from the most proximal margin of the elbow to the proximal margin of the most proximal palmar tubercle); third finger length (measured from the proximal margin of the most proximal subarticular tubercle to the distal digit tip); third finger width (measured at the level of the most distal subarticular tubercle); third finger-tip width (maximal width of the distalmost finger); thigh length (measured from the cloaca to the distalmost knee); crus length (measured from the most proximal knee to the most distal tibiotarsal joint); inner metatarsal tubercle length (maximum length); outer metatarsal tubercle length (maximum length); third toe length (measured from the proximal margin of the most proximal subarticular tubercle to the distal digit tip); third toe width (measured at the level of the most distal subarticular tubercle); third toe-tip width (maximal width of the distalmost toe); and extent of pedal webbing (measured from the distal tip of the third toe to the most proximal lateral margin of the toe where intradigital webbing is attached). Typically, intradigital webbing is described using the formula developed by Heyer (1967, 1997). However, this description does not lend itself to analyses using multivariate statistics. We instead opted for a novel approach by taking a measurement in which greater values indicate less webbing; thus, for example, if a principal components axis loads strongly and positively on this variable it means that species with large positive scores have less webbing than those with negative scores. For some taxa (e.g., K. kalingensis), the outer metatarsal tubercle is essentially absent and thus results in a measurement of nil. Because all mean scores were subsequently natural log-transformed, we added the value of 1 to all mean scores before further calculations. We were thus able to include information from all species for all measurements even though the outer metatarsal tubercle is absent in some species. All measurement data used for these analyses are deposited in Dryad (doi:10.5061/dryad.dj342).

RECONSTRUCTION OF ANCESTRAL ECOTYPES
We evaluated patterns of change in ecotype class across the phylogeny and determine the ecotype class of the MRCA of the Philippines clade. Because of uncertainties in the phylogenetic estimate, we used Bayesian methods for reconstructing ancestral ecotype states across Kaloula phylogeny. In diversification analyses, we used the pruned specieslevel time-calibrated post-burn-in trees from our BEAST analysis to reconstruct ancestral states and analyze patterns of morphological evolution. We integrated over uncertainties in topologies using Bayesian mutational mapping as implemented in SIMMAP 1.5.2 (Bollback 2006), which employs a stochastic model of character change. SIMMAP analyses treated character states as unordered and utilized an equal prior for the bias parameter, a gamma distribution for the rate parameter prior (α = 1.25, β = 0.25, k = 60), 10 sample replicates, and 10,000 draws from the prior distribution in order to have sufficient sampling of the posterior distribution for bias and rates. We evaluated the robustness of our inference to priors for the gamma distribution by conducting similar analyses in which we varied α, k, and θ; because the results of these various analyses were qualitatively similar (data not shown), we interpret our inference as robust.

PHYLOGENETIC RELATIONSHIPS
Of the 2567 sites included in the phylogenetic analyses, 1469 are variable and 1145 are parsimony informative. Below, we discuss patterns of relationships based on maximum likelihood (ML) and Bayesian analyses (including the maximum clade credibility tree There is strong support (PP = 1.00; ML NBS = 96%) for a clade containing divergent and geographically circumscribed populations typically referred to K. baleata, but which we consider a complex of species, some of which are obviously morphologically distinct; we refer to this as the "baleata clade." A new morphologically distinct species within the baleata clade that occurs on Palawan is the only Philippine species of Kaloula not derived from the primary Philippine radiation and, thus, represents a second invasion of the archipelago. Relationships of the baleata clade remain unclear.
In the ML tree, the baleata clade is sister to the endemic Philippines clade (PP = 0.57; ML NBS = 44%) whereas in the MCCT tree, it is sister to the mainland species K. mediolineata (PP = 0.16; ML NBS = 24%). The clade containing K. mediolineata, the baleata clade, and the endemic Philippine clade receives high support (PP = 1.00; ML NBS = 96%). There is also strong support for the relationship of K. verrucosa (China) and K. pulchra (a widespread Asian species) forming successively branching taxa sister to the remaining members of Kaloula (Fig. 1).
Most species, subspecies, and candidate species are recovered as monophyletic with high support ( Table 2). The relationships of one specimen of K. rigida from the type locality remains unclear, but otherwise K. ridiga and K. walteri are each recovered as monophyletic with strong support. A number of well-supported clades exist within the Philippines, though relationships between these clades remain unclear. Kaloula rigida and K. walteri are resolved as sister taxa (PP = 1.00; ML NBS = 100%). There is strong support for the monophyly of an undescribed species from Samar and Leyte islands.

INTRASPECIFIC PATTERNS OF GEOGRAPHIC DIVERSIFICATION
Our intraspecific sampling within several species allows for a characterization of divergence within these taxa (Table 2). For reference to the values reported for other species, the mean p-distance across Kaloula (excluding K. taprobanica) is 7.5%. Two widespread species, K. pulchra (throughout Southeast Asia) and K. picta (throughout the Philippines) exhibit nearly identical sequences across their respective broad geographical distributions (K. pulchra: 0.5%; K. picta: 0.2%). Kaloula baleata consists of a highly divergent and morphologically distinct Vietnam population (> 7% divergent from other species in the "baleata" clade) that is sister to a clade comprising moderately divergent lineages from Borneo, the Malay Peninsula, Borneo, Java + Bali, Sulawesi, and Palawan (< 3% from one another).
Divergences among species in the K. kalingensis clade are an order of magnitude greater than those observed for other species of Kaloula ( Specimens identified as K. c. stickeli were genetically identical to those identified as K. c.

meridionalis.
Moderate levels of both geographic structure and genetic divergence were also detected within allopatric populations of the species K. rigida and K. walteri.

ECOTYPE EVOLUTION
Analysis using Bayesian mutational mapping provides an ambiguous perspective on patterns of ecotype evolution. For example, posterior probabilities of each ecotype for the MRCA of the endemic Philippines radiation are all very similar (tree hole frog: 0.33; ground frog: 0.36; scansorial shrub frog: 0.31). Two contrasting patterns of genetic diversity were unexpected. Two species demonstrated almost no geographic-based genetic variation, despite extensive morphological variation and wide geographical ranges. Both K. pulchra (throughout Southeast Asia) and K. picta (throughout the Philippines) exhibit nearly identical sequences across their respective broad distributions and yet both have been the subject of much discussion by taxonomists who anticipated impending taxonomic subdivision of these species (Parker 1934;Inger 1954Inger , 1966. Populations historically referred to Kaloula baleata (type locality: Java) include several putative new species, including a highly divergent and morphologically distinct Vietnam population (> 7% divergent from other species in the baleata clade) and morphologically similar populations from Peninsular Malaysia, Java + Bali, and Sulawesi (< 3% from one another). In addition, the Philippine exemplar of the baleata clade (Kaloula sp. nov. from Palawan) is a highly distinctive, small-sized, tree-hole breeding species. In accordance with their status as distinct and diagnosable allopatric evolutionary lineages (Wiley 1978;de Queiroz 1998de Queiroz ,1999, we consider these allopatric populations of the baleata clade to be distinct species that require formal description.

Our analysis suggests that
Divergences among species in the kalingensis clade are much greater than those observed between other species of Kaloula (Table 2). The kalingensis clade contains two described species (K. kalingensis from northern Luzon, including the type locality in Kalinga Province, and K. kokacii from the Bicol Peninsula of southern Luzon Island, adjacent to the type locality on Catañduanes Island), and two undescribed species from Panay Island and east Luzon Island, respectively.
Kaloula conjuncta exhibits moderate levels of divergence between clades, not all of which correspond to existing taxonomy. Kaloula c. conjuncta from numerous localities on Luzon Island was paraphyletic with respect to a population on Mindoro Island that is morphologically and acoustically distinct (Brown et al. unpubl. data).
Similarly, populations that are morphologically and acoustically diagnosable as K.
conjuncta negrosensis from Panay and Negros Islands did not form a clade exclusive of other subspecies; in the ML estimate, the Panay population of K. c. negrosensis is most closely related to a putative new species from Sibuyan Island, the latter of which is morphologically and acoustically distinct from all other taxa (Brown et al. unpubl. data).
Mindanao specimens phenotypically similar to the type series of K. c. stickeli (holotype FMNH 60786) are nested within K. c. meridionalis and are genetically identical to this morphologically distinct taxon. Kaloula c. stickeli is morphologically intermediate between K. c. meridionalis (a smaller scansorial species with widely expanded digital disks) and K. picta (larger terrestrial species with non-expanded digital disks) and was originally described from northern portions of the Mindanao aggregate island complex (Samar and Leyte Islands; Inger 1954). The possibility of a hybrid origin for K. c. stickeli seems likely. Specimens included in our analysis were collected on Mindanao Island from a large mixed chorus of individuals of K. c. meridionalis, K. picta and those similar to K. c. stickeli. Interspecific amplexus was observed between K. c. meridionalis and K. picta (RMB and JAM, personal observation). Taken together, the field observations and genetic similarity between specimens identified phenotypically as K. c. stickeli and K. c. meridionalis suggest that K. c. stickeli is a morphologically intermediate phenotype that results from natural hybridization between these K. c. meridionalis and K. picta. As a result of these findings, and for the purposes of our comparative analyses, K. c. stickeli was not treated here as a valid taxon. We anticipate that future taxonomic work will place K. c. stickeli in the synonymy of K. c. meridionalis.
In fact, this hybridization among species in the Philippines radiation is not surprising and fits an emerging pattern highlighting hybridization as a common and important part of radiations (Grant and Grant 2002;Seehausen 2004;Wiens et al. 2006;Cristecu et al. 2010).
Finally, moderate levels of both geographic structure and genetic divergence were also detected within allopatric populations of the species K. rigida and K. walteri. Our one sample of K. rigida from the type locality (Baguio City, Luzon Island) is highly divergent from the remaining K. rigida from the forested mountains of northern Luzon and falls out sister to K. walteri plus the remaining K. rigida. Future studies, involving targeted geographical sampling and additional sampling near Baguio City will be necessary to resolve uncertainty in the status of K. rigida.
In summary, several discrepancies between existing taxonomy and observed distribution of phylogenetic diversity will require an eventual comprehensive taxonomic review of species diversity-an undertaking beyond the scope of this paper. We note that the recognition of the above undescribed species will result in ~ 40% increase in known diversity in Kaloula (see Brown and Diesmos 2002, Brown 2007, and Brown et al. 2008 for review of rapidly changing taxonomy in Philippine anurans). This increase in diversity is unsurprising given recent similar studies of taxonomic diversity and cryptic species in Southeast Asian amphibians (Brown et al. 2008;Stuart and Bain 2008;Brown and Stuart 2012).