Using Ecological Niche Models and Niche Analyses to Understand Speciation Patterns: The Case of Sister Neotropical Orchid Bees

The role of past connections between the two major South American forested biomes on current species distribution has been recognized a long time ago. Climatic oscillations that further separated these biomes have promoted parapatric speciation, in which many species had their continuous distribution split, giving rise to different but related species (i.e., different potential distributions and realized niche features). The distribution of many sister species of orchid bees follow this pattern. Here, using ecological niche models and niche analyses, we (1) tested the role of ecological niche differentiation on the divergence between sister orchid-bees (genera Eulaema and Eufriesea) from the Amazon and Atlantic forests, and (2) highlighted interesting areas for new surveys. Amazonian species occupied different realized niches than their Atlantic sister species. Conversely, species of sympatric but distantly related Eulaema bees occupied similar realized niches. Amazonian species had a wide potential distribution in South America, whereas Atlantic Forest species were more limited to the eastern coast of the continent. Additionally, we identified several areas in need of future surveys. Our results show that the realized niche of Atlantic-Amazonian sister species of orchid bees, which have been previously treated as allopatric populations of three species, had limited niche overlap and similarity. These findings agree with their current taxonomy, which treats each of those populations as distinct valid species.


Introduction
ENMs). By using these methods, it is possible not only to evaluate niche dynamics of a species, but also to determine the most important environmental predictors affecting their potential distribution [24,43,44]. By assuming niche conservatism between sister species and using ENMs, we expect that the potential distribution of one species would reciprocally predict climatically-suitable areas for the other. If the reverse is true, then niches of both species have already evolved into different environmental requirements after allopatric speciation. Another framework to assess the realized niches of different species, while accounting for inherent biases in occurrence data was proposed by Warren et al. [45] and Broennimann et al. [46]. These methods measure niche overlap between each sister species, comparing it with a randomized distribution. Overlaps larger than expected at random indicate that the sister species occupy similar environmental spaces. Although such methods still consider only a portion of the species' realized niches, they control for unbalanced sampling efforts [46]. Along with environmental input data, these methods also require species occurrences. However, the distributions of insect species is still so fragmentary (the so called Wallacean shortfall; [47,48]) that the assessment of insect biodiversity and the effectiveness conservation actions involving these organisms is deeply affected, especially at broad spatial scales. Nonetheless, ENMs are important methods to address these issues [47,48], as well as to inform suitable areas for future field surveys [49].
Given the paleofloristic link between Atlantic and Amazonian rainforests and the controversy related to the taxonomy of Atlantic-Amazonian orchid bees, the main question to be answered in this study regards how much have changes in environmental conditions driven divergence between these lineages. Therefore, we used niche analyses and ENMs to quantify and define the niche features of sister orchid-bee species with allopatric distributions, measuring their niche overlap in an attempt to test whether or not they suffered niche evolution. Additionally, given the overall shortness of knowledge on insect species distribution on the Neotropical region, we used ENMs to predict their potential distributions in both biomes, highlighting relevant areas to future surveys. We expect that realized environmental niches of allopatric and related pairs of Amazonian and Atlantic species from the genera Eulaema and Eufriesea would show significant divergence with little overlap. Additionally, we also expect that the pairs of allopatric but unrelated Eulaema species would have low niche overlaps. Conversely, we expect that the pairs of sympatric and unrelated Eulaema species, which remained under the same environmental conditions during their whole evolutionary history, would show a larger environmental niche overlap when compared to pairs of allopatric and related Atlantic-Amazonian sister species. Finally, if the niches of sister species are conserved, they are expected to predict the potential distributions of one another.

