Insect surveys were conducted in four regions in Morocco as follows (Figure 1 and Supplementary Table 1). (1) The Casablanca-Settat (CS) region in the north-center of Morocco, a semi-arid region dominated by cereal and legume crops. (2) The Rabat-Sale-Kenitra region (RSK) near Maamora forest (25 km from Kenitra city), characterized by sub-humid climate and small farms producing vegetables and avocado. (3) The Fes-Meknes region (FM), a mountainous zone with semi-arid climate located at the foot of the Middle Atlas in Sefrou province. This latter area is unique in Morocco for having Rosaceae orchards such as apple, plum, and cherry in addition to other crops like vegetables and cereals. (4) The Draa-Tafilalet (DT) arid region, including sites located between 4 and 20 km from Errich city and 8 and 20 km from Errfoud. Thus, sites were located between mountainous and oasis areas and characterized by apple orchard, alfalfa, some vegetables, and date palm (Ministry of Agriculture, Fisheries, Rural Development, Water, and Forests, 2020).

We considered nine main crops: faba bean, pumpkin, zucchini, melon, apple, cherry, tomato, eggplant, and okra (Supplementary Table 1). For each main crop in each region, we selected eight fields of 300 m² (30 m × 10 m) (apple and cherry orchards excepted, see below). Five of the eight fields were under the Farming with Alternative Pollinators management approach (FAP) (Christmann and Aw-Hassan, 2012; Christmann et al., 2017). Under the FAP management, the main crop was planted in the central 75% of the field area and other crops (Marketable Habitat Enhancement Plants, MHEP) were planted at the edge of the fields. The three remaining fields were monoculture of the main crop. For apple and cherry, four orchards per crop (i.e., two regular orchards and two orchards with MHEP) were selected. For each orchard, two rows with seven trees were marked. All crops were planted (or orchards were surveyed) in 2018 and replicated in 2019, except pumpkin and faba bean in some regions (Supplementary Table 1). Crop planting date followed local farmer’s schedule in each region. Most crops were planted in different regions at the same time except for zucchini. Following the farmers’ choice, this was planted in spring in DT and FM regions and in summer in RSK and CS regions.

Overall, we completed 37 trials with 22 crops (including main crops and MHEP, see Supplementary Table 1), in a total of 263 fields.

For each trial, three or four insect sampling sessions were conducted depending on the flowering period of the main crop. Two types of sampling methods were used in each session: (1). Transect walks, including observations, and sweep net (about 30–35 cm in diameter) and vacuum collection on flowers (Popic et al., 2013), and (2). Trapping with colored pan traps (white, blue, and yellow) (Westphal et al., 2008).

Four transects were completed in each sampling session as follows: two inside the main crop for 10 min (5 min each), a third in the outer 25% of the field (i.e., in main crop or MHEP, for 5 min in 2018 and 10 min in 2019) and a fourth in the flowering plants (wildflowers and other crops) in the field margin outside the crop (for 5 min). The insect survey in orchards (apple and cherry) consisted of walking alongside trees of the two plots during 10 min, and visits in MHEP were assessed within a 10 min transect. All insects visiting the crops’ flowers were collected or noted if they could not be collected. Pan trapping was done by using two sets of plastic bowls with three colors (blue, white, and yellow). Pan traps were filled with water and a drop of liquid soap. The three sets of pan traps were put in each field in each sampling event for 2 days (one set inside the field and the other outside the field). Floral visitors were identified visually to species level for the bee species Apis mellifera, Bombus terrestris, and Xylocopa pubescens, all which can be unambiguously determined in the field.

Specimens collected from transects and pan traps were prepared and labeled for the following insect groups: Hymenoptera, Aculeata (i.e., bees and wasps), Lepidoptera (i.e., butterflies), and Diptera, Syrphidae (i.e., hoverflies) (Ollerton et al., 2011). Individuals were identified in the laboratory to genus level for bees (Michez et al., 2019) then sent to expert taxonomists for identification to species level (see “Acknowledgments” section). The remaining groups were identified to the lowest taxonomic level possible in the laboratory (Borror and White, 1991).

