- Research article
- Open Access
Divergence at the edges: peripatric isolation in the montane spiny throated reed frog complex
BMC Evolutionary Biology volume 15, Article number: 128 (2015)
Peripatric speciation and peripheral isolation have uncertain importance in species accumulation, and are largely overshadowed by assumed dominance of allopatric modes of speciation. Understanding the role of different speciation mechanisms within biodiversity hotspots is central to understanding the generation of biological diversity. Here, we use a phylogeographic analysis of the spiny-throated reed frogs and examine sister pairings with unbalanced current distributional ranges for characteristics of peripatric speciation. We further investigate whether forest/grassland mosaic adapted species are more likely created through peripatric speciation due to instability of this habitat type.
We reconstructed a multi-locus molecular phylogeny of spiny-throated reed frogs which we then combined with comparative morphologic data to delimit species and analyze historical demographic change; identifying three new species. Three potential peripatric speciation events were identified along with one case of allopatric speciation. Peripatric speciation is supported through uneven potential and realized distributions and uneven population size estimates based on field collections. An associated climate shift was observed in most potentially peripatric splits. Morphological variation was highest in sexually dimorphic traits such as body size and gular shape, but this variation was not limited to peripatric species pairs as hypothesized. The potentially allopatric species pair showed no niche shifts and equivalent effective population sizes, ruling out peripatry in that speciation event. Two major ecological niche shifts were recovered within this radiation, possibly as adaptations to occupy areas of grassland that became more prevalent in the last 5 million years. Restricted and fluctuating grassland mosaics within forests might promote peripatric speciation in the Eastern Arc Biodiversity Hotspot (EABH).
In our case study, peripatric speciation appears to be an important driver of diversity within the EABH biodiversity hotspot, implying it could be a significant speciation mechanism in highly fragmented ecosystems. Extensive peripatric speciation in this montane archipelago may explain the abundance of discrete lineages within the limited area of the EABH, as inferred in remote island archipelagos. Future phylogenetic studies incorporating demographic and spatial analyses will clarify the role of peripatric speciation in creating biodiversity hotspots.
Isolation has been historically considered the most frequent mechanism to limit gene flow, promoting divergence of lineages leading to speciation [1, 2]. The two dominant proposed forms of isolation have been vicariant allopatric (dichopatric) and peripatric, while sympatric and parapatric speciation appear less common  (see below). The classical model of vicariant allopatry starts with an ancestral distribution, which is split into large, disjointed populations ; hereafter referred to as “vicariant allopatry” or simply “allopatry”. Numerous examples have been cited demonstrating the process of allopatry in driving speciation (summarized in [3, 4]), and the signal of vicariant barriers limiting gene flow has been shown in many comparative phylogeographic studies [5–8]. In peripatric speciation, a small population is isolated at the edge of a larger population and rapidly diversifies, often crossing a previous barrier to dispersal and/or occupying a newly available ecological niche/habitat type. Over short periods of evolutionary time, peripatric events leave different evolutionary footprints from allopatry in terms of range, historical demography, and degree of ecological and morphological divergence . Parapatric speciation events also involve a shift in habitat requirements or ecological niche, though unlike in peripatric speciation, there are no physical barriers separating divergent lineages. Sympatric speciation, where species diverge within the same environment without spatial or ecological barriers, has rarely if ever been observed . New analytical approaches using multiple types of data enable testing for these processes in greater detail [10–13].
Different speciation processes produce predictable outcomes in the molecular, morphological, and geographic relationships among lineages. In contrast to the stable population sizes associated with allopatric speciation, a peripatric speciation event yields genetic variation consistent with bottlenecks and low founder population numbers . Peripatric species are also more likely to be ecologically divergent, as they have crossed a barrier to dispersal that previously limited occupation of new habitat at the edge of the ancestral range [3, 15]. Allopatric species, in contrast, more often retain their ecological niche due to the similarity of habitats isolated by vicariant barriers [16–18]. Morphologically, divergence in physical traits, especially those shaped by sexual selection, are predicted to be especially affected in peripatric populations both due to drift in small populations and as reinforcement between sister lineages in close proximity [19–22]. Quantifying molecular, ecological, and morphological metrics among sister lineages can therefore distinguish between these two modes (allopatry and peripatry) in species formation.
Across the African continent, speciation patterns are largely consistent with fragmentation of formally continuous forests through fluctuating climate [23–28]. Such patterns occur in the Eastern Afromontane Biodiversity Hotspot (EABH). Here, isolation of forest endemics on refugial mountain tops appears to be a common driver of exceptional levels of endemism and species diversity (e.g., [26, 29–31]). However, in part due to incomplete phylogenetic and phylogeographic sampling and imprecise distribution models, little attention has been paid to whether regional-scale fragmentation (vicariance/allopatric) or local-scale expansions or niche shifts (peripatric) are the dominant speciation process within this sky-island system.
We evaluate an insular montane radiation of spiny-throated reed frogs within the EABH and assess regional-scale and local-scale speciation mechanisms. Deciphering these two processes for individual radiations is imperative to understand diversification in tropical ecosystems, as both niche conservatism and niche divergence appear to play a major role in amphibian speciation [32, 33]. Establishing the speciation processes involved in the spiny-throated reed frog radiation allows for direct comparisons with other EABH endemic radiations as well as linking these processes to montane and island archipelagos around the world.
Within the EABH, the arrangement of montane areas of varying size and distance from each other provides an ideal setting in which to compare mechanisms of divergence. The mountains of the EABH system are dominated today by two forest vegetation types: rainforest found approximately between 700-1800 m or lower depending upon the mountain block and aspect (referred to as “rainforest” hereafter), and a montane forest/grassland mosaic (referred to as “mosaic” hereafter) above ca. 1800 m in the southern EABH . The extent and distribution of these biomes throughout the EABH has varied through time, with a gradual shift towards increased grasslands from 4 mya onward . The areas of high montane forest/grassland mosaic have historically expanded and contracted and at times likely vanished or become restricted to extremely narrow refugia as forests expanded . Lower elevation rainforest habitats fluctuated to a lesser extent, and were likely persistently stable on montane slopes despite habitat cycling in the adjacent lowland savannahs (classically outlined by ; see [38–41]). The larger fluctuations, smaller potential distributions, and the relative newness of mosaic habitat compared to mid-elevation rainforests, may make the mosaic habitat biome more susceptible to peripatry mechanisms than the more persistent rainforest biome. Though other biome-specific patterns of diversification are known [42, 43], whether peripatric speciation could also be biome-specific is less well documented .
In this study we evaluate the African spiny-throated reed frog complex from the EABH of Tanzania, Malawi, and Mozambique (Fig. 1) and quantitatively assess the isolation mechanisms among sister-lineages. The spiny-throated complex is comprised of six taxa (H. burgessi, H. davenporti, H. spinigularis, H. tanneri, H. minutissimus, and H. ukwiva) , with species occupying localities often restricted to either montane mosaic or rainforest habitats [26, 46–49]. By using molecular, morphological and ecological indicators, we make a first assessment of whether allopatric or peripatric speciation processes are supported and whether peripatric speciation is facilitated by the fluctuating montane forest/grassland mosaic biome. These results highlight a possible speciation syndrome within biodiversity hotspots composed of a variety of habitat types where variation in persistence may generate endemicity.
Summary statistics on the loci used in this study are shown in Table 1. All molecular sequences are deposited along with voucher numbers in GenBank and are available along with all specimen information in Additional file 1: Table S1.
Species tree methods all agreed on a single reconstruction of the evolutionary relationships within this group, though support for the sister relationship of H. ukwiva and H. minutissimus was weak or unresolved in some methods (Fig. 2). All methods show monophyletic relationships supporting six lineages: the previously described species of H. tanneri and H. minutissimus, a new species from the Rubeho Mountains (H. ukwiva), three distinct lineages (H. spinigularis, H. davenporti, and H. burgessi) which once were lumped into a single species (H. spinigularis) (see ). Hereafter the clade previously refered to as H. spinigularis will be refered to in quotes (“spinigularis” clade) and the newly defined species names will be used. RAxML genetrees are shown in Additional file 2: Figure S1. Specimens previously identified as H. spinigularis from the Udzungwa Mountains  were included in this analysis and all were definitively genetically and morphologically H. minutissimus, confirming earlier estimates that only H. minutissimus is found in the Udzungwa Mountain block.
Species delimitation methods (Bayesian General Mixed Yule Coalescent - bGMYC, Bayes factor delimitation - BFD, Bayesian Phylogenetics and Phylogeography - BPP) were largely in agreement in defining lineages within this radiation except for some variation in support of additional distinct lineages in addition to the six described here and in  (Fig. 2; Table 2 and Additional file 2: Table S2). In each case of additionally distinct lineages associated with the named six (Udzungwa Scarp locality of H. minutissimus shown in dark green for bGMYC, Mozambique (Mt. Namuli) specimen for H. spinigularis shown in yellow for BPP, Uluguru Mountains locality of H. burgessi shown in pink for bGMYC), a distinct locality with comparatively low sampling within a large distribution is involved. Additional sampling will be necessary to determine whether additional taxonomic units are supported. The bGMYC probability threshold above which individuals are considered conspecific was set at 0.5 (equivalent to selecting the posterior mean of the analysis) to avoid over lumping or splitting species (as outlined in the program manual). BFD support for the six-species model shown in Fig. 2 was positive (Path sampling: BayesFactor of 5.4, Stepping Stone BF of 5.8 over the next best model), as was BPP3 (Posterior probability = 0.83410).
Effective population size, population bottlenecks and expansions
No populations showed signs of ancestral bottlenecks or expansions, though the narrowly distributed species (H. tanneri, H. ukwiva, H. davenporti; distributions in Table 3) did not have sufficient sample sizes to analyze Extended Bayesian Skyline Plots (EBSP). In EBSP plots, most populations showed potential current bottlenecks, which may be related to recent range reductions of habitats due to climatic and/or anthropogenic events (Fig. 3). Only the population of H. burgessi in the East Usambara Mountains does not appear to show a current reduction, though this is also the population with greatest sampling. EBSP plots of all individuals from deeper nodes did not show bottlenecks or expansions except when including all individuals from all species, which shows an expansion beginning around 8 MYA that peaks around 4 MYA. This coincides with the split between H. tanneri and the H. spinigularis/burgessi/davenporti clade.
BPP estimates of Theta (θ) show essentially stable mean and 95 % confidence intervals for all splits within the tree (e.g., θ from median gamma prior across nodes: ~ 0.003 (~0.0007, 0.007)). However, estimates from the limited-distribution species are based on 1–3 individuals, and thus are not informative. Stable θ estimates of the species with larger distributions (and corresponding divergence nodes) agree with EBSP population estimates that these populations have been stable through time without marked reductions or expansions (Additional file 2: Figure S2).
Measurements for all traits and species summary statistics are presented in . In males, aspects of gular shape (height, width) contributed most to distinguishing species in PC1 (48 %, p = 0.001) (Fig. 4, Additional file 2: Table S3), though overall size (snout-urostyle length, SUL) was also significantly involved. PC2, though contributing much less to distance between the groups, separates H. tanneri and H. spinigularis from all other species based on head length and gular flap height. Both gular shape estimates and SUL were significant in single trait ANOVA estimates (Table 4). For females, head-width, SUL, and leg measurements were the primary distinguishing characteristics in PC1 (62 %, p = 0.01). In ANOVA analyses of these traits controlled for size, only SUL was significantly different between species (Table 4, ). Presence and distribution of gular and ventral spines are different between all species  and is a diagnostic tool for separating species (see ).
Environmental niche divergence, timing of shifts, and range sizes
PCA analysis of full and reduced bioclim variables (19 vs. 7) showed equivalent patterns, and thus only the reduced dataset is discussed. Two major clades of species were identified along PC1 marking the split between a rainforest adapted group (H. burgessi, H. tanneri, and H. spinigularis) and a mosaic-adapted group (H. ukwiva, H. minutissimus, and H. davenporti) (non-overlapping 95 % confidence intervals; Fig. 5, Additional file 2: Table S4). PC1 (53 % explained variance) is dominated by the temperature and precipitation during the coldest and driest times of year. These habitat syndromes, when mapped onto the phylogenetic tree, show two major adaptive habitat shifts to high elevation grasslands (Fig. 5), presuming that the ancestor was forest-adapted in line with paleoclimatic estimates of conditions ~ 4 MYA , before a general shift to grasslands at 1.86 MYA .
In testing potentially peripatric sister species pairs for niche shifts, H. davenporti falls outside of the PCA confidence ellipse probability threshold for H. burgessi into different niche space (Fig. 5). The presence of H. davenporti and H. burgessi in two divergent areas of niche space (mosaic vs. rainforest) implies a major niche shift at the time of speciation. Hyperolius ukwiva localities were outside of the PCA confidence ellipse probability threshold for H. minutissimus (Fig. 5), but only separated by PC2. The appropriate comparison for H. tanneri in a phylogenetic context is between H. tanneri and the ancestor of H. burgessi, H. spinigularis, and H. davenporti. To approximate their ancestral niche, we combined distributional data from these three species into a lumped distribution, which we then compared to H. tanneri in a PCA analysis as above. No separation was seen between H. tanneri and the combined niche space of the “spinigularis” group (H. burgessi, H. spinigularis, and H. davenporti) (Additional file 2: Figure S3).
The two cases of adaptation to grasslands (basal split and H. burgessi/H. davenporti) occurred at differing times (Fig. 5, ~5 MYA and ~ 1 MYA) and therefore a single event cannot have caused these shifts. Both events coincide with periods of increased climatic cycling and a gradual increase in coverage from the forest/grassland mosaic from the previously extensively distributed rainforests of this region. Geographic range estimates (Table 3), underline range asymmetries between sister species .
New lineages within the spiny-throated reed frog clade
Prior to this study and our complementary taxonomic revision , only three species were recognized within the spiny-throated reed frog radiation within the highlands of East Africa: H. spinigularis, H. tanneri, and H. minutissimus. Our field surveys in isolated mountain ranges and associated molecular analyses and morphological assessments of members within this group have documented three new lineages within this radiation: H. burgessi, H. davenporti, and H. ukwiva, each corresponding to a different mountain block or blocks within the EABH. In addition to morphological diagnoses, species delimitation methods based on molecular data were in general agreement in the recognition of each of these species. Deeper population-level sampling should help clarify the levels of divergence within and between species. Nodal support was strong in most methods, though *BEAST showed weak support at basal nodes. Increased sampling of additional molecular markers may help clarify relationships in these deeply divergent lineages.
Peripatric speciation or not
Three species within this radiation are known from single localities in the mountains of the EABH (H. tanneri, H. davenporti, H. ukwiva), raising the question of whether these lineages are the result of peripatric speciation events when isolated from their broadly distributed sister species. Though evidence is limited, many cases of speciation within the spiny-throated reed frogs are consistent with peripatric speciation. In support of peripatric origins, the three species with limited distributions all have smaller realized and potential distributions than their broadly distributed sister lineages (Table 3). Uneven population sizes and population bottlenecks could not be explicitly tested in peripatric species as field observations showed very low population densities and genetic material was limited, yet the rarity of populations and low population sizes of the limited-distribution species were apparent from repeated surveys across localities . In contrast, their sister lineages showed much larger distributions with large and stable population sizes through time. Minimally, the lack of bottleneck patterns from the well sampled, broadly distributed species and limited distribution of sister lineages fits with our a priori predictions for peripatric speciation.
We expected peripatric lineages would differ in climate/habitat, and this was found in 2/3 pairs. Though a habitat shift is not necessary for peripatric speciation, these shifts are less likely in vicariant speciation. The expectation of accentuated divergence in sexually-dimorphic traits in peripatric lineages compared to allopatric was not altogether supported, as all species differed in SUL, gular shape, and the arrangement of gular spines which we assume to be involved in some aspect of sexual selection. There is, however, a complete loss of gular spines in H. tanneri differing substantially from its sister group. This rapid morphological change in H. tanneri, potentially the result of sexual selection, is what would be predicted for peripatric speciation.
Within the allopatric populations of both H. burgessi and H. spinigularis, we find conserved ecological niche/environmental tolerance despite population distances greater than most peripatric species pairs, and no evidence of population bottlenecks. As most species within this radiation are known from single mountain blocks, we were concerned that morphological variation would be high even within species on different mountain blocks (H. spinigularis and H. burgessi have distributions spanning multiple mountain blocks). However, this was not the case in H. burgessi, where morphology was conserved despite substantial divergence times (ca. 0.3-1 mya) between different blocks. Sufficient morphological data are not yet available for the newly discovered Mozambique population of H. spinigularis to determine if they are similarly conserved in relation to Malawian specimens.
Distinguishing between allopatric and peripatric speciation can be difficult, as these processes anchor a continuum of progressively asymmetric isolation. The precise point at which a process ceases to be allopatric and becomes peripatric is therefore blurred. The cases of mixed support in the indicators we examined might therefore be a consequence of the extent to which some of the events were truly “peripatric” or “allopatric” events. Likewise, the indicators themselves (current range sizes, potential ecological niche shifts, and accelerated morphological diversification) may yield ambiguous results when evaluating species pairs that are deeply divergent. With the current data available on the spiny-throated reed frogs, it seems possible that peripatric speciation plays a role in generating biodiversity within the forests and grasslands of the EABH. In systems such as the EABH, which contain a high concentration of old endemics with limited and fragmented distributions in isolated highlands , a community approach may be necessary to clarify the prevalence of these speciation mechanisms in order to more fully understand the formation of its rich biota.
Biome specific speciation mechanisms
Based on the comparisons outlined here from spiny-throated hyperoliids found across these habitats – it seems that peripatric processes may be more pronounced in higher elevation species compared to those in more persistent rainforests (e.g., ). If there were a correlation between the less persistent montane forest/grassland mosaic sites and peripatric speciation, we would expect a common biogeographic pattern repeated across other species. Little data from other amphibians of peripatry has been reported worldwide (e.g., [41, 54, 55]), however, and information is particularly lacking for amphibians within the EABH. Similarly little is known about the phylogenetic relationships and demographic history among EABH high-altitude species and whether patterns exist that might be consistent with peripatry. Evidence in analogous systems, however, suggests that such a distinct boundary in biomes – as in rainforest and grasslands - might be important for generating specific speciation mechanisms. For example, ecotones across Savanna-Forest mosaics of Central Africa been suggested to be important in generating rainforest diversity through peri – and para – patric processes ; however, the generality of such patterns has not be quantified. The role of distinct habitat boundaries promoting speciation processes clearly deserves further investigation, potentially such situations might favor fundamentally different speciation processes compared to other areas (such as vicariance within lower elevation forest habitats). Likewise, there are also clear conservation implications that arise from these considerations. If high elevation species consistently exhibit limited population numbers and historical population genetic bottlenecks, they may be at increased risk if exposed to global anthropogenic pressures such as habitat destruction, fragmentation, and climate change.
Molecular, morphological, range, and niche analyses provide mixed support for cases of peripatric speciation in spiny-throated reed frogs. The highest support is in species pairs exhibiting extreme potential and realized distributional asymmetry, niche shifts, and morphological divergence driven by sexually selected traits: H. burgessi vs. H. davenporti; H. minutissimus vs. H. ukwiva. The older break between H. tanneri and the “spinigularis” clade of H. spinigularis/burgessi/davenporti is less clear, though the limited distribution and unique loss of gular spines (in this clade) in H. tanneri both support the potential for a peripatric split across the 10 km gap separating their current distributions. Though low survey numbers and extremely limited range size in potentially peripatric populations limited our ability to employ demographic reconstructions to confirm peripatric speciation from a signature of founding bottlenecks within the narrowly-distributed species, these conditions themselves support that these are small and divergent populations on the periphery of sister lineages with larger geographic distributions. The emergence of the highland forest/grassland mosaic ecosystems during the same historical interval of this radiation may have facilitated a peripatric mode of speciation as these novel habitats created dispersal opportunities within the larger ranges of rainforest adapted species. Whether the high elevation species that occur along the edge of forest and grassland areas are more amenable to peripatric speciation than rainforest-restricted species awaits additional assessment across more taxonomic groups. However, biome specific peripatric speciation process appears to be supported in some EABH hyperoliids.
We sampled 40 individuals (Hyperolius spinigularis = 8, H. tanneri = 3, H. burgessi = 21, H. ukwiva = 1, H. minutissimus = 6, H. davenporti = 1) from nine mountain blocks, covering the entire range of the spiny-throated reed frog complex (Fig. 1, Sample sizes, Museum numbers and GPS coordinates in Additional file 1: Table S1). Outgroup members of the genus Hyperolius and the family Hyperoliidae were also collected from the field. We collected full voucher specimens, tissue samples, photographs, and microhabitat information. All specimens for this study were collected in accordance with animal ethics guidelines established in the institutions of authors (including Field Museum of Natural History, Science Museo of Trento, and University of Basel).
We extracted total DNA from liver and leg muscle tissue of freshly collected specimens preserved in 95 % ethanol or 20 % DMSO buffer  using the PUREGENE DNA Purification Kit protocol (Qiagen, Valencia, CA). We sequenced one mitochondrial and three nuclear loci: (1) mitochondrial ND2 gene and flanking tRNAs (NADH dehydrogenase subunit 2: 1144 bp), (2) POMC (Pro-opiomelanocortin: exon 629 bp), (3) C-myc (cellular myelocytomatosis proto-oncogene: exons and intron 1335 bp), and (4) Rag-1 (recombination activating gene: exon 1282 bp). Primers are the same as in Lawson (2010) with the addition of Rag-1 primers designed for this study (Table 5, ). Extraction, amplification, sequencing, and cloning of alleles follow Lawson (2010) .
Alignments were performed using MUSCLE  as described in . Alignments were unambiguous in ND2, POMC, and Rag-1, and easily aligned in C-myc despite the indel-prone intron region. All protein-coding regions were translated into amino acids to check for errors. Substitution rates for each locus are described in Lawson (2010) : ND2 0.00957/lineage/my ; C-myc 0.0006334/lineage/my ; POMC 0.000721/lineage/my ; and Rag-1 0.00042/lineage/my , and used for approximate estimates of diversification events. One nuclear allele from each individual (either cloned or phased as outlined above) was randomly discarded to create a 1:1 mtDNA:nuDNA dataset. Summary information for each locus shown in Table 1.
Phylogenetic reconstruction and divergence time estimation
Sequences from  and 10 additional East African Hyperolius and Afrixalus were included in initial analyses to assess monophyly of the spiny-throated reed frogs (Additional files) including an additional proposed close relative, H. parkeri , using phylogenetic methods described below. A number of datasets were used for initial assessment of phylogenetic relationships: species tree (single representative per species and all outgroups), population/species tree (all spiny-throated specimens (ingroup) and all outgroups), and population-level analysis of only spiny-throated specimens. Monophyly of the ingroup individuals was strongly supported in all cases. Final Bayesian inference (BI) and maximum likelihood (ML) phylogenetic analyses were then completed with the full ingroup dataset (all spiny-throated specimens). RAxML (ML) and BPP (BI) analysis included a single arbitrary close outgroup (H. mitchelli), while the BEAST (BI) and *BEAST (BI) rate-calibrated trees were completed without an outgroup for improved precision of branch length estimates. RAxML (version 7.0.4)  analyses used the rapid hill climbing algorithm and the GTRGAMMA substitution model  partitioned by gene and codon. BEAST (versions 1.8 and 2.1.3)  and *BEAST (starBEAST) , were partitioned by locus and codon position (SRD06 model) . BPP3 [69, 70] parameters listed below in species delimitation. For both ML and BI analyses, model parameters were independently optimized for each partition. All analyses were run for 10 million generations, and run twice to ensure stability (results not combined). The first 10 % of total generations were discarded as burnin for both convergence and tree estimates. Convergence was investigated using Tracer (version 1.6)  through a visual inspection of adequate mixing and ESS estimates >200. The maximum clade credibility tree was calculated for BEAST and *BEAST trees using TreeAnnotator in BEAST. ML node support was evaluated by non-parametric bootstrapping  with 1000 replicates performed in RAxML. The Majority Rule Consensus tree from BPP3 is reported from Evolver in PAML 4 .
BEAST analysis was run with a coalescent, constant size tree-prior and a strict molecular clock (as recommended for recent population-level analyses). *BEAST and BPP3 analysis used individuals from the same mountain block as discrete units (after confirmation of monophyly from individual-based tree constructions) with the exception of a single unit containing individuals from the East Usambara Mountains and Nguru Mountains which were not mutually monophyletic. A Birth-Death tree prior was used in *BEAST due to a priori hypotheses that this lineage has likely undergone local population extinctions through its history.
The program BPP3 [69, 70] was used to jointly infer the species tree and species limits. Default priors and settings were used: gamma prior = G(2, 1000) with mean 2/2000 = 0.001 for population size parameters (s, θ); gamma prior G(2, 1000) for the species tree root age (τ0); all other divergence time parameters are assigned the Dirichlet prior  equation two). Analyses were also run with a very large (1,10) and a very small mean (1, 10000) for θ and τ0 to ensure stability of results. A heredity file (G(4,5)) and locus rate file were incorporated to account for the combined mitochondrial and nuclear datasets and associated mutation rates. Each analysis was run twice to confirm consistency between runs. Species delimitation was assessed with rjMCMC algorithm 0, e = 2.5.
To establish species boundaries within the spiny-throated reed frogs, a Bayesian General Mixed Yule-Coalescent (bGMYC) model was implemented in the package bGMYC  in R v. 3.0.3  using 100 random trees from the concatenated BEAST. Simulations were set at 50,000 generations with 40,000 burn-in, sampling every 100th generation. The upper threshold for the number of species was set at 40, to match the number of tips on the tree.
Bayes Factor species Delimitation (BFD; ) was used to compare three alternative species scenarios using stepping stone and pathsampling analysis in *BEAST summarized in Table 2. Based on the depth of nodes and morphological divergence, H. tanneri, H. minutissimus, and H. ukwiva are all considered separate species in all analyses. All *BEAST parameters are as above. These approaches were used to define units for analyzing molecular, morphological and ecological parameters among lineages.
Effective population size
BPP3 was used to estimate coalescent-scaled population sizes (θ = 4Neμ) and time of divergence (τ = μ t) throughout this radiation on the fixed species tree recovered from all phylogenetic analyses. The parameters used are the same as in the joint species delimitation and species tree analysis.
To test for signs of demographic expansion or contraction we implemented multilocus extended Bayesian skyline plots (EBSP) in BEAST v. 1.8 on populations with n > 5 (Mulanje Massif, Malawi; East Usambara Mountains, Tanzania; Nguru Mountains, Tanzania; Udzungwa Mountains, Tanzania) and all individuals within the “spinigularis” clade to calculate population size through time and the probability of bottlenecks or expansion. All EBSP analyses used a strict clock. Operators were modified according to author recommendations and analyses were run for 10 million generations to obtain adequate ESS values.
To determine morphological distinctiveness among members of the species cluster, we took commonly-used limb and cranial measurements for 17 traits , including leg bone (e.g. femur, tibia) and foot lengths and head width (see below). Many of these traits are often assumed to be adaptive, as they are correlated with food source (e.g., head width), perching (e.g., foot length), and other life history traits. If species are morphologically divergent, it is possible that competition or adaptation to different environmental conditions reinforced separation within this radiation. Morphological measurements made in this analysis include: Snout-Urostyle Length (SUL), Head Width (HW), Head Length Diagonal from corner of mouth (HLD), Head Length Diagonal from jawbone end (HLDJ), Nostril-Snout (NS), Inter-narial (IN), Eye to Nostril (EN), Eye Distance (EE), Inter-orbital (IO), Tibiafibula Length (TL), Thigh Length (THL), Tibiale Fibulare Length (TFL), Foot Length (FL), Forelimb Length (FLL), Hand Length (HL), Width of Gular Flap (WGF) Height of Gular Flap (HGF). We also scored specimens for traits associated with sexual selection including presence and distribution of gular spines (which are only present during the breeding season for species that possess them) and color patterns. Morphological measurements were taken of 114 mature specimens using Mitutoyo Absolute Digimatic Calipers (CD-6"C) (Sample sizes and Measurements in supplementary materials). In order to assess the morphological distinctness of these species, we conducted Principal Component analyses on log-transformed data of males and females separately with all variables centered and scaled in R  (95 % confidence ellipse probability threshold from ggbiplot package). To further evaluate whether species were significantly different, we did a Permutational Multivariate Analysis of Variance Using Distance Matrices using the adonis function from the vegan  R package.
Kruskal–Wallis and ANOVA were used to assess variable divergence (standardized by body size) between species for use as diagnostic characters in R.
Niche similarity of all species was assessed by Principle Component Analysis (PCA) using bioclim variables associated with GPS coordinates (Additinal file 1: Table S1) in the MASS and ggbiplot packages in R (95 % confidence ellipse probability threshold). This method allows at least a preliminary view of habitat similarity for severely range-restricted taxa such as the members of this complex. We evaluated both a full bioclim dataset (all 19 standard bioclim variables) and a reduced dataset with Pearson’s correlation coefficients below 0.7: Mean Diurnal Range, Temperature Seasonality, Temperature Annual Range, Mean Temperature of Coldest Quarter, Precipitation of Wettest Month, Precipitation Seasonality, Precipitation of Driest Quarter, Precipitation of Warmest Quarter, Precipitation of Coldest Quarter (30 arc-second resolution, WorldClim database ). In reconstructing ecological shifts across the phylogenetic tree, we rely on the assumption that a species niche remains stable through time after divergence. Though this appears to be widely true [16, 81, 82], there is increased uncertainty associated with assessing extremely locally adapted lineages were spatial autocorrelation might have a more pronounced effect.
Range estimates and field observations of population densities
Species’ range estimates were calculated using standardized approaches conducted on species being evaluated for the IUCN Red List of Threatened Species™ . Estimates of population size were made using molecular techniques when sufficient numbers of individuals were available (above). For those species with very few individuals ever found (H. tanneri, H. ukwiva, H. davenporti), estimates were qualitatively assessed from sampling effort over multiple field seasons (Lawson, Menegon, Loader, personal observations, ) in comparison to all known localities from other spiny-throated reed frogs (H. spinigularis, H. burgessi, H. minutissimus).
Mayr E. Animal Species and Evolution. Cambridge, MA: Harvard University Press; 1963.
Lomolino MV, Riddle BR, Whittaker RJ, Brown JH. Biogeography. 4th ed. Sunderland, MA: Sinauer Associates; 2010.
Coyne JA, Orr HA. Speciation. Sinauer Associates Sunderland, MA: Sunderland, MA; 2004.
Zink RM. The Geography of Speciation: Case Studies from Birds. Evo Edu Outreach. 2012;5:541–6.
Hickerson MJ, Stahl EA, Lessios HA. Test for Simultaneous Divergence Using Approximate Bayesian Computation. Evolution. 2006;60:2435–53.
Soltis DE, Morris AB, McLachlan JS, Manos PS, Soltis PS. Comparative phylogeography of unglaciated eastern North America. Mol Ecol. 2006;15:4261–93.
Carstens BC, Richards CL. Integrating Coalescent and Ecological Niche Modeling in Comparative Phylogeography. Evolution. 2007;61:1439–54.
Victoriano PF, Ortiz JC, Benavides E, Adams BJ, Sites Jr JW. Comparative phylogeography of codistributed species of Chilean Liolaemus (Squamata: Tropiduridae) from the central-southern Andean range. Mol Ecol. 2008;17:2397–416.
Hickerson MJ, Meyer CP. Testing comparative phylogeographic models of marine vicariance and dispersal using a hierarchical Bayesian approach. BMC Evol Biol. 2008;8:322.
Kozak KH, Graham CH, Wiens JJ. Integrating GIS-based environmental data into evolutionary biology. Trends Ecol Evol. 2008;23:141–8.
Knowles LL. Statistical phylogeography. Annu Rev Ecol Evol Syst. 2009;40:593–612.
Freckleton RP, Jetz W. Space versus phylogeny: disentangling phylogenetic and spatial signals in comparative data. Proc R Soc B. 2009;276:21–30.
Pigot AL, Phillimore AB, Owens IPF, Orme CDL. The Shape and Temporal Dynamics of Phylogenetic Trees Arising from Geographic Speciation. Syst Biol. 2010;59:660–73.
Losos J, Glor R. Phylogenetic comparative methods and the geography of speciation. Trends Ecol Evol. 2003;18:220–7.
Blair C, Heckman KL, Russell AL, Yoder AD. Multilocus coalescent analyses reveal the demographic history and speciation patterns of mouse lemur sister species. BMC Evol Biol. 2014;14:57.
Peterson AT, Soberón J, Sánchez-Cordero V. Conservatism of Ecological Niches in Evolutionary Time. Science. 1999;285:1265–7.
Schluter D. The Ecology of Adaptive Radiation. New York: Oxford University Press; 2000.
Dieckmann U, Doebeli M, Metz J, Tautz D. Adaptive Speciation. Cambridge, UK: Cambridge Univ. Press; 2004.
Kaneshiro KY. The dynamics of sexual selection and founder effects in speciation formation. In: Kaneshiro KY, Giddings L, Anderson W, editors. Genetics, speciation, and the founder principle. Oxford, UK: Oxford University Press; 1989. p. 279–96.
Panhuis TM, Butlin R, Zuk M, Tregenza T. Sexual selection and speciation. Trends Ecol Evol. 2001;16:364–71.
Ödeen A, Florin A-B. Sexual selection and peripatric speciation: the Kaneshiro model revisited. J Evol Biol. 2002;15:301–6.
Carson HL, Templeton AR. Genetic Revolutions in Relation to Speciation Phenomena: The Founding of New Populations. Annu Rev Ecol Evol Syst. 1984;15:97–131.
Roy MS. Recent diversification in African greenbuls (Pycnonotidae: Andropadus) supports a montane speciation model. Proc R Soc Lond B. 1997;264:1337–44.
Linder HP. The radiation of the Cape flora, southern Africa. Biol Rev. 2003;78:597–638.
Bowie RCK, Fjeldså J, Hackett SJ, Bates JM, Crowe TM. Coalescent models reveal the relative roles of ancestral polymorphism, vicariance, and dispersal in shaping phylogeographical structure of an African montane forest robin. Mol Phylo Evol. 2006;38:171–88.
Lawson LP. The discordance of diversification: evolution in the tropical-montane frogs of the Eastern Arc Mountains of Tanzania. Mol Ecol. 2010;19:4046–60.
Linder HP, de Klerk HM, Born J, Burgess ND, Fjeldså J, Rahbek C. The partitioning of Africa: statistically defined biogeographical regions in sub-Saharan Africa. J Biogeogr. 2012;39:1189–205.
Demos TC, Kerbis Peterhans JC, Agwanda B, Hickerson MJ. Uncovering cryptic diversity and refugial persistence among small mammal lineages across the Eastern Afromontane biodiversity hotspot. Mol Phylo Evol. 2014;71:41–54.
Ceccarelli FS, Menegon M, Tolley KA, Tilbury CR, Gower DJ, Laserna MH, Kasahun R, Rodriguez-Prieto A, Hagmann R, Loader SP. Evolutionary relationships, species delimitation and biogeography of Eastern Afromontane horned chameleons (Chamaeleonidae: Trioceros). Mol Phylo Evol. 2014;80:125–36.
Poynton JC, Loader SP, Sherratt E, Clarke BT. Amphibian Diversity in East African Biodiversity Hotspots: Altitudinal and latitudinal Patterns. Biodivers Conserv. 2007;16:1103–18.
Fjeldså J, Bowie RCK, Rahbek C. The Role of Mountain Ranges in the Diversification of Birds. Annu Rev Ecol Evol Syst. 2012;43:249–65.
Hua X, Wiens JJ. Latitudinal Variation in Speciation Mechanisms in Frogs. Evolution. 2010;64:429–43.
Kozak KH, Wiens JJ. Climatic zonation drives latitudinal variation in speciation mechanisms. Proc R Soc B. 2007;274:2995–3003.
Lovett JC. Eastern Arc moist forest flora. In: Biogeography and ecology of the rain forests of eastern Africa. Cambridge, UK: Cambridge University Press; 1993. p. 33–55.
De Menocal PB. Plio-Pleistocene African Climate. Science. 1995;270:53–9.
Meadows ME, Linder HP. A Palaeoecological Perspective on the Origin of Afromontane Grasslands. J Biogeogr. 1993;20:345–55.
Hamilton AC. The quaternary history of African forests: its relevance to conservation. Afr J Ecol. 1981;19:1–6.
Lovett JC, Marchant R, Taplin J, Küper W. The oldest rainforests in Africa: stability or resilience for survival and diversity. In: Purvis A, Gittleman JL, Brooks TM, editors. Phylogeny and conservation. Cambridge, UK: Cambridge University Press; 2005. p. 115–25.
Marchant R, Mumbi C, Behera S, Yamagata T. The Indian Ocean dipole – the unsung driver of climatic variability in East Africa. Afr J Ecol. 2007;45:4–16.
Mumbi CT, Marchant R, Hooghiemstra H, Wooller MJ. Late Quaternary vegetation reconstruction from the Eastern Arc Mountains, Tanzania. Quat Res. 2008;69:326–41.
Loader SP, Sara Ceccarelli F, Menegon M, Howell KM, Kassahun R, Mengistu AA, Saber SA, Gebresenbet F, de Sá R, Davenport TRB, Larson JG, Müller H, Wilkinson M, Gower DJ. Persistence and stability of Eastern Afromontane forests: evidence from brevicipitid frogs. J Biogeogr. 2014;41:1781–92.
Crisp M, Cook L, Steane D. Radiation of the Australian flora: what can comparisons of molecular phylogenies across multiple taxa tell us about the evolution of diversity in present–day communities? Phil Trans R Soc Lond B. 2004;359:1551–71.
Cantalapiedra JL, Hernández Fernández M, Morales J. Biomic Specialization and Speciation Rates in Ruminants (Cetartiodactyla, Mammalia): A Test of the Resource-Use Hypothesis at the Global Scale. PLoS One. 2011;6, e28749.
Johnson SD, Linder HP. Speciation and extinction in the Greater Cape Floristic Region. In: Allsopp N, Colville JF, Verboom GA, editors. Fynbos: Ecology, Evolution, and Conservation of a Megadiverse Region. Oxford, UK: Oxford University Press; 2014. p. 119.
Loader SP, Lawson LP, Portik DM, Menegon M. Three new species of spiny throated reed frogs (Anura: Hyperoliidae) from Evergreen forests of Tanzania. in press. BMC Res Notes. 8(1):167
Schiøtz A. Two New Hyperolius (Anura) from Tanzania. Steenstrupia. 1982;8:269–76.
Schiøtz A. Treefrogs of Africa. Chimairath ed. Germany: Frankfurt am Main; 1999.
Schiøtz A, Westergaard MM. Notes on some Hyperolius (Anura: Hyperoliidae) from Tanzania, with supplementary information on two recently described species. Steenstrupia. 2000;25:1–9.
Portik DM, Mulungu EA, Sequeira D, McEntee JP. Herpetological surveys of the Serra Jeci and Namuli massifs, Mozambique, and an annotated checklist of the southern Afromontane archipelago. Herpetol Rev. 2013;44:394–406.
Cerling TE, Bowman JR, O’Neil JR. An isotopic study of a fluvial-lacustrine sequence: The Plio-Pleistocene koobi fora sequence, East Africa. Palaeogeogr Palaeoclimatol Palaeoecol. 1988;63:335–56.
Trauth MH, Larrasoaña JC, Mudelsee M. Trends, rhythms and events in Plio-Pleistocene African climate. Quat Sci Rev. 2009;28:399–411.
IUCN SSC Amphibian Specialist Group 2013. Hyperolius spinigularis/Hyperolius minutissimus. The IUCN Red List of Threatened Species. Version 2014.3. <www.iucnredlist.org> Downloaded on 01 June 2015.
Rovero F, Menegon M, Fjeldså J, Collett L, Doggart N, Leonard C, Norton G, Owen N, Perkin A, Spitale D, Ahrends A, Burgess ND. Targeted vertebrate surveys enhance the faunal importance and improve explanatory models within the Eastern Arc Mountains of Kenya and Tanzania. Diversity Distrib. 2014;20:1438–49.
Rovito SM. Lineage divergence and speciation in the Web-toed Salamanders (Plethodontidae: Hydromantes) of the Sierra Nevada, California. Mol Ecol. 2010;19:4554–71.
Rovito SM, Wake DB, Papenfuss TJ, Parra-Olea G, Muñoz-Alonso A, Vásquez-Almazán CR. Species formation and geographical range evolution in a genus of Central American cloud forest salamanders (Dendrotriton). J Biogeogr. 2012;39:1251–65.
Smith TB, Wayne RK, Girman DJ, Bruford MW. A Role for Ecotones in Generating Rainforest Biodiversity. Science. 1997;276:1855–7.
Seutin G, White BN, Boag PT. Preservation of avian blood and tissue samples for DNA analyses. Can J Zool. 1991;69:82–90.
Hoegg S, Vences M, Brinkmann H, Meyer A. Phylogeny and Comparative Substitution Rates of Frogs Inferred from Sequences of Three Nuclear Genes. Mol Biol Evol. 2004;21:1188–200.
Wiens JJ, Fetzner JW, Parkinson CL, Reeder TW. Hylid frog phylogeny and sampling strategies for speciose clades. Sys Biol. 2005;54:719–748.
Edgar RC. MUSCLE: Multiple Sequence Alignment with High Accuracy and High Throughput. Nucl Acids Res. 2004;32:1792–7.
Crawford AJ. Relative Rates of Nucleotide Substitution in Frogs. J Mol Evol. 2003;57:636–41.
Evans BJ, McGuire JA, Brown RM, Andayani N, Supriatna J. A coalescent framework for comparing alternative models of population structure with genetic data: evolution of Celebes toads. Biol Lett. 2008;4:430–3.
Pickersgill M. Frog Search, Results of Expeditions to Southern and Eastern Africa. Chimairath ed. Germany: Frankfurt am Main; 2007.
Stamatakis A. RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models. Bioinformatics. 2006;22:2688–90.
Stamatakis A, Blagojevic F, Nikolopoulos DS, Antonopoulos CD. Exploring New Search Algorithms and Hardware for Phylogenetics: RAxML Meets the IBM Cell. J VLSI Sign Process Syst Sign Im. 2007;48:271–86.
Drummond AJ, Rambaut A. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol. 2007;7:1–8.
Heled J, Drummond AJ. Bayesian inference of species trees from multilocus data. Mol Biol Evol. 2010;27:570–80.
Shapiro B, Rambaut A, Drummond AJ. Choosing Appropriate Substitution Models for the Phylogenetic Analysis of Protein-Coding Sequences. Mol Biol Evol. 2006;23:7–9.
Rannala B, Yang Z. Bayes estimation of species divergence times and ancestral population sizes using DNA sequences from multiple loci. Genetics. 2003;164:1645–56.
Yang Z, Rannala B. Bayesian species delimitation using multilocus sequence data. PNAS. 2010;107:9264–9.
Rambaut A, Drummond AJ. Tracer ver. 1.5. Program available at 2007; http://tree.bio.ed.ac.uk/software/tracer/.
Felsenstein J. Phylogenies and the comparative method. Am Nat. 1985;125:1–15.
Yang Z. PAML 4: Phylogenetic Analysis by Maximum Likelihood. Mol Biol Evol. 2007;24:1586–91.
Reid NM, Carstens BC. Phylogenetic estimation error can decrease the accuracy of species delimitation: a Bayesian implementation of the general mixed Yule-coalescent model. BMC Evol Biol. 2012;12:196.
R Development Core Team. R: A language and environment for statistical computing. Austria: R Foundation for Statistical Computing Vienna; 2011. 3.
Grummer JA, Bryson RW, Reeder TW. Species Delimitation Using Bayes Factors: Simulations and Application to the Sceloporus scalaris Species Group (Squamata: Phrynosomatidae). Syst Biol. 2014;63:119–33.
Moen DS, Wiens JJ. Phylogenetic Evidence for Competitively Driven Divergence: Body-Size Evolution in Caribbean Treefrogs (hylidae: Osteopilus). Evolution. 2009;63:195–214.
Wickham H. ggplot2: Elegant Graphics for Data Analysis. New York: Springer; 2009.
Oksanen J, Blanchet FG, Kindt R, Legendre P, Minchin PR, O’Hara RB, Simpson GL, Solymos P, Stevens MHH, Wagner H. Vegan: Community Ecology Package. R Package Version 2.0-10. 2013.
Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005;25:1965–78.
Wiens JJ, Ackerly DD, Allen AP, Anacker BL, Buckley LB, Cornell HV, Damschen EI, Jonathan Davies T, Grytnes J-A, Harrison SP, Hawkins BA, Holt RD, McCain CM, Stephens PR. Niche conservatism as an emerging principle in ecology and conservation biology. Ecol Lett. 2010;13:1310–24.
Martínez-Meyer E, Townsend Peterson A, Hargrove WW. Ecological niches as stable distributional constraints on mammal species, with implications for Pleistocene extinctions and climate change projections for biodiversity. Global Ecol Biogeogr. 2004;13:305–14.
We would like to thank J. Miller, H. Müller, D. Blackburn, D. Portik, E. Tonelli, C. Barratt, T. Price, J. Luedtke, & L. Coals for tissues, data and comments on early versions of this manuscript. Also thanks to M. Munissi, E. Mulungu, Nyson, L. Mahler, and B. Zimkus for help in the field. K. Howell, D. Moyer, L. Mazibuko, and J. Poynton provided extensive knowledge of African Amphibians and collecting logistics in Tanzania and Malawi. P. Kaliba was instrumental for fieldwork in Malawi. Special thanks to TAWIRI, COSTEC, Tanzanian Wildlife and Forestry departments for research permits (RCA 2006, 2007-72-Na-2006-19, RCA 2001–272; RCA 2007–153, RCA 2009-306-NA-2009-201, 2011-239-NA-2011-82), and Museums of Malawi for Malawian permits. Thanks to K. Feldheim and the Pritzker lab for help on molecular work and for creating the taq used in this study. Lawson’s PhD work was funded through a GAANN educational grant to the Committee on Evolutionary Biology, University of Chicago. This project was funded through National Science Foundation DDEP grant (NSF 0841641), Hinds Fund of the University of Chicago, and grants through the Council on Africa (Field Museum of Natural History). SPL was funded to conduct surveys and lab work by the following institutes: University of Basel, Swiss National Science Foundation (grant number 31003A-133067), SCNAT, Freiwillige Akademische Gesellschaft Basel, and the Council on Africa (Field Museum of Natural History). Sponsors were not involved in collection, analysis, or interpretation of data, or writing and submitting this manuscript. Finally, we thank the editor (John Wiens) and reviewers for helpful and constructive criticisms that improved the manuscript.
The authors declare that they have no competing interests.
LPL conceived of the study, field collected frogs, sequenced DNA, carried out molecular genetics analyses, morphological analyses, and spatial/niche analyses, and drafted the manuscript. JMB advised on analyses and supported molecular work. MM field collected frogs, calculated range sizes and potential distributions, and advised on Udzungwa identifications in the literature. SPL field collected frogs, sequenced DNA, measured frogs, and with LPL drafted the manuscript. All authors read and approved the final manuscript.
List of specimens analyzed including locality and genbank accession numbers for all loci.
Maximum Clade Credibility Gene Trees from BEAST. Posterior probabilities shown on relevant nodes. Figure S2. Species tree with BPP3 posterior theta values shown for species splits and current populations: Mean (95 % interval) with three theta priors. Top (large mean) = (1, 10). Middle = (2, 1000). Bottom (very small mean) = (1, 10000). Populations IDs match Table 2. Figure S3. PCA analysis of bioclim variables with all members of the “spinigularis” clade (H. burgessi, H. davenporti, H. spinigularis) combined to contrast from H. tanneri. Hyperolius tanneri specimens do not both fall outside of the 95 % confidence interval for the H. spinigularis members, thus no habitat shift is inferred. Table S2. Summary of BPP species delimitation results from 10000 delimitations. Populations were defined as: all individuals from East Usambara and Nguru Mountains (Eung), Uluguru Mountains (Ul), Southern Highlands (SH), West Usambara Mountains (Wu), Udzungwa Mountains (Ud), Rubeho Mountains (Ru), Mulanje Massif (Mu), Mount Namuli (Na), and the outgroup Hyperolius mitchelli (Mit). Table S3. PCA analyses of morphology for Females and Males. Table S4. PCA from reduced set of Bioclim variables. Table S5. jModeltest scores for each locus. Top 8 models from each analysis shown. 24 models considered to facilitate utilization in BEAST. Selected models are in bold and outlined.