Honey bee food resources under threat from climate change

Domain: chemistry
Published: 2025-12-30
Editorial outcome: minor_revision

Abstract
Plant-pollinator interactions are essential for plant productivity but face growing threats from climate change, including vegetation loss and mismatches in flowering. Yet, the consequences for bee food resources remain poorly understood at continental scales. Here, we analyse 2500 samples collected by honey bees (Apis mellifera) between May and August 2023 from 310 locations across Europe using ITS2 metabarcoding. We derive climatic response curves of floral resources and assess exceedance risks of interaction loss under projected climate scenarios. Our findings reveal that rising temperatures and reduced precipitation decrease the diversity of foraging resources across Europe, pushing many plants beyond critical limits. When both warming and drying coincide, the potential for resilience through temporal or spatial buffering is strongly constrained. These declines pose serious risks to bee nutrition, ecosystem functioning, and food security. Our study underscores the urgency of mitigating climate change to preserve vital plant-pollinator systems and the services they sustain.

Introduction
Plant-pollinator relationships are fundamental biotic interactions in ecology due to their key role in plant productivity, food webs, and human well-being. Animal pollinators, principally insects, contribute an estimated $235–$577 billion annually to the global economy by providing pollination services to an estimated 87.5% of the world’s flowering plants and 75% of major global food crops. However, without suitable and abundant plants supplying pollen and nectar, many pollinator populations will decline along with their vital services.

An increase in temperature has already been reported to adversely impact vegetation composition and is known to cause flowering and spatial shifts, which may lead to interaction mismatches between plants and pollinators. The decrease in floral diversity is a key factor driving bee diversity loss, and evidence correlates bee pollinators and wild plants as being tightly linked and in parallel decline. Several studies project that climate change will affect the survival of many plant species. Depending on the climate change scenario, 3-21% of all endemic plant species in Europe could face extinction by 2050. Under the most severe scenario, 22% of total plant species would become critically endangered and 2% extinct by 2080. However, the impact of climate change varies by region, with plant diversity projected to increase in Northern Europe while declining in the Mediterranean region. Furthermore, by the end of the 21st century, many plant species are expected to disappear from over half of their current geographic ranges in Europe, with distribution shifts northwards and uphill or suffer local or global extinction. These changes will severely affect the availability of food resources to pollinators. A decline in floral diversity may limit essential nutrients for pollinators, affecting their overall health and reproductive success. Reduced access to diverse pollen sources can lead to nutritional deficiencies, weakened immune systems, and increased vulnerability to pathogens and parasites, as well as other environmental stressors. Additionally, shifts in floral availability may alter foraging behaviour, forcing pollinators to expand their foraging ranges. These disruptions may ultimately contribute to population declines, reduce pollination efficiency, and negatively impact both wild plant communities and agricultural yields for pollinator-dependent crops.

Our study investigates the impacts of climate change on the food resources of the honey bee Apis mellifera, the dominant global pollinator with the broadest geographical and temporal range in Europe. Although A. mellifera is not in decline, we select it due to its ubiquity, generalised foraging, and critical role as a pollinator in agroecosystems. This species’ extensive spatial and temporal activity makes it a useful system for evaluating large-scale patterns in floral availability under climate change, especially in managed landscapes. The major goal of our study was to quantify the expected impact of climate change on honey bee foraging across different regions of Europe and assess the extent to which resilience may be achieved through spatial or temporal shifts in floral availability. To achieve this, we conduct a systematic survey of honey bee-collected pollen at unprecedented spatial and temporal scales. We analyse 2500 pollen samples collected by honey bees biweekly, for nine sampling rounds, between May and August 2023, from 310 apiaries distributed across all 27 EU countries. We identify the pollen composition of these samples using DNA metabarcoding of the internal transcribed spacer 2 (ITS2). We structure our analysis around five core objectives: (1) assess spatial and temporal patterns of pollen diversity, composition and crop contribution across Europe in relation to climatic gradients, to establish a baseline for evaluating potential impacts of climate change on floral resource availability; (2) quantify historical variation in temperature and precipitation at the local study sites to contextualize and scale future projections; (3) determine the realised climatic niches of foraged pollen plant species by estimating empirical climatic niche response curves to assess their sensitivity to climate change; (4) quantify potential floral resource loss under different climate projection scenarios; and (5) evaluate the resistance and resilience of plant-pollinator interactions by spatial or seasonal compensation through alternative floral resources. Together, these objectives provide an integrated framework for understanding how climate-driven environmental change may affect floral resource dynamics and, consequently, honey bee foraging potential across Europe.

Results
Pollendiversity,composition,andcropcontribution
As a baseline for understanding how floral resource availability may shift under future climate scenarios, we first assessed how honey bee-collected pollen diversity and composition varied across Europe in relation to key climatic gradients. We found that temperature and precipitation gradients, as well as latitude, had significant effects on pollen Effective Species Richness (eH), with a mid-latitude peak particularly evident in late summer. However, these effects had limited explanatory power (quadratic linear model with latitude and month as predictors; eH log-transformed: F(3,2493)=90.27, p<0.001, adjusted R²=0.10; Fig. 1A, B, Supplementary Fig. 1). Seasonal models confirmed that this mid-latitude peak was most pronounced in July and August.

The proportion of pollen from crops per sample was generally low, with an overall average of 19% (20% standard deviation, Fig. 1C), but with peaks during mass-flowering events when crops contributed up to a mean of 50%. These peaks were primarily driven by taxa such as Brassica spp., Helianthus annuus, and Zea mays identified among the 25 most abundant pollen sources (Supplementary Fig. 2). Even during such episodic crop dependencies, honey bees maintained a comparably diverse diet at all sites throughout the season, underscoring the importance of wild floral resources (Fig. 1A).

Pollen composition varied markedly across Europe and over the sampling season, reflecting shifts in the identity of plant species foraged upon (visualised via NMDS; Supplementary Figs. 3 and 4). To quantify how this compositional turnover related to environmental gradients, we used Generalised Dissimilarity Modelling (GDM), which explained 48.47% of the total variance in pollen composition (pseudo-F(3, 48 205)=0.52, p<0.001; Supplementary Fig. 1), with contributions of 27.53% from temperature (ΣI-spline=1.82), 14.14% from geography (ΣI-spline=0.94), and 6.79% from precipitation (ΣI-spline=0.44). These results indicate that floral resource composition is strongly structured by climatic and spatial factors.

Impact of regional warming
To estimate the potential vulnerability of current floral resources under future warming, we next examined regional temperature trends and their effects on the climatic suitability of plant species foraged by honey bees at observed sites. Warming and changes in precipitation can lead to several impacts on bee resources: firstly, they can alter the local plant composition; secondly, they can disrupt plant phenology, affecting key events such as flowering timing; and thirdly, they can impact the quantity and quality of resources. These effects can lead to mismatches in plant-pollinator interactions, changing pollen’s temporal and spatial availability. We were thus interested in how future scenarios might pose risks for the currently available pollen resource under climate change and continuing historical trends. Global warming is projected to reach 2°C by 2040 to 2060, even under aggressive mitigation scenarios. However, such global trajectories can mask substantial regional variation, which is critical for forecasting. We therefore included an analysis of historical warming trends at our actual sampling sites (Supplementary Fig. 5A and 5B) to ground our projections in empirically observed, region-specific climate change and contextualise the effect size of projections. Based on the results of these analyses, sites in Southern Europe have, on average, warmed by 1.8°C over the past 50 years, Central Europe by 1.6°C, and Northern Europe by 0.2°C (Supplementary Fig. 5A and 5B). A previous study with regional forecasts predicts a 2°C increase in Northern Europe by 2040 and Central and Southern Europe by 2030, ahead of global averages.

To study the impacts of rising temperatures and decreasing precipitation on bee pollen resources, we first estimated temperature and precipitation response curves for each plant taxon in the dataset based on their sampling time, location, and corresponding climatic values (for more details, see Methods). We then incrementally adjusted the mean temperature (0 to +5°C) and precipitation (0 to −50 mm) for each location and time. Using these updated climatic conditions, we assessed whether conditions exceeded the percentile thresholds of each taxon’s empirical response curve and thereby whether the conditions remained suitable or became unsuitable. Our results revealed a severe impact of temperature increases on current bee food resources (both wild plants and crops) across Europe, with a gradual reduction in impact from Southern to Northern Europe (Fig. 2; Linear model F(5, 1189)=100.3, p<0.001, adjusted R²=0.30; Supplementary Table 1). Southern European wild plants are most affected, with Mediterranean islands such as Cyprus, Crete, and Malta facing significant risks even with just 1 °C of warming, in which about half of the current bee resources become at risk (Fig. 2, Supplementary Fig. 6A). This trend worsens for higher temperatures. In comparison, Northern Europe is less affected up to a 2 °C rise, although with a notable increase in risk under more severe warming scenarios.

These findings align with Bakkenes et al., who noted that Southern European plant species are more sensitive to warming than Northern ones. Seasonal variations revealed, in general, a higher risk during summer months compared to spring (Fig. 3A, Supplementary Fig. 7A). However, this effect did not persist at 5°C warming, when a high proportion of plants became at risk throughout the season. Temperatures warmer than their critical limits negatively influence the growth and reproduction of species, contributing to biodiversity loss. Projections suggest an increased risk for European species and the possible extinction of 3-21% of endemic plants by 2050. Our empirical data on temperature niches support this alarming notion of the sensitivity of bee pollen resources, with even higher proportions at risk on local scales.

Impact of changing regional precipitation patterns
In parallel, we assessed how changes in precipitation, both observed and projected, might alter the climatic suitability of pollen-providing plants at observed sites, again focusing on species-level risk under plausible climate futures. Historically, precipitation patterns have also changed, and it is projected that precipitation will severely change across most of Europe with climate change (Supplementary Fig. 5C and 5D). For example, Vautard et al. predict a 10% reduction in annual mean precipitation for Southern Europe and an over 10% increase for Northern Europe for the 2 °C warming scenario. However, the effects of precipitation on bee pollen resources are more complex than those of temperature, and annual mean values are not a good proxy for local water availability for plants: the increased intensity and frequency of extreme precipitation events influence annual means to a large degree, but the water will not necessarily be available for plants over extended periods, and direct harm can be caused by flooding. Therefore, we focused on the impacts of a gradual decrease in mean precipitation at a local scale, excluding short-term extreme rainfall events. Based on species-specific climatic response curves, reduced precipitation markedly increased risks for local plants across all latitudes, with the most severe impacts in Southern Europe but also notable effects in Northern and Baltic regions (Fig. 2; linear model F(5, 1138)=56.9, p<0.001, adjusted R² = 0.20; Supplementary Table 1). Northern countries were particularly sensitive to spring rainfall losses, while Southern countries were more affected by summer rainfall deficits (Fig. 3B, Supplementary Fig 7B). This pattern aligns with the Copernicus Climate Change Service’s 2023 report, which predicts increased summer precipitation in Northern Europe and decreased precipitation in Central and Southern Europe for the 2 °C warming scenario. Our results also indicate that Southern Europe, especially the Algarve region (latitude 37°) in Portugal, is already experiencing drought stress with a visible impact on pollen resources, whereas the impact at latitude 46° was less pronounced compared to the rest of Europe (Fig. 2, Supplementary Fig. 6B).

Resistance and resilience
Climate change impacts on a local scale can, however, alter plant species distributions and phenology in various ways, so a high risk for local populations may not necessarily reduce overall plant resources. Affected species may, for example, be replaced by other taxa that can provide suitable resources for honey bees by immigration or changes in flowering phenology. We thus also estimated the proportion of at-risk plants for precipitation and temperature in a combined analysis that could potentially be replaced by species better adapted to the new climatic conditions available in the wider region and/or other seasons. This assessment of resistance (non-risked species) and resilience (availability of replacement species) indicated that small changes in temperature and precipitation might be mitigated by the broader regional species pool (Fig. 4, Supplementary Fig. 8), but the combined effects of temperature increase and precipitation reduction limit this potential. Severe changes, particularly with this combination, create conditions that show low resistance and resilience, i.e. local populations are highly affected and are difficult to fill with replacement taxa. This is likely to lead to significant losses in plant diversity and honey bee food resources (Fig. 4).

Our results are limited to the pollen collected, which may not represent the full resources spectrum at the study sites. Consequently, our findings might underestimate overall resilience. Conversely, unpreferred plants included in the analyses contribute equally to resilience calculations as preferred ones, despite potentially lower nutritional quality or abundance. This could lead to an overestimation of resilience. Therefore, uncertainties in both directions must be considered when interpreting our results. Despite these limitations, we can confidently state that climate change will disrupt local pollination networks particularly when temperature increases are coupled with precipitation losses, with severe implications for plants, honey bee nutrition, and competition with wild pollinators. The limited buffering capacity under combined stress conditions suggests that for many sites, resource losses may not be easily compensated, posing a challenge for pollinator persistence and network stability.

