Molecular Ecology (2015) 24, 1172–1187

doi: 10.1111/mec.13105

Post-fragmentation population structure in a cooperative breeding Afrotropical cloud forest bird: emergence of a source-sink population network M . H U S E M A N N , * 1 L . C O U S S E A U , † 1 T . C A L L E N S , † E . M A T T H Y S E N , ‡ C . V A N G E S T E L , †§ C . H A L L M A N N ¶ and L . L E N S † * * *General Zoology, Institute of Biology, Martin-Luther University Halle-Wittenberg, Halle (Saale), Germany, †Terrestrial Ecology Unit, Department of Biology, Ghent University, Ghent, Belgium, ‡Evolutionary Ecology Group, Department of Biology, University of Antwerp, Antwerp, Belgium, §Entomology Department, Royal Belgian Institute of Natural Sciences, Brussels, Belgium, ¶Sovon, Dutch Centre for Field Ornithology, PO Box 6521, 6503 GA, Nijmegen, Netherlands, **Ornithology Section, Department of Zoology, National Museums of Kenya, Nairobi, Kenya

Abstract The impact of demographic parameters on the genetic population structure and viability of organisms is a long-standing issue in the study of fragmented populations. Demographic and genetic tools are now readily available to estimate census and effective population sizes and migration and gene flow rates with increasing precision. Here we analysed the demography and genetic population structure over a recent 15year time span in five remnant populations of Cabanis’s greenbul (Phyllastrephus cabanisi), a cooperative breeding bird in a severely fragmented cloud forest habitat. Contrary to our expectation, genetic admixture and effective population sizes slightly increased, rather than decreased between our two sampling periods. In spite of small effective population sizes in tiny forest remnants, none of the populations showed evidence of a recent population bottleneck. Approximate Bayesian modelling, however, suggested that differentiation of the populations coincided at least partially with an episode of habitat fragmentation. The ratio of meta-Ne to meta-Nc was relatively low for birds, which is expected for cooperative breeding species, while Ne/Nc ratios strongly varied among local populations. While the overall trend of increasing population sizes and genetic admixture may suggest that Cabanis’s greenbuls increasingly cope with fragmentation, the time period over which these trends were documented is rather short relative to the average longevity of tropical species. Furthermore, the critically low Nc in the small forest remnants keep the species prone to demographic and environmental stochasticity, and it remains open if, and to what extent, its cooperative breeding behaviour helps to buffer such effects. Keywords: approximate Bayesian computation, capture–recapture, census population size, demography, effective population size, gene flow, MARK, STRUCTURE, Taita Hills Received 1 December 2014; revision received 29 January 2015; accepted 4 February 2015

Introduction Population size and genetic connectivity largely determine the vulnerability of populations to demographic and genetic stochasticity (Lande 1988; Ellstrand & Elam Correspondence: Martin Husemann, Fax: +49-345-5527428; E-mail: [email protected] 1 MH and LC contributed equally.

1993; Frankham 1995; Young & Clarke 2000; Amos & Balmford 2001) and are therefore at the core of most conservation biological studies (Richards et al. 2003). The effective size of a population (Ne; Wright 1931), in particular, constitutes a key parameter as it estimates the rate at which genetic variance erodes due to stochastic events in small populations (Crow & Kimura 1970) and influences the rate of evolution (Lanfear et al.

© 2015 John Wiley & Sons Ltd

F R A G M E N T A T I O N G E N E T I C S I N C A B A N I S ’ S G R E E N B U L 1173 2014). If each individual of a population contributes equally to the gene pool of the next generation (ideal population sensu Whitlock & McCauley 1999), then Ne equals the actual number of individuals in the population (census size Nc). However, natural populations almost invariably fluctuate in (census) size (Frankham 1995; Vucetich et al. 1997), are biased towards one sex (Frankham 1995) and/or show high variance in reproductive success (Nunney 1996; Storz et al. 2001), and as such, they do not adhere to idealized population settings. Ne estimates are therefore generally much smaller than Nc values (Caballero 1994; Leberg 2005). In birds, Ne/Nc ratios have typically been shown (or predicted) to range between 0.1 and 0.7 (Nunney 1993; Frankham 1995; Palstra & Ruzzante 2008; see also Ardren & Kapuscinski 2003; Watts et al. 2007; Palstra & Ruzzante 2008 for variation in Ne/Nc ratios within single species). Studies of Ne/Nc ratios have proven particularly informative when evaluating population effects of rapid habitat change, such as the progressing fragmentation of formerly continuous habitat blocks into small, isolated remnants. Habitat fragmentation can affect the genetic population structure of species both directly and indirectly, that is by restricting gene flow, reducing population sizes, and/or increasing genetic drift and inbreeding (reviewed in Frankham 1995; Frankham et al. 2002; Luikart et al. 2010). While few compelling conclusions can be drawn from Ne or Nc estimates separately, low Ne/Nc ratios are generally assumed to provide strong evidence for increased sensitivity to genetic stochasticity (Palstra & Ruzzante 2008; Palstra & Fraser 2012). Traditionally, Ne estimation assumed single isolated populations that received no immigrants over the study interval (reviewed in Luikart et al. 2010). However, even in fragmented populations, dispersal and gene flow are often not negligible and can substantially bias estimates of Ne if not properly being accounted for (Wang & Whitlock 2003; Luikart et al. 2010; Waples & England 2011). A variety of estimators are now available, based on population samples taken at either single or multiple time points that relax the assumption of strict population isolation (Palstra & Ruzzante 2008; Luikart et al. 2010 and references therein). Single-sample methods are based either on linkage disequilibrium of markers or on sibship analysis, whereas temporal genetic methods use the premise that temporal variance in neutral genetic allele frequencies is inversely proportional to Ne (e.g. Barker 2011 for summary). While originally developed for species with discrete generations, these methods have subsequently been extended to cases where generations overlap (Palstra & Ruzzante 2008; Luikart et al. 2010 and references therein). However, the evaluation of multiple methods is necessary as each algorithm estimates a different type of population © 2015 John Wiley & Sons Ltd

size (i.e. inbreeding effective size (most single-sample methods), variance effective size (most temporal methods), Barker 2011) and may perform differently depending on the population structure, gene flow, population size and the sampling effort (e.g. Waples 2005; Luikart et al. 2010; Barker 2011; Holleley et al. 2014). Single-sample methods are now considered to be equally, or even more, reliable than temporal estimators but may in some cases underestimate Ne (Luikart et al. 2010; Barker 2011; Holleley et al. 2014). They are believed to be relatively accurate provided that random, mixed-aged samples include a number of consecutive age classes, approximately equal to a generation length, and populations are not very large (Jorde & Ryman 2007; Luikart et al. 2010; Waples & Do 2010). Given the debated question on the best estimation procedure for Ne, integrating multiple methods may yield the best estimates to judge the conservation status of a species (Waples 2005). One factor that may shape the ratio of effective to census population size is the social system a species lives in (Frankham 1995). For instance, in cooperative breeders where one or more individuals assist a breeding pair with rearing their offspring, the ratio is thought to be only half compared to breeding without helpers (Komdeur & Deerenberg 1997). Despite decades of research, no overall consensus has been reached on the evolution and maintenance of facultative cooperative breeding (Hatchwell 2009; Feeney et al. 2013; Riehl 2013; McDonald 2014) and on the extent to which it may allow organisms to adaptively respond to large-scale environmental change (Walters et al. 2004; Blackmore et al. 2011). While cooperative species share common ecological and life history characteristics such as pronounced philopatry and habitat specialization that may render them more vulnerable to habitat fragmentation and isolation (Walters et al. 2004), cooperative breeding can also be considered as the ‘best-of-a-bad-job’ under constrained environmental conditions such as shortage of breeding vacancies or high costs of dispersal (Hatchwell & Komdeur 2000). As such, it could represent a strategy to dampen adverse effects of habitat fragmentation by which cooperative species, to some extent, may be more resistant than noncooperative species with similar life history traits (Walters et al. 2004). A long-term study of eight forest birds inhabiting the highly fragmented Taita cloud forest archipelago of SE Kenya (ongoing since 1994) revealed sensitivity to both forest degradation and isolation in the cooperatively breeding Cabanis’s greenbul (Phyllastrephus cabanisi). While the species is not (yet) considered at risk within its Central and East African distribution range (Least Concern, IUCN 2014), it was shown to be consistently absent in nearly half of the Taita forest remnants, despite apparent habitat suitability (Lens et al. 2002), and to

