Modeling potential distribution of newly recorded ant, Brachyponera nigrita using Maxent under climate change in Pothwar region, Pakistan

Climate change has been discussed as to exert shifts in geographical range of plants, animals or insect species by increasing, reducing or shifting its appropriate climatic habitat. Globally, Pakistan has been ranked at 5th position on the list of countries most vulnerable to climate change in 2020. Climate change has resulted in the losses of biodiversity and alteration in ecosystem as a result of depletion of natural habitats of species in Pakistan as well as in the world. Ants have been regarded as indicators of environmental change and ecosystem processes. Brachyponera nigrita (Emery, 1895) was reported for the first time from Pakistan (Pothwar region). Objective of our studies was to model geographic distribution of newly recorded ant species, B. nigrita based on two representative concentration pathways (RCP) (RCP4.5 and RCP8.5) for 2050s using maximum entropy model (Maxent) in Pakistan. In modeling procedure, 21occurrence records and 8 variables namely Bio4 (Temperature seasonality), Bio8 (Mean temperature of wettest quarter), Bio10 (Mean temperature of warmest quarter), Bio12 (Annual precipitation), Bio13 (Precipitation of wettest month), Bio15 (Precipitation seasonality), Bio17 (Precipitation of driest quarter) and Bio18 (Precipitation of warmest quarter) were used to determine the current and future distributions. Performance of the model was evaluated using AUC (area under curves) values, partial ROC, omission rates (E = 5%) and AICc (Model complexity).The results showed the average AUC value of the model was 0.930, which indicated that the accuracy of the model was excellent. The jackknife test also showed that Bio4, Bio18, Bio17 and Bio15 contributed 98% for the prediction of potential distribution of the species as compared to all other variables. Maxent results indicated that distribution area of B. nigrita under future predicted bioclimatics 2050 (RCP 4.5 and RCP8.5) would be increased in various localities of Pakistan as compared to its current distribution. In Pothwar region, moderately suitable and highly suitable areas of this species would increase by 505.932321km2and 572.118421km2as compared to current distribution under 2050 (RCP 4.5), while under 2050 (RCP 8.5), there would be an increase of 6427.2576km2and 3765.140493km2 respectively in moderately suitable and highly suitable areas of B. nigrita. This species was associated with termites, collembolans and larval stages of different insects. White eggs, creamy white pupae and many workers of this species were observed in a variety of habitats. Unknown nesting ecology, species identification characters supported with micrographs has been given which will help researchers for further ecological studies.


