Long‐term (1979–2019) dynamics of protected orchid bees in Panama

Plants and pollinators are linked but their dynamics are scarcely known. Chemical monitoring of male “orchid bees” at two sites revealed 75% of species were stable or increasing. Forest bees of 33 species, with live sighting at Pipeline Road (PR), and trapping on Barro Colorado Island (BCI), included 132,000 individuals. No species or community changes occurred in strong El Niño‐Southern Oscillation (ENSO) climate events, which lasted 145 total months during more than 70 bee generations. Parasite and host bees fluctuated in unison. A few very common species, adapted to relatively disturbed habitat, diminished over time at PR, while sightings and traps revealed stable abundance ranks in Euglossa 40 years on BCI. Bee abundance and biomass were stable but 50% of species had few records. Orchid bees appear more stable in older forest, they were evidently insulated from ENSO disturbance, and probably benefit from the increasing abundance of flowering lianas and vines.


| INTRODUCTION
Our understanding of long-term population dynamics relies on data from designated census points over several decades (Didham et al., 2020;Terborgh, 1989). In large nature preserves, a consensus from such data corroborates reasons for conserving forest and wildland residents, including mutualists like pollinators (Nichols & Williams, 2006;Tepedino & Portman, 2021). We asked whether a protected mature forest environment maintains its pollinators. Despite a plethora of discussion, there is almost no information on native bees in this setting (Dicks et al., 2021;Herrera, 2019;Murray, Kuhlmann, & Potts, 2009;Winfree, Griswold, & Kremen, 2007). The consequences of greater or lesser seed production on plant fitness (e.g., Xiao et al., 2017) is also a significant but neglected correlate of pollinator service in a natural habitat (Thomson, 2019) and remains unaddressed here. Using rigorous and continuous data collection, still rare in this line of research (Bonebrake, Christensen, Boggs, & Ehrlich, 2010) and unavailable for this length of time in other studies on pollinators or tropical bees that we are aware of, we focus on pollinator availability and fluctuation as we assess the impact of major cyclic climate events on populations.
Although many species are perennial in tropical forests, large seasonal and yearly changes in insect species abundance occur, and they appear similar to temperate zone patterns (Ackerman, 1983;Wolda, 1978). Available long-term studies on bees or pollinators (Cane, 2021;Soroye, Newbold, & Kerr, 2020;Van Klink et al., 2021) seldom include or analyze global climatic disturbances that influence temperature and drought, or heavy rains-notably the "El Niño-Southern Oscillation" (ENSO). Simulation study often neglects to examine dynamics in either pollinators or plants (Imbach et al., 2017). That general approach, due to a lack of long-term data, must replicate contemporary samples in space (Faleiro, Nemésio, & Loyola, 2018;Nemésio & Silveira, 2006), rather than over time. For the latter, a few investigators observe that natural population trajectories can reverse or change significantly over multiple decades, but remain unexplained (Condit, Pérez, Lao, Aguilar, & Hubbell, 2017;Freed & Cann, 2010;Macgregor, Williams, Bell, & Thomas, 2020;Thomson, 2019). Because pollinator declines are appreciated almost exclusively using historic records (e.g., collection or field sighting) in anthropogenic or disturbed settings (Burkle, Marlin, & Knight, 2013;Cane, 2021;Gegear, Heath, & Ryder, 2021;Zattara & Aizen, 2021), we see the necessity of questioning the conservation and management outlook for large reserves. We tested the impact of nine strong to very strong ENSO climatic events during 1979-2019 using forest bees in Panama. The resident "orchid bee" species, very long-tongued mainly solitary bees with no food stores or management of any kind, were studied throughout this period. The males come quickly to chemical scent baits and do not dwell in established nests but instead shelter in the wild (Roubik & Hanson, 2004). We also compared two census methods simultaneously, and related our data to a pioneer study, in the same location, of 1979. If our data confirm pollinator or bee decline indicated in other studies, then a radical departure in the expectation for forest preserves, their management, and the activity of pollinators is underscored. Instead, we found interesting changes in species abundances, much stability, and some increases, which we examine here.

