Forecasting Divergence: Climate-Driven Habitat Shifts in North American Odonates Depend on Functional Groups

1. Introduction

Future anthropogenic climate change [1] will almost inevitably lead to an accelerated loss in biodiversity [2,3]. However, the effects of climate change will vary across species [4,5], making detailed analyses in specific geographic contexts necessary. Meanwhile, the exact magnitude of future climate change will vary significantly based on the choices made for global socio-economic development [2], including greenhouse gas emissions and land use patterns [6]. The Shared Socio-economic Pathways (SSPs), five qualitative narratives describing potential human development trajectories in the 21st century [7,8], outline distinct environmental scenarios [7], providing quantified data to model potential future changes in global biodiversity [7,8,9].

Arthropods are an important component of global biodiversity. Ubiquitously distributed and frequently observed, arthropods provide various important services in ecosystems, such as pollination and decomposition [10]. In future climate change scenarios, arthropods of the same species will experience climate change impacts very differently across distinct habitats [11]. Analyzing future spatial distribution of arthropods is thus relevant for understanding potential changes in their ecosystem services and overall biodiversity.

Among arthropods, odonates, which includes dragonflies and damselflies, has been described as a model taxon [12] for studying the effects of climate change because they are sensitive to both changes in air temperature, which can directly affect water temperature and influence odonate distribution [13], and changes in precipitation patterns, which can affect the structure and organization of odonate microhabitat, such as vernal ponds [14]. Odonates are also often considered a flagship group strongly associated with ecosystem health [15], monitoring and modeling their distributions can potentially reveal wider patterns for freshwater ecosystems even other arthropods. Moreover, as a charismatic taxa [5], abundant historical and current data for odonate distribution is readily available from both scientific research and amateur observation [12].

Odonates possess distinct functional traits that can contribute to wider ecosystem services [12]: as predators, they contribute to pest regulation; as aquatic species sensitive to water quality, they are indicators of freshwater ecosystem health [16]. The diversity of functional traits can contribute to the provision of a wider range of ecosystem services by directly influencing community assembly and function [17,18]. While previous studies have mainly focused on taxonomic classifications [19], this study characterizes species into functional groups based on functional traits. Compared to taxonomic groups, this slightly more flexible approach allows more meaningful interpretations, particularly on how specific communities respond to and influence the environment [3,17].

To understand how odonate functional groups will respond differently to climate change, I built MaxEnt species distribution models to analyze changes to their distributions across future years and SSPs. Species distribution models use algorithms to apply the presence data of a species on selected environmental variables to predict its geographic distribution [20]. Species distribution models have been employed mainly to guide conservation and restoration efforts [20], but could also be helpful when projecting future scenarios under climate change [5]. When an appropriate algorithm is employed [21], this more holistic method could facilitate addressing the complex, nonlinear changes in diversity and distribution [22]. My study thus conducts species distribution modeling on selected functional groups of odonates to predict changes in their habitats in varying future scenarios.

2. Methods

2.1. Defining the Scopes

Following Woods and McGarvey [23], I categorized odonate species into seven functional groups based on 11 functional traits distributed across four categories: life history, dispersal, morphology, and ecology. The seven functional groups were represented by seven specific genera [23] – Argia, Aeshna, Calopteryx, Hetaerina, Libellula, Sympetrum, and Plathemis – based on record abundance. For further details on the functional group classifications, please refer to Woods and McGarvey [23].

I chose the continent of North America, where there are currently only a few existing attempts to apply species distribution models to odonates [5], as the geographic scope. Based on data accessibility, my study roughly considered North America as all land within the range of 180^∘W to 50^∘W and 5^∘N to 85^∘N. This definition included all land masses and islands around the North American continent, except for parts of Greenland.

My study modeled odonate distributions in 2020, 2050, 2075, and 2100. The year 2020 served as a baseline for the construction of the models before they were applied to future scenarios due to the availability of data on environmental variables. The subsequent 25-year intervals were a balance between the aim of modeling gradual changes in detail and the reality of computational limitations.

