BEES Potential Niche Modeling Distribution and Wing Geometric Morphometrics of Apis mellifera in the Brazilian Pantanal

The natural world varies across space with a mosaic of habitat fragments, and for any given species, some of these fragments are suitable and others are not, with the range of physical conditions where species can persist being called the fundamental niche (Pianka, 1981). While the fundamental niche presents all the characteristics of the n-dimensional space essential for a species in the absence of other organisms, the realized niche is the portion of the fundamental niche Abstract The ecological niche models can be important for biogeographic patterns and processes and geometric morphometrics involves identifying changes that have occurred and comparing them to other specimens from different places and/or environmental conditions, assessing whether the environment is influencing such change. The present work aimed to verify the potential model of distribution for Apis mellifera and analyze if there is variation in the geometric morphometrics in wing venation in the Pantanal. We followed the hypothesis that there is variation in the geometric morphometrics of wings and that the geographically closest groups are more similar. For niche modeling, 44 geographical points and 19 bioclimatic variables were used. For morphometrics, twenty-two anatomical landmarks were plotted at the intersection of the veins. The X and Y coordinates were standardized through Procrustes superimposition, and PCA and MANOVA tests were performed. The predictive model indicated that the center of the Pantanal plain shows the greater probability of occurrence for the species. The most important bioclimatic variables were: average temperature in the rainiest quarter (84%) and average annual temperature (72%). Morphometric analyzes indicate that there was variation between the most distant geographic points. The slight variation between some closely located points in the Pantanal can be related to individual reflections of colonies from other points, since the species has great dispersion capacity. Thus, the distribution of A. mellifera in the Pantanal is possibly related to temperature also accompanied by human occupation and the geometric morphometrics of its wings reflecting aspects of dispersion and population dynamics in the Brazilian Pantanal. Sociobiology An international journal on social insects