We explored the records from both FAP and control fields. We assume that MHEP did not affect the crop pollinator community in FAP fields based on previous research in the same sites. This study has demonstrated using visitation data and pollen load analysis that MHEPs affected neither the crop pollinator’s abundance and richness nor the crop pollinator’s pollen diet (Sentil et al., 2021, 2022).

Initially, we classified insects recorded in this study into five insect groups: wild bees, honeybees, non-bee Hymenoptera, Syrphidae, and butterflies. The mean and the standard error of the relative abundance of each insect group was calculated using (1) the data of all insects recorded in agroecosystems and (2) the data of insects recorded visiting studied crops (22 different crops). We used the relative abundance for each insect group because the sampling effort was not standardized between the different agroecosystems and crops.

For the following analyses, we focused on wild bees. Data from all fields, all sampling dates and both years were pooled to generate a species accumulation curve using the “vegan” package in R. A heatmap was constructed using the “pheatmap” package to visualize interaction between bee genera and crops using the relative abundance of each genus visiting each crop.

To compare regions, we performed the analysis for each of the main crops which were studied in more than one region (zucchini, pumpkin, faba bean, and eggplant). We used data from the two transects sampled inside the main crop and only three sampling events for each crop in each year to remain with coherent sampling effort across all crops and regions. Differences in bee species richness and abundance were analyzed using a generalized linear mixed model “glmer.nb” in the “lme4” R package (Bates et al., 2015). The data distribution used was a negative binomial and we fitted a random intercept model. We used “region” as fixed factor, “year” as crossed random factor and “sampling” as nested random factor within “field.” Then we performed a pairwise comparison between crops and regions using the “glht” function in the “multcomp” R package (Hothorn et al., 2008).

Bee community composition across different regions was analyzed by Non-metric Multidimensional Scaling (NMDS) employing the function “metaMDS” in the R package “vegan.” NMDS was followed by statistical analyses: PERMANOVA (Permutational Multivariate Analysis of variance) and SIMPER (Similarity Percentage Analysis) to identify those species that contributed the most to the observed differences between pollinator communities. Indicator species were also identified using Indicator Value with the “multipatt” function in the package “indicspecies” in R.

A total of 53,361 insects were collected or observed in all agroecosystems surveyed in 2018 and 2019. We recorded a total of 37,317 individuals within crop fields by direct observation and hand netting, of them 3,831 individuals vising wildflowers and additional 12,213 by using pan traps (Supplementary Figure 1). Honeybees and wild bees were the most abundant groups with 37 ± 2% and 45 ± 2% (mean ± SE; N = 263 fields), respectively. Among the other groups, non-bee Hymenoptera represented 12 ± 1%, Syrphidae with 4 ± 0.5% and Lepidoptera with 2 ± 0.5% (Supplementary Table 2). Among the insects visiting target crops, we counted a total of 24,553 honeybees which represented a mean of 49 ± 2%, and 8,400 wild bees that consisted of 33 ± 2% of the mean bee abundance per field (Supplementary Table 2).

Species accumulation curve showed that the number of bee species visiting crops did not reach an asymptote. Basing on the extrapolated values, we sampled as many as 88% of the expected wild bee species richness (Figure 2A and Supplementary Table 3).

A total of 213 species from 39 genera were identified as crop visitors (22% of the national richness) (Figure 2). The most abundant wild bee family visiting crops was the Halictidae (n = 3,940; 47% of the wild bees) followed by Apidae (n = 1,567; 19%), Andrenidae (n = 1,530; 18%), Colletidae (n = 432; 5%), and Megachilidae (n = 368; 4%). The most abundant genus was Lasioglossum (n = 2,046; 29%), followed by Andrena (n = 1,302; 19%), Xylocopa (n = 528; 8%), Amegilla (n = 450; 6%), Nomioides (n = 436; 6%), and Hylaeus (n = 391; 6%). Sixteen genera had an abundance less than 10 individuals and represented less than 1% of the total number of bee genera visiting crops. The most species-rich genus visiting crops was Andrena (53 species), followed by Lasioglossum (35 species), Anthophora (17 species), Hylaeus (16 species), and Eucera (12 species) (Supplementary Table 4).