The Shared Socio-economic Pathways (SSPs) provide narratives that help define the scope of future variability. To capture a wider range of estimates for future odonate distributions while balancing computational demands and limited availability of future predictive data, I selected three SSPs: SSP1-2.6, in which the world gradually shifts to sustainable development, achieving relatively low levels of climate change; SSP2-4.5, in which the world follows a steady, typical development pathway; and SSP5-8.5, in which development is accelerated by increased use of fossil fuel, leading to more drastic effects of climate change.

2.2. Retrieving and Cleaning Data

Presence-only data for North American odonates was retrieved from the Odonata Central dataset [24] hosted on the Global Biodiversity Information Facility (GBIF) platform. The download on April 26, 2025 with 248,533 records [25] was filtered to include only North American records during the period 2015–2025 to align with the defined geographic and temporal scopes, so that historical or geographically irrelevant records do not interfere with the analysis. The data were cleaned using the package pandas (Python, version 2.2.3) [26]. Since our study would explicitly focus on the seven representative genera in North America, all records irrelevant to the selected genera were removed. Since the performance of species distribution models tend to decrease at low species occurrences [21], among all 129 species across the seven selected genera, 41 of them with less than 30 occurrence records were removed. A maximum of 5 species with the most occurrence records were selected from the remaining species in each of the seven genera; notably, there were only 2 existing species in the genus Plathemis, and only 3 Hetaerina species had enough records to be included in the analysis. After unnecessary entries were removed, the remaining 59,307 records, covering a total of 30 species under the seven genera as described in Table 1, were used as presence points in the model.

To capture a range of variables that may influence odonate species distribution, I included climate, human population density, and land use changes in my analysis. Climate has been the primary determinant of species distribution at a continental spatial scale [27]. I retrieved future predictions of 19 bioclimatic variables based on the output of the HadGEM3-GC31-LL global climate model from WorldClim [28]. These predictions, with a spatial resolution of 2.5 $\min$, were presented as averages over 20-year periods ranging from 2021–2040 to 2081–2100 under three SSP scenarios. I cropped the dataset to the defined geographic scope and rescaled it to a resolution of 1 $\mathrm{k}\mathrm{m}$ using the package terra (R, version 1.8-5) [29].

Human population density could also help account for future anthropogenic influences on the habitats of odonates [5]. Future predictions of population density at 5-year intervals with a spatial resolution of 1 $\mathrm{k}\mathrm{m}$ were retrieved from [30], then cleaned to the defined geographic scope using the package terra (R, version 1.8-5) [29].

Land use changes, which have rarely been included in existing species distribution models [5], may similarly influence the habitats of odonates. Future predictions of land use at 5-year intervals with a spatial resolution of 1 $\mathrm{k}\mathrm{m}$ were retrieved from [31], transposed to ensure that the x and y axes match longitude and latitude, respectively, using the package xarray (Python, version 2025.1.1) [32], then cleaned to the defined geographic scope using the package terra (R, version 1.8-5) [29]. All land uses for irrigated crops were combined in a single new “Irrigated Crops Cumulative” category, and rainfed crops in “Rainfed Crops Cumulative”.

Selected topographic factors, such as slope, could also influence drainage basin characteristics and thus the distribution of odonates [5]. Unlike other factors, topographic data would likely have few changes throughout the defined temporal scope and were therefore modeled consistently across all the selected years. Slope and terrain roughness data with a spatial resolution of 1 $\mathrm{k}\mathrm{m}$ was retrieved from [33], then cleaned to the defined geographic scope using the package terra (R, version 1.8-5) [29].

2.3. Selecting Background Points and Environmental Variables

Since the Odonata Central database provided presence-only data, which contains information on observed species presence but not absence [34], pseudo-absence background points were randomly sampled to facilitate the construction of species distribution models [21]. Increasing the number of background points could capture a more comprehensive sample of most potential combinations of environmental factors, allowing for a better performance, but would simultaneously bring a higher computational burden [21]. Nine preliminary MaxEnt models were fitted to different numbers of background points (256 to 65,536). All retrieved variables were included without further selection. The species Aeshna sedula was selected for the preliminary model since it had 1,866 records in total, close to the mean number of records, 1,976.9, of all 30 selected species. To balance performance and runtime, 8,192 background points were used in all subsequent models.