Introduction
Climate Change has been recognized as obvious and emerging global phenomenon. It's adverse effects on biological entities have also been reported throughout the world [1][2][3][4][5][6]. The Fifth Assessment Report recently presented by IPCC (Intergovernmental Panel on Climate Change) emphasized on the reduction of anthropogenic climate change, produced by greenhouse gas emission [7]. The main reason for undertaking such types of protocol is to protect biodiversity from the impacts of global warming [8][9][10]. Impacts of climate change has been reported for humans as well as biodiversity, including disease [11,12] and pest outbreaks [13], shift in species distribution, occurrence and Phenology [14]. The condition is predicted to become spectacular as the magnitude of environmental change accelerates [15]. This situation is expected to be more adverse in developing nations [16]. Climate change has direct effect on the survival, distribution, development [17,18] and abundance of living organisms [19][20][21]. Arthropods are ideal organisms for determining biodiversity. Among insects, members of family Formicidae are significant target group to examine biodiversity [22]. Ants are abundant, ubiquitous and also recognized to be very sensitive to change in temperature and humidity as they affect survival and foraging activity [23].
The taxonomic status of Brachyponera Emery has complex history, is presently synonymized as Pachycondyla, very large genus [24,25]. It comprises species from Oriental, Australian and Ethiopian regions. Oriental region includes a few species from tropics as limited taxonomic work has been done there. Genus Brachyponera is reasonable genus with 19 valid species and 5 subspecies, remarkable among other Ponerines with distinct size and dimorphism among workers and queens. Some species exhibit invasive behavior [26]. Workers of genus Brachyponera are small, solitary, scavengers and predatory in behavior. Most of the species construct their nests in rotten wood or soil [27]. This genus is distinct from other Ponerinae species due to the presence of pests and invasive species having dangerous sting.
ENM (Ecological niche modeling), similar mechanisms of habitat suitability and species distribution modeling have gained recognition as they enable estimation of possible geographic species distribution under climate change scenarios [28]. The methods provide information about the environmental variables related to the ecological niche inhibition on species related to future distribution approximated by climate models [29]. These models result in coarse resolution correlative models of ecological niche that predict the potential geographical distribution of species under current and future climatic conditions. Maxent has been used in a number of studies involving prediction of different ant species in different parts of the world. Senula et al. [30] applied this model to predict the habitat suitability of ants belonging to the genus Trachymyrmex and Mycetomoellerius in the United States. Song et al. [31] modeled Solenopsis invicta (Hymenoptera: Formicidae) using maxent for current and future projected climatic conditions in China. Koch et al. [32] implemented maxent and determined areas of highest suitability for Gracilidris pombero Wild & Cuezzo, 2006.
A limited work has been done on the taxonomy of ants in Pakistan. Recently, Rasheed et al. [33]. recorded 103 ant species from Pakistan. Prior to the present studies, no species of genus Brachyponera has been reported from Pakistan. To the best of our knowledge, to date no species of insects particularly ants have been modeled for projected climatic conditions in Pakistan. Changing climatic conditions would have numerous impacts on the ecological interactions of the species specifically which are affected by environmental factors. Keeping in view this situation, newly recorded ant species, B. nigrita was modeled for current and future climatic conditions in Pakistan.

Ant's collection and identification
Ants were collected from different localities and diverse habitats like leaf litter, under stone, tree barks, soil surface, crop plants, grasses, ornamentals etc. Nests of ants were photographed with the help of professional camera. Ants were collected by using aspirator, net sweeping, handpicking etc and killed using killing bottle (Potassium Cyanide). The killed specimens were preserved in 75% alcohol. Preserved specimens were identified using stage microscope and identification keys by Bingham [34] and Bharti [35]. Micrographs of identified specimens were prepared using tri-nocular microscope (Leica MS5 microscope) attached with camera (Amscope 18 megapixel).

Survey and data collection
About thirty multisite surveys were conducted in four districts of Pothwar region including Capital territory (Islamabad) for the collection of ants in various habitats ranging from grasslands, mountains, forest areas etc. Every sampling site was surveyed 500 ×500 metrs for collection of ants from trees, flowers, ant nests, plants etc. As a result of these surveys, B. nigrita was collected from 29 sites. Data was recorded in the form of geographical coordinates of B. nigrita occurrences with the help of GPS device (Garmin e Trix 10). Data collection technique used in present study was the presence-only data collection method as used in various studies for species distribution modeling [36][37][38][39]. The R package (The spThin) was used for spatial thinning of occurrence points [40]. In order to reduce biases related to spatial autocorrelation, occurrence data were rarefied in a pattern that a pair of points was not closer than three kilometers. This process resulted in a total of 21 occurrence points of B. nigrita. Out of these 17 occurrences were processed in ArcGIS software (Version 10.4.1) and a buffer 70 km was made around all occurrence points for setting model calibration area (M) (Fig 1). M (Calibration area) area is an area which is considered an accessible area for a species over time as described by Barve et al. [41]. Remaining 4 occurrences were used as independent data.