1174 M . H U S E M A N N E T A L . suffer a comparatively strong loss of genetic connectivity over the past decades (Callens et al. 2011). Mark–recapture data (Lens et al. 2002) and homing experiments (Aben et al. 2012, 2014) showed low to moderate rates of mobility and among-fragment dispersal, while genetic data (Vangestel et al. 2013) further showed that juveniles and males tend to exhibit a ‘stay-and-foray’ dispersal strategy while adult females follow a ‘depart-and-search’ strategy within large forest patches. However, such sexand age-related variation was no longer apparent at the between-fragments scale, probably reflecting high costs of dispersal through the inhospitable landscape (Vangestel et al. 2013). At the forest patch level, P. cabanisi showed strong sensitivity to forest degradation, based on a comparison of historic (measured on museum specimen; prefragmentation) and current (measured on field captures; post-fragmentation) levels of fluctuating asymmetry in tarsus length (Lens et al. 1999). Unlike sympatric forest birds in the study area, P. cabanisi forages and breeds in small family groups (see further), and group living species are thought to exceed the carrying capacity of small and disturbed habitat remnants more easily and pay higher costs of biotic interactions near habitat edges and in the landscape matrix than pair-living ones (Van Houtan et al. 2007). Despite the fact that P. cabanisi currently ranks among the best studied forest birds in sub-Saharan Africa, the meta-population dynamics of tropical cooperative breeders in heterogeneous landscapes remain poorly known; yet, such information is essential when aiming to evaluate their long-term viability under persisting habitat change. To assess the vulnerability to rapid genetic erosion and predict temporal changes in genetic meta-population structure under progressing habitat fragmentation, we here analyse capture–mark–recapture histories and genetic data from all five surviving P. cabanisi populations in the rapidly degrading Taita forest archipelago during 2–5-year periods separated by a 6-year time interval. We expect Ne/Nc ratios to be smaller than ratios from those commonly measured in pair-breeding species as fewer individuals breed. We further predict increased variation in the degree of genetic population structuring and of asymmetrical gene flow, that is with higher rates from large populations in high-quality plots to small populations in strongly degraded fragments, in accordance with an emerging source-sink population network (Gaggiotti & Smouse 1996; Gaggiotti 1996).

Materials and methods Study area and species The Taita Hills (SE Kenya, 03°200 S, 38°150 E) represent the northernmost edge of the Eastern Arc Mountains

(EAM) Biodiversity Hotspot, a chain of ancient mountains that runs from southeastern Kenya to southern Tanzania (Fig. 1, Lovett & Wasser 1993). Indigenous forests are currently scattered over ca. 5076 km² (Platts et al. 2010) and are strongly isolated from other East African mountain forests by vast expanses of low-altitude plains. Despite the fact that vast forest tracts have been lost or are severely degraded and fragmented due to anthropogenic activities (Newmark 1998; Stattersfield et al. 1998), the EAM forests still harbour exceptionally high levels of floral and faunal diversity and endemism (Mittermeier et al. 1998). Within the Taita section of the EAM, indigenous forest cover decreased by ca. 50% between 1955 and 2004 (Pellikka et al. 2009), mainly due to clearance for small subsistence agriculture (Bennun & Njoroge 1999; Myers et al. 2000). The remaining forest archipelago currently covers ca. 470 ha of indigenous forest, fragmented into three larger patches and nine tiny remnants that are located on two mountain isolates (Dabida and Mbololo) separated by a low-altitude valley. Phyllastrephus cabanisi is a medium-sized cooperative breeder that inhabits Central to East African moist forests (Keith et al. 1992; Bennun et al. 1996). It lives in small family groups (breeding pair and 0–3 helpers; D. Van de Loock, unpubl. data), occasionally feeds on fruits, but supplies nestlings with invertebrates and small vertebrates (Keith et al. 1992; Lehouck et al. 2009). While the species has been trapped or observed in all twelve Taita forest fragments, breeding is currently restricted to the three largest (Mbololo (MB) 185 ha; Ngangao (NG) 120 ha; Chawia (CH) 86 ha) and two small (Fururu (FU) 8 ha; Ndiwenyi (ND) 4 ha) patches, while immigrants have unsuccessfully attempted to breed in one additional forest remnant (Macha (MA) 2.5 ha). The Taita meta-population comprises a closed network that is strongly isolated from the nearest P. cabanisi populations in the Chyulu Hills (Kenya, ca. 100 km NW) and Pare Mountains (Tanzania, ca. 100 km SW).

Trapping, measurements and census population sizes Between 1996 and 2010, mist net lines were operated in multiple plots per fragment (depending on fragment size) and were evenly spaced out to sample entire plots. Net positions, net lengths (120 m/plot) and daily trapping efforts (06–18 h) were kept constant between trapping sessions. Time intervals between subsequent ringing sessions varied between 1.0 and 4.6 months, and the number of ringing sessions per fragment ranged between 20 and 32 over the 15-year study period. Sampling effort (expressed as total number of trapping days) equalled 638 during the first period and 883 © 2015 John Wiley & Sons Ltd

F R A G M E N T A T I O N G E N E T I C S I N C A B A N I S ’ S G R E E N B U L 1175

Fig. 1 Sampling map showing the location of the Taita Hills within Kenya; the inlet displays the sampling locations within the Taita Hills.

during the second one. Ageing trapped individuals following Jackson (2005) allowed us to accurately differentiate recently fledged juveniles from older (i.e. fully grown) individuals, which was required for estimating previous generation migration rates (see further). After trapped individuals were aged, marked and measured, they were released at their original site of capture. To allow spatio-temporal data analysis, capture records from the five breeding populations were divided into two equal time periods, that is 1996–2000 (hereafter referred to as ‘first period’; 364 unique individuals) and 2006–2010 (‘second period’; 661 unique individuals), separated by a 6-year interval (Table 1). To avoid pseudoreplication, individuals sampled during the first period and still alive during the second were included in the first period only. The age structure of sampled individuals did not differ between the first and second period. To estimate Nc values for each population, we built capture–recapture encounter histories for all individuals © 2015 John Wiley & Sons Ltd

trapped during the first or second period and applied POPAN models (Schwarz & Arnason 1996) implemented in program MARK 6.0 (White & Burnham 1999) to estimate population- and period-specific recapture rates (p), yearly apparent survival rates (Φ), probabilities of entry into the population per session (pent) and total numbers of individuals present (N). We evaluated time-dependent (~t), linear time-trend (~T) and constant (~1) formulations for Φ, p and pent and additionally evaluated the effects of sampling effort by introducing a linear dependence on the number of sampling days during each session. Models were run for all alternative combinations of parameter formulations and were ranked according to Akaike’s information criterion (AICc) values corrected for small sample size (Hurvich & Tsai 1991). For this study, we analysed model-averaged estimates of N based on AICc weights (i.e. accounting for model uncertainty) only and extrapolated the harmonic mean of each period to the total area of indigenous forest (Nc). Standard errors were

F R A G M E N T A T I O N G E N E T I C S I N C A B A N I S ’ S G R E E N B U L 1175

Fig. 1 Sampling map showing the location of the Taita Hills within Kenya; the inlet displays the sampling locations within the Taita Hills.

during the second one. Ageing trapped individuals following Jackson (2005) allowed us to accurately differentiate recently fledged juveniles from older (i.e. fully grown) individuals, which was required for estimating previous generation migration rates (see further). After trapped individuals were aged, marked and measured, they were released at their original site of capture. To allow spatio-temporal data analysis, capture records from the five breeding populations were divided into two equal time periods, that is 1996–2000 (hereafter referred to as ‘first period’; 364 unique individuals) and 2006–2010 (‘second period’; 661 unique individuals), separated by a 6-year interval (Table 1). To avoid pseudoreplication, individuals sampled during the first period and still alive during the second were included in the first period only. The age structure of sampled individuals did not differ between the first and second period. To estimate Nc values for each population, we built capture–recapture encounter histories for all individuals © 2015 John Wiley & Sons Ltd

trapped during the first or second period and applied POPAN models (Schwarz & Arnason 1996) implemented in program MARK 6.0 (White & Burnham 1999) to estimate population- and period-specific recapture rates (p), yearly apparent survival rates (Φ), probabilities of entry into the population per session (pent) and total numbers of individuals present (N). We evaluated time-dependent (~t), linear time-trend (~T) and constant (~1) formulations for Φ, p and pent and additionally evaluated the effects of sampling effort by introducing a linear dependence on the number of sampling days during each session. Models were run for all alternative combinations of parameter formulations and were ranked according to Akaike’s information criterion (AICc) values corrected for small sample size (Hurvich & Tsai 1991). For this study, we analysed model-averaged estimates of N based on AICc weights (i.e. accounting for model uncertainty) only and extrapolated the harmonic mean of each period to the total area of indigenous forest (Nc). Standard errors were

