Mitochondrial control region I and microsatellite analyses of endangered Philippine hornbill species (Aves; Bucerotidae) detect gene flow between island populations and genetic diversity loss

Background The Visayan Tarictic Hornbill (Penelopides panini) and the Walden’s Hornbill (Aceros waldeni) are two threatened hornbill species endemic to the western islands of the Visayas that constitute - between Luzon and Mindanao - the central island group of the Philippine archipelago. In order to evaluate their genetic diversity and to support efforts towards their conservation, we analyzed genetic variation in ~ 600 base pairs (bp) of the mitochondrial control region I and at 12–19 nuclear microsatellite loci. The sampling covered extant populations, still occurring only on two islands (P. panini: Panay and Negros, A. waldeni: only Panay), and it was augmented with museum specimens of extinct populations from neighboring islands. For comparison, their less endangered (= more abundant) sister taxa, the Luzon Tarictic Hornbill (P. manillae) from the Luzon and Polillo Islands and the Writhed Hornbill (A. leucocephalus) from Mindanao Island, were also included in the study. We reconstructed the population history of the two Penelopides species and assessed the genetic population structure of the remaining wild populations in all four species. Results Mitochondrial and nuclear data concordantly show a clear genetic separation according to the island of origin in both Penelopides species, but also unravel sporadic over-water movements between islands. We found evidence that deforestation in the last century influenced these migratory events. Both classes of markers and the comparison to museum specimens reveal a genetic diversity loss in both Visayan hornbill species, P. panini and A. waldeni, as compared to their more abundant relatives. This might have been caused by local extinction of genetically differentiated populations together with the dramatic decline in the abundance of the extant populations. Conclusions We demonstrated a loss in genetic diversity of P. panini and A. waldeni as compared to their sister taxa P. manillae and A. leucocephalus. Because of the low potential for gene flow and population exchange across islands, saving of the remaining birds of almost extinct local populations - be it in the wild or in captivity - is particularly important to preserve the species’ genetic potential.


