Frontiers of Biogeography Non-overlapping climatic niches and biogeographic barriers explain disjunct distributions of continental Urania moths

Larvae of Urania moths feed exclusively on Omphalea plants, which are widely distributed in the Neotropics. However, the distributions of the two Urania species in this region are disjunct. This distributional pattern could derive from the presence of the Andes, but it could also be related to differences in ecological niches, the presence of negative interactions, or the absence of conditions that can only be observed at a habitat level. We tested whether differences in the ecological niches of continental Urania moths play a role in their disjunct distribution. Using species records and climatic variables, we characterized the ecological niches of Urania moths and their host plants and analyzed the overlap of the moths’ niches. Using ecoregions as a proxy of habitat-level environmental conditions, we explored the role of host plant availability on the moth distributions. Suitable conditions for the species were widespread, with a lack of suitability mostly restricted to the Andean highlands. The two-moth distributions were closely related to that of their host plants. There was medium-high overlap of niche models when available conditions were considered; however, niche overlap was not found to be statistically significant. Our results corroborate the barrier effect of the Andes on the dispersal of these moths, but they also show that niche differences contribute to the disjunct distributions of U. fulgens and U. leilus . Furthermore, other non-climatic factors appear to play a crucial role in the disjunction of the species ranges in areas where overlapping suitable conditions are continuous. Our findings support speciation in Urania moths as allopatric and indicate that their disjunct distributions can be attributed to multiple factors. Other studies exploring the causes of similar distributional patterns should consider that a single factor may not be enough to explain such patterns.


