The experimental apiary, pollen load collection and species

The experimental apiary, which typically contains between 5 and 20 hives, is located in the Labex CEMEB platform ‘terrains d'expériences’ in Montpellier (France), a 4-ha area. The platform is located at the periphery of the city. The flora consists of common herbaceous species (e.g. Lolium perenne L. and Thymus vulgaris L.), trees (Quercus ilex L. and Arbutus unedo L.) and shrubs (Phyllierea angustifolia L. and Salvi rosmarinus L.) commonly found in French Mediterranean areas. In addition, the platform is surrounded by private gardens with ornamental species, which are potentially visited by the bees from our experimental hives. In June, in the Mediterranean area, various species are flowering simultaneously, and the weather is usually sunny and dry, which favours intense foraging activity by honeybee communities.

Forager honeybees with a visible pollen load were captured using a 50 ml Falcon® centrifuge tube on the flowers they were foraging on the experimental hive or closely around. Each Falcon tube was used for only one bee and for a single collection session. At least 10 forager bees have been captured for each flowering species. The bee-containing falcon was placed for 10 min in a refrigerator at 4°C. The cold anesthetized bee was then picked up from the falcon and the two pollen loads were carefully collected with forceps. Forceps were cleaned between the samples. Following pollen load collection, honeybees were released. Each pollen load was then placed off into a well of a 24-well plate and kept at room temperature in the dark. At the end of the collection, each pollen load was analysed (20 pollen loads per species or so), first using the pollen analyser and then using an ultra-fast UV/vis spectrometer (see the ‘Manual light-microscopy palynology’ and ‘Fluorescence spectroscopy on pollen load’ sections). A total of 14 pollen species were identified, and the analysis for each species was completed in under 4 h, for a total of 334 pollen loads. Plant species identification was carried out using the botanical expertise of the experimenter, with confirmation provided by the collaborative mobile application Pl@ntNet (https://plantnet.org/) (Neto-Bradley et al., 2025) and manual light-microscopy palynology. The 14 species collected and studied were: Cirsium arvense L. (Asteraceae), Crepis foetida L. (Asteraceae), Diplotaxis erucoides L. (Brassicaceae), Echium vulgare L. (Boraginaceae), Knautia arvensis L. (Caprifoliaceae), Oenothera speciosa L. (Onagraceae), Prunus domestica L. (Rosaceae), Punica granatum L. (Lythraceae), Q. ilex L. (Fagaceae), Rhus coriaria L. (Anacardiaceae), Rubus fruticosus L. (Rosaceae), Rubus ulmifolius L. (Rosaceae), Scabiosa atropurpurea L. (Caprifoliaceae) and Vitis vinifera L. (Vitaceae).

Manual light-microscopy palynology

Pollen identification was performed via light microscopy, enabling the correlation between specific plant taxa and fluorescence spectroscopy signals. Sixteen pollen loads from floral sources were randomly selected for analysis. Pollen grains were prepared through acetolysis following Erdtman's protocol (Erdtman, 1960). The acetolyzed pollen grains were mounted on glass slides using glycerol jelly as the mounting medium, sealed with coverslips and fixed with clear nail varnish to ensure long-term preservation. Microscopic observations were carried out at ×400 to ×1000 magnification using a bright-field microscope (Leica application suite v.4.13; Leica Microsystems (Switzerland) Ltd., Stereo and Macroscope Systems, Heerbrugg, Switzerland). Species identifications were achieved through morphological comparison with pollen reference databases by a trained palynologist. All analyses were conducted at the Institut des Sciences de l'Evolution de Montpellier (ISEM), University of Montpellier.

Fluorescence spectroscopy on pollen load

Fresh pollen loads were placed in a black 384-well flat-bottom plate (#3575; Corning, Glendale, AZ, USA), one load per well. Emission spectra were collected on a ClarioSTAR (BMG Labtech, Ortenberg, Germany) equipped with Linear Variable Filters for excitation and emission. The focus was optimised on a single well. As loads may vary in size, the optimal focus was calculated at several well positions and the optimal position was determined and applied to the whole plate. The gain was adjusted on the whole plate for each excitation wavelength. Emission spectra were then collected for six excitation wavelengths (365, 390, 395, 460, 530 and 625 nm), with a resolution of 1 nm and bandwidths of 10 nm. The emission spectra began 50 nm after the excitation wavelength. The raw data were treated with custom Python scripts.

The pollen analyser

The pollen analyser, which is dedicated to pollen identification, is a cylinder with internal dimensions of 34 mm in diameter and 21 mm in height, easy to use and portable. It is 3D printed with black (lower part) and grey (upper part) acrylonitrile butadiene styrene (ABS). The internal base of the lower part is covered with a small disposable white paper disc on which the pollen ball to be analysed is placed. The upper part of the box is sealed at the top with a 2-mm-thick glass pane measuring 2 cm by 2 cm, enabling a camera positioned 46 mm above the pane to capture clear images of the pollen balls inside. This camera is an 8-megapixel Picam NoIR v.2.1 (Raspberry Pi Holdings, Cambridge, UK) equipped with an 8-mm lens. The images are taken without an IR filter either in the camera or in the lens. Additionally, around the circumference of the upper part of the box, 12 evenly spaced holes accommodate LEDs positioned at a 45-degree angle towards the centre of the box's lower section. These LEDs are arranged into four families of three, with each group oriented 120 degrees apart to provide uniform illumination of the pollen load. The first family consists of RGB LEDs (T-1 3/4 (5 mm) full colour LED lamp, Kingbright, L-154A4SURKQBDZGW) illuminating at 645, 515 and 460 nm respectively. The remaining three groups consist of T-1 3/4, 5-mm UV LEDs operating at wavelengths of 365, 390 and 395 nm. Both the LEDs and the camera are controlled by a Raspberry Pi 4B, enabling precise synchronisation and automated operation.

For each test, a fresh pollen load is placed at the centre of the box on a white paper disc, and image acquisition is triggered via the Raspberry Pi by a custom Python code. First, the red LEDs are lit for 0.13 s while the camera images the pollen load in 0.03 s. Then the system is idle for 0.2 s. The process is repeated sequentially with the other wavelengths (green, blue, UV 365, 390 and 395 nm). At the end of the sequence, we obtain six images corresponding to the six wavelengths of light. For each set of six images corresponding to a given pollen load, the red channel of the image captured under red LED illumination is segmented to identify the position and contour of the load, thus defining the region of interest (ROI) using an optimal threshold determined by OpenCV's thresholding functions in Python (https://opencv.org/). Subsequently, for each excitation wavelength, the RGB channel intensities are averaged within the ROI. Therefore, for each load, 18 measurement points are collected. These measurement points are encoded in the excitation wavelength-pixel colour format. Pixel colours are abbreviated using the initial letter of each channel: green as ‘g’, blue as ‘b’ and red as ‘r’. For example, the intensity measurement in the red pixel after an excitation wavelength of 365 nm is denoted as 365r.

For each species, the 18 measurements from all analysed pollen loads are averaged and displayed on a radar chart with a logarithmic scale, generated using Excel software.

Principal component analysis (PCA)

To further assess the discriminatory power of spectral data on pollen loads, dimensionality reduction was performed using PCA. The analysis was conducted on either the entire dataset or a subset of samples comprising species from the same family or pollen loads with similar colours. In each case, the raw spectra collected at six wavelengths were organised into a single data matrix, where each row corresponded to a single pollen load and each column represented fluorescence intensity for a given excitation and emission wavelength. PCA was implemented in Python using the scikit-learn package (https://scikit-learn.org/). Before analysis, the data were normalised, and a PCA model retaining 10 principal components was constructed. Across all PCA models presented, the first four principal components consistently explained at least 85% of the total variance. The score plots for PC1 vs PC2 are displayed to visualise the separation of pollen loads.

Supervised classification analysis

Raw signal data acquired from both fluorescence spectroscopy and the pollen analyser were used as input for supervised classification algorithms. All classification experiments were conducted using the scikit-learn library (https://scikit-learn.org/). Before analysis, the data were normalised using the StandardScaler function. To identify the most effective classifier, we tested the following algorithms: Logistic Regression (LR), K-Nearest Neighbors, Naive Bayes, Random Forest and Support Vector Machine (SVM) (Cortes & Vapnik, 1995). For each approach, hyperparameter optimisation was performed using the GridSearchCV function. Given the class imbalance in the datasets, a stratified 10-fold cross-validation was applied at each step of the hyperparameter search to ensure that each fold preserved the class distribution. In a multi-class context, accuracy alone is insufficient for evaluating model performance. Therefore, we computed the following metrics for each class: (1) Recall = TP/(TP + FN), which represents the proportion of true positives (TP) relative to the total actual positives, that is, the sum of true positives and false negatives (FN). (2) Precision = TP/(TP + FP), which represents the proportion of true positives relative to the total predicted positives, that is, the sum of true positives and false positives (FP). (3) F1-score, defined as the weighted harmonic mean of precision and recall, calculated as: F1 = 2 × (Precision × Recall)/(Precision + Recall). It ranges between 0 and 1, where higher values indicate a better balance between precision and recall. (4) Receiver Operating Characteristic – Area Under the Curve – One vs Rest (ROC-AUC-OVR) evaluates a model's ability to correctly distinguish between multiple classes. It measures how well the model ranks positive instances higher than negative ones across all class comparisons, using a ‘one-vs-rest’ (OVR) approach. A higher ROC_AUC_OVR score indicates better class separation, meaning the model assigns higher probabilities to the correct classes with minimal overlap. This metric is critical in our multiclass classification problem, especially given the slight overlap between classes observed in the PCA analysis.

To further understand the contribution of individual features and interpret model decisions in our best classifiers, we leveraged the SHapley Additive exPlanations (SHAP) method (Lundberg & Lee, 2017). This approach quantifies the average marginal contribution of each feature across all possible feature combinations, providing insight into feature importance. Two complementary metrics were employed: the mean absolute SHAP values, which quantify the overall influence of each feature on model predictions, regardless of direction and the mean SHAP values, which capture the average contribution of each feature, taking into account the direction of the effect (positive or negative impact). The 20 most influential features were identified based on their mean SHAP values, highlighting those with the most consistent impact on predictions. Each feature is named by combining the excitation and emission wavelengths in the format excitation_emission. The six excitation wavelengths (365, 390, 395, 460, 530 and 625 nm) are labelled L1–L6, respectively. For example, ‘L1_746’ refers to the fluorescence intensity at 746 nm emission under 365 nm excitation.

Cross validation of species identity

From 7 to 15 June 2021, an average of 20 pollen loads of 14 species were collected manually from honeybees foraging on flowers of known species on several consecutive days, ensuring unambiguous species identification (see Fig. 1). Species identity was confirmed either through morphological analysis of the pollen grain using light microscopy or via the PlantNet application. The 14 species collected, which yield a set of pollen loads that are globally heterogeneous in size and colour, span a wide phylogenetic diversity, though seven of them belong to three families (three Rosaceae, two Asteraceae and two Caprifoliaceae).

Pollen load analysis by fluorescence spectroscopy

Pollen, like most plant organs, is auto-fluorescent upon excitation in a wide range of UV, visible and infrared wavelengths. Several studies have demonstrated that spectral signatures enable pollen identification, but current methods rely on fluorescence spectroscopy of soluble extracts, which are labour-intensive and destructive, especially when the pollen source is honeybee-collected pollen load. By contrast, we hypothesized that pollen spectral signatures could be produced directly from pollen load manually collected from the bees' legs, aiming for a non-destructive approach. On the day of collection, each pollen load was analysed by fluorescence spectroscopy using six excitation wavelengths: 365, 390, 395, 460, 530 and 625 nm. For each species, all spectral data were normalised and averaged (Figs S1–S6).

For each species, the shape of the spectra – that is the number and positions of peaks – is reproducible, regardless of the excitation wavelength. However, for some species such as D. erucoides, the relative fluorescence intensities vary significantly from one load to another (Fig. S7). The reasons for such variability remain unclear. We verified that the variability in fluorescence intensity was not caused by the orientation of the pollen load in the well, the position of the well on the plate or the delay between collection and fluorescence measurement. Similarly, no correlation between pollen load size and fluorescence intensity was found (data not shown). The possibility of pollen bleaching due to light exposure can be excluded, as the pollen loads were kept in the dark immediately after collection. Despite this variability, we observed that each species presents a complex fluorescence spectrum with at least one main peak of fluorescence intensity. For example, the spectrum of P. granatum after excitation at 390 nm shows a single peak near 500 nm, but associated with one or two secondary peaks of fluorescence intensity. The spectrum of K. arvensis after excitation at 365 nm shows a maximum near 425 nm, a clear secondary peak near 445 nm and a final broader peak near 595 nm. The superposition of the average fluorescent spectra of each of the 14 pollen load species for each excitation wavelength is shown in Fig. S8. These panels illustrate that the 14 species have specific emission profiles with distinctive spectral features, including different fluorescence peaks or different ratios of these peaks. For each species, the combination of data from these six spectra forms a spectral signature (Fig. 1a). Spectral signature differences can be visually clear in some cases but subtle in others, requiring a quantitative assessment of spectral specificity.

Principal component analysis

Given the complexity and large number of parameters in our dataset, clustering the data using PCA reduces this dimensionality by linearly transforming the data while capturing the most significant variations, making it easier to visualise clustering patterns and compare species data. The raw spectra of each pollen load, collected from the six excitation wavelengths, were concatenated into a single data matrix. The PCA model was elaborated with 10 principal components (PC), four of which explain 90.1% of the variance and whose loadings are presented in Fig. 1(c). These loading profiles reveal distinct patterns across excitation-emission pairs, meaning different spectral regions contribute in complementary ways to species-level discrimination. Plotting PC1 vs PC2 led to the separation of the pollen load spectral data into clusters corresponding to each pollen species (Fig. 1b). This supports that the six fluorescence spectra constitute a spectral signature for each pollen species, although on the 2D graph some of the clusters partially overlap.

As one limitation of pollen identification by chromatic assessment is to distinguish pollen loads of similar colours, it is interesting to evaluate how our method is able to discriminate such pollen loads. To determine the potential link between the spectral pattern of pollen loads and their colour, we first visually classified our collection of pollen loads according to tone. We obtained five groups: two species with white or very light pollen loads (S. atropurpurea and C. arvense), three species with yellow pollen loads (Q. ilex, D. erucoides and P. granatum), two species with orange pollen loads (R. fruticosus and C. foetida), three species with green pollen loads (R. ulmifolius, P. domestica and O. speciosa) and three species with black or very dark pollen loads (E. vulgare, V. vinifera and R. coriaria). We then superimposed the averaged spectra of the species in each of these tone categories (Fig. S9). Although spectra from the same taxonomic groups exhibit clear divergences that could potentially enable differentiation of pollen loads by species, we generally observe similar spectral patterns for pollen loads of the same colour. This similarity is particularly evident between P. domestica and O. speciosa, and even more pronounced between S. atropurpurea and C. arvense. However, the PCA analysis reveals that within each tone group, species have been distinctly clustered, and even a two principal components representation provides a fairly sufficient degree of separation to accurately identify each species (Fig. S9).

Similarly, the classification of the spectra according to their phylogenetic distances illustrates the resolving taxonomic power of the method. Among the 14 species analysed, the families Asteraceae and Caprifoliaceae are represented by two species, and the family Rosaceae has three members, two of which are of the same genus, Rubus (other families were represented by a single species). The spectra of each of these species are superimposed by family (Fig. 2). In the case of the Asteraceae and Caprifoliaceae, the two species exhibit distinct peaks, allowing pollen load classification by spectral comparison, as confirmed by well-separated PCA clusters in 2D (Fig. 2). In the Rosaceae family, P. domestica and R. ulmifolius have the most similar spectra; yet PCA analysis provides sufficient distinction to reliably differentiate each species, as shown in the two-component representation. Moreover, the two Rubus species display distinct spectral peaks, suggesting the method is able to achieve species-level taxonomic resolution (Fig. 2).

Supervised classification analysis

Having established by fluorescence spectroscopy and PCA analysis that pollen loads exhibit species-specific spectral signatures, we then developed an automated method for classifying pollen loads and reliably predicting their species identity. Since our dataset is already labelled with species identities, we trained several supervised learning models on a dataset derived from 70% of the spectral data obtained from pollen load collection. All classifiers, except Naive Bayes, achieve strong performances with accuracy exceeding 90%. These high levels of accuracy underscore the robustness of spectral signatures as reliable taxonomic markers, highlighting their potential for precise species identification. Among all models, the SVM classifier outperformed the others, achieving the highest classification accuracy – over 95% after hyperparameter optimisation – and the highest ROC_AUC_OVR score of 0.997, demonstrating its ability to effectively separate classes (Fig. 3a). The final optimised model used a linear kernel with C = 0.1 and gamma = ‘scale’. The resulting confusion matrix is presented Fig. 3(b). The vertical axis represents the true species identities, while the horizontal axis corresponds to the predicted species identities. The poorest classification performance was observed for D. erucoides, which was misclassified as O. speciosa and R. ulmifolius. Similarly, Prunus fructicosus and R. ulmifolius, belonging to the same genus, remain difficult to distinguish, as expected from their slightly overlapping clusters observed in the PCA analysis (Fig. 2a). However, their classification accuracy remains above 0.75. To further interpret the model's predictions, we employed the SHapley Additive Explanations (SHAP) method (Lundberg & Lee, 2017) to quantify and rank the wavelengths based on their influence on the model's predictions. A global analysis of the SHAP value distribution (Fig. 3c) indicates that no single feature exhibits a dominant colour, meaning that multiple features contribute broadly across different classes. However, some features still have a stronger positive contribution for specific classes – for example, L1_736, which contributes more to the classification of R. coriaria than to other classes. The highest SHAP values correspond to the excitation_emission wavelength pairs that have the greatest impact on the predictive model. However, these SHAP values remain small, at the 10−3 level, indicating that their overall impact on predictions is subtle. These observations suggest that our model, which relies on a large combination of features to predict species identity, is complex and robust enough to generalise well when incorporating new data. Sorting the 20 features with the highest global impact on the SVM predictive model (Figs 3d, 4a) reveals that the decisive wavelengths of the UV/visible emission spectra are clustered into three main regions: c. 450 nm (415–512 nm), 625 nm (597–658 nm) and 720 nm (699–739 nm). Comparing these 20 highest SHAP values across other classifiers with over 90% accuracy confirms this clustering pattern but with slightly different wavelength intervals, suggesting that other wavelengths across the entire UV/visible spectrum also play a significant role in species classification (Fig. 4a).

Development of a pollen analyser prototype

Based on the previous results, we have therefore developed a new methodology able to analyse and classify pollen loads using a system that is more affordable than a fluorescence spectrometer, easy to use and portable. The general principle of our pollen analyser is to collect the radiation emitted by a pollen load using a camera during a short pulse of exciting light produced by various monochromatic LEDs. Conceptually, our pollen analyser records a subsample of the fluorescence spectrum of a pollen load, but with a resolution of 75–100 nm corresponding to the bandwidth of each camera filter. As one of our primary goals is to create a low-cost system, we chose to leverage the Raspberry technology, combining a Pi computer with a v.2.1 NoIR camera that provides raw image data. The Raspberry v.2.1 NoIR camera uses a Sony IMX219 sensor whose characteristics meet our needs, in particular the spectral sensitivity characteristics covering the three spectral areas identified with SHAP values (Pagnutti et al., 2017; Fig. 4a). However, SHAP values analysis of SVM and LR models highlight the importance of emission wavelengths in the deep red. Consequently, the camera's infrared filter has been removed to allow the capture of infrared radiation. These SHAP values further demonstrate that, among all classifiers achieving accuracy rates over 90%, all the excitation wavelengths employed are critical; even the 465 nm (L1) is the most critical for species classification (data not shown). Consequently, the same six excitation wavelengths are retained for the prototype. The prototype of this analysis tool is 3D printed (Fig. 4b) and controlled via the Raspberry Pi 4, which turns the LED on and off and manages image acquisition for data collection. In practice, a pollen load deposited in the round chamber of our analyser is illuminated sequentially, for less than a second, by quasi-monochromatic light of 365, 390, 395, 460, 530 and 625 nm. Each time, an image is recorded. Each of the six images corresponding to the different wavelengths is segmented and analysed by calculating the average light intensity obtained for each of the three types of pixels within the load ROI: blue, green and red. In this way, each pollen load is characterised by 18 measurement values (Fig. 4c).

The same set of pollen loads from 14 species used for fluorescence spectroscopy analysis was also recorded on the day of collection with the pollen analyser. For each species, the average of all measurement points is calculated and the data are represented with a radar chart using a logarithmic scale (Fig. S10). Interestingly, as for the fluorescence spectroscopy approach, several species present variability for several measurement points. The reasons for these discrepancies are not clear yet. Still, we observe that the tendency to get more variable data is more pronounced when the pollen load colours darken. This phenomenon is maximal with E. vulgare whose pollen loads are the blackest of our samples (Fig. S11). The superposition of the 14 species illustrates that each radar chart forms a signature of the foraged plant (Fig. 5a).

Considering the substantial reduction in the number of recorded values compared to spectroscopic data (1381 features for fluorescence spectroscopy vs 18 for the pollen analyser), it was expected that the variations in signature shape between species were less pronounced. Nevertheless, the corresponding PCA model, which explains > 85% of the variance with four principal components whose main loadings are presented in Fig. 5(c), demonstrates that the data extracted from the pollen analyser are organised into clusters as shown by the score plots of PC1 vs PC2 (Fig. 5b). However, some clusters, such as that of P. granatum, are more diffuse compared to the fluorescence spectroscopy data, and the overlap between clusters is slightly more pronounced. Likewise, the radar charts for species with similar pollen load hues reveal consistent patterns, suggesting that colour information is embedded in the shape of the pollen signature (Fig. S12). Consequently, the PCA analysis of these species shows relatively broader clustering with frequent overlap, indicating that the pollen analyser is less effective at distinguishing pollen loads with similar colours (Fig. S12). By contrast, the superimposition of radar charts for species within the same family highlights the lack of correlation between related species and the shape of the radar chart (Fig. 6). For instance, R. fruticosus and R. ulmifolius within the genus Rubus exhibit distinct patterns, demonstrating the resolving power of our pollen analyser at the species level, which is consistent with results previously obtained using fluorescence spectroscopy analysis.

To enable the automatic identification of the plant source of pollen loads using the pollen analyser, we applied the same machine learning strategies as those used with fluorescence spectroscopy data. Once again, the SVM classifier outperformed the other models. While the ROC_AUC_OVR scores for Random Forest, SVM and LR are all closely aligned c. 0.99, SVM is the only model to exceed 90% accuracy, achieving 91.0% on the test data following hyperparameter optimisation (data not shown). This result closely approximates the accuracy obtained from the model trained on fluorescence spectroscopy data. Analysis of the corresponding confusion matrix reveals that the lowest precision occurs for V. vinifera, due to frequent misclassification of E. vulgare as Vitis (Fig. 7a). Although their radar charts display subtle differences (Fig. 7b), the discrimination between these species, characterised by darker pollen colours, is more challenging due to the high variability in the collected data. Another frequently misclassified species is C. foetida, often confused with R. coriaria and whose spectral signatures are highly similar (Fig. 7c).