Introduction
The natural world varies across space with a mosaic of habitat fragments, and for any given species, some of these fragments are suitable and others are not, with the range of physical conditions where species can persist being called the fundamental niche (Pianka, 1981). While the fundamental niche presents all the characteristics of the n-dimensional space essential for a species in the absence of other organisms, the realized niche is the portion of the fundamental niche actually occupied by that species, considering all interspecific interactions in such environment (Hutchinson, 1978).
The ecological niche is the last distributional unit for a species or subspecies (Grinnel, 1917) that unites populational distributions with their environment, with species distributions reflecting the range of physical conditions that support individual survival (Braunisch et al., 2008). In this context, the development of ecological niche models can help to investigate biogeographical patterns and processes, predicting the geographical distribution of species from sparse occurrence  (Guisan & Thuiller, 2005). Species distribution models are empirical precepts that relate field observations to environmental predictor variables, based on statistically or theoretically derived response surfaces (Guisan & Thuiller, 2005). Modeling plays an important role, as it addresses several issues such as the analysis of an exotic/invasive species (Peterson & Vieglais, 2001;Giovanelli et al., 2008) and implications for conservation (Guisan & Thuiller, 2005).
Animals are free to occupy different habitats, however, environments can change in time and space, and individuals must adjust by changing their behavior, physiology and even their structure (Jungers et al., 1995). The measurements taken from structures where shape is studied contain information about both size and conformation (Richtsmeier et al., 2002). Geometric morphometrics aims to describe and represent the geometry of the conformations studied, clearly describing and locating regions where changes occur and graphically reconstituting such changes in shape (Zelditch et al., 2012).
Geometric morphometrics involves identifying changes in homologous structures and comparing them to other specimens from different places and/or environmental conditions, making it possible to assess whether the environment is influencing such change (Francoy et al., 2011). In the last decades, different applications in biology have used geometric morphometrics (Adams et al., 2004). The studies quantify ontogeny (Carvalho et al., 2011, Komatsu et al., 2018, taxonomy, systematics or evolution of morphological characters (Cardini et al., 2007, Salas-Lopez et al., 2017, sexual dimorphism (Pretorius, 2005;Olivier & Aranda, 2018), ecomorphology, functional and biomechanical issues of biological forms (Jungers et al., 1995) and intraspecific geographic variation (Roggero & Passerin d'Entrèves, 2005;Cardini et al., 2007). Insect wings are very suitable structures for conducting morphometric studies, with many studies indicating geographic variation in wing structure (Hoffmann & Shirriffs, 2002;Hoffmann et al., 2005). Geometric morphometrics is precise enough to measure both greater variations between species and intraspecific variations, identifying colony origins in studies about bees as Apis mellifera (Francoy et al., 2008); an exotic species in Brazil.
African bees Apis mellifera scutellata Lepeletier, were introduced to Brazil in 1956 (Kerr, 1967). About one year later, 26 swarms of this bee with their respective queens, escaped and bred with the other subspecies of European honey bees that were introduced in the 19th century: the Italian A. m. ligustica Spinola, the German A. m. mellifera Linnaeus, the Austrian A. m. carnica Pollmann and A. m. caucasica Gorbachec. With this breeding process, polyhybrid populations called Africanized bees arose, presenting predominantly African bee characteristics, such as great ability to swarm and rusticity, high capacity for defense, adaptation to inhospitable environments and ability to reproduce in a shorter life cycle than the other previously mentioned subspecies (Kerr, 1967). Such characteristics allowed the species to rapidly expand its biomass and significantly increase its population and geographic distribution (Gonçalvez, 1994). Apis mellifera is characterized as having very populous colonies with 5,000 to 180,000 individuals (Lindauer & Kerr, 1960), therefore, it has a high recruitment potential to forage and a high number of bees visiting flowers (Kerr, 1967;Sakagami et al., 1967). The Africanized bee has a high swarm capacity and, due to the favorable conditions it is found in the Neotropical territories, spread throughout the vegetation domains in Brazil, including the Pantanal floodplains.
Considering the characteristics of ecological niche and morphometric adjustments in relation to environmental variations, the present work aimed to verify the potential spatial distribution model of A. mellifera in the Pantanal and to analyze if there is variation in the geometric morphometric of its wing venations, indicating where and how such variations occur and whether they are significant between different locations throughout the Pantanal plain, since the region presents divergent characteristics. As hypotheses, we expected variation in the geometric morphometric of the wings and that the geographically closest groups would be more similar to each other within the Pantanal area.

Study area
The Pantanal, one of the essential South American domains, is the largest continuous floodplain on the planet, with an area of approximately 150,000km 2 . The region's climate is sub-humid tropical, with dry and rainy winters and average annual rainfall between 1,000 and 1,500 mm in the North and between 1,000 and 1,200 mm in the South (Peel et al., 2007). The territory originates in the highlands, with predominance of Cerrado, composed of Brazilian savanna (cerrado stricto sensu), Cerradão, natural fields, floodplains and environments such as freshwater or brackish ponds, rivers and streams (Silva et al., 2000). Furthermore, according to its hydrological, soil and vegetation characteristics, the Pantanal can be divided into 11 distinct sub-regions. The altitude level varies from 80 m to 150 m and was formed by the subduction process of a large area, which occurred simultaneously with the appearance of the Andes Mountains (Mercante et al., 2011).

Ecological niche modeling
For niche modeling, 44 georeferenced A. mellifera occurrence points were used, 19 points where specimens were collected (Aranda & Aoki, 2018) and 25 points of occurrence records through bibliographic search (Table 1). The following Portuguese and English keywords were used for the search: "Apis mellifera", "Pantanal", "European bee"/"abelha europeia", "bee production"/"produção apícola". We used 19 bioclimatic variables from the Worldclim database (http://worldclim.org) to build the niche distribution model. Such data is correlated and allows for the visualization and understanding of the spatial pattern. For this model, the DIVA-GIS georeferencing software (V. 7.5) was used to map and analyze spatial distribution data and elucidate geographic and ecological patterns (Hijmans et al., 2012). We used the BIOCLIM algorithm (Hijmans et al., 2012) which calculates environmental similarity in multidimensional space with spatial resolution of 2.5 arch minutes (~5 km²) to create the predictive map. The result is a georeferenced distribution map that demonstrates the probability of species occurrence in a determined area presenting the basic conditions to maintain such species (Hijmans et al., 2012). The model validation test was performed using the values of the area under the curve (AUC), which is an accuracy measurement used for the prediction that measures the forecasting power based on real species distribution data. The AUC ranges from 0 to 1, where a value of 0.5 indicates that the model is no different than random and a value of 1 perfectly discriminates between presence and absence records (Hijmans et al., 2012). Thus, AUC values above 0.5 denote good accuracy and the closer to 1, the better the model.
The specimens were screened, stored in 70% alcohol and mounted on an entomological pin to remove the wing, being established the removal of the right wing due to the firmness in the position of fixation of the specimens. The left anterior wing of each bee was removed near its insertion base on the thorax with an entomological scalpel, was mounted between blade and cover slip, and photographed with an Opticam®OPT 5.1MP camera coupled to a stereomicroscope (Bel photonics series XTL). Cartesian coordinates in the XY plane for 22 anatomical landmarks were manually plotted at the origin and intersection of the wing veins using OPTHD 3.7 software to construction of matrix data (Fig 1).
In order to observe divergence in wing shape within and between the populations of A. mellifera, the X and Y coordinate data was standardized with Procrustes superimposition to remove the effect of wing size, maintaining the shape in relation to the position of anatomical landmarks only. A Principal Component Analysis (PCA) was performed to identify where and how much was the magnitude of changes in landmarks wings. Thin-plate Splines were performed and warps of digitized x/y coordinates of reference points (landmarks) as a square grid associated deformations using the average shape as a reference, and shown expansion factor (or contraction) of the area around each reference point indicating the degree of local growth. The canculation uses the Jacobian of the deformation, with the color-coded expansions being green for expansion and purple for contraction (Dryden & Mardia, 1998). To test the hypothesis that there is a separation of individuals in relation to their geographically closest groups, the scores of the main axes were tested using PCA components with Multivariate Analysis of Variance (MANOVA) to analyze and test differences between the eight pre-established geographic groups (Sadeghi et al., 2009). The analyses and graphics were performed using the free software Past® 3.17 (Hammer et al., 2001), with an alpha of 0.05 for the tests.

Niche Modeling
The potential geographic distribution of A. mellifera encompasses almost the entire region of the Brazilian Pantanal (Fig 2). The predictive model indicates that the center of the Pantanal plain, is the area with higher probability of A. mellifera presence. In the upper Pantanal (Mato Grosso state), registration points are distributed in areas of low probability of occurrence such as Poconé (dark green), medium probability of occurrence such as Barão de Melgaço, Cáceres and Poconé (light green), and high probability of occurrence as Barão de Melgaço (yellow). In the lower Pantanal (Mato Grosso do Sul state), the central area corresponding to Paiaguás and Nhecolândia presents areas with high probability of occurrence in the Pantanal of Aquidauana (yellow), very high probability of occurrence (orange) in the Pantanal of Miranda-Abobral and Nhecolândia, and excellent probability of occurrence in the Pantanal of Paiaguás (red). The main bioclimatic variables responsible for the A. mellifera model were: average temperature in the rainiest quarter (84%) and average annual temperature (72%) ( Table 2). The accuracy assessment (AUC) of the prediction model indicates that the model has excellent accuracy (AUC = 0.970).

Geometric Morphometrics
The wing landmarks of A. mellifera indivuduals showed differences in the measured positions according to the PCA, where the points with the greatest variations and amplitude were the landmarks present in the X axis region (Mark: 01,04,13,14,18,19 20,21 and 22) of the wing (Fig 3; Table 3).
The central locations of vein intersections show differences on both the X and Y axes, with the warmer colors representing the greatest deformations in the landmark positions (Fig 3). Landmarks number 1 and 22 showed high variation due to the cut at the base of the wing, which ended up making it difficult to mark the points (Fig 3)  Average temperature of the driest quarter ºC 23.8 -24.5 (36%)

19
Precipitation of the coldest quarter mm 86 -108.2 (61%) Table 2. Bioclimatic variables and coefficient of variation (CV) corresponding to importance (%) in predicting the occurrence of Apis mellifera in the Brazilian Pantanal. in morphometric interpretation. We only used them to ensure accurate positioning of the areas for other measurements. Variations in geometric morphometrics between groups were significant, with variation between geographical points (MANOVA -Pillai trace = 3.47; F = 1.78; p <0.05) ( Table 4). Groups G3 and G8 showed the greatest difference, with G3 close to the Serra do Amolar region, which has a slightly higher altitude, and G8 Porto Murtinho, which is far from other points (Axis 1 = 26.59; Axis 2 = 15.74; Fig 4).

Discussion
The map shows a high percentage at high temperatures, with 25.5 ºC to 26.4 ºC and 27.1 ºC to 28.2 ºC. This is due to the topographic alignment arranged in the longitudinal direction, which exhibits well-defined morphological features: the plain to the west and the plateau to the east. This disposition influences the behavior of air masses. On one hand, the depression of the Pantanal plain acts as a corridor for polar air masses during the winter in western-central Brazil. On the other hand, seasons are responsible for the retention of hot and dry air causing its typical muggy weather conditions (Alvare et al., 2014). Temperatures are high throughout the year, with an annual average of 25 ºC, remaining above 26 ºC in autumn and winter, and never going below 30 ºC in spring and summer (Alvare et al., 2014). Vital et al. (2012) report few predictive models in the Pantanal region, due to the ancestral subspecies of Africanized bees and indicate that the best model for the Neotropical region reflects colonization by A. mellifera scutellata, which has a low probability of occupation in the Pantanal. The typical climatic conditions of the Pantanal, with high temperatures and low rainfall in regions where there was better prediction of species occurrence, coincide with high temperatures. However, with a large number of extensive lakes and lagoons (Pantanal region of Nhecolândia and Paiagúas), water availability would not be a limiting factor (Gonçgalves et al., 2011;Mercante et al., 2011). Although high temperatures have been reported as limiting factors for the species (Le Conte & Navajas, 2008), the availability of food resources throughout the year contributes to the species' permanence in the area (Salis et al., 2015).
Registration points are concentrated in transition areas with the best occurrence predictions (High, Very high and Excellent), being more concentrated in regions such as Poconé in the northern and Miranda-Abobral in the southern Pantanal, areas that are typically used for ecotourism and highly populated. Human aspects associated with the use of exotic species, such as honey production with certification from the Pantanal, can cause ecological losses in natural areas that have not yet been measured. Understanding which areas are potentially susceptible to colonization can help mitigate future negative impacts.
The changes found in more distant geographic groups may reflect genetic aspects among the colonies in the Pantanal (Francoy et al., 2011;Combey et al., 2018;Henriques et al., 2020) or even in relation to resource use in case of some  Table 3. Mean values and standard deviation (SD) of the X and Y axes for the anatomical landmarks of the left front wings of Apis mellifera individuals collected in the Brazilian Pantanal. bees Ribeiro et al., 2019). The slight variation between some points in the Pantanal could be related to the fact that some individuals collected were not from the same colony, that is, individual reflections of other distinct groups within the same collection point. A minimum distance of 20 km was established between the collection areas to avoid sampling individuals from one point that may be foraging at another point. However, large bees and those that forage massively can fly long distances in search of food resources (Torne-Noguera et al., 2014), which may have caused similarities in closer locations, highlighting differences in the most distant sampling points only and indicating good dispersion capacity. As expected, the closest geographic points showed greater morphometric similarity. Geographically isolated points, such as the Serra do Amolar region (accessible only by water) with apiaries recently established due to the social actions of the NGO ECOA for income generation, or those with a smaller amount of bee records, such as the Porto Murtinho region, may still represent localities with possible genetic isolation Another possibility could be related to the anthropic management of bees in the Pantanal, which introduces different colonies to the same apiary or variations accumulated over time (Francoy et al., 2009;Francoy et al., 2011). Additionally, migratory beekeeping that has been developed by a few beekeepers could have caused such divergence (Reis & Comastri-Filho, 2003;Francoy et al., 2009, Nunes et al., 2012 due to fast development, adaptation and rusticity. Furthermore, if a colony is removed from one place and installed in another place it could also lead to such effects, what happened when an A. mellifera colony was removed from a school in the Barra do São Lourenço community and relocated to the Dois Corações site, located at Morro Saquarema in the Amolar region (Almeida & Reis, 2016). The removal of nests from natural areas to managed apiaries is a common practice adopted by riverside inhabitants in the Pantanal. Genetic analyzes of A. mellifera populations in the Pantanal may reveal the history and routes of migration and colonization of the areas, since the morphology reflects the genetic aspects and we evidence significant differences in the populations (Miguel et al., 2011).
Thus, the distribution of A. mellifera in the Pantanal is related to temperature in the plain and the geometric morphometrics of wings show significant variation between the sampling points, reflecting aspects of dispersion and population dynamics in the Brazilian Pantanal, mainly in more isolated and distant areas.