Contrasting Patterns of Species Richness and Functional Diversity in Bird Communities of East African Cloud Forest Fragments

Rapid fragmentation and degradation of large undisturbed habitats constitute major threats to biodiversity. Several studies have shown that populations in small and highly isolated habitat patches are prone to strong environmental and demographic stochasticity and increased risk of extinction. Based on community assembly theory, we predict recent rapid forest fragmentation to cause a decline in species and functional guild richness of forest birds combined with a high species turnover among habitat patches, and well defined dominance structures, if competition is the major driver of community assembly. To test these predictions, we analysed species co-occurrence, nestedness, and competitive strength to infer effects of interspecific competition, habitat structure, and species′ traits on the assembly of bird species communities from 12 cloud forest fragments in southern Kenya. Our results do not point to a single ecological driver of variation in species composition. Interspecific competition does not appear to be a major driver of species segregation in small forest patches, while its relative importance appears to be higher in larger ones, which may be indicative for a generic shift from competition-dominated to colonisation-driven community structure with decreasing fragment size. Functional trait diversity was independent of fragment size after controlling for species richness. As fragmentation effects vary among feeding guilds and habitat generalists, in particular, tend to decline in low quality forest patches, we plead for taking species ecology fully into account when predicting tropical community responses to habitat change.


Introduction
Habitat fragmentation has profound and mainly negative effects on the long-term viability of indigenous animal and plant species [1][2], in particular in historically stable ecosystems such as tropical rainforests [3][4][5]. While habitat fragmentation mainly causes a decline in species richness at the regional level [6][7] (but see counterexamples in Schmiegelow et al. 1997 and Eastern Arcs and part of the Eastern Afromontane biodiversity hotspot, has been subject to rapid loss, degradation and fragmentation of pristine habitats over the past decades (further details see material and methods) [33]. Despite this transition, indigenous forest remnants still harbour a typical cloud forest avifauna, including many endemic and endangered forest habitat specialists, but also a large number of habitat generalists that also occur in the non-indigenous landscape matrix.
The Taita bird community has been studied intensively over the past two decades, resulting in well supported knowledge about species richness, abundance, and ecological demands. Furthermore, the land use history of the Taita Hills is very well documented, and this combined information offers a strong framework to study how tropical avian communities are shaped in relation to species and landscape traits. Making use of this information, we here test the following three hypothesis: (i) Avian species density increases with decreasing fragment size, resulting in a high proportion of species surviving in forest fragments relative to intact forests; (ii) Species co-occurrence among small forest remnants is aggregated due to crowding effects; and (iii) Habitat specialists respond more strongly to habitat degradation compared to habitat generalists.

Taita Hills study region
The Taita Hills cover an area of around 250 km 2 and are geographically isolated from other mountain blocks to the south (90 km to the Usambara Mts.) and the north (80 km to the Chyulu Hills) (S1 Appendix). Semiarid plains in either direction constitute a strong dispersal barrier for species that depend on moist and cool cloud forest habitat, and this resulted in high levels of endemicity [34][35]. Degradation and fragmentation of the Taita forests started long before the colonial era, when slopes were cleared for agriculture up to the head of the streams [36]. Large-scale forest loss occurred during railway constructions between 1898 and 1924, while in more recent times, forest cover markedly decreased between 1955 and 2004. Even though half of the original indigenous forest has currently been lost, airborne remote sensing of spatio-temporal changes in forest cover [33] revealed that the total forest cover in the Taita landscape remained about the same between 1955 and 2004, mainly due to planting of exotic trees on rocky, barren or eroded areas, secondary bushlands and abandoned agricultural land. In addition to indigenous forest loss, the remaining patches also decreased in forest quality due to pit-sawing, charcoal manufacturing, firewood collection, pole removal and grazing [33]. Three larger forest fragments (Chawia (90 ha), Ngangao (147 ha) and Mbololo (179 ha)), nine smaller ones (< 15 ha), and several tiny patches of indigenous forest remained embedded in a fine-grained mosaic of human settlements and small-holder cultivation plots [33]. Small forest fragments, in particular, continue to suffer from cattle grazing and other forms of habitat disturbance, while the three larger forest fragments vary in the degree of habitat degradation too, being highest in Chawia forest, intermediate in Ngangao forest, and lowest in Mbololo forest [37].
Land-cover information was derived from airborne true-colour images, converted to orthomosaics at a spatial resolution of 0.5 m [33]. Brightness variations were removed by corrections for light falloff and bidirectional effects using the methods developed by Pellikka (1998) [38] after which frames were mosaicked using the EnsoMOSAIC [39]. The resulting mosaics were orthorectified, projected to Transverse Mercator projection with a Clarke 1880 [40] spheroid and Arc 1960 datum, and resampled to 0.5 m ground resolution. The resulting geometric accuracy was within 2 m as verified in the field using GPS. The land cover model was subsequently ground-truthed, revised and fine-tuned during field visits in 2007 and 2008, confirming the correct remote-sensing classification of large patches of closed-canopy forest, exotic plantations, and non-forested habitat [33]. Based on this land cover model, we calculated the following four landscape characteristics using Fragstats v. 3.3 and ArcView 3.2 (ESRI 2013): (i) indigenous forest patch size, (ii) indigenous forest patch perimeter; (iii) percentage of closedcanopy forest cover within 800 m of each indigenous forest patch, and (iv) patch proximity, a distance-weighted, area-based isolation index (PPI) [40]. We related observed species richness to these landscape characteristics to test the first and third hypothesis.

Bird assessments
Understorey bird community metrics were derived from a long-term (1996 to 2010) bird ringing program using standard mist-netting procedures as described in Karr (1981) [41]. Mist netting was conducted in collaboration with the Ornithology Section of the National Museums of Kenya, Ornithology Section. Permission for bird collection was issued by the National Museums of Kenya. Permits to access the forest fragments were provided by the Kenyan Forest Service. As endangered, Taita endemic bird species were involved in this study, its collection was approved and mainly conducted by members of the ethics committee of the National Museums of Kenya personally (P Njoroge, RK Mulwa, O Kioko). Birds were collected using mist-nets, which were regularly controlled, to prevent any negative effects on trapped birds. This activity was approved by the animal ethics committee of the National Museums of Kenya. Mist-net lines were operated in one to seven plots per fragment (depending on fragment size) and were evenly spaced out to sample entire plots, while net positions, net lengths (120m/plot) and daily trapping efforts (06-18h) were kept constant between trapping sessions. Nets were routinely checked at 30-minute intervals so as to promptly remove, process, and release the birds. Time intervals between subsequent ringing sessions varied between 1.0 and 4.6 months, and the number of ringing sessions per fragment ranged between 20 and 32 over the 15 year study period. While mist nets are regarded as likely the best technique for assessing the relative abundances of tropical understorey birds [41,42], habitat modifications such as removal of canopy trees and clearing of the understorey, in particular, may alter flight height of some species, thereby changing their susceptibility to mist-net capture [43]. To minimize this possible bias in the assessment of species richness, we restricted our analysis to the understorey bird community, i.e. species that are reliably caught in mist nets. Therefore, our data, covering 15 years of bird observation, are believed to be highly appropriate to assess total species richness and as well as relative abundances in the smaller fragments sampled with identical sampling effort.
For each bird species we assembled data on body mass (g), average bill culmen length, depth and width (cm), tarsus, tail, and wing length (cm), and average hand-wing index. The dominant principal component of the three bill measures was used to assess bill characteristics, and the respective dominant eigenvector of the wing, tarsus, and tail measures served as a proxy to the type of locomotion [44]. Following Claramunt et al. (2012) [45], we used the hand-wing index to quantify dispersal ability. Each species was assigned to one of four feeding guilds, insectivore, seed-eater, fruit-nectar feeder, and omnivore. We also assembled dietary (coded into ten categories) and foraging stratum preferences (coded into eight categories), respectively [46]. To reduce dimensionality, we calculated the dominant principal components of the matrices and used these in subsequent analyses. The complete species list of bird species, respective abundances per forest fragment, and species traits are given in the S1 Appendix.

Statistical analysis
Analyses were based on matrices containing species relative abundances with species in rows and forest fragments in columns. We composed matrices for each fragment separately and then grouped the three larger (> 90 ha) and nine smaller (<15 ha) fragments in two additional matrices. We calculated the following three metrics of species co-occurrence. First, we estimated species segregation among fragments (negative species associations) by the Cscore [24,47,48] that is a normalized count of the number of checkerboard submatrices ({{1,0},{0,1}} or {{0,1},{1,0}}). As an auxiliary metric of species segregation that focuses on the pattern of species turnover among sites we used the standard proportional turnover beta P ¼ 1 À alpha gamma beta P ¼ 1 À alpha gamma ; where alpha refers to the average richness per fragment site and gamma tom the total observed species richness [49]. Third, we performed a nestedness analysis to identify gradients in species occurrences and richness across fragments [23] using the NODF (nestedness by overlap and decreasing fill) metric of Almeida-Neto et al. (2008) [48]. For the nestedness analysis, rows were always sorted according to species occurrence totals. Finally, we calculated the functional diversity of each feeding guild based on five measured functional traits (bill characteristics, body size, dispersal, locomotion, stratum) using the functional attribute metric FAD [49,50], a measure of total trait space encompassed by the species of a given community calculated as the sum of the Euclidean distances between species in trait space. For comparability, trait expressions were Z-transformed prior to calculation. Co-occurrence analyses were done using the freely available Fortran application NODF [51] and Turnover [52]. Source code is available from WU by request.
For statistical inference of NODF, C-score, and beta P , we used a null model approach and compared the observed co-occurrence metric scores with those obtained from 1000 null matrices each that were randomized using a null model that resamples the matrix with placement probabilities proportional to observed total abundances of rows and columns (proportional abundance model [53]). This is a conservative null model that has been shown to account well for inherent site differences and unequal species colonization probabilities (the mass effect) that are not directly linked to the pattern of interest [24]. Note that this model is equivalent to a neutral model without dispersal limitation and speciation. To account for richness effects [54], raw FAD scores were compared to a null model in which the trait expressions for each single trait were randomly reshuffled among species. Statistical significance was estimated from the respective tail distributions at the two-sided 5% error level. Additionally, we calculated standardized effect sizes (SES = Obs-Exp) / StDev Exp ; Obs and Exp: observed and expected scores, StDev Exp : standard deviation of expectation). SES scores should have values below -1.96 and above +1.96 at the two-sided 5% error level under the assumption that the respective null distribution is approximately normal. To account for multiple testing, all significance levels were Bonferroni corrected.
Possible competitive interactions among species within the seven feeding guilds were assessed following Ulrich et al. (2014) [31] and Soliveres et al. (2015) [32]. For each guild, we calculated 100,000 random species × species competitive strengths matrices, translated these into a column stochastic transition matrix, and used a Markov chain model to predict relative species abundances from this transition matrix within the 12 fragments. We compared predicted and observed relative species abundances by rank order correlation (r C ) and chose the best fitting competition matrix to assess the maximum impact of interspecific competition on community assembly. High values of r C point to the possibility that interspecific competition is a major driver for observed species distributions while low r C values point towards a minor impact of competition [55].
Using the empirical species-area relationship S = 12.4A 0.21 ;r 2 = 0.44 (Fig 1a) we assessed average species densities S A per ha area by S A = S/A 0.21 , with S reflecting the total species richness and A the fragment area. We used a general linear models with identity link function and normal error structure to relate FAD to the categorical variables feeding guild, functional traits, and remnant size group (large-small), and to species richness as the continuous covariate.

