Climatic Suitability and Distribution Overlap of Sawﬂies (Hymenoptera: Diprionidae) and Threatened Populations of Pinaceae

: Ecological Niche Models (ENM) are tools used to predict suitability, based on climatic variables selected and occurrence data of the target species, and characterize the environmental space. Sawﬂies (Hymenoptera: Diprionidae) are one of the main factors threatening forest health in Mexico, with cyclical population outbreaks and a wide range of hosts. In the present paper, we calculate the climatic niche in Mexico of three diprionids, Neodiprion abietis (Harris), N. omosus Smith, and Zadiprion rohweri (Middleton); the ﬁrst and the latter with recent records in Mexico, and N. omosus with presence in Mexico and Guatemala; contrasting them against the distribution records of host species in the country. The climatic suitability of N. abietis was high in the Sierra Madre Occidental where its hosts, Pinus ponderosa , P. strobiformis , and P. menziesii are distributed. For N. omosus , the environmental suitability was projected towards the Southern Altiplano, where it coincides with a small presence of its hosts P. leiophylla and P. ayacahuite ; however, it was possible to calculate its coincidence with more hosts in other biogeographic provinces. Pinaceae species considered under threat, Abies concolor , P. monophylla , and P. strobiformis , have populations within the environmental suitability of the sawﬂies.


Introduction
In Mexico, conifers (Pinophyta or Coniferophyta) are distributed within various types of vegetation [1], with two centers of diversity and evolution of pines: the Transverse Volcanic Axis (with extensions towards the Sierra Madre Occidental and the Sierra Madre del Sur) and the northeast of the country [2]. The country is considered a secondary center of Pinus L. diversity, with 42% of the species of the genus and high levels of endemic pine species (>50%) [3]. However, at present, various factors such as climate change, population growth, and deforestation due to land-use change threaten the distribution of Pinaceae species [4,5], and the survival of around 30 species from diverse genera, e.g., Abies, Picea, Pinus and Pseudotsuga, is at risk [6]. An additional problem is an infestation by phytophagous species, in particular, those that present periodic population outbreaks and that affect large areas [7].
In forestry, gregarious or solitary larval pine sawflies (Hymenoptera: Diprionidae) can cause problems of economic importance [8] and could be associated with outbreaks in their natural distribution sites, as well as recognized as invasive species [9][10][11].
Neodiprion abietis (Harris) is a species found in North America [8,12,13], with newly confirmed detection in Mexico [14]. Biological records indicate an annual generation, with no larval and adult activity from October to April in Mexico [14]. A strong local preference for certain hosts has been noted as a limitation of their distribution [15], in addition to the seasonal presence of the active states of this sawfly [12].
Neodiprion omosus Smith is a species widely distributed from Mexico to Guatemala (Smith, 1988; own collection data). It is considered to have a polyphagous preference for Pinus species [16]. On the other hand, Zadiprion rohweri (Middleton) is spread between the United States (Arizona, California, Colorado, Nevada, New Mexico, and Utah) and Coahuila, Mexico [17,18], affecting multiple species of pine trees.
In Mexico, climatic niche projections have been mainly constructed for species of regional distribution [1,19], but the study of species with the potential to expand their distributions to new areas or to flora with threatened survival, has not been addressed. For a non-native species, the pathway to reach new geographical areas includes several barriers that must be overcome (geography, captivity/cultivation, survival, reproduction, dispersal, and environment) [20,21]. Neodiprion abietis and Z. rohweri are species with recent records of presence in Mexico [14,18], while N. omosus is endemic to the country [16], which has been confirmed from field observation data. In this study, we estimate the ecological niches of these species as a means to explore the geographical spaces that meet the appropriate environmental conditions for the eventual development of their populations, as well as to estimate the host species that could be threatened.

