Pollen meta-barcoding reveals foraging preferences of honeybees (Apis mellifera L.) along space-time gradient in Japan

The availability of pollen in urban-rural landscapes is an essential factor that influences the population dynamics of insect pollinators. The amount and diversity of pollen play a pivotal role in the foraging ecology of pollinators for their growth and health, but investigations on the spatio-temporal patterns of foraged plants remain rare, especially in cities as neo-ecosystems. Here, we explored the temporal foraging habits of a highly polylectic pollinator (Apis mellifera L.) in a study area, including different landscape classes from rural to urban areas. Mixed-pollen in each month and each location (N = 17) were analysed using DNA meta-barcoding to identify plants visited by honeybees. The results showed that the landscape class (rural, suburban and urban areas) explains spatial variations in the plant composition foraged by honeybees, but not in taxa richness. Furthermore, pollen diversity and plant composition showed a strong seasonal dependence. Furthermore, a higher plant richness and foraged woody taxa was found to occur in spring, which was mainly dominated by the genera Prunus and Acer. In summer and autumn, the genera Trifolium and Plantago of the herbaceous stratum were the most visited plants. The Fabaceae, Rosaceae, Brassicaceae, Plantaginaceae, and Onagraceae plant families were the most frequently observed in all combined samples. The present study contributes to a broader understanding of the ecology and floral preferences of honeybees, on which urban planning can rely to promote biodiversity in cities.

Introduction the non-vegetation (NDVI: from -1 to 0.199) from the vegetation (NDVI: from 0.2 to 1) 117 [61,62]. A majority filter, with a 6 × 6 filter kernel size, from the whitebox package in R [63], 118 was applied to smooth the result and aggregate regions of high uncertainty. Landscape 119 classifications at the site level were performed using demographic data (i.e. number of 120 inhabitants per admin units) and landscape metrics from the lconnect [64] and landscape 121 metrics [65] packages in R. To classify our sites along an urban-rural gradient [44,66], we 122 retained selected data: number of inhabitants per km² (dpop), the integral index of connectivity 123 (IIC) [67], the effective mesh size (MESH) [68], Shannon's evenness index (SHEI) [69], the 124 vegetation cover proportion (veg cover), vegetation patch density (Pd) [70], and the median 125 vegetation class NDVI (NDVI median) ( Table 1). We conducted principal component analysis (PCA) of the landscape dataset to visualise the 130 differences among our study sites. The unsupervised k-means clustering method was applied to 131 delineate the landscape category along the urban-rural gradient into k groups. Before initiating 132 the analysis, the data were standardised using the scale function in R to make variables 133 comparable. As a result, the clustering algorithm was independent of any variable unit. The 134 number of k groups required to be defined as the first step was determined using the elbow 135 method [71]. The k-means partitioning analysis was performed using the k-means function with 136 25 random sets [72] and the factoextra package in R for PCA graphical representations [73].

165
Data processing 166 We used "FASTX Barcode Splitter" from Fastx toolkit, a short-reads pre-processing tool, to 167 extract only the Miseq reads with sequences readings that matched the primers used exactly 168 [75]. Next, the reads were denoised and filtered using Sickle software [76]  classifier was used for taxonomic assignment with the output sequence [79,80]. Following the taxonomy classification step, taxa-abundance data and operational taxonomic unit (OTU) data 177 were applied to the R environment [55]. First, the assignment of all OTUs below the identity 178 threshold of 97% was discarded [46,81]. Next, the number of reads was sorted by genus and 179 sample (i.e. site and collection date), and was then expressed as the ratio between the read count 180 and the sum of read number per sample for each genus. Genera accounting for less than 0.05% 181 of the total number of readings for a single sample were excluded to prevent false positives.

182
Moreover, a sample was removed if it accounted for less than 1000 reads to limit inferences 183 from insufficient sequencing depth [82].

184
Taxonomic analysis 185 Given the poor correlation between the actual pollen proportion per foraged plant taxa and the 186 sequencing read proportions from one locus [36,41], the analyses undertaken were based 187 exclusively on the incidence (i.e. presence/absence binary arrays) dataset. This process prevents 188 the misinterpretation of DNA metabarcoding data [83]. The richness of the pollen samples (i.e.

225
Landscape classification 226 The method of differentiating vegetation from impervious surfaces using NDVI provided 227 convincing results after crosschecking, even in the complex environment of an urban matrix.

228
The two first dimensions of the landscape PCA from landscape variables of our sampling sites 229 described a high percentage of the variance (axis 1 = 81.6% and axis 2 = 11.0%; Fig 2).

230
According to the elbow method ( Fig S2) of the k-means partitioning, we classified our sampling August, and September) (Fig 3). foraged in the urbanised landscape than in the countryside (Fig S4).

298
Effect of sampling period 299 The proportion of foraged woody taxa decreased significantly over the seasons (G = 87.5, p < 300 0.001), with a peak of 46% in April and a low of 10% in September. Over the sampling months, 301 honeybees foraged mainly on alien plant taxa ( Fig 5B). Cultivar taxa were more foraged in noticeable phenological turnover in the pollen composition ( Fig 5D) was observed, with May 306 serving as a transition bridge. This has already been highlighted by the discontinuities in the 307 NMDS ordination (Fig 4). Following this shift, the genus Trifolium spp. was highly dominant 308 in the June and July samples. In addition, the herbaceous genera Plantago and Oenothera spp. Trifolium spp., Rosa spp., and Allium spp. were the only genera that were observed throughout 316 the study period. issues on the landscape variables by increasing the foraging area overlaps among the sites [94].

