Predicting Habitat and Distribution of an Interior Highlands Regional Endemic Winter Stoneﬂy ( Allocapnia mohri ) in Arkansas Using Random Forest Models

: Stoneﬂies are a globally threatened aquatic insect order. In Arkansas, a diverse group of winter stoneﬂy (Capniidae: Allocapnia ) have not been surveyed since the 1980s, likely because species-level identiﬁcation requires the rarely-collected adult form. Allocapnia mohri , a regional endemic, was previously commonly found in mountainous, intermittent streams from the Ouachita Mountains ecoregion north to the Ozark Highlands, but no species distributional models including land use or climate variables exist to our knowledge. We collected adults from 71 stream reaches from the historic Arkansas range from November to April 2020 and 2021. We modeled distributions using random forest (RF) models populated with landscape, climate, and both data to determine which were most predictive of species presence. Correlations between landscape or climate variables and presence were examined using multiple logistic regression. The landscape RF models performed better than the climate or landscape + climate RF models. A. mohri presence sites tended to have a greater elevation, a lower mean July temperature, and a greater percentage of very slow inﬁltration soils in the watershed, compared to absence sites. A. mohri was absent at the Ouachita Mountains sites and may be experiencing a range contraction or migration northward.


Introduction
Stoneflies (Plecoptera) are a highly diverse macroinvertebrate group encompassing more than 1000 species and representing 99 genera and ten families in North America [1]. Though the plecopteran group is expansive, little species-specific data on habitat, feeding habits, and distributions exist, and several existing records need updating to more accurately estimate their population status and to manage their conservation. This information is essential as Plecopterans are estimated to be one of the most endangered groups of insects globally [2,3]. Possible reasons for this include the alteration of plecopteran habitat by anthropogenic influences, such as climate and land use changes [3]. Here we model distributions of an endemic stonefly species to close the gap in species-specific stonefly knowledge, and to investigate whether climate, landscape, or a combination of factors could be affecting spatial distributions of stoneflies.
Several insect species, including the nymphs of some species of stonefly, are reliant on allochthonous resources [4], and the energy transported by these terrestrial materials, along with the energy acquired through autochthonous resources, is the primary source of matter and energy for many streams, particularly low order forested streams [5][6][7]. Shredding stoneflies will consume autumn-shed leaves, then these nymphs are subject Hydrobiology 2023, 2 197 to predation, thereby using the energy acquired through allochthonous sources. [8]. This cycling in of allochthonous material into the stream makes stoneflies a vital component of the stream food web [9]. The reliance on fallen leaves for food adds a habitat requirement of course particulate organic matter (CPOM) for shredder stoneflies in this macroinvertebrate group, and several stonefly species have strict water oxygenation and temperature requirements [10,11], increasing the habitat needs for these sensitive plecopteran nymphs. One such group of threatened plecopterans is the winter stonefly.
Stoneflies require specific environmental conditions, such as stream water temperatures and hydroperiod [12], and these factors will be affected by global climate change [13]. Ephemeroptera, Plecoptera, and Trichoptera (EPT) richness has displayed a decline of 81.6% worldwide, coming to an average loss of 1.9% loss of EPT richness each year. This 1.9% loss in richness mirrors the 1.88 • C increase in water temperature over the 40 years [14]. Several stonefly species have a thermal cue for egg hatching, and if stream temperatures do not reach a thermal low-end cue, nymphs will not hatch, causing extinction of species in certain thermally inadequate streams [12]. The effect of changing hydroperiod can greatly affect stoneflies, as they spend most of their lives as aquatic nymphs. Plecopterans are sensitive to benthic habitat loss, due in part to their need to lay eggs inside refugia within the stream (i.e., aquatic plants and woody debris), and their need to emerge using this refugia [12]. In some streams, stonefly nymphs can retreat to the hyporheic zone in dry periods, and it was found that several stonefly species even have diapausing eggs to outlast the dry stream conditions [15], though for many species and for many systems the hyporheic zone is not as important as refugia. Small stoneflies have been found to aggregate on stones, and larger stoneflies aggregate in leaf patches during high flow, as well as a general movement toward areas with flow in drying summer months [16].
Land-use change is another key driver of declines in stonefly populations [3]. Loss of riparian forest will likely amplify increases in stream water temperatures [17], reducing forested area in catchments, which may reduce stonefly habitat availability [18], and this, along with climate change, may cause reductions in stonefly populations due to their strict thermal requirements [10]. Further, standing crops of coarse benthic organic detritus can be directly related to the secondary production of shredding macroinvertebrates, such as stoneflies [19,20]. So, the loss of forested land in catchments can be strongly linked to a lower standing crop of litter detritus and likely reduced food availability for shredding stoneflies.
Forested land is often converted to agricultural or urban land uses that may also contribute to stonefly population declines and local species extinctions. Anthropogenic activities and land usage around stream systems directly influence the physical and chemical characteristics of the streams [21], including the amounts of sewage inputs, salts from road de-icing, and agricultural fertilizer and pesticides that seep into streams [22][23][24]. Anthropogenic activities can also result in impervious soil surrounding streams, leading to slow infiltration of water [25], causing an increase in speed and volume of storm runoff [26], and an increase in pollutants as a whole, ultimately causing a degradation of the aquatic habitat [27] and a loss of benthic macroinvertebrate diversity [28].
Anthropogenic nutrient inputs often lead to an increase in primary production [29], and this eutrophication can lead to low levels of dissolved oxygen (DO) during certain times of the day [29], thereby killing sensitive species in the stream, such as stoneflies [30]. This increase in nutrients can also stimulate microbial decomposition of leaf litter, the food source of shredders such as stoneflies [31]. While microbial colonization is thought to improve the nutritional quality of leaf litter for shredders, this is not always the case [32][33][34][35]. Further, elevated microbial decomposition may reduce the quantity of coarse detrital carbon available for food in stream ecosystems [36]. These effects can compound over time and be detrimental to stonefly species in the long-term.
Plecopterans are sensitive to these stream impurities and to microhabitat disturbances which makes them an important biological model to indicate water quality [37]. Most of the studies on benthic macroinvertebrates in relation to water quality evaluations have dealt with the richness of the EPT community because these groups are intolerant to pollutants in water bodies, but in the EPT community, Plecopterans are often the most sensitive to changes in water quality [38].
While the common causes of global-scale stonefly declines are known, a finer-scale, regional perspective is needed for conservation of sensitive groups, such as winter stoneflies. Winter stoneflies (Order: Plecoptera) emerge from the nymph stage in the winter months, and they are known as "clean-water insects" due to their need for highly oxygenated water [39]. A genus of winter stoneflies (Family: Capniidae, Genus: Allocapnia) emerges from the nymph stage to a small-winged adult across a wide variety of intermittent to perennial streams on vertical surfaces in and around the water body and lacks stream-tostream vagility [12]. The Arkansas Level III Ecoregions known as the Ozark Highlands, Boston Mountains, Arkansas Valley, and Ouachita Mountains have a high diversity of Allocapnia. During the last regional census in the late 1980 s, seven endemic species of Allocapnia were found in Arkansas [12], several of which have recently been listed as species of greatest conservation need (SGCN) according to the Arkansas Game and Fish Commission's 2019 Arkansas Wildlife Action Plan (AWAP).
A. mohri is a regional endemic that was commonly collected in the last 1980's census [12], making it a good possible candidate for assessing important habitat factors and for distribution modeling. A. mohri populations during the previous census occurred broadly across Strahler stream orders one through five and were also found broadly across sites with differing flow permanence from intermittent sites with a completely dry stream bed or only pools with no flowing water for part of the year to perennial sites [12]. However, they were not detected at perennial sites with a significant groundwater flow [12]. This may be due to groundwater-based streams lacking the necessary temperature cue for egg hatching, whereas runoff dominated streams get cold enough with ice and snow melt entering the stream to cue hatching.
During the last widespread survey of Allocapnia in Arkansas, geospatial climate and land use data for stream catchments were not readily available [12] making it difficult to determine what factors are associated with A. mori distributions in the region. Species distribution models (SDMs) are empirical models relating field presence or abundance data to environmental predictor variables [40]. The use of SDMs has gained popularity in the past 20 years in ecological fields as geospatial mapping data with high resolution has become more readily available to conservation biology, making it possible to gauge the effect of climate, land use, and other environmental changes on the distribution of rare and understudied species [41].
One particular way of creating distribution maps that is gaining popularity is Random Forest (RF) modeling. Random Forest modeling is a group of tree predictors based on classification and regression tree (CART). The trees depend on the value of a random, independently sampled vector, so the method reduces bias [42]. Random Forest makes predictions by machine learning and then utilizes out-of-bag samples for model validation [43]. This method is excellent at minimizing overfitting due to its random nature and is effective at modeling spatial data even though it is not a spatial method. Random Forest modeling is superior to other popular SDM procedures, such as Mahalanobis Typicalities, a method adopted from remote sensing analyses. It is comparable to Maxent, a statistical mechanics approach, when one has a small-to-moderate number of the collected species presence/absence records (n = 38-94), and in species with low dispersal abilities [44]. The average number of records of each species of Allocapnia sampled in the previous census was 39 [12], and Plecoptera have reduced or absent flight capabilities, so Random Forest will likely work well for modeling winter stonefly species.
No SDMs exist for A. mohri to our knowledge. However, it was found to be the most widely distributed species of Allocapnia in Arkansas in the 1980's ( Figure 1A). Though this species used to be the most common through Arkansas, the species range has not been evaluated since the 1990s census. Also, due to the lack of geospatial data, no distribution models including land use and landcover or climate variables have been developed for A. mohri. Thus, our study objective was to recensus the Arkansas population and use the data to spatially model A. mohri distributions in Arkansas using RF models. Landscape, climate, and a combined landscape and climate data set were used to populate our RF model to determine which factors were most predictive of species presence. We then used multiple logistic regression to examine four hypotheses: higher elevations [45]. Hypothesis H2-A. mohri will be more common in areas with high percentages of forested land [46] because leaf litter inputs are an important food source for detritivorous stoneflies [47]. Hypothesis H3-A. mohri will be more common in streams with high rainfall and low water temperatures [46]. Hypothesis H4-A. mohri will be less common in streams with a slow soil infiltration rate, due to their sensitivity to stream impurities [12].
Our resulting RF SDMs are some of the first published for a North American Plecopteran species (but see [48]), and the landscape RF model performed better than the climate or combined RF models. The MLR modeling indicated both landscape and climate variables were important for predicting presence. A. mohri was absent at the Ouachita Mountains sites in 2020 and 2021 possibly indicating range changes [12].

