### 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 [38] 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 analyses

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.

### Laboratory procedures

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)

### Molecular sexing

We sexed individuals with the primer pair P2/P8 using the protocol described in Griffiths et al. [39]; 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, [40]). 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).

### Phylogenetic reconstruction

Molecular phylogenies were estimated using maximum likelihood and Bayesian inference, as implemented in PHYML 3.0 [41] and MRBAYES 3.1.2 [42], respectively. The most appropriate models of nucleotide substitution were determined with DT_MODSEL [43]. 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 [44] 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 [45]. 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 [46], as implemented in DNASP 5.0 [47]. 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 [48] 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 [49] 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 [50] 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}, [51]), as implemented in ARLEQUIN v3.1 [52]. 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 [53].

We used Fu's *Fs* test (1000 replicates) and Ramos-Onsins and Rozas R_{2} statistic [24], as implemented in DNASP 5.0, to detect signatures of demographic change. The significance of the R_{2} statistic was assessed using 1000 coalescent simulations. We also used Bayesian skyline plots [54] 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. [55] using the HKY model and a strict molecular clock. The MCMC simulations were run for 20^{6} iterations, with genealogies and model parameters being sampled every 1000 iterations. The Bayesian skyline plots (BSPs) were visualized in TRACER V.1.5 [44].

### 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. [56]] although its magnitude continues to be debated [57, 58]. Subramanian et al. [59] 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. [57] 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. [61] and the 'standard' avian mtDNA rate of 2.1% divergence per million years. MCMC chains were run for 5*10^{6} steps and were sampled every 1000 steps.

### Population assignment using STRUCTURE

We used STRUCTURE V2 [62] 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*10^{6} iterations (burnin: 2*10^{5} iterations) from K = 1 to K = 5. The number of clusters (populations) was estimated using ΔK [63].

### Multi-locus network and species tree approaches

We used POFAD V1.03 [64] and SPLITSTREE V4.0 [65] 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, [66]) 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[43]. 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 [[68]: 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 [69] to fit the data to a model that included both isolation and migration. IMa estimates six parameters scaled to the neutral mutation rate (μ): θ_{pop1} (4Ne_{pop1}μ), θ_{pop}2 (4Ne_{pop2}μ), θ_{popA} (4Ne_{popA}μ), 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 [70], and a mutation rate of 1.05*10^{-8} substitutions/site/year (s/s/y) for mtDNA (6.1% per million years, [61]), 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 [71]. 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*
_{Usambara}.

### 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 [74]. 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 [75].