Australian Stingless Bee Nest Microbiome

Supplementary material for the paper 16S Amplicon Metabarcoding from Nest Materials of Native Australian Stingless Bees

Authors

Affiliation

Boyd Tarlinton

School of Biology and Environmental Science, QUT

Flavia Carmelina Massaro

School of Biology and Environmental Science, QUT

Caroline Hauxwell

School of Biology and Environmental Science, QUT

Published

November 10, 2022

Introduction

This repository contains the code and documentation required to reproduce the results presented in our paper “16S Amplicon Metabarcoding from Nest Materials of Native Australian Stingless Bees”. We describe how the raw metabarcoding datasets were processed using QIIME 2 and further analysed in R to produce the figure included in our paper. Beta diversity statistics describing microbiome variability associated with the different bee species and material types sampled are calculated, visualised, and used for hypothesis testing.

Data Availability

All data are available from NCBI SRA with study accession SRP405832.

Samples were sequenced in two batches, identified by the last two characters of each sample name.

Processing in `QIIME 2`

Raw sequence data were processed on the QUT HPC Lyra. Batches B1 and B2 were analysed separately using the following scripts: B1, B2.

Setup

Load packages and pre-defined functions.

library(phyloseq)
library(qiime2R)
library(tidyverse)
library(magrittr)
library(symbioteR)
library(vegan)
library(ggforce)
library(hrbrthemes)
library(concaveman)
library(ggh4x)
library(ggtext)
library(svglite)
source("functions.R")

Load Data

We read the .QZA artifacts as phyloseq objects, merge them, and add metadata.

hiveMB_1 <- readQZAFolder("Data/H1")
hiveMB_2 <- readQZAFolder("Data/H2")

hiveMB <- merge_phyloseq(hiveMB_1, hiveMB_2)

sample_data(hiveMB) <- read_tsv("Data/Metadata.tsv") |>
  mutate(Year = str_extract(Date_Harvest, "\\d{4}$")) |> 
  column_to_rownames("ID")

Filter Samples

We remove any mitochondrial or chloroplast ASVs, and any that were unclassified at the kingdom level.

hiveMB <- removeNonbacterial(hiveMB)

Beta Diversity

Perform Ordination

Convert to relative abundance and perform NDMS ordination of Bray-Curtis dissimilarities.

set.seed(42)

hiveMBRel <- hiveMB |> 
  makeObservationsRelative()

ord <- ordinate(hiveMBRel, 
                distance = "bray", 
                method = "NMDS", 
                k = 4, 
                maxit = 1000,
                trace = 0)

ord

Call: metaMDS(comm = veganifyOTU(physeq), distance = distance, k = 4, trace = 0, maxit = 1000)

global Multidimensional Scaling using monoMDS

Data: veganifyOTU(physeq) Distance: bray

Dimensions: 4 Stress: 0.08892112 Stress type 1, weak ties Best solution was repeated 4 times in 20 tries The best solution was from try 13 (random start) Scaling: centring, PC rotation, halfchange scaling Species: expanded scores based on ‘veganifyOTU(physeq)’

Plot Ordination

Plot the NMDS ordination.

custPal <- viridis::turbo(n = 4)[c(1:2, 4)]

ordPlot <- plot_ordination(hiveMBRel, ord, justDF = TRUE) |> 
    ggplot(data = _, 
           aes(x = NMDS2,
               y = NMDS1,
               shape = Material, 
               color = Bees)) +
    geom_point(size = 3, alpha = 0.8) +
    theme_ipsum(base_family = "sans") +
    geom_mark_hull(aes(x = NMDS2, 
                       y = NMDS1,
                       x0 = NMDS2,
                       y0 = NMDS1,
                       fill = Year, 
                       label = Year),
                   inherit.aes = FALSE,
                   label.margin = margin(),
                   con.cap = 0,
                   label.fill = NA,
                   concavity = 3,
                   alpha = 0.2) +
    scale_colour_manual(values = custPal, "Species",
                        labels = c(expression(paste(italic("A. australis"))),
                                   expression(paste(italic("T. carbonaria"))),
                                   expression(paste(italic("T. hockingsi"))))) +
    scale_fill_brewer(palette = "Pastel2", guide = "none") +
    scale_shape_manual(values = c(15, 16, 17, 18)) +
    theme(legend.text.align = 0,
          legend.position = "bottom",
          legend.box = "vertical",
          legend.margin = margin(),
          text = element_text(family = "sans"),
          aspect.ratio = 3.5/4.5) +
    scale_x_continuous(expand = expansion(add = 0.5)) +
    scale_y_continuous(expand = expansion(add = 0.5))