Study Design
Our study area included 73 sample reaches stratified across the Level II Ecoregion Ozark/Ouachita-Appalachian Forests, which lies within the Level I Ecoregion, the Eastern Temperate Forests. The sample reaches were further spread into the following Level III Hypothesis H1-A. mohri will be more likely to occur in mountainous regions and higher elevations [45]. Hypothesis H2-A. mohri will be more common in areas with high percentages of forested land [46] because leaf litter inputs are an important food source for detritivorous stoneflies [47]. Hypothesis H3-A. mohri will be more common in streams with high rainfall and low water temperatures [46]. Hypothesis H4-A. mohri will be less common in streams with a slow soil infiltration rate, due to their sensitivity to stream impurities [12].
Our resulting RF SDMs are some of the first published for a North American Plecopteran species (but see [48]), and the landscape RF model performed better than the climate or combined RF models. The MLR modeling indicated both landscape and climate variables were important for predicting presence. A. mohri was absent at the Ouachita Mountains sites in 2020 and 2021 possibly indicating range changes [12].

Study Design
Our study area included 73 sample reaches stratified across the Level II Ecoregion Ozark/Ouachita-Appalachian Forests, which lies within the Level I Ecoregion, the Eastern Temperate Forests. The sample reaches were further spread into the following Level III Arkansas Ecoregions: Ozark Highlands (n = 15), Boston Mountains (n = 29), Arkansas Valley (n = 13), and Ouachita Mountains (n = 16) ( Figure 1). These Level III ecoregions historically had high stonefly diversity compared to other ecoregions in Arkansas [12]. Sample locations were based on historical studies reporting presence of Allocapnia genera [12]. We also stratified sample sites by Strahler stream order with several representing 1st order (n = 28), 2nd order (n = 15), 3rd order (n = 10), 4th order streams (n = 11), 5th order (n = 8), and 6th order (n = 1) (Table 1). Finally, the hydrologic flow regime was evaluated using the GeoCrawler Google Earth application (application located in Appendix B https://doi.org/10.1002/rra.2838, accessed on 1 January 2023). Due to the high amount of headwater streams sampled, and the stream locale belonging mainly to the Eastern Temperate Forests, most study sites had a high probability of belonging to the intermittent flashy flow regime (n = 46) [49].

Local Site Characterization
We measured the length and five wetted widths along evenly spaced transects for each habitat unit (e.g., riffles, pools, runs) at wadeable sites in August and September 2020-2022. Stream depth was measured at each width transect by taking five measurements across the cross-section of the stream. Each sample stream reach was 200 m in length, allowing all reaches to measure at least six times the wetted width of the stream, and each reach included at least one pool and riffle sequence [50]. Discharge was measured at the base of the reach using the mid-section method [51] at wadeable sites with no U.S. Geological Survey (USGS) gaging station. Otherwise, gaging station discharge was used.
Stream permanence was also estimated according to two separate metrics-one method described by Poulton and Stewart (1991) [12] and the other by Sheldon and Warren (2009) [52] in August and September of 2020-2022, which is generally the driest time of the year. The first method involved giving streams a ranking of A through D depending on the dryness over time, with an A ranking given to streams with a dry bed for part of the year and a D ranking given to streams with a permanent flow through the year along with a significant underground source. The second metric was a stream drying metric wherein we measured the linear extent of visible surface water in the stream reach in marked 50 m sections on each visit, then converted the wetted extent to a proportion, and calculated a time-weighted annual average for each site. This stream drying metric is primarily an index of time without surface flow, whereas the former gives a large-scale idea of the stream permanence. Using both metrics gives a clearer idea of the flow characteristics of each stream.

Landscape-Level Site Characterization
We collected two main types of landscape data: table and raster. Tabular climate and watershed data was collected virtually for all sites via Model My Watershed [53], and it was used in Multiple Logistic Regression (MLR) models (see Statistical analyses and modeling section). Several high-resolution landscape-level rasters were collected for the state of Arkansas to create distribution models of the species (see Statistical analyses and modeling section), including a land use raster (NLCD 2016) [54], a soil variable raster (gNATSGO) [55], and a large-scale environmental data raster (WorldClim) [56].

A. mohri Presence
Adult stonefly sampling started in early December and continued through mid-April during the 2020-2021 and 2021-2022 field seasons to match the emergence periods of A. mohri. Several approaches have been used to determine presence of Allocapnia populations in the past [12,52,57,58]. None of these methods have been evaluated to determine detection probability, which requires multiple samples over the emergence period. We chose to focus on adult collections, and all sample reaches were visited at least nine times across the emergence period to allow for evaluation of detection probability. We searched for adults in all areas of the 200 m reaches but particularly focused searches on compacted leaves and rocks in riffles, tree trunks, woody debris, and leaf litter debris piles on the side of the stream for 20 minutes per reach. All adults were collected by hand with forceps and immediately preserved in 95% ethanol, transported back to the laboratory, sexed, and identified to species using a regional taxonomic key [12]. The species were then placed on a distributional map.

Statistical Analyses and Modeling
We chose Random Forest (RF) to determine the possible distribution of a species of Allocapnia mohri relative to using MaxEnt and Mahalanobis because of its reduced bias with moderate to large sample size and transparent readout [43]. The RF species distribution models in this study were created in ArcGIS [59], which also produced out-of-bag score, mean squared error, Matthews Correlation Coefficient, and F1 statistics to evaluate model accuracy. The out-of-bag (OOB) error is the average error for each tree, calculated using predictions from the trees that do not contain it in their respective bootstrap sample. This score provides information on the expected performance on new, unseen data, and generally, smaller OOB scores correspond with better models. The mean squared error (MSE) is a measure of the prediction accuracy of the model. There are no set acceptable limits for MSE, but in general the lower the MSE, the more accurate the model. F1 scores are a weighted average of the precision and recall of a model, where an F1 score reaches its best value at 1 and worst score at 0. The Matthews Correlation Coefficient scores are based on results for all four confusion matrix categories (true positives, false negatives, true negatives, and false positives), and a higher score corresponds to a better result. It scores proportionally both to the size of positive and negative elements in the dataset, making it a great metric for unbalanced data [60,61].
Three different RF models were created to parse out effects of land cover, soil, and climatic data. Each model was trained with 80% of the collected presence and absence data, and 20% of the data was used to validate and test the models. Model 1 was the "Landscape Model," which included all land use and soil variable rasters along with an elevation raster. Model 2 was the "Climate Model," which included only climatic raster data. Model 3 was the "Combined Model," which included soil, land cover, and climatic variable rasters.
Random forest models are limited in that they do not depict the way that variables relate to one another, and they do not determine if the variables are significantly affecting A. mohri presence. Therefore, they are not appropriate for testing hypotheses about important habitat factors describing the landscape or climate niche of A. mohri in headwater streams. To resolve this issue, the data were evaluated using a multiple logistic regression (MLR). The MLR models were created in R (glm{stats}) [62]. A total of four separate models were created based on the four established hypotheses for this study ( Table 2). The Hypothesis 1 (H1) model evaluated the ecoregion and elevation variables in the study to see if A. mohri was more common in mountainous regions and higher elevations. The Hypothesis 2 (H2) model included land use data to determine if A. mohri was more common in areas with a high percent of forested land. The Hypothesis 3 (H3) model included annual precipitation total and annual mean temperature variables to determine if A. mohri is more common in streams with high rainfall and low water temperatures. The Hypothesis 4 (H4) model included water infiltration rate variables to determine if A. mohri would be less common in streams with impacted soils. All landscape and climate predictor variables were tested for model assumptions including binomial distribution, independence of error terms, linear relationship between logit(Y) and X, and lack of collinearity among predictors. After checking for assumptions, some variables were excluded from analysis due to collinearity issues (correlation among predictor variables), and the variables that lacked logistic restraints and were determined to be of highest biological importance in the model were retained. All models were then compared and ranked using Akaike Information Criterion (AIC) to determine which model was best at predicting A. mohri presence, and Hosmer-Lemeshow goodness-of-fit tests and likelihood-ratio tests were run to compare all models to the global model.

Results
We collected and identified 3463 Allocapnia mohri individuals, with 2519 males and 944 females. This species began emerging around the end of November with the first A. mohri collection being made on 23 November 2020, nearly coinciding with the previously stated early December start of the emergence period for A. mohri [12]. We found A. mohri adults at 18 of our 73 sample reaches. Adults were found across stream orders one through five with no bias towards a particular stream order (Pearson's chi-square = 3.42, p = 0.63, Figure 2A). A. mohri was not found at our 6th order stream site. A. mohri tended to be found in the Boston Mountains and Arkansas Valley ecoregions (Pearson's chi-square = 0.25, p = 0.04, Figure 2B). Adults were also found across all hydrologic permanence categories with no bias towards a particular category (Pearson's chi-square = 1.28, p = 0.73, Figure 2C).
Random Forest models estimated the potential distribution of Allocapnia mohri across Arkansas. The Landscape Model ( Figure 3A) had an out-of-bag score of 0.594, an F1 score of 0.86, and MSE of 6.938, and a MCC of 0.83. Further evaluation of the visual present/absent prediction output of the RF model shows that 11 known present sites fall within the produced predicted bounds and 43 known absent sites correctly fall within the produced predicted absent bounds, meaning the visual prediction output for the Landscape model correctly sorted 54 present and absent sites out of 73 (73.9%). Random Forest models estimated the potential distribution of Allocapnia mohri across Arkansas. The Landscape Model ( Figure 3A) had an out-of-bag score of 0.594, an F1 score of 0.86, and MSE of 6.938, and a MCC of 0.83. Further evaluation of the visual present/absent prediction output of the RF model shows that 11 known present sites fall within the The Climate Model ( Figure 3B) had an out-of-bag score of 0.375, a F1 score of 0.80, a MSE of 11.008, and a MCC of 0.75 (Table 3). Further evaluation of the visual present/absent prediction output of this RF model shows that 15 known present sites fall within the produced predicted bounds and 42 known absent sites fall inside the produced predicted absent bounds, meaning the visual prediction output for the Climate model correctly sorted 55 present and absent sites out of 73 (78.1%).  The Combined Model ( Figure 3C) had an out-of-bag score of 0.406, an F1 score of 0.57, a MSE of 9.631, and a MCC of 0.45 (Table 3). Further evaluation of the visual present/absent prediction output of this RF model shows that 15 known present sites fall within the produced predicted bounds and 38 known absent sites fall inside the produced predicted absent bounds, meaning the visual prediction output for the Landscape model correctly sorted 53 present and absent sites (72.6%).
Next, MLR models were used to test our study hypotheses. The H1 model showed total maximum elevation was significant in predicting A. mohri presence (p < 0.05, chisquare = 0.34), and elevation was positively correlated with species presence ( Figure 4A, Table 4). The H2 model showed that no land use variable was significant in predicting A. mohri presence (p > 0.05; Table 4); however, there was a negative correlation between urban land and species presence. The H3 model showed that the hottest month temperatures were significant in predicting A. mohri presence (p < 0.05, chi-square = 0.34), and higher temperatures were negatively correlated with species presence ( Figure 4B, Table 4). The H4 model was the best model according to AIC ranking. The results from the H4 model revealed that very slow soil infiltration was significant in predicting A. mohri presence (p The Climate Model ( Figure 3B) had an out-of-bag score of 0.375, a F1 score of 0.80, a MSE of 11.008, and a MCC of 0.75 (Table 3). Further evaluation of the visual present/absent prediction output of this RF model shows that 15 known present sites fall within the produced predicted bounds and 42 known absent sites fall inside the produced predicted absent bounds, meaning the visual prediction output for the Climate model correctly sorted 55 present and absent sites out of 73 (78.1%). The Combined Model ( Figure 3C) had an out-of-bag score of 0.406, an F1 score of 0.57, a MSE of 9.631, and a MCC of 0.45 (Table 3). Further evaluation of the visual present/absent prediction output of this RF model shows that 15 known present sites fall within the produced predicted bounds and 38 known absent sites fall inside the produced predicted absent bounds, meaning the visual prediction output for the Landscape model correctly sorted 53 present and absent sites (72.6%).
Next, MLR models were used to test our study hypotheses. The H1 model showed total maximum elevation was significant in predicting A. mohri presence (p < 0.05, chi-square = 0.34), and elevation was positively correlated with species presence ( Figure 4A, Table 4). The H2 model showed that no land use variable was significant in predicting A. mohri presence (p > 0.05; Table 4); however, there was a negative correlation between urban land and species presence. The H3 model showed that the hottest month temperatures were significant in predicting A. mohri presence (p < 0.05, chi-square = 0.34), and higher temperatures were negatively correlated with species presence ( Figure 4B, Table 4). The H4 model was the best model according to AIC ranking. The results from the H4 model revealed that very slow soil infiltration was significant in predicting A. mohri presence (p < 0.01, chi-square = 0.53), and the correlation between A. mohri presence and very slow soil infiltration was positive ( Figure 4C, Table 4). All MLR models were ranked using AIC, and the best MLR was found to be the H4 model (Table 4).
Hydrobiology 2023, 2, FOR PEER REVIEW 10 < 0.01, chi-square = 0.53), and the correlation between A. mohri presence and very slow soil infiltration was positive ( Figure 4C, Table 4). All MLR models were ranked using AIC, and the best MLR was found to be the H4 model (Table 4).

Discussion
Stoneflies are a major component to aquatic diversity [47]. Many species are globally imperiled [3,63], and yet little is known about species-specific habitat requirements and tolerances to changing land use and climate variables. Our study provides one of the first RF species distributional models for a North American plecopteran species. It also provides more detailed habitat and climate sensitivity information and updates distribution data that can aid conservation of a regional endemic winter stonefly species.
A 1980 s census provided some information about A. mohri's habitat [12] that is generally consistent with the present study's findings, but our study adds to their results. We also found A. mohri broadly inhabited stream orders and hydrologic permanence categories. Specifically, A. mohri has been historically found at intermittent lower-order sites that have a completely dry stream bed or only have pools during parts of the year; therefore, they may have desiccation-resistant, diapausing eggs [64]. They also occurred at higher-order perennial sites; we found them at perennial sites that also have a significant underground spring source. A. mohri was historically found in first through fifth-order streams, but not in sixth-order streams. We found this same pattern. However, we only censused one sixth-order stream, which was located in the Ouachita Mountains. So, our study has limited power evaluating species presence in larger rivers.
Curiously, we did not find A. mohri in any of our 16 Ouachita Mountains stream reaches. A. mohri was found in Ouachita Mountains streams in the 1980 s [12]. Further, it was found in the Little Glazypeau Creek and Alum Fork of the Saline River drainages in 2000-2001 [52]. We have study reaches in these watersheds and did not find them despite visiting the reaches nine times during the winter emergence period. Their absence in the Ouachita Mountains ecoregion may be a statistical artefact, but it is possible that their southern range is shifting northward. The Random Forest model using only landscape predictor variables was the best model of A. mohri distribution according to the OOB, MCC, F1, and MSE scores ( Table 3). The Landscape model scored the highest in MCC out of all models and this metric in particular is an appropriate metric for this dataset, being that it works proportionally to the size of both positive and negative elements in the data, and the dataset that was used was uneven, with much more absence values than presence values. The Landscape model had the best accuracy in correctly predicting A. mohri absence when evaluating the incorrectly sorted present and absent sites in the RF prediction output. However, the Landscape model incorrectly sorted 19 present and absent sites out of the total 74 sample sites and incorrectly predicted that A. mohri should be present in the Ouachita Mountains ecoregion.
The Climate model was determined to be less accurate than the Landscape model according to the OOB, MCC, F1, and MSE scores, though it was the most accurate at sorting the present sites in the RF prediction output and had only 16 incorrectly sorted present and absent sites out of 74 sample sites. Though the Climatic model had the lowest OOB out of all three models, those scores test how accurately the model predicts the data it was given, meaning our Climate model did a great job in predicting the sites that we already knew about. The purpose of our distributional models was to predict the distribution of A. mohri outward into the rest of the state, which the Climate RF model does not do as well as the Landscape RF model, according to the other metrics used. The Climate model did correctly predict that A. mohri would not be present in the Ouachita Mountains ecoregion.
A. mohri tended to be present in watersheds with a greater maximum elevation ( Figure 4, Table 3). Elevation is known to be an important component for plecopterans, as stonefly presence increases as elevation increases [12,45,65,66], and this is even true for rare stonefly species, like Beloneuria jamesae [67]. Elevation has also been found to be an important factor in stonefly abundance [45], showing that mountain streams are beneficial for plecopterans as a whole. The link between elevation and species presence could be due to high elevation streams being a lower stream size order (1-3) [68] with riparian vegetation contributing allochthonous inputs to the stream, which shredding macroinvertebrates utilize. The forested canopy also shades the stream and keeps water temperatures cooler compared to areas further downstream [69].
Stoneflies have been shown to be sensitive to stream temperature, and as water temperatures increase, streams become less hospitable to stoneflies [70]. Water temperature can be important to winter stonefly life cycles. For example, diapausing eggs may require a thermal cue to hatch [12]. Streams with average July temperatures reaching 26.4 o C did not tend to have A. mohri present ( Figure 4B). It is possible that these temperatures did not get low enough to cue eggs or nymphs [12]. Mean July and January temperatures increased as one moves from the Ozark Highlands south to the Ouachita Mountains ecoregion ( Figure 5A and B, respectively). The Ouachita Mountains streams had a higher minimum and maximum mean temperature across the year than other ecoregions possibly explaining the absence of A. mohri. Higher elevation streams in the region may be a refuge from warmer stream temperatures [71]. However, the maximum elevation of our stream study reach watersheds was also the lowest in the Ouachita Mountains ecoregion ( Figure 5B). Future studies should endeavor to sample some higher elevation streams in the Ouachita ecoregions to see if they might be A. mohri refuges. If these higher elevation streams act as shelter for temperature sensitive species, such as Allocapnia mohri, they might not be a haven for much longer. As climate change increases the temperatures of the vulnerable cool mountain streams, cold water adapted species like Allocapnia mohri will lose their thermal refuges [71] and eventually be replaced by more thermal generalist species, affecting the diversity of freshwater ecosystems [72][73][74][75].
peratures increase, streams become less hospitable to stoneflies [70]. Water temperature can be important to winter stonefly life cycles. For example, diapausing eggs may require a thermal cue to hatch [12]. Streams with average July temperatures reaching 26.4 o C did not tend to have A. mohri present ( Figure 4B). It is possible that these temperatures did not get low enough to cue eggs or nymphs [12]. Mean July and January temperatures increased as one moves from the Ozark Highlands south to the Ouachita Mountains ecoregion ( Figure 5A and B, respectively). The Ouachita Mountains streams had a higher minimum and maximum mean temperature across the year than other ecoregions possibly explaining the absence of A. mohri. Higher elevation streams in the region may be a refuge from warmer stream temperatures [71]. However, the maximum elevation of our stream study reach watersheds was also the lowest in the Ouachita Mountains ecoregion ( Figure  5B). Future studies should endeavor to sample some higher elevation streams in the Ouachita ecoregions to see if they might be A. mohri refuges. If these higher elevation streams act as shelter for temperature sensitive species, such as Allocapnia mohri, they might not be a haven for much longer. As climate change increases the temperatures of the vulnerable cool mountain streams, cold water adapted species like Allocapnia mohri will lose their thermal refuges [71] and eventually be replaced by more thermal generalist species, affecting the diversity of freshwater ecosystems [72][73][74][75].
Catchment soil infiltration category percentages were included in the best logistic model for predicting A. mohri presence ( Table 2). This result may be due to the inherent soil texture found in each ecoregion. Soil texture refers to the percentage of sand, silt, and clay, which are the major factors affecting water infiltration. Species presence dropped to zero in the Ouachita Mountains ecoregion, which had the highest infiltration values of all visited ecoregions ( Figure 5D). Soil in the Arkansas Valley, Boston Mountains, and Ozark Highlands ecoregions had a higher proportion of clay (a very small-pored soil) compared to the Ouachita Mountains, which has more sandy soil. Therefore, this correlation between very slow infiltration and species presence could just be due to the correlation between soil type and ecoregion.  Catchment soil infiltration category percentages were included in the best logistic model for predicting A. mohri presence ( Table 2). This result may be due to the inherent soil texture found in each ecoregion. Soil texture refers to the percentage of sand, silt, and clay, which are the major factors affecting water infiltration. Species presence dropped to zero in the Ouachita Mountains ecoregion, which had the highest infiltration values of all visited ecoregions ( Figure 5D). Soil in the Arkansas Valley, Boston Mountains, and Ozark Highlands ecoregions had a higher proportion of clay (a very small-pored soil) compared to the Ouachita Mountains, which has more sandy soil. Therefore, this correlation between very slow infiltration and species presence could just be due to the correlation between soil type and ecoregion.

Conclusions
Our resulting RF SDMs are some of the first published for a North American Plecopteran species. Though the landscape RF model performed better than the climate or combined RF models, our study indicates that climate and landscape factors are important habitat variables controlling A. mori distributions, as seen in our resulting MLR modeling, which suggested both landscape and climate variables were important for predicting presence. Even though landscape variables such as maximum elevation and soil infiltration categories led to strong models, these variables have likely not changed much over the last thirty years as A. mohri shifted from present to absent in censuses occurring in the Ouachita Mountains. The Landscape RF model also clearly predicts A. mohri should be present there, but the Climate RF model does not. Additional surveys for adults in the Ouachita Mountains stream reaches within predicted presence watersheds based on the RF model, particularly those with greater maximum elevation, could help determine if these reaches may be a refuge for A. mohri populations.
Climate change and associated warming stream waters can be a significant threat to macroinvertebrate taxa and can be responsible for species losses [3,65,70,76]. Warming stream temperatures may be causing losses of A. mohri in streams at the southern edge of its range (e.g., the Ouachita Mountains), shifting its distributions northwards. Evidence of range shifts due to effects of climate change has been observed in other macroinvertebrate species [75,77], and could be expected to occur for sensitive species in a continually warming climate.