Tissue and morphological data collection
Our morphometric data was drawn from 29 specimens deposited in the Zoological Museum, University of Copenhagen, as well as from notes on an additional 61 specimens deposited at other institutions (Field Museum, National Museum of Natural History, National Museums of Kenya, and Museum Alexander König). All morphological measurements and the scoring of plumage patterns were undertaken by J. Fjeldså: wing-length (flattened-chord, from carpal joint to tip of the longest primary feather) was measured to the nearest 0.1 mm using a wing-rule; weight was measured using a Pesola spring scale to the nearest 0.1 g; bill-length (from the bill tip to the base of the exposed bill on the skull), tail-length (from the insertion point on the pygostyle to the tip of the longest feather) and tarsus-length were measured using Vernier calipers to the nearest 0.1 mm. Plumage colors were matched to known color standards  to objectively define colors across specimens.
We collected tissue from 124 individuals (60 albigula, 26 debilis and 38 rabai, altogether 59 vouchers) that span the entire distributional range of the species (Additional File 5, Figure 1): 96 of the samples were fresh tissues and 28, including 24 debilis samples, were obtained from museum specimens (toe-pads).
Morphological measurements were log-transformed (log x+1) to reduce variance between characters. Differences in univariate measures between taxa were tested using paired t-tests. Principal component analysis (PCA) was used to summarize morphological variation. PCA was performed on individuals and only principal components (PCs) with eigenvalues greater than one were extracted. The factor matrix was rotated using the varimax method to optimize variable loadings. The resulting rotated factor matrix was used to determine which of the original variables were most highly correlated with the PCs.
DNA was extracted from tissue or blood using the Qiagen extraction kit (Valencia, CA) following the manufacturer's protocol. We extracted the DNA from toe-pads in a room dedicated to ancient DNA labwork. We used the same extraction protocols as for the fresh samples, but added 20 μl of dithiothreitol (DTT, 0.1 M) to facilitate the digestion of these tissues. We amplified and sequenced four loci, one mitochondrial (ND2), two autosomal (Beta-fibrinogen intron-5, FGB; GAPDH intron-11, GAPDH) and one Z-linked (BRM intron-15, BRM). PCR-amplification was performed using standard protocols, only the annealing temperature varied (54-60°C). Locations on the chicken genome and primer sequences used are detailed in the Additional File 6. All sequences have been deposited in GenBank (Accession Numbers HQ716721-HQ717145)
We sexed individuals with the primer pair P2/P8 using the protocol described in Griffiths et al. ; sex was deduced based on the number of bands (two for females and one for male). For the museum specimens, sex was inferred from the specimen labels. We could not obtain any PCR-product using the P2/P8 primer pair for 14 individuals that were homozygous at the Z-linked locus we sequenced (BRM). For these 14 individuals, we aimed to PCR-amplify and sequence three further Z-linked loci (CHDZ, ACO1, SPIN1, ). Our success with PCR-amplification was variable; we considered an individual to be a female if no heterozygous position was detected in at least two loci (minimum length of Z-linked data: 834 bp, maximum length of Z-linked data: 2849 bp).
Molecular phylogenies were estimated using maximum likelihood and Bayesian inference, as implemented in PHYML 3.0  and MRBAYES 3.1.2 , respectively. The most appropriate models of nucleotide substitution were determined with DT_MODSEL . Two analyses of four Metropolis-coupled MCMC chains (one cold and three heated) were run for ten million iterations with trees sampled every 100 iterations. The number of iterations discarded before the posterior probabilities were calculated varied among analyses. We checked that the potential scale reduction factor (PSRF) approached 1.0 for all parameters and that the average standard deviation of split frequencies converged towards zero. We also used TRACER v1.5  to ascertain whether our sampling of the posterior distribution had reached a sufficient effective sample size (ESS) for meaningful parameter estimation. We made use of sequences of P. hypochloris (DQ402215) and P. xavieri (DQ402219) as outgroups.
Testing for selection
To test whether selection was acting on the mitochondrial protein-coding gene (ND2), we used the McDonald-Kreitman test . The stop codon was excluded from the analyses, leaving a total of 1038 bp for 110 Tiny Greenbuls. We used sequences of P. hypochloris (DQ402215) and P. xavieri (DQ402219) as outgroups. Significance was assessed using Fischer's exact test at a threshold of 0.05.
We tested for selection acting on the nuclear loci by using the HKA test , as implemented in DNASP 5.0 . We used sequences from one Yellow-Streaked Greenbul Phyllastrephus flavostriatus as an outgroup (S. Lokugalappatti unpubl. data).
Determining the phase of alleles
We used PHASE V2.1.1  to infer the alleles for each nuclear locus. Three runs, using different seed values, were performed and results were compared across runs. We used the recombination model and ran the iterations of the final run 10 times longer than for the other runs. We used a threshold of 0.75 to consider a SNP to be satisfactorily phased  and individuals that did not satisfy this threshold were removed from further analysis.
Population genetic analyses
Haplotype diversity (Hd), nucleotide diversity (π) and Watterson's theta (θ) were estimated with DNASP 5.0. We used TCS 1.21  to reconstruct a 95% statistical parsimony network. Overall genetic structure of populations was investigated with an analysis of molecular variance (AMOVA, 1000 permutations) using pairwise distances (ϕst, ), as implemented in ARLEQUIN v3.1 . In order to test for any correlation between geographic (shortest straight line) and genetic distances (ND2 uncorrected-p) we performed a Mantel test in GENALEX 6 .
We used Fu's Fs test (1000 replicates) and Ramos-Onsins and Rozas R2 statistic , as implemented in DNASP 5.0, to detect signatures of demographic change. The significance of the R2 statistic was assessed using 1000 coalescent simulations. We also used Bayesian skyline plots  on the mitochondrial data set for the albigula and rabai lineages to estimate more complex scenarios of population dynamics. This method is independent of a priori defined demographic models and tree reconstructions, and is thus suitable for taxa with complicated population history. Analyses were run in BEAST 1.5.4.  using the HKY model and a strict molecular clock. The MCMC simulations were run for 206 iterations, with genealogies and model parameters being sampled every 1000 iterations. The Bayesian skyline plots (BSPs) were visualized in TRACER V.1.5 .
mtDNA divergence times
Inferring divergence times within species is a challenging task as no internal fossil calibration points can be used for most species. The existence of a time dependency of substitution rates appears to now be well accepted [e.g. ] although its magnitude continues to be debated [57, 58]. Subramanian et al.  suggested that the time dependency phenomenon could primarily be attributed to non-synonymous substitutions. They estimated the mean rate of evolution at four-fold degenerated sites from complete mtDNA sequences of Adelie Penguins (Pygoscelis adeliae) to be 0.073 (95% HPD: 0.025-0.123 s/s/myr). We applied this mutation rate, and the associated uncertainty, to the four-fold degenerated sites within our mitochondrial data set (192 sites). We used BEAST 1.5.4 [55, 60] with an uncorrelated lognormal molecular clock model, coalescent tree prior with constant population size and a GTR+Γ model of sequence evolution: the same model of sequence evolution used by Subramanian et al.  to estimate the rate. We compared this new molecular rate with a molecular clock rate for ND2 (6.2%/Ma) proposed by Arbogast et al.  and the 'standard' avian mtDNA rate of 2.1% divergence per million years. MCMC chains were run for 5*106 steps and were sampled every 1000 steps.
Population assignment using STRUCTURE
We used STRUCTURE V2  to infer how many populations could be distinguished based on the three nuclear loci. We only included individuals (n = 80) for which: 1) all three nuclear loci could be satisfactorily phased, and 2) sequences of the three nuclear loci were available. We compared the optimal number of populations estimated by STRUCTURE and the probability of each individual being assigned to each population across analyses. We assumed an admixture model with correlated allele frequencies and let alpha vary among populations. We ran 2*106 iterations (burnin: 2*105 iterations) from K = 1 to K = 5. The number of clusters (populations) was estimated using ΔK .
Multi-locus network and species tree approaches
We used POFAD V1.03  and SPLITSTREE V4.0  to build a multi-locus network based on the mitochondrial and nuclear data sets. We used the same set of individuals that were used in the STRUCTURE analyses. We used uncorrected-p distances as input for POFAD and made use of the standardized matrix for network reconstruction.
We used the species tree approach (*BEAST, ) implemented in BEAST 1.5.4 [55, 60], to estimate the relationship among the three subspecies (albigula, debilis and rabai). Species tree approaches implement the coalescent to estimate a species tree based on the individual gene trees; this approach has been shown to outperform the traditional concatenation approaches in that incomplete lineage sorting is explicitly taken into account. We defined each subspecies as a 'species'. We used all individuals for which full phase information was available and assumed a strict molecular clock model for all loci and made used of the best-fit model for each partition, as determined with DT_MODSEL. Thus, each locus had its own model and clock rate specified. We ran the MCMC chains for 100 million iterations.
We made use of the software BAYESIAN PHYLOGENY AND PHYLOGEOGRAPHY V2.0 (BPP v2.0, [67, 68]) to estimate the speciation probability for the case where all three subspecies are considered species. This decision was based on the results from the STRUCTURE analyses, multi-locus network and species tree analyses. The method implemented in BPP v2.0 accommodates the species phylogeny as well as lineage sorting due to ancestral polymorphism. A speciation probability of 1 on a node indicates that every species delimitation model visited by the reverse-jump MCMC algorithm supports the marginal posterior probability inference that two lineages descend from a particular node as 'species'. We consider to speciation probability values greater than 0.95 as strong-support for the occurrence of a speciation event. A gamma prior was used on the population size parameters (θs) and the age of the root in the species tree (τ0), whereas the other divergence time parameters were assigned a Dirichlet prior [: equation 2]. We ran the rjMCMC analyses for 500 000 generations with a burn-in period of 10 000 and different starting seeds. We ran each analysis at least twice. We evaluated the influence of the priors on the posterior probabilities by changing the priors for θ and τ0, assuming either small or large ancestral population size with G set to (2, 2000) and (1, 10), respectively, and shallow or deep divergence with G set to (2, 2000) and (1, 10), respectively.
Coalescent analyses using the Isolation with Migration model
We used the Markov chain Monte Carlo method implemented in the program IMa  to fit the data to a model that included both isolation and migration. IMa estimates six parameters scaled to the neutral mutation rate (μ): θpop1 (4Nepop1μ), θpop2 (4Nepop2μ), θpopA (4NepopAμ), t (T/μ, where T is the time since population divergence in years before present), m1 (2 M/θpop1, where M is the effective number of migrants moving into population 1) and m2 (2 M/θpop2, where M is the effective number of migrants moving into population 2). We defined inheritance scales to reflect the difference in modes of inheritance among the loci used: 0.25 for the mtDNA locus, 0.75 for the Z-linked locus and 1.0 for the two autosomal loci. We used an HKY model of nucleotide substitution for the mtDNA and an infinite-sites model for each of the three nuclear loci. Parameters and genealogies were sampled every 100 steps for the 5 and 10 million step runs. To assess convergence we monitored the extent of autocorrelation, parameter trend lines and the effective sample size (ESS) throughout the run and we also compared the results between three independent runs.
We used a generation time of 1.7 years, which reflects the average for several passerine species , and a mutation rate of 1.05*10-8 substitutions/site/year (s/s/y) for mtDNA (6.1% per million years, ), and thus a per locus rate of 3.17*10-5 substitutions per year. For the Z-linked locus we assumed a mutation rate of 3.6*10-9 s/s/y and for autosomal loci we selected a rate of 3.61*10-9 s/s/y . This translated into per locus rates of: FGB 2.01*10-6 s/l/y, GAPDH 1.17*10-6 s/l/y, and BRM 1.31*10-6 s/l/y. The geometric mean of the combined mtDNA and ncDNA was 3.14*10-6 s/l/y. We tested for intralocus recombination using the GARD and SBP algorithm, as implemented in HYPHY [72, 73].
We performed four pairwise comparisons: 1) albigula versus rabai, 2) albigula versus rabai with individuals collected on Mt Kanga excluded (see results), 3) rabai versus debilis, and 4) albigula
Nguru versus albigula
Bioclimatic data analyses
Current bioclimatic data (Bioclim variables 1-18, 30 arc-seconds resolution = c. 1 × 1 km) were collected from WORLDCLIM http://www.worldclim.org/current and compiled using DIVA-GIS v7.2 . The climatic envelope for each of the sampling points was then extracted. We performed a Principal Component Analysis on the 18 bioclimatic parameters extracted from each sampling point using a correlation matrix approach in R 2.10.1 .