Foraging of Honeybees from Different Ecological Areas Determined through Melissopalynological Analysis and DNA Metabarcoding

1. Introduction

Beekeeping provides humanity with one-third of the food consumed. By pollinating plants, bees provide food for humans and animals and preserve plant biodiversity on Earth. Pollination of crops by bees plays a vital role in plant biodiversity and is the basis for crop production [1]. In recent decades, a decline in bee populations has been observed worldwide. The reasons for this involve many factors: the use of pesticides [2,3], climate change [4], crop landscapes [5,6,7], pathogens and parasites of bees [8,9].

The preservation of the bee population depends on the available food supply and environmental conditions. Bee feeding can be monitored by the composition of honey obtained from plants. The botanical composition of honey depends on geographical location, climate, and agricultural methods [10]. To provide appropriate floral resources to bees, a better understanding of their food preferences is required. Knowing which plants honeybees feed on is critical for understanding the relationship between flower availability, honeybee nutrition, and health [11].

One of the most critical factors influencing bee feeding is the environment. The food supply of honeybees (Apis mellifera) and other pollinators depends entirely on the surrounding landscape’s floral resources (pollen and nectar). Intensification of land use in agricultural landscapes leads to a decrease in the diversity of native plants and a lack of available nectar and pollen, which is reflected in the interaction of plants and pollinators [12]. Monoculture production leads to the formation of “green deserts” with reduced biodiversity. In monoculture landscapes, bee colony growth occurs during flowering and is associated with increased flower availability [13]. However, this is followed by a period of food shortage, which leads to deterioration in the health of bees [5,14,15]. One way to solve this issue could be the preservation of islands of semi-natural or natural landscapes to support bee populations in monoculture fields [5,14].

Another human-modified landscape is the urban environment. According to some researchers, when compared between urban and rural landscapes, bees favor agricultural landscapes for foraging [16,17]. According to other authors, honeybees in agricultural and urban environments use the same resources to obtain nectar [18,19]. Urban or urbanized environments and suburban areas can become a food source for bees due to the landscaping of areas, parks, and gardens with native and ornamental plants. Urban beekeeping has recently attracted worldwide attention to greening as a beneficial method for conserving bees and promoting urban greening, an important measure to counter climate warming [20].

Unmanaged forests are considered the natural habitat of the western honeybee Apis mellifera and are believed to be significant sources of pollen and nectar. However, the supply of resources in managed forests is spatially and temporally limited. Rutschmann et al. note that in deciduous and coniferous forests, bees do not have enough food resources, and they more often forage for food in meadows and arable lands [21].

Bees’ foraging preferences vary depending on flower shape, color, scent, and nectar. It was reported that of the available flowering plants, honey bees used only 11% of plant genera to collect nectar or pollen [11,22]. According to several researchers, bees are good pollinators for native plants. Of the plants that bees use during the year, 54% are native, while ornamental plants serve them more as a source of pollen [23,24]. Horticultural plants are rarely abundant in honey samples, but they may be necessary to increase forage pollen’s nutritional diversity [23]. Bees may have biased visitation of specific plant genera because they provide relatively large amounts of nectar or pollen compared to other plants [19].

To study the feeding preferences of bees, the choice of methods for determining the floral composition of honey is essential. For many years, the primary method for assessing the floral composition of honey was the micromorphological analysis of pollen using light microscopy or melissopalynology. This method requires considerable time, labor costs, and exceptional knowledge in palynology [25]. In addition, interpretation of the results can be difficult because the pollen of some species is difficult to distinguish [26]. The development of next-generation sequencing (NGS) methods has made it possible to use them in various areas of scientific research, including studying the interaction of living organisms with the environment [27,28,29].