Multicollinearity between the bioclimatic, land use, human population density, and topographic variables was assessed with variable inflation factor (VIF) test using the package usdm (R, version 2.1-7) [35]. I adopted a relatively lenient VIF threshold of 25, reflecting the large sample size and the high variance of the independent variable [36]. After 18 variables were removed, the multicollinearity of the remaining 22 environmental variables (Table 2) satisfied the selected threshold.

I conducted a second round of screening to select the main environmental variables [37,38]. For each of the 30 species (Table 1), an initial MaxEnt species distribution model with presence records and background points was fitted to all 22 remaining variables using the maxent function in the package dismo (R, version 1.3-16) [39]. All variables that did not meet a permutation importance threshold of 2.5 were dropped for the respective species; the remaining variables were included in the final models. The number of selected variables ranged from 3 for Sympetrum ambiguum to 10 for Argia fumipennis, Hetaerina vulnerata, Plathemis lydia, and Sympetrum corruptum.

2.4. Constructing Final Species Distribution Models

MaxEnt (maximum entropy) [40] is considered an efficient and well-performing option for species distribution models [21]. MaxEnt is a machine-learning technique that estimates probability distributions from incomplete information that are closest to uniform, calculated under the constraint that the expected value of each environmental variable under this estimated distribution matches the empirical average from the species occurrence data [40]. It was also specifically designed for handling presence-only data [40], suitable for the Odonata Central database [24].

The 8,192 background points were randomly sampled from the defined geographic scope, then combined with the existing presence-only data for each of the 30 species. The combined data was randomly divided through an 80:20 train-test split using the folds function in the package predicts (R, version 0.1-17) [41]. The larger training set was used to construct a final MaxEnt model based on the environmental variables selected through the initial models for the baseline year of 2020 using the maxent function in the R package dismo (R, version 1.3-16) [39]. The smaller testing set was used to evaluate the final model through the area under the receiver operating characteristic curve (${\mathrm{AUC}}_{\mathrm{ROC}}$), a number within the range $(0,1)$, using the pa_evaluate function in the package predicts (R, version 0.1-17) [41]. A random prediction would result in ${\mathrm{AUC}}_{\mathrm{ROC}}=0.5$, while ${\mathrm{AUC}}_{\mathrm{ROC}}\ge 0.75$ would indicate that the model was well-fitted [42]. The 30 models for the species had a mean testing ${\mathrm{AUC}}_{\mathrm{ROC}}=0.968$ with a standard deviation of $0.017$, suggesting that all models were well-fitted.

The models were then extrapolated to each of the year-SSP combinations in the defined temporal and variability scopes. The output for each estimation was given in the form of a continuous estimate of suitability from 0 to 1; a higher value indicates greater suitability. The average estimate of suitability throughout the defined geographic scope was recorded as the overall habitat suitability.

A threshold was used to convert the model’s continuous estimate into a binary prediction. The maximum sum of specificity sensitivity (maxSSS) threshold would find the highest value of the sum of specificity (true presences) and sensitivity (true absences) [43]. This threshold has been proven effective in previous studies [44,45] and was especially suitable for this scenario due to the abundance of presence data points [44]; though traditionally used on presence-absence data, it has demonstrated good performance on presence-only data [43]. In the resulting binary prediction, range size was defined as the area of all 1$\mathrm{k}\mathrm{m}$ grids predicted as the habitat. The centroid, a coordinate pair, was calculated as the mean coordinate of all grids marked as the predicted habitat. Range size and centroid were recorded for each species under each year-SSP combination.

In total, each model returned three key outputs for analysis. Overall habitat suitability was the mean of the direct predictive output of the species distribution model – a value between 0 and 1 – under the defined geographic scope. Range size was the total area of a species’ predicted habitat obtained by applying the maximum sensitivity specificity sum (maxSSS) threshold to the direct output of the species distribution model. Centroid change was the Euclidean distance between the geometric centroid of the predicted habitat of a species in a specific future scenario and the geometric centroid of the habitat in the baseline year of 2020.

