Genetic diversity and population structure of the Saharan honey bee Apis mellifera sahariensis from southeastern Morocco: introgression assessment and implications for conservation

This study assessed the conservation status of the Saharan honey bee (Apis mellifera sahariensis) from southeastern Morocco using 12 microsatellite loci to examine genetic diversity and hybridization with other subspecies. Samples from 148 colonies were clustered into seven populations representing the expected distribution of A. m. intermissa and A. m. sahariensis, and reference samples from two European subspecies, A. m. carnica and A. m. mellifera, were included. Moroccan honey bees showed higher genetic diversity than European reference samples, and genetic structure analysis revealed two distinct clusters in Morocco separated by the High Atlas Mountains (FST = 0.05). Although high rates of hybridization with A. m. intermissa jeopardize the genetic integrity of the Saharan honey bee, no evidence of introgression was detected from the European reference subspecies. Additionally, we found that the probability of assignment to Saharan subspecies decreased with increasing human management intensity and precipitation. These findings are important for developing a conservation strategy for the Saharan honey bee in Morocco.


INTRODUCTION
The importance of the western honey bee (Apis mellifera L.) as a substantial insect pollinator is well established and hugely emphasized (Hung et al. 2018;Papa et al. 2022). Several studies have proved the significant contribution of this species to humans' food security by pollination services in commercial crops and many natural habitats (Klatt et al. 2014). Natively found in Africa, Europe, and western Asia, honey bees are currently subdivided into at least 30 different subspecies (Ruttner 1988;Engel 1999;Requier et al. 2019), which have been assembled into seven evolutionary lineages: Africa (A, L, and U lineages), Europe (M and C lineages), and western Asia (O and Y lineages) (Dogantzis et al. 2021). Unfortunately, this diversity is lately confronted with multiple stressors, 1 3 31 Page 2 of 15 which have led to the registration of high levels of colony losses worldwide each year (Gray et al. 2019(Gray et al. , 2020(Gray et al. , 2022. One of the possible causes of these losses is the introgressive hybridization with non-native subspecies, where foreign queen importation and migratory beekeeping are considered to play a key role in enhancing this phenomenon (De la Rúa et al. 2009). In fact, the hybridization of local populations of honey bees adapted to local flora and climate with allochthonous honey bees is not without consequences for the survival of such populations (De la Rúa et al. 2009;Büchler et al. 2014).
There is a growing concern that introgressive hybridization will undermine the genetic integrity of indigenous taxa to the level of extinctions (Ellis et al. 2018). In fact, hybridization and introgression are considered problematic for the conservation of honey bee subspecies as these two processes can lead to a combination of issues threatening their existence (De la . In fact, the modernization of the beekeeping sector and beekeeping practices such as the importation of foreign queens, large scale queen breeding, and the regular movement of colonies, aiming at increasing and diversifying honey production, has heavily impacted the integrity of the genetic pool of locally adapted subspecies and causes genetic pollution of these varieties by introgression (De la Rúa et al. 2009).
In Morocco, there should be present two welldifferentiated subspecies: the Tellian honey bee Apis mellifera intermissa (Maa 1953) and the Saharan honey bee Apis mellifera sahariensis (Baldensperger 1932) (Ruttner 1988). Although the morphometric study of Cornuet et al. (1988) suggested the existence of a third subspecies in the extreme north, Apis mellifera major (Ruttner 1975), there is little evidence to support it . With its remarkable red-yellow coloration, the Saharan honey bee is natively found inhabiting southeastern Morocco's oases (Ruttner 1988). It is separated from the other African subspecies by the Sahara Desert and from the Tellian honey bee by the High Atlas Mountain. The important barrier of the High Atlas separates the Sahara from the rest of Morocco, thus isolating and, at the same time, preventing the local subspecies of the honey bee (A. m. sahariensis) from mixing with those of the North (A. m. intermissa). Unfortunately, several incidences during the 1980s have led to the importation of a large number of colonies, mainly from the subspecies A. m. intermissa, into the natural range of the Saharan honey bee (Aglagane et al. 2022). Currently, with the modernization of the beekeeping sector in Northwest Africa and the practice of transhumance, this subspecies is assumed to be at an advanced status of hybridization (Cornuet et al. 1988;De la Rúa et al. 2007;Chahbar et al. 2013;Loucif-Ayad et al. 2015).
Microsatellite loci or simple sequence repeats (SSR) are an excellent genetic marker for studying genetic diversity, structure, and hybridization in honey bee subspecies (Chahbar et al. 2013;Loucif-Ayad et al. 2015;Eimanifar et al. 2020;Tanasković et al. 2022) because they are codominant with Mendelian inheritance (biparental transmission), highly polymorphic (individual variability), multi-allelic, and neutral (not subject to natural selection) (Estoup et al. 1995). In this study, we have evaluated, using SSR, the genetic diversity and population structure of honey bee populations that were sampled across the natural range of the Saharan honey bee A. m. sahariensis. This geographic race has gained very limited attention in terms of studying its genetic diversity, population structure, and conservation status in Morocco. The only studies dated back to 1988, and 2001, respectively (Cornuet et al. 1988Garnery et al. 1995;Franck et al. 2001), while the most recent one (De la Rúa et al. 2007) was conducted using a limited number of both samples and loci. Since then, the field of apiculture in Morocco has changed drastically especially with the modernization of beekeeping practices, increasing exercise of transhumance, and the importation of foreign queens and colonies. Thus, the mean driver of our study was to provide solid proof that this subspecies is at an advanced state of introgressive hybridization particularly with its immediate neighbor A. m. intermissa and to provide the basic elements needed to start a national wide conservation strategy for the preservation of Moroccan honey bee subspecies. Additionally, since C-lineage honey bees are prevalent in apiculture worldwide, it can be expected that they were imported to Morocco and that traces of hybridization with these C-lineage can also be detected in the genetic makeup of local populations. We hypothesized that hybridization would be more likely in areas with more intensive human management.

Sampling
One worker honey bee per colony was sampled from the nest of 98 A. m. sahariensis colonies between 2019 and 2020. The samples represent four populations from southeastern Moroccan provinces, namely: Errachidia (N = 21 colonies), Ouarzazate (N = 29), Tinghir (N = 27), and Zagora (N = 21) (Figure 1). A detailed justification of the choice of the four populations can be found in Aglagane et al. (2022). For comparison and to test our hypothesis that the Saharan honey bee is at an advanced state of hybridization, 50 worker honey bees representing the Tellian subspecies A. m. intermissa were also collected directly from flowers at 10 localities from two regions located north of the High Atlas Mountains (Marrakech, N = 25 and Tanger, N = 25; each sample represents one colony). To enable the analysis of honey bee samples collected across different regions of Morocco, a more natural division than the boundaries of provinces was required to group them into separate populations. For this purpose, we used K-means clustering based on geographic coordinates (longitude and latitude) and two significant climate features (annual mean temperature and annual precipitation, i.e., WorldClim variables BIO1 and BIO12, respectively, downloaded for Morocco using "geodata" R package ver. 0.5-3). The variables were standardized before analysis. We utilized the Elbow method as a clue to determine the optimal number of groups for our K-means clustering. To balance the level of detail in our analysis with the need for sufficiently large groups, we decided to distinguish seven populations ( Figure 1

DNA extraction
The samples were preserved in 90% ethanol and, after bringing them to the laboratory, stored in the refrigerator (− 20 °C) until DNA extraction, which we performed using the Insect DNA Kit (Omega Bio-Tek, Norcross, GA, USA) following the manufacturer's protocol.

Microsatellite analysis
In the study, we used 13 nuclear SSR loci (Solignac et al. 2003), amplified in two multiplex reactions: multiplex 1: Ap66, Ap43, A7, A88, A113, and Ap55; multiplex 2: Ap243, Ac011, Ap249, Ap090, Ap238, Ap103, and Ap226. To enable simultaneous detection of these loci during capillary electrophoresis, forward primers for these loci were 5′-labelled with fluorescent dyes. Multiplex PCR was performed using the Multiplex PCR Kit (QIAGEN, Hilden, Germany) following the recommended protocol in a final reaction volume of 10 μL (5 μL of 2 × QIAGEN Multiplex Master Mix, 4 μL of primer mix, and 1 μL of template DNA). The PCR cycling started with an initial incubation at 95 °C for 15 min. It was followed by nine touchdown cycles of 94 °C for 30 s, 60 °C (− 0.5 °C per cycle) for 1 min 30 s and 72 °C for 1 min and 24 cycles of 94 °C for 30 s, 55 °C for 1 min 30 s and 72 °C for 1 min. Finally, tubes were incubated at 72 °C for 10 min. To perform amplifications, PTC200 thermocycler was used (BioRad, Hercules, CA, USA). The separation of fragments was performed on an automated sequencer ABI PRISM 3130xl (Applied Biosystems, Foster City, CA, USA) using the internal size standard (LIZ 600; Applied Biosystems). The resulting electropherograms were scored using GeneMarker ver. 2.6.3 (SoftGenetics, State College, PA, USA).

Population genetic analyses
We started our analysis by checking whether the allele frequencies at SSR loci could be biased due to the presence of null alleles or genotyping failure. These issues represent a significant challenge in population genetic studies, as they may lead to the overestimation of population differentiation due to reduced gene diversity. By using INEst ver. 2.153 (Chybicki and Burczyk 2009), we evaluated the possibility of null allele presence and the genotyping failure rates for each of the nine populations separately and also for Tellian and Saharan samples treated as single populations. Posterior distribution was based on 500,000 steps, sampling every 50 steps and discarding the first 50,000 steps as burn-in. Since highly elevated inbreeding levels seem unlikely in honey bees considering polyandry and almost panmictic mating mode at drone congregation areas (Baudry et al. 1998), we assumed that excess of homozygotes in our population could result from null alleles and/or genotyping failure. Thus, we evaluated the following models: n (null alleles), b (genotyping failure), nb (both null alleles and genotyping failure), and null (no null alleles and genotyping failure). We chose the best model based on the deviance information criterion (DIC; Chybicki et al. 2011).
Then, after eliminating possibly biased loci, we calculated several statistics to compare genetic variability in populations: the mean number of alleles (Na), the allelic richness (AR), expected (He), observed heterozygosity (Ho), and fixation index (F IS ), using "hierfstat" ver. 0.5-11 (Goudet 2005). The genetic parameters obtained were compared between Saharan, Tellian, and European reference samples using analysis of variance (ANOVA) using SPSS ver. 22 software. In addition, since allelic richness may underestimate the genetic diversity in the presence of rare alleles, allele accumulation curves over sampling effort were used to illustrate the rate at which new alleles were sampled. The Chao index was used to estimate the expected total number of alleles based on the rarefaction analysis. The analysis was performed in "iNEXT" ver. 3.0.0 (Hsieh et al. 2016).
We tested departure from Hardy-Weinberg equilibrium (HWE) with the exact test based on Monte Carlo permutations of alleles implemented in "pegas" ver. 1.2 (Paradis 2010). The Bonferroni correction was used to adjust for multiple testing with p.adjust function from "stats" package ver. 4.2.2. To measure the extent of divergence between the sampling locations and the four subspecies, we estimated pairwise F ST values following the Weir and Cockerham (1984) method using the "hierfstat" ver. 0.5-11 R package. To obtain 95% confidence intervals, we used a bootstrapping approach over loci.

Population structure and admixture analysis
We used STRU CTU RE ver. 2.3.4 (Pritchard et al. 2000), a Bayesian model-based clustering method, to evaluate the population ancestry and possible gene flow between the populations of honey bees studied. We selected the admixture model because it assumes that individuals may have ancestry from multiple populations. To enhance the accuracy of population assignments, we inferred Lambda (the allele frequencies parameter) separately for each source population, as allele frequencies may differ between honey bee lineages. The initial value of Lambda was set to 1.0, which is a common practice as it allows for flexibility in the estimation process (Porras-Hurtado et al. 2013). The analysis did not include any prior information about the sample origin and was performed using 250,000 burn-in steps and 750,000 Markov Chain Monte Carlo iterations. We ran 16 independent runs for each K value. To infer the optimal K (i.e., the most likely number of ancestral gene pools that contributed to the gene pool of the studied populations), we utilized KFinder ver. 1.0 (Wang 2019). This software computes the parsimony index (PI) along with the widely used ΔK statistic (Evanno et al. 2005). Then, bar plots for clustering modes were produced using "pophelper" ver. 2.3.1 (Francis 2017). Levels of introgression for each sampled population were investigated using the Q-values. This approach uses the Q-index generated by STRU CTU RE for each individual for K = 4. Interpretation of Q-values was conducted following five threshold values of 80, 90, 95, 98, and 99% as described in (Ellis et al. 2018). An individual is then considered a "pure Saharan honey bee" or "pure Tellian honey bee" when the Q value for the A. m. sahariensis or A. m. intermissa cluster equals or exceeds the threshold value of choice.
Due to the restrictive assumptions of the population genetic model used by STRU CTU RE, which are often not met in real populations (such as populations at Hardy-Weinberg equilibrium and linkage equilibrium between loci), we additionally utilized discriminant analysis of principal components (DAPC) (Jombart et al. 2010) to examine population subdivision. The "adegenet" ver. 2.1.10 package (Jombart 2008) was used for these analyses, which were performed on several datasets (1) all Moroccan populations and European reference subspecies, (2) only the seven Moroccan populations separately, (3) selected "pure" individuals (Q ≥ 0.90) of the four subspecies (two Moroccan and two reference European). DAPC was run using principal components representing 90% of the variation and discriminant functions representing n-1 groups. Finally, scatter plots were produced for visual inspection of the clusters. The function adegenet::find.clusters() was first used to infer the number of groups and group membership of individuals within study population.

Relationship between assignment probabilities and human management
We conducted a regression analysis to explore the relationship between honey bee introgression and human management. To model the assignment probability, which is a proportion between 0 and 1, we used a generalized linear model with a binomial family (Zuur et al. 2009). The human footprint index (HFI) was utilized as an explanatory variable to estimate the intensity of human management in each sampled location. HFI is a global dataset that measures the impact of human activity on the terrestrial landscape based on factors such as population density, land use, and infrastructure development (Sanderson et al. 2002). To account for the potential impact of climate on the relationship between assignment probabilities and HFI, we included bioclimatic variables BIO1 and BIO12 (mean annual precipitation and mean annual temperature, respectively, Fick and Hijmans 2017) in the model in addition to HFI. By incorporating these variables, we were able to investigate potential interactions and confounding effects between human management and climate on honey bee introgression. We used the extract function from "raster" ver. 3.6-20 package to obtain the values of the independent variables for each sampling location. To facilitate coefficient comparison, we standardized the explanatory variables measured on different scales, transforming them to a mean of zero and a standard deviation of one. To explore the potential influence of the explanatory variables, we used the dredge() function from "MuMIn" ver. 1.47.5 package (Bartoń 2023) to test all possible combinations of the three variables. We ranked the models using the Akaike information criterion (AICc), and the best-supported models with delta AICc < 2 were averaged to obtain parameter estimates.

Genetic diversity
The INEst analysis with the BIC criterion pointed to the n model (model with null alleles) as the best model explaining homozygote excess in all Moroccan populations. On the contrary, the best model for two European reference populations (A. m. carnica and A. m. mellifera) was the null model, i.e., the model without null alleles and/or genotyping failures, indicating that in this case, neither the null alleles nor the amplification failures were present in a significant proportion.
According to INEst results, the frequency of null alleles in Moroccan samples were negligible for all loci except for locus A7, for which a significant proportion of null alleles has been estimated for some populations: . Since null alleles could confuse population structure analyses, we omitted locus A7 and kept 12 loci for subsequent genetic diversity and population structure analyses.
The overall genetic diversity parameters are presented in Table I. Moroccan honey bees had a higher level of genetic diversity than the European subspecies used as reference samples. For Moroccan populations and 12 loci analyzed, the average allelic richness per locus ranged from 7.02 (North Morocco) to 8.57 (West Sahara), while it was equal to 5.48 and 3.95 in A. m. mellifera and A. m. carnica, respectively (Table I). A similar trend was present in all other genetic parameters showing higher values in Saharan populations over the other reference subspecies. Total allelic richness extrapolated with the Chao estimator showed an increase in the overall number of new SSR alleles sampled based on the number of genotyped individuals, with Moroccan subspecies showing higher rates over European reference samples (Figure 2).
Populations of the Saharan honey bee from southeastern Morocco displayed higher genetic diversity than the Tellian populations. However, the results of the ANOVA test showed no significant difference between these two subspecies. On the contrary, both Moroccan subspecies exhibited higher and significant (p < 0.05) differences when their genetic parameters (Ho, He, and AR) were compared with the European reference subspecies. The exact test revealed no departures from Hardy-Weinberg equilibrium in any of the studied populations.
To evaluate the genetically differentiated populations within our samples dataset, multilocus average F ST values were estimated among the Saharan and Tellian populations and among Moroccan and reference samples of the two European subspecies. Within the Moroccan populations, pairwise F ST values ranged between 0.00 and 0.04, indicating little to moderate genetic differentiation among these populations (Figure 3). Surprisingly, two Tellian populations, West Morocco and North Morocco, showed the highest pairwise F ST among Moroccan populations (0.04). However, this could be attributed to strong introgression in Page 7 of 15 31 West Morocco. When populations were grouped by subspecies and hybrids were omitted from the analysis, we found a moderate differentiation between A. m. sahariensis and A. m. intermissa (0.05, p < 0.001). In addition, both Moroccan subspecies showed strong differentiation from A. m. carnica and A. m. mellifera (pairwise F ST values from 0.24 to 0.37).

Population structure
The STRU CTU RE analysis indicated that the most likely number of clusters for the studied populations was four (K = 4), as supported by both the highest parsimony index and ΔK values (Figure 4). At K = 4, the Moroccan honey bee populations were assigned to two clusters. One cluster predominantly comprised samples from West and North Morocco, which could be interpreted as A.m. intermissa, while the other cluster was more common in samples from Saharan populations and identified as A.m. sahariensis (Figure 4). On the other hand, each of the two European reference subspecies has clustered separately. Levels of introgression into the Saharan honey bees from the Tellian and European reference subspecies were estimated with Q-values, and the results are shown in Table II). The choice of the most stringent threshold (99%) reported no pure Saharan honey bee colonies, while only five colonies (5.10% of total samples) were considered pure when applying a threshold of 98% (Table II). It is worth mentioning that four Tellian samples were identified in the natural range of the Saharan honey bee when considering a threshold of 95%. On the other hand, no significant introgression of C and M lineages was detected in Moroccan populations.
Discriminant analysis of principal components revealed that three to four clusters were the best model describing the population structure of our dataset ( Figure 5). The three clusters correspond to a deep division into evolutionary lineages A, C, and M. The four clusters obtained represent subspecies A. m. carnica, A. m. intermissa, A. m. sahariensis, and A. m. mellifera. The two European reference subspecies have clustered separately from each other and Moroccan samples; however, these latter showed highly overlapping groups ( Figure 5). Results of inferred membership based on the find. cluster function for K = 4 have shown that within the natural range of the Saharan subspecies, 64% of sampled colonies could be assigned to A. m. sahariensis ( Figure 6). Unfortunately, the method was not able to identify hybrids and assigned individuals to either A. m. sahariensis or A. m. intermissa subspecies. According to both STRU CTU RE and DAPC, the West and South Sahara populations showed the least introgression ( Figure 6).

Relationship between assignment probabilities and human management
Based on the Akaike information criterion, the model that contained HFI and precipitation was the highest-ranked model with a weight of 0.397. The model with HFI, temperature, and precipitation was the next best model with a weight of 0.368, followed by the model with HFI and temperature, which had a weight of 0.187. Since all three models had a delta value of less than 1.5, they were considered equivalent (Table III). The combined weight of these three models was 0.952, indicating that the remaining five models, which were based on different combinations of variables, had poorer performance. (Table III) presents the AIC-weighted modelaveraged parameter estimates obtained from the top-performing models. The averaged parameter estimates (β) were − 0.719 (p = 0.01) for HFI, − 0.605 (p = 0.04) for BIO_12, and 0.203 (p = 0.3, ns) for BIO_1. Therefore, the analysis revealed a significant negative effect of HFI and precipitation on assignment probability, indicating that honey bee populations in drier areas with lower human management intensity were more likely to belong to the putative A. m. sahariensis genetic cluster (Figure 7).

DISCUSSION
As expected, genetic diversity in Moroccan populations was significantly higher than in European reference subspecies represented by A. m. mellifera (Lineage M) and A. m. carnica (Lineage C). The same result was previously reported from numerous published papers showing high diversity among African honey bees over European ones (Franck et al. 2001;Loucif-Ayad et al. 2015;Miguel et al. 2015;Techer et al. 2016Techer et al. , 2017. African honey bee populations exhibit two exceptional characteristics -pronounced migratory behavior and higher swarming tendency -suggesting the presence of larger effective population sizes, hence the higher genetic diversity (Franck et al. 1998(Franck et al. , 2001. Our results also suggest higher but not significant genetic diversity when comparing southeastern populations (Saharan honey bee) with populations north of the High Atlas mountain (Tellian honey bee). The same trend was also reported in other publications dealing with the same subspecies (Franck et al. 2001;De la Rúa et al. 2007;Chahbar et al. 2013).
The analysis of population differentiation based on multilocus pairwise F ST has revealed a moderate difference between Saharan and Tellian populations (Fst ≈ 0.05 when hybrids were omitted from the computation). Our findings support  those of Chahbar et al. (2013), who also reported an F ST value of 0.05. However, if all collected individuals are included in the analysis, the F ST estimates are lower (usually, in the range of 0.01-0.02), which is a clear indication of the progressing homogenization of honey bees in the region. In contrast, the two Moroccan subspecies displayed high levels of distinctiveness from the two European reference subspecies, which was also observed for the Algerian populations (Loucif-Ayad et al. 2015). West Sahara population showed the highest differentiation from populations of the Tellian region (0.03-0.04), which indicates it as the least population affected by hybridization as was also confirmed by STRU CTU RE and DAPC.
Population structure analysis has proved the presence of four clusters in our dataset. This result was supported with both methods used (STRU CTU RE and DAPC). Analysis of the STRU CTU  7. Visualization of the highest-ranking model for predicting the probability of assignment to A. m. sahariensis based on A human footprint index (HFI) and B annual precipitation (bio_12). The solid lines represent the fitted relationship between each predictor and the probability of occurrence, and the shaded areas represent the 95% confidence intervals.
RE results indicated that when the number of clusters was set to K = 3, we identified three distinct groups separating the Moroccan populations (lineage A) together apart from European reference samples (lineages C and M). This result corroborates the previously identified lineage relationship (Arias and Sheppard 1996). When K = 4, two distinctive clusters emerge among Moroccan samples. The first one predominated north of the High Atlas mountain (the Tellian region), while the share of second cluster increased towards southeastern Morocco (the Saharan region). The same separation was previously identified for the Algerian honey bees proving a south/north differentiation (Loucif-Ayad et al. 2015). We suggest that the two clusters may correspond to the subspecies A. m. intermissa and A. m. sahariensis, which agrees with the proposed earlier distribution of the two subspecies. It is not possible to use morphological data to confirm this finding because the reference data from Oberursel Data Bank differ markedly from the current population from Morocco (Aglagane et al. 2022). We detected no significant admixture of M and C lineages into A lineage from Morocco, which might be unexpected given that reports have stated that foreign importations have occurred, especially during the French colonization of Morocco (Estoup et al. 1995). Previous mtDNA analysis also failed to confirm the presence of European mitotypes in this part of northern Africa , which is consistent with our results. On the other hand, Estoup et al. (1995) reported European microsatellite gene alleles in a population from Tiznit (A. m. intermissa).
As for the levels of introgression between the Saharan and the Tellian subspecies, results are quite alarming, indicating extremely high levels of admixture. Based on the Q thresholds established at 0.80, only 43 (43.88%) out of 98 colonies were clustered as "pure Saharan" against 55 colonies identified as either "introgressed" or "Tellian." The number of "pure Saharan colonies" decreases as the Q-thresholds go more stringent. This situation was similar to that of Ellis et al. (2018) in which they found 15 pure A. m. mellifera colonies out of 43 at Q = 0.8, ending up with only 4 pure colonies at Q = 0.99. These findings were also corroborated when using the find.clusters function of the adegenet package, where 64% of colonies in Sahara were assigned to A. m. sahariensis. These results strongly highlight the artificial sympatry of Moroccan subspecies as a major threat and highpoints an urgent demand for a conservation management strategy to stop further hybridization and to start a rehabilitation plan for the Saharan subspecies. It is also interesting to note that there are signs of introgression from the Saharan honey bee into A. m. intermissa especially in the population of West Morocco, i.e., in the Marrakech region (although not as pronounced as in the Saharan region). This situation, even though it may increase the withinsubspecies genetic diversity by the introduction of new alleles (Harpur et al. 2012), will eventually decrease the among-subspecies differentiation leading thus to a large-scale homogenization of the two geographic races and, inevitably, the loss of local adaptations (De la ). In particular, A. m. sahariensis is known for its adaptive capacity to support extreme and harsh climatic conditions with temperatures ranging from − 8 to 50 °C (Ruttner 1988). This adaptability is enormously crucial to the survival of honey bee populations in southeastern Morocco, especially with the current rapid climate change in the area. Thus, failure to take appropriate measures to protect local populations of honey bees may lead to the loss of such adaptive traits long-shaped by natural selection .

CONCLUSION
The current study analyzed the genetic diversity and population structure of the Saharan honey bee A. m. sahariensis from southeastern Morocco using nuclear SSR loci. The results revealed the higher genetic diversity of the Saharan honey bee over the Tellian (although not significant) and the European reference (highly significant) subspecies. The population structure of Moroccan honey bee populations indicated the presence of two distinct clusters representing A. m. sahariensis and A. m. intermissa, south and north of the High Atlas Mountains, respectively. As for the levels of introgression, we proved that the Saharan honey bee is at an advanced state of hybridization, especially with its immediate neighbor, A. m. intermissa. This situation indicates the pressing demand for a conservation strategy for this remarkable geographic race. Our study highlights that both human management and climate (precipitation) could influence chances of survival of the local Saharan honey bee.