Results
We recorded a total of 17,520 individuals from 69 bird species in the 12 Taita forest fragments (Table 1), of which 36 species belonged to the insectivorous guild and 12 species to the frugi/ nectarivorous guild (Table 2). Species richness increased moderately with area (Fig 1a) and perimeter (r = 0.61, P = 0.03). The best predictor of fragment richness was the total number of individuals caught (Fig 1b). Richness did not significantly vary with fragment isolation (Pearson r = 0.27, P > 0.3) and increased weakly with forest cover within the matrix (Pearson r = 0.51, P = 0.09). In fragments below 15 ha, species richness was independent of fragment  area (Fig 1a, r = 0.17, P > 0.5) but tended to be positively related to fragment perimeter, albeit not statistically significant (r = 0.61, P = 0.08). Most species rich were the Ngangao (39 species) and Chawia (38 species) forest fragments, while the large and predominately pristine Mbololo fragment was comparatively poor in species richness (32 species). Of the smaller fragments, Fururu and Macha were most species rich (30 and 32 species, respectively). Average species density was independent of fragment isolation (Fig 2a), forest cover (Fig  2b), and fragment perimeter (Fig 2c). In the smaller fragments, species density increased with total abundance (Fig 2d). In the larger fragments, species density was at an intermediate level compared to the smaller fragments (Fig 2d). A comparison of the relative abundances between the smaller and the larger fragments revealed a significant shift in relative abundance between the two fragment types. Of the 22 species with relative abundances below 0.001 in the larger fragments, 18 achieved higher relative abundances in the smaller fragments (not shown). This shift was accompanied by a sharp decline of five species (Columba larvata, Phyllastrephus placidus, Phylloscopus ruficapillus, Turdus olivaceus, Zoothera gurneyi) in the smaller fragments.
A separate analysis based on forest specialist species only (Table 1) also revealed an increase in richness with fragment size (Fig 1a) and abundance (Fig 1b), however, less strong so compared to an analysis with all species included (Fig 1). The proportion of forest specialists decreased with fragment area (r = -0.32, P = 0.32) and abundance (r = -0.73, P < 0.01), while specialist species density was not significantly related to habitat isolation (Fig 2a), forest cover (Fig 2b), fragment perimeter (Fig 2c), or specialist abundances (Fig 2d).
A comparison of species richness between the three larger (total area 416 ha, 56 species) and nine smaller fragments (total area 46 ha, 55 species) ( Table 1) revealed a loss of 14 species (25%) and a gain of 13 species (23.6%) in the latter. The larger fragments contained 17 species of forest specialists, the smaller ones 14 species. Despite the fact that the nine smaller fragments comprised only 11% of the area of the larger fragments, total species richness decreased by 2% (1 species) only. The prevalence of negative nestedness and positive beta P SES scores (Table 2) further indicated spatial turnover in species composition among the larger and smaller fragments, however, we did not find direct evidence that the latter was caused by competitive interactions. Indeed, correlations between the competitive interaction matrix and the observed Table 2 distribution of abundances (Table 2) were on average weak and explained at most 50% of variance in abundance. Finally, we compared functional diversity between feeding guilds and between larger and smaller fragments (Fig 3a, Table 3). Except for omnivores, functional attribute diversity significantly differed between feeding guilds and was always larger in the large fragments (Fig  3a). We observed the same pattern for forest specialists (Fig 3b) although the differences were statistically not significant due to the small number of species. Linear modelling   Table 3. General linear mixed modelling detected significant (*: parametric P(F) < 0.05; ***: P < 0.001) differences of functional attribute diversity (FAD) among feeding guilds, and guild×functional trait combinations when using raw FAD as dependent variable, but not when using SES transformed values (traits reshuffling null model accounting for differences in species richness between fragments (Table 3) revealed that these differences were mainly caused by co-variation with species richness. Although FAD significantly differed between feeding guilds and guild specific differences in species traits (Table 2), FAD did not vary with fragment size (Table 3). FAD also did not significantly vary with any of the categorical variables (guild, trait, size) when using standardized effect scores (Table 3).

Discussion
Within an island biogeographic framework [56], species richness is predicted to vary positively with patch size and negatively with patch isolation. Along these lines, Canale et al. (2012) [57] reported a strong decline in species richness in tropical forest fragments of less than 10 ha (equivalent to a lowered species density), approximately the upper size limit of our small forest fragments. Hanski et al. (2013) [58] predicted such a decline while demonstrating that SARs that only account for area effects tend to overestimate species richness in fragmented landscapes if fragments are highly isolated, and they proposed an extension of the power function SAR that downsizes species richness in small fragments. In turn, a number of empirical studies reported a small island effect [12] where species density in very small islands or habitat patches increased [59][60]. Our results do not support either prediction. Species richness in the smallest fragments did not negatively deviate from the observed SAR (Fig 1a) although we observed a tendency to independence of area below 10 ha. However, richness was very closely related to the number of individuals indicating that habitat capacity is probably the main trigger of species richness. Contrary to Hanski et al. (2013) [58], area corrected richness (species density) was not linked to fragment isolation (Fig 2a).
The counter-intuitive finding that species richness in small fragments was only 2% lower than in large ones, is in line with other studies that reported patch connectivity to be of higher importance than patch size [61], in particular for generalist species with a broad ecological amplitude that can easily cross the landscape matrix and (re)colonise small forest patches [62,63]. This generalist-focused explanation is corroborated by the fact that specialist species increased less in richness in the larger fragments (Fig 1). Bird mobility in heterogeneous landscapes tends to vary among feeding guilds [64][65][66], with small understorey insectivores often showing the highest sensitivity to fragment isolation [67]. The resulting high species turnover (i.e. partial replacement of sedentary, forest-restricted specialists by mobile, matrix-tolerant generalists in small, degraded forest fragments) might be responsible for the weak SAR observed in our study.
Previous studies showed that bird assemblages strongly vary in composition with habitat area and structure (such as the degree of fragmentation), and predicted higher species richness in large, connected forest patches [68][69] [62] argued that loss of tropical species is mainly driven by habitat destruction, rather than fragmentation. According to these authors, species with broad ecological amplitude can readily exploit the landscape matrix in which forest fragments are embedded and may even gain from habitat fragmentation. The high species turnover between isolated forest fragments (independent of fragment size) observed in our study is in line with these predictions.
Various studies have reported increasing mobility with decreasing habitat integrity [72][73][74]. Yet, Price (2006) [75] showed that frugivorous species conducted shorter, rather than longer, movements in fragmented forests, despite their high intrinsic mobility. Based on earlier studies of mobility [76][77] and gene flow [78] in the same study area, bird species currently surviving in the Taita Hills forest also appear to differ in metapopulation dynamics. Results of this study add to the effect of functional connectivity [79] and support the need to define species-(or guild-) specific fragmentation thresholds [80] in conservation. Indeed, even though we found species richness to be only marginally affected by patch area, there are documented cases of local species extinctions in small Taita forest fragments [42]. In this respect, our results are also relevant for the ongoing SLOSS debate [81][82]. Although each large fragment was at least twice as large as all smaller fragments together, they hosted at most 75% of the total species richness of the small fragments ( Table 2). As such, our findings corroborate the theoretical expectation of Lasky & Keitt (2013) [26] that networks of small habitat remnants would be able to maintain higher total species richness than homogenous habitat blocks of the same size. Yet, contrary to these authors, we did not find reduced fragment (alpha) diversities within each of the smaller fragments. We interpret these findings as evidence that tropical forest birds respond to increased habitat fragmentation by higher mobility [78], or that there is a debt between former habitat destruction and ongoing local extinction [83].
We were surprised to see that functional trait identity was not linked to the pattern of species-co-occurrence. Species interactions are mediated by traits and classical community assembly theory predicts co-existing species to show low levels of trait similarity [84]. Yet, we did not find strong evidence that the bird communities in our study fragments are shaped by competitive interactions (Table 2), which might explain why trait expression was not significantly linked to co-occurrence. This finding contrasts to recent evidence by Bregman et al. (2015) [16] who reported increased levels of interspecific competition within bird guilds at decreasing fragment size. However, these authors used indirect evidence within the community assembly framework that assumes that overdispersion of phylogenetic and functional traits is a consequence of competitive interactions (Darwin's competition-relatedness hypothesis, reviewed in Allan et al. 2013, Götzenberger et al. 2012) [85][86]. However, Cahill et al. (2008) [87] found little evidence for this assertion and it is now well known that overdispersion (spatial segregation) might stem from different processes, particularly from filter effects within heterogeneous landscapes [88] and even from dispersal limited neutral community assembly [89]. Here we used a direct way to assess whether any set of competitive strength relationships between species is able to predict observed abundance distributions. For the smaller fragments this relation was highest in omnivores where competition explained at most 13% of variance in species relative abundances (r C = 0.36, Table 2) comparing to 50% in the larger fragments. Apparently species density in the smaller fragments yet did not reach the threshold for intense competitive effects and thus is of rather marginal importance for community assembly.
Total trait space, as expressed by FAD, increases with species richness [52]. However, various authors found homogenization effects in fragmented landscapes where smaller fragments are devoid of habitat specialists and regionally rare species [90][91][92]. Consequently, this selective pattern of species extinction should translate into a reduced effective functional diversity that is the degree of FAD after accounting for richness effects. Rather than observing such a pattern, FAD remained constant after correction for richness differences (Table 3). While small tropical forest fragments are hence able to maintain a high effective functional diversity, communities will inevitably collapse when fragment areas become exceedingly small, which raises the question about a minimal tropical forest fragment size. SES-transformed FAD scores of the three smallest fragments (Mwachora, Kichuchenyi, Wundanyi, Table 1) were indeed smaller (average SES FAD = -0.26±0.14) than those of the three large fragments (average SES FAD = 0.13±0.15), although not statistically significant. A one ha area is probably at the lower boundary for a functioning understorey bird community in the Taita forest archipelago.

Conclusion
In conclusion, results of this study do not point to a single ecological driver of the observed variation in avian species composition among indigenous fragments of the Taita forest archipelago. Interspecific competition does not appear to be a major driver of species segregation, since in small forest fragments no single competitive strength hierarchy was able to predict at least a major part of the observed species abundances distributions. In larger fragments, however, the relative importance of competition might be higher. We did not find strong evidence for habitat filtering either, and community-wide species co-occurrences and joint occurrences between pairs of species were neither nested nor segregated, as would be expected if competition would be the main driver of species occurrences. The frequency of pairwise species segregation was even much below the level expected under random association. However, the large differences in species density and species richness between the two smallest fragments highlight that local peculiarities might heavily constrain species richness, irrespective of subsequent patterns of co-occurrences. Such variable, and partly opposing, responses of single guilds to fragmentation suggest that it is vital to take species ecology into consideration when predicting community-wide responses to habitat change in tropical forests.