Currently, the DNA metabarcoding method has become an alternative for studying the biodiversity of the botanical composition in honey. Metabarcoding of plant loci allows the simultaneous sequencing and identification of multiple taxa and is an effective tool for describing ecological interactions between plants and pollinators [30]. Plant DNA sequencing can be used not only to determine the origin of honey over large geographic areas [31] but also to determine the adulteration of honey [32]. The method has higher sensitivity and resolution in identifying plant species than microscopic analysis [10], and may allow the identification of taxa that cannot be distinguished using light microscopy due to their low abundance in samples [33]. Despite obvious advantages, such as many barcode sequences allowing the identification of multiple species in a single reaction, the method does not accurately infer the relative abundance of pollen types due to possible quantitative errors encountered during DNA extraction: amplification and sequencing. Metabarcoding remains a high-cost method for several laboratories and requires appropriate conditions, equipment, and specialist training. In addition, one of the significant disadvantages of metabarcoding is the imperfection of databases for identification [33]. However, the reduced time and scalability make pollen metabarcoding a promising complementary tool for monitoring plant taxa [34], and laboratory costs are recouped by reduced human effort and reduced dependence on specialized taxonomic knowledge.

Choosing universal markers capable of distinguishing closely related plant taxa is particularly important when carrying out metabarcoding. The following genetic markers are most often used to determine the taxonomic composition of plants: plastid genes of ribulose bisphosphate carboxylase (rbcL), maturase K (matK), chloroplast marker trnL and ITS (mainly ITS2) [35]. Although the rbcL barcode is effective in identifying the botanical origin of honey with 99 to 100 percent confidence and has a high level of universality, in most studies, the use of the rbcL plastid region together with the nuclear ribosomal marker ITS2 strikes a balance between universality and discriminability [36]. Although ITS is thought to detect more species than rbcL, combining both markers provides more reliable evidence of geographic origin [15,31,34,37,38].

The purpose of this study was to compare bee foranging using melisopalynological analysis and double DNA metabarcoding in different ecological areas of Kazakhstan: an intensive farming, a reserved, and an urbanized. Studying bee foranging is necessary to understand the use of food resources and plant preferences in different ecological, climatic, and anthropogenic zones.

2. Materials and Methods

2.1. Sampling

Honey samples were selected from three ecological areas in Kazakhstan. These regions differed in geographical location, climatic conditions, environmental conditions, and landscape (Figure 1).

Intensive farming area - located in the northern part of the Karaganda region of Central Kazakhstan (49°48’“N 73°07’“ E). The climate of this area is sharply continental, with hot, moderate summers and cold winters with little snow. The main part of this area is occupied by grass-wormwood steppe on dark chestnut and chestnut soils. This is the main area of rain-fed agriculture and virgin land plowing. The nomadic apiary in this area was located near buckwheat, mustard, and sunflower fields. Periodically, after the flowering of the primary plants, the apiary wandered near the floodplains of the rivers.

The reserved area was located in the zone of the State Natural Mikhailovsky Nature Reserve, located in the North of Kazakhstan (53°29’ “N 61°39’“ E) in the forest-steppe zone of the Kostanay region. The reserve was created in 1967 to preserve the area’s hunting and commercial species of animals, such as forest birds. The territory’s climate is sharply continental, with harsh winters and hot summers. The reserve is located in a zone of mixed forests bordering steppes. The stationary apiary was located in the middle of the forest. In this area, the beekeeper planted plants around the apiary to feed the bees.

Urbanized area, the largest agglomeration of Kazakhstan - the zone of Almaty with its suburbs (43°15’ “N 76°54’“ E). This area is located in the extreme southeast of Kazakhstan at the foot of the Trans-Ili Alatau mountains, the northwestern part of the Tian Shan mountain range. Almaty’s climate is much milder due to relatively high temperatures in winter. Due to favorable climatic conditions, this region’s vegetation is rich in plant diversity. Grains, melons, and other cultivated plants grow at the foot of the mountains. Stationary apiaries from which honey was collected were located in summer cottages between the urban and suburban areas.

Forty-two honey samples were collected for analysis in autumn in 2023, including 15 samples from the intensive farming area, 12 from the reserved, and 15 from the urbanized area. All samples were collected sterilely in plastic containers and transported to the laboratory, where they were stored at 4° C.

2.2. Melissopalynological Analysis