2.5. Analysis

Changes between each of the predicted future scenarios and the baseline scenario of the respective species were calculated based on the raw outputs of the MaxEnt models. Percent change from baseline was calculated for overall habitat suitability and range size, while Euclidean distance from baseline was calculated for centroid.

A linear model was fit to the log-transformed response variables using lm function from stats (R, version 4.4.3) [46], with year as an independent variable alongside SSP and genus as interacting terms. Logarithmic transformation was applied to account for the non-negative, right-skewed nature of the data. Since the models on overall habitat suitability and range size did not meet all assumptions of residual normality, we performed non-parametric bootstrapping using the package boots (R, version 1.3-31) [47] to evaluate the stability of coefficient estimates and confidence intervals. Most of the resulting confidence intervals of the coefficients were within acceptable ranges to validate the linear models, though slight deviations for the interacting factors of SSP2:Hetaerina and SSP2:Plathemis, presumably due to the lower number of selected species in the respective genera, suggests that these particular results should be interpreted with more caution.

Marginal means for each genus and each SSP scenario were then calculated using the R package emmeans [48], allowing direct comparison across genera and SSPs. Statistical significance was evaluated through a non-parametric assessment on the 95% confidence intervals of the marginal means [49]. Furthermore, the permutation importance of each environmental variable, calculated as the normalized percent decrease in training AUC after random permutation of its values [50], was retrieved from the final models to assess which environmental variables were more significant at influencing odonate habitats; variables excluded from the initial models, and thus unused in the final models, were considered to have a permutation importance of 0.

3. Results

3.1. Overall Habitat Suitability Changes

All seven functional groups showed an increase in overall habitat suitability, which reflects whether the environment across North America is generally more or less conducive to the survival of a species (Figure 1), with six groups being statistically significant. Aeshna ($\overline{x}=+28.830\%$, ${\mathrm{CI}}_{95\%}=[+22.902\%,+35.044\%]$) will experience the largest increase. Sympetrum ($\overline{x}=+9.178\%$, ${\mathrm{CI}}_{95\%}=[+4.155\%,+14.444\%]$), Hetaerina ($\overline{x}=+8.527\%$, ${\mathrm{CI}}_{95\%}=[+2.123\%,+15.331\%]$), and Calopteryx ($\overline{x}=+8.607\%$, ${\mathrm{CI}}_{95\%}=[+3.609\%,+13.845\%]$) shared similar and overlapping predictions. The exception is the functional group represented by Libellula ($\overline{x}=+4.441\%$, ${\mathrm{CI}}_{95\%}=[-0.365\%,+9.478\%]$), whose increase in overall habitat suitability is not statistically significant, suggesting a potential future risk of competitive disadvantage, vulnerability, and the need for conservation focus compared to the other species.

The effects of SSP scenarios varied across genera (Figure 2). SSP1, which describes more sustainable development, displayed the smallest differences between functional groups, consistently leading to slight increases ranging from $+3.923\%$ (Sympetrum) to $+14.805\%$ (Aeshna). By contrast, SSP5, which describes an accelerated climate change, led to the largest differences between predicted changes and current suitability, resulting in a considerable increase in overall habitat suitability for functional groups such as those represented by Aeshna ($+53.438\%$) and Hetaerina ($+38.027\%$). In some functional groups, however, the effects of SSP was not significant, such as those represented by Calopteryx (SSP1 $+10.295\%$, SSP5 $+10.194\%$) and Libellula (SSP1 $+3.607\%$, SSP5 $+3.330\%$).

3.2. Range Size Changes

Range size change across functional groups, which reflects whether more or less area is considered favorable for the genera, similarly increased (Figure 3). All functional groups but Libellula ($\overline{x}=-2.969\%$, ${\mathrm{CI}}_{95\%}=[-8.250\%,+2.616\%]$) were predicted to have a greater range size, from Sympetrum ($\overline{x}=+0.342\%$, ${\mathrm{CI}}_{95\%}=[-5.119\%,+6.117\%]$) to Aeshna ($\overline{x}=+26.998\%$, ${\mathrm{CI}}_{95\%}=[+20.086\%,+34.308\%]$). For three of the seven functional groups (Libellula, Sympetrum, and Calopteryx), the predicted increase in range size was not statistically significant.