F R A G M E N T A T I O N G E N E T I C S I N C A B A N I S ’ S G R E E N B U L 1177 Waples & Do 2008). A random mating model and a critical allele frequency value of 0.05 were applied. A second point estimate was generated with ONESAMP 1.2 (Tallmon et al. 2004). Here we used the minimum and maximum population sizes of 2 and 500, respectively. Both estimators have been suggested to yield reliable results when 10 or more unlinked loci and over 30 individuals are available (Luikart et al. 2010), yet might underestimate Ne in the face of null alleles (Sved et al. 2013) or population substructure (Holleley et al. 2014). Due to the presence of null alleles and low polymorphism in one marker, however, we were only able to include 9 markers in the analysis, thereby potentially limiting its statistical power. Third, we used two temporal approaches, the moment-based method as implemented in NeEstimator and a method based on allele frequency shifts performed in TEMPOFS (Jorde & Ryman 2007). Both programs were run with sampling plan 2, and 95% confidence intervals were derived from jackknifing across all loci. While more refined temporal approaches have recently become available (see Barker 2011; Luikart et al. 2010 for review), our data violated some of the underlying assumptions. Furthermore, as the long generation times of tropical species with slow life histories typically hamper sufficient spacing between temporal sampling periods (Waples 2005; Johnstone et al. 2012), we based our final estimate for Ne on the harmonic mean of the single-sample point estimates only. As the latter provide an estimate of Ne in the parental generation (Waples 2005) and Nc estimates correspond to the current population, it was only possible to derive population-specific Ne/Nc ratios for the first period, that is by using Ne from the second period and Nc from the first one. An estimate of meta-population Ne was generated under the island model as described by Gomez-Uchida et al. (2013) using the harmonic means of the Ne estimates of all populations and a global FST estimate calculated with ARLEQUIN 3.5 (Excoffier & Lischer 2010) (see below).

(Meta)-population structure and demographic history We performed multiple analyses of molecular variance (AMOVA) in ARLEQUIN to partition the genetic variance among populations, between both sampling periods, and within periods. Next, we used ARLEQUIN to calculate global fixation indices and pairwise FST estimates among populations. Genetic population structure during the first and second period was inferred from Bayesian admixture models implemented in STRUCTURE 2.3.4 (Pritchard et al. 2000). STRUCTURE probabilistically assigns individuals into genetic clusters (K) depending on their multilocus genotypes that minimize deviation from Hardy–Weinberg equilibrium and linkage © 2015 John Wiley & Sons Ltd

disequilibrium within clusters. We used the admixture model with correlated allele frequencies allowing gene flow between populations, and K was restricted between 1 and 10. A total of 100 independent iterations (each 150 000 sweeps long and discarding the first 50 000 sweeps as burn-in) were run for each K. Following Evanno et al. (2005), we selected the ‘true’ number of K using the second-order rate of change in the log probability of the data between successive K values (DK) as modelled with STRUCTURE HARVESTER (Earl & von Holdt 2012). Cluster assignments across the 10 replicate runs for each K were aligned with CLUMPP 1.1.2 (Jakobsson & Rosenberg 2007), and results were visualized with DISTRUCT 1.1 (Rosenberg 2004). Meta-population dynamics were further assessed by estimating previous generation migration rates during each period with BIMr 1.0 (Faubet & Gaggiotti 2008). BIMr uses the multilocus genetic disequilibrium in migrants and their recent descendants to infer the proportion of immigrants in a given population, and as such, it relaxes the assumption of Hardy–Weinberg equilibrium (Faubet & Gaggiotti 2008). We ran 10 replicates to ensure convergence of MCMC. Each replicate started with 20 short pilot runs of 1000 iterations followed by 100 000 iterations discarded as burn-in. We then ran 5 000 000 iterations from which samples were drawn every 100 iterations for a total sample size of 50 000 samples for each replicate. We used the model with correlated allele frequencies (F-model) that takes the population admixture that may have taken place before the last generation of migration into account; this procedure is believed to improve migration estimates (Faubet & Gaggiotti 2008). We selected parameter estimates from the run with the lowest Dassign (i.e. the probability of a particular assignment given a migration rate) to facilitate convergence (Faubet et al. 2007). The 95% highest posterior density intervals (HPDI) were used to test significance of asymmetrical pairwise migration rates among the five populations and between both periods. Because BIMr assumes that sampling has taken place before migration, consistency of our results was checked by re-analysing DNA samples from 222 individually marked nestlings and recent fledglings that were still present in their natal population upon sampling (CH=98, NG=95, MB=10, FU=16, ND=3; samples only available for the second period). Finally, the effective number of migrants from population j into i during the previous generation was calculated by multiplying the immigration rate mij (i.e. probability that alleles in population i came from population j during the previous generation) with the effective size Ne of population i (Tables 1 and 2). Deviation from mutation–drift equilibrium was tested for each temporal sample (period) from each population

1178 M . H U S E M A N N E T A L . Table 2 Estimates of effective population size and the ratio of effective to census population size. Ne estimates were derived as the harmonic means of two single-sample estimators and two temporal approaches. Given are the point estimates and 95% confidence intervals. The ratios of Ne to Nc were calculated from the harmonic means of Ne from the second sampling period and the Nc estimates from the first period. Meta-Ne estimates for the first and second period of 172 and 218 individuals were estimated from the point estimates of Ne and the global FST as derived from ANOVA (meta-Nc for the first period and second period are 1680 and 1900); the global ratio of meta-Ne to meta-Nc equals 0.13

Population

NeLD (second period)

NeABC (second period)

Moment-based NeEstimator

TempoFS

Harmonic means (and ranges) of single-sample point estimates (second period)

Chawia Fururu Mbololo Ndiwenyi Ngangao

75.3 22.5 53.7 4.9 76.7

52.0 27.7 39.7 9.0 71.6

116.5 34.7 119.2 26.4 98.3

68 7 40 26 138

61.521 24.83 45.65 6.35 74.06

(61.2–95.1) (17.5–30.3) (39.5–79.0) (3.4–7.8) (61.5–98.7)

(28.9–105.6) (19.9–46.3) (28.3–69.4) (7.1–13.2) (48.6–154.1)

(53.8–348.0) (13.3–285.5) (43.7–3181.4) (7–∞) (46.9–268.7)

with BOTTLENECK 1.2.02 (Piry et al. 1999). This program generates the distribution of heterozygosity expected under mutation–drift equilibrium for each locus and population (or cluster, see higher), based on the assumption that recently bottlenecked populations exhibit a faster reduction in the number of alleles than in the level of heterozygosity, resulting in an excess of heterozygotes. Because the mutation model underlying the microsatellite markers is unknown, data were analysed under both the two-phase model (TPM) and stepwisemutation model (SMM) (Di Rienzo et al. 1994; Jarne & Lagoda 1996; Luikart & Cornuet 1998; Piry et al. 1999). When employing TPMs, combinations of 95% singlestep mutations and 5% multistep mutations were used, with a variance of 30 among multiple-step mutations (104 replications) (Piry et al. 1999). Expected values were compared with observed heterozygosity levels as calculated from observed allele frequencies (sensu Nei et al. 1975; Nei 1987), and Wilcoxon’s signed-rank tests were used to test for significant heterozygote excess. The demographic history of the three largest populations (CH, NG, MB) was further inferred through approximate Bayesian computation (ABC) analysis of posterior distributions of population divergence times and population size changes, implemented in DIYABC 2.0.4 (Cornuet et al. 2014). In this procedure, the likelihood distribution of the Bayes formula is replaced by a measure of similarity between observed and simulated data that are summarized by a set of summary statistics (Cornuet et al. 2008). As the precise generation time of P. cabanisi is not known, we compared two historical models differing in the assumed number of generations between both sampling periods (t0 and t1). We thereby assumed a generation time of 6 years (resulting in a between-sampling interval of one generation) in the first scenario and of 3 years (i.e. two generation sampling

(29–∞) (14–∞) (15–∞) (6–∞) (72–1780)

(28.9–105.6) (17.5–46.3) (28.3–79.0) (3.4–13.2) (48.6–154.1)

Nc (first period)

Ne/Nc

424 12 845 13 386

0.14 2.02 0.05 0.43 0.19

(332–516) (8–16) (641–1049) (7–19) (304–468)

interval) in the second scenario. Based on the current topography of the Taita Hills, we constrained both scenarios such that populations CH and NG diverged from an ancestral population A1 at time t2, while this ancestral population and population MB diverged from another ancestral population A2 backward in time at t2 + t3 (Fig. 2). As a prior for Ne, we used a narrow normal distribution centred at the harmonic mean of Ne estimates per period (Tables 1 and 2). All other priors are listed in Table 3. Statistics used to summarize the genetic variation within and between the three populations and two periods comprised the mean number of

