INTRODUCTION
The distribution of many highland Neotropical bird species is naturally fragmented as a consequence of past climatic changes (Haffer Reference HAFFER1974, Vuilleumier Reference VUILLEUMIER1970, Weir Reference WEIR2006). The glacial–interglacial cycles of the late Pleistocene, and in particular the climatic shifts after this epoch, had a drastic impact on dispersal, speciation and distribution of the avifauna in the Neotropical highlands. Climatic changes fragmented the more continuous distribution of highland species into a series of disjunct populations restricted to the upper elevations of mountain ranges (Barrantes Reference BARRANTES2009, Fjeldså Reference FJELDSÅ, Balslev and Luteyn1992, Fjeldså & Krabbe Reference FJELDSÅ and KRABBE1990, Vuilleumier Reference VUILLEUMIER1970, Weir Reference WEIR2006). By being isolated by topographic discontinuities, species numbers in many of these highland habitat islands (hereafter called highland islands) are affected by a series of factors, in some cases similar to those affecting insular communities, such as island area, distance from a source of species, distance from other islands, and the interaction of biological processes with the geological dynamics of islands (MacArthur & Wilson Reference MACARTHUR and WILSON1967, Nores Reference NORES1995, Rosenzweig Reference ROSENZWEIG1995, Simberloff & Abele Reference SIMBERLOFF and ABELE1982, Vuilleumier Reference VUILLEUMIER1970, Vuilleumier & Simberloff Reference VUILLEUMIER and SIMBERLOFF1980, Whittaker et al. Reference WHITTAKER, TRIANTIS and LADLE2008). Dynamic processes such as local extinction and colonization are expected to be responsible for non-random patterns of species composition in a large number of these isolated continental assemblages (Patterson & Atmar Reference PATTERSON and ATMAR1986, Patterson & Brown Reference PATTERSON and BROWN1991). Often smaller assemblages contain successive subsets of the species in larger assemblages following a nested pattern (Patterson Reference PATTERSON1990, Patterson & Atmar Reference PATTERSON and ATMAR1986).
Glacial cycles during the late Pleistocene favoured the dispersal of birds from South and North America to the recently formed highland mountains in southern Central America (Barrantes Reference BARRANTES2009, Hackett Reference HACKETT1995, Haffer Reference HAFFER1974, Lowell et al. Reference LOWELL, HEUSSER, ANDERSEN, MORENO, HAUSER, SCHLUCHTER, MARCHANT and DENTON1995, Pérez-Emán Reference PÉREZ-EMÁN2005, Prance Reference PRANCE1982, Stiles Reference STILES and Janzen1983a, Winker & Pruett Reference WINKER and PRUETT2006). The subsequent switch to a more tropical climate after the Pleistocene, and probably during other interglacial periods during this epoch, caused a contraction of the distribution of highland vegetation in these mountains (Hooghiemstra et al. Reference HOOGHIEMSTRA, CLEEF, NOLDUS and KAPPELLE1992, Islebe et al. Reference ISLEBE, HOOGHIEMSTRA and VAN DER BORG1995, Reference ISLEBE, HOOGHIEMSTRA and VAN VEER1996) – and likely also of the avifauna associated with this vegetation – to the upper elevations of high mountains (Barrantes Reference BARRANTES2009). The reduction of the original distribution of highland birds had two major effects. First, it increased the isolation of the highland avifauna of Costa Rica and western Panama from similar South and North American avifauna. Second, it caused a fragmentation within the highlands of Costa Rica and western Panama into highland islands, separated by geographic and climatic barriers (Hackett Reference HACKETT1995, Puebla-Olivares et al. Reference PUEBLA-OLIVARES, BONACCORSO, ESPINOSA DE LOS MONTEROS, OMLAND, LLORENTE-BOUSQUETS, PETERSON and NAVARRO-SIGÜENZA2008).
The effect of this isolation of Costa Rican montane areas has long been recognized on the basis of morphological differences across populations of several avian species (Barrantes & Sánchez Reference BARRANTES and SÁNCHEZ2000, Slud Reference SLUD1964, Stiles Reference STILES1983b, Reference STILES1985). However, the effect of area of available habitat on the number of highland endemic species (e.g. the species–area relationship predicted by the Island Biogeography Model) and the effect of past reduction of available habitat on the composition of species assemblages has never been tested using data from Costa Rica's highland islands. To this end, we first tested the hypothesis that area predicts the number of endemic and non-endemic bird species present on highland islands. We then tested the hypothesis that the reduction in highland vegetation that occurred during the last interglacial periods affected the composition of endemic birds on highland islands (e.g. nested distribution vs. independent subsample of species). Finally, we hypothesized that differential species extinction due to reduction of island area had a larger effect on the composition of species on each highland island than dispersal limitation due to barriers between islands.
HIGHLAND ENDEMIC AVIFAUNA
The Costa Rican–western Panama highland endemic avifauna, as defined by Wolf (Reference WOLF1976), consists of those species that occur in the Talamanca mountains above 2400 m. The habitat used by these endemic species includes upper montane forest, subalpine forest, elfin forest and paramo (Barrantes Reference BARRANTES, Kappelle and Horn2005, Barrantes & Loiselle Reference BARRANTES and LOISELLE2002, Stiles Reference STILES and Janzen1983a, Reference STILES1985). Endemic species are also present in the other three Costa Rican mountain ranges (Central, Tilarán and Guanacaste). In the Talamanca mountain range, the lower limit of some endemic species is below 2400 m (Stiles & Skutch Reference STILES and SKUTCH1989), but in all cases their abundance is higher above this altitude (Wolf Reference WOLF1976). Populations of most highland endemic species are naturally fragmented by geographic discontinuities, such as mountain passes and watersheds (Barrantes Reference BARRANTES2009, Chavarría-Pizarro et al. Reference CHAVARRÍA-PIZARRO, GUTIÉRREZ-ESPELETA, FUCHS and BARRANTES2010).
METHODS
Field data collection
We conducted bird surveys from December 1993 through January 2005 in highland forests of the four main mountain ranges in Costa Rica to estimate species richness: four sites in the Talamanca mountain range, three in the Central Volcanic mountain range, one in the Tilarán mountain range, and two in the Guanacaste mountain range (Figure 1). Sites included the highest peaks as well as the northern and southern limits of each mountain range, except in the Guanacaste mountain range where volcanic activity precluded sampling in a large portion of highland vegetation along its southern boundary. We visited each site from 7 to 20 times and during each visit we walked into the forest along available trails identifying all birds heard or seen. We additionally deployed 5–8 mist nets (12 × 2 m) during each visit to capture secretive species that can be difficult to see or hear. The number of visits and surveying time was initially planned to be proportional to the area of available habitat for endemic highland species. However, we extended the surveys at each field site until no new species were recorded in at least three consecutive visits. Presence–absence of species obtained with both survey methods was combined to obtain species richness for each site. For the analyses we combined all sites in the Talamancas since we did not detect any difference in species richness among sites (data not shown) and because there are no geographic discontinuities along this mountain range that are expected to prevent bird movement (Gómez Reference GÓMEZ1986). We did not sample birds in the highlands of western Panama; however the bird list reported by Ridgely & Gwynne (Reference RIDGELY and 1989) for those highlands is identical to our list from the southernmost extreme of Costa Rica.
Figure 1. Study sites (indicated by black arrows) surveyed for this study along the four Costa Rican mountain ranges.
Statistical analysis
We assessed the relationship between area of highland islands and species richness for endemic and non-endemic species separately using the linear regression for the island biogeography model: ln S = ln c + z ln A; where S is the number of species estimated, A is the area, and c and z are fitted constants that depend on the studied system (MacArthur & Wilson Reference MACARTHUR and WILSON1967). Observed species number for each highland island was compared with species richness estimated by the model using a chi-square goodness-of-fit test. The area of each highland island (Table 1) was obtained from Barrantes & Loiselle (Reference BARRANTES and LOISELLE2002).
Table 1. Highland sites visited in each Costa Rican mountain range and respective estimated areas (ha). Areas of all Talamanca sites were combined due to absence of geographic discontinuities.
We also determined whether endemic species composition across highland islands follows a nested distribution, or if each island contains an independent subset of species of the entire highland endemic avifauna. Island assemblages are considered to be nested when species found in assemblages from smaller islands are a subset of the species found on larger islands. We assessed nestedness by measuring deviation of the observed assemblage of species from a distribution of randomly constructed assemblages. According to Patterson & Atmar (Reference PATTERSON and ATMAR1986), deviation from perfect nestedness occurs when a given species is absent from richer assemblages summed over all species, so that nested assemblages have less deviation from perfect nestedness. In this study nestedness was assessed using the NODF (nestedness based on overlap and decreasing fill) metric proposed by Almeida-Neto et al. (Reference ALMEIDA-NETO, GUIMARÃES, GUIMARÃES, LOYOLA and ULRICH2008), which improves the result of the temperature metric proposed by Atmar & Patterson (Reference ATMAR and PATTERSON1993). With this metric it is possible to decompose total nestedness (entire presence–absence matrix) into two components: the effect of sites (i.e. rows) and the effect of species (i.e. columns). The NODF metric was calculated for rows (i.e. sites) and for the entire matrix using the Vegan package (version 1.17; http://cran.r-project.org), implemented in the R Statistical Language (version 2.10; http://www.R-project.org). The statistical significance of the calculated metrics was assessed by contrasting the observed value to NODF estimates from 10 000 simulated communities. We created simulated species assemblages using null models with the r1 constraint (Wright et al. Reference WRIGHT, PATTERSON, MIKKELSON, CUTLER and ATMAR1998). The r1 model draws species for each site with probabilities based on the species’ marginal frequencies, until the site's original richness is attained. Normal probability scores (two-tailed z-scores) were used to compare observed values to the distribution of simulated temperature values. Results were comparable when other null models (r00, r2, sensu Wright et al. Reference WRIGHT, PATTERSON, MIKKELSON, CUTLER and ATMAR1998) were tested. This analysis was not applied to species that are not highland endemics, since distribution of these species is not fragmented throughout most of the country.
To infer the causality of nestedness we conducted the statistical analysis described by Lomolino (Reference LOMOLINO1996). This analysis assesses probability of two potential causal mechanisms for nestedness in insular communities: differential migration due to increased distance or barriers from a species source, and differential extinction due to a reduction in available habitat (Cutler Reference CUTLER1991, Lomolino Reference LOMOLINO1996, Patterson Reference PATTERSON1990, Patterson & Atmar Reference PATTERSON and ATMAR1986). We calculated the statistics defined by Lomolino (Reference LOMOLINO1996):%PN = 100 × (R−D)/R, where %PN = % of perfect nestedness, R = mean number of departures from random simulations and D = number of departures in the order matrix. The significance of the D statistic was obtained by comparing the calculated value to a distribution obtained from 1000 random permutations of the original matrix, using a function written in the R language (available upon request). Lomolino's statistics were calculated for two presence/absence species matrices. The first matrix ordered sites by area of available habitat with Talamanca as the largest site and Cacao volcano as the smallest (Figure 1). The second matrix ordered sites based on isolation with Talamanca as the most isolated and Barva volcano the least. Degree of isolation was defined as the vertical distance between the lower limit of the distribution of highland endemic species and the highest point of the mountain passes (or watersheds) between two adjacent highland islands. The first matrix tests the hypothesis that nestedness is caused by differential extinction, while the second considers differential migration as the causal factor of nestedness. We additionally tested nestedness causality using a different second matrix. This matrix ordered the sites by isolation caused by contraction of highland vegetation after the Pleistocene. Those mountains that first became separated as highland vegetation retreated were considered for this analysis more isolated than other mountains that became isolated later (Barrantes Reference BARRANTES2009).
RESULTS
We recorded a total of 74 bird species, 36 highland endemic and 38 non-endemic species in the highlands of Costa Rica (Appendix 1). The number of species increased for endemics and non-endemics from Cacao volcano (Guanacaste mountain range) to the Talamanca mountain range (Figure 1). The change in species number was primarily correlated with the amount of available habitat area in each locality. For highland endemic species, changes in area explained 77% of the variance in species number across localities (F 1,5 = 16.8, P = 0.009, c = −0.136, z = 0.358; Figure 2). We obtained similar results for non-endemic species with 75% of the variance in species number across localities explained by a change in area (F 1,5 = 16.8, P = 0.009, c = 1.291, z = 0.232). The slope was steeper for endemic species, though it did not differ significantly from that of non-endemics (t 10 = 0.01, P = 0.98), nor from the slope (z = 0.335) calculated by Vuilleumier & Simberloff (Reference VUILLEUMIER and SIMBERLOFF1980) for species present in 23 paramo patches in South America (t 26 = 0.001, P = 0.99, z = 0.334).
Figure 2. Relationship between area of available habitat and number of highland bird species (species number and area were ln-transformed). Number of endemic (a) and non-endemic (b) bird species sampled at different highland islands on the four mountain ranges in Costa Rica. Dotted lines show 95% confidence intervals for the regression line.
The difference between observed species number and the number predicted by the linear model (species–area relationship) varied across localities. In the Talamancas, the number of species was significantly lower than that predicted by the model (χ21 = 11.6, P = 0.0006), whereas in Barva and Poás the number of species was higher than the number predicted by the model (χ21 = 8.45, P = 0.0034; and χ21 = 10.3, P = 0.0013 respectively). At all other localities, species numbers did not deviate significantly from values predicted by the model.
The entire community of highland endemics showed a strongly nested pattern (NODF = 78.0; z = 2.25; P < 0.002). Sites also showed a similar pattern (Nsites = 95.2, z = 3.78, P < 0.0001; Figure 3), indicating that the observed distribution of endemic species across the highland localities is congruent with a nested subset of communities. In this case, the Talamanca mountain range, with the largest area of available habitat, contained all highland endemic species, whereas smaller highland islands toward the north-west contained smaller subsets of this avifauna (Table 1). Nestedness in this highland mountain system is more likely explained by differential extinction (D = 12, R = 44.5, P = 0.01,%PN = 73.0) than by differential migration due to current topographic barriers (D = 39, R = 44.3, P = 0.368,%PN = 12.0), or differential migration due to isolation caused by contraction of highland vegetation (D = 33, R = 44.2, P = 0.233,%PN = 25.4).
Figure 3. Nestedness plot with binary presence–absence of species in the ordinate and study sites in the abscissa. Species presence is shown as a shaded block, while absence is in white. The maximum incidence of species is shown on the upper left corner of the matrix and the dashed curved line is the fill line that indicates the distribution of species in a perfectly nested matrix.
DISCUSSION
Highland endemic birds present in the highlands of Costa Rica and western Panama have a stronger affinity with Andean species and to a lesser degree with Nearctic species (Barrantes Reference BARRANTES2009). This suggests that ancestral bird species probably dispersed from South and North America to Central America since the late Pliocene (Hackett Reference HACKETT1995, Haffer Reference HAFFER1974, Smith & Klicka Reference SMITH and KLICKA2010, Stiles Reference STILES and Janzen1983a, Winker & Pruett Reference WINKER and PRUETT2006). Ancestral population dispersal was likely favoured during glacial peaks, when vegetation currently typical of the highlands extended to lower elevations, but dispersal routes were later severed during interglacial periods. As climate became gradually more tropical during interglacial periods (e.g. after the Pleistocene), arrival of new species was reduced, isolating the Central American highlands from the Nearctic region and from South American mountains, promoting speciation in Costa Rican and western Panamanian mountains (Barrantes Reference BARRANTES2009, Hackett Reference HACKETT1995, Haffer Reference HAFFER1974, Winker & Pruett Reference WINKER and PRUETT2006).
The number of highland endemic and non-endemic bird species present in the Costa Rican highlands is significantly affected by the area of highland islands. Despite our small sample size, dictated by the geographic characteristics of the system studied, the effect of area on the number of species is comparable with that obtained by Vuilleumier & Simberloff (Reference VUILLEUMIER and SIMBERLOFF1980) in a much more extensive geographic area in South American highlands. They found nearly the same effect of area on the number of paramo species present in 23 patches, indicating that similar processes (e.g. extinction imposed by habitat reduction over interglacial periods during and after the Pleistocene, or differential migration due to the effects of barriers) may have similar effects on independent bird communities. However, our results should be viewed with caution, since each point (locality) in a small sample size, particularly extreme values (e.g. Cacao, Barva, Poás), may strongly affect the significance of the species–area relationship.
Within Costa Rican mountain ranges, the assemblages of highland endemic species across disjunct highland islands follow a spatially nested distribution. Endemic species found in the most depauperate isolated areas located in the north-west part of the country (e.g. Guanacaste and Tilarán, Figure 1) are present in progressively more species-rich islands. These islands have progressively larger areas towards the south-east part of the country, with the Talamanca mountain range containing all highland endemic species. This nested pattern in these assemblages of highland bird species is more likely the result of differential extinction due to a reduction of available habitat, rather than a consequence of limited migration.
Our results are also congruent with the hypothesis proposed by Barrantes (Reference BARRANTES2009) to explain the general distribution of endemic species in the highlands of Costa Rica and western Panama. This hypothesis is based on the reconstruction of highland vegetation using fossil pollen from the Dryas stadial (c. 18 000 ya; Hooghiemstra et al. Reference HOOGHIEMSTRA, CLEEF, NOLDUS and KAPPELLE1992, Islebe et al. Reference ISLEBE, HOOGHIEMSTRA and VAN DER BORG1995, Reference ISLEBE, HOOGHIEMSTRA and VAN VEER1996), and proposes that highland vegetation and highland avifauna had a continuous distribution from western Panama to the northernmost extreme mountains in Costa Rica during glacial periods. During interglacial periods, as the climate turned progressively more tropical, the distributions of highland bird species contracted and became confined to the upper elevations of mountain ranges as highland vegetation withdrew to higher altitudes where climatic conditions were more favourable. The progressive reduction in available habitat probably also led to a progressive increase in the extinction rate of the birds associated with highland vegetation. Extinction was apparently higher on isolated mountain peaks with relatively low elevation, where available habitat for endemic birds covers only a small area around the summit. Our analyses strongly suggest that changes in habitat availability produced by climatic shifts in the late (and after) Pleistocene was an important factor shaping current bird species assemblages in southern Central America.
Local extinction has been proposed to be a primary factor shaping nested patterns of species distribution in naturally isolated areas, such as archipelagos and isolated areas on mountain tops in different regions worldwide (Brown Reference BROWN1971, Martínez-Morales Reference MARTÍNEZ-MORALES2005, Mayr & Diamond Reference MAYR and DIAMOND2001, Patterson & Brown Reference PATTERSON and BROWN1991, Worthen Reference WORTHEN1996, Wright et al. Reference WRIGHT, PATTERSON, MIKKELSON, CUTLER and ATMAR1998). The evidence obtained in this study supports the hypothesis of extinction, rather than restricted migration, to explain the current pattern of nested species assemblages in the highlands of Costa Rica. This mountain system represents a unique opportunity to study the effects of isolation on morphological, genetic and behavioural (e.g. song) divergence, and which may reflect incipient speciation in a large number of species confined to small geographic areas (Barrantes & Sánchez Reference BARRANTES and SÁNCHEZ2000, Chavarría et al. Reference CHAVARRÍA-PIZARRO, GUTIÉRREZ-ESPELETA, FUCHS and BARRANTES2010). Our results are of particular interest as present global climate changes may have an impact on the vegetation structures in highland mountains and further reduce the available habitat for endemic highland bird species (Pounds et al. Reference POUNDS, FOGDEN and CAMPBELL1999).
ACKNOWLEDGEMENTS
We are grateful to Johel Chaves-Campos, Jessica Eberhard, Jeffrey Ross-Ibarra and an anonymous reviewer for constructive comments that greatly improved the manuscript, Julio Sánchez for collaboration in the field, and Federico Valverde for logistical support. Financial support was provided by Vicerrectoría de Investigación, Universidad de Costa Rica.
Appendix 1. Distribution range between Talamanca and the northernmost extreme of endemic and non-endemic bird species sampled at different mountains on the highlands of Costa Rica. Nomenclature and order of species follow the American Ornithologists’ Union (1998) and supplements. Talam, Talamanca; Mirav, Miravalles; Montev, Monteverde.