Varying effects of SSP on range size across functional groups also resembled the results on overall habitat suitability (Figure 4). SSP1 generally displayed the smallest differences between functional groups, leading to slight changes ranging from $-2.957\%$ (Sympetrum) to $+20.119\%$ (Aeshna). SSP5 led to the largest differences between functional groups, with Aeshna increasing ($+50.672\%$), and Libellula decreasing ($-10.116\%$).

3.3. Centroid Change

Centroid change, or how far a species has migrated from its original habitat, varied significantly among species (Figure 5). The functional group represented by Aeshna was predicted to experience the greatest centroid change ($\overline{x}=390.781$, ${\mathrm{CI}}_{95\%}=[320.402,476.619]$), with a statistically significant difference compared to all the other functional groups, followed by a more moderate change experienced by the functional group represented by Sympetrum ($\overline{x}=[212.138$, ${\mathrm{CI}}_{95\%}=173.932,258.735]$). The functional groups represented by Libellula, Plathemis, and Argia were the least affected with a centroid change of $<100$ $\mathrm{k}\mathrm{m}$; they also had the narrowest confidence intervals, suggesting prediction robustness.

The effects of SSP on centroid were consistent across functional groups (Figure 6). SSP5 led to the highest change in centroid, with a change of $1,008.547$ $\mathrm{k}\mathrm{m}$ for Aeshna. SSP2, the intermediate scenario, typically led to a moderate change, ranging from 48.069% (Libellula) to 71.075% (Hetaerina) compared to the changes shown in SSP5 scenarios. SSP1 resulted in the lowest changes overall, such as merely $45.127$ $\mathrm{k}\mathrm{m}$ for Plathemis, again suggesting how this more sustainable scenario could bring more certainty and assurance to future conservation efforts.

Table 3. Summary of marginal means and confidence intervals, grouped by genus, for overall habitat suitability (adjusted by $-100\%$), range size (adjusted by $-100\%$), and range centroid distance change.

Genus	General Suitability (%)		Range Size (%)		Centroid Distance ($\mathrm{km}$)
	$\overline{x}$	${\mathrm{CI}}_{95\%}$	$\overline{x}$	${\mathrm{CI}}_{95\%}$	$\overline{x}$	${\mathrm{CI}}_{95\%}$
Aeshna	$+28.830$	$[+22.902,+35.044]$	$+26.998$	$[+20.086,+34.308]$	$390.781$	$[320.402,476.619]$
Argia	$+13.406$	$[+8.188,+18.876]$	$+10.328$	$[+4.323,+16.678]$	$88.692$	$[72.719,108.174]$
Calopteryx	$+8.607$	$[+3.609,+13.845]$	$+4.550$	$[-1.141,+10.567]$	$113.658$	$[93.188,138.624]$
Hetaerina	$+8.527$	$[+2.123,+15.331]$	$+7.743$	$[+0.233,+15.815]$	$136.880$	$[105.927,176.877]$
Libellula	$+4.441$	$[-0.365,+9.478]$	$-2.969$	$[-8.250,+2.616]$	$92.167$	$[75.568,112.413]$
Plathemis	$+16.842$	$[+8.456,+25.877]$	$+11.122$	$[+1.712,+21.402]$	$86.727$	$[63.358,118.717]$
Sympetrum	$+9.178$	$[+4.155,+14.444]$	$+0.342$	$[-5.119,+6.117]$	$212.138$	$[173.932,258.735]$

Table 3. Summary of marginal means and confidence intervals, grouped by genus, for overall habitat suitability (adjusted by $-100\%$), range size (adjusted by $-100\%$), and range centroid distance change.

Genus	General Suitability (%)		Range Size (%)		Centroid Distance ($\mathrm{km}$)
	$\overline{x}$	${\mathrm{CI}}_{95\%}$	$\overline{x}$	${\mathrm{CI}}_{95\%}$	$\overline{x}$	${\mathrm{CI}}_{95\%}$
