Genetic Differentiation in Insular Lowland Rainforests: Insights from Historical Demographic Patterns in Philippine Birds

Phylogeographic studies of Philippine birds support that deep genetic structure occurs across continuous lowland forests within islands, despite the lack of obvious contemporary isolation mechanisms. To examine the pattern and tempo of diversification within Philippine island forests, and test if common mechanisms are responsible for observed differentiation, we focused on three co-distributed lowland bird taxa endemic to Greater Luzon and Greater Negros-Panay: Blue-headed Fantail (Rhipidura cyaniceps), White-browed Shama (Copsychus luzoniensis), and Lemon-throated Leaf-Warbler (Phylloscopus cebuensis). Each species has two described subspecies within Greater Luzon, and a single described subspecies on Greater Negros/Panay. Each of the three focal species showed a common geographic pattern of two monophyletic groups in Greater Luzon sister to a third monophyletic group found in Greater Negros-Panay, suggesting that common or similar biogeographic processes may have produced similar distributions. However, studied species displayed variable levels of mitochondrial DNA differentiation between clades, and genetic differentiation within Luzon was not necessarily concordant with described subspecies boundaries. Population genetic parameters for the three species suggested both rapid population growth from small numbers and geographic expansion across Luzon Island. Estimates of the timing of population expansion further supported that these events occurred asynchronously throughout the Pleistocene in the focal species, demanding particular explanations for differentiation, and support that co-distribution may be secondarily congruent.


Introduction
Widespread lowland rainforest bird species are of great interest for studying phylogeographic structure, which links genetics and geography [1,2,3,4]. In regions where lowland rainforest has fewer putative geographic/ecological barriers, species may be expected to be widespread, with limited phylogeographic structure [5,6,7]. However, many bird species exhibit strong genetic differentiation across seemingly continuous lowland forests, suggesting that current or past landscape features restrict gene flow between populations [8,9,10,11,12,13,14,15]. In the naturally fragmented insular landscapes of Southeast Asia, widespread lowland rainforest species often show extensive phylogeographic structure associated with repeated cycles of connection and isolation of island complexes across the region during Quaternary climatic changes [16,17,18,19,20]. However, the possibility of isolation and genetic differentiation within continuous lowland forests has received relatively little attention compared to isolation across marine barriers (but see [21,22] and references therein).
The Philippine archipelago provides a natural model to test for genetic differentiation in continuous habitats within islands. With the exception of Palawan and some of its offshore islands which were apparently united to the Sunda Shelf (reviewed in [46]), the Philippine archipelago is oceanic in origin (reviewed in [47]). After its origin by complex geological activity, the archipelago was subjected to climatic and sea-level changes in the Pleistocene [48,49] that produced cycles of isolation and aggregation of islands into larger landmasses, providing opportunities for both dispersal and isolation [48]. These Pleistocene Aggregate Island Complexes, or PAICs [49] were apparently never joined to other PAICs, because they were separated by deep water channels (>120 m, [48]).
The dominant paradigm for explaining the high diversity levels in the Philippines has been based on the arrangement of exposed land among aggregated islands during Pleistocene sealevel changes (reviewed in [32]), predicting that most speciation events among PAICs are due to vicariance [16]. Recent work, however, has shown that the PAIC paradigm may not explain many observed diversity patterns (reviewed in [32]), and that dispersal or within-PAIC differentiation may have played a more important role than expected. We tested the prevalence of the PAIC paradigm in three primarily lowland rainforest passerine birds: Blue-headed Fantail (Rhipidura cyaniceps, Rhipiduridae), White-browed Shama (Copsychus luzoniensis, Turdidae) and Lemon-throated Leaf-Warbler (Phylloscopus cebuensis, Phylloscopidae). These taxa are codistributed within the Greater Luzon PAIC, which includes the present-day islands of Luzon, Polillo, Catanduanes, and Marinduque; and in the Greater Negros-Panay PAIC, which includes Negros, Panay, Cebu, Ticao, and Masbate islands. Each species is common in forested habitats on these islands, and they are regularly found together in mixed species flocks. Both inter-and intra-island plumage differences have been documented in these species, as suggested by recognized subspecies throughout their range (Fig 1, [29,31]). In each species, multiple described subspecies are endemic to Greater Luzon (two in R. cyaniceps and P. cebuensis, three in C. luzoniensis), and single subspecies are endemic to Greater Negros-Panay.
We used tools from phylogenetics and phylogeography to estimate historic demographical patterns in these co-distributed Philippine passerines. Given the extensive effects of Pleistocene climatic changes in Southeast Asia [48,50,51,52,53], we may expect substantial changes in their geographic distribution, and consequently, that some populations may have experienced demographic changes. However, the effects of Pleistocene climatic changes in the oceanic Philippines remain little studied; thus, phylogeographic approaches may be of great utility in proposing primary hypotheses about the evolutionary history of the Philippine biota. Results were used for discriminating between different scenarios likely responsible for producing phylogeographic structure, such as the PAIC paradigm [48]. Based on a previous study that documented within-PAIC differentiation [54], we tested the following hierarchical biogeographic hypotheses for the co-distributed species in this study (Fig 2): a null hypothesis (H 0 ) corresponding to the classical PAIC paradigm, in which genetic differentiation is expected to be partitioned between clades restricted to PAICs; alternative hypotheses involving genetic differentiation as a result of different colonization events, resulting in unrelated clades in a single PAIC (H 1 ), and within-island genetic differentiation, in which sister clades are distributed in a single PAIC (H 2 ).