N e2 t2+t3 N e1

t2 N eCh_t1

N eNg_t1

NeMb_t1

t1 (Period1) N eCh_t0

NeNg_t0

NeMb_t0 t0 (Period 2)

Chawia

Ngangao

Mbololo

Dabida

Fig. 2 Demographic scenarios used for approximate Bayesian computation analysis (scenario 1: one generation between t0 and t1; scenario 2: two generations between t0 and t1) with Ne1=NeCh_t1+NeNg_t1 and Ne2=Ne1+NeMb_t1. Prior and posterior distributions of the model parameters are shown in Table 4. © 2015 John Wiley & Sons Ltd

F R A G M E N T A T I O N G E N E T I C S I N C A B A N I S ’ S G R E E N B U L 1177 Waples & Do 2008). A random mating model and a critical allele frequency value of 0.05 were applied. A second point estimate was generated with ONESAMP 1.2 (Tallmon et al. 2004). Here we used the minimum and maximum population sizes of 2 and 500, respectively. Both estimators have been suggested to yield reliable results when 10 or more unlinked loci and over 30 individuals are available (Luikart et al. 2010), yet might underestimate Ne in the face of null alleles (Sved et al. 2013) or population substructure (Holleley et al. 2014). Due to the presence of null alleles and low polymorphism in one marker, however, we were only able to include 9 markers in the analysis, thereby potentially limiting its statistical power. Third, we used two temporal approaches, the moment-based method as implemented in NeEstimator and a method based on allele frequency shifts performed in TEMPOFS (Jorde & Ryman 2007). Both programs were run with sampling plan 2, and 95% confidence intervals were derived from jackknifing across all loci. While more refined temporal approaches have recently become available (see Barker 2011; Luikart et al. 2010 for review), our data violated some of the underlying assumptions. Furthermore, as the long generation times of tropical species with slow life histories typically hamper sufficient spacing between temporal sampling periods (Waples 2005; Johnstone et al. 2012), we based our final estimate for Ne on the harmonic mean of the single-sample point estimates only. As the latter provide an estimate of Ne in the parental generation (Waples 2005) and Nc estimates correspond to the current population, it was only possible to derive population-specific Ne/Nc ratios for the first period, that is by using Ne from the second period and Nc from the first one. An estimate of meta-population Ne was generated under the island model as described by Gomez-Uchida et al. (2013) using the harmonic means of the Ne estimates of all populations and a global FST estimate calculated with ARLEQUIN 3.5 (Excoffier & Lischer 2010) (see below).

(Meta)-population structure and demographic history We performed multiple analyses of molecular variance (AMOVA) in ARLEQUIN to partition the genetic variance among populations, between both sampling periods, and within periods. Next, we used ARLEQUIN to calculate global fixation indices and pairwise FST estimates among populations. Genetic population structure during the first and second period was inferred from Bayesian admixture models implemented in STRUCTURE 2.3.4 (Pritchard et al. 2000). STRUCTURE probabilistically assigns individuals into genetic clusters (K) depending on their multilocus genotypes that minimize deviation from Hardy–Weinberg equilibrium and linkage © 2015 John Wiley & Sons Ltd

disequilibrium within clusters. We used the admixture model with correlated allele frequencies allowing gene flow between populations, and K was restricted between 1 and 10. A total of 100 independent iterations (each 150 000 sweeps long and discarding the first 50 000 sweeps as burn-in) were run for each K. Following Evanno et al. (2005), we selected the ‘true’ number of K using the second-order rate of change in the log probability of the data between successive K values (DK) as modelled with STRUCTURE HARVESTER (Earl & von Holdt 2012). Cluster assignments across the 10 replicate runs for each K were aligned with CLUMPP 1.1.2 (Jakobsson & Rosenberg 2007), and results were visualized with DISTRUCT 1.1 (Rosenberg 2004). Meta-population dynamics were further assessed by estimating previous generation migration rates during each period with BIMr 1.0 (Faubet & Gaggiotti 2008). BIMr uses the multilocus genetic disequilibrium in migrants and their recent descendants to infer the proportion of immigrants in a given population, and as such, it relaxes the assumption of Hardy–Weinberg equilibrium (Faubet & Gaggiotti 2008). We ran 10 replicates to ensure convergence of MCMC. Each replicate started with 20 short pilot runs of 1000 iterations followed by 100 000 iterations discarded as burn-in. We then ran 5 000 000 iterations from which samples were drawn every 100 iterations for a total sample size of 50 000 samples for each replicate. We used the model with correlated allele frequencies (F-model) that takes the population admixture that may have taken place before the last generation of migration into account; this procedure is believed to improve migration estimates (Faubet & Gaggiotti 2008). We selected parameter estimates from the run with the lowest Dassign (i.e. the probability of a particular assignment given a migration rate) to facilitate convergence (Faubet et al. 2007). The 95% highest posterior density intervals (HPDI) were used to test significance of asymmetrical pairwise migration rates among the five populations and between both periods. Because BIMr assumes that sampling has taken place before migration, consistency of our results was checked by re-analysing DNA samples from 222 individually marked nestlings and recent fledglings that were still present in their natal population upon sampling (CH=98, NG=95, MB=10, FU=16, ND=3; samples only available for the second period). Finally, the effective number of migrants from population j into i during the previous generation was calculated by multiplying the immigration rate mij (i.e. probability that alleles in population i came from population j during the previous generation) with the effective size Ne of population i (Tables 1 and 2). Deviation from mutation–drift equilibrium was tested for each temporal sample (period) from each population

1180 M . H U S E M A N N E T A L . K = 3 (i.e. CH/MB/NG-FU-ND) for both periods (Fig. 3, Appendix S2, Supporting information). During the second period, the level of genetic admixture slightly increased for each value of K, most strongly so from population MB to CH and to cluster NG-FU-ND (Fig. 3). Overall, BIMr analyses yielded consistent migration estimates across all replicates while migration rates inferred from nestlings/fledglings concurred with those derived from adult genotypes (Fig. 4, Appendix S3, Supporting information). However, immigration rates in FU appeared to be incongruent as opposing scenarios emerged among replicates in period 1 and between adult and nestling/fledgling genotypes in period 2. In the first period, the most optimal model assigned zero immigrants in FU. Likewise, in period 2, no immigrants were detected for FU when considering nestling/fledgling genotypes only, whereas migration rates based on adult birds indicated 50% immigrants in FU. Overall, migration rates during the second period were highly congruent with those during the first period as indicated by the extensive overlap of both HPD intervals. Yet, there was weak support for a decrease in migration from FU to CH and a slight increase in migration from MB and ND to CH, and from all fragments to FU over time. None of the populations (or population clusters) showed evidence of heterozygosity excess in either period (Wilcoxon’s signed-rank test; TPM and SMM: all P > 0.05). However, during the second period, populations NG and FU showed significant heterozygosity deficiency for the stepwise-mutation model compared to levels expected under mutation–drift equilibrium (Wilcoxon’s signed-rank tests; two-tailed probabilities; NG: TPM P = 0.129; SMM P = 0.010; FU: TPM P = 0.129; SMM P = 0.010). This deficit was significant

Period 1

Table 4 Test statistics, variance components and percentage variation explained as obtained from an analysis of molecular variance with predefined populations and periods (see text for details)

Source of variation Among populations Among periods within populations Within periods Total

Sum of squares

Variance components

Percentage variation

97.907 21.774

0.10903 0.01427

3.55223 0.46496

2975.633 3095.314

2.94617 3.06948

95.98281