Bio-climatic predictors
A set of 19 bio-climatic predictor variables were obtained from official source (WorldClim website; Global Climate Data) for B. nigrita species distribution and suitability analysis. They were raster layers with a spatial resolution of 30 arc seconds. These variables are often used in many ecological and biogeographic studies for modeling species distribution, niche modeling and habitat suitability of different organisms [42,43]. They are availed from various sources like WorldClim database, CliMond database etc.
These bioclimatic variables have been approved to indicate general patterns of precipitation and temperature, extremity and seasonality of temperature. These layers were converted according to model calibration area (M) by extract by mask function using ArcGIS software. These were then converted to Ascii format. The multicollinearity between 19 bio-climatic predictors was determined using Pearson correlation coefficient (r) with help of ENM Tools software v1.3.1.,cross-correlation value set at (r) � 0.90 in order to remove strongly correlated predictors [44][45][46][47]. (Table 1). The 8 variables were selected based on their potential relevance to explain the distribution of species. The Maxent (Version 3.3.4k) [48] was run along with occurrence points, the results of jackknife test were used to determine dominant environmental variables that determined the species potential distribution. After screening of climatic variables, finally eight variables selected for the model run; Bio4 (Temperature seasonality), Bio8 (Mean temperature of wettest quarter), Bio10 (Mean temperature of warmest quarter), Bio12 (Annual precipitation), Bio13 (Precipitation of wettest month), Bio15 (Precipitation seasonality), Bio17 (Precipitation of driest quarter) and Bio18 (Precipitation of warmest quarter).

Settings of the Maxent model using R package Kuenm
Maxent was run using R package Kuenm [49] with the aim to produce preliminary models with occurrence data and environmental predictors. Maxent has been used for modeling different species in various parts of the world and is a renowned tool implemented for species distribution modeling. R package Kuenm automates ENM processes; help the users to manage model complexity and testing of data with a number of parameters. Occurrence data (17 points) was divided randomly into 5 testing and 12 training for model calibration and internal testing. We already set aside 4 occurrences as independent data for independent model testing. Three sets of selected 8 variables were also randomly prepared as Set1 (Bio4, Bio10, Bio13 and Bio18), Set 2 (Bio8, Bio12, Bio15 and Bio17) and Set 3 (combination of all eight selected variables).These three sets were named as M variables. Kuenm_ceval function of Kuenm was used for the preparation of preliminary models. Calibration of these models was done for the evaluation of best possible combination of parameters to be selected for running Maxent in order to prepare best model for our species. Candidate models were prepared using all combination of 17 value of regularization\ multiplier and 29 feature class combinations as used by Cobos et al. [49]. Later on evaluation of these models was carried out on the basis of partial ROC (with 500 iterations), omission rates (E = 5%) and AICc (Model complexity).Then from these model sets, best model was selected with delta AICc values � 2. kuenm_ceval function was used for evaluation and selection of best candidate model. Similarly kuenm_mod and kuenm_feval function of the package were used for evaluation and selection of final models. Maxent was run with 500 iterations, 10 random replicate analyses, all other options were of default settings of kuenm package.

Transfer of calibration area model results
After the selection of best model for calibration area, model results were transferred to Pothwar and Pakistan region for current and future scenarios with RCPs 4.5 and 8.5 for the year of 2050 based on 5th assessment report of IPCC. For this purpose parallel bioclimatic data layers were obtained for three GCM (HadGEM2-AO, HadGEM2-ES and BCC-CSM1-1). These predictor layers were prepared with same procedure as given above as G variables and model results were projected to two regions and two future scenarios. Finally we explored 3 GCMs 1 time period 2 RCPs = 6 future transfers of models. Results of Maxent were imported in ArcGIS software and further processed for the preparation of maps. Raster layers (Pothwar region only) were projected to WGS 1984 UTM zone 42 North and reclassified for the measurement of geographical area of the species. Area calculation was performed only for Pothwar region (area of interest).

Development of binary models
Binary models were developed on the basis of relationship of calibration points and raw model output of maxent. Raw outputs were imported to ArcGIS software along with calibration points. Then with extract values to point's function of ArcGIS, we decided threshold values; the highest thresholds corresponding to 0%, 5%, and 10% omission error rates as used by Ashraf et al. [50]. Change in high, moderate, and low suitability areas was calculated with same criteria for current and all future scenarios. The set threshold were as follows threshold values from0-0.18586 (E = 0-5%) as low suitability area, 0.18586-0.22912 threshold value (E = 5-10%) as moderate suitability area and above 0.22912 threshold value as highly suitable area. Future-transfer maps were also prepared with the same criteria as done for present day.