The heatmap revealed that most of the 22 crops were visited by dominant bee genera like Lasioglossum, Andrena, Xylocopa, and Amegilla, while some rare genera like Pseudapis, Lithurgus, and Melecta were visiting only one crop (Figure 3).

The family Fabaceae represented by faba bean, alfalfa, sweat pea, green bean, and grass pea, was visited mostly by the bee genera Anthophora, Amegilla, Eucera, Megachile, and Xylocopa, but also by some species of Andrena (such as A. abjecta, A. flavipes, and A. fulvicornis) and Lasioglossum. The Cucurbitaceae (zucchini, pumpkin, melon, armenian cucumber, and cucumber) was mostly visited by Lasioglossum, Ceratina, and X. pubescens. Crops of the family Solanaceae, including tomato and eggplant, had X. pubescens as the main visitor. However, we also collected bees from other genera visiting these crops like Amegilla, Halictus, and Lasioglossum but with a lower abundance. The Apiaceae (coriander and anise) and Brassicaceae (canola and arugula) attracted a high number of bees, these crops were visited mainly by Andrena, Hylaeus, and Lasioglossum. Okra (Malvaceae) was visited by a high abundance of Amegilla, Ceylalictus, Lasioglossum, and X. pubescens. In sunflower fields, the only crop which represented the Asteraceae family, we found that the most abundant species was X. pubescens, followed by the family Halictidae (genera Halictus, Lasioglossum, and Nomioides), but also with some bees of the genera Amegilla and Bombus. The family Lamiaceae represented by chia was frequented mostly by Ceratina cucurbitina, and two species of Osmia: O. caerulescens and O. tricornis. Finally, the wild bee species visiting apple flowers included two species of Andrena (A. flavipes, and A. propinqua), three species of Lasioglossum (L. lucidulum, L. malachurum, L. algericolellum), B. terrestris, Eucera nigrilabris, and Anthophora fulvitarsis, while cherry was mostly visited by Lasioglossum (L. algericolellum) and B. terrestris (Figure 3 and Supplementary Table 6).

Accounting to how the abundance and species richness of wild bees visiting each crop differed between regions, the GLMM analysis revealed significant differences for zucchini, faba bean, and eggplant (all p < 0.05, Supplementary Table 7). Zucchini and faba bean in Draa-Tafilalet hosted the highest mean bee abundance (13.05 ± 3.2 individuals for zucchini and 10.08 ± 1.92 individuals for faba bean) and mean bee species richness (3.25 ± 0.54 for zucchini and 5.17 ± 0.65 for faba bean) followed by Fes-Meknes for zucchini (5.36 ± 1.03 individuals; 2.22 ± 0.31 species), and Casablanca-Settat for faba bean (4.49 ± 0.69 individuals, 2.20 ± 0.26 species).

For eggplant, the mean bee abundance and species richness were found to be higher in Rabat-Sale-Kenitra (1.23 ± 0.27 individuals: 0.65 ± 0.12 species) than in Casablanca-Settat region (0.36 ± 0.11 individuals, 0.31 ± 0.09 species) (Figure 4). However, no statistical differences were found between regions in the abundance and species richness of wild bees visiting pumpkin and apple fields (Supplementary Table 7).

Wild bee composition differed between regions for all studied crops (PERMANOVA test: all p < 0.005). Based on pairwise analyses, we found differences in bee composition between all pairs of regions for each crop except between the regions Rabat-sale-Kenitra and Casablanca-Settat for wild bees visiting zucchini (p = 0.15), and between Rabat-Sale-Kenitra and Fes-Meknes for zucchini and pumpkin (p = 0.08 and p = 1, respectively) (Figure 5).