Methods
Pollensampling
Pollen samples were collected simultaneously every two weeks between May 1st and August 27th 2023, from 310 stationary apiaries distributed across all 27 European Union (EU) countries. These pollen samples were voluntarily collected and provided by participating beekeepers from their own managed honey bee colonies, solely for scientific research and under written consent. Sampling was organised through national coordinators in each EU member state, who served as contact points between the INSIGNIA-EU consortium [https://www.insignia-bee.eu] and participating beekeeper citizen scientists and ensured adherence to standardised sampling procedures. Sampling was performed by activating standard pollen traps for approximately 24 hours on two hives per apiary during routine beekeeping practice. As domesticated and managed Apis mellifera colonies and their hive products are not subject to wildlife collection regulations or national access restrictions within the EU, no governmental sampling permits were required. No animals were harmed during pollen collection. All publicly shared data were anonymised in accordance with GDPR and with the agreement of the contributing beekeepers. This included removal of beekeeper identifiers, addresses, and any other personal information, as well as truncation of geographical coordinates in public metadata. Anonymised consent documentation for participating citizen scientist beekeepers is available from the authors upon request.

About 5g of pollen per apiary (an equal homogeneous amount from each hive) was gathered and stored in capped vials together with 12g of silica gel, following Quaresma et al. Of the 2790 expected pollen samples, only 2500 were analysed. The shortfall of 290 samples was due either to PCR amplification failure, logistic issues, or the absence of collected pollen.

DNA metabarcoding
Samples were homogenised in a magnetic stirrer by adding 4ml of sterile ultrapure water to 2g of silica-dried pollen. A volume of 200µl of the homogenate solution (~50mg of pollen per sample) was placed in 1.5ml tubes and centrifuged at maximum speed for 3min. After centrifugation, the supernatant was removed. Before DNA extraction, 300µl of CTAB buffer was added to the 1.5mL tubes containing the stored pollen. After thorough mixing, pollen samples were transferred to 2.0ml screw cap tubes containing a custom-made mixture of zirconia beads of varying sizes (100µm: 80mg; 200µm: 38mg; 400µm: 85mg; 800µm: 97mg; 3mm: 2 beads) to ensure efficient rupture of pollen exines across a wide range of grain sizes. Pollen samples were then ground in a Precellys 24 tissue homogeniser (Bertin Instruments) for three cycles at 6200 rpm for 5s. DNA extraction was then carried out using the Maxwell® RSC Pure Food GMO and Authentication Kit (Promega) with the Maxwell® automatic extractor (Promega), following the manufacturer’s instructions. Positive and negative controls were included in each extraction batch.

DNA metabarcoding was performed using the dual-indexing approach following Sickel et al. for the internal transcribed spacer 2 (ITS2) nuclear DNA fragment. The first PCR was carried out in triplicate in a 10µl total volume containing 5µl of Q5 High-Fidelity 2X Master Mix (New England Biolabs), 0.5μl of primers ITS-S2F47 (5’-ATGCGATACTTGGTGTGAAT-3’) and ITS-S4R48 (5’-TCCTCCGCTTATTGATATGC-3’) at 10μM, and 1μl of DNA. The PCR cycling profile consisted of 98°C for 3min, 35 cycles of 98°C for 10s, 52°C for 30s, and 72°C for 40s, and a final extension of 72°C for 2min.

Amplicons were then subjected to a second PCR to incorporate unique indexes into each triplicate separately. This indexing reaction was performed in a total volume of 14µl, containing 7μl of KAPA HiFi HotStart ReadyMix PCR Kit (Kapa Biosystems), 1.4μl of each oligonucleotide (1μM), and 2.8μl of each first-round PCR amplicon, previously diluted 1:1. Thermal cycling conditions were 95°C for 3min, followed by 8 cycles of 95°C for 30s, 55°C for 30s, 72°C for 30s, and a final extension of 72°C for 5min. Indexed amplicons were normalised using SequalPrepTM Normalisation Kit 96 wells (Invitrogen) to ensure equalised libraries. After normalisation, all the samples from each PCR plate were combined into one pool.

The concentration of each pool was assessed using the Epoch Microplate Spectrophotometer (BioTek Instruments). The pools were then combined equimolarly, purified using 0.75× reversible immobilisation paramagnetic beads (Agencourt AMPure XP), and eluted in 50μl. Final pools were quality-controlled on a TapeStation 2200 using the HSD5000 kit (Agilent Technologies) and a SYBR green quantitative PCR assay using the KAPA Library Quantification kit (Kapa Biosystems) with the QuantStudio™ 5 Real-Time PCR (ThermoFisher Scientific). The sequencing libraries were diluted to 2nM, spiked with 10% Illumina-generated PhiX control library, and then sequenced on the Illumina MiSeq using 2×250 cycles v2 chemistry according to the manufacturer’s instructions.

Pools were de-multiplexed in the BaseSpace Sequence Hub based on unique indexes. Raw sequencing data were deposited publicly at NCBI SRA under BioProject PRJNA1198597. Data were processed according to Leonhardt et al. and the pipeline available at https://github.com/chiras/metabarcoding_pipeline. Raw sequencing reads were merged, quality filtered (length <200bp and > 500bp, no ambiguous base pairs), dereplicated, denoised, and chimaeras were removed to generate amplicon sequence variants (ASVs). ASVs were directly classified at the species level using global alignments with VSEARCH v2.15.2, a threshold of at least 97% sequence similarity (top-hit selection), as in Sickel et al. Unclassified reads were subjected to SINTAX hierarchical classification, with a threshold of 90%. The ITS2 referenced databases used in both taxonomic classification steps were created by Quaresma et al. Taxa that could not be confidently distinguished at the species but at genus level, i.e., due to recent evolutionary divergence and frequent hybridisation resulting in low interspecific variability in the ITS2 region (e.g., Brassica napus, B. rapa, and B. oleracea), were treated at the genus level only.

Climatedata
Mean monthly temperature and monthly accumulated precipitation were extracted from the high-resolution (0.5°) gridded dataset CRU TS4.08, provided by the Centre for Climate Research at the University of East Anglia. Climate data was only available monthly, resulting in two sampling points per site sharing the same climate metadata.

Statisticalanalysis
Data were analysed using the R v4.4.1 software within the RStudio v2024.04.2+764 environment using the R packages phyloseq v1.48.0, tidyr v1.3.1, speedyseq v0.5.3.9021, dplyr v1.1.4, viridis v0.6.5, bipartite v2.20, ggplot2 v3.5.1, ggsci v3.2.0, gghighlight v0.4.1, rnaturalearth v1.0.1, sf v1.0.16, geosphere v1.5.18, gdm v1.55.0.9.1, scales v1.3.0, foreach v1.5.2, doParallel v1.0.17, doSNOW v1.0.20, progress v1.2.3, patchwork v1.3.0.9000. The R code is available at https://github.com/chiras/HoneyBee-ResistanceResilience. Processed data and metadata of samples are publicly available at this GitHub repository and as well under https://doi.org/10.5281/zenodo.17578272. Metabarcoding data were transformed to relative read abundances (RRAs). Taxa contributing less than 1% to a given sample were excluded to reduce noise and focus on ecologically meaningful contributors to the pollen diet. Multiple ASVs from the same plant species were aggregated into a single taxon to avoid artificial inflation of diversity metrics.

Pollendiversity: For each sample (site x time point), we estimated Shannon entropy (H’) and Effective Species Richness (eH), calculated the crop contribution to diets, and visualised spatial and temporal trends across Europe. This allowed us to examine how floral resource diversity, crop type, and composition vary geographically and seasonally, providing a quantitative baseline (Objective 1) against which projected climate-driven changes can be evaluated. To assess large-scale spatial patterns, we fitted an exploratory linear model relating H′ to latitude, including latitude as a quadratic term to capture a non-linear gradient, as hypothesised mid-latitude effect. Model diagnostics showed no evidence of heteroscedasticity (Breusch–Pagan p=0.24) and approximate normality of residuals (Shapiro–Wilk W=0.98), supported by visual inspection. We also fitted month-stratified models to evaluate seasonal variation in the latitude–diversity relationship. eH and crop contribution were used to visualise patterns across Europe, latitudes, and months.

A non-metric multidimensional scaling (NMDS) using Bray-Curtis dissimilarities was conducted to visualise compositional variation across sites and time points in an unconstrained ordination space. The ordination was performed with five dimensions (k=5), selected based on minimising stress (final stress=0.098). To quantify the contribution of environmental variables, we applied Generalised Dissimilarity Modelling (GDM), which models compositional turnover as a function of geographic and climatic gradients. GDM allowed us to estimate the proportion of variation in pollen composition explained by site, temperature, and precipitation. To identify dominant pollen sources in the dataset, we calculated the 25 most abundant plant taxa across all samples by summing RRAs. These were visualised as stacked bar plots to highlight the relative contribution of key wild and cultivated taxa. This summary aimed to contextualise floral resource composition across Europe (Objective 1), particularly in relation to the reliance on crop species.

Historical developments of temperature and precipitation at the sites: To contextualise the magnitude of projected climate change scenarios (up to +5°C and –50 mm), we examined historical changes in temperature and precipitation at each of the 310 sampling sites over the past 50 years. This was done to ensure that our future scenarios are ecologically realistic and grounded in empirical trends (Objective 2). While large-scale climate change is well established, confirming site-specific warming and drying trajectories is essential for interpreting regional foraging risks. To this end, we applied Wilcoxon paired tests comparing monthly temperature and precipitation between 1970–1975 and 2018–2023 for each site and calendar month. These non-parametric tests were chosen to accommodate deviations from normality and to provide robust, site-level statistical confirmation of historical climate change. Summaries by latitude and season are shown in Supplementary Fig. 5.

Climatic response profiling of plant taxa: Our dataset spanned much of Europe but was irregularly sampled and not spatially continuous. We therefore restricted risk assessments to species at sites where they were directly observed, applying an empirical exceedance analysis. Modelling approaches (e.g., SDMs) were not pursued because they typically rely on spatially continuous covariate layers and denser, raster-like coverage to extrapolate across unsampled areas, which are conditions not met here with our point-local data. In preparation for such risk assessment, we estimated climatic response profiles for each plant taxon by pairing the empirical association between pollen abundance and site and date-level temperature and precipitation (Objective 3). RRA of each species was paired with the corresponding climate values at each sampling location and date. These values were then used to construct non-parametric kernel density estimates (KDEs), yielding species-specific distribution-free approximations of their realised climatic niche envelopes. These KDE profiles reflect the abundance-weighted climatic associations of each species, without assuming a parametric distribution or functional form. These profiles were then used in an exceedance analysis to assess risk under future climate scenarios (see below). To ensure robustness, only species with data from >5 unique sites and that occurred in at least 10 different samples were retained. As illustrative examples, temperature and precipitation response curves are shown for the 25 most abundant taxa in Supplementary Fig. 10.

Exceedance-based floral resource risk assessment: For each site and sampling time point, we calculated the percentile-based risk intervals of the site-specific temperature and precipitation within the empirical KDE-based climatic response curve of each plant taxon (Objective 4). This analysis was applied using both current and projected climate conditions. Projections included incremental warming scenarios (0 to +5°C in 0.5°C increments) and precipitation reduction scenarios (0 to −50 mm per month in −5 mm steps), for each site and time point. Plants were classified into risk categories based on these percentile thresholds: for temperature, values ≤90th percentile threshold were considered low risk, >90th moderate risk, >95th high risk, and >99th critical risk; for precipitation: values ≥10th percentile threshold were low risk, <10th moderate risk, <5th high risk, and <1st critical risk. We then aggregated the RRAs of plants in each risk category per site and time point. This exceedance-based approach allowed us to estimate the proportion of floral resources that may experience climate conditions outside the observed range of their realised niches at each site and date. This method is transparent, data-driven, and conservation-relevant. Conceptually, it aligns with frameworks such as the IUCN Red List, which applies threshold-based criteria to assess ecological stress and extinction risk. These proportions were summarised by latitude to examine broad geographic trends. Means per latitude were visualised using a composite bar plot, showing the effects of increasing temperatures and reduced precipitation across different latitudes. We were further interested in how resource risks change over the sampling season, so we created graphs showing the proportion of pollen resources classified as at least “moderate risk” using the classification described above along the temperature and precipitation gradient at different latitudes. Graphs are shown in the main text only for a subset of scenarios, while all scenarios are visualised in the Supplementary Information. Finally, we fitted a multifactorial linear model to test how the proportion of floral resources at risk was affected by climate change scenario, latitude, sampling month, and their interactions. This allowed us to evaluate whether projected risks varied systematically with geography and season.

Resistance and resilience: To better understand the potential for ecological persistence or recovery under climate change, we quantified resistance and resilience of floral resource communities in response to combined temperature and precipitation stress (Objective 5). We conducted a combined exceedance analysis of precipitation and temperature by calculating, for each sample and scenario combination, the proportion of species classified as at risk according to the above definitions due to either precipitation reduction or temperature increase. Resistance was defined as the proportion of non-risk taxa in the community, computed as:

Resistance=1 − S risk / (S no-risk + S risk)

where S risk is the number of taxa in the sample classified as at least “moderate risk” under the given scenario, and S no-risk is the number classified as “low risk”.

To assess resilience potential, we identified replacement taxa that could hypothetically compensate for taxa at risk (S risk), assuming range shifts or phenological substitution across space or time. A taxon was considered a suitable replacement for a given at-risk species if it met the following criteria: (1) occurred within a 500 km radius of the focal sample (i.e., geographically proximal); (2) was classified as “low risk” for both temperature and precipitation under the same climate scenario; (3) belonged to the same plant family, serving as a proxy for morphological and nutritional similarity and for maintaining phylogenetic representation; (4) had not already been selected as a replacement in that sample, to avoid redundancy and maintain taxonomic diversity. For each plant taxon in S risk, we searched for suitable replacement candidates across the full dataset. Resilience was then defined as the proportion of at-risk taxa for which at least one suitable replacement taxon was identified, calculated as:

Resilience= S replacement / S risk

where S replacement is the number of at-risk taxa for which replacements were available under the criteria above.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The article studies the effects of climate change on honey bee-plant interactions. This is done using a huge dataset of pollen samples that were collected by honey bees in a host of European countries. The pollen were identified using DNA metabarcoding. The authors applied various data analyses to conclude that climate change will negatively affect foraging resources of honey bees. A R-based pipeline accompanies the article, which I greatly appreciate as it facilitated the review work greatly (despite the somewhat unhandy sequential nature of said pipeline).\nIn the introduction, the authors state that they selected honey bees as their model organism, which is fine. However, the statement “this makes it an ideal species for assessing the broader impacts of climate change on plant-pollinator interactions” is overly general. I would like to point out (having some, but very limited, experience with community datasets of pollinator-plant interactions) that most bee species exhibit wildly different (foraging) behavior from honey bees. Thus, any general statements as to the broader impacts of the study should be limited to social, domestic, bees. Honey bees are not a good model organism for social wild bees, and not at all for solitary wild bees. That does not take away from the importance of the study, honey bee pollination is critical for society, and understanding how this might change in the future is thus also vital.\nI am a statistical ecologists specializing in statistical methods for multivariate analysis, and so I will mostly focus on the data analyses, while (mostly) refraining from the ecological conclusions from the authors. Another reviewer is more suitable to comment on those aspects of the article. The authors should be aware that there is an ongoing shift in statistical community ecology in moving away from “traditional” methods used for data analyses, and more towards model-based approaches. In this debate, I am on the model-based “side”, which leads me to being overly critical on the motivation for and execution of the data analyses in the article.\nIn the “main text” the authors note that “to assess potential risks to pollination networks, they examined climatic niche distributions of the plant species in the pollen samples, while estimating the likelihood of interaction losses under climate change scenarios”. When diving further into the methods, the authors categorize their analyses in 1) pollen diversity, 2) historical developments of temperature and precipitation, 3) plant temperature and precipitation niches, 4) food resource risk assessment, and 5) resistance and resilience. These analyses by all accounts in line with most what most of community ecology does. Unfortunately, most of community ecology poorly understands the motivation for applying certain analyses, and from a statistical standpoint, most of the analyses in community ecology have little merit. I continue my review by the same categorization that the authors followed.\nI note that most of the analyses are insufficiently motivated with respect to the goal that the authors outline in the main text: what questions, subgoals, or otherwise do these analyses answer? I urge the authors to carefully consider how each analysis or method contributes to the goal of the paper, and elaborating this in the article. Some of the details necessary to understand what is done in the analyses are actually covered in the “results” section, while the “methods” section omits large amounts of detail.\n1) Pollen diversity\nL491-2: “Effective species richness was estimated plotted on a map”, why? Which goal does this serve or which question does it answer?\nL493: “We fitted a quadratic model of effective species richness explained by latitude”, why? What kind of quadratic model? What assumptions were checked to make sure that the results were valid? Does the model predict effective species richness well or did the authors not consider its predictive accuracy at all? I suspect the latter.\nL493-495: “A NMDS was conduced to determine the effect of temperature and precipitation” I fail to see how an NMDS contributes to this effect: NMDS is an unconstrained ordination, and only post-hoc can the axes (which are unitless and rotationally variant) be “tested” for the “effect” of precipitation and temperature. Here, “testing” means to check if the scores of the axes (i.e., in an ordination plot) are structured in some way: increasing in temperature from left to right or decreasing in precipitation from top to bottom, for example, but in no way can the authors from this conclude that pollen diversity is affected by these environmental variables.. The authors fail to clarify the number of axes in the NMDS (usually determined in terms of stress), while the two used axes might very poorly represent the authors’ data. Quantifying axis relevance in NMDS is not a trivial endeavour (unlike in other classical ordination methods that come with accompanying eigenvalues), though it is generally done by looking at the “stress” measure that NMDS minimizes. Regardless, a more careful consideration of the number of axes would not have made this a more convincing analysis. Afterall, NMDS is not a method for estimating effect of environmental variables, this is done using constrained ordination methods (RDA or CCA) or more advanced statistical methods such as vector GLMs, reduced rank vector GLMs (i.e., model-based constrained ordination) or generalized linear latent variable models (model-based ordination and JSDMs generally).\n- L495-497: “GDM was used to estimate how much of the pollen composition was explained by site, temperature, and precipitation”. Here, the authors are admitting using a tool to get a metric out, which (by all accounts) is considered by statisticians to have little merit for scientific importance. If the authors so wished, the same metric could have been calculated from their post-hoc analysis of NMDS scores, but instead the authors applied another method, which could have been entirely avoided.\n- L497-498: Barplots were made for summarizing the most abundant pollen. Why? What goal do these serve?\n2) Historical developments\nThe authors used a Wilcoxon paired test to test in the mean temperature and precipitation of the two time windows has changed. Personally, I think this is a well accepted fact in science, and thus again do not see the added value of this test to the study.\n3) Plant temperature and precipitation niches\nL508-509: “Temperature data were rounded to 0.5°C, and\nprecipitation data were square-root transformed due to a skewed distribution”. Why did the authors round the temperature data, and why is a skewed distribution for the precipitation data a problem? I worry that the authors do not understand the assumptions that are made with their analyses, how to check the consequences of these assumptions, and consequences (or lack thereof) of violating certain assumptions (such as normality).\nL509-512: “Temperature and precipitation niche distributions for each plant species were estimated as normal distributions given the mean and standard deviation of temperature and precipitation, respectively, weighted by the abundance data over all samples.”\nAll details with respect to the analyses are omitted: one can only understand what the authors have done by digging into the accompanying github repository, because generally the article omits so much detail that it would be impossible to understand, let alone to reproduce, without the accompanying code.\nFrom the text, “estimate” would indicate to me that the authors fitted a regression of some sort (the reference to transformed skewed data would imply a linear regression), but the code makes clear that the authors do not fit a model at all. What the authors do takes some time to figure out, as the code is not extensively commented throughout. The authors calculated the sample mean and sample variance of the relative abundances of plant species, but how these constitute valid estimates of “distribution” is entirely unclear to me (i.e., a statistical ecologists that frequently works for methods of species distributions).\nThe only product in the article in terms of “distribution” that I can find is Figure 1A, and Extended data Figure 10. Figure 1A plots effective richness on a map, and thus does not constitute a predicted distribution of any kind.\nExtended Data Figure 10 is also not a “distribution”, but it is a response curve, despite the fact that the authors call it a distribution. In Extended Dat Figure 10 the sample mean and sample variance of summed relative abundances are used to plot a Gaussian kernel, without any regard for the accuracy of the response curve (since it is not a distribution, while response curve might be an apt name for what is shown here), and how it is relevant for the original data. The authors cannot make statements about this, as no model has been fitted whatsoever.\nThis can be seen as an attempt at a linear regression, where temperature and precipitation are separately used as response variables, and “accumulated species” is used as a categorical variable, thus estimating the mean of temperature and precipitation per species, but while allowing for heterogeneity of variances per “accumulated species”. Making that connection exposes multiple deficits with the authors’ approach: 1) the code could be much more concise, 2) a normality assumption for precipitation is not realistic: precipitation cannot be negative, 3) a skewed error distribution is not a problem, the authors are only interested in prediction, and violation of the normality assumption only affects the statistical uncertainty in a linear regression (and the authors are not quantifying statistical uncertainty at all). In fact, by using the Wilcoxon signed-rank test on temperature and precipitation the authors already account for the non-normality effect on test results, 4) the authors are not analysing their relative abundances: they are analysing temperature and precipitation while stating that they are analysing the species data, i.e., the authors are confusing the relationship of y and x: temperature and precipitation affects the abundance of species, and the distribution of species is predicted using temperature and precipitation. Here, the authors are using species abundances to calculate the (mean) temperature or precipitation of sites where species were present, which effectively ignores nearly all information in the pollen data. Although the x and y relationship can be reversed without consequence when there is only one variable impacting the relative abundance of pollen at sites in a linear regression, this is not true when we assume that there is more than one variable affecting the pollen, as is there case here. Further, there are most likely confounders that affect pollen distribution, which are also not accounted for in this analysis. In summary, the collinearity of precipitation and temperature is very likely to impact the results, and separate linear regression cannot be fitted to the data while expecting accurate results.\n4) Food resource risk assessment\nAs noted, the calculation of distributions is inadequate, and 4) is based on 3), so that 4) too needs a considerable rework. Further, note that in L517-518 the authors write “the quantile was calculated”, but there is no such thing as “a quantile” as a quantile defined is a range, the word does not take meaning until that range is defined. This is in contrast to percentiles.\nIn this part of the analysis, the authors fit a linear model and perform model selection using AIC, after which statistical significance is assessed. I have two notes on this: 1) stepwise model selection has considerable issues that the authors should address, details on this can be found widely in the (ecological) literature. The debate is not “hot” anymore, and AIC model selection has largely been put on the back-burner due to its widely recognized flaws, 2) model selection by information criteria cannot be combined with the interpretation of statistical significance. Omitting variables by information criteria naturally leads to better fitting models, and thus artificially increased statistical significance (inflated Type I error). This is commonly referred to as “p-hacking”, and is considered a statistical sin. The two paradigms: information criteria and testing (exploratory and confirmatory), are not unifiable.\n5) Resistance and resilience\nThis calculation depends on the previous analyses with considerable flaws, and thus also needs a rework. At the same time, the equation the authors provide in the text is not referenced, so determined whether or not this is a sensible estimate of resistance (or resilience) is entirely unclear.\nOverall, I conclude that the analyses in the article lack both statistical support and ecological motivation. I would have personally reanalysed the data, but due to the sequential nature of the pipeline, and the lack of intermediate products stored, I could not find sufficient time to do the analysis. My recommendation is for the authors to completely redo their analyses, as it can be done considerably more succinctly, and because I expect various (if not all) of the results to change. I recommend that the authors dive into the world of multivariate statistical models, as at least three, if not four, of their analysis can be made redundant by fitting a single vector GLM (as in the mvabund R-package) or single GLLVM (as in the gllvm R-package). The compositional nature of the data can be addressed by appropriately selecting the model structure: the total read count should be included as an offset, and the response distribution should be one of Poisson, negative-binomial or zero-inflated versions thereof. I would be happy to provide the authors with more details for a reanalysis if needed.\nFor now, I have refrained from extensive commenting on the text. That is not by lack of interest – the study is interesting and I would like to see an improved version. However, I expect the present version to still undergo significant changes when the analyses have been improved, so that commenting on (e.g.,) the results is not relevant at this stage.\nBert van der Veen, 14-03-2025\n(Remarks on code availability)\nThe code is available, and can be run without (much) further intervention, on which the authors have done a reasonable good job.\nThe code is not well commented throughout, so that working through the code, and understanding what the authors did, took a considerable amount of time, and was challenging. Ideally, this should be improved. Having said that, it would not have been possible to understand what the authors did in the article, without reviewing the code.\nNote that in the first \"Umbrella\" script the authors incorrectly (for anyone else accessing the code) set a working directory.\nThe sequential nature of the pipeline, without storage of and reference to intermediate results, makes reanalysis extremely difficult.\nThe pipeline would function better (for reanalysis) if each analysis was stand-alone.", "reviewer_2": "The manuscript presented to me is very well written. It was a pleasure to read it, and I found it difficult to point out any major shortcomings. The topic is very significant, and the research is original and important. The methods are suitable, although honeybees are not the best representative of bees generally, no other species would allow for such a whole-season analysis and a sufficient number of samples. The Authors are aware of this and discuss their result accordingly. I only missed one piece of information about the apiaries. Did the Authors use only stationary apiaries for their sampling? Moving apiaries into or close proximity to mass flowering and honey-giving crops is a common practice and may affect the results.\nThe description of the results is very detailed, and the Authors did a good job discussing the eventual shortcomings of the climate scenario results. I think it was even a bit too good of a job, and part of these descriptions/interpretations of the received results should be moved to the Discussion, which is short and lacking these interpretations.\nI also found a few minor points, which could be improved. In Fig. 1A. the sampling points are somewhat outside the map. I understand that it is because some samples were coming from islands like Lesvos and Cyprus? Can you change the map underneath to show not just the EU countries but simply all the countries? It would make it easier for the reader to locate the most southeastern sample locations correctly.\nFig2. Can you divide it into two figures? One showing crops, the other wild plants? Although I understand that it shows the total impact, it is difficult to differentiate visually in this form.\n(Remarks on code availability)\nYes, I have checked the URL however, I don't feel competent to review it.", "reviewer_3": "This is a research paper reporting findings from a 4 month study in 2023 across a very impressive number and range of sites in Europe where the authors collected honeybee pollen from aviaries at each of these sites simultaneously every 2 weeks over a 17 week period (presumably). The findings indicate the range of flowering plant species that honeybees visit over May to August at 310 locations, and relate the number to environmental variables. The results provide evidence of the reduced diversity of flowering plants that honeybees are feeding on, and this appears to be related to increased temperature and reduced rainfall. There is no data however on any negative effect on the honeybees, nor any data on forage available at the time of pollen collection. Given the limited range of flowering plants honeybees feed from (seemingly flagged in this paper though framed a little differently), the findings might have major implications for other flowering species and other visiting insects, particularly in light of the fact that we know pollinators are in decline.\nI would suggest a re-framing of the use of honeybees as a model, given I believe they do have a preference for a small number of highly abundant flowering plants, as flagged in this study.\nThere is some lack of clarity and minor suggestions for changes below.\nLine 23+ Abstract is missing elements that would be useful, re the dataset, timing, to make it clear that this was one year. This is modelling based on scenarios.\nLine 38 Some odd phrasing, e.g. how the word ecology is used in the first sentence, I also question putting value on different ecological interactions.\nLine 66 requires a citation. Honeybees are generalists but by focusing on this group you miss the deeper flowers that long tongued groups depend on, and the species that have few inflorescences which I think needs mentioning. Managed honeybees collect pollen which can be very easily located and collected making this a useful group to use to investigate. What’s more, your findings further suggest that honeybees feed on a small number of presumbaly abundant flowering plants. This to me might make the unknown element, the impact of the predicted changes in your analysis on flower visiting insect populations alarming as you mention in the discussion. But flags the importance of monitoring impacts on other species, in addition to an understanding of what is available at these locations at the time.\nLine 85 it isn’t clear to me whether latitude and temp and/or rain are correlated which is somewhat critical for the interpretation?\nLine 74 2500 in total, simultaneously every 2 weeks from 310 apiaries, one apiary per location? May-August in 2023 is that then 8 sampling dates? Which is 2480 samples. This could just clarified, and could be clearer in the intro and I think should be stated in the abstract at least the timing as this is important.\nLine 206 - I’m not certain these results (species richness 5-10) indicates a diverse range of collected pollen, some summary stats would be useful to clarify this and to clarify the proportions of the top flowering plant groups that honeybees are feeding on in the results.\nLine 428 why not a specific date?\nLine 484 R studio is a platform for analysing data using R software, I think this needs to be referred to appropriately.\nLine 482 this indicates pseudo-replication, it’s unclear whether this was accounted for in modelling?\n(Remarks on code availability)", "reviewer_4": "The manuscript under review, “Honey bee food resources under threat from climate change,” provides valuable insights into the restoration of plant–animal interactions amidst rapid environmental changes and increasing human interference. It effectively conveys the urgency of monitoring and maintaining the balance of species interactions. The manuscript is generally well-written, with a clear and coherent presentation. However, I have a few concerns and suggestions that I believe could strengthen the study:\nSampling Period: The sampling timeframe appears narrow for drawing conclusions about climate impacts. I recommend extending the study over at least 1–2 additional years to capture interannual variability.\nMethodological Comments and Suggestions:\nLine 92: It would be beneficial to include a main figure showing plant abundance data stratified by location.\nLines 205–207: Could the authors comment on the occurrence of low-abundance pollen taxa and whether they might result from secondary contamination?\nIt would be valuable to include a comparison across multiple bee species to assess species-specific preferences in pollen foraging.\nLine 438: Please clarify how pollen samples were transferred into 2 ml tubes. Also, include a reference for the zirconia bead size-sorting procedure.\nLine 445: Consider placing the Maxwell® 440 RSC PureFood kit information here, and include product details in parentheses.\nLine 450: The strategy of maintaining equal DNA concentrations for each sample at the start of PCR is commendable. Did this approach apply to all samples, including those with low concentrations?\nLine 454: Please elaborate on how subsamples were selected for indexing. Was pooling involved, and if so, what was the basis? The “1:1 dilution of the first amplicons” requires further clarification.\nPooling Strategy: How were pools equimolarly combined? Please provide methodological detail.\nLine 473: The use of 97% sequence similarity—do the authors believe this threshold is sufficient for species-level resolution?\nLines 475–477: The method used to merge undetermined taxa with Brassica is unclear. Please elaborate.\nLine 472: Which database(s) were used for taxonomic assignment?\nOverall, this manuscript addresses an important and timely topic. With these clarifications and additions, it could make a strong contribution to our understanding of pollinator–plant dynamics under climate stress.\n(Remarks on code availability)"}, "author_response": "Dear Reviewers,\nWe sincerely appreciate the valuable and constructive feedback provided on our\nmanuscript. We have carefully addressed all comments, as detailed below, and agree that\nthe revisions have strengthened the manuscript.\nReviewer #1 (Remarks to the Author):\nThe article studies the effects of climate change on honey bee-plant interactions. This is\ndone using a huge dataset of pollen samples that were collected by honey bees in a host\nof European countries. The pollen were identified using DNA metabarcoding. The authors\napplied various data analyses to conclude that climate change will negatively affect\nforaging resources of honey bees. A R-based pipeline accompanies the article, which I\ngreatly appreciate as it facilitated the review work greatly (despite the somewhat unhandy\nsequential nature of said pipeline).\n>> We thank the reviewer for their thoughtful and constructive comments, which have\nsignificantly improved the manuscript. We also appreciate the recognition of the R-based\npipeline and acknowledge the limitations of its sequential structure. We have now\nclarified in the code and repository documentation that key intermediate results are saved\nduring runtime. However, certain steps, such as the generation of plant-specific density\nfunctions, require sequential execution and recall in subsequent parts, due to their\ndependency on in-memory structures, and cannot be modularised.\nIn the introduction, the authors state that they selected honey bees as their model\norganism, which is fine. However, the statement “this makes it an ideal species for\nassessing the broader impacts of climate change on plant-pollinator interactions” is\noverly general. I would like to point out (having some, but very limited, experience with\ncommunity datasets of pollinator-plant interactions) that most bee species exhibit wildly\ndifferent (foraging) behavior from honey bees. Thus, any general statements as to the\nbroader impacts of the study should be limited to social, domestic, bees. Honey bees are\nnot a good model organism for social wild bees, and not at all for solitary wild bees. That\ndoes not take away from the importance of the study, honey bee pollination is critical for\nsociety, and understanding how this might change in the future is thus also vital.\n>> We thank the reviewer for this insightful comment. We agree that honey bees differ\nsubstantially in foraging behaviour from wild and solitary bees due to traits such as\ncolony-level recruitment and generalist foraging. Our study does not aim to generalise\npollinator behaviour but to assess how climate change may affect floral resources\navailable to Apis mellifera, and the resulting nutritional implications. While species-\nspecific responses will differ, honey bees forage on a broad range of plants, many of\nwhich are shared with other pollinators. Thus, the patterns of floral resource shifts we\nreport remain relevant for agroecosystem dynamics and may have indirect implications\nfor other bee taxa. We have clarified this scope in the Introduction and Discussion and\nremoved the \"model organism\".\nI am a statistical ecologists specializing in statistical methods for multivariate analysis,\nand so I will mostly focus on the data analyses, while (mostly) refraining from the\necological conclusions from the authors. Another reviewer is more suitable to comment\non those aspects of the article. The authors should be aware that there is an ongoing shift\nin statistical community ecology in moving away from “traditional” methods used for\ndata analyses, and more towards model-based approaches. In this debate, I am on the\nmodel-based “side”, which leads me to being overly critical on the motivation for and\nexecution of the data analyses in the article.\n>> We thank the reviewer for their thoughtful remarks and for outlining their perspective\non statistical methodology in community ecology. We acknowledge the increasing role of\nmodel-based approaches and appreciate their utility. Our study does in fact employ\nseveral model-based methods, including Generalised Dissimilarity Modelling (GDM)\nand multifactorial linear models. That said, we agree that the rationale for specific\nanalyses was not always fully articulated, partly due to prior formatting constraints\nimposed by a different journal (see reply to reviewer 2), where methods were detached\nand significantly shortened. We have now revised the manuscript to clarify the motivation\nfor each analytical step.\nWhile we recognize the value of methodological reflection, we respectfully note that\ngeneral commentary on broad statistical trends, particularly when not accompanied by\nspecific model suggestions or references (e.g., statistical vs. mechanistic frameworks,\nJSDMs, or vector GLMs), may be better suited to an own perspective piece rather than a\nmanuscript review. We remain confident that our analytical framework is appropriate to\nthe questions posed and the data structure at hand, and we are grateful to the reviewer for\nhelping us strengthen the presentation of our approach.\nIn the “main text” the authors note that “to assess potential risks to pollination networks,\nthey examined climatic niche distributions of the plant species in the pollen samples,\nwhile estimating the likelihood of interaction losses under climate change scenarios”.\nWhen diving further into the methods, the authors categorize their analyses in 1) pollen\ndiversity, 2) historical developments of temperature and precipitation, 3) plant\ntemperature and precipitation niches, 4) food resource risk assessment, and 5) resistance\nand resilience. These analyses by all accounts in line with most what most of community\necology does. Unfortunately, most of community ecology poorly understands the\nmotivation for applying certain analyses, and from a statistical standpoint, most of the\nanalyses in community ecology have little merit. I continue my review by the same\ncategorization that the authors followed.\n>> While we acknowledge the comment regarding broader trends in community ecology,\nwe emphasize that our analyses were selected based on their relevance to our research\nquestions and are supported by established methods. We have now clarified our\nmotivations by clearer objectives and more explicitly in the revised manuscript, especially\nat the end of the introduction.\n>> With respect to the reviewer’s general statement that “most of community ecology\npoorly understands the motivation for applying certain analyses,” we respectfully note\nthat such broad characterisations fall outside the scope of this manuscript and do not\ncontribute constructively to the review of our work\nWe also note that none of the other reviewers raised concerns about the suitability or\njustification of the analyses. Nevertheless, we are grateful for the opportunity to improve\nthe choice, clarity and framing of our methods and hope that the revisions address the\nreviewer’s concerns.\nI note that most of the analyses are insufficiently motivated with respect to the goal that\nthe authors outline in the main text: what questions, subgoals, or otherwise do these\nanalyses answer? I urge the authors to carefully consider how each analysis or method\ncontributes to the goal of the paper, and elaborating this in the article. Some of the details\nnecessary to understand what is done in the analyses are actually covered in the “results”\nsection, while the “methods” section omits large amounts of detail.\n>> We agree that the connection between our overarching aim and the individual analyses\ncould be more clearly articulated. In response, we have revised the Introduction to\nexplicitly state the five analytical objectives and added brief linking statements at the start\nof each relevant subsection to clarify how each analysis contributes to addressing the\noverall research question. We agree, these revisions improve the logical flow and clarify\nthe purpose of each analytical step. Thank you for pointing this out!\n1) Pollen diversity\nL491-2: “Effective species richness was estimated and plotted on a map”, why? Which\ngoal does this serve or which question does it answer?\n>> We clarified that mapping effective species richness serves to address Objective 1 by\nproviding a baseline of floral resource diversity across Europe, against which climate-\ndriven changes can be evaluated.\nL493: “We fitted a quadratic model of effective species richness explained by latitude”,\nwhy? What kind of quadratic model? What assumptions were checked to make sure that\nthe results were valid? Does the model predict effective species richness well or did the\nauthors not consider its predictive accuracy at all? I suspect the latter.\n>> We thank the reviewer for this important point. The quadratic model was used for\nexploratory purposes, specifically to examine broad spatial patterns in floral resource\ndiversity (Objective 1), including the hypothesised mid-latitude peak. We have now made\nthis motivation more explicit in the text. As the reviewer correctly notes, this analysis was\nnot intended for prediction, but to describe the geographic structure of diversity across\nsites and seasons. While we understand that the reviewer may prioritise predictive\napproaches, we consider this valuable, particularly given that bee foraging diversity has\nnot previously been assessed at this scale. The final model showed no evidence of\nheteroscedasticity (Breusch–Pagan p = 0.24) and approximate normality (Shapiro–Wilk\nW = 0.98), supported by residual diagnostics. As noted in the Results, the model’s\nexplanatory power was low (adjusted R² = 0.10), which is consistent with its exploratory\npurpose.\nL493-495: “A NMDS was conduced to determine the effect of temperature and\nprecipitation” I fail to see how an NMDS contributes to this effect: NMDS is an\nunconstrained ordination, and only post-hoc can the axes (which are unitless and\nrotationally variant) be “tested” for the “effect” of precipitation and temperature. Here,\n“testing” means to check if the scores of the axes (i.e., in an ordination plot) are\nstructured in some way: increasing in temperature from left to right or decreasing in\nprecipitation from top to bottom, for example, but in no way can the authors from this\nconclude that pollen diversity is affected by these environmental variables.. The authors\nfail to clarify the number of axes in the NMDS (usually determined in terms of stress),\nwhile the two used axes might very poorly represent the authors’ data. Quantifying axis\nrelevance in NMDS is not a trivial endeavour (unlike in other classical ordination\nmethods that come with accompanying eigenvalues), though it is generally done by\nlooking at the “stress” measure that NMDS minimizes. Regardless, a more careful\nconsideration of the number of axes would not have made this a more convincing analysis.\nAfterall, NMDS is not a method for estimating effect of environmental variables, this is\ndone using constrained ordination methods (RDA or CCA) or more advanced statistical\nmethods such as vector GLMs, reduced rank vector GLMs (i.e., model-based constrained\nordination) or generalized linear latent variable models (model-based ordination and\nJSDMs generally).\n>> We thank the reviewer for this detailed and helpful critique. We fully agree that\nNMDS is an unconstrained ordination method and does not allow for direct inference\nabout environmental drivers. We acknowledge that our original phrasing was imprecise,\nand have revised the manuscript to clarify that NMDS was used solely to visualise\ncompositional dissimilarity among samples, not to assess environmental effects. The\nstress value (0.098) is now included in the Supplementary Figure legend, and the number\nof NMDS axes (k = 5) is clarified in the Methods.\nL495-497: “GDM was used to estimate how much of the pollen composition was\nexplained by site, temperature, and precipitation”. Here, the authors are admitting using\na tool to get a metric out, which (by all accounts) is considered by statisticians to have\nlittle merit for scientific importance. If the authors so wished, the same metric could have\nbeen calculated from their post-hoc analysis of NMDS scores, but instead the authors\napplied another method, which could have been entirely avoided.\n>> Regarding the use of Generalized Dissimilarity Modelling (GDM), we respectfully\ndisagree with the characterisation that it lacks scientific merit. GDM is a well-established\nmodel-based approach designed to quantify spatial and environmental turnover in species\ncomposition. Unlike NMDS, it allows explicit modelling of dissimilarity as a nonlinear\nfunction of geographic and environmental gradients, and is widely used in biodiversity,\nmacroecology, and biogeography. We selected GDM precisely because it accommodates\ncomplex species–environment relationships without relying on assumptions of linearity\nor normality, and because it is well suited to large-scale, high-dimensional community\ndatasets such as ours. We believe this makes it more appropriate than either constrained\nordination or regression on NMDS scores, given the structure and scope of our data. We\nhave revised the manuscript to make the distinction between NMDS and GDM clearer.\nL497-498: Barplots were made for summarizing the most abundant pollen. Why? What\ngoal do these serve?\n>> We thank the reviewer for raising this point. The barplots and summaries of the most\nabundant pollen taxa support Objective 1, which is to characterise the composition and\ndiversity of floral resources collected by honey bees across Europe. Identifying dominant\ntaxa (both wild and cultivated) provides essential ecological context for interpreting\npatterns of diversity, potential vulnerability, and seasonal dynamics. This compositional\noverview also informs Objectives 4 and 5, which assess climate-driven risks and the\npotential for resilience. We have now clarified this rationale in the Methods and Results.\nWe also wish to emphasise the value of such large-scale descriptive data, which remain\nrare despite their importance for understanding foraging ecology at continental scale.\n2) Historical developments\nThe authors used a Wilcoxon paired test to test in the mean temperature and precipitation\nof the two time windows has changed. Personally, I think this is a well accepted fact in\nscience, and thus again do not see the added value of this test to the study.\n>> We thank the reviewer for this comment. While overall climate warming is well\nestablished, regional and even more local patterns can vary significantly. The included\nthe Wilcoxon paired tests were to confirm that such changes are detectable at the scale of\nour actual sampling locations. This ensures that the historical context of our study is\ngrounded in site-specific climate trajectories, which is important for interpreting scenario-\nbased projections. In this context, we now adapted the range of simulations for\nprecipitation to reflect a better range of historical developments on local scale. We have\nclarified this rationale in the revised Methods. We also refer the reviewer to the\naccompanying figure and caption for further context.\n3) Plant temperature and precipitation niches\nL508-509: “Temperature data were rounded to 0.5°C, and precipitation data were\nsquare-root transformed due to a skewed distribution”. Why did the authors round the\ntemperature data, and why is a skewed distribution for the precipitation data a problem?\nI worry that the authors do not understand the assumptions that are made with their\nanalyses, how to check the consequences of these assumptions, and consequences (or lack\nthereof) of violating certain assumptions (such as normality).\n>> The reviewer is correct that skewness is not inherently problematic unless parametric\nassumptions are made. Our initial approach relied on estimating species-specific climatic\nmeans and standard deviations, which assumes approximately normal distributions. This\napproach followed established works, as e.g.:\nKühsel, S., Blüthgen, N. High diversity stabilizes the thermal resilience of\npollinator communities in intensively managed grasslands. Nature\nCommunications 6, 7989 (2015). https://doi.org/10.1038/ncomms8989,\nespecially the framework illustrated in Figure 1.\nRounding temperature values to 0.5 °C and applying a square-root transformation to\nprecipitation were intended to reduce noise and skewness in support of this parametric\nframework.\nHowever, we recognise that normality may not hold across all species, and our initial,\nmanual validation focused only on the most abundant taxa. In response, we have now\nreplaced the parametric approach with a non-parametric kernel density estimation (KDE)\nmethod. This allows us now to estimate plant–climate associations without assuming\nnormality and increases the robustness of our climate risk estimates. In fact, the results\nand interpretation did not change. New changes are reflected in the revised Methods\nsection and implemented in the updated analysis pipeline.\nL509-512: “Temperature and precipitation niche distributions for each plant species\nwere estimated as normal distributions given the mean and standard deviation of\ntemperature and precipitation, respectively, weighted by the abundance data over all\nsamples.”\nAll details with respect to the analyses are omitted: one can only understand what the\nauthors have done by digging into the accompanying github repository, because generally\nthe article omits so much detail that it would be impossible to understand, let alone to\nreproduce, without the accompanying code.\nFrom the text, “estimate” would indicate to me that the authors fitted a regression of\nsome sort (the reference to transformed skewed data would imply a linear regression),\nbut the code makes clear that the authors do not fit a model at all. What the authors do\ntakes some time to figure out, as the code is not extensively commented throughout. The\nauthors calculated the sample mean and sample variance of the relative abundances of\nplant species, but how these constitute valid estimates of “distribution” is entirely unclear\nto me (i.e., a statistical ecologists that frequently works for methods of species\ndistributions).\n>> We thank the reviewer for this close and thoughtful reading. We acknowledge that the\noriginal text may have been unclear and thereby led to misunderstandings about both our\ngoals and the methods used. In particular, our use of the term “distribution” was intended\nin a statistical sense (i.e., a response profile over climate space), but we recognise that it\nmay be interpreted as referring to spatial distributions or even species distribution models\n(SDMs), which we do not apply. We have now revised the terminology throughout to\nrefer instead to “climatic response profiles” or “abundance-weighted climatic\nassociations.”\n>> As the reviewer correctly observes, our analysis does not involve modelling species\nabundance or presence as a function of climate via regression. In our original approach,\nwe summarised the climatic conditions associated with each plant taxon by calculating\nabundance-weighted means and standard deviations of temperature and precipitation.\nHowever, we agree that this parametric framework assumes normality and may not\nadequately capture all plants' climatic associations, therefore we switched to the KDE\napproach as detailed in the previous comment response.\n>> These KDE-derived response profiles are then used to assess species-specific climate\nsensitivity by evaluating whether projected site-level temperature and precipitation values\nfall outside the observed climatic envelope of each species. This is not a predictive SDM,\nbut a data-driven exceedance analysis designed to highlight the potential for mismatch\nbetween future local climates and observed foraging niches. We have substantially revised\nthe Methods section to clarify each step of this procedure and ensure that it is\nunderstandable and reproducible without referring to the source code.\nThe only product in the article in terms of “distribution” that I can find is Figure 1A, and\nExtended data Figure 10. Figure 1A plots effective richness on a map, and thus does not\nconstitute a predicted distribution of any kind.\n>> We have clarified that Figure 1A shows observed, not predicted, effective species\nrichness across sites.\nExtended Data Figure 10 is also not a “distribution”, but it is a response curve, despite\nthe fact that the authors call it a distribution.\n>> We agree that the term “distribution” may lead to confusion with species distribution\nmodels (SDMs) rather than statistical representations. We now refer to the curves shown\nin Supplementary Figure 10 as empirical response curves not to be confused with fitted\nor theoretical distributions. Terminology has been updated throughout the manuscript to\nimprove clarity.\nIn Extended Dat Figure 10 the sample mean and sample variance of summed relative\nabundances are used to plot a Gaussian kernel, without any regard for the accuracy of\nthe response curve (since it is not a distribution, while response curve might be an apt\nname for what is shown here), and how it is relevant for the original data. The authors\ncannot make statements about this, as no model has been fitted whatsoever.\n>> We thank the reviewer for this comment and appreciate the opportunity to clarify. In\nthe original version, Supplementary Figure 10 displayed Gaussian curves based on\nabundance-weighted means and standard deviations, following a parametric framework\nused in previous works (e.g. Kühsel & Blüthgen 2015, Nat. Commun.). However, we\nagree that this approach imposes assumptions that may not hold valid across all taxa, was\nonly verified manually for a subset, and does not constitute a model fit. In response to the\nreviewer’s concern, and in line with above mentioned changes, we have replaced the\nGaussian curves following with non-parametric kernel density estimates (KDEs), which\nmore faithfully represent the empirical climate–abundance associations of each species\nwithout assuming normality. These KDEs directly reflect the response curves used in the\nrisk analysis, and the figure is intended solely for descriptive illustration of the 25 most\nabundant taxa. We have updated the figure, caption, and main text.\nThis can be seen as an attempt at a linear regression, where temperature and precipitation\nare separately used as response variables, and “accumulated species” is used as a\ncategorical variable, thus estimating the mean of temperature and precipitation per\nspecies, but while allowing for heterogeneity of variances per “accumulated species”.\nMaking that connection exposes multiple deficits with the authors’ approach: 1) the code\ncould be much more concise,\n>> We have revised the relevant sections in the code to improve clarity and conciseness.\nWe apologize that there were still remnants and variables of older workflows not used\nanymore and might have caused confusions later. Such lines were removed.\n2) a normality assumption for precipitation is not realistic: precipitation cannot be\nnegative,\n>> We agree and have replaced this method with a non-parametric, kernel density-based\napproach that does not assume any specific distributional form. This change improves\nrobustness, especially for skewed or bounded variables like precipitation.\n3) a skewed error distribution is not a problem, the authors are only interested in\nprediction, and violation of the normality assumption only affects the statistical\nuncertainty in a linear regression (and the authors are not quantifying statistical\nuncertainty at all). In fact, by using the Wilcoxon signed-rank test on temperature and\nprecipitation the authors already account for the non-normality effect on test results,\n>> We agree and note that our revised kernel density-based approach no longer relies on\ndistributional assumptions. To clarify, the Wilcoxon test was used only for historical\nclimate trends at sampling locations, not for plant response data.\n4) the authors are not analysing their relative abundances: they are analysing\ntemperature and precipitation while stating that they are analysing the species data, i.e.,\nthe authors are confusing the relationship of y and x: temperature and precipitation\naffects the abundance of species, and the distribution of species is predicted using\ntemperature and precipitation. Here, the authors are using species abundances to\ncalculate the (mean) temperature or precipitation of sites where species were present,\nwhich effectively ignores nearly all information in the pollen data. Although the x and y\nrelationship can be reversed without consequence when there is only one variable\nimpacting the relative abundance of pollen at sites in a linear regression, this is not true\nwhen we assume that there is more than one variable affecting the pollen, as is there case\nhere.\n>> We thank the reviewer for this detailed comment and appreciate the opportunity to\nclarify. We believe this reflects a misinterpretation of our analytical framework. We are\nnot treating temperature or precipitation as response variables, nor do we model species\nabundance or distribution using regression techniques. Our approach is as follows:\n• For each plant species, we construct an empirical climatic response curve (no\nlonger referred to as a statistical “distribution”), now based on kernel density\nestimation. These curves represent the species’ realised climatic niche,\nsummarising the temperature and precipitation conditions under which the species\nwas observed, weighted by its relative abundance.\n• For each site and climate scenario, we compare projected temperature and\nprecipitation to these species-specific curves to assess whether future conditions\nfall outside the observed climate envelope. Species are then assigned to risk\ncategories based on exceedance thresholds (e.g., >90th, >95th, >99th percentile).\n• We then aggregate the relative abundances of species falling into each risk\ncategory at a site to estimate the proportion of floral resources potentially\nvulnerable to climate change.\nThis is a threshold-based exceedance analysis, not a regression model. Our focus is on\nthe species historically foraged at each site and whether projected climatic conditions are\nlikely to exceed their empirically observed environmental ranges. This logic aligns with\nwidely used conservation risk assessment frameworks. For example, the IUCN Red List\nclassification system relies on threshold-based criteria to assess extinction risk. These\nthresholds are not derived from complex models but are applied transparently and directly\nto empirical data to inform conservation decisions. Our method follows a similar\nrationale: it uses observed climatic bounds for each species and quantifies site-level risk\nby assessing the proportion of the community likely to fall outside these bounds under\nfuture scenarios. This ensures ecological interpretability and relevance, particularly in\npolicy and conservation contexts.\nWe reworked the text to make this clearer.\nFurther, there are most likely confounders that affect pollen distribution, which are also\nnot accounted for in this analysis. In summary, the collinearity of precipitation and\ntemperature is very likely to impact the results, and separate linear regression cannot be\nfitted to the data while expecting accurate results.\n>> We agree that temperature and precipitation may covary, and that additional\nenvironmental or ecological factors likely influence pollen composition and plant\noccurrence. We are aware that ecological responses are complex and multifactorial.\nTherefore, the following resistance-resilience analysis already jointly considers\ntemperature and precipitation scenarios to evaluate potential buffering through temporal\nor spatial resource shifts. Such covariation was already noted in the manuscript, especially\nin the resilience section, and we now further emphasise that our species–climate\nassociations do not capture all potentially relevant drivers. However, our approach avoids\nthe pitfalls of multicollinearity in regression-based analyses by not relying on regression\nat all. Instead, we apply a (now non-parametric) abundance-weighted climate profiling\nmethod and a threshold-based exceedance analysis.\nWe also believe that the structure of our dataset supports this approach: our pollen samples\nare geographically extensive but spatially sparse and not rasterized. Our analysis focuses\nonly on whether a species is likely to experience climatic conditions beyond its\nempirically observed envelope at the sites where it is currently found, where we can\nassume constancy in other covariates that would need to be incorporated as variables into\npredictive modelling beyond the local sites. Attempting to model spatial presence or\nabundance, such as through SDMs, would be unwarranted given the absence of critical\ncovariates (e.g., land use, soils, agricultural practices) that would substantially affect\nextrapolation beyond our observed sites. In this context, a simple, abundance-weighted\nclimate profiling approach combined with a threshold-based exceedance analysis, both\ngrounded in empirical observations at local sites, is more defensible than attempting to\npredict beyond the data.\nWe hope these rationales and clarifications are now clearly reflected in the revised text.\n4) Food resource risk assessment\nAs noted, the calculation of distributions is inadequate, and 4) is based on 3), so that 4)\ntoo needs a considerable rework.\n>> As mentioned before, we have revised our method for estimating climatic niches of\nspecies to a kernel density estimation (KDE) to empirically model each species’ realized\nniche along temperature and precipitation gradients. KDE avoids assumptions of\nunimodality and symmetry, and allows us to better reflect skewed or multimodal\ndistributions that are common in ecological data. This is now implemented also in this\nsection.\nFurther, note that in L517-518 the authors write “the quantile was calculated”, but there\nis no such thing as “a quantile” as a quantile defined is a range, the word does not take\nmeaning until that range is defined. This is in contrast to percentiles.\n>> We appreciate the precision in the reviewer’s comment and agree that our original\nphrasing led to confusion. We have revised occurrences of this issue in the manuscript to\n“percentile thresholds” or “percentile-based risk intervals”.\nIn this part of the analysis, the authors fit a linear model and perform model selection\nusing AIC, after which statistical significance is assessed. I have two notes on this: 1)\nstepwise model selection has considerable issues that the authors should address,\ndetails on this can be found widely in the (ecological) literature. The debate is not “hot”\nanymore, and AIC model selection has largely been put on the back-burner due to its\nwidely recognized flaws, 2) model selection by information criteria cannot be combined\nwith the interpretation of statistical significance. Omitting variables by information\ncriteria naturally leads to better fitting models, and thus artificially increased statistical\nsignificance (inflated Type I error). This is commonly referred to as “p-hacking”, and is\nconsidered a statistical sin. The two paradigms: information criteria and testing\n(exploratory and confirmatory), are not unifiable.\n>> We agree with the concerns and have entirely removed the stepwise model selection.\nIn most cases, this did not change final models.\n5) Resistance and resilience\nThis calculation depends on the previous analyses with considerable flaws, and thus also\nneeds a rework. At the same time, the equation the authors provide in the text is not\nreferenced, so determined whether or not this is a sensible estimate of resistance (or\nresilience) is entirely unclear.\n>> The rework is done given the Kernel based approach now. The resistance and\nresilience calculations are based on transparent and ecologically grounded definitions that\nquantify the proportion of currently observed taxa likely to persist or be replaced under\nchanging climatic conditions. While simple in formulation, these metrics reflect common\nconservation-oriented concepts. As the equations are straightforward proportions derived\nfrom clearly stated assumptions, we believe external referencing is unnecessary, but we\nnow clarify that these are custom definitions tailored to our risk assessment framework.\nWe however included references that underline the general concepts of resistance and\nresilience.\nOverall, I conclude that the analyses in the article lack both statistical support and\necological motivation. I would have personally reanalysed the data, but due to the\nsequential nature of the pipeline, and the lack of intermediate products stored, I could not\nfind sufficient time to do the analysis.\n>> We respectfully disagree with the reviewer’s assessment. Many of the concerns raised\nappear to stem from misunderstanding our framework (exceedance analysis vs.\nregression) and earlier ambiguities in our presentation, which we have now addressed\nthrough clearer narrative structure, improved terminology, and more transparent\ndocumentation of the analytical steps. In the revised manuscript, we have clarified the\necological motivation and statistical rationale for each part of the workflow. Several of\nsuggestions of Reviewer #1 have now been implemented. None of these steps were\nquestioned by the other reviewers.\nThe full analysis pipeline is reproducible and follows open science best practices. All\nscripts and data are openly available. Intermediate outputs were in fact stored at most\nstages of the analysis, and while these files were not reloaded in subsequent scripts, we\nhave now ensured they are consistently named and clearly organised in a dedicated\nintermediate.data folder for transparency and ease of reuse.\nMy recommendation is for the authors to completely redo their analyses, as it can be done\nconsiderably more succinctly, and because I expect various (if not all) of the results to\nchange.\n>> We have restructured the entire analytical workflow, now using a kernel density\nestimation framework to derive species-specific climatic envelopes. This revision affects\nall downstream analyses. Importantly, despite this substantial methodological change, the\nmain results and interpretations remained consistent. If anything, the updated results\nindicate more pronounced risk levels than previously estimated, suggesting that our\noriginal analysis was conservative. This consistency reinforces the robustness of our\nfindings and the validity of our general approach.\nI recommend that the authors dive into the world of multivariate statistical models, as at\nleast three, if not four, of their analysis can be made redundant by fitting a single vector\nGLM (as in the mvabund R-package) or single GLLVM (as in the gllvm R-package). The\ncompositional nature of the data can be addressed by appropriately selecting the model\nstructure: the total read count should be included as an offset, and the response\ndistribution should be one of Poisson, negative-binomial or zero-inflated versions thereof.\nI would be happy to provide the authors with more details for a reanalysis if needed.\n>> We appreciate the reviewer’s suggestion to explore multivariate frameworks such as\nvector GLMs (e.g., via the mvabund package) and generalized linear latent variable\nmodels (GLLVMs). These tools are indeed powerful for specific types of inference,\nparticularly when the goal is to model multivariate responses.\nHowever, our study is conceptually distinct. We do not aim to fit multivariate species-\nenvironment regression models or predict abundance patterns. Instead, we aim to assess\nspecies-level climatic risk by estimating the likelihood of future environmental conditions\nexceeding the empirically observed climatic envelopes of currently foraged plant taxa.\nThis approach prioritizes ecological interpretability and transparency over complexity.\nFrom a statistical perspective, methods like GLLVMs require estimating large numbers\nof parameters across all species simultaneously. Given our dataset (>1,100 taxa across\n>2,000 samples), these models would face significant challenges related to computational\npower, overfitting, and identifiability, particularly for rare taxa. Moreover, they introduce\nassumptions about distributions, link functions, and latent structure that are neither\ntestable in our context nor aligned with our aims.\nBy contrast, our revised framework uses empirical, abundance-weighted kernel density\nestimates to characterise species’ realised climatic niches, followed by a percentile-based\nexceedance analysis under climate scenarios. This approach is non-parametric,\ntransparent, and grounded in observed data. It also aligns well with conservation-focused\nrisk assessment frameworks, such as the IUCN Red List, which similarly applies defined\nthreshold criteria (e.g., ≥75% loss of a species' population or range experiencing a\nspecified level of decline) to assess levels of threat. Like our method, IUCN assessments\nprioritise interpretability, data sufficiency, and operational simplicity in support of real-\nworld decision-making. For these reasons, both conceptual and statistical, we chose to\nretain our current approach, which is well aligned with our study aims and the structure\nof the data.\nFor now, I have refrained from extensive commenting on the text. That is not by lack of\ninterest – the study is interesting and I would like to see an improved version. However, I\nexpect the present version to still undergo significant changes when the analyses have\nbeen improved, so that commenting on (e.g.,) the results is not relevant at this stage.\n>> We understand the reviewer’s decision, though we note that the updated analyses,\nnow based on a non-parametric framework, do not violate any potential distributional\nassumptions questioned. The results remain stable and support the same ecological\nconclusions, now with increased statistical rigour and greater clarity in exposition.\nWe believe that part of the reviewer’s critique reflects a misalignment of methodological\nperspective. Our goal is not to model compositional variation via multivariate regression,\nbut to identify ecological risk through exceedance of empirically observed species-\nspecific climatic distributions. This is a different question, requiring a different statistical\nlens. We hope the revised manuscript now makes this distinction clear and demonstrates\nthat our methods are well suited to our data structure and ecological objectives.\nBert van der Veen, 14-03-2025\n>> Thank you for your thorought review and valuable suggestions.\nReviewer #1 (Remarks on code availability):\nThe code is available, and can be run without (much) further intervention, on which the\nauthors have done a reasonable good job.\n>> Thank you\nThe code is not well commented throughout, so that working through the code, and\nunderstanding what the authors did, took a considerable amount of time, and was\nchallenging. Ideally, this should be improved. Having said that, it would not have been\npossible to understand what the authors did in the article, without reviewing the code.\n>> Thank you. We have reworked the entire code in reply to this comment. We included\nmore comments now, also removed parts of the code that were not used anymore and may\nhave been misleading and raised concerns.\nNote that in the first \"Umbrella\" script the authors incorrectly (for anyone else accessing\nthe code) set a working directory. The sequential nature of the pipeline, without storage\nof and reference to intermediate results, makes reanalysis extremely difficult.\n>> We thank the reviewer for pointing this out. In the original code, intermediate results\nwere in fact stored, but all were placed in the base directory. We have now reorganized\nthe pipeline to include a dedicated folder (intermediate.data) for storing these files, and\nwe have improved the naming conventions for greater clarity and reproducibility. One\nmajor computational step is now omitted if intermediate results are present, and data\nloaded from file, not recalculated (was previously still executed, although intermediate\nfiles were present).\nWe also note that some steps in the pipeline, such as kernel density estimation, produce\nR objects (e.g., functions or distributions) that cannot be meaningfully serialized or\nreloaded across scripts without custom handling. For these components, the current\ndesign, where objects are passed between sequential scripts, is necessary for\nreproducibility and model integrity. We have added explanatory comments in the code to\nmake this logic clearer for external users, and declare which explicit variables are needed\ndownstream. We further specified which dependencies need to be run for individual parts,\nsince in fact not all rely on all others.\nThe pipeline would function better (for reanalysis) if each analysis was stand-alone\n>> As noted above, making each script fully stand-alone is unfortunately not feasible due\nto the sequential nature of the workflow and dependencies between steps, particularly for\ncomplex objects such as kernel density functions and preprocessed matrices. However, as\nwith any R-based workflow, users may store and reload intermediate objects (e.g., via\n.RData files) at any point. This provides flexibility for reanalysis or adaptation without\nrestructuring the analytical logic.\nReviewer #2 (Remarks to the Author):\nThe manuscript presented to me is very well written. It was a pleasure to read it, and I\nfound it difficult to point out any major shortcomings. The topic is very significant, and\nthe research is original and important. The methods are suitable, although honeybees are\nnot the best representative of bees generally, no other species would allow for such a\nwhole-season analysis and a sufficient number of samples. The Authors are aware of this\nand discuss their result accordingly.\n>> Thank you for your positive appreciation and valuable comments, which helped us\nimprove the manuscript and enhance its clarity. We also discuss differences between\nhoney bees and other bees now in more detail, also due to comments of other Reviewers.\nI only missed one piece of information about the apiaries. Did the Authors use only\nstationary apiaries for their sampling? Moving apiaries into or close proximity to mass\nflowering and honey-giving crops is a common practice and may affect the results.\n>> Yes, all apiaries were stationary, and we have now included this information in the\nMethods section. Please see lines 283-285.\nThe description of the results is very detailed, and the Authors did a good job discussing\nthe eventual shortcomings of the climate scenario results. I think it was even a bit too\ngood of a job, and part of these descriptions/interpretations of the received results should\nbe moved to the Discussion, which is short and lacking these interpretations.\n>> Thank you for your comment. You are absolutely right, the two sections are\nunbalanced, due to much of the discussion being embedded within the Results. Before\nsubmitting to Nature Communications, we initially submitted to Nature which required a\ndifferent manuscript format. During the reformatting process for Nature\nCommunications, this imbalance was introduced, and we apologize for that. Since Nature\nCommunications allows for a combined Results and Discussion section, we opted to\nmerge the two sections in order to preserve as much as possible the structure and content\nof the version already evaluated by the reviewers.\nI also found a few minor points, which could be improved. In Fig. 1A. the sampling points\nare somewhat outside the map. I understand that it is because some samples were coming\nfrom islands like Lesvos and Cyprus? Can you change the map underneath to show not\njust the EU countries but simply all the countries? It would make it easier for the reader\nto locate the most southeastern sample locations correctly.\n>> We appreciate the reviewer’s observation. Rather than expanding the map to include\nall countries, which would not resolve the issue due to projection constraints, we clarified\nthe sampling locations by explicitly naming the relevant islands in the figure caption. The\nrevised caption now reads: “The map displays the spatial distribution of pollen samples\ncollected across the 27 EU countries, including Cyprus and Malta, as well as the French\nisland of Corsica and the Greek islands of Chios, Crete, Lesvos, and Samos, with colour\nrepresenting Effective Species Richness (eH)...”\nFig2. Can you divide it into two figures? One showing crops, the other wild plants?\nAlthough I understand that it shows the total impact, it is difficult to differentiate visually\nin this form.\n>> Thank you for this valuable suggestion. We have updated Figure 2 and legend\naccordingly. The revised figure retains a unified layout but now mirrors crop and wild\ntaxa around zero: crop-related impacts are shown on the left (negative axis), and wild\nplant impacts on the right (positive axis). This preserves the total view while making the\ndistinction between crop and wild components clearer and more intuitive. Also the Figure\nwas updated using color-blind friendly colors according to journal policy.\nReviewer #2 (Remarks on code availability):\nYes, I have checked the URL however, I don't feel competent to review it.\nReviewer #3 (Remarks to the Author):\nThis is a research paper reporting findings from a 4 month study in 2023 across a very\nimpressive number and range of sites in Europe where the authors collected honeybee\npollen from aviaries at each of these sites simultaneously every 2 weeks over a 17 week\nperiod (presumably). The findings indicate the range of flowering plant species that\nhoneybees visit over May to August at 310 locations, and relate the number to\nenvironmental variables. The results provide evidence of the reduced diversity of\nflowering plants that honeybees are feeding on, and this appears to be related to increased\ntemperature and reduced rainfall. There is no data however on any negative effect on the\nhoneybees, nor any data on forage available at the time of pollen collection. Given the\nlimited range of flowering plants honeybees feed from (seemingly flagged in this paper\nthough framed a little differently), the findings might have major implications for other\nflowering species and other visiting insects, particularly in light of the fact that we know\npollinators are in decline.\n>> Thank you for your positive appreciation and valuable comments, which helped us\nimprove the manuscript and enhance its clarity.\nI would suggest a re-framing of the use of honeybees as a model, given I believe they do\nhave a preference for a small number of highly abundant flowering plants, as flagged in\nthis study.\n>> Thank you for this observation. We agree that honeybees commonly forage on a\nlimited set of highly abundant flowering plants. However, we would like to clarify that\nthis behavior reflects their generalist foraging strategy rather than a true preference in the\necological sense. In ecological terms, a \"preference\" typically implies selection for a\nresource that is not the most abundant, i.e., choosing disproportionately compared to\navailability. Honeybees, by contrast, tend to track and exploit the most abundant floral\nresources, which makes them effective indicators of overall floral availability in the\nlandscape. That said, we agree with your suggestion to reframe the notion of honeybees\nas a \"model\" in this context, and we have revised the relevant section of the manuscript\naccordingly to better reflect this nuance.\nThere is some lack of clarity and minor suggestions for changes below.\nLine 23+: Abstract is missing elements that would be useful, re the dataset, timing, to\nmake it clear that this was one year. This is modelling based on scenarios.\n>> We added a sentence in the abstract that summarises the sampling. We hope it\naddresses your comment adequately. Please see lines 26-28.\nLine 38: Some odd phrasing, e.g. how the word ecology is used in the first sentence, I also\nquestion putting value on different ecological interactions.\n>> We fully agree with this comment and thank you for pointing it out. We have revised\nthe sentence to address your concern and hope it is now clearer. Please see lines 38-39.\nLine 66: requires a citation.\n>> We added Reference 2 and 19 that support the statement.\nHoneybees are generalists but by focusing on this group you miss the deeper flowers that\nlong tongued groups depend on, and the species that have few inflorescences which I\nthink needs mentioning. Managed honeybees collect pollen which can be very easily\nlocated and collected making this a useful group to use to investigate. What’s more, your\nfindings further suggest that honeybees feed on a small number of presumbaly abundant\nflowering plants. This to me might make the unknown element, the impact of the predicted\nchanges in your analysis on flower visiting insect populations alarming as you mention\nin the discussion. But flags the importance of monitoring impacts on other species, in\naddition to an understanding of what is available at these locations at the time.\n>> We totally agree that the potential wild pollinator impact is alarming and important to\ninvestigate. We expanded and highlighted this now more in the discussion.\nLine 85: it isn’t clear to me whether latitude and temp and/or rain are correlated which\nis somewhat critical for the interpretation?\n>> We thank the reviewer for this insightful observation. In our dataset, temperature and\nlatitude are indeed correlated, as well as precipitation. We agree that this covariation is\nimportant when interpreting spatial trends in pollen diversity and climate sensitivity. In\nour view, this also reinforces the ecological motivation of the study: plant communities,\nand therefore the pollen collected by honey bees, are already naturally structured by\nclimatic gradients such as temperature and precipitation, which are themselves partly\norganised along latitude. This inherent covariation is not a confounder to be eliminated,\nbut a fundamental ecological axis of variation that underlies foraging patterns across\nEurope. We leverage this variation with our sampling design to estimate climate niches.\nTo reduce misattribution of effects, we included latitude and month of sampling in our\nlinear models rather than temperature and precipitation, which are confounded with\ngeography. In contrast, in the dissimilarity models we partitioned the individual\ncontributions of geography, temperature, and precipitation.\nLine 74: 2500 in total, simultaneously every 2 weeks from 310 apiaries, one apiary per\nlocation? May-August in 2023 is that then 8 sampling dates? Which is 2480 samples. This\ncould just clarified, and could be clearer in the intro and I think should be stated in the\nabstract at least the timing as this is important.\n>> The number of sampling rounds during the four-month period was nine, not eight.\nTherefore, the expected number of samples would be 2,790 (instead of 2,480). Of these,\nwe did not obtain data for 290 samples, either due to PCR amplification failure or samples\ncould not be collected for logistic reasons. We have added these details to both the\nIntroduction (lines 75-78) and the Methods section (lines 283-285), as requested.\nLine 206: I’m not certain these results (species richness 5-10) indicates a diverse range\nof collected pollen, some summary stats would be useful to clarify this and to clarify the\nproportions of the top flowering plant groups that honeybees are feeding on in the results.\n>> We thank the reviewer for this observation. We would like to clarify that the values\nrefer to Effective Species Richness (eH), also known as Hill number 1, which incorporates\nboth species richness and evenness of the pollen taxa, as noted in the figure legend. To\naid interpretation, we have added a supplementary figure to the GitHub repository\ndisplaying the untransformed species richness (i.e., species counts without evenness\nweighting). Further, a conservative 1% relative abundance threshold was applied at the\nsample level, based on positive control experiments, to exclude spurious low-frequency\ntaxa. This might also reduce the richness, but we consider this important to reduce the\nlikelihood of introducing contamination into the analyses. We hope this clarifies the range\nand diversity of the collected pollen.\nMay June July August\n25\n20\n15\n10\n5\n0\n40 50 60 40 50 60 40 50 60 40 50 60\nLatitude\nLine 428: why not a specific date?\n>> We added the specific dates (please see lines 283-285), as requested. Thank you for\nthe suggestion.\nLine 484: R studio is a platform for analysing data using R software, I think this needs to\nbe referred to appropriately.\n>> We changed the sentence to accommodate your comment; Please see lines 348.\nLine 482: this indicates pseudo-replication, it’s unclear whether this was accounted for\nin modelling?\n>> We thank the reviewer for raising this important concern. Indeed, climate data were\nonly available at monthly resolution per site, so all samples within a given site and month\nnecessarily share the same climate values. While this could introduce pseudo-replication,\nit is important to note that our niche estimation is based on weighted abundance data and\nis not implemented as a statistical model with replicated error structure; thus, correction\nfor pseudo-replication is not applicable at that stage. In the final linear models, data were\nre-aggregated to a single value per site and month, i.e. back to the level at which true\noriginal data exists.\nReviewer #4 (Remarks to the Author):\nThe manuscript under review, “Honey bee food resources under threat from climate\nchange,” provides valuable insights into the restoration of plant–animal interactions\namidst rapid environmental changes and increasing human interference. It effectively\nconveys the urgency of monitoring and maintaining the balance of species interactions.\nThe manuscript is generally well-written, with a clear and coherent presentation.\nHowever, I have a few concerns and suggestions that I believe could strengthen the study:\nSampling Period: The sampling timeframe appears narrow for drawing conclusions\nabout climate impacts. I recommend extending the study over at least 1–2 additional years\nto capture interannual variability.\n>> We agree that interannual variability is important when investigating long-term\necological trends. However, the inclusion of over 2,500 samples from 27 countries,\ncollected simultaneously over several months, provides robust spatial and seasonal\nssenhciR\nseicepS\nLatitude\n354045505560\nreplication that compensates for the lack of multi-year data. This design enables us to\nassess plant climatic responses and climate-driven risk patterns across a broad range of\nreal-world conditions — which is the primary aim of the study. While extending the study\nacross multiple years would certainly strengthen assessments of interannual consistency\nand climatic stability, this was beyond the scope of the current project. The study was\nfunded for one season by the European Commission (DG Environment) and relied heavily\non the motivation and coordination of a large network of citizen scientists. Given the\nunprecedented scale, geographic breadth, and standardization of the sampling effort, we\nbelieve the data are well-suited to the study’s objectives and provide meaningful,\ngeneralizable insights into climate-related risks to floral resources.\nMethodological Comments and Suggestions:\nLine 92: It would be beneficial to include a main figure showing plant abundance data\nstratified by location.\n>> We appreciate the suggestion. While we agree that spatially stratified abundance data\nare informative, the scale of the dataset, comprising over 310 locations and up to nine\nsampling rounds per site, would render such a figure overly dense and difficult to interpret\nwithin the scope of the main manuscript. Moreover, such a figure would not directly serve\nthe core research questions addressed here, i.e. climate change risks. That said, the full\nraw data per location will be made publicly available upon publication, and relative\nabundance summaries will be released by the European Commission. We are also\npreparing a follow-up manuscript focused on spatial patterns in floral resource\ncomposition, in which location-specific relative abundance data will be presented in\ngreater detail.\nLines 205–207: Could the authors comment on the occurrence of low-abundance pollen\ntaxa and whether they might result from secondary contamination?\n>> We thank the reviewer for raising this point. While secondary contamination cannot\nbe entirely ruled out, we took several precautions to minimize its impact and ensure data\nreliability. First, the bioinformatic pipeline includes multiple quality-control steps,\nincluding the removal of chimeric sequences, fungal taxa, and other artifacts. Second, a\nconservative 1% relative abundance threshold was applied at the sample level, based on\npositive control experiments, to exclude spurious low-frequency taxa. Third, for all niche-\nbased analyses, we restricted the dataset to plant species that occurred in at least ten\nindependent samples, as noted in the Methods. Finally, most of the analyses rely on\naggregated abundance values across broader risk or spatial categories, which reduces the\ninfluence of occasional low-abundance taxa. We acknowledge that some low-abundance\nrecords, though above the filtering threshold, may reflect scout foraging behavior or\nincidental collection of wind-dispersed pollen. However, given the multi-layered filtering\nand aggregation steps, we do not expect these occurrences to meaningfully affect the\nresults. For clarity and cohesion, we have not expanded this discussion in the main text,\nbut would be happy to elaborate further in the Supplementary Information if the editors\ndeem it appropriate.\nIt would be valuable to include a comparison across multiple bee species to assess\nspecies-specific preferences in pollen foraging.\n>> We agree that comparing pollen foraging across multiple bee species could yield\nvaluable insights into species-specific floral preferences. However, the focus of this study\nwas on honeybees, the most widely distributed generalist pollinator, and on assessing the\nclimate-related risks to its floral resources at local scales over an entire continent.\nIncluding other bee species, particularly wild bees, would introduce substantial variability\nin geographic distribution, behaviour, floral specialization, phenology, and habitat\nrequirements. Given the diversity of over 2,000 bee species in Europe, applying a\nstandardized sampling protocol across 27 countries would be logistically and ecologically\nunfeasible. We also believe that our spatially and seasonally broad approach would not\nbe easily transferable to other species, and thus results regarding climate-induced risks\nwould not be directly comparable. Moreover, the study already represents an\nunprecedented effort in terms of scale, coordination, and sampling detail, and was\nconstrained by time, funding, and available manpower. That said, we fully agree that\nspecies-specific comparisons represent a valuable direction for future research,\nparticularly in more localized or targeted studies that can accommodate the ecological\ndiversity of wild bees. We added more details pointed towards such as future perspectives\nin the discussion.\nLine 438: Please clarify how pollen samples were transferred into 2 ml tubes. Also,\ninclude a reference for the zirconia bead size-sorting procedure.\n>> Thank you for pointing this out. We revised the relevant sentences to incorporate your\nsuggestions. Please see lines 295-297 for details on how pollen samples were transferred\ninto 2 mL tubes, and lines 298-299 for the zirconia bead size-sorting procedure. Although\nthis information was originally described in the Guideline for Apicultural Citizen Science\nto Apply the Honey Bee Colony for Bio-monitoring of the Environment (van der Steen,\nPinto, & Quaresma et al.), we have included it here to facilitate reproducibility for readers\ninterested in following the protocol.\nLine 445: Consider placing the Maxwell® 440 RSC PureFood kit information here, and\ninclude product details in parentheses.\n>> Done (please see lines 301-304).\nLine 450: The strategy of maintaining equal DNA concentrations for each sample at the\nstart of PCR is commendable. Did this approach apply to all samples, including those\nwith low concentrations?\n>> We previously used to diluted the DNA concentration to 10ng/uL, before the 1st PCR\nstage, when working with smaller sample sets. Inadvertently, a reference to this procedure\nwas left in the Methods section when copying from earlier protocols—our apologies for\nthe oversight. In the present study, due to the large number of samples, we were unable\nto homogenize DNA concentrations. However, prior tests indicated that this lack of DNA\nconcentration dilution does not impact the results. Moreover, the paper “Standard\nMethods for Pollen Research” (https://doi.org/10.1080/00218839.2021.1948240),\nspecifically Section 8.1, does not recommend DNA dilution as a required step. We\nremoved that piece of information from the Methods section.\nLine 454: Please elaborate on how subsamples were selected for indexing. Was pooling\ninvolved, and if so, what was the basis? The “1:1 dilution of the first amplicons” requires\nfurther clarification.\n>> We have revised the sentence to improve clarity; Please see lines 312-314.\nPooling Strategy: How were pools equimolarly combined? Please provide\nmethodological detail.\n>> We added “The concentration of each pool was assessed using the Qubit Fluorometric\nQuantification method (ThermoFisher Scientific)”. Please see line 320-321.\nLine 473: The use of 97% sequence similarity—do the authors believe this threshold is\nsufficient for species-level resolution?\n>> We apologize for the lack of clarity — this should have read “at least 97% sequence\nsimilarity (top hit).” The final assignment was always based on the top hit, which in most\ncases was ≥99–100% identical. The sentence has been rephrased for accuracy. The 97%\nthreshold is a commonly applied and well-supported cutoff for ITS2-based plant\nmetabarcoding, particularly in bee-collected pollen studies. While we acknowledge that\nITS2 resolution varies among plant taxa, this threshold generally achieves good species-\nlevel discrimination for angiosperms and follows established guidelines (e.g., Keller et\nal. 2021; Sickel et al. 2015; Leonhardt et al. 2022). See below however the reply to the\nnext comment that we also found taxa where this was not valid and we therefore manually\nreduced the taxonomic assignment level.\nLines 475–477: The method used to merge undetermined taxa with Brassica is unclear.\nPlease elaborate.\n>> We apologize for the confusion caused by unclear wording. We did not merge\nundetermined taxa with Brassica. Rather, during manual verification and plausibility\nchecks, we found that several closely related Brassica species (e.g., B. napus, B. rapa,\nB. oleracea, and hybrids) could not be reliably distinguished based on ITS2 sequences due\nto insufficient divergence in the reference database. Therefore, we chose to reduce the\ntaxonomic resolution for these and treat them collectively as a Brassica species complex,\nto avoid false precision in species-level assignments. We have rephrased the manuscript\ntext accordingly for clarity (Line 338-341).\nLine 472: Which database(s) were used for taxonomic assignment?\n>> The reference databases used herein were constructed in our previous publication\n(“Semi-automated sequence curation for reliable reference datasets in ITS2 vascular plant\nDNA (meta-) barcoding”, by Quaresma et al. 2024\n(https://www.nature.com/articles/s41597-024-02962-5) and were mentioned in line 475\n(line 337-338 in the revised version).\nOverall, this manuscript addresses an important and timely topic. With these clarifications\nand additions, it could make a strong contribution to our understanding of pollinator–\nplant dynamics under climate stress.\n>> Thank you for your valuable comments, which helped us improve the manuscript and\nenhance its clarity.\nResponse to Reviewers\nReviewer #1 (Remarks to the Author):\nI agree with the authors that various of my concerns have been addressed by careful\nclarification of the authors' methodologies.\n> We appreciate your acknowledgement.\nThe context I attempted to convey in my review, which the authors found misplaced and\ncounter constructive, seems largely be lost on the authors, which is disappointing. I would\nhave liked to see a finer appreciation for the nuances of applied statistical methodologies\nand statistical best practices reflected in the authors their work, as it would have improved\nthe quality of the manuscript considerably in my perspective. However, I do not think it\ncrucial for final publication in nature communications.\n> Thank you again for this final comment. We would like to reiterate that we believe our\napplied methods are well suited to the structure and nature of the data and the study's\nresearch questions. We nevertheless appreciate the reviewer's perspective and consider\nthese comments valuable impulses for future work. We also thank the reviewer for\nacknowledgement of the manuscript's readiness for publication.\nFirst and foremost, the authors have provided carefully drafted replies to all my comments,\nwhich I appreciate. Some of my concerns the authors shared, some of which they did not,\nwhich is their prerogative and requires no further addressing. My main concern, that various\nanalyses were not clearly connected to goals or aims of the study, has been carefully\naddressed, and the manuscript carefully revised, which has greatly improved the readability\nand transparency of the study.\n> Thank you very much for this positive assessment.\nReviewer #1 (Remarks on code availability):\nThe authors have chosen to keep their pipeline in a sequential fashion, which I believe to be\na poor choice. It does not facilitate reproduction or review, albeit the authors believe their\nstudy to be in line with open science initiatives. Be that as it may, in line with open science\ninitiative does not correspond to following \"best practices\", but as the authors have\nconveyed unwillingness to change this aspect of their study, there is nothing left to improve\nin the code.\n> It might have gone unnoticed, but we reworked the code in the previous revision so that\nonly a few components remained sequential. Intermediate files are now also stored to\nfacilitate reproducibility, and overall code readability has been much improved. Only\nsections that are inherently dependent remain sequential, reflecting the sequential\ncharacter of the programming language used, and we clearly defined dependencies. We are\nuncertain why this was evaluated as not being in line with “best practices,” since the revised\ncode adheres to common standards for transparency, documentation, and reproducibility.\nReviewer #2 (Remarks to the Author):\nThe Authors did a good job improving their manuscript. All points raised were corrected or\nclarified. The manuscript now reads much better and the results are clear.\n> Thank you very much for your approval and kind evaluation.\nReviewer #3 (Remarks to the Author):\nIt looks to me like the authors have done an excellent, thorough job of addressing the\ncomments and amending the manuscript.\n> Thank you very much for your positive feedback and approval.\nReviewer #4 (Remarks to the Author):\nDear authors,\nThank you for the current revision. The manuscript is suitable for publication now. I have no\nadditional comments.\n> Thank you very much for your approval\nReviewer #4 (Remarks on code availability):\nDear authors,\nThank you for the current revision. The manuscript is suitable for publication now. I have no\nadditional comments.\n> Thank you very much for your approval"}, {"round": 2, "reviewer_comments": {"reviewer_1": "I agree with the authors that various of my concerns have been addressed by careful clarification of the authors' methodologies. The context I attempted to convey in my review, which the authors found misplaced and counter constructive, seems largely be lost on the authors, which is disappointing. I would have liked to see a finer appreciation for the nuances of applied statistical methodologies and statistical best practices reflected in the authors their work, as it would have improved the quality of the manuscript considerably in my perspective. However, I do not think it crucial for final publication in nature communications.\nFirst and foremost, the authors have provided carefully drafted replies to all my comments, which I appreciate. Some of my concerns the authors shared, some of which they did not, which is their prerogative and requires no further addressing. My main concern, that various analyses were not clearly connected to goals or aims of the study, has been carefully addressed, and the manuscript carefully revised, which has greatly improved the readability and transparency of the study.", "reviewer_2": "The Authors did a good job improving their manuscript. All points raised were corrected or clarified. The manuscript now reads much better and the results are clear.", "reviewer_3": "It looks to me like the authors have done an excellent, thorough job of addressing the comments and amending the manuscript.", "reviewer_4": "Dear authors,\nThank you for the current revision. The manuscript is suitable for publication now. I have no additional comments.\n\nDear authors,\nThank you for the current revision. The manuscript is suitable for publication now. I have no additional comments."}, "author_response": ""}]

Decision letter