Identification characters of Brachyponera nigrita (Emery, 1895) (Fig 2)
Description (worker) head. Head in full-face view clearly longer than wide with posterior margin straight and feebly convex sides; mandibles with series of small-median tooth, apical one acute and longest; clypeus narrow and slightly emarginated posteriorly; eyes of small size, maximum diameter (0.16 mm); head in full-face view eyes located below the mid length of head and just parallel to antennal insertion; antennal scape hardly reaching posterior head margin. Mesosoma. Promesonotal outline slightly indented; mesonotal grove distinct; propodeal spiracles rounded-ovoid; propodeal dorsum in profile making well defined declivity. Sculpture. Cephalic surface with minute punctations; clypeus with faint vestiges; mandibles with minute sculpture throughout; dorsum of pro-mesonotum smooth, metanotum weakly punctured; petiole slightly reticulate; gaster tergites smooth and shining but weakly reticulate in profile. Pilosity. Dorsum of head with small, white pubescent allover; antennae with abundant setae on funicular segment; clypeus with 4-5 long protrusive setae; gaster with scattered erect and suberect setae in profile. Color. Body uniformly dark brown-blackish. Nesting ecology (Fig 3) B. nigrita was collected from Northern and forest areas of Himalayan foothills of Pakistan, in Pothwar region. This species was collected from variety of habitats. Mostly its nests were found under the stones or under piles of decomposed leaves especially of pine leaves. They were found associated with collembolans in grass roots under the stone. One colony of these ants was observed under the stone with eggs, creamy white pupae along with isopods. Multiple nests under stone were also observed having plant debris in the nest. Termites were also present in some nests of these ants. Lepidopterous larvae have also been observed in a few nests. Many nests were present along the tree trunks under the bark of tree attacked with termite. Some nests were also found around the aphid colonies along the trunks of some plants. The collection sites were characterized by good amount of moisture, relative low temperature and dense pine forest. environmental variables. Model performance was evaluated based on statistical significance (Partial_ROC), omission rates (OR) and the Akaike information criterion corrected for small sample sizes (AICc) [52] (Fig 4A and 4B). Among these only 1251 models were statistically significant and 2 were statistically significant models meeting omission rate and AICc criteria (Tables 2 and 3 and Fig 5). We used average of these two models for final model map preparation. Best model

PLOS ONE
performance was determined at 5% threshold (Omission rates). Map of calibration area (Fig 6) predicted with maxent indicate that this species is currently distributed in Kashmir, Punjab, coastal areas of Sindh and Balochistan and a limited area of KPK provinces of Pakistan.
Percentage contribution of predictor variables. The results of Jackknife test after running Maxent gave relative contribution of each bioclimatic variable for B. nigrita prediction (Fig 7). The main variables for predicting potential distribution of B. nigrita were Bio4 (Temperature seasonality), Bio18 (Precipitation of warmest quarter), Bio17 (Precipitation of driest quarter), Bio15 (Precipitation seasonality), Bio8 (Mean temperature of wettest quarter) and Bio12 (Annual precipitation) for the current geographical distribution with their contributions of 70.6%, 11.7%, 10.7%, 5%, 1.3% and 0.8% respectively (Table 4). Bio12 is the least contributing variable out of these variables. As a whole these four variables i-e (Bio4, Bio18, Bio17 and Bio15) contributed 98% for the prediction of potential distribution of the species as compared to all other variables. Bio4, Bio18 and Bio17 are also important for permutation (Table 4).
Model performance evaluation based on AUC values. When maximum entropy principle of Maxent is applied, it has been concluded that areas under receiver operating characteristic curves (AUC) test is 0.988 for current prediction and future projection [53,54]. Our calculated AUC value shows that the results are significant (0.930) (Fig 4A). AUC value was considered as best and significant on the basis of knowledge based on previous studies [55,56] Current and future habitat suitability of Brachyponera nigrita in Pakistan. The Maxent model provided the species future distribution under hypothetical projected climatic values for 2050 (RCP 4.5 and RCP 8.5) which were compared with bioclimatic values for current species distribution. The results (Figs 8-10) showed that averaged future predictions from bioclimatic Maxent model of 2050 (RCP 4.5) and (RCP 8.5), there would be an increase in moderately suitable area and highly suitable areas as compared to current distribution in Pakistan. In future this species is going to increase its suitable areas especially in KPK, AJK Punjab, Sindh and Balochistan. Its current predicted distribution is in some areas of Punjab, Kashmir, Sindh and Balochistan. This distribution is predicted to be with substantial increase in future in various areas of KPK, AJK Punjab, Sindh and Balochistan.
Current and future hypothetical habitat suitability maps of Pothwar region. The Maxent model provided the species future distribution under hypothetical projected climatic values for 2050 (RCP 4.5 and RCP 8.5) which were compared with bioclimatic values for current