ggsave("Out/Figures/Figure1.svg", ordPlot, device = "svg", bg = "white", 
       width = 6.5, height = 6, dpi = 1200)

Significance Test

Perform a PERMANOVA test.

dist <- phyloseq::distance(hiveMB, method = "bray")
adonis2(dist ~ Bees + Material + Year + Hive_ID, 
                as(sample_data(hiveMB), "data.frame"))

Table 1: PERMANOVA Results
	Df	SumOfSqs	R2	F	Pr(>F)
Bees	2	4.358568	0.1770692	7.916509	0.001
Material	3	3.857593	0.1567168	4.671055	0.001
Year	2	1.071197	0.0435180	1.945625	0.004
Hive_ID	10	3.215226	0.1306203	1.167969	0.095
Residual	44	12.112472	0.4920758	NA	NA
Total	61	24.615056	1.0000000	NA	NA

Taxonomic Abundance

Key Genera

We highlight the abundance of a number of key taxonomic groups referenced in the paper.

abPlot <- hiveMBRel |> 
  tax_glom("Genus") |> 
  psmelt() |> 
  group_by(Genus, Material, Bees) |> 
  mutate(Bees = case_when(Bees == "Austroplebia australis" ~ "A. australis",
                          Bees == "Tetragonula carbonaria" ~ "T. carbonaria",
                          Bees == "Tetragonula hockingsi" ~ "T. hockingsi")) |> 
  summarise(Proportion = mean(Abundance)) |> 
  filter(Genus %in% c("Lactobacillus", "Bombella", "Gilliamella",
                      "Snodgrassella")) |> 
  ggplot(aes(x = Material, y = Genus, fill = Proportion, 
              label = round(Proportion, 2))) +
    geom_tile() +
    geom_text(colour = "white") +
    facet_wrap(~ Bees) +
    theme_ipsum(base_family = "sans") +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
          strip.text = element_text(face = "italic"),
          legend.position = "none")

ggsave("Out/Figures/Figure2.png", abPlot, device = "png", bg = "white", 
       width = 8, height = 4, dpi = 1200)

Gammaproteobacteria

We note the abundance of ASVs from the class Gammaproteobacteria.

hiveMBRel |> 
  tax_glom("Class") |> 
  psmelt() |> 
  group_by(Bees, Material, Class) |>  
  summarise(Proportion = round(mean(Abundance), 2)) |> 
  filter(Class == "Gammaproteobacteria") |> 
  select(!Class)

Table 2: Mean Proportion of Gammaproteobacteria
Bees	Material	Proportion
Austroplebia australis	Food	0.00
Austroplebia australis	Honey	0.50
Austroplebia australis	Pollen	0.53
Austroplebia australis	Propolis	0.00
Tetragonula carbonaria	Food	0.00
Tetragonula carbonaria	Honey	0.11
Tetragonula carbonaria	Pollen	0.04
Tetragonula carbonaria	Propolis	0.00
Tetragonula hockingsi	Food	0.00
Tetragonula hockingsi	Honey	0.26
Tetragonula hockingsi	Pollen	0.10
Tetragonula hockingsi	Propolis	0.00

Export Data

We provide a copy of the phyloseq object that can be easily imported for further analysis in R.

saveRDS(hiveMB, "Out/hiveMB.rds")

Download it here.

Session Info

sessionInfo()

R version 4.2.2 (2022-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19044)

Matrix products: default

locale:
[1] LC_COLLATE=English_Australia.utf8  LC_CTYPE=English_Australia.utf8   
[3] LC_MONETARY=English_Australia.utf8 LC_NUMERIC=C                      
[5] LC_TIME=English_Australia.utf8    

attached base packages:
[1] stats     graphics  grDevices datasets  utils     methods   base     

other attached packages:
 [1] svglite_2.1.0        ggtext_0.1.2         ggh4x_0.2.2         
 [4] concaveman_1.1.0     hrbrthemes_0.8.0     ggforce_0.4.1.9000  
 [7] vegan_2.6-4          lattice_0.20-45      permute_0.9-7       
[10] symbioteR_0.0.0.9000 magrittr_2.0.3       forcats_0.5.2       
[13] stringr_1.4.1        dplyr_1.0.10         purrr_0.3.5         
[16] readr_2.1.3          tidyr_1.2.1          tibble_3.1.8        
[19] ggplot2_3.4.0        tidyverse_1.3.2      qiime2R_0.99.6      
[22] phyloseq_1.42.0     