Distributional data
We gathered distributional data on orchid-bee species from i) literature records (see Supplementary Material for complete list of published papers); ii) online databases [Species Link (http://splink.cria.org.br), Global Biodiversity Information Facility (http://www.gbif.org)]; and iii) personal collection data from one of us (AN). The species studied had at least 20 unique occurrence records (Table S1). We checked all occurrences regarding their taxonomic reliability, based on the known distribution of the species, and excluded dubious or unreliable records. We used Google Earth [50] to acquire proxy geographical information coordinates for records lacking exact geographic coordinates.

Environmental layers, modeling procedures, and evaluation
Species occurrences were overlaid in a grid of cells of size of 5 arc-min ranging from south USA to south South America. Using this same grid and considering all 19 bioclimatic variables from WorldClim (http://www.worldclim.org/), we derived principal components (hereafter PCs) to be used as new environmental layers in our distribution models. We selected seven PCs (<97% of the original climatic variation; Table S2) as our environmental variables. Given the overall biases and uncertain nature of ENMs, different algorithms may produce divergent species potential distributions [51]. Therefore, we generated the species potential distributions using four different algorithms: i) Envelope Score (ENVSC), a quantitative version of BIOCLIM [52,53]; ii) Mahalanobis Distance (MAHAL) [54], iii) GARP with best subsets [55], and iv) Maximum Entropy (MAX) [56,57]. While ENVSC and MAHAL are simpler algorithms that usually need presence data only to produce the species' potential distributions, MAX and GARP are artificial intelligence methods that are generally more complex, and correctly predict the species known occurrences more often [58]. We used the software MaxEnt to run MAX [57], and openModeller Desktop [59] for all the other three modeling algorithms.
We randomly divided all occurrences into ten 70%-30% training-testing subsets. We used the ROC threshold to balance both omission and commission errors while modeling the species distributions [60], to cut the suitability matrices of the modeled species in all modeling algorithms into presence/absence maps. True Skilled Statistics (TSS hereon; [61]), a threshold-dependent statistics, was used to assess model performance. TSS vary from 21 to +1, where negative and around zero values indicate that distributions are no better than random, while values near +1 represent perfect agreement between the observed and the modeled distributions. Acceptable models are those with at least 0.5, and excellent those with TSS around 0.7.

Quantifying the overlap between pairs of orchid-bee species
We used the same 19 bioclimatic variables from WorldClim in Broennimann et al.'s [46] framework, implemented in R [62], to compare niche features of all six species considered. We used a Principal Components Analysis (PCA) calibrated on the entire environmental space of the study background (Broennimann et al.'s [46] PCA-env) to measure the niche overlap of all species. We compared Atlantic-Amazonian species in pairs (El. atleticana versus El. meriana; El. niveofasciata versus El. bombiformis; and Ef. atlantica versus Ef. ornata). This routine compares the environmental conditions available for species within the full study background with its observed occurrences.
At first, this method encompasses both the niche and the density of occurrences of the species. It calculates the available environmental space, defined by the first two axes from the PCA-env, for each study area. This method corrects for potential sampling bias on occurrence records using a smooth kernel density function [46]. Later, the method obtains an observed niche overlap score for each species pair using Schoener's D metric [63], which varies from 0 to 1, representing a gradient of complete dissimilarity to fully overlapping niches. Then, two different randomization routines are used to test the hypotheses of niche evolution. First, the niche equivalency test compares whether the environmental niche overlap is different from random, while maintaining the original sample sizes. Secondly, the niche similarity test compares the niche overlap of one range randomly distributed over its background keeping the other unchanged (1R2), and then does the reciprocal comparison (1r2). We repeated each randomization process 100 times, producing a null distribution of overlap values to which the observed score was compared. If the observed overlap was significantly smaller, both occurrence sets were different, which would imply that the species were occupying distinct segments of environmental space [46]. Once we had only one pair of sister Eufriesea species with a minimum amount of unique occurrences, we did not compare them with other allopatric and unrelated Eufriesea sister species.
Determining priority areas for future surveys of orchid-bee species After we produced 40 presence/absence maps for each species (10 with each algorithm), we produced a mean consensual distribution map for each species with those which achieved TSS values .0.4. Many orchid bees depend on densely forested areas to occur, since those areas are expected to offer them suitable environmental conditions. Therefore, we used the information on Atlantic Forest remnants (http://www.sosma.org.br/) and the worldwide tree cover data [64] to respectively detect areas with at least 16 km 2 for the Atlantic (El. atleticana, El. niveofasciata, and Ef. atlantica) and the Amazonian species (El. bombiformis, El. meriana, and Ef. ornata) that are suitable areas for future surveys. The data on forest remnants in the Amazon and the Atlantic forests were converted to raster files with grid cells with 0.041˚(<4 km), allowing us to detect rainforest remnants with at least 16 km 2 , the same resolution to which the final species consensual maps were downscaled. Later, we identified grid cells where each orchid bee was predicted to occur in rainforest remnants with, at least, this minimum area.

Orchid bees' realized niche comparisons
As predicted, both allopatric and unrelated Eulaema species and all allopatric sister species from Atlantic-Amazonian forests (Eufriesea species included) had small niche overlap. Conversely, sympatric and unrelated Eulaema species had a large niche overlap. Generally, Amazonian orchid bees were associated with higher temperatures and less seasonal rainfall than Atlantic species (Fig. 1A-G; Fig. S1). Nevertheless, Amazonian species had a wider climatic range than Atlantic ones ( Fig. 1A-C, Fig. S1).
Pairs of sympatric ( Fig. 1D-E) and allopatric ( Fig. 1F-G) Eulaema species had a lower niche overlap and a lower niche equivalency. Niche similarity of El. atleticana with its Amazonian sister species, El. meriana, was smaller than expected by chance, but the reverse was not true (Table 2; Fig. 1A). Similarly, pairs of allopatric and unrelated Eulaema species showed smaller overlaps, and the similarity between Atlantic species and their related allopatric counterparts was also small ( Table 2; Fig. 1F-G). Conversely, sympatric but unrelated Eulaema species had high overlap of equivalent niches, but their realized niches were less similar ( Table 2; Fig. 1D-E), even for sympatric Eulaema species.

Potential distributions of orchid bees
Regardless of the pair of allopatric species considered, those inhabiting the Atlantic coast always had higher TSS values than Amazonian ones (Table 3). However, TSS bombiformis (red) and the Atlantic El. niveofasciata (blue). C) Niche overlap between the Amazonian Ef. ornata (red) and the Atlantic Ef. atlantica (blue). D) Niche overlap between the sympatric but unrelated Amazonian Eulamea species, El. bombiformis (red) and El. meriana (blue). E) Niche overlap between the sympatric but Atlantic Eulaema species, El. atleticana (red) and El. niveofasciata (blue). F) Niche overlap between the allopatric and unrelated Amazonian El. meriana (red) and the Atlantic El. niveofascita (blue). G) Niche overlap between the allopatric and unrelated Amazonian El. bombiformis (red) and the Atlantic El. atleticana (blue). The solid red and blue thin lines correspond to 100% of the available (background) environment for each species considered in the analyses. Red and blue shadings surrounded by thick lines correspond to the density of occurrences of each species per grid cell. values for the Eufriesea species were similar, independently of the modeling algorithm considered. In general, TSS values were either acceptable (higher than 0.5) or excellent (higher than 0.7), but the TSS values for El. meriana, Ef. atlantica, and Ef. ornata showed poor prediction rates for ENVSC and MAHAL (Table 3).
While both Amazonian Eulaema species showed wider potential distributions, usually ranging from Central America to central South America, the two Atlantic species showed narrower potential distributions, especially along the Atlantic coast in South America (Fig. 2), similarly to Eufriesea species. Nonetheless, the distribution of the Amazonian Ef. ornata was more constrained than that of the Amazonian species of Eulaema (Fig. 2). In general, Atlantic species did not project suitable areas for Amazonian species and the opposite was also true (Fig. 3). Finally, while allopatric species of Eulaema did not predict one another potential distributions, the sympatric species did.
Interior areas in the state of Bahia and eastern portions of Sergipe and Alagoas (Fig. 3) were priority areas for future surveys of Atlantic species, since they have large remnants of Atlantic Forest. Since the distribution of Ef. atlantica was more constrained than their Atlantic forest relatives, only interior areas in Bahia were indicated as potential for future surveys. Northern Brazil and northern South America (Fig. 3) are priorities areas for future surveys of Amazonian species. Similarly to Ef. atlantica in the Atlantic Forest, the constrained potential distribution of Ef. ornata in the Amazon revealed areas for future surveys only in northern Brazil, Guianas, and Suriname.

