Habitat degradation has over time formed synergy with other factors to contribute to dwindling populations of both fauna and flora by altering their habitats. The disturbance of natural habitats affects the diversity of both vertebrates and invertebrates by altering both feeding and nesting sites for which organisms are known to depend on for survival. Little is known of the extent to which vulnerable habitats could shape the diversity of most indigent pollinators such as African meliponine bee species in tropical ecosystems. This study was conducted to determine how disturbance could shape the natural occurrence of African meliponine bee species in different ecological habitats of Taita Hills, leading to changes in their diversity. A total of four species depicted by the Renyi diversity profile was recorded in five of the six main habitat types surveyed, and a further extrapolation with Shannon index (EH) also predicted the highest species richness of 4.24 in a deciduous habitat type. These meliponine bee species (Hypotrigona gribodoi, Hypotrigona ruspolii, Meliponula ferruginea (black), and Plebeina hildebrandti) were observed to be unevenly distributed across all habitats, further indicating that mixed deciduous habitat was more diverse than acacia-dominated bush lands, grasslands, and exotic forest patches. Geometric morphometrics categorized all four meliponine bee species into two major clusters—cluster 1 (H gribodoi, H ruspolii, M ferruginea (black)) and cluster 2 (P hildebrandti)—and further discriminated populations against the 4 potential habitats they are likely to persist or survive in. Each habitat appeared to consist of a cluster of subpopulations and may possibly reveal ecotypes within the four meliponine populations. This has revealed that unprecedented conversions of natural habitats to agroecosystems are a key driving factor causing increased habitat isolation and vulnerability in this Afromontane region which may potentially distort local assemblages of native pollinators, such as meliponine bee species.
As majority of landscapes in Kenya drastically change due to habitat destruction and land conversion, the composition of pollinator communities is likely altered unknowingly with the creation of pollinator deficits.1 This in turn threatens the normal provision of pollination services, which is critical for the survival of flowering plants.2,3 Pollinators are important species because plants and ecosystems largely depend on pollinators for stability due to the multiple roles they play in maintaining the viability of pollinator-dependent plants, which in turn supports herbivore and carnivore survival within a food chain.4 Meliponine bees are a group of eusocial insects that form part of this niche and they play an important role in the pollination process of plant life, particularly plants in natural and seminatural habitats.5,6 They are also considered to be crucial pollinators in tropical forests78–9 and visit more than 100 plant species in a given habitat.10
Land-use changes commonly take the form of habitat alteration, fragmentation, and isolated forest patches which are brought about by different anthropogenic activities, increasing vulnerability and negatively influencing biodiversity of flora and fauna species in habitats.11
In Kenya, the Eastern Arc Mountain is listed as one of the world’s 30 biodiversity hotspots having some of the richest concentrations of endemic plants and animals on earth.12 The eastern arc mountains runs along the Tanzanian and Kenyan coasts, which now lies between two newly classified hotspots (eastern Afromontane hotspot and the coastal forests). Most of this region is in Tanzania, which begins in the Eastern Arc Mountain and runs along the Rufiji water catchment; however, a narrow gorge close to the Kenyan-Tanzanian border follows this Eastern Arc Mountain, ending its northernmost limits in the Taita Hills of Kenya. The unique habitats of the Eastern Arc Mountains are notably fragmented, leading to rapid habitat loss with consequential effects on both flora and fauna species within key sites to become highly vulnerable to extinction. Agricultural encroachment, timber extraction, and charcoal production are listed as the greatest threats to the survival of most flora and fauna species. Taita Hills possess a high level of endemic fauna and flora,13,14 but ironically they are one of the most degraded areas in the Eastern Arc Mountains, having lost about 99% of its original cloud forest during the past 50 years.151617–18 Some plant species unique to this region include the African violet (Streptocarpus teitensis) which is restricted to a small patch in Ngangao forest, Ceropegia verticillata, Chassalia discolor taitensis, Coffea fadenii, Impatiens englerii teitensis, Impatiens teitensis, and Zimmermannia ovatoa. Also, some endangered endemic bird species include the Taita thrush (Turdus helleri), Taita apalis (Apalis thoracica fuscigularis), and Taita white eye (Zosterops poliogaster). Some other notable endemic amphibians in this hill include the Taita reed frog (Hyperolius viridiflavus) and the forest gecko (Cnemaspis dickersonii).17,19 Ecologically, habitat features are important in predicting the diversity of species and their population size as plants and animals are highly dependent on the quality of their habitats.5,20 The fragmentation of natural and seminatural habitats is regarded as a major threat to biodiversity,11 having negative effects on ecological processes such as primary productivity, population recovery from disturbance, interspecific competition, community structure, and fluxes of energy and nutrients. This can have important ecological consequences at the population, community, and ecosystem levels, and in most cases such effects are comparable in magnitude to natural disasters. However, it is not clear how strongly these apply in nature, as studies to date have been biased toward manipulations of species diversity in vulnerable habitats, and little is known about the ecological interactions of other factors, such as forest fragment size, level of fragment isolation and level of degradation, which may influence ecological processes for native bee species. Recent studies on African meliponine bee species have indicated that these bees are strongly associated with indigenous forested areas for both nesting and foraging requirements.2122–23 African meliponine bees are reported to be one of the many invertebrates mostly affected by forest degeneration caused majorly by anthropogenic activities.24,25 Recent studies on the ecology of African meliponine bee species in countries such as Uganda23,26 and Kenya27 has mentioned the importance of intact and undisturbed habitats as a key driving factor for meliponine bees to thrive, but the extent to which these groups of pollinators are distributed in vulnerable habitats caused by increasing habitat isolation in tropical regions has not been confirmed.
Materials and Methods
Study area and sampling method
This study was conducted in 2 locations, namely, the lowlands and the highlands, with 3 sites selected in each area—lowlands (Mwatate, Msau, and Mugama) and highlands (Mwachora forest, Chawia forest, and Kishenyi) in Taita Hills (Figure 1).
Taita Hills comprise 2 distinct microclimates which are found in lowlands and highlands, respectively. The lowland is mostly characterized by dry and hot climatic conditions, with Mwatate, Msau, and Mugama being grassland habitats characterized by sparsely dispersed number of indigenous tree species, whereas the highland is characterized by wet and cold climatic conditions, with Chawia forest, Mwachora forest, and Kishenyi being cloud forests characterized by mixed indigenous and exotic tree species. The study sites were chosen based on various features, such as forest fragment size, level of forest fragment isolation, forest fragment age, and level of degradation. The lowlands lie along an altitude of ~600 to 1000 m a.s.l. with severely disturbed forest fragments comprising deciduous tree species, and agroforestry is practiced on an extensive scale. The highlands lie along an altitude of ~1200 to 2200 m a.s.l with more relatively protected forest fragment patches comprising majorly of an uneven distribution of exotic tree species. Both areas are unique as they represent a mixture of indigenous and exotic vegetation which could provide potential nesting and foraging habitats for meliponine bee species. In both study locations, meliponine bee species were sampled using standardized transect walk method28 from the months of March to September 2014 (combining both the long rainy season and dry season). In each of the two study locations, nesting colonies of the 4 meliponine bee species, namely, H gribodoi, M ferruginea (black), H ruspolii and Plebeina hildebrandti, were surveyed following a successive gradient. In each study site (25 ha), 20 linear transects of 250 m × 20 m each were established using a Global Positioning System (GPS) receiver to mark coordinates with relation to habitat type. Meliponine bees were sampled using the conventional complementary method, belt transect (direct observation of nesting colonies synonymous to a visual census),29 and data such as nesting site/substrate, GPS coordinates of nest, and names of nesting trees were recorded.
Field surveys were performed during the sunny days to facilitate viewing of foraging bees exiting their colonies. Nest inspections were performed on every substrate having the likelihood of accommodating nests, such as trees, termite mounds, and the ground.30,31 Specimens from different colonies were preserved for morphological identification and genetic characterization to confirm species identity. The number of meliponine bee species and their colonies observed per transect in the different habitats were recorded.
Specimen identification by wing morphometrics and DNA barcoding
Representative specimens from each of the colonies were examined by the biosystematics unit of the International Centre of Insect Physiology and Ecology (icipe), Nairobi, Kenya, and tentatively identified as H gribodoi, H ruspolii, M ferruginea (black) and P hildebrandti based on external morphology.
Specimens comprising approximately 20 foragers were sourced as a representative from each feral colony. The right forewing of each forager was removed and placed between a 35-mm microscope glass slide and cover slip. Each individual wing was captured based on morphometric characters with a digital camera connected to a stereomicroscope.32,33 Wing images were captured and further created in JPEG format, with 1 TPS file created from the image files using tpsUtil software (version 1.49). Approximately 8 homologous points of correspondence were plotted at specified junctions of the wing venation using tpsDig2 software version,3435–36 with one single TPS file grouping each of the processed wings. The remaining collected specimens were deposited at the biometrics unit of icipe, Duduville Campus, Nairobi.
A total of 36 individuals were selected from this pool of morphologically identified specimens and their genomic DNA extracted using a guided protocol.37 The COI region was selected and used based on its demonstrated ability in resolving generic relationships within arthropods species.3839–40
Polymerase chain reaction conditions were optimized and followed with an initial denaturation step at 96°C for 2 minutes, followed by 35 cycles at 96°C of denaturation for 30 seconds. Then, an annealing cycle at 50°C for 30 seconds and elongation step at 72°C for 1 minute followed by an initial and final extension step at 72°C for 10 minutes were performed. A prestained agarose gel (1.5%) with ethidium bromide was used to visualize the polymerase chain reaction (PCR)–amplified products. A run time of 45 minutes was used to fully separate the bands and then visualized with a UV transilluminator. A total volume of 10 µL of PCR product was digested with exonuclease II and shrimp alkaline phosphatase for 15 minutes at 37°C prior to sequencing, essentially to remove any residual primers and deoxynucleotide triphosphates. Bidirectional sequencing of the PCR products was outsourced to inqaba biotec, South Africa. Specimen sequences for COI gene were aligned using the Geneious (version 8.1) software program,41,42 and an appropriate model of sequence evolution was determined using the model with the lowest information criterion. A maximum likelihood phylogenetic tree was then generated using a General Time Reversible (GTR) model in phyML and Gamma model in Mr Bayes.42,43 Assessment of branch support was done with 1000 bootstrap replicates to generate a neighbor joining (NJ) tree and estimate the confidence relations in the NJ tree.44 The comparisons of nucleotide sequences of H gribodoi, H ruspolii, M ferruginea (black) and P hildebrandti were performed by alignment with Liotrigona madecassa (accession number: HQ012823) which served as the closest related outgroup using the BLASTX (National Center for Biotechnology Information). A maximum likelihood phylogenetic tree was then generated using a GTR model in phyML and Gamma model in Mr Bayes.43,44
Analysis of variance (ANOVA) using R statistical package was used to compute the significant effect of habitat type on species abundance. A nonlinear regression model, such as the species accumulation curve, was used to estimate the number of meliponine bee species represented in the whole surveyed area.45 The species accumulation curves were used to estimate species richness and rank abundance of meliponine bee species across varying habitats types.46 Biodiversity indices (species richness, abundance, and Shannon index) were computed using the Biodiversity R package47 installed in R software. Species richness, species diversity (using Shannon index and Renyi diversity profiles) and the proportion of habitat type with most abundant meliponine bee species were computed using Renyi diversity profiles.48 Similarity index was also used to derive dendrograms that establish similarities between habitats types in terms of species composition49 (Table 1).
Summary of surveyed habitat types in Taita Hills of Kenya.
MorphoJ software (version 1.03)50 was used to create Cartesian coordinates of the 8 targeted landmarks which were then Procrustes aligned to determine existing shape variations among the different species. The data points were subjected to principal component analysis (PCA); canonical variate analysis (CVA), discriminant function analysis (DFA), Procrustes ANOVA, and regression analyses were performed to further delineate the different bee species.
After all characters were measured, comparisons between the two study sites (highlands and lowlands) were performed using ANOVA and the Tukey test for a posterior comparison among means. Differences in wing venation between the two sites by means of a contingency G test were performed, and then a PCA using a correlation matrix was performed on all log-transformed metric characters.36 Colony principal component scores (PCs) were obtained by multiplying the character coefficients by their mean value for each colony. Colony PCs from both sites were compared by means of ANOVA and were plotted orthogonally against the axes of components to obtain a comparative spatial distribution of all species within the two study sites.
Meliponine species composition
A total of four different meliponine bee species were identified on the basis of morphological differences to species level.51 All four species were recorded in all habitats sampled in the lowlands as H gribodoi was the most dominant species across the two study sites, accounting for 58.3% of all recorded nests, followed by M ferruginea (32.0%). The proportion of P hildebrandti (7.4%) and H ruspolii (2.3%) species revealed that they were the least distributed species; this indicates that all four species were highly dispersed in distribution but in uneven patterns across all sampled habitats, showing a rank abundance of (1) H gribodoi, (2) M ferruginea, (3) P hildebrandti and (4) H ruspolii, respectively (Table 2).
Rank abundance of bee species recorded in sampled locations of Taita Hills.
Meliponine species nest abundance
Varying ranges of colonies of meliponine bee species within five of the six sampled habitat types (indigenous and exotic forests, grasslands, woodlands, and bush lands) revealed an unequal distribution (Figure 2A). Species nest abundance was skewed to the left particularly in dispersed habitats of the lowlands characterized by mixed deciduous woodlands (MDW), Acacia-dominated bush lands (ADBL), and grasslands (GL), signifying a normal distribution; however, no nest was recorded in the mixed highland forest of Chawia and was excluded. The range of nest abundance was uniform in lowland habitats characterized by acacia trees, GL, and MDW compared with both indigenous cloud forests and exotic forests. These habitats—MDW, ADBL and GL—had significant numbers of unfragmented sites which explained the positive effect on nest abundance (P = .003) compared with highland habitats (indigenous mixed forests [IMF], exotic forest patches [EFP]) which had high numbers of fragmented habitats (Figure 2B), thus revealing a distinct preference between the two main habitats (highlands and lowlands).
Meliponine bee species richness and diversity
The Renyi diversity profile recorded a total richness of 4 species across the five main habitat types sampled, and a further extrapolation with Shannon index (Evenness) also predicted a total species richness of 4.24 (Table 3). The species accumulation curve climaxed at ~80 sampling points for a total species richness of 4.24 (Figure 3). The Renyi profiles indicate that habitats characterized by MDW are more diverse and provide potential nesting sites than GL, ADBL, EFP, and IMF habitats in descending order (Table 3). The diversity profiles of EFP and IMF could not be adequately ordered, as their profile curves frequently overlapped (Figure 4). At the α = 0 scale, IMF habitat overlapped EFP habitat; at the α = 1 scale (Shannon index), species diversity was ranked in sequential order: MDW > ADBL > GL > EFP > IMF; at the α = 2 scale (Simpson index), species diversity showed the same pattern for the three dispersed habitats in the lowlands. Shannon diversity extrapolation for each habitat predicted more species in MDW and GL than for other sampled habitats. Meliponine bee species were grouped according to similarity in habitat types and preferred nesting substrates and four distinct groups were distinguished according to the evenness in species distribution. Group A (EH = ∝) comprised EFP and IMF habitats, group B (EH = 1) comprised ABDL habitats, group C (EH = 1) comprised GL habitats, and group D (EH = 0.5) consisted of only MDW (Figure 5A), whereas three commonly used nesting substrates were categorized into tree (T), ground (G), and homestead (H) (Figure 5B). Approximately 140 nests which comprised the 4 species were found nesting in trees, whereas less than 27 nests which comprised three species were found in homesteads, with the least number of nests found nesting underground and comprising only one species.
Diversity indices for each habitat sampled.
DNA sequences of the COI region were edited and aligned with the program Geneious (version 8.1) software program,37 and an appropriate model of sequence evolution was determined with the least information criterion (Figure 6). Generally, the PCA and CVA recorded significant differences among all the species and among habitats. All 12 factors of eigenvalues were found to be less than 1, which accounted for 99.17% of data variability. Graphical representation of CVA (Figure 7) scores shows a clear differentiation of species within all habitats sampled, and the DFA also revealed significant differences within populations from the different habitats with values of P < .0001. In general, 99.59% of all specimens were correctly classified according to the respective habitats; with H gribodoi populations accounting for 93.96%, M ferruginea (black) for 3.57%, and H ruspolii for 2.05%, whereas P hildebrandti recorded the least DFA (0.40%).
Our study on meliponine bee assemblages in this biodiverse hotspot provides the first documentation of their natural occurrence in these vulnerable habitats, demonstrating their resilience to survive in disturbed ecosystems. This documentation is consistent with other ecological studies that reported low bee abundance and species richness with increasing agricultural intensity from a wide variety of agroecosystems.525354–55 Although studies have identified habitat loss1,56,57 arising from human activity as the key factor driving declines of native species worldwide, the synergistic effect of fragmentation on the reduction of meta-population networks was not clearly demonstrated. Our surveys showed that there was a distinct difference in the nest abundance of meliponine bee species with relation to habitat type and preferred nesting substrates, giving a clear evidence of low distribution and diversity in the highlands compared with the lowlands which may have possibly resulted from the conversion of natural habitats to agricultural dominated habitats, which is the primary form of land-use change and the largest cause of native habitat loss and fragmentation.27,58,59 This further confirms the dominance of agroecosystems worldwide to increasingly restrict bee populations existing between the interface of agricultural and natural habitats or solely within agricultural areas, as is currently observed in this hotspot. Currently, there is no consensus on how bee communities could potentially respond to biological features in these fragmented isolations as empirical studies have revealed a range of responses restricted to fragment size.22,59,60 Species richness recorded in both sites indicates a divergence from a normal pattern of diversity of meliponine bee species assemblages within this hotspot, although this is comparatively different from the species recorded in Kakamega forest in Kenya61; it unmistakably signifies the negative effects of habitat fragmentation in predicting the diversity of bee species within an ecosystem. Generally, agroecosystems that contain a mixture of seminatural habitats throughout any particular landscape can maintain significant levels of bee diversity and abundance,6263–64 even at regional scales,6566–67 as demonstrated by the majority of H gribodoi species which was more dominant and naturally occurred in all habitat types but at variable proportions, which may be attributed to its high plasticity in nesting in varied habitat types. Other studies15 affirm that land-use intensity and proximity to seminatural habitats best explained bee species richness across landscapes, but loss of bee species richness was not solely the result of declines within such habitats but also increased homogenization of plant community composition between and within habitats, which could be a contributing factor acting in synergy with land-use intensity.
Our surveys have indicated that habitats characterized by mixed decidious woodland/acacia tree–dominated habitats presents itself as a much preferred habitat for nesting as profile curves indicated that more species could be identified with increased sampling sites in such unique habitats and on more tree species as preferred nesting hosts. This is confirmed by similar studies29,68,69 on bee pollinators and plant pollinator interactions in fragmented landscapes. The 4 species recorded directly from sampling are close to the JEvenness extrapolated predicted value of 4.24. The species accumulation curve indicated ~80 sampling points as adequate to recover at least four species in such vulnerable habitats.
We further demonstrated that geometric morphometric analyses could successfully segregate all the four meliponine bee species, distinctly grouping them into two clusters—cluster 1 (H gribodoi, H ruspolii, and M ferruginea (black)) and cluster 2 (P hildebrandti)—and successfully discriminating populations against four different habitats in Taita Hills. Each habitat appeared to consist of cluster of subpopulations and may possibly reveal ecotypes within the 4 meliponine populations; studies conducted by Owen70 also demonstrated that this tool can also been able to fully discriminate bumble bee species Bombus terrestris into possible ecotypes.
A major reason for this clustering of species would be the superficial resemblance of the 3 species belonging to cluster 1 (H gribodoi, H ruspolii, and M.ferruginea (black)) with regard to morphological similarities in forewing characters (open submarginal cells, anterior region of the submarginal cross-vein faintly visible, and nondistinct veins) and cluster 2 (P hildebrandti) which has distinct marginal cells, closed submarginal cells, and distinct veins. Also, characteristic vegetation types and climatic conditions each habitat appeared to have may have ultimately altered morphological characters for greater survival in such habitats.
The results of a PCA on the morphological measurements corroborated with molecular results, revealing the specimens clustering in four different clades (H gribodoi, H ruspolii, M ferruginea (black), and P hildebrandti, respectively. This conclusively shows that integrating DNA barcoding with morphometrics can help in segregating species that have high levels of similarities, ie, Hypotrigona spp. As there is a distinct lack of data on how pollinator communities disassemble, predictions arising from the recent proliferation of simulation studies based on networks of pollinator interactions with their nesting habitats are a valuable source of future testable hypotheses.
In summary, our study provides the first documentation of meliponine bee assemblages within such vulnerable habitats of this biodiverse hotspot and further reveals higher species diversity in certain habitats characterized with deciduous tree species that are indigenous to such habitat. Similar trends with respect to habitat composition revealed higher variation within sites that had more density with tree species indigenous to these habitats than continuous forested landscapes of the highlands mostly dominated by exotic tree species, implying higher heterogeneity within indigenous vegetation of the lowlands than in exotic forested landscapes of the highlands. This further confirms the immediate need to conserve tree species that are indigenous to these habitats, due to the high tendencies of logging them for charcoal production. Conservation of natural landscapes could prevent environmental threats to the existence of these bee species and reduce the gradual extinction of these tree species which are used as preferred nesting substrates in such vulnerable habitats.
The authors would like to thank OWSD/TWAS for providing B.O.B with a postgraduate fellowship, Behavioural and Chemical Ecology Unit (BCEU), icipe, and for providing professional and logistic support; Daisy Salifu of the statistics department for guidance on appropriate statistical tools; Josephat Bukhebi of biosystematics unit of icipe for identification of the species; and finally, the University of Nairobi. The authors are particularly grateful to Mzee Mwakolomba of Tsavo Ecotourism Association, Taita Hills, for facilitating field work and logistics.
 Six peer reviewers contributed to the peer review report. Reviewers’ reports totaled 1241 words, excluding any confidential comments to the academic editor.
 Financial disclosure The author(s) received no financial support for the research, authorship, and/or publication of this article.
 Conflicts of interest The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
 BOB and FA conceived and designed the experiments and contributed to the writing of the manuscript; BOB analyzed the data and wrote the first draft of the manuscript; BOB, FA, PNN, and LI agree with manuscript results and conclusions, jointly developed the structure and arguments for the paper, and made critical revisions and approved the final version. All authors reviewed and approved the final manuscript.
 As a requirement of publication, author(s) have provided to the publisher signed confirmation of compliance with legal and ethical obligations including, but not limited to, the following: authorship and contributorship, conflicts of interest, privacy and confidentiality, and (where applicable) protection of human and animal research subjects. The authors have read and confirmed their agreement with the ICMJE authorship and conflict of interest criteria. The authors have also confirmed that this article is unique and not under consideration or published in any other publication, and that they have permission from rights holders to reproduce any copyrighted material. Any disclosures are made in this section.