loaded via a namespace (and not attached):
  [1] readxl_1.4.1           backports_1.4.1        Hmisc_4.7-1           
  [4] systemfonts_1.0.4      plyr_1.8.7             igraph_1.3.5          
  [7] splines_4.2.2          GenomeInfoDb_1.34.2    digest_0.6.30         
 [10] foreach_1.5.2          htmltools_0.5.3        viridis_0.6.2         
 [13] qpcR_1.4-1             fansi_1.0.3            checkmate_2.1.0       
 [16] googlesheets4_1.0.1    cluster_2.1.4          tzdb_0.3.0            
 [19] Biostrings_2.66.0      modelr_0.1.9           extrafont_0.18        
 [22] vroom_1.6.0            extrafontdb_1.0        timechange_0.1.1      
 [25] jpeg_0.1-9             colorspace_2.0-3       rvest_1.0.3           
 [28] textshaping_0.3.6      haven_2.5.1            xfun_0.34             
 [31] crayon_1.5.2           RCurl_1.98-1.9         jsonlite_1.8.3        
 [34] survival_3.4-0         iterators_1.0.14       ape_5.6-2             
 [37] glue_1.6.2             polyclip_1.10-4        gtable_0.3.1          
 [40] gargle_1.2.1           zlibbioc_1.44.0        XVector_0.38.0        
 [43] V8_4.2.2               Rhdf5lib_1.20.0        Rttf2pt1_1.3.11       
 [46] BiocGenerics_0.44.0    DEoptimR_1.0-11        scales_1.2.1          
 [49] DBI_1.1.3              Rcpp_1.0.9             viridisLite_0.4.1     
 [52] gridtext_0.1.5         htmlTable_2.4.1        bit_4.0.4             
 [55] foreign_0.8-83         Formula_1.2-4          stats4_4.2.2          
 [58] DT_0.26                truncnorm_1.0-8        htmlwidgets_1.5.4     
 [61] httr_1.4.4             RColorBrewer_1.1-3     ellipsis_0.3.2        
 [64] pkgconfig_2.0.3        NADA_1.6-1.1           farver_2.1.1          
 [67] nnet_7.3-18            dbplyr_2.2.1           deldir_1.0-6          
 [70] utf8_1.2.2             labeling_0.4.2         tidyselect_1.2.0      
 [73] rlang_1.0.6            reshape2_1.4.4         munsell_0.5.0         
 [76] cellranger_1.1.0       tools_4.2.2            cli_3.4.1             
 [79] generics_0.1.3         ade4_1.7-20            broom_1.0.1           
 [82] evaluate_0.18          biomformat_1.26.0      fastmap_1.1.0         
 [85] ragg_1.2.4             yaml_2.3.6             bit64_4.0.5           
 [88] knitr_1.40             fs_1.5.2               robustbase_0.95-0     
 [91] rgl_0.110.2            nlme_3.1-160           xml2_1.3.3            
 [94] compiler_4.2.2         rstudioapi_0.14        curl_4.3.3            
 [97] png_0.1-7              zCompositions_1.4.0-1  reprex_2.0.2          
[100] tweenr_2.0.2           stringi_1.7.8          highr_0.9             
[103] gdtools_0.2.4          Matrix_1.5-1           multtest_2.54.0       
[106] vctrs_0.5.0            pillar_1.8.1           lifecycle_1.0.3       
[109] rhdf5filters_1.10.0    data.table_1.14.4      bitops_1.0-7          
[112] R6_2.5.1               latticeExtra_0.6-30    renv_0.16.0           
[115] gridExtra_2.3          IRanges_2.32.0         codetools_0.2-18      
[118] MASS_7.3-58.1          assertthat_0.2.1       rhdf5_2.42.0          
[121] minpack.lm_1.2-2       withr_2.5.0            S4Vectors_0.36.0      
[124] GenomeInfoDbData_1.2.9 mgcv_1.8-41            parallel_4.2.2        
[127] hms_1.1.2              grid_4.2.2             rpart_4.1.19          
[130] rmarkdown_2.17         googledrive_2.0.0      Biobase_2.58.0        
[133] lubridate_1.9.0        base64enc_0.1-3        interp_1.1-3

Introduction

Data Availability

Processing in QIIME 2

Setup

Load Data

Filter Samples

Beta Diversity

Perform Ordination

Plot Ordination

Significance Test

Taxonomic Abundance

Key Genera

Gammaproteobacteria

Export Data

Session Info

Processing in `QIIME 2`