Background
The Philippine archipelago consists of more than 7,000 islands, mostly of volcanic origin, with a well known geological and plate tectonic history [1]. According to various authors e.g. [2][3][4][5][6], during the Pleistocene glaciations, the Philippine islands were less fragmented than today; groups of islands formed composite Pleistocene Aggregate Island Complexes (PAICs) [5], that were isolated from one another by deep-water channels. Apart from Greater Palawan, which, biogeographically, does not belong to the oceanic region of the Philippines [4], there were three major PAICs: Greater Luzon, Greater Negros-Panay, and Greater Mindanao [3]. The geological history has had a great impact on the evolution of the Philippines' fauna and flora [7] resulting in an extremely high rate of over 57% endemism in the major faunal groups [8], with many lineages limited to single islands or PAICs [4]. Among Philippine vertebrates, this is especially true for reptiles, amphibians, and mammals. In contrast, birds show, due to their volant nature, relatively low levels of endemism. Here, dispersal may be thus more an important driver of species diversity than vicariance [9].
Because of their species richness, the Philippines belong to the 25 biodiversity hotspots on earth, but are simultaneously the one with the lowest percentage (3%) of remaining primary vegetation [10]. According to satellite data from 1987, forest cover has declined to 23.8%, but varies considerably across the archipelago [11]. For example, Luzon retained 24% forest cover and Mindanao 29%; in the Western Visayas, Panay retained 8% and Negros 4%. Even these low numbers are probably overestimates, as only a proportion comprises closedcanopy forest [11,12]. In 2002, estimates for the whole archipelago of 21.7% and 24.4% (depending on source) forest cover reflect some slowing down in the rate of deforestation [13,14]. The remaining forest is of critical importance for the survival of many threatened species, including the hole-nesting hornbills (Bucerotidae), which, due to their size, depend on primary or old secondary growth for their reproduction.
The distribution of hornbills in the Philippines largely corresponds to the composition of the PAICs. These birds are a prominent faunal element of the Archipelago such that, based on the occurrence of the different hornbill species, Kemp [15] divided the Philippines into six zoogeographical regions, the Luzonian, the Mindoroan, the Panay-Negrosian, the Samarian the Mindanaoean, and the Suluan, corresponding to the three major and two smaller PAICs (according to Heaney [2], the Samarian and the Mindanaoean formed one PAIC).
The focal species of this study, the Visayan Tarictic Hornbill (Penelopides panini) and the Walden's Hornbill (Aceros waldeni), occur in the Panay-Negrosian region (Greater Negros-Panay, Figure 1), the most threatened faunal region of the Philippines [16]. For comparison, their less endangered (= more abundant) respective sister taxa, the Luzon Tarictic Hornbill (P. manillae) from the Luzonian (Greater Luzon) and the Writhed Hornbill (A. leucocephalus) from the Mindanaoean (Greater Mindanao) are included in this study [15].
A. waldeni and P. panini are endemic to the western islands of the Visayas, which constitute -between Luzon and Mindanao -the central island group of the Philippine archipelago. Due to deforestation and hunting, their populations decreased rapidly [12]. P. panini is known from Panay (including offshore islands), Guimaras, Negros, Masbate, and Ticao Island (here with a second recognized subspecies P. p. ticaensis [17,18]). Today, it is still fairly common on Panay, less so on Negros, and extinct, or nearly so, on all other islands ( Figure 1A) [19][20][21], EC pers. obs. Recently, a remnant population of three birds was found in Milagros, Masbate, but this very small native forest patch is slowly being converted into agricultural land, rendering the survival probability of this remnant population very low [22]. P. panini is listed as endangered by BirdLife International [12].
A. waldeni used to occur on three islands, Panay, Guimaras, and Negros ( Figure 1B) [18], but has survived with a breeding population of substantial size only on Panay, in the forest of the Central Panay Mountain Range [23]. It is extinct on Guimaras and may be extinct on Negros as well, as no recent sightings have been reported [20], EC pers. obs. It is listed as critically endangered by BirdLife International [12].
Captive breeding of both threatened Visayan hornbill species is very rare and incipiently. The world's only captive population of A. waldeni (22 individuals) is held in a breeding facility on Panay. In 2005, for the first time, a total of three offspring were successfully bred by two pairs [24]. P. panini have bred there successfully since 2003 [24]. Former captive reproduction success within the genus Penelopides could have also include P. panini, but due to uncertainty about species level in this genus [15,18], the different species most probably had been hybridized in captivity [25].
The less endangered P. manillae is classified into two subspecies, separated morphologically by size, the width of the tail band and the glossiness of the upperparts [15,18,26]; P. m. manillae occurs on Luzon and smaller adjacent islands, such as Marinduque and Catanduanes, P. m. subniger on Polillo and Patnanungan Islands ( Figure 1A) [18,[26][27][28]. BirdLife International [12] lists the species in the category 'least concern'. The less endangered A. leucocephalus occurs on Mindanao and its satellite islands ( Figure 1B) [18]. Its conservation status is poorly known, but it is considered locally fairly common in suitable habitat and is listed as near threatened by BirdLife International [12].
In order to support conservation efforts for Philippine hornbill species, an assessment of the genetic population structure of the remaining wild populations is a basic requirement. Not many population genetic studies have so far been carried out in the Philippines; this is especially true for birds. Patterns of colonization and diversification in Passeriformes have been studied in a phylogeographic context by Jones and Kennedy [7,29] and Silva-Iturriza et al. [30], and at higher taxonomic levels by Oliveros and Moyle [9] and Sánchez-González and Moyle [31]. Here, we reconstruct the phylogeographic history of the two Penelopides species in order to relate it to the geological history of the Philippines. We also give estimates on the genetic diversity both in the endangered species (A. waldeni and P. panini) and their less threatened sister taxa (A. leucocephalus and P. manillae) in order to test, whether the decline in population sizes in the endangered species is also reflected in lower genetic diversity.

Sampling and DNA extraction
One drop of blood was taken from captive hornbills and stored in Queen's Lysis Buffer [32]. The number of sampled individuals per population, their origin and abbreviation used in the following text and figures are given in Table 1. In addition, 15 feather samples of P. panini and 7 feather samples of A. waldeni were collected from molted feathers in nest holes of wild birds on Panay. Furthermore, 24 museum samples from The Natural History Museum, Tring (U.K.; BMNH) and the Museum für Naturkunde, Berlin (Germany; ZMB) -all dating back to the 19 th century -were included in the study (Table 2).
DNA extraction of blood samples was performed using the DNeasy Tissue Kit (Qiagen) according to the manufacturer's instructions for blood samples. Feather samples were digested with 15 μl proteinase K (20mg/ ml) in the presence of 40 μl DTT and 5 μl Carrier RNA (1 μg/μl) in 450 μl buffer solution (100 mM Tris, 10mM EDTA, 100mM NaCl, 0.1% SDS). After incubation at 55°C over night, we added again 15 μl proteinase K (20mg/ml) and incubated the mix for one more hour. Afterwards, we purified the samples according to the manufacturer's instructions of the DNeasy Tissue Kit (Qiagen), but increased EtOH to 400 μl and the washing solutions to 600 μl. DNA was eluted twice in 50 μl of the buffer solution provided with the kit. DNA extraction from toe-pad samples of museum specimens were performed using the DNeasy Tissue Kit (Qiagen), according to the manufacturer's instructions, adding 5 μl Carrier-RNA (1 μg/μl). In all steps, solution volumes were doubled. DNA was eluted twice, initially in 50 μl and then in 30 μl of the kit buffer solution. DNA extraction from museum samples was carried out in a dedicated ancient DNA laboratory, physically separated from the laboratory where contemporary DNA was handled.

mtDNA amplifying and sequencing
The mitochondrial genome of Philippine hornbills contains two control regions of which we amplified the 5'end of the first control region (CRI) as described for fragment 3 in Sammler et al. [33]. Products were sequenced with primer AcePen_Glu-for to obtain the targeted first~650 base pairs.
To verify the authenticity of amplificates from historical samples, all the mtDNA PCRs of museum samples were repeated by another person in a separate laboratory in a plant genetics department, where animal DNA had never been handled. Specifically, all repeated PCR/ sequencing analyses of the first fragment yielded identical results in 1 st and 2 nd PCR. For the second fragment, 1 st and 2 nd PCR also yielded identical results, except for two occasions with a single bp mismatch. Scarcity of DNA precluded a third PCR, but these two mismatches translate into an per base pair error rate of about 0.014%, a measure in the range of known Taq error rates. There was no indication of an elevated error rate in the historical samples.
Sequences were aligned in BioEdit version 7.0.5.3 [34]. Mitochondrial haplotypes were defined based on a 591-593 bp data partition. Distance matrices and statistical parsimony haplotype networks were constructed using TCS 1.21 with default parameter settings, but a fixed connection limit at 100 steps [35].

Microsatellite genotyping
All blood samples were genotyped at 12-19 previously published polymorphic microsatellite loci (Table 4). PCR and fragment size analysis of listed loci were carried out as described in Sammler et al. [36]. Annealing temperatures of Bbi and Bubi primers were set to 55°C.

mtDNA data analysis
Standard measures of genetic diversity (haplotype diversity, HD [38]; nucleotide diversity, ND [39]) as well as of genetic divergence among island populations (fixation index, F ST [40]), were calculated as implemented in the software package ARLEQUIN version 3.5.1.2 [41]. When more than two populations were compared (i.e., P. panini), we performed an analysis of molecular variance (AMOVA) with the same software. Significance of F ST / AMOVA results was evaluated by permutation analysis (1,000 permutations). Diagnostic mutations of each island population were identified visually in the alignment. ARLEQUIN was also used to assess the fit of the demographic expansion model [42] to hornbill genetic variations (each population singularly). Mismatch distributions were used to fit the model of sudden demographic expansion. Goodness of fit was assessed by the sum of square deviations (SSD) between the observed and the expected mismatch and its significance (P-SSD) was determined by a parametric bootstrap with 10,000 replicates. Mismatch distributions were not calculated for A. leucocephalus because of the small sample size of the single population of this species included in the study.

Microsatellite data analysis
The program FSTAT 2.9.3.2 [43] was used to calculate allelic richness (AR). F ST values were calculated with ARLEQUIN. Significance of these values was evaluated  [45]. Genetic subdivision was evaluated estimating the likelihood and sample composition of independent runs of subgroups (K = 1-5), assuming an admixture model, with a burn-in length of 100,000 and a data collection period of 900,000 iterations. To check for convergence of the Markov Chain parameters, ten replicate runs for each value of K were performed. The value of K best suited for the data was determined according to Evanno et al. [46] and Earl and von Holdt [47]. Demographic signatures of bottlenecks were investigated using the BOTTLENECK 1.2.02 software [48] under three different possible evolutionary models, i.e., Infinite Allele Model (IAM), Stepwise Mutation Model (SMM), and Two-Phase Model (TPM). Estimations of significant deviations from the null hypothesis of mutation-drift equilibrium were based on 10,000 replications.

mtDNA diversity and divergence
Penelopides panini  (Figure 2A). The AMOVA revealed a highly significant overall divergence among populations (F ST =0.572, p<0.001), such that 57% of the mtDNA variation was apportioned to divergence among populations. All single populations had significantly diverged from each others (pairwise F ST between 0.284 and 0.664; p<0.01; see below for details). Three individuals originating from Negros (PpNe11, 40, 55) carry two unique haplotypes that differ by 18-24 substitutions from all other individuals coming from the same area. Both haplotypes are genetically more similar to haplotypes originating from Panay/Guimaras (4 and 5 substitutions away from the respective closest haplotype, Figure 2A). Apart from these three individuals, the two recent populations of Panay and Negros are clearly separated by at least 17 mutational steps in the network and seven diagnostic mutations in the alignment. The fixation coefficient (F ST ) between these populations is 0.665 (p<0.001), but rises to 0.765 (p<0.001) when calculated excluding the three individuals exhibiting deviant haplotypes (Pp11/Pp40 and Pp55). Haplotype diversities range from 0.900 (Negros) to 1.000 (Masbate). The Negros population has a nucleotide diversity higher than that found on Panay (ND = 0.015 vs. 0.010, Table 5), but when the two exceptionally divergent haplotypes (PpNe11/40, PpNe55) are removed from computations, ND drops to 0.004. The highest number of base substitutions between two individuals of the species is 28, corresponding to a nucleotide divergence of 4.7%.
The museum specimens from Guimaras and Masbate are all represented by unique haplotypes that cluster closer Pp_GT_22 P. panini, P. manillae [36] to haplotypes from Panay than to those from Negros ( Figure 2A). There are no diagnostic mutations in the alignment which would clearly separate them from the Panay population. However, individuals of those two populations do not cluster randomly between the Panay haplotypes, but form separate subclusters (Figure 2A)   No haplotype is dominant in frequency ( Figure 2B). Apart from individual Pmm28, both subspecies, P. m. manillae und P. m. subniger, are separated by at least 18 base substitutions in the network and can be distinguished from one another by 7 diagnostic mutations in the alignment. The F ST between the two subspecies is 0.502 (0.529 when excluding Pmm28; both values are significant, p<0.001). P. m. subniger has a lower haplotype and nucleotide diversity than P. m. manillae (HD=0.910, ND = 0.007 vs. HD=0.966, ND=0.037, Table 5). The highest number of base substitutions between two individuals of the species is 70 (11.8% sequence divergence).
In the network, P. m. manillae forms three haplogroups (I, II, III, Figure 2B) separated by large genetic gaps. The first (I) is formed by four museum specimens and Pmm17. Two of the museum specimens (ZMB 9641, ZMB 21772) come from the Manila area, the other two (ZMB 2000.22208, ZMB 32616) are labeled "Culion", an island north of Palawan. Seventeen substitutions separate haplogroup I from the second one (II). This is formed by most of the individuals sampled, including two museum specimens (ZMB 21774, ZMB 2000.22207), which also come from the Manila area. The third haplogroup (III), formed by Pmm24/25, is separated by at least 52 substitutions from haplogroup II. This genetic divergence is by far larger than that between the two morphologically differentiated subspecies P. m. manillae and P. m. subniger ( Figure 2B).
Similar to what was observed in P. panini, we found also in P. manillae one individual (Pmm28) carrying a haplotype that does not fit with its morphological identification. Pmm28 is assigned morphologically to P. m. manillae, but its haplotype clusters closer to those of P. m. subniger (4-11 base substitutions), while there are at least 21 substitutions between Pmm28 and the closest P. m. manillae haplotype. Mismatch distributions never deviate from the purely demographic expansion model (0.156 ≤ P-SDD ≤ 0.371).
With nine individuals, one haplotype is the most abundant; all others are only represented by four or less individuals ( Figure 2C). mtDNA diversity is lower in Guimaras (HD=0.700, ND=0.005) than in Panay (HD=0.897, ND=0.011). The highest number of base substitutions between two individuals is 16 (2.7% sequence divergence). The five museum specimens originating from Guimaras carry three unique haplotypes that cluster among the Panay haplotypes. There are no diagnostic mutations in the alignment which would separate the individuals from Panay and Guimaras; the fixation index is 0.084 (p=0.063). The haplotype of the single museum specimen originating from Negros (BMNH 1897.5.13.471) is separated from the next Panay haplotype by 9 mutational steps in the network ( Figure 2C) and 4 diagnostic mutations in the alignment (Table 5). Mismatch distributions never deviate from the purely demographic expansion model (0.393 ≤ P-SDD ≤ 0.373).

Aceros leucocephalus
Based on 592 bp sequence in 4 specimens, we found 41 polymorphic sites (all transitions). All individuals bear unique haplotypes ( Figure 2D) [GenBank: JX273814-JX273817]. HD is 1.000 and ND is 0.038, both values however associated with high uncertainty (as reflected in SD values), due to small sample size. The highest number of base substitutions between two individuals is 30 (5.1% sequence divergence).

Microsatellite diversity and divergence
The characteristics of the 14 Pp loci and their usability for population genetics analyses are presented by Sammler et al. [36]. Six of the eight analyzed Bbi and Bubi loci [37] were polymorphic, at least in the genus Penelopides. Their number of alleles per locus varied from 1 to 8 in A. waldeni, 1 to 6 in A. leucocephalus, 2 to 13 in P. panini, and 2 to 19 in P. manillae. We confirmed that no linkage disequilibrium could be found between any pair of loci.

Penelopides panini
STRUCTURE results of independent runs with K ranging between 1 and 5 supported two as the most likely number of clusters. The clear division into two groups ( Figure 3A) corresponds to the two populations from Panay and Negros. Intriguingly, the individuals carrying mtDNA haplotypes deviant from their region of origin (PpNe11/40, PpNe55; Figure 2A) are inconspicuous with regard to their microsatellite genotype. Interestingly, in two individuals from Panay (PpPa134 and PpPa78), about 25% of their genotypes is assigned to the Negros cluster.
The number of private alleles and the allelic richness over all loci is higher in the Panay than in the Negros population (PA = 2.11 vs. 0.61; AR = 5.07 vs. 4.78, Table 5) while this pattern is reversed for the observed and expected heterozygosities ( Table 5). In the Negros population no loci deviate from the HWE while two loci are significantly out of HWE in the Panay population due to a lack of heterozygotes (data not shown). The fixation index between these two populations is 0.085 (p<0.001). For the Panay population, BOTTLENECK found a normal L-shaped allele frequency distribution as expected under mutation-drift equilibrium, suggesting that the effective population size has remained constant in the past [48]. For the Negros population, however, a shifted mode is indicative of a recent bottleneck.

Penelopides manillae
STRUCTURE was run with K = 1-5 and supports three as the most likely number of genetic clusters in this species. Individuals morphologically determined as P. m. subniger build a first cluster, and individuals of P. m. manillae build a second cluster, except Pmm24 and Pmm25, which group on their own ( Figure 3B). Pmm24 and Pmm25 also had extremely divergent mitochondrial haplotypes ( Figures 2B and 3). Four individuals (Pms1, Pms9, Pms6, Pms3) can be assigned partially both to the subniger-and the manillae-cluster. The individuals Pmm17 (haplogroup I of P. m. manillae; Figure 2B) and Pmm28, clustering closer to P. m. subniger in the mtDNA network ( Figure 2B), do not stand out at the microsatellite genotypes, but are genetically similar to the other individuals assigned to P. m. manillae on morphological grounds.
The number of private alleles and the allelic richness over all loci is higher in P. m. manillae than in P. m. subniger (PA = 3.789 vs. 0.737; AR = 5.20 vs. 3.53, Table 5). H O is nearly the same in the two subspecies while H E is slightly higher in P. m. manillae than in P. m. subniger (Table 5); one locus deviates from HWE in either subspecies due to heterozygote deficiency (data not shown). The fixation index between these two populations is 0.050 (p<0.001). For P. m. manillae, computations with BOTTLENECK found a normal L-shaped allele frequency distribution as expected under mutationdrift equilibrium (constancy in the past population effective size) [48], whereas for P. m. subniger a shifted mode indicates a recent bottleneck.

Aceros waldeni
Analyses with STRUCTURE did not find any subpopulation structure in the sampling on Panay (K = 1 got best support; data not shown). H E is slightly higher than H O in the species (Table 5) and one locus shows deviations from HWE due to a lack of heterozygotes (data not shown). Computations with BOTTLENECK found a normal Lshaped allele frequency distribution as expected under mutation-drift equilibrium [48].

Penelopides panini
Whereas the microsatellite data for P. panini show a clear division between the two populations from Panay and Negros (the sole exceptions being PpPa134 and PpPa78, Figure 3A), mtDNA revealed that three individuals from Negros (PpNe11, 40, 55) carry two haplotypes that are genetically more similar to haplotypes found on Panay, Guimaras and Masbate than they are to the remaining haplotypes found on Negros (Figure 2A). If we do not consider these two deviant haplotypes and the microsatellite genotypes of PpPa134 and PpPa78, the overall clear-cut subdivision in the genetic architecture of the species would be largely consistent with the geological history of the area. We hypothesize that substantial gene flow between Panay/Guimaras and Negros ceased at the end of the last Pleistocene ice age when the Panay-Negrosian PAIC broke down into separated islands due to eustatic changes in the sea level. While such a Pleistocene sea level-related model of speciation is rejected in favor of an earlier, Pliocene diversification for, e.g., rodents of the genus Apomys [50] or for the fruit bat Haplonycteris fischeri [51], the latter author suggests that sea-level changes may have played an important role for more recent divergence at the within PAICs level, such as among Negros and Panay. The shortest distance between the two islands today is about 15 km but, due to the presence of Guimaras and several other intervening islets, the longest open water distance is just 5 km. Tarictic hornbills are not strong flyers, and group-living, territorial species often make poor candidates for long-distance dispersal, even to offshore islands [52]. However, the three birds PpNe11, 40, 55 presumably descended from at least two female ancestors, which had recently, i.e., some generations ago, migrated from Panay or Guimaras to Negros. Since their haplotypes cluster relatively close to haplotypes from both islands (Figure 2A), the exact origin cannot be ascertained. Due to the general maternal inheritance of the mtDNA in vertebrates including birds (although, in rare occasions, paternal leakage has been found [53]), the haplotypes of these putative immigrants had persisted in the Negros population (Figure 2A), while most of their Panay/Guimaras ancestry has vanished at the microsatellite level due to the bi-parental inheritance of these markers ( Figure 3A).
In the last 100 years, both Panay and Negros suffered from severe deforestation; Guimaras became completely denuded [11,14]. This complete deforestation might have forced birds from Guimaras to abandon the island and to seek still extant forested areas on Negros and Panay. These islands would then still host a (partial) mitochondrial heritage of Guimaras, although the maternal lineages carried by all 19 th century Guimaras museum specimens could not be found in our sampling. These Guimaras lineages have probably gone extinct. However, in this regard, a Guimaras origin of PpNe11/40, 55 seems the most plausible.
Today, the shortest distance between the remaining forest patches of Panay and Negros is about 100 km, a barrier that has become increasingly difficult for forest species to cross. As a consequence, we assume that gene flow between the two areas is currently strongly reduced or completely interrupted. The clear subdivision in the microsatellite data ( Figure 3A) supports this hypothesis, the sole exceptions being PpPa134 and PpPa78.
Interestingly, for another forest bird species, the Philippine bulbul Hypsipetes philippinus, despite a small sample size of only three individuals from Negros, mtDNA suggested multiple independent colonization events for this island and a lack of a direct link between populations from Panay and Negros [30]. These authors, however, did not analyze any nuclear marker. A similar situation, again based only on mtDNA, was found for four out of the five passerine species studied by Jones and Kennedy [7,29]. However, in a fifth species, the Island Verditer-Flycatcher Eumyias panayensis, the populations of Negros and Panay are clearly separated and form sister clades [7]. Generally, dispersal across islands appears much more infrequent in our hornbill species than in the more vagile passerines.
Whereas the population of P. panini on Panay shows no indication of a recent bottleneck, the one on Negros does. This comes along with lower values of genetic diversity in the Negros population (Table 5). Negros retained half of the forest cover of Panay (4% vs. 8%) [11]; in addition, the forest cover on Panay (located mainly in the Central Panay Mountain Range) is almost continuous, while it is quite fragmented on Negros. We are aware that our sample size for this latter population is rather small for a precise estimate of its genetic diversity. We do not have access to 30 individuals, which is recommended for bottleneck analyses by Luikart et al. [54]. Small sample sizes are likely to miss alleles at low frequency [55] and might cause allele frequency distributions to resemble those of a bottlenecked population. Nonetheless, both our genetic analysis and the ecological factors are consistent with a declining and genetically depauperate population on Negros.
Interestingly, not only the historic museum specimens from Guimaras, but also the three from Masbate cluster close to the Panay haplotypes (Figure 2A). The Negros population is thus the only one clearly separated from the other insular populations (except for three individuals: PpNe11, 40,57). Unfortunately, we could not genotype the museum specimens at the microsatellite loci; mtDNA, however, suggests genetic connectedness among Guimaras, Panay and Masbate. In the case of Guimaras, this is not surprising; only a narrow water channel of about 2 km separates nowadays this island from Panay ( Figure 1). Obviously, this channel had not acted as such a strong geographic barrier as the sea between Panay/Guimaras and Negros. In the case of Masbate, it is noteworthy that the over water distance between Panay and Masbate is at least 35 km and thus longer than that between Negros and Panay. We hypothesize that the small islets lying in a row between Panay and Masbate could have acted as stepping stones facilitating dispersal.
Similarly to what was found for the Guimaras specimens, the Masbate specimens also carry unique haplotypes. Thus, Masbate also hosts an independent (yet closely related to Panay) lineage, which still survives with a few last individuals [19,22] or might have already got lost, if the population went extinct since the last census.

Penelopides manillae
Both mtDNA and microsatellites indicate a separation between the two recognized subspecies P. m. manillae and P. m. subniger. However, divergence within one subspecies (manillae) is comparable to divergence between them (see Figure 2B and Table 5). In the haplotype network ( Figure 2B), Pmm28 is an outlier. Although morphologically assigned to P. m. manillae, it clusters within the P. m. subniger haplogroup. The microsatellite profile of Pmm28 (coefficient of membership to manillae of about 90%; Figure 3B) indicates that this bird descends from a female ancestor that migrated from Polillo to Luzon generations ago.
Haplogroup I in the network ( Figure 2B) is formed by Pmm17 and four museum specimens. Two of these museum specimens (ZMB 2000.22208, ZMB 32616) are labeled "Culion", an island north of Palawan. For lack of any historical evidence for an occurrence of P. manillae on Palawan and its satellite islands, which, biogeographically, even do not belong to the Philippines, but to Borneo [4], we cannot rule out the possibility of these specimens having been transferred to Culion by humans. The correct geographic origin of these two individuals is thus considered uncertain. The other two (ZMB 9641, ZMB 21772) come from or around Manila. Two further museum specimens coming also from the Manila area (ZMB 21774, ZMB 2000.22207) cluster in haplogroup II. This implies that in the Manila area two divergent mitochondrial lineages coexist. Such a subdivision is not evident at the microsatellite level, though ( Figure 3B). We cannot rule out the possibility that these results might be biased by our relatively small sample size that could have resulted in missing potentially extant intermediate haplotypes. On the other hand, given the strikingly contrasting pattern yielded by microsatellites, we tend to favor a scenario where two distinct mitochondrial lineages have coexisted for an extended period of time, either because of ancient larger population size or originating from two formerly separated populations within Luzon which later met and merged in the area around Manila.
In contrast, haplogroup III (Pmm24/25) stands out in both the mitochondrial and nuclear data sets, suggesting the occurrence of a further distinct genetic lineage on Luzon. P. manillae has a very wide distribution. Although encompassing islands such as Marinduque and Catanduanes which are as far away from the main island of Luzon as the Polillo Islands are, no further subspecies for those islands have been identified morphologically. Our data show that the lineage represented in our sampling by Pmm24/25 is genetically differentiated from both P. m. manillae (to which it was assigned in the first place) and P. m. subniger. More importantly, it is as distant from P. m. manillae as the latter subspecies is from P. m. subniger. While the instances of deviant mt haplotypes with inconspicuous microsatellite genotypes are indicative of historical gene flow (see above), the microsatellite genotypic compositions of Pms1, Pms9, Pms6, Pms3 might suggest ongoing gene flow, as these specimens presumably represent direct hybrids between a P. m. subniger female (inferred from the mtDNA type) and a P. m. manillae male. Although these samples originate from a non-breeding facility, we cannot exclude with certainty that these hybrids have been produced in captivity [25], such that our inference about naturally occurring gene flow and hybridization should be considered preliminary.
Apart from these hybrids, we hypothesize that P. m. manillae and P. m. subniger became separated when the sea level rose at the end of the Pleistocene and the Polillo Islands became disconnected from the main island of Luzon. While these islands were land-bridge connected to southern rather than central Luzon [2], they are nowadays about 20 km away from central Luzon. In spite of deforestation on Luzon, the coastline facing the Polillo Islands is still afforested, and migration between the islands would be still possible. Although tarictic hornbills are not renowned for long-distance dispersal, even to nearby offshore islands [52], migration events become the more likely, the more the pressure on the Polillo and Patnanungan populations increases due to logging, mining, and conversion of primary forest into coconut plantations [28]. Thus, gene flow might have been directed from the small Polillo Islands to the large island of Luzon by individuals attempting to escape the increasing human pressure. Such a scenario would explain the observed phylogeographic pattern; this, in turn, would be at odds with the theory of island biogeography that presumes large islands to serve as source for dispersal towards smaller islands [56,57].
The small population size of P. m. subniger as compared to that of P. m. manillae might be the result of the small size of the Polillo Islands combined with the loss of adequate habitat. This would explain the inferred recent bottleneck and the low values of genetic diversity of P. m. subniger (Table 5). Again, it should not be overlooked that our sample size for the subspecies is rather small.

Aceros waldeni
In the haplotype network ( Figure 2C), all three haplotypes of the five individuals from Guimaras cluster among the Panay haplotypes. Unfortunately, we do not have microsatellite data for these museum specimens; however, mtDNA (Table 5) suggests a certain degree of genetic connectedness between these two islands. The water channel between Guimaras and Panay is nowadays about 2 km wide and obviously did not act as an insurmountable barrier to gene flow. Between these two populations of A. waldeni, we did not find any major phylogeographic break. Nevertheless, each of the three haplotypes is unique. Thus, this island hosted independent (yet closely related) mitochondrial variants that probably got lost when the population went extinct.
However, despite having analyzed only a single museum specimen from Negros (BMNH 1897.5.13.471), similar to what was observed in P. panini, this individual indicates a separation between the Negros and the Panay/Guimaras population ( Figure 2C). The same processes described for P. panini above may thus also have had an impact on A. waldeni. The Negros population is believed to be extinct or nearly so, as no recent sightings have been reported [EC pers. obs.]. Thus, this differentiated mitochondrial lineage has already gone extinct or is at least threatened with extinction, if any remaining hornbills are still alive.

Conclusions
P. manillae is genetically more structured at the mtDNA level than P. panini. The latter is undergoing a very rapid continuing decline, associated with genetic depletion. Appreciable populations survived only on Panay and Negros. These two populations are genetically differentiated, as are those extinct or almost extinct from Guimaras and Masbate. This means, at least in the case of Guimaras, that part of the original genetic diversity of the species has already got lost.
At the microsatellite level, we found evidences of a recent bottleneck for P. panini from Negros, but also for P. m. subniger (Polillo Islands). While the latter population might be small in size because the islands it colonizes are also small, Negros is a large island and its population shrank solely because of loss of adequate habitat.
For A. leucocephalus, our data suggest that this species is (or was) genetically highly structured. A. waldeni has very rapidly decreased in number [12]. It survives with a substantial population size only on Panay, i.e., in the forest of the Central Panay Mountain Range [23]. A. waldeni generally exhibits lower genetic diversity values, both at microsatellites and mtDNA, than its less endangered sister species A. leucocephalus. Nevertheless, no sign for a recent bottleneck was found for A. waldeni. For the extinct population on Guimaras we could show that an independent (yet closely related) lineage probably has already got lost. From what we found in P. panini,we are inclined to hypothesize that the extinct or almost extinct population of A. waldeni on Negros is most likely genetically differentiated from that still extant on Panay. This hypothesis is further supported by the mt haplotype of the single museum specimen from Negros.
To sum up, the populations of the analyzed hornbill species are generally structured by the geological history of the area, although we found a few indications for across-island dispersal. By comparing patterns of genetic divergence and variability in both endangered Visayan hornbills species to that of their (less endangered) sister taxa, we revealed lower genetic diversity associated with their dramatic population decline and the extinction of genetically differentiated populations. Conservation efforts maintaining the surviving wild birds of virtually extinct populations are therefore particularly important, as well as the preservation of their genetic potential in captivity. This applies especially to the potentially still existing populations of A. waldeni on Negros and P. panini on Masbate.