Aeshna	$+28.830$	$[+22.902,+35.044]$	$+26.998$	$[+20.086,+34.308]$	$390.781$	$[320.402,476.619]$
Argia	$+13.406$	$[+8.188,+18.876]$	$+10.328$	$[+4.323,+16.678]$	$88.692$	$[72.719,108.174]$
Calopteryx	$+8.607$	$[+3.609,+13.845]$	$+4.550$	$[-1.141,+10.567]$	$113.658$	$[93.188,138.624]$
Hetaerina	$+8.527$	$[+2.123,+15.331]$	$+7.743$	$[+0.233,+15.815]$	$136.880$	$[105.927,176.877]$
Libellula	$+4.441$	$[-0.365,+9.478]$	$-2.969$	$[-8.250,+2.616]$	$92.167$	$[75.568,112.413]$
Plathemis	$+16.842$	$[+8.456,+25.877]$	$+11.122$	$[+1.712,+21.402]$	$86.727$	$[63.358,118.717]$
Sympetrum	$+9.178$	$[+4.155,+14.444]$	$+0.342$	$[-5.119,+6.117]$	$212.138$	$[173.932,258.735]$

3.4. Environmental Variable Permutation Importance

Permutation importance of environmental variables were collected from the final models to evaluate their influence on odonate distributions (Figure 7). The results suggested that bioclimatic variables, especially temperature-related variables like Annual Mean Temperature (BIO1, $\overline{x}=49.500\%$, $\mathrm{SE}=3.935\%$) were a particularly strong determinant of odonate distributions, followed by Annual Mean Precipitation (BIO12, $\overline{x}=13.931\%$, $\mathrm{SE}=2.600\%$) and Temperature Annual Range (BIO7, $\overline{x}=11.287\%$, $\mathrm{SE}=1.351\%$). Selected land use types like Needleleaf Evergreen Tree Boreal (PFT2, $\overline{x}=5.938\%$, $\mathrm{SE}=1.060\%$) and Human Population Density (HPD, $\overline{x}=5.267\%$, $\mathrm{SE}=1.235\%$) also had relatively significant contributions. However, most land use types, as well as topographic slope ($\overline{x}=2.677\%$, $\mathrm{SE}=0.550\%$), had much smaller influence overall, suggesting that they had low relevance to odonate species distribution.

4. Discussion

The results of MaxEnt and linear models revealed that the aggregate habitat of all odonate functional groups will shift to the north and expand into smaller pockets, with the greatest changes observed under SSP5. Among the environmental predictor variables, bioclimatic variables such as Annual Mean Temperature (BIO1) showed the most significant impact. Significant differences were observed between functional groups, indicating a functional approach to species distribution modeling under global change scenarios may be an effective way to understand potential changes to ecosystem services.

The statistically significant predictions on distance shift indicate that all odonate species will migrate to the north (polewards) in the coming decades. This kind of habitat shift may be feasible for odonates since North America does not have major east-west geographic constraints, such as the Alps in Europe, that would prevent such movement [12]. However, the habitat shift of odonates may trigger chain reactions leading to habitat shifts of multiple related species in a shared ecosystem [51], creating new challenges and opportunities for species conservation of and beyond odonates.

Notably, stronger habitat changes compared to the baseline scenario and greater differences between functional groups are observed under SSP5, the trajectory with more severe global warming. All species are expected to shift further to the north. Species whose overall suitabilities are predicted to increase will have even larger range sizes, while species with lower overall suitabilities will likely require more conservation attention. This highlights the uncertainty that a less sustainable pathway might bring to future biodiversity and conservation efforts, underscoring the necessity to reduce carbon emissions to meet the 2 degrees climate target [52].

As the habitat of species shifts to the north in future decades, many will encounter changes in landforms that split a species into relatively isolated smaller fragments, potentially leading to a loss in genetic diversity within a species and eventually higher extinction risks [53]. This is supported by the model’s predictions that range size is predicted to increase slower than overall habitat suitability. As a case study, the predicted distribution changes for Calopteryx angustipennis, or Appalachian Jewelwing, in 2100 under SSP5 shows a drastic shift to the north from its baseline distribution in southern United States to enter Canada (Figure 8). Compared to the current large cluster in the southeastern United States, the species will likely split into four smaller fragments, occupying Nova Scotia, the tri-state area, the Great Lakes, and the peaks of the Appalachian Mountains.