343
This was not tested in the present study because we considered that there was no spatial 344 influence on foraging preferences between each sampling colony due to the high variations in 345 their foraging behaviour [95,96], even at the local scale for colonies of the same apiary [97]. in its infancy, particularly in the data processing of imagery classification by deep learning 357 [100]. However, despite these spatial limitations, the approach led to satisfactory classification, 358 which paves the way for further investigations.

359
The use of DNA meta-barcoding allowed for the identification of the pollen composition of the 360 samples. This methods yields to the identification of 307 taxa, which is higher than previous 361 studies owing to a greater sampling effort among sites and seasons [46,47,52]. This reflects a 362 great diversity of plant resources and highlights foraging patterns, regardless of the landscape.

363
Through our molecular analysis, we gauged potential biases that may arise from PCR 364 amplification, sequencing errors, or arbitrary choices in the bioinformatics pipeline [101,102].
a metagenomics approach against light microscopy. They recommended using a multi-locus 367 meta-barcoding strategy in combining plastid loci, such as the matK or trnL markers with ITS2, 368 to improve the correlation between the number of reads and classic palynology identification 369 [103,104]. However, in the present study, we used one pair of ITS1 primers [74], which led to 370 the consideration of all the generated datasets on their incidences due to the poor correlations 371 between the number of sequence reads and visually identified pollen grains [36,41,105]. insects.

377
The landscape type did not influence forage plant richness, as reported in previous studies 378 [49,51]. However, considering the spatial variation of taxa composition among the sites of our 379 area of interest (i.e. beta diversity) [106], our results showed a spatial structure of foraged plant 380 communities in countryside, suburban, and urban environments (i.e. by merging the urban and 381 urban centre areas) (Fig 4). It is likely that honeybee colonies modify their foraging preferences 382 due to the high prevalence of unattractive ornamental flowers in urban landscapes [14]. The 383 urban matrix also offers smaller spread patches and less dense floral resources, which contribute 384 to the foraging change of workers [47,97]. This shift in the prospected flora is also coupled to 385 fulfil the nutritional demand with a diverse and complementary floral diet [107,108]. Therefore, 386 it shows the importance of taking beta diversity and not only the local richness into account to 387 understand the community structure of foraged plants throughout space and time scales to 388 protect local biodiversity and support conservation planning [53,109]. Despite the landscape 389 structure of the foraged plant community, the trait-based analysis revealed no significant 390 preference for honeybee colonies, except for herbaceous plants in rural landscapes (Fig 5A).
In all the samples, 35 plant taxa were shared throughout all the landscapes studied, 392 corresponding to 45% of all read counts. The top three plant families were Fabaceae, Rosaceae, 393 and Brassicaceae in the samples (Fig S4). The presence of grasses (Poaceae) may seem 394 surprising among the most frequent families in the samples (Fig S4), given that these plants are 395 considered unsuitable for A. mellifera resource needs [110]. In view of their dominance, it 396 seems unlikely that anemophilic pollen is the result of contamination by pollen blown from 397 flowers onto the body of the bees, such as the rice paddy field in Japan [111]. Other recent 398 studies have shown that pollinators (i.e. bees and syrphids), particularly Apis bees, interact with 399 wind-pollinated plant species for their nutrient or nesting requirements [112]. This means that 400 these grass species could be implemented in rural and urban management recommendations.

401
In contrast to the minor variations in landscape diversity, we observed a strong seasonal effect 402 on plant richness (Fig 3), the foraged plant community, and traits (Fig 4 and 5B). We observed 403 higher plant richness and foraged woody taxa in spring than in the other two seasons (Fig 3 and   404 5B), as reported previously [82,113,114]. The genera Prunus spp. and Acer spp. (Fig 5D) 405 dominate the foraged woody stratum during this season, as these taxa constitute an important 406 source of pollen for the early season [82,103,114,115]. Indeed, bee breads with high proportions 407 of both genera were positively correlated with a higher protein content [107]. Particularly, the 408 complexity and high range of foraged plants is known to be beneficial to the "nutritional value" 409 of these bee bread stocks, and thus to honeybee immunity [28]. After the spring period, the 410 proportion of foraged herbaceous strata gradually substituted woody taxa to reach 411 approximately 90% of herbaceous foraged taxa in September, in agreement with previous 412 studies [82,114]. This growing herbaceous stratum is mainly dominated by Trifolium spp. and

413
Plantago spp. which, because of their long flowering period or floral consistency, explains the 414 lower richness of taxa foraged in summer and autumn [114,116,117]. Moreover, clover species 415 (i.e. Trifolium spp.) are highly ubiquitous in grasslands, such as meadows for rural areas or parks and gardens for urban areas [114,118] and may contribute to the concentration effect of 417 amino acid content [107]. Concerning the temporal dynamics of the native traits (Fig 5B) August by a rupture of the dominant flower prevalence, which may correspond to the seasonal 424 dearth of floral resources (Fig 4 and 5B) [114,121]. From the perspective of the honeybees, 425 they mitigate this effect by increasing their foraging range, requiring extra effort for sometimes 426 worthless rewards [121]. Therefore, in an urban greening context, it would be relevant to ensure 427 the diversity and abundance of floral resources near apiaries by targeting seasonal gaps in pollen 428 collection [11]. Finally, the study of other co-variables, such as brood monitoring or estimating 429 the pollen collection/reserve of each colony, could be used to compare the conditions of each 430 sampled colony and improve our understanding of the foraging patterns of the colonies [122].

431
The availability of floral resources needs to be sufficient to host both domesticated honeybees 432 and local wild pollinators. The percentage of impervious surfaces plays a major role in the 433 pollination of biodiversity [123]. Therefore, populations of managed honeybees must be 434 regulated to ensure wild pollinator populations are not adversely affected [18,124] Table S1. Climatic data of the study sites