Introduction
There are several combined factors that determine the regions of the Earth that a species occupies. For any given species, the literature mentions that these factors include suitable climate and habitat, availability of food, lack of natural enemies and competitors, and the capacity to reach favorable localities (Good 1931, Cain 1944, Darwin 1859, Udvardy 1969, Soberón and Peterson 2005. Darwin (1859) stressed that the environment alone cannot explain geographic distributions of species because favorable environmental conditions often occur outside of a species range. Another factor affecting the distribution of a species is habitat suitability, which is understood to be the availability of certain structural conditions in the landscape that are necessary for survival or reproduction. For instance, species distribution may be affected by the availability of nesting places or the presence of conditions for the establishment of propagules (Andrewartha and Birch 1954). Because habitat-level characteristics can only be appreciated on a finer scale than climatic conditions, and information is often limited (Lindenmayer et al. 1991), these factors are usually disregarded in macroecological analyses of species distributions (Soberón 2007).
Many examples of favorable climatic conditions that occur beyond a species distribution exist in nature (non-equilibrium with the environment; Peterson 2003, Svenning and Skov 2004, Seliger et al. 2020); however, quantitative analyses exploring this phenomenon are limited. To explore how distributions out of equilibrium with environmental conditions could derive from multiple factors, we present an example of two moth species with apparently similar niches but allopatric geographic distributions. Our analysis is quantitative, and the methods are novel in the sense that overlap is measured using a model of niche, in environmental space, and considering available environmental conditions for each species.
The neotropical genus Urania includes four species of diurnal moths (U. boisduvalii, U. leilus, U. fulgens, and U. sloanus; Smith 1991, Nazari et al. 2016); however, U. solanus has become extinct due to habitat degradation (Lees and Smith 1991). These moths feed during their larval stages exclusively on plants of the genus Omphalea (Smith 1983, Lees and Smith 1991, Smith 1991, Smith 1992, Nuñez-Penichet and Barro 2020 with no reported preference for a particular species of these plants. The distribution of Urania species has been described in several observational studies (Williams 1937, Williams 1958, Young 1970, Smith 1972, Smith 1983, Lees and Smith 1991, Smith 1991, Meerman and Boomsma 1997, Barro and Rodríguez 2005, Murillo-Hiller 2008. Recently, Nuñez-Penichet et al. (2019) used an ecological niche modeling approach to characterize the potential distribution and migratory routes of U. boisduvalii, which is endemic to Cuba, but no similar studies exist for the two continental species of this genus.
Given the close relationship between the moths and their host plants, Smith (1983) hypothesized that the distribution of the Uraniidae species must exactly coincide with the distribution of plants of the genus Omphalea. In the continental Neotropics, the occurrences of the host plants of Urania are widely distributed from Veracruz, Mexico to Bolivia (Lees and Smith 1991) in a broad variety of habitats like humid lowlands (e.g., swamp forest, along rivers, on the edges of sandy beaches), costal shrubs, as well as limestone hills and other karstic areas (Smith 1992). However, the distributions of the two Urania species in this area are disjunct. U. fulgens is distributed from Veracruz to the north of Colombia (with a subspecies, U. f. poeyi, endemic to Cuba) and U. leilus is distributed from the south of Colombia to Bolivia (Lees and Smith 1991). This allopatric distribution is intriguing because the availability of the host plants does not seem to be a limiting factor for the establishment of either of these species.
Despite the intriguing disjunction in the distribution of these species, no clear explanation has been proposed. Only Smith (1972) and Nazari et al. (2016) mentioned the Andes as a potential barrier for the dispersion of the two species of moths. Besides this topographic barrier, other types of barriers (at a macroecological scale) could be limiting the distribution of these two species. For instance, another potential explanation for the disjunct distribution of continental Urania moths could be a difference in their ecological niches. Low similarity in the preferences of environmental conditions of these two species could explain why they have non-overlapping distributions. One approach to address this question is by analyzing the niche overlap between these two moth species. This method has been used previously, among other applications (i) to explore the overlap between a species native and invasive range (Banerjee et al. 2019); (ii) to analyze niche overlap between co-occurring native and exotic species (Pascual-Rico et al. 2020);and (iii) to complement species taxonomic distinction by considering the niche overlap between sister taxa (Zhang et al. 2014).
In order to understand the causes for the disjunct distributions of the moths, in this study, we aimed to: (1) characterize the ecological niches of the continental Urania moths and its host plants; (2) evaluate whether the niches of these two moth species are different; and (3) explore landscape-level factors that may contribute to maintaining the disjunct distribution of these moths. We used ellipsoidal envelope models (Farber and Kadmon 2003) to characterize the ecological niches of Urania moths and Omphalea plants. We projected these models to current and past scenarios in environmental and geographic spaces, to evaluate how strong the barrier effect of Andean highlands is. To test whether differences in the ecological niches of the moths exist, we measured niche overlap using a novel approach based on ellipsoid envelopes (similis Qiao et al. 2016) that takes into account the very heterogeneous distribution of climates available for species within their ranges (a point often disregarded). We also searched for landscape-level barriers in the area where the moth distributions become disjunct by considering a layer of terrestrial ecoregions.

Data
The occurrence data of U. fulgens, U. leilus, and Omphalea species were downloaded from the Global Biodiversity Information Facility database (GBIF; https://www.gbif.org/). We cleaned the data by removing duplicates and records with uncertain or missing coordinates, as previously described by Cobos et al. (2018). We also removed records included at distances <5 km from another locality (spatial thinning; Anderson 2012) to avoid model overfitting derived from multiple records in the same pixel in the raster variables. After these steps, we obtained 180 records for all species of Omphalea genus (of the 268 original records), 108 for U. fulgens (of the 393 original records), and 34 (of the 72 original records) for U. leilus (Fig. 1). All species of Omphalea were joined to model a niche that characterized the environments in which all the neotropical host plant species were present (Nuñez-Penichet et al. 2016, Nuñez-Penichet et al. 2019. Host plant species from other tropical areas were not included because they were not relevant to this study. Data cleaning and spatial thinning were undertaken in R 3.6.2 (R Core Team 2019) using the package ellipsenm (Cobos et al. 2020; available at https://github.com/marlonecobos/ellipsenm).
We used bioclimatic variables at a 2.5' resolution (~4 km) for current and past scenarios (mid-Holocene and Last Glacial Maximum) from the WorldClim database version 1.4 (available at https://www.worldclim.org/). We excluded the variables mean temperature of wettest quarter, mean temperature of driest quarter, precipitation of warmest quarter, and precipitation of coldest quarter, because they may contain artifacts as a result of combining temperature and precipitation information (Escobar et al. 2014). The variables were masked to a hypothesis of accessible area for each species (Fig. 1). We tested the correlation of the variables and selected the ones with values of r ≤0.8 and that were more biologically relevant for the species of interest, based on our knowledge of their natural history (Table 1; Simões et al. 2020). Variables from past scenarios corresponded to one General Circulation Model (GCM), the Community Climate System Model (CCSM).
When characterizing ecological niches and estimating potential distributions, it is important to consider the area that has been accessible to the studied species (M; Soberón and Peterson 2005, Barve et al. 2011). These accessible regions restrict the set of conditions used for modeling, which is important especially for correlative approaches (Peterson 2011). Although our modeling approach (ellipsoid envelopes) does not need background information (see Ecological niche modeling below), accessible areas were necessary to perform a test of statistical significance when measuring niche overlap. We determined accessible areas for our three species based on known occurrences and dispersal abilities. The highest points of the Andes served as limits of accessible areas for the moths where their ranges adjoin, because no record of one or the other species exists on the side of the mountain range that is opposite to its known distributional area (Fig. 1).

Ecological niche modeling
We created ellipsoidal envelope ecological niche models (Farber and Kadmon 2003) for the three taxa of interest, the species of Omphalea (all records Mahalanobis distances (MD) were calculated from the centroid to each environmental combination in the area of interest. Values of suitability (S) were obtained by finding a multivariate normal distribution from MD; S ~ exp(-0.5 × MD) and assuming that the highest values of suitability were closer to the centroid (Maguire 1973). The limits of the ellipsoids were found using a chi-squared cumulative probability distribution on the MD (Etherington 2019) and by imposing a value to exclude extreme data (e.g., 0.95 excludes 5% of occurrences with the longest MDthe extremes). Finding the limits of the ellipsoid allows us to distinguish suitable from non-suitable areas, and excludes extreme environmental values that may have been derived from sampling biases or other artifacts during the georeferencing processes. We used ellipsoids as niche models because it has been hypothesized that they may yield a closer approximation to a fundamental niche estimate than other popular methods that fit arbitrary shapes to the observation points without much concern about its biological interpretation (Drake 2015, Jiménez et al. 2019. Furthermore, ellipsoids require considerably fewer assumptions and decisions regarding parameters compared with other popular methods, like Maxent (Merow et al. 2013). Finally, the volume of the fitted ellipsoids can be easily obtained, and it is a measure of niche breadth.
Our modeling approach consisted of using the clean-thinned occurrence data and the variables selected for each species to fit ellipsoids that exclude a percentage of marginal data (e.g., 5%, as described above). Because the data always contain biases and no perfect setting exists in modeling exercises, we aimed to explicitly include variability. To accomplish this, we produced 10 ellipsoid envelopes based on distinct subsamples of 75% of the occurrence data (replicates). For each species, our final models were obtained by averaging the centroids and the covariance matrices of the 10 replicates produced (we checked for the average of covariance matrices to be positive definite). For the sake of comparison, geographic projections of all models were undertaken based on the areas accessible to the host plants (Fig. 1). In addition, values of prevalence were measured in geography (proportion of total suitable area) and in environmental space (proportion of the multidimensional space identified to be suitable considering only distinct combinations of environmental conditions). Ecological niche modeling was performed using the ellipsenm R package.
We represented ecological niches and suitability levels in environmental and geographic space. Geographic projections were binarized (U. fulgens: 0 and 1; U. leilus: 0 and 2; where 1 and 2 represent suitable values) and then summed to represent how suitable areas for each species alone (values of 1 or 2) and coincident suitable areas (values of 3) were distributed. The threshold for suitability was obtained using the approach described above, aiming to exclude the 5% of the data with the most extreme values. Model projections to geography and other calculations described above were performed for current and past climatic scenarios.

Ecological niche overlap
To compare the niches of the two moths, we used a test based on ellipsoid envelopes (see above, Ecological niche modeling), in which we measured the amount of overlap between ellipsoids corresponding to the species. To measure overlap, points of an environmental background are checked to be inside or outside the ellipsoid of species A, B, or both. We considered two types of overlap for our analyses. The first type, full-overlap, measures overlap of the ellipsoids considering background data (1,000,000 points) that fills uniformly the multidimensional space containing the two ellipsoids (Fig. 2). Here, points are extracted randomly from such uniform distribution and checked to determine where they are located considering the two ellipsoidal niches. In the second type, backgroundoverlap, measurements were performed using only the environments that are accessible to the species (which may or may not fill all the volume of the two ellipsoids). In the later type, points are extracted from this non-uniform, existing background, and checked to be inside of one, the other, or both ellipsoids (Fig. 2). The value of overlap is then given by the proportion of total points in the intersection of two ellipsoids (A and B): overlap = J = |A ∩ B| / |A ⋃ B| (bars denote To test the statistical significance of the value found using the background-overlap approach, we measured the overlap between ellipsoid pairs that were created with randomly sampled points from the background of each species (the N for random samples is the number of occurrences per species). We repeated this process 1,000 times and compared these results with the observed value of overlap. In this context, the null hypothesis is that the two ellipsoids fitted to actual observations overlap at least as much as random-data ellipsoids. To reject this hypothesis, the observed value of overlap must be as extreme or more extreme than the lower confidence limit (5%) of the overlap values derived from comparing random-data ellipsoids derived from each species background. The significance test was not performed for the results of the full-overlap approach (because the uniform distribution of background points used in this approach does not consider limitations derived from accessible environmental conditions), but we illustrated it for comparison purposes. These analyses were implemented for the three scenarios considered in our analyses, in which the species niches were characterized using only current data, but available environmental conditions were modified by changing the scenario when testing overlap and statistical significance. We performed these analyses using the routines implemented in the ellipsenm R package. The code used to perform ecological niche modeling and overlap analyses is available at https://github.com/claununez/Continental_Urania.

Landscape-level barrier identification
In addition, we analyzed the ecoregions present in our study area, especially the ones in which records of the host plants exist and the ones that are close to where the moth ranges adjoin. We did this to detect potential landscape-level barriers that may be limiting the distribution of the moths. An ecoregion was considered as used if it contained one or more records of Omphalea species. The layer of terrestrial ecoregions used was downloaded from the site of the World Wide Fund for Nature (https://www. worldwildlife.org/publications/terrestrial-ecoregionsof-the-world; Olson et al. 2001).

Results
The geographic projections of ecological niches of all species, for current and past scenarios, were similar, showing widespread suitable areas across the Neotropics (Fig. 3, Figs. S1-S2). Unsuitable areas for all the species were found in the north and central parts of Mexico, the north and central areas of the Andean mountain range, and the southernmost parts of Brazil. The southernmost part of Mexico, a large portion of Guatemala, Panama, and Mato Grosso and Bahia (Brazil), were also identified to be unsuitable for U. leilus (Fig. 3, Figs. S1-S2). Geographic patterns of suitability changed in past scenarios, displaying a decrease in highly suitable areas, especially for U. leilus and Omphalea species. A reduction of suitable areas for U. fulgens was more pronounced in the mid-Holocene, whereas these areas were even smaller in the Last Glacial Maximum for U. leilus and Omphalea species (Table 2).  In regard to environmental space, the ellipsoidal niche of U. fulgens had a larger hyper-volume than the ellipsoid of U. leilus (13.09 × 10 10°C .mm 2 and 45.23 × 10 7°C .mm 2 , respectively). The niches of both moths were found to be associated with high temperatures and precipitations levels, although U. fulgens was found to be present in regions with higher precipitations, regarding the wettest months (Fig. 3, Figs. S1-S2). The niche for Omphalea spp. had a wide range in annual precipitation, although it was located in regions of low precipitation seasonality and high annual mean temperatures (Fig. 3, Figs. S1-S2).
The suitable areas for both species of moths that coincided in geographic space were also widespread; 95% of the suitable area of U. leilus overlapped with U. fulgens, and this value was higher in the mid-Holocene and Last Glacial Maximum scenarios (~99% for both). In contrast, the suitable areas for U. fulgens that coincided with those of its sister were less widespread during the past scenarios (70%, 65%, and 63% for current, mid-Holocene, and Last Glacial Maximum conditions, respectively; Fig. 4, Fig. 3S). In regard to the entire area of projection for our models (i.e., the accessible area of Omphalea spp.), 49%, 40%, and 41% of it was found to be suitable for both species of moths under current, mid-Holocene, and Last Glacial Maximum conditions, respectively. In environmental space, the three scenarios, conditions used by U. fulgens included most of the climatic characteristics of the areas in which U. leilus was present (Fig. 4,  Fig. 3S). Under current conditions, the value of full overlap between the ellipsoidal models of the two moths was low (J = 0.13); however, the overlap based on the available background for both moth species was considerably higher (J = 0.62). Nonetheless, based on the significance test, the null hypothesis-that the niches overlap-was rejected (p = 0.003; Fig. 4). Therefore, the niches were considered to be different. Similar results from background-overlap analyses were found for the mid-Holocene and Last Glacial Maximum scenarios (J = 0.67, p = 0.011; and J = 0.63, p = 0.016, respectively; Fig. 3S).
When we analyzed the overlap in the region where ranges of U. fulgens and U. leilus adjoin, we identified three possible connection points between the suitable areas for these moths (Caja Seca on Maracaibo lake, Venezuela; Cúcuta, Colombia; and close to El Placer in the National Park Cordillera Natural de los Picachos, Colombia; Fig. 5). The environmental conditions in those areas were found to be suitable only for U. fulgens in the current scenario, but some of them were found to be suitable for both moth species in past scenarios (Fig. 5).
The ecoregions that we found to be occupied by Omphalea spp. included mostly moist forests, with a few exceptions, such as pine-oak forests in Central America and a few other smaller ecoregions in eastern Brazil. In the region where the ranges of the moths adjoin, Omphalea species were only present in moist forests. A clear barrier of non-used ecoregions was observed across this region, including the areas where the three possible connections of suitable conditions exist (Fig. 5).

Discussion
We found that, although it was clear that U. leilus had a narrower niche than its sister species, the environments in current and pasts scenarios, used by U. fulgens and U. leilus, and consequently their niches, were relatively similar. However, although the overlap value found between the two niches appears large (J = 0.62, 0.67, and 0.63, when using available environments in current, mid-Holocene, and Last Glacial Maximum scenarios, respectively), it was significantly lower than random expectations (Fig. 4,  Fig. S3). Therefore, the niches of the two moths do not overlap. The apparently high overlap value found for the two niches concurs with the niche conservatism theory (Peterson et al. 1999, Soberón and Miller 2009, Peterson 2011, because the ancestral populations of these species apparently separated recently (≤2.7 Mya; Figure 5. Close up of the Andean region where distribution ranges of Urania fulgens and U. leilus adjoin. Representation of the coincident suitable areas between the two moths is presented, highlighting the possible causes for their disjunct distribution. Nazari et al. 2016). However, the absence of statistically significance for the overlap found may be due to a difference in the combination of environmental factors in their respective distribution ranges. Despite the medium-high similarities of conditions used by the two species, the disparity in the portion of environments used considering what is available, supports the hypothesis that their distributions are disjoint because of niche differences. However, the presence of either geographical or some other type of dispersal barriers has certainly contributed to maintain these moth ranges separated, especially in areas where connectivity of suitable conditions is present.
Observations of individuals of the Urania genus have been reported in the literature at elevations ranging from 0 -1,200 m ( Rodríguez 2005, Murillo-Hiller 2008), with the highest elevation recorded at 2,430, according to data from GBIF (Antioquia, Colombia). These observations suggest that even for these moths, which are strong fliers that can move long distances (Dudley and DeVries 1990, Smith 1992, Murillo-Hiller 2008, the Andes mountain range may represent a dispersal barrier. This also coincides with reports by Smith (1972) as well as Nazari et al. (2016). These last authors mention that a South American lineage of Urania may have dispersed to Central America after the formation of the Isthmus of Panama (~15-11 Mya) and that the uplift of the Andes ~5-2.7 Mya divided the ancestral distribution range of these moths. Similarly, our models support these former observations because the high Andean regions were found to be non-suitable for the host plants as well as the two species of Urania (Fig. 3). These high grounds were also identified as unsuitable for Omphalea species at a landscape-level when compared with the occupied ecoregions of these plants. This may explain the absence of Urania moths in the few connected areas that were detected between the niches of these two moth species in the region where their ranges adjoin (Fig. 5). Our results also showed that suitable areas for the host plants coincided with the most suitable areas for the moths (Fig. 3), which further agrees with previous suggestions that the presence of these plants determines the distribution of Urania (Smith 1983, Lees and Smith 1991, Smith 1992, Nuñez-Penichet et al. 2016). However, not all the areas with Omphalea plants were found to be occupied by Urania moths at any given time (Smith 1983) because, it is hypothesized, the Andean and the non-used ecoregion barriers prevent U. fulgens from exploring the accessible areas of U. leilus and vice versa.
The niche and geographic projection of suitable areas for U. fulgens were broader than those of its sister species in all considered scenarios (Fig. 3, Figs. S1-S2), which could be related to the use of more heterogeneous environments by U. fulgens. Because of this broader niche, a large portion of South America located beyond the Andes mountain range was found to be suitable for U. fulgens, even though this species has not been recorded in this area. In contrast, suitable areas for U. leilus were closer to its actual geographic distribution in South America, and included only a few areas in Central America (Fig. 3, Figs. S1-S2).
Suitable areas for U. leilus included mostly localities where rainforests were the dominant vegetation, thus the absence of suitability in Panama was unexpected (Fig. 3, Figs. S1-S2). The areas of Brazil that we found to be non-suitable for this species were a part of the Mato Grosso Seasonal Forest, the Caatinga, and the tropical savanna of the Cerrado, which are known to be arid with very marked seasonal rain (Klink and Machado 2005). The environmental characteristics in these areas were completely different from those present in the Amazonian region and in the known localities for this moth. The high climatic suitability that we found for U. leilus in a large part of the Cerrado could be associated with higher levels of precipitation compared with other parts of this ecosystem (Fig. 3, Figs. S1-S2; Soley-Guardia et al. 2014). Therefore, precipitation may play a crucial role in determining the distribution limits of U. leilus in its accessible areas. In contrast, projections of the niche of U. fulgens in South America indicated that this species may also be able to tolerate lower levels of precipitation than its sister species.
Our results from niche overlap analysis based on the species backgrounds and statistical significance tests assume that: (i) niches can be represented by ellipsoids, and (ii) measuring overlap requires consideration of the actual distribution of accessible environments. Although the last point has been almost universally ignored throughout the niche literature, except for a few reports (Hirzel et al. 2002, Broennimann et al. 2012, it is extremely relevant. As shown in our results, overlap values obtained with assumption of uniform distributions of combined environmental factors (full overlap) can be very different from those obtained based on the actual availability of environmental data. This point was first raised by Colwell and Futuyma (1971), but it has seldom been taken into account in the ecological niche modeling literature. In the terminology of Peterson et al. (2011), we are comparing existing niches, which means the actual existing environmental combinations inside the models of "fundamental niches" (the ellipsoids). Although using ellipsoids as a model for niche is not new (e.g., Osorio-Olvera et al., 2020), including information about the highly non-uniform environmental spaces at the time of testing for overlap is a rather recent practice (Qiao et al. 2016). In fact, the test of statistical significance of the observed overlap that is based on available environmental conditions is new, and only implemented by Cobos et al. (2020) in the R package ellipsenm. As our overlap analyses were performed in environmental space, we were able to test overlap using niches characterized with current data but exploring different scenarios of accessible environments (current and two past scenarios). This is a new practice that gave us the ability to test if distinct climatic configurations would affect the results of overlap analyses. We hope that with this example, we contribute to illustrate how informative it is to consider the shape of environmental space at the time of assessing ecological niches similarities.
In summary, the combination of the Andes acting as a geographic barrier, the absence of overlap between the moth niches, and the presence of a landscape-level barrier due to unsuitable ecoregions for the host plants (Fig. 5), may explain the disjunct distribution of U. fulgens and U. leilus. Other factors at a microhabitat level and biotic interactions should be explored in future studies because they may also influence the distribution patterns of these moths.

Supplementary Material
The following materials are available as part of the online article at https://escholarship.org/uc/fb Figure S1. Predictions of ecological niche models for Urania fulgens, U. leilus, and Omphalea species, under the mid-Holocene conditions. Figure S2. Predictions of ecological niche models for Urania fulgens, U. leilus, and Omphalea species, under the Last Glacial Maximum conditions. Figure S3. Coincident suitable areas of Urania fulgens and U.leilus in past scenarios.