for both mutation models when FU, NG and ND were clustered (Wilcoxon’s signed-rank test; TPM P = 0.014; SMM P = 0.010). Based on our ABC analysis, scenario 2 (i.e. generation time of 3 years) was selected as the better of the two models (posterior probability = 1 for both methods of model comparison). Estimates of divergence times for this scenario indicated that populations CH and NG diverged 48 [15; 381] years ago (mode, 95% HPD), while their ancestral population diverged 90 [45; 1,386] years ago from population MB (Table 3, Fig. 2). As DIYABC does not allow for the implementation of gene flow in the model and our BIMr analyses do suggest the presence of gene flow, these estimates need to be treated with caution. Yet, considering the current lack of alternatives, they provide the best possible approximation.

Discussion Despite its protected status, the Eastern Afromontane biodiversity hotspot is under constant threat (Newmark 1998; Stattersfield et al. 1998). The Taita Hills that are home to several endemic species of birds and other

Period 2

K=2

K=3

K=4

Fig. 3 Structure plots for five P. cabanisi populations sampled during two non-overlapping periods (period 1: 1996–2000; period 2: 2006–2010) as inferred from Bayesian genetic clustering for K = 2, 3 and 4 (with K = 3 evaluated as the best model). Each bar represents an individual partitioned according to his admixture proportions from genetic clusters. © 2015 John Wiley & Sons Ltd

F R A G M E N T A T I O N G E N E T I C S I N C A B A N I S ’ S G R E E N B U L 1181

Period 1

Period 2 5

6

1

5

2

2

1

13

1

3

3

22

6 15

11

2 6

13

Fig. 4 Previous generation effective number of migrants Nei*mij for period 1 and period 2, with Nei the harmonic mean estimates from Table 1 and mij the BIMr estimates derived from adult genotypes (replicates of each period with the lowest Dassign) (Appendix S3, Supporting infomation).

taxa, in particular, have suffered under strong anthropogenic pressure (Pellikka et al. 2009). Severe forest fragmentation during the 1960s was earlier shown to have caused strong population fragmentation in several Taita cloud forest specialists, as inferred from spatial patterns of population persistence (Lens et al. 2002) and temporal shifts in population connectivity (Callens et al. 2011). Building on these results, we here present a detailed demographic and genetic study of all five remnant populations of P. cabanisi over a 15-year post-fragmentation interval. We found small effective population sizes potentially resulting in strong stochastic effects, dependence of the ratio of Ne and Nc on the size of the population, a very recent divergence of the larger populations following deforestation as revealed by ABC models, and a meta-population structure with source and sink subpopulations. The effective meta-population size was under the critical population size limit of 500– 5000 individuals often suggested for vertebrate taxa (Lynch & Lande 1998; Franklin & Frankham 1998; Traill © 2015 John Wiley & Sons Ltd

et al. 2007; Jamieson & Allendorf 2012; but see Frankham et al. 2014 for critical notes on these numbers), and some populations appear prone to local extinction due to outward migration that is not fully compensated by inward migration. While the overall Ne to Nc ratio of 0.13 for the meta-population size is relatively low in comparison with other bird species, our temporal data suggest that some populations may have slightly increasing, rather than decreasing, census and effective sizes, and both the census and effective metapopulation sizes increased too. Yet, the time period over which these trends were documented is still short relative to the average longevity of tropical species, and it remains to be seen whether the species will be able to thrive despite its strong population signatures of forest fragmentation and isolation. In particular, given the critically low census population sizes in the small forest remnants, it probably remains prone to demographic and environmental stochasticity and might therefore not be viable in the long term.

1182 M . H U S E M A N N E T A L .

Post-fragmentation population structure Comparing the genetic structure and connectivity between both periods suggests that populations remained relatively stable during the time frame of our study. Structure analyses detected three clusters for both periods and migration estimates changed only slightly. Bayesian analyses further suggested a very recent divergence of the greenbul populations starting about 90 years ago (~1920), with a first split between the Dabida fragment (CH, NG) and MB and a second split ~50 years ago (~1960) between forest fragments CH and NG. As this second recent divergence corresponds to a period of intense deforestation starting in the 1950s, indigenous forest fragmentation likely contributed to population divergence in this region. However, as stated above, these time estimates have to be treated with caution due to the limitation of the method in considering gene flow. Despite the significant differentiation into three genetic clusters, gene flow appears to be relatively high, subpopulations remain relatively well connected, and small forest fragments seem to act as essential stepping stones connecting populations of the larger forest patches within a meta-population network. Asymmetrical rates of gene flow between patches further suggest that source-sink dynamics are becoming established, with the population of the most intact forest fragment (MB) serving as main source for the other large populations. However, the apparent ongoing, or even increased, functional connectivity in this study does not necessarily warrant long-term population sustainability as multi-species studies have shown habitat fragmentation induced population declines despite unaltered gene flow rates (Ford 2011; Amos et al. 2014). In addition, the apparent high rates of outward migration in the small forest remnants may be an artefact of the relatively low samples sizes and therefore should be treated with caution. Moreover, Bayesian analyses (on which the migration estimates were based) assume drift-migration equilibrium that is rarely present in natural populations (Faubet & Gaggiotti 2008); hence, migration estimates need to be interpreted as relative (rather than absolute) measures.

Demographics and population sizes Despite Ne estimates as low as 4 (ND) and 7 (FU) individuals during 1996–2000 (the generation previous to this interval) that are well below the thresholds of 50 and 500 individuals generally assumed to minimize the effects of genetic inbreeding and drift (Lande 1988; Franklin & Frankham 1998; Lynch & Lande 1998), none of the populations showed evidence for a recent

population bottleneck (heterozygote excess) in eitherperiod, despite a potentially high power to detect it after few generations (Cornuet & Luikart 1996). In contrast, two populations (NG and FU) showed lower levels of heterozygosity than expected under an equilibrium SSM during the second period. While such a heterozygote deficit might result from inadequate sample sizes or the presence of null alleles (Cornuet & Luikart 1996), sample sizes were comparable among populations and there was no indication of null alleles (see MICROCHECKER results). Hence, we believe the observed deficit to have resulted from an influx of new alleles from genetically distinct populations (in population NG and FU; see above) as recent immigrants can quickly increase the number of rare alleles without substantially affecting the level of heterozygosity mimicking an increase in population size (Cornuet & Luikart 1996). This scenario was supported by the corresponding increase in Nc in NG and Ne in both populations from the first to the second period and by the presence of inward migration in NG and FU (Fig. 4). Apart from population MB, estimates of Ne increased in all populations during the second period (Table 1), probably as a result of increased gene flow among them. Counter-intuitively, predation rates on P. cabanisi eggs and nestlings were significantly lower, and reproductive success significantly higher, in poor quality fragments and closer to fragment edges (inverse edge effect), possibly due to the decrease in nest predators in these highly disturbed areas (Spanhove et al. 2013). Ne estimates, however, even when based on spatiotemporal genetic data and multiple statistical estimators, such as in this study, remain prone to sampling variation, in particular in populations with mixed age structure, overlapping generations, strong gene flow or when few generations separate temporal samples (Wang & Whitlock 2003; Waples & Yokota 2007; Palstra & Ruzzante 2008; Luikart et al. 2010; Tallmon et al. 2010; Sved et al. 2013; Holleley et al. 2014). In our study, temporal samples were only spaced out for one generation, which may cause Ne estimates to vary in diverse and complex ways (i.e. depending on patterns of population-specific survivorship and age classes sampled; Luikart et al. 2010). Therefore, we based our harmonic mean Ne estimates on single-sample estimators, despite the fact that adult survival rates did not differ between populations nor periods (Table 1), and equal random mixtures of all age classes were sampled during both time periods. Given the fact that we sampled all age classes, we believe that the basic assumptions of single-sample estimators were met. Furthermore, single-sample and temporal estimates of Ne yielded comparable estimates, suggesting that potential biases associated with small sample sizes were probably small. We therefore © 2015 John Wiley & Sons Ltd

F R A G M E N T A T I O N G E N E T I C S I N C A B A N I S ’ S G R E E N B U L 1183 consider our Ne estimates fairly reliable, despite the relatively low number of genetic markers applied that inevitably reduced the analytical power. Similar to the increase in population size in most populations, both estimates of meta-population size (effective and census) increased during the second period as well. Together with the increased levels of gene flow, this may indicate that the species is slowly adapting to – or coping with – the fragmentation of its habitat and/or that connectivity has recently improved. Data of sympatric bird species from the same forest archipelago earlier showed that the effects of fragmentation strongly depend on the dispersal ability of a species (Lens et al. 2002), allowing species with better dispersal abilities to maintain reasonably high levels of genetic diversity despite population isolation (e.g. Callens et al. 2011). Therefore, despite the potentially increasing population size and gene flow, the low dispersal ability of P. cabanisi (Lens et al. 2002) is in agreement with the still low effective meta-population size (218 individuals) estimated in this study, which is far below the lower limit of 500–1000 individuals considered to buffer negative effects of drift and inbreeding (Franklin & Frankham 1998; Jamieson & Allendorf 2012).

