Floral hosts of leaf-cutter bees (Megachilidae) in a biodiversity hotspot revealed by pollen DNA metabarcoding of historic specimens

prior to further analyses. Each plant family identified by RDP classifier was treated as an opera- tional taxonomic unit (OTU). Data were handled in two different ways. Biases during DNA extraction of pollen and DNA amplification may lead to a skew in the abundances of certain taxa being identified [2, 31]. To circumvent this, plant family data were converted to absence/ presence counts for each plant family for each sample in a biome, referred to as detection counts. These counts were summed across samples for each plant family to obtain plant family prevalence within a biome across all samples. Data were also analysed as read counts per sam- ple. OTU tables were generated for read and detection counts for family abundances. The Sickel et al. (2015) ITS2 database does not represent South African plant diversity ade- quately to make species-level classifications [26]. In order to get an estimation of angiosperm species abundance ranges in different bee species, sequence reads from each sample were clus- tered into OTUs using the cluster_otus command in usearch v8.0.1517 [32]. All reads per sam- ple were then assigned to OTUs using usearch_global alignment, and taxonomy added to the OTU sequences with the utax command using the utax-compatible ITS2 database from [3]. During the taxonomy assignment, confidence scores were not assigned, and a raw score cut- off of 10 was chosen for a species assignment to be included. OTU tables were again created for both read counts and detection counts for species abundances, disregarding species not

from angiosperm families. Rarefaction curves were drawn for plant family and species assign- ments for all bee species using vegan v. 2.3.4 [33] in R v. 3.2.4 [34]. Taxonomic assignments were checked for local family occurrence against the Plants of Southern Africa (POSA) data- base [35]. For the Succulent Karoo biome, the search was confined to the Western Cape, Northern Cape and Cape Floristic Region (as defined in the database), for the Savanna biome the remaining regions were selected as search criteria and for the widespread species all the South African provinces were selected. The POSA database, the vegetation map of South Africa [36], and a quarter degree grid square (QDGS) shape file were used to calculate the fre- quency of occurrences per family within each biome in South Africa. Spatial analyses were car- ried out using the rgdal package v.1.2–4 in R v. 3.2.4. A list of recorded flower visits (from field observation) was also generated for each bee species using available literature [37–39] and the Catalogue of Afrotropical Bees, accessed through the online Global Biodiversity Information Facility (GBIF, www.gbif.org/dataset/da38f103-4410-43d1-b716-ea6b1b92bbac [40]) for comparison. To compare bee taxa and to assess the effect of time since collection on the mean number of taxa detectable in pollen loads, generalised linear models that incorporated a Poisson distribu- tion and log link function were used in SPSS 23.0 (IBM Corp., Armonk, New York). Separate models were run with angiosperm species and families as response variables. To account for

statistical non-independence among bees sampled from the same localities, collection locality was treated as a subject in generalised estimating equations that used an exchangeable correla- tion matrix. Bee species was treated as a fixed factor and time was treated as a covariate. Signif- icance was assessed using Wald statistics, and post-hoc comparisons among means were carried out using the Dunn-Sidak procedure. Marginal (model-adjusted) means were obtained by back-transformation from the log scale, which also results in asymmetrical standard errors. To assess whether bee species carried different assemblages of pollen, a similarity matrix for square-root transformed data on detection counts were calculated using the Bray–Curtis method and then plotted in two-dimensions with non-metric multi-dimensional scaling (NMDS) using Past 3.14 [41]. The significance of differences in pollen assemblages among bee species was assessed using ANOSIM, a non-parametric permutation procedure based on the similarity matrix underlying the ordination. Observed R-values were compared with the distri- bution of R-values generated by up to 10,000 random permutations of the sample labels in order to assess statistical significance. PLOS ONE Floral hosts of Megachilidae revealed by metabarcoding PLOS ONE | https://doi.org/10.1371/journal.pone.0244973 January 21, 2021 5 / 16