Preprint
Article

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

Altmetrics

Downloads

79

Views

36

Comments

0

A peer-reviewed article of this preprint also exists.

Submitted:

29 July 2024

Posted:

30 July 2024

You are already at the latest version

Alerts
Abstract
The environment significantly impacts the life of bees and their feeding. This study aimed to investigate bee foraging using melissopalynological analysis and DNA metabarcoding in ​​intensive farming, reserved, and urbanized areas. The highest alpha diversity was observed in the reserved and in the intensive farming areas. The urbanized area has less diversity. In the intensive farming area predominated: Sinapis, Helianthus, and Fagopyrum; in the reserved area: Melilotus, Helianthus, and Brassica. In the urbanized area garden plants: radish (Raphanus sativus) and cucumber (Cucumis hystrix) and agricultural plants: soybean (Glycine max) and melon (Cucumis melo) were often found. The most significant agreement was between the rbcL and the melissopalynological analysis. The ITS2 revealed equal matches with both rbcL and melissopalynology but this marker missed or underestimated some genera. Trifolium pretense and Brassica nigra, were identified simultaneously by three analyses. Species Convolvulus arvensis, Melilotus officinalis, Echium vulgare, Brassica rapa, Helianthus divaricatus, and Onobrychis viciifolia were found in all ecological areas. Imperfect databases put some limits for identification of some taxa using metabarcoding. Further research and expansion of plant databases is needed. Studying the food preferences of bees in different environmental conditions and landscapes is necessary to develop measures to preserve their population.
Keywords: 
Subject: Environmental and Earth Sciences  -   Ecology

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 cm3 of distilled water and left to swell for 2-3 hours (mixture 1). To 70 cm3 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 × 2500 10 ,
where N represents the number of pollen counted on the slide.
The percentage of occurrence was calculated using the formula:
Occurrence percentage = T o t a l   n u m b e r   o f   p o l l e n   g r a i n s   o f   a   p a r t i c u l a r   s p e c i e s T o t a l   n u m b e r   o f   p o l l e n   g r a i n s   o b s e r v e d × 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.

Supplementary Materials

Table 2. The most common plant taxa determined by melissopalynology. Table 3. Plant taxa at the species level identified using DNA metabarcoding with ITS2 and rbcL markers.

Author Contributions

Conceptualization, D.S., D.A. and S.P.; methodology, D.S., A.K. and T.P.; data curation, S.P.; analysis, D.S.; A.K.; K.A.; A.A.; resources, A.K.;, T.N.; writing—original draft preparation, D.S.; writing—review and editing, D.A., A.K. and S.P.; supervision, D.A; project administration, D.S. All authors have read and agreed to the published.

Funding

This research is funded by the Science Committee of the Ministry of Science and Higher Education of the Republic of Kazakhstan (Grant No. AP19676294).

Data Availability Statement

The data presented in this study are publicly available at the NCBI Sequence Read Archive (Bioproject ID: PRJNA1126018). Other data are contained within the article and Supplementary Materials.

Acknowledgments

Acknowledgments: We would like to thank the beekeepers of Kazakhstan for their assistance in completing the work. We express special gratitude to the beekeepers: Kuleshov Viktor Alexandrovich, Nuradil Gabit Nuradiluly, Bogdanova Elena Nikolaevna, Oleinikov Igor, Molzhigitov Yermek.

Conflicts of Interest

The authors declare that there are no conflicts of interest.

