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

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 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, 58]). Extraction, amplification, sequencing, and cloning of alleles follow Lawson (2010)
26].

Table 5. Primers use in this study

Alignments were performed using MUSCLE 60] as described in 26]. 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) 26]: ND2 0.00957/lineage/my 61]; C-myc 0.0006334/lineage/my 26]; POMC 0.000721/lineage/my 26]; and Rag-1 0.00042/lineage/my 62], 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 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. parkeri63], 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)
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 spiny-throated 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: 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 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).