Discussion
Brachyponera nigrita is distinct from B. luteipes and B. jerdoni in having longer joints of flagellum of antennae rather than broad [34]. Collected specimen were compared with published description of Bingham [21] and found similar. During collection, it was observed that workers were foraging at algae affected tree in forest. Schmidt and Shattuck [57] reviewed that Ponerine ants prefer to make their nests in soil, leaf litter or decaying wood. As it is clear from our results too that nests of B. nigrita were observed from same type of habitats. Additionally,  our studies indicates that this species have some sort of association with other organisms like isopods, collembolan, termites etc. This indicates that this species may be predatory in nature too as discussed by Schmidt and Shattuck [57]. Its nesting habit in soil show that it likes low temperature and high humidity, which is also evident from our results of modeling too. It will move to low temperature areas in future. Results of our model showed that geographical distribution of the species is dominantly influenced by four bioclimatic variables namely Bio4, Bio18, Bio17 and Bio15. Bio4 and Bio18 are also important for permutation. From these results of Maxent model, it can be concluded that B. nigrita would be increased under the influence of future climate change. Maxent results for current species distribution map reveal that B. nigrita is predicted to be increased in its distribution pattern in future under future predicted bioclimatics 2050 (RCP 4.5 and RCP8.5) in Pothwar region (Figs 8-10). Currently its distribution is limited in Attock and Chakwal region but in future its suitable areas are going to be increased as a whole in both districts in addition to increase in Rawalpindi and Islamabad. Highly suitable area of this species will also increase in Jhelum district of Pothwar region. Our studies are in line with other studies done in various parts of the world on ant's distribution. It has been proved that ant's distribution is largely limited by climate so is considered as most important factor globally affecting their distribution [58][59][60]. In various studies, species distribution models have been used to evaluate the potential impacts of future climate on habitat suitability and ant's distribution. Kwon et al. [61] modeled the future distribution range of different ant species using generalized additive models and described that 12 species would shrink their suitable habitats, whereas five species would expand their distribution. Bertelsmeier et al. [62] and Bertelsmeier et al. [63] modeled 15 of the most invasive ant species for future and found that 8 would experience decrease in suitable areas and 5 would expand their distribution under future climate. According to retrospective analysis done by Needleman et al. [43] fire ants are expected to have an expanded range distribution due to predicted warming trends.

Conclusions
B. nigrita prefer to make nest in soil, under stones, under the bark of tree etc. It may have some association with other organisms living in its nest. White eggs, pupae and adults may be further studied for estimating its associations with other organisms in the nests. There is a need of lab based studies for determination of its ecological role in various habitats. It may be concluded from the future predicted maps by Maxent model, future climate change would influence the potential distribution of B. nigrita. There will be a total increase in highly suitable habitat of the species in various localities of Pothwar region as well as different provinces of Pakistan with future bioclimatics. Its increased distribution will ultimately affect the total ecology of different organisms in various habitats in Pakistan. So it is recommended that species distribution modeling and GIS tools should be used to predict the impacts of future climatic conditions on various organisms as they are being used worldwide to cope with changing climatic conditions.