Discussion
Here, we showed that while allopatric orchid bees had a small niche overlap, an indicative that their niches adapted to different environmental conditions, sympatric species had very similar environmental niches. On the same direction, our ENMs showed that while the pairs of sympatric species predicted suitable climatic conditions for one another, the same did not happen with the pairs of allopatric species: ENMs for Amazonian species did not predict suitable areas for their counterparts in the Atlantic Forest and vice versa. Our results agree with previous studies on niche conservatism of sister species (e.g., [24,65,66]). Specifically, we found that sympatric but unrelated Eulaema species correctly predicted some portion of the environmental space of the other parapatric congeneric species, which suggests these species likely experienced similar environmental pressures throughout their evolutionary history. Conversely, we found that allopatric but closely-related orchid bees in the Amazon and Atlantic Forest did not predict suitable areas for one another, what challenges both theory and previous studies on niche conservatism [24,65,66]. Such patterns also hold if we consider the cross-comparisons of the realized niches of allopatric, sympatric, and sister allopatric orchid-bee species. Realized Niches among Orchid Bees The contrasting environmental niches found for the orchid bees from both biomes could be due to varying competition and parasitism rates by other orchid bees [67] or differences in micro-habitat preference in those biomes. Nonetheless, given the spatial scale of our study three other factors seem to be more likely. First, environmental conditions may have been more constant for Amazonian species than in Atlantic Forest. Consequently, Amazonian lineages would not need to adapt to new climatic conditions. As a result, species would have to develop physiological adaptations to cope with the expanding/contracting environmental gradient on the Atlantic coast. Second, natural characteristics of the two biomes, such as a larger area and narrower climate variation in the Amazon could facilitate the spread of Amazonian species. Also, differences in latitudinal and altitudinal variation between the Amazon and the Atlantic forest may to be key elements to explain the range and niche differences. The Atlantic Forest spreads over 27 degrees in latitude [68] and ranges from the sea level up to 2,700 m, whereas the Amazon Basin rarely reaches beyond 1,000 m [69]. Lastly, longitudinal variation is also important to explain such differences, since the dry forests of the inner parts of the Atlantic forest are climatically different from those closer to the coast [70]. Therefore, past cycles of forest expansion and contraction coupled with contrasting environmental conditions in both Amazonian and Atlantic biomes may have forced each subpopulation inhabiting these areas to adapt to different environmental conditions, which later shaped their niches and distributions, promoting their subsequent speciation.
The importance of changing environmental conditions for the divergent speciation of the Atlantic-Amazonian species following forest expansions and contractions were also found by Ló pez-Uribe et al. [71], at least for two of the species analyzed (El. meriana and El. bombiformis). However, it is noteworthy that those authors, without further explanation, considered both Amazonian and Atlantic populations of El. meriana and El. bombiformis as belonging to the same taxonomic entity, what led them to group all occurrences to generate their ecological niche model, contrarily to the established taxonomy classification that considers them as four different taxa [29,30,39,72,73], as also supported by our data. Moreover, their own molecular data suggest that Atlantic populations of both Eulaema species seem to be monophyletic. Even though these bees have high dispersal rates [74][75][76], Ló pez-Uribe et al. [71] showed that their distribution may have been severely affected by South American paleoclimatic instability, which differed between Amazonian and Atlantic forests and supports the niche differences we found for the orchid bees considered in this study. Future studies such as that could be useful to assess how orchid bees responded to past climate changes in South America.
One criticism about niche shifts between related taxa is the use of biased species occurrences data to study realized niche and infer processes driving their fundamental niche. Since species distribution is the intersection between biotic, abiotic, and historically reached habitats, much of the fundamental niche may still not be accessed [42,77,78]. Additionally, sampling bias may also prevent the assessment of all reachable areas [42]. Finally, the assumptions of ENMs techniques are subjected to criticisms, since key biological processes (e.g., migration, demographic features, species competition, and predation) are implicit, while a more explicit treatment would be more efficient to assess the species' fundamental niche [79][80][81].
Although our analyses were based on the species realized niches, we believe that our results may help understand speciation patterns of these species. Our results agree with previous studies using the same methods to discuss niche conservatism of exotic species [82,83] and the division of a single lineage into several ones [84][85][86][87]. Thus, given the idiosyncratic environmental niches, our results support the splitting of Amazon and Atlantic populations into different taxonomic entities [28][29][30][31][32][33][34][35][36][37][38]88].
The distributions of insect species in the tropics suffers from knowledge gaps [47,48], specially taxonomic identity (Linnean shortfall) and spatial distribution (Wallacean shortfall). Orchid bees are no exception [89], especially regarding taxonomic issues [90]. Nonetheless, ENMs are important tools to assess the distribution of ''data deficient'' species [48], allowing them to be included in broad-scale conservation studies. Because orchid bees depend on humid forested areas [26,27], some have been reported as threatened [91] or extinct in fragmented areas in the Atlantic Forest [28,33]. Therefore, integrating ENMs and information on original vegetation remnants are urging. Given this scenario, future surveys in the Amazonian and Atlantic Forest fragments highlighted by our analyses may add important information on species distribution and their natural history.

