Introduction

The biodiversity crisis that has been documented in detail for several taxa (^{Dirzo et al., 2014}), includes insects of the order Odonata (dragonflies and damselflies) with 1 in 10 odonate species being threatened (^{Clausnitzer et al., 2009}). This despite their fairly high dispersal ability and flexibility to occupy different habitats (^{Stoks & Córdoba-Aguilar, 2012}). Similar to other aquatic insects, deforestation and susceptibility of anthropogenic change for their freshwater habitats (especially lentic waters; ^{Clausnitzer et al., 2009}) seem major drivers underlying threat (^{Dijsktra et al., 2014}).

Argia damselflies are a highly speciose genus with up to 144 species (^{Caesar & Wenzel, 2009}; ^{Garrison, 1994}), with a large occurrence throughout the American continent. Some generalities of their ecology are that they use both lentic and lotic freshwaters, and occur in different habitats that comprise tropical and non-tropical areas (^{Garrison, 1994}). Argia species are mostly small-sized odonates with males of some species being either territorial (^{Guillermo-Ferreira & Del Claro, 2012}) or nonterritorial (^{Sánchez-Guillén et al., 2014}). The International Union for Conservation of Nature (from now on, IUCN) has incorporated 36 species of Argia in the Red List, of which only A. sabino appears threatened (“vulnerable” category; ^{IUCN, 2017}). This reduced number implies that further efforts are badly needed to uncover more Argia species. One criteria to be included as part of IUCN assessment is the extent of occurrence (criterion B; ^{IUCN, 2017}) which can be gathered using ecological niche models (^{Guisan & Zimmermann, 2000}; ^{Gouveia et al., 2011}).

With the aim of providing biological information of odonates that have not been included in the IUCN Red List, in this work we have assessed the distribution area of 6 Argia species: Argia cuprea (Hagen, 1861), A. funcki (Selys, 1854), A. garrisoni Daigle, 1991, A. harknessi Calvert, 1899, A. munda Calvert, 1902, and A. rhoadsi Calvert, 1902. Except for A. garrisoni which appears as least concern (^{IUCN, 2017}), the remaining 5 species have not been assessed.

Materials and methods