The preparation of honey for melissopalynological analysis was carried out as described previously with some modifications [39]. Ten grams of honey were weighed and transferred into a 50 ml test tube. Twenty ml of distilled water (20–40 °C) were added to dissolve the honey, which was centrifuged at 4680 rpm for 15 minutes. The supernatant was removed. Twenty ml of distilled water was added to the sediment and centrifuged again for 5 minutes at 4680 rpm. The water was drained, the test tube was turned at an angle of 45° and dried on filter paper.

Glycerin gelatin was prepared as follows: 10 g of gelatin was poured into 60 cm³ of distilled water and left to swell for 2-3 hours (mixture 1). To 70 cm³ of glycerol, one μl of a solution of basic Ziehl fuchsin was added (produced by the Research Center for Pharmacotherapy (RCF), St. Petersburg, Russia) (Mixture 2). Mixture 2 was added to mixture 1, stirred, and heated in a water bath until a homogeneous mass was formed.

Preparation of a smear. The sediment was thoroughly mixed using a dispenser with a replaceable tip, transferred to a glass slide preheated to 40 °C, and evenly distributed over an area of 22 * 22 mm.

The glass with the sediment was heated at a temperature not exceeding 40 °C until the sediment was completely dried. A drop of glycerin gelatin was applied to a cover glass heated to a temperature of 40 °C and distributed crosswise along the diagonals. The cover glass was slowly lowered onto the dried sediment (to avoid the appearance of air bubbles). To ensure uniform glycerin gelatin distribution and optimal pollen swelling, the preparation was heated for 5 minutes at a temperature not exceeding 40 °C. After the glycerin gelatin hardened, the preparation was examined under a microscope.

Light microscopy. Pollen grains were viewed and identified under an Olympus BX43 light microscope (Olympus Corporation, Japan) with ADFImage Capture software.

Study of the pollen spectrum. Pollen grains were identified using the Pollen Grain Atlas [40]. Pollen counts and percentages of pollen types were calculated based on the frequency of pollen grains in different honey samples [39].

Pollen identification and counting were carried out in groups of 100 along five parallel, equidistant lines from one edge of the cover slip to the other. Square coverslip = 50 views per side. The cover glass area = L x L = 50 x 50 = 2500.

The following equation was used to obtain the total pollen count per slide:

Total pollen count per slide = $N\times {\displaystyle \frac{2500}{10}}$ ,

where N represents the number of pollen counted on the slide.

The percentage of occurrence was calculated using the formula:

Occurrence percentage = ${\displaystyle \frac{Total number of pollen grains of a particular species}{Total number of pollen grains observed}}\times 100$

2.3. DNA Extraction and Amplification

The preparation of pollen grains for DNA extraction was carried out as in previous studies [32]. Twelve and a half grams were taken from each honey sample in triplicate. Twenty-five ml of ultrapure water (UPW) was added and placed in a water bath at 50 °C for 15 minutes (until the honey dissolved) and centrifuged at 7000 g for 40 minutes. The supernatant was taken to the bottom mark. Pellets were collected from three test tubes and joined into one. Thirty ml of water were added to the sediment and centrifuged at 7000 g for 20 minutes. The supernatant was removed and the pellet was used for DNA extraction. Next, 200 μl of CD1 buffer from the DNeasy Plant Pro Kit (Qiagen, Hilden, Germany) was added to the pellet and vortexed for 10 minutes using ceramic beads on a VORTEX Genius 3 instrument at maximum speed. Forty µL of Proteinase K was added to each sample and incubated at 56 °C for 1 hour, then transferred to a new 1.5 μl tube and 200 µL of CD2 buffer were added. Further, all manipulations were carried out according to the protocol provided by the kit manufacturer.

To determine the botanical composition of honey from different ecological zones, we used double metabarcoding with ITS2 and rbcL markers. Amplification was carried out across two target regions, ITS2 and rbcL (Table 1). Poppy pollen was taken as a positive control. PCR-grade water was negative quality control for both PCR assays. Negative and positive controls were amplified and sequenced along with honey samples.