References

  1. Khalifa, A.M.; Elshafiey, E.H.; Shetaia, A.A.; El-Wahed, A.A.; Algethami, A.F.; Musharraf, S.G.; AlAjmi, M.F.; Zhao, C.; Masry, S.H.D.; Abdel-Daim, M.M.; et al. Overview of Bee Pollination and Its Economic Value for Crop Production. Insects 2021, 12, 688. [Google Scholar] [CrossRef] [PubMed]
  2. Belsky, J.; Joshi, N.K. Effects of Fungicide and Herbicide Chemical Exposure on Apis and Non-Apis Bees in Agricultural Landscape. Front. Environ. Sci. 2020, 8, 81. [Google Scholar] [CrossRef]
  3. Zewdie, A.; Amssalu, B.; Alemayehu, G.; Asaminew, T. Evaluating the Impact of Commonly Used Pesticides on Honeybees (Apis mellifera) in North Gonder of Amhara Region, Ethiopia. Journal of Toxicology 2023, 2634158, 13. [Google Scholar]
  4. Landaverde, R.; Rodriguez, M.T.; Parrella, J.A. Honey Production and Climate Change: Beekeepers’ Perceptions, Farm Adaptation Strategies, and Information Needs. Insects 2023, 14, 493. [Google Scholar] [CrossRef] [PubMed]
  5. Dolezal, A.G.; St. Clair, A.L.; Zhang, G.; Toth, A.L.; O’Neal, M.E. Native habitat mitigates feast–famine conditions faced by honey bees in an agricultural landscape. PNAS 2019, 116, 25147–25155. [Google Scholar] [CrossRef] [PubMed]
  6. Machado, T.; Viana, B.F.; da Silva, C.I.; Boscolo, D. How landscape composition affects pollen collection by stingless bees? Landscape Ecol. 2020, 35, 747–759. [Google Scholar] [CrossRef]
  7. Millard, J.; Outhwaite, C.L.; Kinnersley, R.; Freeman, R.; Gregory, R.D.; Adedoja, O.; Gavini, S.; Kioko, E.; Kuhlmann, M.; Ollerton, J.; et al. Global effects of land-use intensity on local pollinator biodiversity. Nat Commun. 2021, 12, 2902. [Google Scholar] [CrossRef] [PubMed]
  8. De Landa, F.G.; Alberoni, D.; Baffoni, L.; et al. The gut microbiome of solitary bees is mainly affected by pathogen assemblage and partially by land use. Environmental Microbiome 2023, 18, 38. [Google Scholar]
  9. Insolia, L.; Molinari, R.; Rogers, S.R.; et al. Honey bee colony loss linked to parasites, pesticides and extreme weather across the United States. Sci. Rep. 2022, 12, 20787. [Google Scholar] [CrossRef]
  10. Laha, R.C.; de Mandal, S.; Ralte, L.; et al. Meta-barcoding in combination with palynological inference is a potent diagnostic marker for honey floral composition. AMB Expr. 2017, 7, 132. [Google Scholar] [CrossRef]
  11. de Vere, N.; Jones, E.L.; Gilmore, T.; Moscrop, J.; Lowe, A.; Smith, D.; Herarty, M.J.; Creer, S.; Ford, C.R. Using DNA metabarcoding to investigate honey bee foraging reveals limited flower use despite high floral availability. Scientific RepoRts 2017, 7, 42838. [Google Scholar] [CrossRef] [PubMed]
  12. Millard, J.; Outhwaite, C.L.; Kinnersley, R.; Freeman, R.; Gregory, R.D.; Adedoja, O.; Gavini, S.; Kioko, E.; Kuhlmann, M.; Ollerton, J.; et al. Global effects of land-use intensity on local pollinator biodiversity. Nat Commun. 2021, 12, 2902. [Google Scholar] [CrossRef] [PubMed]
  13. Silliman, M.R.; Schürch, R.; Malone, S.; Taylor, S.V.; Couvillon, M.J. Row crop fields provide mid-summer forage for honey bees. Ecology and Evolution 2022, 12. [Google Scholar] [CrossRef] [PubMed]
  14. St Clair, A.L.; Zhang, G.; Dolezal, A.G.; O’Neal, M.E.; Toth, A.L. ; Diversified Farming in a Monoculture Landscape: Effects on Honey Bee Health and Wild Bee Communities. Environmental Entomology 2020, 49, 753–764. [Google Scholar] [CrossRef] [PubMed]
  15. Querejeta, M.; Marchal, L.; Pfeiffer, P.; Roncoroni, M.; Bretagnolle, V.; Gaba, S.; Boyer, S. ; Environmental variables and species traits as drivers of wild bee pollination in intensive agroecosystems—A metabarcoding approach. Environmental DNA. 2023, 00, 1–14. [Google Scholar] [CrossRef]
  16. Sponsler, B.D.; Matcham, E.G.; Lin, C.-H.; Lanterman, L.J.; Johnson, R.M. ; Spatial and taxonomic patterns of honey bee foraging: A choice test between urban and agricultural landscapes. Journal of Urban Ecology 2017, 1–7. [Google Scholar] [CrossRef]
  17. Couvillon, M.J.; Schürch, R.; Ratnieks, F.L.W. Dancing bees communicate a foraging preference for rural lands in high level agri-environment schemes. Curr. Biol. 2014, 24, 1212–1215. [Google Scholar] [CrossRef] [PubMed]
  18. McMinn-Sauder, H.; Lin, C.-H.; Eaton, T.; Johnson, R. ;. A Comparison of Springtime Pollen and Nectar Foraging in Honey Bees Kept in Urban and Agricultural Environments. Front. Sustain. Food Syst. 2022, 6, 825137. [Google Scholar] [CrossRef]
  19. Lucek, K.; Galli, A.; Gurten, S.; et al. Metabarcoding of honey to assess differences in plant-pollinator interactions between urban and non-urban sites. Apidologie 2019, 50, 317–329. [Google Scholar] [CrossRef]
  20. Tanaka, K.; Nozaki, A.; et al. Using pollen DNA metabarcoding to profle nectar sources of urban beekeeping in Kōtō-ku, Tokyo. BMC Res Notes 2020, 13, 515. [Google Scholar] [CrossRef]
  21. Rutschmann, B.; Kohl, P.L.; Steffan-Dewenter, I. Foraging distances, habitat preferences and seasonal colony performance of honeybees in Central European forest landscapes. Journal of Applied Ecology 2023, 60, 1056–1066. [Google Scholar] [CrossRef]
  22. Leponiemi, M.; Freitak, D.; Moreno-Torres, M.; et al. Honeybees’ foraging choices for nectar and pollen revealed by DNA metabarcoding. Sci Rep 2023, 13, 14753. [Google Scholar] [CrossRef] [PubMed]
  23. Park, B.; Nieh, J.C. Seasonal trends in honey bee pollen foraging revealed through DNA barcoding of bee-collected pollen. Insect. Soc. 2017, 64, 425–437. [Google Scholar] [CrossRef]
  24. Jones, L.; Lowe, A.; Ford, R.C.; Christie, L.; Creer, S.; de Vere, N. Temporal Patterns of Honeybee Foraging in a Diverse Floral Landscape Revealed Using Pollen DNA Metabarcoding of Honey. Integrative and Comparative Biology 2022, 62, 199–210. [Google Scholar] [CrossRef] [PubMed]
  25. Bell, K.L.; de Vere, N.; Keller, A.; Richardson, R.T.; Gous, A.; Burgess, K.S.; Brosi, B.J. Pollen DNA barcoding: Current applications and future prospects. Genome 2016, 59, 629–640. [Google Scholar] [CrossRef]
  26. Bänsch, S.; Tscharntke, T.; Wünschiers, R.; Netter, L.; Brenig, B.; Gabriel, D.; Westphal, C. Using ITS2 metabarcoding and microscopy to analyse shifts in pollen diets of honey bees and bumble bees along a mass-flowering crop gradient. Mol Ecol. 2020, 29, 5003–5018. [Google Scholar] [CrossRef]
  27. Ptaszyn ́ska, A.A.; Latoch, P.; Hurd, P.J.; Polaszek, A.; Michalska-Madej, J.; Grochowalski, Ł.; Strapagiel, D.; Gnat, S.; Załuski, D.; Gancarz, M.; et al. Amplicon Sequencing of Variable 16S rRNA from Bacteria and ITS2 Regions from Fungi and Plants, Reveals Honeybee Susceptibility to Diseases Results from Their Forage Availability under Anthropogenic Landscapes. Pathogens 2021, 10, 381. [Google Scholar] [CrossRef] [PubMed]
  28. Wirta, H.; Abrego, N.; Miller, K.; Roslin, T.; Vesterinen, E. DNA traces the origin of honey by identifying plants, bacteria and fungi. Scientific Reports 2021, 11, 4798. [Google Scholar] [CrossRef]
  29. Daugaliyeva, A.; Daugaliyeva, S.; Ashanin, A.; Beltramo, C.; Mamyrova, L.; Yessembekova, Z.; Peletto, S. Prokaryotic Diversity of Ruminal Content and Its Relationship with Methane Emissions in Cattle from Kazakhstan. Life 2022, 12, 1911. [Google Scholar] [CrossRef]
  30. Lowe, A.; Jones, L.; Witter, L.; Creer, S.; de Vere, N. Using DNA Metabarcoding to Identify Floral Visitation by Pollinators. Diversity 2022, 14, 236. [Google Scholar] [CrossRef]
  31. Khansaritoreh, E.; Salmaki, Y.; Ramezani, E.; Azirani, T.A.; Keller, A.; Neumann, K.; Alizadeh, K.; Zarre, S.; Beckh, G.; Behling, H. Employing DNA metabarcoding to determine the geographical origin of honey. Heliyon 2020, 6, e05596. [Google Scholar] [CrossRef] [PubMed]
  32. Chiara, B.; Cerutti, F.; Brusa, F.; Mogliotti, P.; Garrone, A.; Squadrone, S.; Acutis, P.L.; Peletto, S. Exploring the botanical composition of polyfloral and monofloral honeys through DNA metabarcoding. Food Control. 2021, 128. [Google Scholar] [CrossRef]
  33. Smart, M.D.; Cornman, R.S.; Iwanowicz, D.D.; McDermott-Kubeczko, M.; Pettis, J.S.; Spivak, M.S.; Otto, C.R.V. A Comparison of Honey Bee-Collected Pollen From Working Agricultural Lands Using Light Microscopy and ITS Metabarcoding. Environmental Entomology 2017, 46, 38–49. [Google Scholar] [CrossRef]
  34. Milla, L.; Schmidt-Lebuhn, A.; Bovill, J.; Encinas-Viso, F. Monitoring of honey bee floral resources with pollen DNA metabarcoding as a complementary tool to vegetation surveys. Ecological Solutions and Evidence 2022, 3, e12120. [Google Scholar] [CrossRef]
  35. Sickel, W.; Ankenbrand, M.J.; Grimmer, G.; Holzschuh, A.; Härtel, S.; Lanzen, J. ; Increased efficiency in identifying mixed pollen samples by meta - barcoding with a dual - indexing approach. BMC Ecol. 2015, 15, 1–9. [Google Scholar] [CrossRef]
  36. Murthy, K.M.; Khandayataray, P.; Ralte, L.; Laha, R.; Samal, D. Identification of plant species in multi flower honey by using Ribulose-Bisphosphate carboxylase gene (RBCL) coding region as barcode marker, Mizoram, Northeast India: An Indo: Burma hotspot region. Journal of Entomology and Zoology Studies. 2019, 7, 1475–1483. [Google Scholar]
  37. Bell, K.L.; Loeffler, V.M.; Brosi, B.J. An rbcL reference library to aid in the identification of plant species mixtures by DNA metabarcoding. Applications in Plant Sciences. 2017, 5, 1600110. [Google Scholar] [CrossRef]
  38. Saravanan, M.; Mohanapriya, G.; Laha, R.; Sathishkumar, R. DNA barcoding detects floral origin of Indian honey samples. Genome. 2019, 62, 341–348. [Google Scholar] [CrossRef]
  39. von Der Ohe, W.; Oddo, L.P.; Piana, M.L.; Morlot, M.; Martin, P. Harmonized methods of melissopalynology. Apidologie 2004, 35, S18–S25. [Google Scholar] [CrossRef]
  40. Karpovich, I.V.; Drebezgina, E.S.; Elovikova, E.A.; Legotkina, G.I.; Zubova, E.N.; Kuzjaev, R.Z.; Hismatullin, R.G. Atlas pyl’cevyh zeren/Pollen atlas/. Ekaterinburg Ural’skij rabochij 2015, 320. http://www.apiworld.ru/1453835538.html.
  41. Chen, S.; Yao, H.; Han, J.; Liu, C.; Song, J.; Shi, L.; et al. Validation of the ITS2 Region as a Novel DNA Barcode for Identifying Medicinal Plant Species. PLoS ONE 2010, 5, e8613. [Google Scholar] [CrossRef] [PubMed]
  42. https://github.com/apallavicini/bc4q; https://github.com/apallavicini/PLANiTS.
  43. Banchi, E.; Ametrano, C.G.; Greco, S.; et al. PLANiTS: a curated sequence reference dataset for plant ITS DNA metabarcoding. Database 2020, 2020, baz155. [Google Scholar] [CrossRef] [PubMed]
  44. de Melo Moura, C.C.; Setyaningsih, C.A.; Li, K.; et al. Biomonitoring via DNA metabarcoding and light microscopy of bee pollen in rainforest transformation landscapes of Sumatra. BMC Ecol Evo 2022, 22, 51. [Google Scholar] [CrossRef] [PubMed]
  45. Martin, S.G.; Hautier, L.; Mingeot, D.; Dubois, B. How reliable is metabarcoding for pollen identification? An evaluation of different taxonomic assignment strategies by cross-validation. PeerJ 2024, 12, e16567. [Google Scholar] [CrossRef] [PubMed]
  46. Noël, G.; Mestrez, A.; Lejeune, P.; Francis, F.; Miwa, M.; Uehara, K.; Nagase, A. Pollen meta-barcoding reveals foraging preferences of honeybees (Apis mellifera L.) along space-time gradient in Japan. bioRxiv, preprint. 2021. [Google Scholar]
