Introduction
Among birds, Psittacidae is a family of major conservation concern. Nearly half of the 344 species in the world (Snyder et al. Reference Snyder, McGowan, Gilardi and Grajal2000) are considered threatened under the IUCN Red List criteria, including nine already ‘Extinct in the Wild’ due high levels of habitat loss and trapping for illegal trade (BirdLife International 2015).
Although research on the ecology and behaviour of psittacids has increased in the last decade, there is an urgent need to determine trends in distribution at larger scales (Martin et al. Reference Martin, Perrin, Boyes, Abebe, Annorbah, Asamoah, Bizimana, K, Bunbury, Brouwer, Diop, Ewnetu, Fotso, Garteh, Hall, Holbech, Madindou, Maisels, Mokoko, Mulwa, Reuleaux, Symes, Tamungang, Taylor, Valle, Waltert and Wondafrash2014). The few monitoring programmes for this family are usually focused on one, often globally threatened species, or are limited to narrow geographic areas (but see Hille Reference Hille, Wiedenfeld, Lezama-Lopez, Brightsmith and Patten2014). Status assessment of range and population trends for psittacids are usually based on static representation of species distribution (range maps; Snyder et al. Reference Snyder, McGowan, Gilardi and Grajal2000). However, range maps could provide a misleading interpretation of trends because they lack estimates of uncertainty in under- or over-prediction and assume homogeneous probability of occurrence across the range (Peterson et al. Reference Peterson, Soberón, Pearson, Anderson, Martínez-Meyer, Nakamura and Araújo2011). Data on presence and absence are required in order to track population changes, but most studies are limited in time and resources, and many apparent absences are in fact lack of detections (Kéry & Schmidt Reference Kéry and Schmidt2008). Standardised survey methods and appropriate analytical tools allow for the inclusion of heterogeneous detection probabilities, and the calculation of reliable estimates of presence or absence and their associated uncertainty (MacKenzie et al. Reference MacKenzie, Nichols, Lachman, Droege, Royle and Langtimm2002). Occupancy models have been widely used for monitoring avian populations (Baumgardt et al. Reference Baumgardt, Sauder and Nicholson2014). These models use repeated detection and non-detection data (detection histories) at each location to jointly estimate the probability of presence (Ψ), and the probability of detection (p; MacKenzie et al. Reference MacKenzie, Nichols, Lachman, Droege, Royle and Langtimm2002). Indeed, the most common application of detectability estimates is to determining whether a species is, in fact, present at a given site when not detected, but fewer case studies provides evaluations of the level of confidence that can be placed in a particular non-detection observation (Wintle et al. Reference Wintle, Walshe, Parris and McCarthy2011). The accuracy with which a non-detection could be interpreted as a true absence, may have direct implications in our ability to confidently interpret a current species distribution and hence, our capacity to monitor temporal and geographical changes (Garrard et al. Reference Garrard, Bekessy, McCarthy and Wintle2014).
In Venezuela, as in several other Neotropical countries, the high diversity of psittacids is combined with an increasing rate of land transformation (Rodríguez et al. Reference Rodríguez and Giraldo Herández2010), illegal wildlife trade (Sánchez-Mercado et al. Reference Sánchez-Mercado, Asmüssen, Rodríguez, Moran, Cardozo-Urdaneta and Morales2020), and limited resources for monitoring and conservation effort (Rodríguez Reference Rodríguez2014). The IUCN reports declining regional trends for 34 of the 50 Psittacidae species occurring in the country (Table 1), and six species are already under some threat category in the Venezuelan Red Data Book (Rodríguez et al. Reference Rodríguez, Rojas-Suárez, García-Rawlins and Rojas-Suárez2015). However, a more detailed national assessment requires the evaluation of current status and trends for the whole family. As a first step, a national bird monitoring programme was carried out in 2010 as a component of the Neotropical Biodiversity Mapping Initiative (NeoMaps; Rodríguez et al. Reference Rodríguez, Rodríguez, Ferrer-Paris and Sánchez-Mercado2012, Ferrer-Paris et al. Reference Ferrer-Paris, Rodríguez, Good, Sánchez-Mercado, Rodríguez-Clark, Rodríguez and Solís2013). This programme provides an important source of detection and non-detection records, ideal for fitting occupancy models allowing the establishment of a baseline to evaluate temporal and spatial changes in distributions (Ferrer-Paris et al. Reference Ferrer-Paris, Sánchez-Mercado, Rodríguez-Clark, Rodríguez and Rodríguez2014, Berkunsky et al. Reference Berkunsky, Simoy, Cepeda, Marinelli, Kacoliris, Daniele, Cortelezzi, Díaz-Luque, Friedman and Aramburú2015). A previous analysis using NeoMaps data showed that even the most widespread psittacids from the Amazona genus could be experiencing negative changes in their distribution (Ferrer-Paris et al. Reference Ferrer-Paris, Sánchez-Mercado, Rodríguez-Clark, Rodríguez and Rodríguez2014). Here, we used the NeoMaps data available for all psittacid species in Venezuela to provide a complete description of the geographic distributions of most species occurring in the country. We used NeoMaps detection/non-detection records to fit occupancy models and evaluated the uncertainty and reliability of the resulting predictions for presence probabilities, and the suitability of the data to explain the lack of detections across the survey sites. Furthermore, we provide recommendations for improving future surveys to monitoring distribution changes in the country and any other Neotropical country.
Table 1. Psittacidae species reported for Venezuela. Distribution description, conservation categories, and population trend for each species according IUCN is shown. The number of NeoMaps sampling sites overlapping with the expected distribution of the species according to the available range maps from BirdLife is shown. The number of sites with detections, the ratio of current/expected detections and the number of detections for each Venezuelan psittacid species reported in Global Biodiversity Facility (GBIF) in 2010 is shown. The total number of sites sampled was 1,350.
Methods
Study species
A total of 50 species from 19 genera of psittacids occur in Venezuela, although the presence of Sun Parakeet Aratinga solstitialis has only been confirmed in the disputed territory of Guyana Esequiba (Hilty Reference Hilty2003, Rojas-Suárez pers. comm.). Seven species are endemic or almost endemic to Venezuela (Amazona barbadensis, Amazona bodini, Nannopsittaca panychlora, Pyrrhura egregia, P. emma, P. hoematotis and P. rhodocephala). Fourteen further species have restricted distribution in the country (Table 1; Hilty Reference Hilty2003).
Range maps for all these species were obtained from BirdLife (BirdLife International 2008) and clipped to the region between 0–13°N and 59–73ºW which includes all of Venezuela and neighbouring regions. Within this polygon we retrieved 21,860 presence records for all 50 psittacid species from the Global Biodiversity Information Facility (GBIF Occurrence Download http://doi.org/10.15468/dl.ofmi8y, 10 June 2016).
Field survey
The NeoMaps bird survey was performed between March and April 2010 by a team composed of seven expert ornithologists and several field assistants (methods fully described in Rodríguez et al. Reference Rodríguez, Rodríguez, Ferrer-Paris and Sánchez-Mercado2012). The sampling universe consisted of 170 half-degree cells defined in the Venezuelan Biodiversity Grid, which cover over half of the country, but do not include the southern forest regions (Figure 1). Twenty- seven cells were selected using a stratified sampling design based on environmental and biogeographical variables.
Figure 1. Sampling universe consisted in 170 half-degree cells defined in the Venezuelan Biodiversity Grid. Numbers indicate NeoMaps’ cells code visited by survey teams in 2010.
Standardised field sampling protocols for birds were implemented along a 40-km roadside transect within each cell. Two surveys were performed during two consecutive days in each transect: on the first day, 3-min point counts were performed at 50 stops, 800 m apart. On the second day, cumulative species lists were recorded at a selection of 10 stops sampled for 9 min each, divided into three consecutive 3-min periods. Total sampling effort was 108 hours of bird surveys (Rodríguez et al. Reference Rodríguez, Rodríguez, Ferrer-Paris and Sánchez-Mercado2012).
For this analysis we built detection histories for each psittacid species recorded by NeoMaps. We considered each stop as a ‘site’ (i; 1,350 sites, 50 stops across 27 transects), and each timed survey period of 3 min as a ‘observation’ (j), with duration d = 3 min. For the first day survey, detections were recorded as ‘1’ and lack of detections as ‘0’. For the cumulative list of the second day the detection history was filled with ‘0’ until the first detection, and with null values (N) afterward. Thus, valid detection histories for the second day are 1NN, 01N, 001 and 000, or NNN if the site was not visited on the second day (Ferrer-Paris et al. Reference Ferrer-Paris, Rodríguez, Good, Sánchez-Mercado, Rodríguez-Clark, Rodríguez and Solís2013). Time of day was used as an observation covariate. Sites covariates were extracted from the spatial location of each site.
Site covariates
Ranges of psittacid species are often described in terms of elevation, aridity and vegetation cover (Hilty Reference Hilty2003). Taking this into consideration, we searched for site covariates that could describe the environmental conditions during sampling and decided to use time series of remotely sensed data (Kerr et al. Reference Kerr, Southwood and Cihlar2001). In order to obtain representative data on climatic and vegetation conditions at the time of the survey, we matched the location and date of each observations with time-series of environmental variables derived from the Moderate Resolution Radio Spectrometer (MODIS) sensors in Terra-Satellites and queried them using the global MODIS Subsetting Tool (Land Processes Distributed Active Archive Center (LP DAAC 2014), and the Climate Hazards Group InfraRed Precipitation with Station data archive (CHIRPS version 2.0; Funk et al. Reference Funk, Peterson, Landsfeld, Pedreros, Verdin, Shukla, Husak, Rowland, Harrison, Hoell and Michaelsen2015). We calculated the representative value of the variable for the year prior to the sampling time (approx. March 2009–March 2010).
We considered that the effect of elevation could be properly described by the annual mean temperature, thus we used the Land Surface Temperature with Daily Cycle (LST; MOD11A2, version 6, 1-km spatial resolution) as a measure of local temperature during day time (Wan et al. Reference Wan, Hook and Hulley2015). We used the annual mean value of the Enhanced Vegetation Index (EVI; MOD13Q1, version 5, 250 m resolution) as proxy for vegetation cover, because it measures the chlorophyll concentration across all vegetation components (Didan Reference Didan2015). We used total annual precipitation and total annual potential evapotranspiration as proxies for water balance. We used the CHIRPS precipitation data (PREC; version 2, 1-km resolution) and the Potential Evapo-transpiration (PET; MOD16A2, version 6, 1-km resolution; Running et al. Reference Running, Mu and Zhao2017).
Occupancy models
We used a single-season occupancy model based on zero-inflated binomial models (MacKenzie et al. Reference MacKenzie, Nichols, Royle, Pollock, Bailey and Heines2006) to estimate the probability of occurrence for species detected in the surveys (Ψ). The occupancy state (zi) of site i was modelled as zi ~ Bernoulli (Ψi), while the observation process was modelled as yij|zi ~ Bernoulli (zi * pij) in which pij represented site and occasion specific detection probability. Covariates of Ψi (site covariates) and pij (observation covariates) were modelled using the logit link (Fiske and Chandler Reference Fiske and Chandler2011).
We fitted eight models representing different combinations of covariates for probability of detection and probability of occurrence. First, we considered models with constant probability of detection (p(.)), and others that assumed detection changed linearly with time of the day in hours (p(h)). Regarding probability of occurrence, we defined a null model with constant probability (Ψ(.)) and alternative models considering the effects of vegetation (Ψ(V), using second degree polynomials of mean EVI), climatic (Ψ(C), second degree polynomials of mean LST and total PREC and PET), or both vegetation and climatic covariates (Ψ(VC)). The models were fitted with data from sampling regions that had at least one detection or that overlapped with the expected distribution of the species according to the available range maps from BirdLife (BirdLife International 2008) and GBIF presence records for 2010 (GBIF 2018).
We evaluated the individual performance of each model using the corrected Akaike Information Criterion (AICc; Burnham and Anderson Reference Burnham and Anderson2002). Then we used the model with the best performance for each species to explain the lack of detections across the survey sites. For the sites without detections, we calculated the conditional probability of occurrence given that the species was not detected (MacKenzie et al. Reference MacKenzie, Nichols, Royle, Pollock, Bailey and Heines2006). This probability (ΨCONDL), considers two components: whether sampling effort was enough to detect the species at least once, conditional on its presence (p* = 1 - Prod (1 - p)), and the unconditional probability of occurrence given the values of the site covariates (Ψ’). We used the unmarked, raster, and AICcmodavg packages of R to fit the models (Fiske and Chandler Reference Fiske and Chandler2011).
We visualised the spatial distribution of the unconditional probability of occurrence (Ψ’) for the whole country for the species with more than 15 detections, based on the model with the highest support for each species (Table 2) and values of the vegetation and climatic covariates. We used the predict function of unmarked package (Fiske and Chandler Reference Fiske and Chandler2011) and a raster stack of predictive variables at a resolution of 1 km.
Table 2. Top performing occupancy models for 13 psittacid species with at least one detection during NeoMaps surveys. The sampling size used to fit each model is shown as the total number of sites within species range sampled during NeoMaps surveys, as well as the number of sites where each species was detected is indicated (detections). AICc = corrected Akaike Information Criterion. ΔAICc = the difference between the AIC for the i th model and the lowest AIC among all the models. AICw = relative weight from the differences in values of AICc. LL = 2log likelihood. The model with the best performance by species is in bold.
Results
NeoMaps sampling in 2010 detected 26 of the 50 species of psittacids present in Venezuela (Table 1). The most detected species were Brown-throated Parakeet Eupsittula pertinax (190 detections), Yellow-crowned Amazon Amazona ochrocephala (143 detections) and Orange-winged Amazon Amazona amazonica (97). For six species, NeoMaps sampling provided more presence records than GBIF data for the year 2010, including two detections of Orange-cheeked Parrot Pyrilia barrabandi and three detections for Maroon-tailed Parakeet Pyrrhura melanura (Table 1).
NeoMaps also provided detections outside the BirdLife distribution ranges for six species (Amazona amazonica, Amazona farinosa, Ara militaris, Diopsittaca nobilis, Orthopsittaca manilata, and Pionus menstruus. The ratio between actual and expected detections was usually lower than 10%, except for Amazona amazonica, A. ochrocephala, Chestnut-fronted Macaw Ara severus and E. pertinax, with values between 11% and 15% (Table 1). Detailed methods and results are shown in Appendix S1 in the online supplementary material.
Model fitting
The number of models fitted to each species was limited due to non-convergence or unrealistic estimates of coefficients. For the 12 species with less than five detections one or two models could be fitted. The four species with more detections also had several candidate models: Eupsittula pertinax (eight models), Orange-chinned Parakeet Amazona amazonica, Brotogeris jugularis and Blue-headed Parrot Pionus menstruus (six) and Green-rumped Parrotlet Forpus passerines (five), the rest of the species had three or four models fitted (Table 2).
Although we were able to fit models for 25 out of 26 detected species, we discarded the models for 12 species with less than five detections due to obvious over-fitting in probabilities of presence or detection. Among the 13 remaining species with reliable models, for five of them the model with the lowest AICc considered constant detectability (Table 2). For Scarlet Macaw Ara macao, Amazona ochrocephala and Eupsittula pertinax the model suggested constant high probability of detection (< 0.3; Figure 2a), while for Amazona farinosa and Pyrrhura picta, constant low probability (> 0.2; Figure 2a). For White-eyed Parakeet Psittacara leucophthalmus, Black-headed Parrot Pionites melanocephala, Military Macaw Ara militaris, Yellow-shouldered Amazon Amazonas barbadensis, Blue-headed Parrot and Orange-winged Amazon Amazona amazonica the probability of detection was low with important variation across the time, but for Brotogeris jugularis the detection was also variable but higher (Figure 2a).
Figure 2. Spatial prediction of the (unconditional) probability of occurrence for the whole country based on the model with highest support for each species (Table 2) and the values of the vegetation and climatic covariates. Darker colours indicate higher probabilities. a) Amazona amazonica; b) Amazona farinosa, c) Amazona ochrocephala; d) Brotogeris jugularis; e) Eupsittula pertinax; f) Forpus passerinus; g) Pionus menstruus.
For three species, the models with constant probability of occurrence (p(h)ψ(.)) had the lowest AICc. Models including vegetation covariates either assuming constant detectability (p(.)ψ(V)) or not (p(h)ψ(V)) had the lowest AICc for two species, Ara militaris and Painted Parakeet Pyrrhura picta. However due the low number of detections (> 10) we were not able to perform spatial prediction of unconditional probability of occurrence for these species.
Models including both climatic and vegetation covariates either assuming constant detectability (p(.)ψ(VC)) or not (p(h)ψ(VC)), were selected for five species (Table 2). The spatial prediction of unconditional probability of occurrence for Amazona amazonica; Figure 2a) showed a widespread distribution, with the highest values in the most forested and humid part of the country, in the east (south and north of the Orinoco river) and in the west (the Maracaibo Lake basin). For Southern Mealy Amazon Amazona farinose; Figure 2b) the highest probabilities of occurrence were in the dry areas along the coast, north of the country. Yellow-crowned Amazon Amazona ochrocephala; Figure 2c) and Brown-throated Parakeet Eupsittula pertinax; Figure 2e), show a widespread distribution across the country, only excluded from the Venezuelan Andes (Figure 2c,e). However, for E. pertinax, higher probabilities were predicted in the central flood plains, and the north-west, characterised by high temperature and lower vegetation cover (Figure 2e), while for A. ochrocephala higher probabilities were focused on more humid areas with moderate forest cover.
Models including only climatic covariates had the lowest AICc for two species, Orange-chinned Parakeet Brotogeris jugularis and Green-rumped Parrotlet Forpus passerinus. B. jugularis showed a restricted distribution focused on the western part of the country, in the Maracaibo Lake basin (Figure 2d), while F. passerinus, Figure 2f) had the highest probabilities in the dry areas across the coastal north.
Conditional probability of occurrence
For most species a great proportion of sites which lacked detections had low conditional probabilities of occurrence (ΨCONDL > 0.2). For Ara macao, Amazona ochrocephala, Brotogeris jugularis and Eupsittula pertinax the model estimated extremely high probabilities of detection (Figure 3a) and low probabilities of occurrence (Figure 3b), suggesting that the species are reliably detected were they are present. For Forpus passerinus, Pionus menstrus (Blue-headed Parrot), and Amazona amazonica the situation was the opposite, suggesting that the species is often present but seldom detected (Figure 3a,b).
Figure 3. Model predictions. a) Detection probability (p*) = Sampling effort required to detect the species at least once conditional on its presence. b) Conditional probability of presence given that the species was not detected (ΨCONDL).
Discussion
The observed declines in psittacid populations across Venezuela highlights the need for monitoring programmes that can reliably detect occurrence in a cost-effective and logistically feasible manner. The first national bird monitoring program in Venezuela developed by NeoMaps is an important step to achieve this aim by providing confident detection and non-detection records for 92% of psittacid species in the country (48 species with sites sampled). With this data we were able to fit occupancy models for 50% of Venezuelan psittacid species, which provide: 1) more reliable data on species absences, 2) better understanding of the importance of factors affecting psittacid occurrence, and 3) improving sampling strategy to get more confident occurrence probabilities.
Reliability of species absence
For highly mobile species such as parrots, the lack of detections is likely to be a combination of insufficient sampling effort and true absences during the time of survey (Ferrer-Paris et al. Reference Ferrer-Paris, Sánchez-Mercado, Rodríguez-Clark, Rodríguez and Rodríguez2014). For at least nine of the missing species, the lack of detection is most probably due to low sampling effort within their distribution (less than four transects, 200 sampling sites). For Aratinga solstitialis and Spectacled Parrotlet Forpus conspicillatus, no single NeoMaps survey locality matched the expected distribution (Table 2). For eight further species the expected distribution overlapped with survey localities along a single transect (less than 50 sampling sites; Amazona autumnalis, Ara militaris, Brotogeris cyanoptera, Hapalopsittaca amazonina, Pionus chalcopterus, Pyrilia barrabandi, Pyrrhura caeruleiceps and Pyrrhura melanura). All these species require targetted sampling in order to monitor their populations in the future. The most notable “absences” are those of Red-fan Parrot Deroptyus accipitrinus, Caica Parrot Pyrilia caica and Dusky Parrot Pionus fuscus, which were expected in six or seven transects (more than 300 sampling sites) south of the Orinoco river, but were never detected. For these species the presence records in GBIF during 2010 were also scarce (> 8; Table 2). Several species of Neotropical psittacids perform seasonal movements following availability of food as well as adapting to novel foods in modified environments (Juniper and Park Reference Juniper and Park1998). So in those cases of non-detections in spite of the high sampling effort, taking into account seasonal movements during the sampling design would likely improve the detection probability. Evidence from Río Manu in Peru shown a three-fold decline in the number of large macaws encountered during the dry season compared with the rainy season, which coincides with a sharp decline in plant energy production of the forest during the dry season (Renton Reference Renton2002).
For the species with low probability of detections (Figure 3a), we were able to provide useful insights to improve occurrence predictions by analysing the components of conditional probability of occurrence in those sites where the species was expected to occur, but was not detected. For example for three species (Forpus passerinus, Pionus menstrus, and Amazona amazonica), the probability of detection taking into account the current sampling effort was too low (> 0.2) in areas where environmental conditions would imply high probabilities of presence (Figure 3a,b). This suggests that sampling effort should be increased to generate reliable estimations of occurrence. Evidence from parrot communities in northern Bolivia suggest that detection probability significantly improves with a larger sampling effort (nine weeks; Berkunsky et al. Reference Berkunsky, Simoy, Cepeda, Marinelli, Kacoliris, Daniele, Cortelezzi, Díaz-Luque, Friedman and Aramburú2015). Recent empirical studies suggest that the detectability rate is significantly correlated with sighting frequency, and less conspicuous psittacids may require longer observation session in order to register the presence of some species (Rodrigues et al. 2012).
Factors affecting species occurrence
Factors related to vegetation, climatic conditions or both were important in explaining distributions of 10 species (Table 2; Figure 2). We used remotely sensed data to better represent the vegetation and climatic conditions prior to NeoMaps’ sampling period, which were heavily influenced by one episode of the ‘El Niño-Southern Oscillation’ (ENSO; May 2009 and April 2010), resulting in a severe drought, diminished water bodies and overall drier vegetation throughout the country, which in turn could have affected the probability of detection or occurrence of some species (Lentino and Portas Reference Lentino, Portas, Morales, Novo, Bigio, Luy and Rojas-Suárez1994, Hilty Reference Hilty2003). A low number of detections and subsequent low probability of occurrence of some widespread species in the Orinoco floodplains (llanos) could thus be explained by this extended drought and its consequences on vegetation growth and resource availability (Ferrer-Paris et al. Reference Ferrer-Paris, Rodríguez, Good, Sánchez-Mercado, Rodríguez-Clark, Rodríguez and Solís2013).
We were able to predict spatial distribution for seven species (Figure 2), which reflects in general the expected distribution for these species, but for Amazona farinosa (Figure 2b) the model also predicts high probabilities in the north-central region (Falcón and Lara states) where the species is absent, probably due to biogeographic constrains or the influence of additional variables not included in our model.
The low number of detections of several species often resulted in the selection of null and constant models, with low predictive power. However, for three species (Brotogeris jugularis, Amazona ochroceohala and Eupsittula pertinax; Figure 3a) with high probability of detection, better models could be obtained by including other predictive variables. This is similar to other studies where, for example, the distribution of food resources, improved the model predictions for both specialist and generalist parrot species in the Brazilian cerrado (De Araújo et al. Reference De Araújo, Marcondes-Machado and Costa2014). Topographic variables like slope and orientation also determined the distribution pattern of Mexican psittacid species, because they account for complexity in the landscape, especially at local scales (Plasencia-Vázquez et al. Reference Plasencia-Vázquez, Escalona-Segura and Esparza-Olguín2014).
Improving sampling strategy
Our study indicates that important changes in the sampling design and modeling approach are necessary to improve occurrence predictions in those sites where a species is expected to occur but is not detected. Given the low ratio of actual/expected detections (~10%), we suggest that more effective parrot monitoring programmes require: a) increasing sampling effort to improve estimates of probability of occurrence for all psittacids in Venezuela. This implies adding more sampling localities and days, and optimising survey time to those periods of the day when detection probabilities are higher; b) implementing additional surveys per year to improve estimates of seasonal patterns; c) including variables related with temporal use of resources and habitat heterogeneity in the survey and models, and d) alternatively, combining records from systematic surveys with other sources of data (collections, literature, GBIF) which could significantly increase the sample size allowing fitting of more informative models (Ferrer-Paris et al. Reference Ferrer-Paris, Sánchez-Mercado, Rodríguez-Clark, Rodríguez and Rodríguez2014).
Supplementary Materials
To view supplementary material for this article, please visit http://dx.doi.org/10.1017/S0959270919000522.
Acknowledgements
Funds for this research were provided by the Instituto Venezolano de Investigaciones Científicas (Proyect number 1071). NeoMaps data were originally curated by G. A. Rodríguez; E. Blanco and E. Goncalves helped curating GBIF records.