Conclusions
In this study, based in SDMs and analyses of the species' realized niches in South America, we have shown that sister Atlantic vs. Amazonian orchid bees, namely El. atleticana vs. El. meriana, El. bombifromis vs. El. niveofasciata, and Ef. atlantica vs. Ef. ornata, had different potential distributions and realized niche features. Nonetheless, we also showed that sympatric species have more similar environmental niches than that observed for the pairs of sister Atlantic-Amazonian orchid bees. These results support the view that the orchid bee sister lineages from Amazonian and Atlantic forests are currently under divergent selection, and consequently should be treated as different species. Additionally, we also highlighted important areas to be considered in future surveys in both biomes. We believe that similar approaches to those we used here are important to unveil both the biogeographic history of the South American fauna, as well as to optimize and guide future biological surveys of these species in the continent. Finally, as already developed for El. meriana and El. bombiformis, we recommend that future studies should also employ phylogeographic analyses of these (and other) orchid-bee species so we could better understand how past range expansion and contraction events in the continent may have originated the faunal and floral similarities currently found in those biomes. Figure S1. Results obtained from the PCA-env approach developed by Broennimann et al. (2012). A) PCA-env obtained with WorldClim's 19 raw environmental variables, considering the whole extent of South and Central Americas, used as the background in the analysis. PC percentages refer to the amount of variation explained by each PCA axis. B) and C) Contributions of each environmental variable to the first and second PCA axes, respectively. doi:10.1371/journal.pone.0113246.s001 (TIF)