Obtaining Presence Records
Records were obtained from field trips, entomological collections, scientific literature, and biodiversity databases with free access online. Material from field collection was preserved and mounted with standard practices. Sawflies were determined using original descriptions and keys and compared with digital images from BOLDSYSTEMS [22]. Voucher specimens were deposited in the Forest Insect Collection of the Pabellón Experimental Field, belonging to the Instituto Nacional de Investigaciones Forestales, Agrícolas y Pecuarias (Aguascalientes, Mexico) (Supplementary Material, Table S1). Databases were assembled as follows. For each species, to reduce bias resulting from spatial clustering of occurrences and to improve the performance of the models [23], the 'spThin' package [24] was used with Rstudio ® ver. 3.3.
[25] and points with less than 5 km between each other were eliminated. All occurrence records from biodiversity databases were checked for credibility, taxonomy, and nomenclature [26], then filtered removing: (1) records with low coordinate precision, (2) duplicate data, (3) records without geographic coordinates, and (4) those outside its historical range [27,28].
Neodiprion abietis. The total number of georeferenced occurrence records was 67: 48 from open access biodiversity databases (GBIF, Naturalist) [29,30], 16 records from scientific literature, and 3 field collection records (Supplementary Material, Table S2). Fifty-one occurrence records passed the spatial correlation filter. Subsequently, the occurrence records were divided into training points (to perform the model) and evaluation points. Given that most of the occurrence records were gathered from databases, 37 occurrence records (72.54%) were used as training points, and 14 (27.45%) as evaluation points (Supplementary Material, Table S2).
Neodiprion omosus. The total number of occurrence records was 51, 47 of which were derived from field collections and four from scientific literature (Supplementary Material, Table S2). After spatial correlation was performed, 27 were kept for further analysis. Since most of the occurrence records were derived from field collection data, 30% (8) were separated for evaluation points and 70% (19) for training points, selecting them at random.
Zadiprion rohweri. The database had 22 occurrence records in total, from which one occurrence record corresponded to field collection, and 21 were obtained from scientific literature (Supplementary Material, Table S2). From this, 17 occurrence records passed the filter, with 30% (5) separated for evaluation points and 70% (12) for training points; these were all selected at random as most came from scientific literature collections.
The distribution of the host plants was obtained from the GBIF and Southeast Regional Network of Expertise and Collections [31] databases. Occurrence records were cleaned and filtered as mentioned above. The list of selected hosts is based on [12,16,18,[32][33][34][35] and included in Table 1.

Environmental Variables
Nineteen bioclimatic variables were downloaded from www.worldclim.org (accessed on 31 May 2022) with a spatial resolution of 30 s; the area of projection is Mexico (1:4,000,000). Two analyses of ecological niche models (ENM) were run to create candidate models for each species. The first aimed to obtain the jackknife values to know the contribution of each bioclimatic variable and then remove all correlated variables (0.80, −0.80) based on Spearman correlation performed in Past ® ver. 4.07 [36]. Then, in the second analysis, to statistically evaluate the highest possible combination of environmental and Maxent parameterization variables, the candidate models were created from two folders (set 01 and set 02). For both sets, bioclimatic variables without correlation were included. In set 01, only those that contribute to the calculation of the model were included, while for set 02 all the variables were included, regardless of their contribution (see Tables 2 and 3, and  Supplementary Material, Table S2). Table 2. Parameter values used to define ecological niche models of three species of sawflies (Hymenoptera: Diprionidae) in Mexico. M = Model; multiplier regularization values = 1, 2, 3; F = features classes -linear (l), quadratic (q), product (p), threshold (t), hinge (h).

Model Calibration
Ecological Niche Models (ENM) applications always involve the transfer of models to areas of interest. Environmental data from one region to another may include new environmental conditions [39], so calibrating the models is critical. The species calibration area (M) representing areas that have been accessible to the species during its biogeographic history [40,41] was defined by world ecoregions [42] characterized by their biodiversity: endemism of species, the rarity of high rates, species richness, unusual ecological or evolutionary phenomena, and the global rarity of their habitat type (e.g., Mediterranean forests, temperate and scrub forests). The selection of ecoregions was subject to the presence of the species, as well as the suitability of forest areas, according to the BAM (favorable biotic conditions, favorable abiotic conditions, and accessibility of the species to this set of favorable biotic and abiotic conditions) diagram for understanding species' likely distributions [43,44]. Environmental information was extracted from the training points in QGIS ® ver. 3.16 Hannover [45].

Development of Candidate Models
The candidate models were run with the Rstudio ® ver. 3.3 interface, Maxent ver. 3.4 [46] and the Kuenm package, which automates important calibration and evaluation steps in ENM [47], following the recommendation to evaluate the best potential combination of parameters (features classes and regularization multiplier) to select the most appropriate model [48]. The kuenm_cal function of the Kuenm package [47] allows for the creation of different models with different parameters. In this study, the multiplier regularization values were 1, 2, and 3. The levels of this parameter reduce or increase the over adjustment of the models. Regularization multiplier levels 2 to 4 times higher than 1 produce a substantially lower over adjustment, low multiplier regularization values tend to result in overly complex models [49]. In turn, they were combined with all possible compositions of feature classes: linear (L), quadratic (Q), product (P), threshold (T), and hinge (H) with a total of 174 candidate models for each set per species.

Model Evaluation
The evaluation was performed in Rstudio ® ver. 3.3 with the kuenm_feval function of the Kuenm package [49]. We evaluated 174 models for each species distributed in two sets. This package evaluates all candidate models and selects one or more models. Model selection is based on importance, predictability, and complexity, in that order of priority. Models are first filtered to detect those that are statistically significant; the omission rate criterion applies to this small set of models; finally, among the significant and lowemission candidate models, those with AICc delta values of less than two [49], lower omission rate, and significance value (pROC) ≤ 0.05 -through the kuenm_feval function, evaluating the relationship between the omission error for independent points and the

Model Extrapolation Analysis
The Multivariate Environmental Similarity Surface (MESS) analysis was performed with the Maxent software ® ver. 3.4.1. (Steven J. Phillips, New York, NY, USA) [51]. This analysis compares the environmental similarity of the native area from which M (accessible area for the species) was delimited and identifies areas where one or more environmental variables are outside the training range. Negative values indicate a new climate, and magnitude indicates the degree to which a point is outside the range of its predictors [52].

Parameterization and Environmental Variables That Defined the Models
The results of model evaluation for each species are shown in Table 2. The ecological niche was defined by a combination of bioclimatic variables particular to each species of sawfly (Table 3). In general, there were few contributing bioclimatic variables (3)(4) necessary to calculate the ecological niche, and for each model, the great majority of the contribution came from only two variables (Table 3).

Environmental Suitability
Environmental suitability in Mexico varied among the species studied. For N. abietis the ecological niche was estimated mainly along the Sierra Madre Occidental, until reaching the limit with the Volcanic Axis, where the environmental suitability of the central zone to the east of the latter biogeographic province decreases. From this point, the calculated niche extends north along the Sierra Madre Oriental and south along the southern Madre Sierra del Sur. Medium to almost high environmental suitability was present in the external limits of the provinces Soconusco and Altos de Chiapas. Medium environmental suitability was predicted along the northern Gulf of Mexico ( Figure 1A).
Neodiprion omosus showed lower environmental suitability compared to N. abietis. Furthermore, N. omosus exhibited disjointed suitability. The main area that met its environmental requirements was located in the Southern Altiplano, mainly towards the Sierra Madre Occidental and north of the Volcanic Axis. Medium to higher suitability was extrapolated in areas that run irregularly along with the Sierra Madre Oriental, south-west of the Volcanic Axis and center-south of Oaxaca, with smaller areas with medium to lower suitability calculated south of the Gulf of Mexico and south-west of Petén. In the province of Baja California, it was possible to project good environmental suitability towards the west ( Figure 1B).
Of the three species studied, Z. rohweri showed the most area with environmental suitability focused on two large spaces. The first of these extended across the province of California, the central and northern part of Baja California as well as with a tendency to decrease environmental suitability in the province of Del Cabo. A larger area was calculated in the center of the country, extending from the Pacific Coast to the Sierra Madre Oriental, and from northwest of the Volcanic Axis to the south of the Northern Altiplano. Smaller regions with good environmental suitability were estimated east of the Volcanic Axis and southwest of Oaxaca, and with lower suitability in central and northern Veracruz, and northwest of Petén ( Figure 1C).

Distribution of Host Species
The environmental suitability model calculated for N. abietis coincides very well with the distribution of its host plants registered in Mexico, mainly with P. ponderosa, P. strobiformis, and P. menziesii (Figure 2A,B). The distribution of Abies concolor is restricted to small populations in North Mexico [54] and is not illustrated in this manuscript, due to protection issues [6,54], and to minimize the risk of species extinction [55]; however, it coincided with low environmental suitability for N. abietis. The records of P. strobiformis and P. menziesii located more to the south coincide equally with the environmental suitability of N. abietis in the Sierra Madre del Sur and around Los Altos de Chiapas.

Multivariate Environmental Similarity Surface (MESS) Analysis
For N. abietis, a bioclimatic analogy between its native distribution and the model was projected in the biogeographic provinces of California, the northwestern Northern Highlands,  The projection of environmental suitability of N. omosus in the biogeographic province of the Altiplano Sur gathered the largest calculated area, but did not correspond well to records of host plants (P. leiophylla and P. ayacahuite). However, for the Sierra Madre Occidental, the host plant distribution suggests the possibility of a wider distribution of the sawfly than previously recorded. For the Volcanic Axis, there is a coincidence between the calculated niche and the distributions of P. leiophylla, P. ayacahuite and P. lawsonii. For the Sierra Madre Oriental, the calculated niche coincides with populations of P. patula, a species recorded from the north of this province south to the biogeographic province of Oaxaca, and in adjacent parts of the Sierra Madre del Sur. In these two provinces, the ecological niche of N. omosus coincides with at least three hosts: P. ayacahuite, P. leiophylla and P. patula. On the contrary, although environmental suitability was calculated for the biogeographic provinces present in the peninsulas of Baja California (California, Baja California, Del Cabo) and Yucatan (Yucatan, Petén) and east of the Gulf of Mexico, there are no records of the species here included as hosts ( Figure 2C-E).
Zadiprion rohweri showed environmental suitability that coincides with the total distribution of P. monophylla populations (not illustrated due to protection issues) [6,50]. It also coincides with the distribution of P. cembroides in the part between the Southern Altiplano, the Sierra Madre Occidental, and the Sierra Madre Oriental. The populations of P. cembroides present in Del Cabo escape the calculated environmental suitability for Z. rohweri, as do the populations present in the border areas of Mexico with the United States. All records of P. edulis in the country (Southern Highlands, Volcanic Axis, Sierra Madre Occidental) match areas of high suitability for Z. rohweri ( Figure 2F).

Multivariate Environmental Similarity Surface (MESS) Analysis
For N. abietis, a bioclimatic analogy between its native distribution and the model was projected in the biogeographic provinces of California, the northwestern Northern Highlands, and along the Sierra Madre Occidental); therefore, the predicted reliability is more appropriate because they do not present new combinations of environmental variables for the species. The MESS analysis confirms that the risk of extrapolation is low in the areas predicted for the calculated niche model of N. abietis, since the non-analogous areas were projected in the Gulf of Mexico and Pacific Coast biogeographic provinces ( Figure 3A).
For N. omosus, the situation was different. The MESS analysis shows several areas where this sawfly will face novel climates (northern areas of the country), as indicated by its negative values (until −63.118), and represent the strict extrapolation of the analysis. The response of the species to the conditions present in the biogeographic province of Oaxaca and the Sierra Madre Oriental should also be analyzed as they are less analogous to the ecological niche of the species ( Figure 3B).
In the case of Z. rohweri, MESS did not reveal areas of coincidence with the projected ecological niche model, which shows that the risk of extrapolation is low ( Figure 3C). calculated niche model of N. abietis, since the non-analogous areas were projected in the Gulf of Mexico and Pacific Coast biogeographic provinces ( Figure 3A).