Taxon sampling
The three focal species are endemic to the Greater Luzon + Greater Negros-Panay PAICs. In addition to these islands R. cyaniceps is also found on Tablas Island (Romblon Island Group). Neither R. cyaniceps nor C. luzoniensis are present on Cebu (Fig 1). Of the total samples, 40 corresponded to Copsychus luzoniensis, 28 to Phylloscopus cebuensis and 50 to Rhipidura cyaniceps (S1 Table, see Acknowledgements). To ensure thorough assessment of differentiation patterns within Greater Luzon, we sequenced multiple individuals from as many localities as possible. A Gratuitous Permit (GP) to conduct research and collect specimens in all sampling localities was issued by Mundita S. Lim, Director of the Biodiversity Management Bureau (BMB), Republic of the Philippines. Field studies did not involve endangered or protected species. Birds were captured using mist-nets. Nets were checked every hour, with birds immediately released if not needed for future study. Birds were euthanized via thoracic compression or isoflurane open-drop. This project operated under the University of Kansas Animal Care and Use Committee (IACUC approval AUS no. 174-01), issued to R.G.M. at the University of Kansas.
Because species-level relationships of R. cyaniceps and C. luzoniensis have been studied previously, outgroup choice was straightforward. Species selected as outgroups for C. luzoniensis were Copsychus niger from Palawan Island, Philippines, and Copsychus malabaricus, a widely distributed species in southeastern Asia [19,20]. Outgroup selection for R.cyaniceps relied on two recent studies [55,56], showing that the three traditionally recognized endemic Philippine Rhipidura (including cyaniceps) are a monophyletic group, thus we included samples of the other two endemic Philippine species. Finally, because the monophyly of P. cebuensis has been questioned, with some authors lumping this species and P. olivaceus into a single taxon [57,58], we included samples of the latter, which is likely the sister species of P. cebuensis [59]. We also included samples of P. ijimae and P. coronatus as additional outgroups [60].  Heaney (1985), study areas are in dark grey. Dotted circles represent sampling localities (See S1 Table for details). Maps show distribution for each taxa, and the subspecies described (Dickinson et al. 1991). Dashed lines represent putative borders for currently accepted subspecies in Luzon. Intra-Island Differentiation in Philippine Birds DNA sequencing DNA was extracted from frozen tissue using Proteinase K digestion procedures following the manufacturer's protocols (DNeasy; Qiagen, http://www.qiagen.com/). Markers selected for this work have been used widely in bird systematics and biogeography, and include the entire second subunit (ND2, 1041 bp) and third subunit (ND3, 351 bp) of nicotinamide adenine dinucleotide dehydrogenase. These markers were amplified with the external primers L5215 and H6313 for ND2 [61] and L10755 and H11151 for ND3 [62], as well as the ND2 internal primers 487L (provided by C. Oliveros C., unpublished), ND2-SWH [63], and H5766 [61].
Genomic DNA was amplified using 5-primeTaq DNA polymerase under standard PCR thermocycling protocols and visualized in agarose gels stained with ethidium bromide. Resulting products were cleaned with ExoSAPIT (GE Healthcare Corp.) and the purified products were cycle-sequenced with ABI Prism BigDye v3.1 terminator chemistry. Cycle-sequenced products were purified with ethanol precipitation, and sequenced on an ABI 3730 automated sequencer. Sequences were aligned using MAFFT [69], as implemented in Geneious 7.0.2 [70]. Alignments for each gene were further inspected by eye, and insertions and deletions (indels) were adjusted as necessary. Heterozygous positions in the nuclear introns, determined by double-peaks of similar height, were coded following the IUPAC ambiguity codes.