| Study sites and field methods
We investigated a lowland tropical pollinator community during four decades in protected forests, one a 16 km 2 island and one close by on the mainland. Our field methods, established in 1980 (Roubik & Ackerman, 1987) have been widely employed for bee studies in tropical America (Ackerman, 1983;Roubik, 2004). Our study area's core forest is approximately 15 Â 50 km (GoogleEarth, 2020 image), of which 280 km 2 is protected in Barro Colorado Nature Monument and nearby Soberanía National Park (Supporting Information; Methods).
The studies were along Pipeline Road (PR) and on trails of Barro Colorado Island (BCI). The study area, classified as moist seasonally dry forest with 2,600 mm annual rainfall, 9 N latitude, 80 W longitude, has been protected since 1914, when the Panama Canal was first operating. It contains secondary and old forest types. Our bee monitoring methods were developed in a year-long study of wild euglossine "orchid bees" on BCI during 1979-1980(Ackerman, 1983, using 16 chemical attractants and nondestructive live sighting. We used the three most effective baits at PR and continued in monthly censuses 1979-2018. Using the permanent forest study plot (Condit et al., 2017) near the hilltop center of BCI, our additional monitoring 2009-2019 used 10 McPhail traps (BioBest Group NV: Belgium), placed 1.3 m aboveground, at 300 m intervals along five trails. Traps were employed for 1 week, baited with 7 ml reagent grade 1-8 cineole (Merck: Germany) mixed in 100 ml of vehicle radiator coolant, then bees were taken to the laboratory for identification (taxonomic guides in Supporting Information). Like the live-sighting method, the positions of attractant baits were fixed during our study. Bees were generally monitored monthly or quarterly.

| Community statistics
Our analyses of abundance among all but the rarely sampled bee species, and trends in long-term data that combined all species, used nonparametric statistics, leastsquares regression, Pearson correlations and multivariate methods developed from vegetation analysis. Rank clock plots of dominant species, species richness, species turnover, mean rank shifts and rate of community change were also evaluated (Supporting Information; Methods).

| ENSO climate anomaly impact
Detailed ENSO weather data (Oceanic Niño Index [ONI], www.cpc.ncpe.noaa.gov) allowed analysis of 9 total El Niño or La Niña "strong to very strong" intervals, which spanned 22 years and included 145 months, 58 of which we had data from at PR ( Table 2, Table S1). Relatively warm or cool years in the tropical Pacific can result in substantially greater or less than normal precipitation and T A B L E 1 Total counts of male bees using 16 chemical baits on BCI (1979BCI ( -1980, at Pipeline Road, Soberania National Park (PR, 1979(PR, -2018 using three maximally effective baits (live counts), and at BCI (2009-2019) using the most attractive chemical, in 10 replicate traps; Eg. = Euglossa, Ef. = Eufriesea, El. = Eulaema, Ex. = Exaerete; abundances deemed common, uncommon or rare (with no species analyses here) follow natural breaks in abundance, i.e. for BCI at dressleri and nigrita, then at pulchra and villosiventris, at PR, championi and flammea. BCI 2009BCI 1979-1980PR 1979-2018 Genus

| RESULTS
Data from our medium-to long-term sampling of the 4 native orchid bee genera (Euglossa, Eulaema, Eufriesea, and the parasitic bee Exaerete) comprise 257 time-series points, each with 8-28 bee species included (means = 22 and 13, respectively, PR and BCI). There were 596 months in the combined studies; therefore, slightly less than half the total months spanning our studies had data. An average sample on PR was 342 individuals (range 53-923), and the mean of total bees trapped at BCI, using 10 traps, was 1,093 (range 114-2,750). We identified all, including 44 species and 132,000 individuals (Table 1). The same species of bees at PR were previously censused for a year using 16 attractants on BCI and as expected, annual bee abundance peaks were followed by pronounced lows in late wet season. Our trap methodology was ineffective; however, censusing the largerbodied bees, Exaerete, Eufriesea, and Eulaema, because it was designed to prevent too many bees from entering (see Section 4). Total abundance and biomass were stable, both at PR and BCI (Figures 1 and 2), with about 75% of all species stable or increasing (Figures S1-S3). The 19 BCI Euglossa registered the same relative abundance rankings in 1979-1980 and 2009-2019 (Wilcoxon z = À1.489, p = .136; Table 1). The most common Eulaema species and Eg. imperialis (the largest Euglossa) were among the top-10 bees in abundance both on BCI (1979-1980, top 3, 4, and 5 species) and at PR (top 1, 5, and 6 species). Seasonality in species abundance varied little ( Figure S4). Even though no difference was found in the two-ranked species abundance lists that we compared, separated by 30-40 years in 19 Euglossa at BCI, this was not the case for the PR-more recent BCI comparison (z = À2.296, p = .021) or for the orchid bees of PR (1979PR ( -2018 and BCI (1979BCI ( -1980, z = À4.443, p < .0001. In the long term, PR displayed the positive directional change of an unstable community (R 2 = .009, F 1,701 = 6.45, p = .011), but did not in the shorter-term at BCI (R 2 = .014, F 1,43 = 0.60, p = .44, Figure S5). The directional change may be interpreted as the decline in Eg. imperialis, the most common species-the "Imperial bee", Eg. imperialis, which totaled 26-48% of all bees in all three censuses (Table 1, Figure 2). Although it declined at PR, its abundance was stable during the last half of the series (F tests, p < .001, p = .294). Two of the most abundant Eulaema on BCI (Ackerman, 1983) and PR also declined significantly at PR ( Figure S2) and further account for the community trend. Other, larger Eulaema increased. Rank clock plots suggested that in the long-term (PR) Eg. imperialis was declining but in recent years populations experienced a slight increase ( Figure S4). Regression analysis also showed that the Imperial bee either declined or increased, comparing two halves of the study at PR, while seven other Euglossa steadily increased and five declined (Figures 2, S1, and S3). Euglossa dissimula appeared to be stable despite much variation, whereas Eg. despecta increased slightly after 1979 ( Figure S1).
Community dynamics during the study period ( Figure S5) indicate there were small differences in species richness over time, slightly decreasing, although 2018 richness was higher than in 1979. There was a large peak of species turnover in 1998. Exclusion of Imperial bee data in our analyses, or dividing the almost 40-year duration of study into two parts, changed the long-term dynamics of all Euglossa from stable to increasing at PR, but had no effect on stability of the total euglossine community there; neither were affected on BCI (Figure 2). These parallel study periods at our two sites were closely similar for Euglossa spp. Remarkably, in total over 98% of individual bees consisted of only half the total species monitored in either study, one with all encountered orchid bees, and one counting only one of four genera. Many species were always in low abundance. Further, the more common species tended to increase while less common species decreased (Fisher test, p = .005, Table 2), with notable exceptions.
The larger bees, Eulaema and Eufriesea, and Exaerete, their parasite, had sufficient live sighting data for an analysis for PR. Of 11 species, 3 increased, 2 declined, and 6 were stable. We found a significant positive relationship between Exaerete and certain Eulaema and Eufriesea, the potential hosts of similar size (r = .143-.266, p = .0-.05; Table S3). The larger Eulaema (bombiformis, meriana) increased, along with their probable brood parasite Exaerete frontalis, as the smaller El. nigrita and El. cingulata declined. At PR, there was a considerably greater proportion of El. bombiformis, the largest orchid bee, compared to the original live sighting at BCI (Table 1). One relatively abundant Eufriesea, Ef. corusca, did not appear in either of our studies because it comes to vanillin attractant, which we did not employ.
The linear regression results for Euglossa, with 22 species monitored on PR and 19 at BCI, showed 7 increasing, 6 declining and 9 stable in the former, and 2 declining and 17 stable (2 increased slightly) in the latter (F tests, p > .05; Figures S1 and S3). One, Eg. dressleri, declined in both studies, while Eg. hansoni (p = .06) and Eg. variabilis (p = .09) increased on BCI. The BCI data (2009-2019) lacked Eg. crassipunctata, a common species at PR and on BCI (Table 1; Ackerman, 1983;Roubik & Ackerman, 1987). Overall, 26 and 11% of Euglossa species declined at PR or BCI, respectively. Eight rare species, including Euglossa, among the 42 species followed on PR, did not appear in the latter half of the 40-year PR study (also see Roubik, 2001).
Within 145 ENSO "climatic event" months (a strong or very strong La Niña or El Niño, 9 total events), our 58 census results were no different from those of 130 neutral climate periods at PR (p = .57, Table 2). None of the 33 "core" species at PR displayed a correlation between abundance and climate anomalies (Table S1). At both study sites, years categorized according to ONI events had no bearing on faunal composition of bee communities (adonis tests, PR: F 2,35 = 0.023, p = .967; BCI: F 2,7 = 0.407, p = .812; Figure S6).

| DISCUSSION
As found in other studies using permanent plots for a similar length of time (trees on BCI, Condit et al., 2017; nesting bees and wildflowers in the north temperate zone, Cane, 2021) nature is often not in equilibrium. Our findings contradict those of Faleiro et al. (2018). Their data were based on sampling in many sites for abundance differences among 37 orchid bee species in Brazilian Amazonia, some associated with disturbed habitats (see also below). Predictions were presented from a climate modeling approach that uses contemporary species presence/absence data. General decline, with some T A B L E 2 ENSO impact (3-month running averages given by the ONI), bee abundance, and trend tests for orchid bees (Euglossa, Eulaema, Exaerete, Eufriesea) monitored on BCI (Nature Monument) and Pipeline Road (Parque Nacional Soberanía) Panama, during 1979-2019 exceptions, was envisioned, accompanying climate change. We assume that the climate changes during our Panama studies encompass most of those seen over the time of orchid bee residence-including the largest changes during the Neogene-but we may be mistaken. Did our census data contradict presence/absence predictions for particular reasons? We remain largely ignorant of the evolutionary history of the bees we studied and the communities in which they perform. Cases of fossil studies cover far greater timescales (Botta, Dahl-Jensen, Rahbek, Svensson, & Nogués-Bravo, 2019;Clark & McLachlan, 2003), but our study systems do not allow such an extensive view. Although a strong El Niño dry season was positively correlated with orchid bee abundance at PR during the first 240 months of our study (Roubik, 2001), we found no apparent longer-term responses in the bee community or among individual species to the ONI events as a whole. Our census data included almost half the months in a total of 50 years of study. The ENSO events arose during 145 of those 596 months at the PR site, and the corresponding 58 monthly bee census results did not predictably vary to any extent that we could detect. In addition, the observation that seeds are more abundant on BCI in El Niño years (Wright & Calderon, 2006), thus more abundant flowers may then influence bee abundance, lacks confirmation here. We are only beginning to discover which resources are of major importance to most bees. We suggest that neither presence/absence data nor species dynamics in long-term studies are conclusive for judging species or community resilience in natural habitats and native communities. Both evolutionary and adaptive potentials remain to be deciphered. Nonetheless, a significant advantage for core species within protected reserves, on a scale comparable to foraging ranges of the species in question, further examined below, seems to be the inescapable conclusion of the present work.
Bee populations were mostly on the rise or stable, 74% at PR (39 years), 79% on BCI (11 years). A plausible explanation is the preservation and increase of major orchid bee resources. The increasing abundance of vine and liana flowers in the Neotropics (Schnitzer, Bongers, Burnham, & Putz, 2014) a prominent orchid bee resource (Dressler, 1982;Ramirez, Dressler, & Ospina, 2002;Roubik & Moreno, 2021), should support orchid bee population increases and stability. The primary floral resources required by orchid bees-sheltered nectar in long, tubular flowers, or sheltered pollen in poricidal anthers-and floral resin, may well be relatively less affected by rain or drought. Their phenology, however, is not a variable that we can analyze here. Euglossines differ from other common tropical bees, which hoard food in their nests and often use easily foraged, open flowers that have maximum abundance during dry season and early wet season (Roubik, 1989).
Remarkably, because they constituted less than 2% of the total individuals censused, 50% of orchid bee species were scarcely detected annually at our sites. We suggest that many studies that fail to encounter certain species have a significant challenge in verifying that their sampling is adequate or sufficient. We were nonetheless able to accumulate data to make inferences about the population dynamics of scarce bees ( Table 2) and found that less abundant bees were more likely to decline. Yet some of the most abundant species also decreased over time. Our two study sites were dissimilar in having a proportionately greater old forest area in BCI compared to the PR area, which, considering the several to many kilometer flight capability of orchid bees (Pokorny, Loose, Dyker, Quezada-Eu an, & Eltz, 2015;Wikelski et al., 2010), may be an underlying cause of differences, although bees on the island certainly access the mainland. A migratory lifestyle for euglossine males, but not their females-whose sheltered nests in the wild are difficult to encounter and very little studied-is conceivable. Male dispersal covers 50 km for a medium-sized Euglossa in the Yucatan, Mexico lowland (Pokorny et al., 2015). Male euglossines, which have no nest, can potentially find sheltered or favorable habitats within Central Panama over large distances. There are varied elevations and habitats, but we have no information on whether orchid bees use multiple nests or habitually migrate. The general result for mean ranks analysis on PR indicated some instability at the community level, as certain abundant bees became relatively less common. In contrast, the subset of species studied on BCI, the Euglossa, had rank abundances no different after 40 years, on the forest island.
That a dozen species were dominant, while 42% were uncommon among bees censused at both sites, suggests competing hypotheses. Extreme abundance of Imperial bees does not fit the expected lognormal distribution of species abundance in such diverse natural communities (Chisholm & Pacala, 2010). If bee behavior at the baits influenced the skewed numbers among species, then the time spent collecting a fragrance bait, or even its detection, may explain differences among species counts. In addition, a given bee species came to either 1, 2, or all 3 baits at PR, so species were sampled unevenly. However, with a fixed protocol, such variables are controlled by sampling rigor. Removal of individuals at baits during intervals of very high abundance does not diminish species counts, and a high (negative) correlation exists between counts and the time of first arrival (Roubik, 2001). Orchid bee abundance in the forest, as suggested by other studies, is indicated by the males at baits (Ackerman, 1983;Janzen, DeVries, Higgins, & Kimsey, 1982), and unless the sex ratio changes, their numbers should be reliable monitoring indicators. More importantly, because Ackerman's (1983) year-long sampling with 16 chemical attractants was very similar to our BCI data, the results with a single attractant and considering the same Euglossa species suggest live sightings agree with trapping data, and suggest the forest community was quite stable. In contrast, the less skewed dominance of the Imperial bee during BCI's initial survey (Table 1) and the design of traps so that not all large bees were collected, but smaller bees enter, demands further qualification and limits generalization based on the genus Euglossa alone. We purposely limited trapping of the three abundant largebodied genera, which live long and pollinate many flowers (Dressler, 1982). The Euglossa as a whole were increasing moderately after data of the single dominant species, Imperial bees, were excluded for PR. And on the contrary, on BCI no trend was found in either case (Figure 2). Absence of a trend does not signify stability in itself, but the absence of a trend on BCI compared to PR in the abundance of the same Euglossa species, with similar annual variation, strongly suggests the island experienced less change than the mainland.
Conservation efforts seek an expected result, but in studies such as ours, questions still remain about the status of succession in the habitat (Clark & McLachlan, 2003;Clements, 1936), the methods used to follow its residents (Bonebrake et al., 2010;White, 2019), and here, the multifaceted natural histories of bees and plants (Dressler, 1982;Ramirez et al., 2002). Tree mortality and growth in the 50-ha BCI plot has stabilized after drought in the 1980s, and the forest is changing (Condit et al., 2017). We also see orchid bees in a particular biological paradigm and time frame. The conservation of orchid bees, because they are not managed nor producers of commercial products, yet obviously visit and interact with many flowering plants, requires reasonably intact, mature natural habitat. In this status, during a 40-year interval, we found no decline or susceptibility to cyclic climatic changes, although suggested in much shorter studies (Ramirez, Hern andez, Link, & L opez-Uribe, 2015;Vega-Hidalgo et al., 2020). The magnitude of seasonal fluctuations in abundance was large, but we observed little directional change over time.

| CONCLUSION
Eight rare species in the first 20 years of study on PR did not appear later, and Eg. dressleri was declining at both our study areas, which might mean they ordinarily centered in higher elevations or different rainfall regimes (Roubik & Ackerman, 1987). Nevertheless, none of the core species in our study declined or increased during periods of a given ENSO climate event, usually a year or longer. Those drought or high rainfall periods (Table S1) greatly exceeded an average orchid bee generation time of 2 months (Dressler, 1982) and allow our study to assess changes in survival or reproduction during over 70 bee generations that were monitored during such conditions. As to whether "wild" bees are increasing, declining or maintaining their populations, and whether "climate change" is driving those populations in any particular direction seems now, on the available evidence, largely unknown for non-Anthropocene habitat. Our study suggests great overall stability during decades. Some tropical forests still maintain their organismal inventories, and may be more amenable to partial restoration or rescues than commonly believed. Modeling research efforts, although they may show what certain conditions portend, are not based on the needed deep history, but rather upon what can be garnered from contemporary environments and very few detailed studies. Studies like ours, which include highly adapted species assemblages, may still be improved upon. We need innovative ways to compare stable core species assemblages, their remnants, or ancient assemblages, with those that we see today.
We also hypothesize that community succession is occurring, whereby some abundant species decline and others take their places, exemplified with Euglossa and particularly the Imperial bee, Eg. imperialis. Some Eulaema and Imperial bees depend on abandoned nest cavities (e.g., from Atta ants) that are more common in disturbed sites, and those species range from Brazil to Mexico (Faleiro et al., 2018;Roubik & Hanson, 2004;Silva, Macêdo, Ascher, & DeMarco Jr., 2015). They remove competitive pressure as they decline in older forest and provide opportunities for other species. Although orchid bee communities in the large forest we studied were not declining, we still lack basic information to propose specific reasons for their dynamics. and José Alejandro Ramírez Silva made field collections and laboratory identifications.

ETHICS STATEMENT
Recognition of live orchid bee males was done in the field by DWR and sampling was carried out using traps on BCI by coauthors.