Discussion
In ENM, all the "inferences made from the models are only as good as the models themselves" [56]. Model selection is then a careful step and must be taken following rigorous criteria [46]. Several metrics for evaluating the model performance have been proposed [57], but some of them are inadequate (i.e., the area under the curve of the receiver operating characteristic [58]). Researchers in this area must select models with statistical significance, high predictive ability, and low model complexity [35,40,59]. Here, the ecological niche models were characterized by bioclimatic variables selected by rigorous statistical criteria. Although there is no consensus on the optimal number of variables that define an ecological niche, those calculated here meet the condition of parsimony: less complex models can give the most logical explanation of the environmental suitability of the species of interest, and can estimate its potential distribution (but not the distribution per se) [60]. However, the low number of variables in the definition of the models (three for N. abietis: precipitation of wettest month, precipitation seasonality, and mean temperature of the coldest quarter; three for N. omosus: isothermality, precipitation seasonality, minimum temperature of coldest month, and temperature annual range; and four for Z. rohweri: mean temperature of wettest quarter, temperature seasonality, and precipitation of driest month) can be considered as a sign of climate vulnerability, since dependence on one or two variables could make the species vulnerable to sudden changes over time [61].
The ecological niche of N. abietis is delimited by two bioclimatic variables associated with precipitation: precipitation of the wettest month and precipitation seasonality, which together explain 99.8% of the model, similar to that recorded for the white coffee stem borer Monochamus leiconotus (P.) [62]. Precipitation of the wettest month has been identified in previous studies as an important determinant in the distributions of many species and is associated with environmental suitability for phytophagous insects in particular [63,64]. Precipitation seasonality is the variation in annual monthly precipitation [38], with a deep effect on water quality and availability [65]. This bioclimatic variable has a marked impact on other insect species (i.e., Tryoza erytreae Del Guercio, Anthophora curta Provancher, A. squammulosa Dours) [66,67], indicating its association to environments with contrasting wet and dry seasons. The current distribution of this diprionid indicates there is a disjointed distribution, possibly influenced by sampling bias, with an abundance of records in Canada and the United States [8,12,13], but few in Mexico [14]. However, according to the estimated climatic niche, it may be more widely distributed, and future collection efforts could confirm this prediction.
The ecological niche of N. omosus was determined mainly by temperature isothermality, and precipitation seasonality, which together amounted to 94.1% of the information (Table 3). Isothermality quantifies the range of temperature oscillation from day to night with respect to the annual oscillation from summer to winter [38], indicating low tolerance of this sawfly to sites with large temperature fluctuations between daylight and nighttime hours [68].
In contrast, for Z. rohweri, the most important bioclimatic variables are related to temperature (mean temperature of the wettest quarter and temperature seasonality), which together make up 88.1% of the calculated niche model ( Table 3). As in the case of the studied species of Neodiprion, it seems to be true that each species of Zadiprion has a particular set of bioclimatic variables to define the ecological niche [i.e., annual mean temperature was the best predictor for Z. jeffreyi Smith [10]). It is necessary to carry out a greater number of studies to understand the influence of these variables on the biological cycle of sawflies. As with N. abietis, researchers should look for this species in the areas identified as suitable.
In general, when studying the environmental suitability of insect pest and their hosts, there is a tendency for an adequate match between both niches [69]. In the case of the sawflies (Hymenoptera: Diprionidae) here studied, it seems that there is a correspondence between each modeled niche and the observed distribution of the hosts. These diprionids are generally polyphagous species with wide ranges, where they are often exposed to different host species [16,18], these factors have an important influence on fewer defoliation rates [70]. However, the discrepancies found between the ecological niche of N. omosus and the complete distribution of its host plants indicate that their projections should be taken with caution. Additional future distribution occurrences and host association information will reveal a better estimate of its ecological niche.
Sawflies are a problem for forest stands when sudden outbreaks occur, associated with climatic factors unfavorable to plant species, due to conditions of the trees themselves [70], and/or by physical aspects of the affected sites [71]. For Mexico, the calculated ecological niche of three species of sawflies allows us to estimate that there are new areas where the combination of environmental conditions and the presence of host plants might favor the development of population outbreaks that could be harmful to populations of Pinaceae species. It is recommended to establish preventive monitoring schemes for the early detection of these diprionid species. Particular attention must be paid to A. concolor, P. monophylla, and P. strobiformis, threatened species [6,72] whose protection also benefits associated flora and fauna [73].
Finally, the environmental variables used in this study represent current conditions and can only be used for short-term predictions and management decisions, or for conducting intensive sampling efforts.