SIMPER analysis showed that X. pubescens and halictid bees, especially Lasioglossum species, were the major contributors to the dissimilarity between regions for zucchini and pumpkin, whereas species of Eucera, Anthophora, Amegilla, and Xylocopa contributed to the compositional differences between regions for faba bean. Comparison between the two regions where eggplant was sampled revealed that only four species had a cumulative contribution of more than 70% to the dissimilarity index (X. pubescens with 53%, Anthophora sp. with 8%, Amegilla quadrifasciata with 7% and Halictus fulvipes with 5%) (Supplementary Table 8).

Based on the indicator species analysis, 14 bee species visiting zucchini were strongly associated with one or several regions. Seven species were associated with Draa-Tafilalet, two with Fes-Meknes and only one species was associated with Rabat-Sale-Kenitra. Four other species were found in more than one region: Lasioglossum villosulum, Lasioglossum leucozonium, and Lasioglossum interruptum in two regions and Lasioglossum malachurum in three regions (Tables 1, 2). Among the bee species visiting pumpkin, only four species were significantly associated with one region (Table 1) whereas Lasioglossum villosulum and Xylocopa pubescens were common and widespread species across regions (Table 2). For faba bean, six species were associated with Draa-Tafilaltet, two species were found to be indicators of the regions Rabat-Sale-Kenitra and Casablanca-Settat (Table 1) and only one species was associated with both Casablanca-Settat and Rabat-Sale-Kenitra (Table 2).

This study shows that agroecosystem in Morocco and particularly crop plants host high relative abundance of wild bee species, i.e., 45% in all agroecosystems, and 33% of all crop visitors, compared to that found in other countries. For example, in a study across five continents, wild bees accounted for 23% of crop visitors (Rader et al., 2016) and in the United States they accounted for 26% of the visitors (Reilly et al., 2020).

Furthermore, we recorded 213 wild bee species as crop visitors, i.e., 22% of national species richness. This number is high compared to the estimated global species richness of bees currently known to be crop visitors, 785 bee species, i.e., 4% of global species richness (Kleijn et al., 2015). The high species richness of wild bees found on crops further supports Morocco’s status as a hotspot of bee diversity (Michener, 1979; Lhomme et al., 2020).

Our results are in accordance with many studies linking bee species and crops. We found that the most commonly represented genera in our study sites, like Lasioglossum, Andrena, and Xylocopa, were visiting multiple crops in agreement with the findings of Hutchinson et al. (2021).

Our results are also in line with studies showing a link between the floral traits of crops and their bee visitors (Garibaldi et al., 2015). We found that long tongued bees like Eucera, Anthophora, and Amegilla were more abundant at Fabaceae flowers that are characterized by flowers with a deep corolla tube (Aouar-sadli et al., 2008; Shebl and Farag, 2015; Marzinzig et al., 2018). Likewise, Solanaceae crops that typically require buzz pollination for fertilization were found to be visited mostly by X. pubescens and A. quadrifaciata, the two species which represent groups of efficient buzz pollinators (Gemmill-herren and Ochieng, 2008; Jayasinghe et al., 2017; Toni et al., 2020; Udayakumar et al., 2021).

Unlike other crops, we found a relatively small number of species and individuals visiting apple and cherry orchards. This could be explained by our sampling of these crops during cold and rainy weather conditions, when only a few species that are more tolerant to suboptimal conditions were active (including the genera Anthophora, Bombus, Eucera, and Lasioglossum). A similar pattern of floral visitation was found along a climatic gradient from Spain to the Netherlands with a growing dominance of the genera Bombus and Lasioglossum (Weekers et al., 2022). The higher abundance of Anthophora and Eucera in our study sites is typical to countries with a warmer Mediterranean climate, such as Morocco and Spain (Patiny and Michez, 2007; Weekers et al., 2022). Particularly in cherry orchards, we recorded mostly bees of the genus Lasioglossum while Andrena were less frequent, contrary to other studies in Germany and Belgium where these bees were the most abundant (Holzschuh et al., 2012; Eeraerts et al., 2020). Additionally, as woody Rosaceae are not a common element of the native Moroccan flora, the Moroccan Andrena community lacks many of the Eurasian species typically associated with these flowering plants, potentially explaining their low abundance in our surveys.