Of all the environmental variables used to predict future habitat changes, the high permutation importance of all temperature variables lends further support to the projected impacts of climate on odonates [5,12]. On the other hand, slope and many land use types had low permutation importance. Future efforts at modeling species distribution should focus more on climatic variables, but can potentially keep the few most relevant land use types for odonates in the region, such as Needleleaf Evergreen Tree Boreal (PFT2) in North America, as supplementary determinants.

The models and analyses have some limitations worth further consideration. Geographically, the isolation between the eastern and western sides of North America can lead to overprediction of species distributions. Since major north-south features in North America, such as the Rocky Mountains and the Great Plains, can prevent east-to-west migration, a species from one side of North America would be unable to colonize a region on the other side even if the climate and land use patterns are similar. For instance, though C. angustipennis is currently only found in the southeastern states [24], the baseline model falsely suggested an additional tiny piece of C. angustipennis habitat in central Mexico, where similar environmental variables are present (Figure 8). The model could thus potentially overpredict overall habitat suitability and range size and shift all range centroids to the west. However, this consideration applies to all predictions, including the baseline year and in future scenarios. Since the final analyses compare the future changes to the baseline year, the generalizable results still reflect relatively accurate changes throughout the years. Additionally, the uneven representation across genera, from 2 for Plathemis to 5 for the other genera, resulted in stronger uncertainties in functional groups represented by less species, limiting the strength of interpretation.

Despite potential limitations, this study’s results can help to design conservation strategies. Of the seven representative genera included in the study, species with negative predicted range size changes and higher habitat fragmentation should be considered more vulnerable and ought to be prioritized in conservation efforts. Such conservation efforts can be applied to other genera in the same functional group. Alternatively, the more vulnerable functional groups identified by my study can be primary targets for additional species distribution modeling to further specify the variations within the functional groups.

Though the study was limited to North America, it would be important to project the methods and findings of the study to a global scope. Species distribution modeling has only been applied to less than a quarter of all discovered odonate species [5]. Since the current study revealed that responses to climate change could vary significantly among functional groups, more modeling attempts would be necessary to discover patterns for different functional groups in other spatial contexts. Notably, the quality of presence data in other spatial contexts can be vulnerable to sample selection bias, in which more presence observations are recorded in easily accessible areas [55], which can be mitigated by extracting such bias from the presence-only dataset to be intentionally introduced when selecting background points [55]. Meanwhile, the study can also include other potential environmental variables for more robust predictions. For instance, data describing river characteristics such as water chemistry [5] and agrochemical usage [56] could impact earlier stages of odonate life cycle, especially naiads, and thus influence the distribution of adult odonates [5]: potential anthropogenic damages to water quality, for example, might reduce future habitats for odonates. Beyond modeling with environmental variables, some studies have highlighted the interactive effects of flora and vegetation on odonate abundance [57]; however, others have found weaker correlations between plant diversity and odonate distribution, instead recommending the current approach[58]. Expanding the scope of the study could help identify the causes of nonlinear variations in predicted biodiversity changes, facilitating the conservation of important practical and cultural values in diverse ecosystems [3].

By modeling the distribution of North American odonate functional groups using MaxEnt, this study took a novel approach to categorizing species for species distribution modeling and provided insights into the future distribution of odonates under climate change scenarios. In particular, most functional groups will experience northward range expansions, increased habitat suitability, and greater fragmentation, with effects amplified under SSP5 and Libellula experiencing potential range contractions. These findings emphasize the need for targeted conservation strategies tailored to each functional group under different climate scenarios, especially in functional groups predicted to suffer from habitat fragmentation. Future research can consider adapting the functional-group-based approach to other taxonomic groups and geographic regions to provide critical insights for conservation planning that considers both species vulnerability and cascading effects on ecosystem services.