PCR amplification was performed in a Nexus gradient thermal cycler (Eppendorf AG, Hamburg, Germany). The total volume of the PCR mixture was 25 μl: GoTaq^® G2 Hot Start Colorless Master Mix, 2X (Promega Corporation, USA) – 12.5 µl; 1.2 μl of each primer (10 μM); 5 µl DNA (about 20 ng/µl) and 5.1 µl PCR pure water. The PCR product was checked in a 2% agarose gel and on an Agilent 2100 Bioanalyzer (Agilent Technologies, Waldbronn, Germany). No PCR product was observed in the extraction and amplification controls. Next, the PCR product was purified using the Agencourt AMPure kit (Beckman Coulter Inc. Beverly, Massachusetts, USA). The addition of Illumina Nextera indices was carried out in the PCR step. In the second stage of PCR, the reaction mixture consisted of 12.5 μl GoTaq^® G2 Hot Start Colorless Master Mix, 2X (Promega Corporation, USA), 2.5 μl index primers i5 and i7, 2.5 μl PCR product from the first step and five μl PCR water. Amplification program: one cycle at 95°C for three minutes, then eight cycles of amplification at 95°C for 30 seconds, one cycle at 55°C for 30 seconds, one cycle at 72°C for 30 seconds, and one cycle at 72°C for five minutes. The indexed PCR product was also purified using the Agencourt AMPure purification kit. The concentration and size of the PCR product were tested on an Agilent 2100 Bioanalyzer using the Agilent DNA 1000 Kit (Agilent Technologies, Waldbronn, Germany). The PCR product was quantified at each step using a Qubit 4.0 Fluorometer (Thermo Fisher Scientific) and the Qubit™ dsDNA HS Assay Kit (Life Technologies, Oregon, USA). The finished libraries were diluted to a 4 nM concentration and combined into a shared pool. The library was denatured to 10 pM concentration, and 15% PhiX was added. Sequencing was performed on a MiSeq Illumina instrument using the MiSeq Reagent Kit v3 - 600 cycles (Illumina corp., USA).

2.4. Bioinformatics and Data Analysis

Data obtained from amplicon sequencing were analyzed using the CLC Genomics Workbench software v. 24.0.1 (Qiagen). The raw FASTQ data underwent analysis using the “Data QC and OTUs Clustering” and “Estimate Alpha and Beta Diversity” workflow tools of the Microbial Genomics Module. Paired-end reads were joined and trimmed for low-quality score (Qscore < 0.05), nucleotide ambiguity (max of 2 nucleotides allowed), adapter sequences, and length. Duplicate sequences were merged and aligned against the PLANiTS database (Bianchi et al., 2020) and the reference database that Bell et al. (2017), developed for ITS2 and rbcL markers, respectively. Chimeric reads were removed, and taxonomy was assigned, creating OTU tables. OTU clustering was performed, setting a 97% similarity threshold for both databases and allowing the creation of new OTUs using a taxonomy similarity percentage of 80. The OTU sequences that were not taxonomically assigned based on the reference databases were submitted to Blast analyses versus the NCBI Nucleotide database. Taxa profiles of the negative and positive control were examined to assess procedure correctness and rule out cross-contamination, then they were not considered further. As rarefaction curves did not reach a plateau in all samples, comparisons were made at a sequencing depth of 20000 reads. OTU tables generated bar plots at the genus level for ITS2 and rbcL markers. Alpha diversity (diversity within groups using total number) and beta diversity (diversity between groups using the Bray-Curtis method and principal coordinate analysis, PcoA) were assessed, considering samples according to the ecological area where honey was collected (intensive farming, reserved and Urbanized areas). Statistical support for alpha diversity was determined using the Kruskal–Wallis test, while the PERMANOVA test was applied for beta diversity based on the Bray-Curtis matrix. OTU tables were used to perform a generalized linear model test of differential abundance based on ecological zones. The ‘Differential Abundance Analysis’ tool in CLC Genomic Workbench performs a TMM normalization to ensure sample compatibility by adjusting library sizes.