Phylogenetic Analyses
Phylogenetic relationships were reconstructed from the concatenated dataset for each taxon through Maximum Likelihood analysis (ML) as implemented RAxML 7.0.3 [71], which allows for different sequence evolution models to be incorporated into the analysis; nodal support was assessed via non-parametric bootstrapping [72] with 1000 replicates. We also used Bayesian Inference (BI) on the complete dataset for each species using MrBayes 3.2. [73]. Each dataset was partitioned by gene and codon positions for the nuclear intron and mitochondrial genes respectively [74,75]. The Akaike Information Criterion (AIC), as implemented in MrModeltest [76], was used to determine the best substitution model for each partition. BI was implemented for 10 7 generations and sampled every 500 generations. Stationarity of the MCMC chains was assessed in Tracer v1.5.0 [77], after which the first 30% generations were discarded as initial burn-in. All remaining trees in the summary were used to produce a single 50% majority-rule consensus tree. In order to ensure that examination of tree space was appropriate, topological convergence was assessed in the on-line application AWTY [78] by using the compare function, which plots posterior probabilities of all splits for paired MCMC runs. Inspection for stationarity revealed that parameter and topological space were searched thoroughly.

Population genetic parameters
All population parameters were estimated from the mitochondrial dataset (ND3 and ND2) for each species. Genetic diversity was assessed using indices of haplotype diversity (Hd) and nucelotide diversity (π), with samples grouped by present-day island boundaries (Fig 1).
Following the phylogenetic results, Luzon samples were further divided according to the obtained clades, corresponding to northern Luzon and the Bicol Peninsula. However, the geographic structure in these clades showed the Zambales region (in western Luzon) samples grouped with the northern Luzon clade in both R. cyaniceps and P. cebuensis, but with the Bicol Peninsula clade in C. luzoniensis. These arrangements suggest particular dynamic biotic processes in the Zambales region, for which we conducted separated analyses for the populations of the three species. As a relative measure of divergence between populations, we estimated Nei's genetic distance values (Dxy), using DNAsp v5 [79] with a Jukes-Cantor correction [80], between groups of populations identified in the phylogenetic results.
Genetic structure in the three species was explored using three-way AMOVA and Fst. For the AMOVA analysis, samples were arranged in groups corresponding to the clades found in the phylogenetic analyses. The significance of the AMOVA results was assessed through 10,000 non-parametric permutations. The parameter Fst was calculated through pairwise differences between haplotypes; significance of the Fst parameter was also assessed with 10,000 permutations. Finally, as an additional measure of gene flow, we calculated the number of migrants per generation (Nm). Because some studies have shown that populations of R. cyaniceps and C. luzoniensis in the Greater Negros-Panay PAIC and the Romblon Island Group are evolutionarily distinct from those in Greater Luzon [20,55,81], AMOVA analyses were repeated using the same parameters as described above but excluding populations in these island groups. All statistics were calculated in Arlequin ver. 3.5.1.3 [82]. Interpretation of Fst and Nm values followed the guidelines in Hartl and Clark [83].
We tested whether populations of the three species had experienced demographic changes by calculating Fu's F s statistic [84], which indicates whether individual populations are evolving according to the Wright-Fisher model. Significance of Fu's F s was determined using a p-value of 0.02 as suggested by Fu [84]. We also calculated Tajima's D statistic [85], another measure of the selective neutrality of markers in a population. Significance of Fu's F s and Tajima's D were calculated by constructing 1000 coalescent simulations in DnaSP v5 [79].
Population history was further inferred by plotting mismatch distributions [86,87] and calculating their significance using Ramos-Onsins and Rozas R 2 , which is better suited for small sample sizes [88]. R 2 significance was estimated through 1000 coalescent simulations for the different clades (and the Zambales subpopulation) of each species in DnaSP v5 [79]. For clades showing significant population growth, the parameter Tau (τ) was used to calculate the time t of potential step-wise expansion from a relatively small, but constant population to a large population of size Ɵ 1 over t generations in the past, with t = τ/2u, where τ = age of expansion (in mutational units) and u = 2μk, where μ = mutation rate and k = the length of the sequence [87]. We used mutation rates of 1 x 10 −9 substitutions/site/year (s/s/y), which is supported by a large dataset of studies on passerine birds [89], 2.7 x 10 −9 s/s/y, which is a ND2 specific mutation rate calculated from mockingbirds in the Galapagos Islands [90], and a faster rate of 4 X10 -9 s/s/y, which accounts for uncertainties about mitochondrial substitution rates [91,92]. Some researchers have suggested that the estimation of τ according to Rogers [93] often leads to an underestimation of the age of expansion, due to the omission of heterogeneity in mutation rates among sites [94]. Thus, we calculated τ values by applying a least-squares approach, as implemented in Arlequin ver. 3.5.1.3 [82]. Confidence intervals for τ were calculated using 3000 parametric bootstrap replicates [94], as implemented in Arlequin ver. 3.5.1.3 [82].
Finally, because it has been suggested that bifurcating trees may not always fully represent intraspecific phylogenies due to the coexistence of ancestral and derived haplotypes in a given sample [95], we also estimated Median-joining networks using a median-joining method [96] in Networks 4.6.0.0 (http://www.fluxus-engineering.com), assigning equal weights to all variable sites and with default values for the epsilon parameter (Ɛ = 0).

Phylogenetics
In all, 1967 characters were included in the complete dataset for R. cyaniceps, 2198 characters for C. luzoniensis, and 2429 characters for P. cebuensis. Sequences for R. cyaniceps, C. luzoniensis, and P. cebuensis are deposited in GenBank (S1 Table). Sequence characteristics for each gene partition and the selected models of evolution are provided in S2 Table. Phylograms from the mitochondrial DNA (mtDNA) dataset were identical to those obtained from the complete dataset in R. cyaniceps (Fig 3, see [55]), as well as for C. luzoniensis (Fig 4) and P. cebuensis (Fig 5). Nuclear phylograms (nDNA) for each species showed two sister clades including all of the Greater Luzon PAIC samples and Greater Negros-Panay PAIC (Negros and Panay islands) samples. However, no structure in the nDNA datasets was apparent within Greater Luzon, likely due to the fourfold higher effective population size, longer coalescence times, male biased dispersal, and slower mutation and rates of nDNA markers, consistent with expectations from coalescent theory [97].
ML and BI results were largely congruent in all three cases (Figs 3, 4, and 5). All analyses showed each of the three taxa as monophyletic. In each case, there were three relatively wellsupported clades: a Western Visayas clade (Greater Negros-Panay PAIC), and two clades from Greater Luzon, corresponding to northern and southern parts of the island. Samples from the Zambales region in western Luzon were embedded within the northern Luzon clade in R. cyaniceps and P. cebuensis, but with the Bicol Peninsula clade (including Catanduanes Island) in  C. luzoniensis. One additional monophyletic grouping included the Tablas island samples of R. cyaniceps, as previously reported [55].

Phylogeography
Genetic structure (Table 1) analyses revealed high values of differentiation between PAICs. High values for Nei's corrected distance (Dxy) were found for comparisons between Greater Luzon PAIC and Greater Negros-Panay PAIC for R. cyaniceps and C. luzoniensis, consistent with recent systematic treatments in which these populations have been recognized as full species (R. albiventris and C. superciliaris, respectively, [20,55,81,98]). Divergence also was found within PAIC boundaries, as in Greater Luzon PAIC, where values of genetic flow and genetic differentiation are also high [82]. High levels of haplotype diversity (Hd) and low levels of nucleotide diversity (π) were observed in the populations of the three species (Table 2), suggesting rapid population growth according to Category 2 in Grant and Bowen [99].
Mismatch distributions (Fig 6) and Ramos-Onsins and Rosas´R 2 values (Table 2) showed similar patterns with respect to geographic areas across species. In northern Luzon, mismatch distributions for R. cyaniceps and P. cebuensis showed a unimodal pattern, suggesting population expansion; R 2 values were significant in both species. In C. luzoniensis, a ragged pattern is apparent, suggesting that populations may be at demographic equilibrium. This same unimodal pattern is also apparent in Bicol Peninsula populations of R. cyaniceps and C. luzoniensis, in which significant R 2 values were obtained for both species, suggesting recent population expansion in the area. For P. cebuensis, the Bicol Peninsula clade showed a multimodal pattern, suggesting persistence through small and isolated populations or a population bottleneck [100]. The Zambales populations of both R. cyaniceps and C. luzoniensis (small sample size in P. cebuensis prevented analysis) showed a bimodal pattern, suggesting either that for both species, there has been a relatively constant effective population size in the past [100], or probably, an admixture of distinct populations, which may be expected as both of these subpopulations are part of larger genetic groups located either in Northern Luzon and the Bicol Peninsula, respectively. Finally, the Greater Negros-Panay in R. cyaniceps and C. luzoniensis showed a multimodal pattern, suggesting either population bottleneck or persistence of small and isolated populations, which seems highly probable as this PAIC is presently partitioned in smaller islands fragments, in comparison to Greater Luzon PAIC. However, these results are only preliminary and should be taken with caution, given the sample size and geographic coverage, which prevented conclusive demographic inferences, especially for P. cebuensis.
Three-way AMOVA results (Table 3) suggested that the greatest variation in genetic structure of the three species is found among groups, which correspond to clades obtained in the  [83] were obtained in all comparisons involving populations from Greater Negros-Panay PAIC and populations in Greater Luzon PAIC (Table 2), as predicted by the PAIC paradigm [49]. However, high and significant Fst values (>0.72) were also obtained for within-PAIC comparisons, particularly involving Northern Luzon and the Bicol Peninsula populations in the three species. Values of Fst from populations in Zambales and Northern Luzon showed no significant genetic structure in R. cyaniceps and P. cebuensis, but were moderate (0.21) for C. luzoniensis populations from the Zambales and Bicol Peninsula. Gene flow was low [83] among most populations within the three focal species (Table 1), suggesting less than one migrant per generation. However, gene flow was more prevalent between populations of R. cyaniceps in northern Luzon and the Zambales (4.11 migrants per generation), and between populations of C. luzoniensis from Zambales and Bicol (1 migrant per generation). Population expansion in Northern Luzon was supported by significantly negative values obtained for Fu's F s [84] in both R. cyaniceps and P. cebuensis (Table 2). Remaining populations may have experienced population growth; however, small sample sizes may influence statistical power to detect them. Tajima's D values [84] corroborated the results obtained with the Fu's F statistic, also suggesting population expansion or stabilizing selection [99]. Additionally, the Bicol Peninsula clade of R. cyaniceps showed significant negative values, suggesting also demographic expansion in this population. Although estimated population expansion dates are heavily dependent on the mutation rates used, all of them and their respective confidence intervals fall within the Pleistocene (Table 4). Median population expansion times using the standard 2% [89] rate support early

Intra-Island Differentiation in Philippine Birds
Pleistocene for R. cyaniceps (1.5 mya) and early to middle Pleistocene for P. cebuensis (0.8 mya). In the Bicol Peninsula, estimates of population expansion strongly suggest two different events. For C. luzoniensis, estimates suggest mid to early Pleistocene events, whereas populations of R. cyaniceps may have expanded in the late Pleistocene. Mitochondrial haplotype networks (Figs 3, 4, and 5) for the three species reflected the same population structure revealed by the phylogenetic analyses, in which high (in R. cyaniceps and C. luzoniensis, 40 and 41 mutational steps respectively) to low (3 in P. cebuensis) numbers of mutational steps separated Greater Negros-Panay populations from Luzon populations. Luzon birds were grouped in two different haplotype clusters in R. cyaniceps and C. luzoniensis, in which 19 and 7 mutational steps, respectively, separated the Bicol Peninsula clades from Northern Luzon clades. Haplotype networks in the three species showed a star-like pattern, suggesting historical demographic expansion after isolation [3]. Apparent signs of secondary contact for C. luzoniensis were observed in the Mingan Mountains in eastern Luzon, where analyses showed representatives of north and south haplotype groups. One of the haplotypes  nested with the Bicol peninsula samples, the other nested with Northern Luzon samples (Fig  4).

Discussion
The observation that endemic species' distributions in the Philippines corresponded to biogeographic subprovinces produced a paradigm that dominated biogeographic and taxonomic research for over 25 years, in which repeated cycles of connection and isolation within different PAICs occurred in concert with climatic shifts and sea level change in the Pleistocene, promoting biotic evolution in the archipelago (reviewed in [33]). Although the PAIC paradigm provides a functional explanation for biotic evolution in the Philippine archipelago, molecular analyses have suggested that the this paradigm only explains a portion of diversity patterns in the archipelago, highlighting mixed models as potentially better explanations [18,33,39,54]. Genetic differentiation in the three lowland passerine species in this study also corresponds to a mixed model. Consistent with the PAIC paradigm (H 0 ), phylogenetic trees of the three species showed monophyletic sister groups in Greater Luzon and Greater Negros-Panay PAIC. The pattern departs from strict PAIC expectation within Luzon, with substantial structure discovered in this clade in all three species (H 2 ). Populations in each of the species showed levels of genetic divergence higher than those proposed for the recognition of full species [11,63,101,102], except for P. cebuensis. Genetic differentiation within Luzon is particularly evident in R. cyaniceps and C. luzoniensis, where high Fst values are coupled with a low estimated number of migrants, suggesting the existence of a past barrier to gene flow within the island. Additionally, the association of samples from the Zambales region either with northern Luzon (in R. cyaniceps and P. cebuensis), or the Bicol Peninsula, as in C. luzoniensis suggests that this region in western Luzon has had a dynamic biotic history, in which lowlands may have connected at different times with Northern Luzon or the Bicol Peninsula lowland rainforests.
Following the premises of the PAIC paradigm, most attention on bird differentiation in the Philippines has been devoted to patterns among PAICs [29]. In contrast, within-island differentiation has received relatively little attention [59,103,104,105]. Phylogeographic results in this study showed that genetic differentiation also has occurred within the limits of Greater Luzon PAIC, as in the Zambales region, where phylogenetic patterns suggest a dynamic history, with populations of different species having close relationships to populations on opposite ends of the island.
All of the estimated dates and confidence intervals for population expansion suggest that these events occurred throughout the Pleistocene, spanning several proposed glacial cycles [106], corresponding with models of rainforests expansion and contraction during glacial periods in southeastern Asia [52]. Evidence for demographic expansion in populations of the three species support a scenario in which isolation occurred in different refugia located in Northern Luzon for R. cyaniceps and P. cebuensis and in the Bicol Peninsula for C. luzoniensis and R. cyaniceps. In the case of C. luzoniensis, population expansion is further supported by secondary contact between the two distinct phylogroups in the Mingan Mountains, in eastern Luzon.
Considered together, intra-PAIC differentiation patterns and evidence of demographic expansion suggest that populations of some species were fragmented during at least some Pleistocene climatic cycles [49,59]. A recent study, however, suggested that lowland rainforest connectivity in the Philippines increased during the Pleistocene Last Glacial (LGM, [53]). Although forest may have been widespread in the LGM, differentiation and even population expansion in the focal species likely predated the onset of that period in the three, thus suggesting the main effect of lowland rainforest connectivity was probably biotic redistribution. Although forest connectivity may have increased during glacial periods in some regions [53], the structure and communities of these lowland rainforests may have been different from present-day communities [107,108], rendering the habitats unsuitable, which may support that populations of the three focal species in Luzon have maintained their differentiation due to ecological vicariance [43,55]. Environmental niche modelling research has shown that response to environmental change is species-specific [109,110] and that modifications in some of the abiotic variables may influence distributional patterns, which in turn may have promoted genetic divergence [54]. Also, despite forest connectivity, genetic distances and Fst values for the focal species closely resemble the patterns detected in montane forest taxa [59,103,105], for which geographic isolation has played a key role (reviewed in [111]).
Taken together, phylogenetic and phylogeographic patterns in these lowland rainforest endemic birds point to a more dynamic evolution of the Philippine archipelago biota than predicted by the PAIC paradigm. For the three species in this study, a model incorporating the effects of sea level change and climatic changes seems adequate to explain historical demographic patterns. Although sea level changes may have functioned as barriers promoting divergence between PAIC populations, habitat isolation due to climatic changes may have restricted (and perhaps reinforced in different glacial cycles) gene flow in populations within the same PAIC [42,54]. Results in this study and other studies have revealed that diversity in the Philippines is not only structured between PAICs, but also within PAICs. This result is consistent with recent studies in mammals and reptiles (e.g. [39,42,46]), and birds [54,54,59,101,103,105,112].
Our study adds to the growing body of evidence showing that speciation and genetic differentiation may occur within a single island, probably as a consequence of the Pleistocene climatic changes. Recent studies for birds (reviewed in [21]) and mammals in Borneo (e.g. [113]) shown that speciation and differentiation may be more common than expected, either in landbridge islands in the Sunda Shelf, or in oceanic islands, as in the Philippine archipelago [36, 37, 38, 39, 40, 42, 43, 46, 104, 105, this study], where lowland differentiation in birds is widespread within the two largest Philippine islands of Luzon [54,104] and Mindanao [54], yet it remains undocumented in smaller islands.

Taxonomy
Results in this and other studies have repeatedly suggested that the current taxonomy for the Philippine biota does not accurately reflect diversity in the archipelago, which means that current species diversity estimates are better interpreted as conservative [114,115]. Birds have been long considered as a group with a relatively good taxonomic knowledge; however, the application of the biological species concept in allopatric contexts, such as the Philippine archipelago, has been controversial because of difficult inferences about reproductive isolation between geographically isolated populations [116].
Phylogenetic patterns showed two main clades: one grouped all of the Greater Luzon samples; the other grouped all of the Greater Negros-Panay samples. High genetic distances and Fst values indicate substantial differentiation between these populations in all focal species, with the exception of P. cebuensis (1.3% average). Populations of the other two taxa in Greater Negros-Panay both show deep genetic differentiation from Luzon populations and diagnostic characters that allow them to be recognized as full species [20,55].
Additional genetic variation was also discovered within the bounds of a single PAIC. In Greater Luzon PAIC, phylogenetic analyses showed two well-supported monophyletic sister groups. This arrangement apparently agrees with current taxonomy, as there are two recognized subspecies for each taxon [29,57]; however, geographic ranges between recognized subspecies and phylogeographic groups do not match. In the case of R. cyaniceps, the two recognized subspecies cyaniceps and pinicola were included in the same clade, suggesting that morphological variation may be clinal [117] and even ecological, as pinicola is mainly restricted to pine forests [29,117]. The second clade includes all of our samples from the Bicol Peninsula. An average genetic differentiation of 3% (range 2.9-3.1%), and high Fst values (average 0.9, range 0.85-0.95) between the two groups suggests a long period of genetic isolation. Disagreement between current taxonomy and phylogeographic groups was also found in C. luzoniensis. Samples from Luzon Island were included in two sister clades, with populations from northern Luzon sister to those from western, central, and southern Luzon (Bicol Peninsula), and Catanduanes Island. Deep genetic divergence (average 2.9%, range 2.4-3.3%) and high Fst values (average 0.72, range 0.71-0.73) also suggest long isolation and consequent differentiation. Finally, the only case where taxonomy apparently agrees with our work is in P. cebuensis. The two recognized subspecies luzoniensis and sorsogonensis seem to match the monophyletic groups obtained; this subespecific arrangement is supported by a low genetic differentiation (average 0.3%, range 0.3-0.4%) but high Fst values (average 0.86, range 0.7-1).
Throughout the northern Philippines, molecular studies in birds have found values of genetic divergence ranging from 2.7% to 13.8% between populations of the same taxon, suggesting that species boundaries in a number of avian species should be revised [11,19,55,105,118,119]. These values bracket the genetic divergence found in the two clades of R. cyaniceps and C. luzoniensis found in Luzon Island. It has been suggested that evidence from a single dataset may not be enough, underscoring the need for additional evidence such as morphological and song characters [120,121]. In the focal species, morphological differentiation is evident when comparing Greater Luzon and Greater Panay Negros lineages, as each has diagnostic characters. However, this task is complicated when comparing lineages within Luzon, as they look almost identical at first sight. However, a cursory inspection in C. luzoniensis revealed that birds from Bicol and the Polillo islands differ from the northern Luzon birds in the size of the white spots on the undertail (a trait that allowed the recognition of the Polillo taxon parvimaculatus), suggesting that diagnostic characters are present. In R. cyaniceps, birds from the Bicol Peninsula apparently have darker plumage [116]. This situation is reversed in P. cebuensis, which has the lowest genetic divergence but clear diagnosable characters; subspecies in Luzon differ in the amount and brightness of yellow in the throat and undertail coverts, being brighter in the Bicol Peninsula sorsogonensis [122].
Species delimitation may be controversial due to different philosophies and species concept applicability [123]. This may be even more complicated when genetic divergence is not accompanied by clear morphological difference and when this divergence has occurred within the same island e. g. [124]. However, genetic differentiation may occur without corresponding morphological differentiation [125], and recent work in the Philippine archipelago has demonstrated that speciation has occurred within a single island [36,37,38,39,40,42,43,46,104,105]. Whatever philosophy is applied, genetic differentiation and speciation has occurred within a single PAIC, thus deviating from classic PAIC differentiation expectations.
Supporting Information S1