Conclusions
The ecological niches of three sawflies (Hymenoptera: Diprionidae), Neodriprion abietis, N. omosus and Zadiprion rohweri in Mexico were calculated with a low number of bioclimatic variables (3-4) but strong statistical support. The combination of modeled ecological niches and observed host distributions reveals that dioprionid sawflies have the potential to increase their distribution in places where they are not recorded yet. Endangered host plants either have low environmental suitability (Abies concolor for N. abietis), allowing them to escape from potential damage, or occur in places with the ideal conditions for a sawfly geographical expansion (P. strobiformis for N. abietis; P. monophylla for Z. rohweri), thus increasing their vulnerability.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/f13071067/s1, Table S1: Specimens deposited of Neodiprion abietis and Neodiprion omosus in the Forest Insect Collection of the Pabellon Experimental Field, belonging to the Instituto Nacional de Investigaciones Forestales, Agricolas y Pecuarias (Aguascalientes, Mexico), and used in the present study. Table S2: Distribution points used for the ecological niche calculation of three sawflies with forestry importance in Mexico; Table S3: Bioclimatic variables, regularization multiplier levels, features classes used for each sawfly species for the calculation of ecological niches.

Data Availability Statement:
The data presented in this study are available in Supplementary material in Tables S1-S3.

Conflicts of Interest:
The authors declare no conflict of interest.