3. Results

3.1. Melissopalynological Analysis

A total of 42 samples from three ecological areas were studied by melissopalynological analysis (15 honey samples from the intensive agriculture zone, 12 samples from the reserved area, and 15 from the urbanized area). Light microscopy detected a total of 10203 pollen morphotypes in all samples. Four thousand fifty-five morphotypes were identified at the family level, 3364 at the genus level, and 2784 at the species level; 953 pollen were not identified.

According to melissopalynological analysis the predominant plants in the intensive farming area were Fagopyrum esculentum Moench (53%), followed by Helianthus annuus L. (17%), Urtica dioica L (7%). In the reserved area, the most numerous species was Helianthus annuus L. (24%), Fagopyrum esculentum Moench, Echium vulagare L. and Brassica napus L. (16%) each accounted for 17%. In the urbanized area, the predominant taxa were: Onobrychis arenaria (Kit.) DC. (19%), Brassica napus L. (17%), Melilotus albus Medik. (13%), Erysimum cheiranthoides L. (12%) (Figure 2 and Figure 3. Table 2. Suplementary file).

3.1. Metabarcoding

Forty-two honey samples were examined using metabarcoding using ITS2 and rbcL markers. One sample (A8) did not amplify for ITS2 for technical reason, and it was excluded from analysis.

After Illumina MiSeq sequencing and paired-end merging of the forward and reverse reads, we obtained an average of 399429 raw reads for ITS2 and 499559 raw reads for rbcL. The average fragment length was 239 bp for ITS2 and 346 bp for rbcL.

ITS2 metabarcoding revealed 233 taxa at the family level, 229 at the genus level, and 138 at the species level in all samples with the marker. rbcL metabarcoding identified a total of 326 plant taxa at the family level, 201 at the genus level, and 274 at the species level. Considering the three regions under study, there were 117 taxa at the genus level in the intensive farming area; 102 plant genera in the reserve area; and 173 in the urbanized area.

Figure 4 shows the taxonomic composition at the level of plant genera in honey samples from three different ecological areas of Kazakhstan.

DNA metabarcoding with ITS2 primers showed that, among the most frequently occurring plant genera, the following predominated in the intensive farming area: Sinapis (54%), Helianthus (14%), Picris (6%), Fagopyrum (5%), and Melilotus (4%). According to data with rbcL primers, the predominant species were Helianthus (19%), Fagopyrum (16%), Picris (12%), Linaria (11%), and Brassica (9%).

In the reserved area, the most common, according to ITS2 data, were Melilotus (26%), Helianthus (21%), Brassica (17%), Berteroa (6%), Filipendula (5%), Onobrychis (5%). According to rbcL data, the dominant plant genera in the reserved area were Helianthus (22%), Melilotus (19%), Brassica (15%), Solidago (12%), Erodium (7%), Lonicera (6%).

In the urbanized area, according to data obtained with the ITS2 marker, the most common genus were Raphanus (17%), Helianthus (16%), Reseda (12%), Trifolium (9%), Melilotus (9%). The rbcL marker showed that in the urbanized area, the dominant plant taxa were Cucumis (25%), Melilotus (11%), Helianthus (9%), Gossypium (9%), Peganum (8%), Malva (6%).

The most frequently occurring plant species for both markers are shown in Table 3 (Supplementary file 2). As can be seen from the table, the following plant species predominated in the intensive farming area: Sinapis alba (72%), Helianthus divaricatus (20%), and Fagopyrum esculentum (16%). In the reserved area, the most common were: Helianthus divaricatus (24%), Filipendula vulgaris (20%), Melilotus albus (20%), Echium vulgare (17%), Melilotus officinalis (16%). The dominant plant species in the urbanized area were Raphanus sativus (31%), Trifolium pratense (16%), and Cucumis melo (15%).

Despite the fact that the number of species-level taxa was considerable, especially with a primer to the rbcL gene, we took taxa at the genus level for analysis. This is because existing databases may not contain reference sequences of native plant species.

