Divergence at the edges: peripatric isolation in the montane spiny throated reed frog complex

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.


Background
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 [3] (see below). The classical model of vicariant allopatry starts with an ancestral distribution, which is split into large, disjointed populations [2]; 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][6][7][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 [9]. 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 [3]. New analytical approaches using multiple types of data enable testing for these processes in greater detail [10][11][12][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 [14]. 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][17][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][20][21][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][24][25][26][27][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][30][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 spinythroated reed frogs within the EABH and assess regionalscale 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 [34]. 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 [35]. 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 [36]. 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 [37]; see [38][39][40][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 [44].
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 spinythroated complex is comprised of six taxa (H. burgessi, H. davenporti, H. spinigularis, H. tanneri, H. minutissimus, and H. ukwiva) [45], with species occupying localities often restricted to either montane mosaic or rainforest habitats [26,[46][47][48][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

Phylogenetic relationships
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 [45]). 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 [48] 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
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 [45] (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

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).

Morphological divergence
Measurements for all traits and species summary statistics are presented in [45]. 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, headwidth, 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, [45]). Presence and distribution of gular and ventral spines are different between all species [45] and is a diagnostic tool for separating species (see [45]).
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 [50], before a general shift to grasslands at 1.86 MYA [51].
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 [52].

Discussion
New lineages within the spiny-throated reed frog clade Prior to this study and our complementary taxonomic revision [45], 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 [53]. 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 sexuallydimorphic 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 spinythroated 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 [34], 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 spinythroated hyperoliids found across these habitatsit seems that peripatric processes may be more pronounced in higher elevation species compared to those in more persistent rainforests (e.g., [41]). 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 biomesas 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 periand parapatric processes [56]; 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.

Conclusions
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.

Molecular data
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 [57] using the PURE-GENE DNA Purification Kit protocol (Qiagen, Valencia, CA). We sequenced one mitochondrial and three nuclear loci: (1) (Table 5, [58]). Extraction, amplification, sequencing, and cloning of alleles follow Lawson (2010) [26].

Phylogenetic reconstruction and divergence time estimation
Sequences from [26] 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 [63], 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) ratecalibrated trees were completed without an outgroup for improved precision of branch length estimates. RAxML (version 7.0.4) [64] analyses used the rapid hill climbing algorithm and the GTRGAMMA substitution model [65] partitioned by gene and codon. BEAST (versions 1.8 and 2.1.3) [66] and *BEAST (starBEAST) [67], were partitioned by locus and codon position (SRD06 model) [68]. 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) [71] 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 [72] with 1000 replicates performed in RAxML. The Majority Rule Consensus tree from BPP3 is reported from Evolver in PAML 4 [73]. 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.

Species delimitation
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 [70] 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 spinythroated reed frogs, a Bayesian General Mixed Yule-Coalescent (bGMYC) model was implemented in the package bGMYC [74] in R v. 3.0.3 [75] 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; [76]) 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.

Morphological divergence
To determine morphological distinctiveness among members of the species cluster, we took commonly-used limb and cranial measurements for 17 traits [77], 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: 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 [78] (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 [79] 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 divergence
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 [80]). 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™ [52]. 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, [53]) in comparison to all known localities from other spiny-throated reed frogs (H. spinigularis, H. burgessi, H. minutissimus).