Figure 1. Map showing different places of collection honey samples from ecological zones of Kazakhstan.
Figure 1. Map showing different places of collection honey samples from ecological zones of Kazakhstan.
Preprints 113598 g001
Figure 2. The plant composition of honey from different ecological areas was determined using melissopalynology. A – Representation of the main plant taxa by areas; B - intensive farming area; C - reserved area; D - urbanized area (in percent).
Figure 2. The plant composition of honey from different ecological areas was determined using melissopalynology. A – Representation of the main plant taxa by areas; B - intensive farming area; C - reserved area; D - urbanized area (in percent).
Preprints 113598 g002
Figure 3. The most common plant species in honey samples from different ecological areas. Microscope magnification x400.
Figure 3. The most common plant species in honey samples from different ecological areas. Microscope magnification x400.
Preprints 113598 g003
Figure 4. Relative abundance of the most common plant genera found in honey from different ecological areas. A – plant composition of honey with ITS2; B - plant composition of honey with rbcL. The abundance of taxa is indicated as a percentage of reads. The height of the bars represents the relative frequency of each plant genus in samples from each study region. Colors are used to differentiate each plant genus. The legend indicates the 17 most abundant plant genera.
Figure 4. Relative abundance of the most common plant genera found in honey from different ecological areas. A – plant composition of honey with ITS2; B - plant composition of honey with rbcL. The abundance of taxa is indicated as a percentage of reads. The height of the bars represents the relative frequency of each plant genus in samples from each study region. Colors are used to differentiate each plant genus. The legend indicates the 17 most abundant plant genera.
Preprints 113598 g004
Figure 5. Alpha Diversity (Total Number). Estimates of plant community diversity metrics detected in honey samples using OTU assignments obtained using both metabarcoding markers (ITS2 on the left, rbcl on the right). The colors indicate three types of ecological areas (intensive farming area, reserved area and urbanized area). Alpha diversity shows differentiation between all zone types (p > 0.05) for all methods used in this study. The results of the Kruskal-Wallis test are indicated in the upper corner of each graph.
Figure 5. Alpha Diversity (Total Number). Estimates of plant community diversity metrics detected in honey samples using OTU assignments obtained using both metabarcoding markers (ITS2 on the left, rbcl on the right). The colors indicate three types of ecological areas (intensive farming area, reserved area and urbanized area). Alpha diversity shows differentiation between all zone types (p > 0.05) for all methods used in this study. The results of the Kruskal-Wallis test are indicated in the upper corner of each graph.
Preprints 113598 g005
Figure 6. - Principal components analysis of read counts for amplicon sequence variants belonging to different ecological areas, calculated using two genetic markers: point = ITS2; triangle = rbcL. The color of the figure outline indicates honey samples from different ecological areas: green – Intensive farming area; blue – reserved area; pink – urbanized area. The axes represent the first and second principal components, the percentage deviation of which is indicated in parentheses.
Figure 6. - Principal components analysis of read counts for amplicon sequence variants belonging to different ecological areas, calculated using two genetic markers: point = ITS2; triangle = rbcL. The color of the figure outline indicates honey samples from different ecological areas: green – Intensive farming area; blue – reserved area; pink – urbanized area. The axes represent the first and second principal components, the percentage deviation of which is indicated in parentheses.
Preprints 113598 g006aPreprints 113598 g006b
Figure 7. Venn diagrams showing the number of matches by species of the most common plants from all honey samples. Legends indicate the plant species. The number of matches for one species is indicated in parentheses.
Figure 7. Venn diagrams showing the number of matches by species of the most common plants from all honey samples. Legends indicate the plant species. The number of matches for one species is indicated in parentheses.
Preprints 113598 g007
Table 1. Primers and amplification mode.
Table 1. Primers and amplification mode.
Primer name Nucleotide sequence Amplification mode and product length in nucleotide pairs Source
rbcLaf + adaptor TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGATGTCACCACAAACAGAGACTAAAGC Thermal cycling conditions were:
95 °C for 2 minutes;
95 °C for 30 seconds, 50 °C for 1 minute 30 seconds, 72 °C for 40 seconds (35 cycles);
72 °C for 5 minutes,
30 °C for 10 seconds.
~704 bp.
de Vere, N. et al. Using DNA metabarcoding to investigate honey bee foraging reveals limited flower use despite high floral availability. Sci. Rep. 7, 42838; doi: 10.1038/srep42838 (2017). [11]
rbcLr506 + adaptor GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGAGGGGACGACCATACTTGTTCA
ITS2 S2F TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGATGCGATACTTGGTGTGAAT Thermal cycling conditions were: 94℃ 5 min.
94℃ 30 sec, 56℃ 30 sec, 72℃ 45 sec, 40 cycles.72℃ 10 min
~350-500 (size range) in bp.
Chen S, Yao H, Han J, Liu C, Song J, Shi L, et al. (2010) Validation of the ITS2 Region as a Novel DNA Barcode for Identifying Medicinal Plant Species. PLoS ONE 5(1): e8613. https://doi.org/10.1371/journal.pone.0008613 [41]
ITS2 S3R GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGACGCTTCTCCAGACTACAAT
Note: the Illumina adapter area is highlighted in bold.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

Subscribe

© 2024 MDPI (Basel, Switzerland) unless otherwise stated