Ne to Nc ratio Based on Komdeur & Deerenberg (1997) and references therein, we hypothesized low Ne/Nc ratios in P. cabanisi relative to other ecologically related species, due to its cooperative breeding behaviour. The global Ne/Nc estimate of 0.13 falls within the range of 0.10–0.14 reported for the majority of vertebrate taxa (Frankham 1995; Palstra & Ruzzante 2008), but is at the low end of what has been reported for birds in the same studies. For instance, in a range of bird species encompassing Darwin’s finches, owls, hawks and cranes, most estimates were above 0.4 (Frankham 1995 and references therein), while in common species such as house sparrows (Passer domesticus) and great tits (Parus major), estimates were even above 0.6 (Nunney & Elam 1994) . In contrast, rare and endangered species generally show much lower ratios, often below 0.1 (e.g. 0.03 in little spotted kiwi Apteryx owenii, Ramstad et al. 2013; 0.09 in jabiru stork Jabiru mycteria, Lopes et al. 2013). Given this wide range in avian Ne/Nc ranges, we cannot determine whether the low ratio in P. cabanisi reflects its cooperative mating system in the study area, other mechanisms related to its threat status, or a combination of both. Ne/Nc ratios of single populations strongly varied, that is between 0.05 (MB) and 2.02 (FU). Such strong variation, and such high ratios in the small fragments in particular, may have several nonexclusive reasons: (i) in small populations, the strength of sexual selection may decrease and more individuals may contribute to the © 2015 John Wiley & Sons Ltd

gene pool (Price 1998), hence leading to smaller group sizes and increasing effective population sizes. While P. cabanisi group sizes did not consistently vary among fragments in the Taita archipelago (L. Lens unpubl. data), possible differences in social group structure (e.g. level of extra-pair or extra-group matings) cannot be ruled out; (ii) effective population sizes may be elevated in small sink populations through new genetic diversity introduced by immigrants, while census population sizes reflect the carrying capacity of the habitat; and (iii) in very small populations, population sizes are prone to strong stochasticity and fluctuations, and high Ne/Nc ratios may therefore partially reflect sampling artefacts. Irrespective of the underlying mechanism involved, the high Ne/Nc ratios observed in smaller fragments, together with the inverse relationship between Nc and the standardized variance in family size commonly observed in small populations (Hedrick 2005), are assumed to compensate genetic effects typically associated with low Nc (reviewed in Palstra & Ruzzante 2008). In salmonid fishes, for example, such ‘genetic compensation’ mechanism was associated with particular breeding properties that increase fertilization success at low breeding densities (e.g. Jones & Hutchings 2001, 2002). Whether, and to what extent, such mechanisms may also operate in fragmented populations of cooperative breeders, such as P. cabanisi, for instance, through reduced numbers of nonbreeding helpers under low K values remains a topic of further research. Spatiotemporal analysis of genetic and demographic data over a recent 15-year time span did not provide evidence for increased genetic erosion and isolation in fragmented populations of P. cabanisi, but rather showed slightly increased levels of genetic admixture and effective population sizes. It is currently unclear which ecological or evolutionary factors may have triggered such population responses. As part of a series of conservation initiatives, over 150 000 indigenous tree seedlings have been raised in community-owned tree nurseries and planted within indigenous forest remnants and the Taita landscape matrix since 2006 (Githiru et al. 2011). While it is tempting to believe that these actions underlie the observed changes in genetic population structure in P. cabanisi, this is unlikely given the short time frame involved and a possible genetic lag phase before changes would manifest in the population genetic structure. Also, given the very small census population sizes in the small forest remnants, this species (and ecologically equivalent ones) probably remains prone to demographic and environmental stochasticity. Further study is therefore required to examine whether the particular mating system of cooperative breeding species increases or decreases their vulnerability against stochastic processes, and to what extent strong selective

1184 M . H U S E M A N N E T A L . pressures may ultimately change evolutionary trajectories in fragmented populations.

Acknowledgements We thank T. Spanhove, V. Lehouck, M. Chovu, S. Karimi, T. Brooks, D. Gitau, T. Imboma, J. Kageche and P. Kariuki for field assistance, T. Schenck, A. Van Vlaslaer, B. Bytebier and M. Githiru for logistic laboratory and field support and coordination and R.C.K. Bowie for sharing unpublished (at the time of testing) sequences of the primer set Pfi04 and Pfl54. Fieldwork and genetic analyses were funded by research grants G.0149.09 and 3G030813 of the Research Foundation Flanders (FWO) and through contacts facilitated by FWO research community WO.037.10.

References Aben J, Adriaensen F, Thijs K et al. (2012) Effects of matrix composition and configuration on forest bird movements in a fragmented Afromontane biodiversity hotspot. Animal Conservation, 15, 658–668. Aben J, Strubbe D, Adriaensen F et al. (2014) Simple individual-based models effectively represent Afrotropical forest bird movement in complex landscapes: implications for improved management recommendations. Journal of Applied Ecology, 51, 693–702. Amos W, Balmford A (2001) When does conservation genetics matter? Heredity, 87, 257–265. Amos N, Harrisson KA, Radford JQ et al. (2014) Species- and sex-specific connectivity effects of habitat fragmentation in a suite of woodland birds. Ecology, 95, 1556–1568. Ardren WR, Kapuscinski AR (2003) Demographic and genetic estimates of effective population size (N-e) reveals genetic compensation in steelhead trout. Molecular Ecology, 12, 35–49. Barker JSF (2011) Effective population size of natural populations of Drosophila buzzatii, with a comparative evaluation of nine methods of estimation. Molecular Ecology, 20, 4452–4471. Bennun L, Njoroge P (1999) Important Bird Areas in Kenya. EANHS, Nairobi. Bennun L, Dranzoa C, Pomeroy D (1996) The forest birds of Kenya and Uganda. Journal of the East Africa Natural History Society, 85, 23–48. Blackmore CJ, Peakall R, Heinsohn R (2011) The absence of sex-biased dispersal in the cooperatively breeding greycrowned babbler. Journal of Animal Ecology, 80, 69–78. Caballero A (1994) Developments in the prediction of effective population size. Heredity, 73, 657–679. Callens T, Galbusera P, Matthysen E et al. (2011) Genetic signature of population fragmentation varies with mobility in seven bird species of a fragmented Kenyan cloud forest. Molecular Ecology, 20, 1829–1844. Cornuet JM, Luikart G (1996) Description and power analysis of two tests for detecting recent population bottlenecks from allele frequency data. Genetics, 144, 2001–2014. Cornuet JM, Santos F, Beaumont MA et al. (2008) Inferring population history with DIY ABC: a userfriendly approach to approximate Bayesian computation. Bioinformatics, 24, 2713–2719.