3.1. Differences between Foraging Preferences of Bees in Different Ecological Zones: Alpha and Beta Diversity

High alpha diversity was observed in the reserved and in the intensive farming areas. In the urbanized area, diversity based on the total number of OTUs is lower, as confirmed by both markers (Figure 5).

As can be seen from Figure 5, alpha diversity was significantly lower in the urbanized area than in the others (p value = 0.00001). The rbcL marker showed greater species diversity within groups.

To assess beta diversity between honey samples from different zones, principal component analysis (PCoA) based on OTU-level Bray-Curtis distances was used (Figure 6). PCA analysis confirmed significant differences between the three studied ecological zones (Bonferroni p ≤ 0.05).

3.1. Comparison of Metabarcoding and Melissopalynology Data

A comparison of agreements in identifying dominant plant species across the three analyses is presented in a Venn diagram (Figure 7).

As can be seen from Figure 7, of the most common plant species, the greatest agreement is observed between the rbcL marker and palynological analysis. While the ITS2 marker revealed almost the same number of matches with both rbcL and palynology. And only two species, Trifolium pretense and Brassica nigra, were identified simultaneously by three analyses. Species such as Convolvulus arvensis, Melilotus officinalis, Echium vulgare, Brassica rapa, Helianthus divaricatus, Onobrychis viciifolia were not only identified by different methods, but also occurred in different ecological areas.

4. Discussion

This study determined the foraging features of bees from different ecological areas. To more thoroughly determine the botanical composition of honey, along with traditional melissopalynological analysis, plant DNA metabarcoding was used at two loci: ITS2 and rbcL. Previous studies reported that the ITS2 marker was superior to rbcL in identifying plant taxa, and errors in identification at the genus level were low [31]. Even with identity rates reaching 100%, Bell et al. obtained fewer high-quality sequence reads from rbcL than ITS2. The authors attribute this to the fact that the rbcL fragment is longer and has less overlap between forward and reverse reads, resulting in more reads being removed during filtering [37]. At the same time, Carneiro et al., using a ITS2 dataset, found fewer plant families in comparation with the rbcL [44]. To obtain better results when assessing the biodiversity of the honey’s plant composition, several authors recommend using multple loci [31,36,44].

In our study, the higher number of plant taxa in honey was obtained using the rbcL marker (326 taxa), while 233 taxa were identified using ITS2, and 219 taxa were identified using melissopalynological analysis. When comparing the three analyses considering the most represented plants, most matches were between the rbcL marker and the melissopalynological analysis. The ITS2 marker revealed almost the same number of matches with rbcL and melissopalynology, but this marker missed or underestimated some relevant genera (e.g., Sinapis). Only two species, Trifolium pretense and Brassica nigra, were identified by all three analyses. Species such as Convolvulus arvensis, Melilotus officinalis, Echium vulgare, Brassica rapa, Helianthus divaricatus, and Onobrychis viciifolia were not only identified by different methods but also occurred in different ecological areas.

Despite several advantages of NGS, we also used the traditional method of melissopalynology to determine the plant composition of honey. Studies by Laha et al. showed that palynological data revealed many plant species that could not be identified using NGS. The authors attribute this to the incompleteness of the plant databases [10]. Smart et al. noted that although DNA sequencing analysis identified a more significant number of taxa, this method is limited by the availability of relevant sequences for comparison in reference databases. Some taxa were incorrect and indicated gaps in the database within the phylogenetic lineage [33]. Aligning local species sequences with a local reference database instead of a global database increased accuracy for ITS2 from 63% to 75%, and for rbcL from 46% to 58% [45]. The results of our research indicate that plant taxa databases do not always contain information about local plants. Therefore, reliability was not always correct at the species level, especially with the ITS2 primers. More plant taxa were detected using the rbcL marker, especially at the species level.