We found that bee abundance, species richness and bee community composition were significantly different between regions for several of the studied crops.

There are several possible factors that can account for this variability between regions. Firstly, the studied regions are different in terms of the landscape composition and crop diversity that surround the surveyed fields. Studies have repeatedly shown that closer proximity of crop fields to semi-natural areas has a positive impact on wild bee diversity via augmentation of nesting resources in depauperate agricultural landscape (Ricketts et al., 2008; Kennedy et al., 2013; Chatterjee et al., 2020; Geeraert et al., 2020). In addition, higher crop diversity in agroecosystems may increase bee density (Raderschall et al., 2021). Secondly, the variability of bee diversity and community composition might be related to the climatic variability among regions in our study (sub humid, semi-arid, and arid regions). Bee species richness was found to be higher in arid and warm regions compared to wet regions (Michener, 2007; Patiny et al., 2009). This finding is in line with our result in Draa-Tafilalet (an arid region) where bee abundance and species richness in faba bean and zucchini fields were higher compared to other regions. Finally, seasonal differences throughout the year (Osorio-canadas et al., 2016). Previous studies have shown differences in bee abundance and species richness between dry and rainy seasons (Oertli et al., 2005; Samnegård et al., 2015; Balfour et al., 2018). These results may explain the differences that we found for zucchini in Draa-Tafilalet and Fes-Meknes where the crop was surveyed in early summer and was associated with a higher number of bee species as compared to Casablanca-Settat and Rabat-Sale-Kenitra, the two sites which were sampled in autumn.

Overall, 26 species were identified as indicator species of the regional crop-visiting bee community. The identified indicator species are among the most abundant species in agroecosystems (Supplementary Table 5) and are widespread in Morocco (Lhomme et al., 2020). Identifying indicator species by region and by crop could form the basis for selecting conservation priority species (Bladt et al., 2008). We found for example that only Lasioglossum species (three species) are indicators in Rabat-Sale-Kenitra and Fes-Meknes for zucchini compared to Draa-Tafilalet, in which seven species from six genera were found: Lasioglossum, Xylocopa, Amegilla, Halictus, and Nomioides. This suggests that different conservation strategies may be required between regions given the identified differences in the agriculturally relevant bee faunas.

This study provides basic information about wild bees in Moroccan agroecosystems and shows the importance of these habitats in hosting a high diversity of wild bees. Identifying the bee species associated with each crop and their abundance in different regions could help to identify conservation strategies for wild bees in Moroccan agroecosystems. For example, Lasioglossum and Andrena were the most abundant and species rich bee genera recorded in this study. These ground nesting bees could be at risk from intensive agriculture because of deep tillage (Ullmann et al., 2016; Christmann, 2022), soil compaction by heavy machinery, accumulation of pesticides as well as the combined effects of these factors. Reducing machinery or adopting no tillage in farms could conserve ground nesting bees (Shuler et al., 2005). Assessment of pesticide impact and risk to solitary ground-nesting bees as demonstrated by Christmann (2022), may be particularly important for the protection of pollinators in developing countries with high pesticide use (Dicks et al., 2021). For faba bean, cherry and apple that bloom in early season, seeding adequate complementary flowering crops (Christmann and Aw-Hassan, 2012; Sentil et al., 2021), wildflowers strips (Garibaldi et al., 2011; Sutter et al., 2017; Muñoz et al., 2021) or hedgerows (Morandin and Kremen, 2013) may help to attract early flying bees like the genera Eucera, Anthophora, and Bombus and help to conserve these bees in the agricultural landscape.

This first comprehensive study on the bee fauna associated with crops in Morocco highlights the relatively few knowledge available and the urgent need for future research in this country including: (a) research on the value of pollination services (Blaauw and Isaacs, 2014; Christmann et al., 2017, 2021, 2022; Anougmar, 2021; Sabbahi, 2021), (b) applied research to enhance pollinator diversity, and (c) assessment of different global change drivers (i.e., climate change) and their influence on pollinator distribution and diversity (Scheper et al., 2013).