Cornuet JM, Pudlo P, Veyssier J et al. (2014) DIYABC v2.0: a software to make approximate Bayesian computation inferences about population history using single nucleotide polymorphism. DNA sequence and microsatellite data. Bioinformatics, 30, 1187–1189. Crow JF, Kimura M (1970) An Introduction to Population Genetics Theory. Harper & Row, New York. Di Rienzo A, Peterson AC, Garza JC, Valdes AM, Slatkin M, Freimer NB (1994) Mutational processes of simple-sequence repeat loci in human-populations. Proceedings of the National Academy of Sciences of the United States of America, 91, 3166–3170. Do C, Waples RS, Peel D, Macbeth GM, Tillett BJ, Ovenden JR (2013) NEESTIMATOR V2: reimplementation of software for the estimation of contemporary effective population size (Ne) from genetic data. Molecular Ecology Resources, 1, 209– 214. Earl DA, von Holdt BM (2012) STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conservation Genetics Resources, 4, 359–361. Ellstrand NC, Elam DR (1993) Population genetic consequences of small population size - implications for plant conservation. Annual Review of Ecology and Systematics, 24, 217–242. England PR, Cornuet JM, Berthier P, Tallmon DA, Luikart G (2006) Estimating effective population size from linkage disequilibrium: Severe bias in small samples. Conservation Genetics, 7, 303–308. Evanno G, Regnaut S, Goudet J (2005) Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Molecular Ecology, 14, 2611–2620. Excoffier L, Lischer HE (2010) ARLEQUIN suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Molecular Ecology Resources, 10, 564–567. Excoffier L, Laval G, Schneider S (2005) ARLEQUIN (version 3.0): an integrated software package for population genetics analysis. Evolutionary Bioinformatics Online, 2005, 1. Faubet P, Gaggiotti OE (2008) A new Bayesian method to identify the environmental factors that influence recent migration. Genetics, 164, 1567–1587. Faubet P, Waples RS, Gaggiotti OE (2007) Evaluating the performance of a multilocus Bayesian method for the estimation of migration rates. Molecular Ecology, 16, 1149–1166. Feeney WE, Medina I, Somveille M et al. (2013) Brood parasitism and the evolution of cooperative breeding in birds. Science, 342, 1506–1508. Ford HA (2011) The causes of decline of birds of eucalypt woodlands: advances in our knowledge over the last 10 years. EMU, 111, 1–9. Frankham R (1995) Effective population-size adult-population size ratios in wildlife - a review. Genetical Research, 66, 95– 107. Frankham R, Ballou JD, Briscoe DA (2002) Introduction to Conservation Genetics. Cambridge University Press, Cambridge. Frankham R, Bradshaw CJA, Brook BW (2014) Genetics in conservation management: revised recommendations for the 50/ 500 rules, Red List criteria and population viability analyses. Biological Conservation, 170, 56–63. Franklin IR, Frankham R (1998) How large must populations be to retain evolutionary potential? Animal Conservation, 1, 69–70.

© 2015 John Wiley & Sons Ltd

F R A G M E N T A T I O N G E N E T I C S I N C A B A N I S ’ S G R E E N B U L 1185 Gaggiotti OE, Smouse PE (1996) Stochastic migration and maintenance of genetic variation in sink populations. The American Naturalist, 147, 919–945. Gaggiotti OE (1996) Population genetic models of source sink metapopulations. Theoretical Population Biology, 50, 178–208. Githiru M, Lens L, Adriaensen F, Mwang’ombe J, Matthysen E (2011) Using science to guide conservation: from landscape modeling to increased connectivity in the Taita Hills, SE Kenya. Journal of Nature Conservation, 19, 263–268. Goldstein DB, Ruiz Linares A, Cavalli-Sforza LL, Feldman MW (1995) An evaluation of genetic distances for use with microsatellite loci. Genetics, 139, 463–471. Gomez-Uchida D, Palstra FP, Knight TW, Ruzzante DE (2013) Contemporary effective population and metapopulation size (Ne and meta-Ne): comparison among three salmonids inhabiting a fragmented system and differing in gene flow and its asymmetries. Ecology and Evolution, 3, 569–580. Goudet J (1995) FSTAT (version 1.2): a computer program to calculate F-statistics. Journal of Heredity, 86, 485–486. Hatchwell BJ (2009) The evolution of cooperative breeding in birds: kinship, dispersal and life history. Philosophical Transactions of the Royal Society B-Biological Sciences, 364, 3217–3227. Hatchwell BJ, Komdeur J (2000) Ecological constraints, life history traits and the evolution of cooperative breeding. Animal Behaviour, 59, 1079–1086. Hedrick PW (2005) Large variance in reproductive success and the Ne/N ratio. Evolution, 59, 1596–1599. Holleley CE, Nichols RA, Whitehead MR, Adamack AT, Gunn MR, Sherwin WB (2014) Testing single-sample estimators of effective population size in genetically structured populations. Conservation Genetics, 15, 23–35. Hurvich CM, Tsai CL (1991) Bias of the corrected AIC criterion for underfitted regression and time-series models. Biometrika, 78, 499–509. IUCN (2014) The IUCN Red List of Threatened Species. Version 2014.3. Downloaded on 28 November 2014. Jackson C (2005) Ageing Afrotropical birds in the hand: a revised new system. Afring News, 34, 60–67. Jakobsson M, Rosenberg NA (2007) CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics, 23, 1801–1806. Jamieson IG, Allendorf FW (2012) How does the 50/500 rule apply to MVPs? Trends in Ecology and Evolution, 27, 578–584. Jarne P, Lagoda PJL (1996) Microsatellites, from molecules to populations and back. Trends in Ecology & Evolution, 11, 424– 429. Johnstone DL, O’Connell MF, Palstra FP, Ruzzante DE (2012) Mature male pall contribution to the effective size of an anadromous Atlantic salmon (Salmo salar) population over 30 years. Molecular Ecology, 22, 2394–2407. Jones MW, Hutchings JA (2001) The influence of male parr body size and mate competition on fertilization success and effective population size in Atlantic salmon. Heredity, 86, 675–684. Jones MW, Hutchings JA (2002) Individual variation in Atlantic salmon fertilization success: implications for effective population size. Ecological Applications, 12, 184–193. Jorde PE, Ryman N (2007) Unbiased estimator for genetic drift and effective population size. Genetics, 177, 927–935.

© 2015 John Wiley & Sons Ltd

Keith S, Urban EK, Fry CH (1992) The birds of Africa. Vol. 4: Broadbills to Chats, Academic Press, London. Komdeur J, Deerenberg C (1997) The importance of social behavior studies for conservation. In: Behavioral Approaches to Conservation in the Wild (eds Clemmons JR, Buchholz R), pp. 262–276. Cambridge University Press, Cambridge, UK. Lande R (1988) Genetics and demography in biological conservation. Science, 241, 1455–1460. Lanfear R, Kokko H, Eyre-Walker A (2014) Population size and the rate of evolution. Trends in Ecology & Evolution, 29, 33–41. Leberg P (2005) Genetic approaches for estimating the effective size of populations. Journal of Wildlife Management, 69, 1385– 1399. Lehouck V, Spanhove T, Vangestel C, Cordeiro NJ, Lens L (2009) Does landscape structure affect resource tracking by avian frugivores in a fragmented Afrotropical forest? Ecography, 32, 789–799. Lens L, Van Dongen S, Wilder CM, Brooks TM, Matthysen E (1999) Fluctuating asymmetry increases with habitat disturbance in seven bird species of a fragmented afrotropical forest. Proceedings of the Royal Society B - Biological Sciences, 266, 1241–1246. Lens L, Van Dongen S, Norris K, Githiru M, Matthysen E (2002) Avian persistence in fragmented rainforest. Science, 298, 1236–1238. Lopes IF, Mino CI, Rocha CD, Oliveira DMM, Del Lama SN (2013) Inferred kinship patterns reveal low levels of extrapair paternity in the endangered Neotropical Jabiru Stork (Jabiru mycteria, Aves: Ciconiiformes). Genetica, 141, 195–203. Lovett JC, Wasser SK (1993) Biogeography and Ecology of the Rain Forests of Eastern Africa. Cambridge University Press, Cambridge. Luikart G, Cornuet JM (1998) Empirical evaluation of a test for identifying recently bottlenecked populations from allele frequency data. Conservation Biology, 12, 228–237. Luikart G, Ryman N, Tallmon DA, Schwartz MK, Allendorf FW (2010) Estimation of census and effective population sizes: the increasing usefulness of DNA-based approaches. Conservation Genetics, 11, 355–373. Lynch M, Lande R (1998) The critical effective size for a genetically secure population. Animal Conservation, 1, 70–72. McDonald PG (2014) Cooperative breeding beyond kinship: why else do helpers help? EMU, 114, 91–96. Mittermeier RA, Myers N, Thomsen JB, da Fonseca GAB, Olivieri S (1998) Biodiversity hotspots and major tropical wilderness areas: approaches to setting conservation priorities. Conservation Biology, 12, 516–520. Myers N, Mittermeier RA, Mittermeier CG, da Fonseca GAB, Kent J (2000) Biodiversity hotspots for conservation priorities. Nature, 403, 853–858. Nei M (1987) Molecular Evolutionary Genetics. Columbia University Press, New York. Nei M, Maruyama T, Chakraborty R (1975) Bottleneck effect and genetic variability in populations. Evolution, 29, 1–10. Newmark WD (1998) Forest area, fragmentation, and loss in the Eastern Arc Mountains: implications for the conservation of biological diversity. Journal of East African Natural History, 87, 29–36. Nunney L (1993) The influence of mating system and overlapping generations on effective population size. Evolution, 47, 1329–1341.