Melissopalynological data more accurately determined plant taxa at the species level, especially in the reserved area and the intensive farming area. To determine plant taxa, the Atlas of pollen grains was used, based on plants growing in a natural area similar in climatic and geographical conditions to the regions of Northern and Central Kazakhstan. However, DNA metabarcoding of honey made it possible to discover more plant taxa, especially in the urbanized area, where there is an abundance of plants growing in the country’s southern regions, which are not included in the Pollen Plant Atlas. In the urbanized area, some specific plants were more accurately identified using metabarcoding: Peganum harmala, Cucumis melo subsp. Melo, Glycine max, Nitraria sphaerocarpa, Saussurea elegans.

Landscape class (rural, suburban, and urban) explains spatial variation in the composition of plants fed by honey bees but not in taxa richness [46]. In our case, the highest alpha diversity was observed in the reserved and in the intensive farming areas. The urbanized area had less diversity based on the total number of OTUs than the others (p value = 0.00001).

The predominant genera in the intensive farming area were Sinapis (54%), Helianthus (14%), Linaria (11%), Brassica (9%), Picris (6%), Fagopyrum (5%), Melilotus (4%). These data were concordant for both markers, with slight deviations in percentage, and with the data of melissopalynological analysis. In this area, apiaries are located near fields with monoculture plants: buckwheat, mustard, and sunflower. Previous studies reported that bees prefer these plants during the flowering in monoculture zones. However, after this period, there is a lack of food, which leads to a deterioration in the nutrition of bees, their productivity, and health [5,27]. In such cases, an essential factor is the correct management of the apiary, aimed at finding available food resources. In our case, apiaries migrated from monoculture fields to river floodplains. Therefore, in addition to monoculture plants, an abundance of other plants was observed in the intensive farming area.

In the reserved area, the stationary apiary was surrounded by deciduous and coniferous forests. It has been observed that in forested landscapes with high forest cover, bees are limited in pollen, and they more often forage for food in meadows and arable lands [21]. In this area, the beekeeper planted plants around the apiary, including those brought from other places to feed the bees. The research results show that the bees used the available plant species around the apiary. In the reserved area, the following plant taxa prevailed: Melilotus (26%), Helianthus (22%), Brassica (17%), Solidago (12%), Erodium (7%), Lonicera (6%), Berteroa (6%), Filipendula (5%), Onobrychis (5%). The biodiversity of plants for feeding bees was greatest in the reserved area. This may be due to the fact that the beekeeper’s management allowed the bees to use additionally available resources near the apiary.

According to metabarcoding data in urbanized area the following plant genera are predominated: Cucumis (25%), Raphanus (17%), Helianthus (16%), Reseda (12%), Trifolium (9%), Melilotus (9%), Gossypium (9%), Peganum (8%), Malva (6%). At the same time, a large share was accounted for by garden plants, such as radish (Raphanus sativus), cucumber (Cucumis hystrix), as well as agricultural plants, such as soybean (Glycine max), melon (Cucumis melo).

5. Conclusions

The study of bee foraging in different ecological areas is necessary to understand the feeding preferences of bees and the development of measures to preserve the bee population. Our research shows that feeding bees depends on the availability of plants, the surrounding landscape, and climatic and geographical conditions. Human management of apiaries has an essential influence on the feeding of bees. This made it possible to diversify the monotonous diet of bees in the intensive farming area and provide food for bees in forest in reserved area. An urbanized environment additionally provides crops for bee nutrition and can serve as one of the ways to preserve bee populations.

To determine the botanical composition of honey, it is advisable to couple melissopalynology and metabarcoding. The melissapalynological method was more accurate in identifying plant taxa, especially at the species level, but it relies on the availability of local reference Atlas. However, honey DNA metabarcoding has revealed a higher number of plant taxa. Specifically, more plant taxa were detected using the rbcL marker than ITS2, especially at the species level. The imperfection of databases, particularly for native plants, makes accurate identification using metabarcoding sometimes difficult. Therefore, further research and expansion of the plant database is necessary.

The study of bee feeding using DNA metabarcoding in three ecological areas was carried out for the first time in Kazakhstan. Our findings may have important implications for the development of beekeeping in the study area.