Presences of the 6 Argia species were compiled from literature records, confirmed records in OdonataCentral (http://www.odonatacentral.org), Conabio records (www.conabio.gob.mx), GBIF records (www.gbif.org), and published papers, thesis and data provided by experts in odonates. All data were checked carefully for geographic accuracy which, when in doubt, included asking specialists (e.g., Rosser Garrison). Notice that most records were gathered from OdonataCentral, which is a repository of odonate record information. We are confident that the biases of identification should be minimal because these data were provided by experts for these animals. Still, we limited our analyses to the 6 species indicated above as these were the ones we were fully confident in regards to the quality of their records.

Distribution models were built only for species with more than 30 localities. Thus, the final data set included 1585 unique presences of the 6 species, which were predominantly adult specimens: 759 records for A. cuprea, 261 for A. funcki, 255 for A. garrisoni, 185 for A. harknessi, 86 for A. munda and 39 for A. rhoadsi. Although these records spanned the last 40 years, most came from the last 15 years. The database is available upon request.

We used the WorldClim 1.4 (www.worldclim.org) data set (^{Hijmans et al., 2005}) at 0.041666669 degrees (or 5 km) of resolution to set predicting bioclimatic variables for our models. To establish a set of uncorrelated climatic variables, we intersected the variables with target group points, and with 10,000 points randomly selected in the extension of the study area (M). We eliminated some variables with an exploratory data analysis of contribution of variables using jackknife and Pearson correlation analyses (i.e., any value > 0.7). Thus, we selected variables with low correlation and high contribution to reduce the parametrization of the models. After this, the final data set included uncorrelated variables which had more biological importance for Argia species, and more contribution according to the jackknife analysis. Such variables were: mean diurnal temperature range (bio 02), temperature seasonality (bio 04), mean temperature of warmest quarter (bio 10), precipitation of wettest month (bio 13), precipitation of driest month (bio 14), and precipitation seasonality (bio 15). Our study area included North America where our study species occur.

Our study area was North and Central America where members of the genus Argia occur (^{Caesar & Wenzel, 2009}; ^{Garrison, 1994}), covering land between the latitudes 53.00 to 0.00N, and the longitudes -130.00 to -55.00W. To choose the best background, preliminary species distribution models were generated with Maxent 3.3.3k (^{Phillips et al., 2006}) with target group points, 10,000 points randomly selected in the extension of the study area (M), and with special extent delineating M for each particular species with ecoregions (World Wildlife Fund; www.worldwildlife.org/) for North America, and biogeographical provinces (Conabio) for Mexico. Models were constructed setting several parameters to default (‘Auto features’, convergence = 10-5, maximum number of iterations = 500). However, we used random seed (with a 30 test percentage), 10 replicates, removed duplicate records, ran bootstrap replicated type, with no extrapolation and no clamping. All this to find which combination of settings and variables generated the best outcomes (highest area under the curve, or AUC) while minimizing the number of model parameters, as well as producing ‘closed’, bell-shaped response curves guaranteeing model calibration and transferability (^{Elith et al., 2010}). The best background was 10,000 points randomly selected in the extension of the study area.

Final models were built with Biomod (Biodiversity Modelling) package in R software. This package is a platform for predicting species’ distribution, including the ability to model the distribution using various techniques and test patterns (^{Thuiller et al., 2009}). When using Biomod, we trained models using 4 widely used algorithms: maximum entropy (Maxent), random forest (RF), generalized boosting methods (GBM) and multivariate adaptive regression splines (MARS). These models have shown good performance in terms of predictive power (^{Broennimann et al., 2012}; ^{Pliscoff & Fuentes-Castillo, 2011}; ^{Reiss et al., 2011}). From individual models obtained with these different algorithms, we generated a “consensus model”. Such model combination is the best logistic compromise to avoid either overfitting and overpredicting (^{Merow et al., 2014}). In other words, this reduces biases and limitations of using only individual models. 70% of data was used for training, and 30% for validation with 10 replicates. The final validation of models was performed with TSS (True Skills Statistics), average net rate of successful prediction for sites of presence and absence (^{Liu et al., 2009}), ranging from -1 to 1, where the more positive values indicate a higher degree of accuracy and discrimination model (^{Allouche et al., 2006}) (Table 1). It is noteworthy that the result of these models is not the area that species occupies absolutely, because these models do not consider population dynamics, dispersibility, interactions with other species and human impacts. However, these models can make right transferences where species can be potentially found given their environmental conditions. This is based on the assumption that the distribution known of species provides enough information to characterize its environmental requirements.

A total of 420 models were generated, whose performance was assessed by means of the AUC and TSS statistic, while minimizing the number of model parameters. The best presence/absence models using the “10 percentile training presence” are presented. This threshold was used because we prefer to err on the side of caution accepting that a 10% of our presences could be problematic (for a similar rationale, see ^{Sánchez-Guillén et al., 2013}). The best models of current climatic conditions of species were used to generate projections into the geography.

Results

Our results indicated that A. garrisoni showed a distribution area of less than 20,000 km² (see below). Since a distribution area smaller than 20,000 km² is one of those criteria that allows categorizing whether any species may be threatened (^{IUCN, 2017}), we took a step further to characterize the habitat of this species for any field assessment and/or future protection effort (for a similar approach, see ^{Cuevas-Yáñez et al., 2015}). For this characterization, all types of land use corresponding to the habitat of the species were selected and trimmed using both the soil cover layer and the vegetation cover of the Instituto Nacional de Estadística y Geografía (Inegi) (Union of Layers, Series V: ^{Inegi, 2011}), removing, for example, urbanized or arid zones that do not correspond to A. garrisoni’s habitat, rendering a “clean” layer with this selection. From the overlap of this layer on the distribution map of the species obtained in the Biomod package, a reduced distribution area was provided compared to the initial distribution, but more adjusted to the actual location of the species. The map shown in the results is the one obtained after making the land use cut-off corresponding to the habitat for this species. We finally obtained the area in km² by calculating the number of cells to a size of 0.041666669, with ArcMap® software.

Ecological niche models predicted a considerably large distribution for all species except for A. garrisoni which had 9,991 km² (Fig. 1; Table 1). Our trimming exercise to predict a more suitable habitat for this species, rendered a more reduced area of 8,038 km²

Figure 1 Distribution estimates and maps generated according to ecological niche modeling for 6 Argia species. 

Discussion

Five out of 6 Argia species that we modelled here showed fairly large distribution areas (over 70,000 km²). This large area estimates for these 5 species coincides with those gathered for other Argia species such as A. anceps, A. extranea, A. oenea, and A. tezpi that mainly occur in central Mexico (^{Nava-Bolaños et al., 2017}). It seems that the niche properties of the entire genus are highly conserved as assessed by their large degree of area overlapping (i.e., sympatry; ^{Nava-Bolaños et al., 2017}). Such similar niche properties imply frequent copulations with heterospecifics. This is the case of A. rhoadsi and A. munda which are found mainly in northeastern Mexico and whose distribution extends from the southern part of Texas to the lower part of San Luis Potosí (^{Dunkle, 2004}), inhabiting lagoons, streams and riversides (^{Abbott, 2001}). Both species can be frequently found in heterospecific matings (A. Córdoba-Aguilar, unpub. obs.). As a matter of fact, possibly the radiation of the entire genus may not necessarily be related to their ecological pressures but to sexual selection forces as has been proposed for other damselflies (^{Wellenreuther & Sánchez-Guillén, 2016}). This non-adaptive radiation may render species not to be habitat-specific. Perhaps the exceptions are those species that are at risk such as A. sabino (^{IUCN, 2017}) and A. garrisoni.

Our results can be used to fulfill current IUCN Red List knowledge. As indicated before, no risk assessment has been carried out for A. cuprea, A. funcki, A. harknessi, A. munda and A. rhoadsi. Given the large area of occurrence estimation, these 5 species should be placed in a low risk category by not meeting the criteria of a reduced extension of occurrence. However, this is not the case of A. garrisoni, which is an endemic species from Mexico (^{González-Soriano & Novelo-Gutiérrez, 2014}). Despite the reduced distribution estimated by our methodology, this species was abundant in the Mexican state of San Luis Potosí in 1999, in the roots of citrus along rivers and perched in banks of partially dry sludge (R. Garrison pers. comm.). Ecological niche models suggested that this species would occur in the following biogeographic provinces: Sierra Madre Oriental, Oaxaca, and Gulf of México (^{Conabio, 1997}). However, since 1999 there have been no more recent records for this species. According to the IUCN criteria, the fulfillment of even a single criterion, in this case the B1, allows the change from least concern to vulnerable which is when an area of less than 20,000 km² is estimated (^{IUCN, 2017}). Even when one may see such criteria as too vague, restricted distribution areas correlate well with the probability of extinction in odonates (the smaller the area, the more likely that a species will disappear; ^{Korkeamäki & Suhonen, 2002}), so our risk concerns for A. garrisoni are likely to be founded.