1186 M . H U S E M A N N E T A L . Nunney L (1996) The influence of variation in female fecundity on effective population size. Biological Journal of the Linnean Society, 59, 411–425. Nunney L, Elam DR (1994) Estimating the effective populationsize of conserved populations. Conservation Biology, 8, 175–184. Palstra FP, Fraser DJ (2012) Effective/census population size ratio estimation: a compendium and appraisal. Ecology and Evolution, 2, 2357–2365. Palstra FP, Ruzzante DE (2008) Genetic estimates of contemporary effective population size: what can they tell us about the importance of genetic stochasticity for wild population persistence? Molecular Ecology, 17, 3428–3447. Pellikka PKE, Lotjonen M, Sijander M, Lens L (2009) Airborne remote sensing of spatiotemporal change (1955-2004) in indigenous and exotic forest cover in the Taita Hills, Kenya. International Journal of Applied Earth Observation and Geoinformation, 11, 221–232. Piry S, Luikart G, Cornuet JM (1999) BOTTLENECK: a computer program for detecting recent reductions in the effective population size using allele frequency data. Journal of Heredity, 90, 502–503. Platts PJ, Ahrends A, Gereau RE et al. (2010) Can distribution models help refine inventory-based estimates of conservation priority? A case study in the eastern arc forests of Tanzania and Kenya. Diversity and Distributions, 16, 628–642. Plummer M, Best N, Cowles K, Vines K (2006) CODA: Convergence Diagnosis and Output Analysis for MCMC. R News, 6, 7–11. Price T (1998) Sexual selection and natural selection in bird speciation. Philosophical Transactions of the Royal Society of London, B, 353, 251–260. Pritchard JK, Stephens M, Donnelly P (2000) Inference of population structure using multilocus genotype data. Genetics, 155, 945–959. R Core Team (2014) R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. URL http://www.R-project.org/. Ramstad KM, Colbourne RM, Robertson HA, Allendorf FW, Daugherty CH (2013) Genetic consequences of a century of protection: serial founder events and survival of the little spotted kiwi (Apteryx owenii). Proceedings of the Royal Society B, 280, 20130576. Raymond M, Rousset F (1995) GENEPOP (version-1.2) - Population genetics software for exact tests and ecumenicism. Journal of Heredity, 86, 248–249. Richards CM, Emery SN, McCauley DE (2003) Genetic and demographic dynamics of small populations of Silene latifolia. Heredity, 90, 181–186. Riehl C (2013) Evolutionary routes to non-kin cooperative breeding in birds. Proceedings of the Royal Society of London Series B-Biological Sciences, 280, 20132245. Rosenberg NA (2004) DISTRUCT: a program for the graphical display of population structure. Molecular Ecology Notes, 4, 137–138. Rousset F (2008) GENEPOP ‘007: a complete re-implementation of the genepop software for Windows and Linux. Molecular Ecology Resources, 8, 103–106. Schwarz CJ, Arnason AN (1996) A general methodology for the analysis of capture-recapture experiments in open populations. Biometrics, 52, 860–873. Spanhove T, Callens T, Jallmann CA, Pellikka P, Lens L (2013) Nest predation in Afrotropical forest fragments shaped by inverse edge effects, timing of nest initiation and vegetation structure. Journal of Ornithology, 155, 411–420.

Stattersfield AJ, Crosby MJ, Long AJ, Wege DC (1998) Endemic Bird Areas of the World: Priorities for Biodiversity Conservation. BirdLife International, Cambridge. Storz JF, Bhat HR, Kunz TH (2001) Genetic consequences of polygyny and social structure in an Indian fruit bat, Cynopterus sphinx. Ii. Variance in male mating success and effective population size. Evolution, 55, 1224–1232. Sved JA, Cameron EC, Gilchrist AS (2013) Estimating effective population size from linkage disequilibrium between unlinked loci: theory and application to fruit fly outbreak populations. PLoS ONE, 8, e69078. Tallmon DA, Luikart G, Beaumont MA (2004) Comparative evaluation of a new effective population size estimator based on approximate Bayesian computation. Genetics, 167, 977–988. Tallmon DA, Gregovich D, Waples RS et al. (2010) When are genetic methods useful for estimating contemporary abundance and detecting population trends? Molecular Ecology Resources, 10, 684–692. Traill LW, Bradshaw CJA, Brook BW (2007) Minimum viable population size: a meta-analysis of 30 years of published estimates. Biological Conservation, 139, 159–166. Van Houtan KS, Pimm SL, Halley JM, Bierregaard RO, Lovejoy TE (2007) Dispersal of Amazonian birds in continuous and fragmented forest. Ecology Letters, 10, 219–229. Van Oosterhout C, Weetman D, Hutchinson WF (2006) Estimation and adjustment of microsatellite null alleles in nonequilibrium populations. Molecular Ecology Notes, 6, 255–256. Vangestel C, Callens T, Vandomme V, Lens L (2013) Sex-biased dispersal at different geographical scales in a cooperative breeder from fragmented rainforest. PLoS ONE, 8, e71624. Vucetich JA, Waite TA, Nunney L (1997) Fluctuating population size and the ratio of effective to census population size. Evolution, 51, 2017–2021. Walsh PS, Metzger DA, Higuchi R (1991) Chelex-100 as a medium for simple extraction of DNA for PCR-based typing from forensic material. BioTechniques, 10, 506–513. Walters JR, Cooper CB, Daniels SJ, Pasinelli G, Schieg K (2004) Conservation biology. In: Ecology and Evolution of Cooperative Breeding in Birds (eds Koenig WD, Dickinson JL),pp. 197–209. Cambridge University Press, Cambridge, UK. Wang JL, Whitlock MC (2003) Estimating effective population size and migration rates from genetic samples over space and time. Genetics, 163, 429–446. Waples RS (2005) Genetic estimates of contemporary effective population size: to what time periods do estimates apply? Molecular Ecology, 14, 3335–3352. Waples RS, Do C (2008) LDNe: A program for estimating effective population size from data on linkage disequilibrium. Molecular Ecology Resources, 8, 753–756. Waples RS, Do C (2010) Linkage disequilibrium estimates of contemporary Ne using highly variable genetic markers: a largely untapped resource for applied conservation and evolution. Evolutionary Applications, 3, 244–262. Waples RS, England PR (2011) Estimating contemporary effective population size on the basis of linkage disequilibrium in the face of migration. Genetics, 189, 633–644. Waples RS, Yokota M (2007) Temporal estimates of effective population size in species with overlapping generations. Genetics, 175, 219–233.

© 2015 John Wiley & Sons Ltd

F R A G M E N T A T I O N G E N E T I C S I N C A B A N I S ’ S G R E E N B U L 1187 Watts PC, Saccheri IJ, Kemp SJ, Thompson DJ (2007) Effective population sizes and migration rates in fragmented populations of an endangered insect (Coenagrion mercuriale: Odonata). Journal of Animal Ecology, 76, 790–800. Weir BS, Cockerham CC (1984) Estimating F-Statistics for the analysis of population-structure. Evolution, 38, 1358–1370. White GC, Burnham KP (1999) Program MARK: survival estimation from populations of marked animals. Bird Study, 46, 120–139. Whitlock MC, McCauley DE (1999) Indirect measures of gene flow and migration: FST not equal 1/(4Nm+1). Heredity, 82, 117–125. Wright S (1931) Evolution in Mendelian populations. Genetics, 16, 97–159. Young AG, Clarke GM (2000) Genetics, Demography and Viability of Fragmented Populations. Cambridge University Press, Cambridge.

T.C. and L.L. collected the data. M.H., L.C. and T.C. analysed the data. M.H., L.C., T.C. and L.L. wrote a first draft of the manuscript and lead writing. All authors read and approved the final version of the manuscript.

© 2015 John Wiley & Sons Ltd

Data accessibility The microsatellite genotype data and the capture– mark–recapture data are available in Dryad: doi:10.5061/dryad.8n5d7.

Supporting information Additional supporting information may be found in the online version of this article. Appendix S1. PCR conditions and DNA specifications of 10 microsatellite markers for P. cabanisi. Appendix S2. Inference of the optimal number of genetic clusters (K) of Cabanis’s greenbul in the Taita Hills for two nonoverlapping time periods (Period 1: 1996-2000; Period 2: 20062010) based on the STRUCTURE algorithm. Appendix S3. Migration rates between five Cabanis’s greenbul (sub-)populations in the Taita Hills.