Research article
Diesel exhaust particles alter gut microbiome and gene expression in the
bumblebee Bombus terrestris
Dimitri Seidenath1, Alfons R. Weig2,
Andreas Mittereder3, Thomas
Hillenbrand3, Dieter Brüggemann3,
Thorsten Opel4, Nico Langhof4,
Marcel Riedl1, Heike Feldhaar1*,
Oliver Otti1,5
1Animal Ecology I, Bayreuth Center of Ecology and
Environmental Research (BayCEER), University of Bayreuth,
Universitätsstrasse 30, 95440 Bayreuth, Germany
2Keylab Genomics and Bioinformatics, Bayreuth Center
of Ecology and Environmental Research (BayCEER), University of Bayreuth,
Universitätsstrasse 30, 95440 Bayreuth, Germany
3Department of Engineering Thermodynamics and
Transport Processes, University of Bayreuth, Germany,
Universitätsstrasse 30, 95440 Bayreuth, Germany
4Department of Ceramic Materials Engineering,
University of Bayreuth, Prof.-Rüdiger-Bormann-Str. 1, 95447 Bayreuth,
Germany
5Applied Zoology, TU Dresden, Zellescher Weg 20b,
01062 Dresden, Germany
ORCID ID Heike Feldhaar 0000-0001-6797-5126
ORCID ID Oliver Otti 0000-0002-2361-9661
ORCID ID Alfons R. Weig 0000-0001-8712-7060
*Corresponding author : Heike Feldhaar
Phone: +49921552645, e-mail:
feldhaar@uni-bayreuth.de
Keywords : transcriptome, air pollution, particulate matter,
insect decline, brake dust, pollinator
Abstract
Insect decline is a major threat for ecosystems around the world as they
provide many important functions, such as pollination or pest control.
Pollution is one of the main reasons for the decline, besides changes in
land use, global warming, and invasive species. While negative impacts
of pesticides are well studied, there is still a lack of knowledge about
the effects of other anthropogenic pollutants, such as airborne
particulate matter, on insects. To address this, we exposed workers of
the bumblebee Bombus terrestris to sublethal doses of diesel
exhaust particles (DEPs) and brake dust, orally or via air. After seven
days, we looked at the composition of the gut microbiome and tracked
changes in gene expression. While there were no changes in the other
treatments, oral DEP exposure significantly altered the structure of the
gut microbiome. In particular, the core bacterium Snodgrassellahad a decreased abundance in the DEP treatment. Similarly, transcriptome
analysis revealed changes in gene expression after oral DEP exposure,
but not in the other treatments. The changes are related to metabolism
and signal transduction which indicates a general stress response. Taken
together, our results suggest potential health effects of DEP exposure
on insects, here shown in bumblebees, as gut dysbiosis may increase the
susceptibility of bumblebees to pathogens, while a general stress
response may lower available energy resources. However, experiments with
multiple stressors and on colony level are needed to provide a more
comprehensive understanding of the impact of DEPs on insects.
Introduction
Global biodiversity loss is one of the major challenges humanity
currently faces (Diaz et al. 2006, Dirzo et al. 2014). Especially the
rapid decline in insects is cause for concern, as they provide or
contribute to many important ecosystem functions such as pollination,
nutrient cycling, pest control, and linking trophic levels (Cardoso et
al. 2020, Noriega et al. 2018). Pollution is one of the major reasons
for the decline besides intensification of land use, climate change, and
invasive species (Milicic et al. 2021, Sanchez-Bayo & Wyckhuys 2019).
Pesticides harm insects on many different levels ranging from subtle
changes in the gut microbiome over behavioral changes to increased
mortality (Desneux et al. 2007, Motta et al. 2018, Ndakidemi et al.
2016). Other anthropogenic pollutants might also contribute to the
observed declines in insects, but their impacts are often less well
studied (Cameron & Sadd 2020, Feldhaar & Otti 2020, Sanchez-Bayo &
Wyckhuys 2019). Airborne particulate matter deriving from traffic or
industrial processes has become ubiquitous in the environment (Gieré &
Querol 2010, Zereini & Wisemann 2010). While the harmful effects on
mammals, in particular humans, have been intensively studied, research
investigating the impact on insects remains scarce (Kim et al. 2015,
Valavanidis et al. 2008). Insects can encounter these pollutants in
various ways, e.g. by foraging in contaminated areas, consuming
contaminated food or direct deposition on the insect’s cuticle (Feldhaar
& Otti 2020, Lukowski et al. 2018, Negri et al. 2015). The airborne
particulate matter might enter an insect´s body via oral ingestion or
the tracheal system (Feldhaar & Otti 2020, Negri et al. 2015). Social
insects might be at an increased risk, as pollutants are transferred to
and stored in their nests, which could lead to a higher exposure to
conspecifics and the brood (Feldhaar & Otti 2020, Hladun et al. 2016).
Vehicle brake dust and diesel exhaust particles (DEPs) are major classes
of airborne particulate matter deriving from traffic released into the
environment (Hamilton & Hartnet 2013, Harrison et al. 2012, Rönkkö &
Timonen 2019). Brake dust particles contain various metals and phenolic
compounds, depending on the brake lining used (Iijima et al. 2007,
Thorpe & Harrison 2008). Exposure of different invertebrate species to
such particles showed mixed effects. Particulate matter contamination in
soil did not affect colony founding in the ant Lasius niger(Seidenath et al. 2021). However, soil-feeding earthworms
(Eisenia fetida ) showed a strongly increased mortality when
exposed to soil spiked with brake dust particles (Holzinger et al.
2022). DEPs have a different composition than brake dust. They are
composed of an elemental carbon core with adsorbed organic compounds,
such as polycyclic aromatic hydrocarbons (PAHs), and traces of metals
and other elements (Greim 2019, Wichmann 2007). Exposure to high doses
of diesel exhaust particles (up to a concentration of 2 g/L) in food
over a period of seven days reduced survival in Bombus terrestrisworkers compared to controls (Hüftlein et al. unpublished).
Many classical ecotoxicology approaches focus on the effect of a
substance on mortality, growth, or reproduction. However, pollutants can
also have more subtle sublethal effects on insects which may have severe
consequences in the long-term (Straub et al. 2020). Direct sublethal
effects include changes in physiology such as stress reactions or
detoxification processes. By interacting with microorganisms inside the
insect’s body, oral exposure to pollutants may indirectly affect insect
health.
Most eukaryotic organisms and their associated microbes form an entity,
the so-called holobiont (Theis et al. 2016, Zilber-Rosenberg &
Rosenberg 2008). In insects, microorganisms can be found in the
digestive tract, the exoskeleton, the hemocoel, or within cells (Douglas
2015). The insect gut microbiome has a range of functions that include
protection from pathogens, detoxification, digestion, and the production
of essential nutrients (Engel & Moran 2013). Social bumblebees
(Bombus spp.) and honeybees (Apis mellifera ) are model
organisms to study gut microbiota as their gut microbiome is rather
simple and highly conserved (Engel et al. 2016, Kwong & Moran 2016,
Zhang & Zheng 2022). A few core bacterial taxa dominate the gut
microbiome of bumblebees: Snodgrassella , Gilliamella ,Schmidhempelia , Bifidobacteriaceae (Bifidobacterium and
Bombiscardovia ) and two clusters within Lactobacillaceae (Hammer et al.
2021, Koch & Schmid-Hempel 2011a, Martinson et al. 2011). While many
functions of the bacterial symbionts in bumblebees have been proposed,
only very few have been demonstrated in experiments (Hammer 2021, Zhang
& Zheng 2022). Resistance to the common trypanosomatid parasiteCrithidia bombi is higher in bumblebees with an intact microbiome
compared to microbiota-free individuals (Koch & Schmid-Hempel 2011b).
Moreover, infection outcomes of C. bombi vary with host
microbiota composition rather than genotype (Koch & Schmid-Hempel
2012). The gut microbiome of bumblebees may be important for
detoxification as microbiota-free individuals had lower survival when
exposed to toxic concentrations of selenate (Rothman et al. 2019).
Examining the effects of anthropogenic pollutants, such as airborne
particulate matter, on the gut microbiome is an important tool for
assessing their risk for insect health (Duperron et al. 2020). Even with
a conserved gut microbiome, the relative abundance of core bacteria and
the presence of other microorganisms will vary with age, diet and
changing environmental parameters (Kwong & Moran 2016, Koch et al.
2012). Different pollutants affect the microbial composition of bee
guts. In honeybee workers, pesticides or antibiotics change the relative
and absolute abundance of core gut microbiota species (DeGrandi-Hoffmann
et al. 2017, Motta et al. 2018, Raymann et al. 2017). An array of
environmental toxicants, such as cadmium, copper, selenate, and hydrogen
peroxide, alter the gut microbiome of Bombus impatiens at
field-realistic concentrations (Rothman et al. 2020). These shifts in
the microbial community may affect bumblebee health. Intestinal
dysbiosis, compositional and functional alteration of the microbiome, is
associated with various diseases and health problems in humans and
vertebrates (De Gruttola et al. 2016, Levy et al. 2017, Shreiner et al.
2015). In insects, dysbiosis negatively affects reproductive fitness,
immunity, and resistance to pathogens (Ami et al. 2010, Daisley et al.
2020, Raymann et al. 2017).
Transcriptome analysis is a sensitive tool to characterize sublethal
effects of potentially harmful substances on a molecular and cellular
level (Prat & Degli-Esposti 2019, Schirmer et al. 2010). Changes in
gene expression help to identify biological processes, such as stress
responses and detoxification processes, at an early stage. Exposure to
different pollutants have been shown to induce changes in gene
expression in several insect species. Mosquitos (Aedes aegypti )
exposed to anthropogenic pollutants (insecticides, PAHs) increased the
expression of genes related to detoxification, respiration and cuticular
proteins (David et al. 2010). Fireflies (Luciola leii ) showed a
similar response when exposed to benzo(a)pyrene, a widespread PAH (Zhang
et al. 2019). In different bee species, the neonicotinoids imidacloprid,
thiamethoxan, and clothianidin induce an upregulation of metabolic,
immune and stress response genes (Aufauvre et al. 2014, Bebane et al.
2019, Christen et al. 2018, Colgan et al. 2019, Gao et al. 2020, Shi et
al. 2017). The expression of genes related to detoxification was higher
in honeybees (A. mellifera ) exposed to heavy metals than in
controls (Al Naggar et al. 2020, Gizaw et al. 2020, Zhang et al. 2018).
In contrast to pesticides, the effects of other environmental
pollutants, such as particulate matter, on gene expression in bees as
well as their gut microbiome are largely unclear. To address this
knowledge gap, we exposed workers of the buff-tailed bumblebeeBombus terrestris to airborne particulate matter deriving from
traffic and investigated changes in the gut microbiome and gene
expression. Bumblebees were fed sugar water spiked with sublethal
concentrations of brake dust or diesel exhaust particles (DEPs). Adding
to this oral exposure, one group of bumblebees was exposed to DEPs via
air to enable potential uptake in the tracheal system. We expect changes
in the composition of the gut microbial community, as previous research
showed changes due to different metals in a closely relatedBombus species (Rothman et al. 2020). Moreover, we expect changes
in the expression of detoxification and metabolic genes, indicating an
increased stress level, as the toxic compounds in the particulate matter
may interfere with bumblebee physiology.
Methods
Bumblebee keeping
Four queenright colonies of B. terrestris were ordered from
Biobest (Westerlo, Belgium) in March 2021. Colonies were kept in a
climate chamber at 26° C and 70 % humidity under a constant, inverted
12:12 h light:dark cycle. Colonies were provided with sugar water (50%
Apiinvert, Südzucker AG, Mannheim, Germany) and pollen (Imkerpur,
Osnabrück, Germany) ad libitum .
Experimental procedure
At the beginning of the experiment, adult workers from the four colonies
were randomly assigned to one of six treatments. Control : fed
with sugar water only (50% Apiinvert) (n=56); Solvent control :
fed with sugar water spiked with 0.02% (v/v) of the emulsifier Tween20
(n=56); Brake dust : fed with sugar water spiked with 0.02%
(v/v) of the emulsifier Tween20 and 0.4 g/l brake dust particles (n=56);DEP :
fed with sugar water spiked with 0.02% (v/v) of the emulsifier Tween20
and 0.4 g/l diesel exhaust particles (n=56); Flight control :
fed with sugar water (50% Apiinvert) and allowed to fly once per day in
a plastic box (7 x 7 x 5 cm, EMSA, Emsdetten, Germany) for 3 minutes
(n=24); DEP flight : fed with sugar water (50% Apiinvert) and
allowed to fly once per day for 3 minutes in a plastic box (7 x 7 x 5
cm, EMSA, Emsdetten, Germany) that contained 1.5 (+-0.1) mg of diesel
exhaust particles (n=24).
The experiment was conducted in a climate chamber at 26° C and 70 %
humidity under a constant 12:12 h light:dark cycle. Bumblebees were kept
in Nicot cages (Nicotplast SAS, Maisod, France) connected to a 12 ml
syringe (B. Braun SE, Melsungen, Germany) with the tip cut off, that
contained 2 ml of the respective feeding solution (ad libitum ).
Every day the syringes were replaced with fresh ones to prevent molding
or bacterial growth in the food. The exposure lasted for seven days. At
the end of the experiment, the animals were frozen at -20° C.
Within a week after the end of the experiment, we randomly selected
twelve (three workers per colony) bumblebees per treatment for
transcriptome analysis (N=72). Additionally, for the control, solvent
control, brake dust and DEP treatment, we randomly selected 20
bumblebees (five workers per colony) for microbiome analysis (N=80),
respectively.
Generation and collection of diesel exhaust particles (DEPs)
DEPs were collected from a four-cylinder diesel engine (OM 651, Daimler
AG, Stuttgart, Germany) during a repeating cycle of transient and
stationary operating points, resembling an inner-city driving scenario
with stop-and-go intervals. The engine was operated on a test bench with
a water-cooled eddy-current brake as previously described in Zöllner
(2019). DEP samples were collected by an electrostatic precipitator
(OekoTube Inside, Mels-Plons, Switzerland). A fast response differential
mobility particulate spectrometer DMS500 (Combustion, Cambridge,
England) was applied to measure sub-micron particle size distributions
of raw exhaust samples. Depending on engine load and speed during the
inner-city cycle, solid particles showed a median diameter between 52.1
± 1.8 nm and 101.9 ± 1.7 nm. DEP composition was characterized by
thermogravimetric analysis (TGA, STA 449 F5 Jupiter, Netzsch-Gerätebau
GmbH, Selb, Germany). A fraction of 72.2 ± 1.1 % of the DEP mass was
attributed to elemental carbon, 23.2 ± 0.9 % w/w to organic fractions
and 4.6 ± 0.7 % w/w to inorganic matter. Quantification of PAHs
revealed concentrations of 444 ppm for pyrene, 220 ppm for phenanthrene,
and 107 ppm for fluoranthene.
The elemental composition of the DEP samples was analyzed by
Inductively-Coupled Plasma Optical Emission Spectrometry (ICP-OES,
Optima 7300 DV, PerkinElmer Inc., Waltham, United States of America) and
interpreted according to Zöllner (2019). It showed fractions of calcium
(1.63 % w/w), zinc (0.53 % w/w) and phosphorus (0.50 % w/w) that can
be traced back to diesel fuel and lubrication oil. Copper (1.03 % w/w),
aluminum (0.02 % w/w) and iron (0.02 % w/w) can be attributed to
abrasion of piston rings, cylinder head and engine block material,
respectively. In addition, small amounts of boron (0.13 % w/w),
magnesium (0.10 % w/w), molybdenum (0.03 % w/w), natrium (0.02 % w/w)
and sulphur (0.17 % w/w) were found.
Generation of brake dust particles
The brake dust particles provided by the Chair of Ceramic Materials
Engineering of the University of Bayreuth are derived by LowMet brake
pads (provided by TMD Friction Holdings GmbH, Leverkusen, Germany) that
were milled for three minutes in a vibrating cup mill with a tungsten
carbide grinding set (Pulverisette 9, Fritsch GmbH, Idar-Oberstein,
Germany). LowMet brake pads are common and representative for passenger
cars in Europe and consist of non-ferrous metals (25 % (w/w)), steel
wool (15 % (w/w)), petrol coke (12 % (w/w)), sulphides (10 % (w/w)),
aluminum oxide (5 % (w/w)), resin (5 % (w/w)), graphite (4 % (w/w)),
mica (4 % (w/w)), silicon carbide (3 % (w/w)), barite (2 % (w/w)),
fibers (2 % (w/w)), and rubber (1 % (w/w)) (Wiaterek 2012). The
particle size distribution of the milled, fine-grained powder was
measured with a laser diffraction particle size analyzer (PSA 1190 LD,
Anton Paar GmbH, Ostfildern-Scharnhausen, Germany). The mean particle
size found was 10.19 ± 4.37 µm (D10 = 0.68 µm (10% of all particles
being smaller in diameter than this size), D50 = 5.76 µm (median
particle size), D90 = 25.87 µm (90% of particles being smaller in
diameter than this size)).
Bumblebee gut microbiome analysis
Prior to dissection bumblebees were defrosted and rinsed in 70%
ethanol, 90% ethanol and twice in ultrapure water. We placed each
bumblebee on an autoclaved square of aluminum foil (5 x 5 cm) and opened
the abdomen with sterilized tweezers and scissor. After carefully
separating the gut from the crop and transferring it to an Eppendorf
tube, we snap-froze the gut in liquid nitrogen. All samples were stored
at -80° C until further processing.
PCR amplification and sequencing of 16S rDNA fragments
Metagenomic DNA of bumblebee gut samples was purified using the
NucleoMag DNA Bacteria kit (Macherey-Nagel, no. 744310, Düren, Germany)
after disruption of samples with 1.4 mm (diam.) ceramic beads (no.
P000912-LYSK0A, Bertin Instruments, Montigny-le-Bretonneux, France) in a
FastPrep-24 bead beating device (MPbio, Irvine, USA) following the
instructions of the manufacturer. The metagenomic DNA was diluted to a
concentration of 5 ng/µl, and 2.5 µl DNA was used to amplify 16S rDNA
fragments using primers 515F-Y (Turner et al. 1999) and 806RB (Apprill
et al. 2015) as described in the 16S Metagenomic Sequencing Library
Preparation protocol (Part # 15044223 Rev. B, www.illumina.com). Sample
libraries were barcoded using the Nextera XT index kit (v2 set A,
www.illumina.com), combined in equimolar amounts, and sequenced on
Illumina’s iSeq-100 platform using a 293 cycle single-end R1 mode.
Demultiplexing of reads was performed by the iSeq-100 local run manager
and sample-specific reads were saved in FastQ format.
Microbiome analysis
Statistical analyses of the microbial data were performed using QIIME2
(Bolyen et al. 2019) and R 4.2.1 (R Core Team 2022). Forward reads of
16S rDNA fragments (R1 reads) were analyzed using the QIIME2 microbiome
analysis package (ver. 2021.11; Bolyen et al. 2019). Unless indicated
otherwise, all analysis tools were used as plugins of the QIIME2
package. The respective parameters used along the analysis steps are
readily accessible by provenance information in the QIIME2 data files
(available as supplemental data). In brief, the following analysis steps
were performed: Demultiplexed reads were trimmed for 16S primer
sequences (plugin cutadapt; Martin 2011), denoised, dereplicated and
chimera-checked (plugin DADA2; Callahan et al. 2016) resulting in
amplified sequence variants (ASVs). Rare ASVs were filtered using the
median frequency of ASVs over all samples. Taxonomic classification of
ASVs was performed (plugin feature-classifier; Bokulich et al. 2018)
using the pre-fitted sklearn-based taxonomy classifiers based on the
SILVA reference database (ver. 138.1; (Quast et al. 2013; Yilmaz et al.
2014). ASVs that could not be taxonomically assigned at any taxonomic
level (‘unassigned’) as well as samples with less than 3,900 reads in
total were removed prior to subsequent analysis steps. Alpha diversity
metrics, such as Shannon diversity index, Faith’s phylogenetic
diversity, Pielou’s eveness, and observed ASVs, were obtained using the
QIIME2’s ‘core-metrics-phylogenetic’ workflow (plugin diversity),
rarefied to 3,900 reads per sample. To find significant differences in
α-diversity we fitted generalized linear models (GLMs) with treatment as
factor. We checked model assumptions using model diagnostic test plots,
i.e., qqplot and residual vs. predicted plot from the packageDHARMa (Hartig 2022). Depending on model assumptions, we then
used Kruskal-Wallis tests or produced F-statistics with the functionAnova () from the package car (Fox & Weisberg 2019) to
calculate p-values for differences between treatments. For significant
treatment effects, we ran pairwise comparisons. In the case of a
significant Kruskal-Wallis test, pairwise comparisons were done using
Dunn’s test for multiple comparisons with Benjamini-Hochberg correction
(package dunn.test (Dinno 2007)). In the case of a significant
ANOVA, pairwise comparisons were made using Tukey HSD post-hoc test with
Benjamini-Hochberg correction from the package multcomp (Hothorn
et al. 2008). Differential abundance of the rarefied data we analyzed
using the package DESeq2 with a negative binomial distribution, a
significance level cutoff of FDR < 0.01, replacement of
outliers turned off, and cooksCutoff turned off (Love et al. 2014).
Compositional differential abundance analysis was performed using Aldex2
(plugin aldex2; Fernandes et al. 2013). Beta diversity of the sparse,
compositional microbiome data were calculated using QIIME2’s plugin
DEICODE which performs a robust Aitchison PCA (Martino et al. 2019).
Significance was tested in a PERMANOVA with 999 permutations followed by
pairwise PERMANOVA with Benjamini-Hochberg (BH) correction for multiple
testing (Anderson 2008). We used the packages qiime2R (Bisanz
2018) and mia (Ernst et al. 2022a) to import and process the
microbiome data in R. Data were arranged using the package tidyr(Wickham & Girlich 2022) and were plotted using the packagesggplot2 (Wickham 2016), ggpubr (Kassambra 2020), andmiaViz (Ernst et al. 2022b).
Transcriptome analysis of whole bumblebee abdomens
Bumblebees were defrosted and rinsed in 70% ethanol, 90% ethanol and
twice in ultrapure water prior to dissection. The abdomen was cut off
with sterile scissors, placed in an Eppendorf tube and snap-frozen in
liquid nitrogen. All samples were stored at -80° C until further
processing.
RNA sequencing
Total RNA was prepared from abdomen samples using the RNeasy Lipid
Tissue kit (Qiagen, no. 74804, Hilden, Germany). RNA-Seq libraries were
constructed from 100 ng RNA using the NEBNext Ultra II Directional
Library Prep Kit for Illumina (New England Biolabs, no. E7760, Ipswich,
USA) in combination with the NEBNext Poly(A) mRNA Magnetic Isolation
Module (New England Biolabs, no. E7490, Ipswich, USA). The samples were
combined at equimolar amounts and sent out for sequencing on an Illumina
device in 150 bp paired-end mode (Genewiz, Leipzig, Germany). A total of
1.470 million reads, corresponding to an average of 19.5 million reads
per sample, were obtained.
Differential expression analysis
RNA-Seq reads were further analyzed using the OmicsBox bioinformatics
platform (v. 2.0.36, www.biobam.com).
Unless indicated otherwise, all tools used for differential expression
analyses are accessible within the OmicsBox platform. RNA-Seq reads were
preprocessed by Trimmomatic (Bolger et al. 2014) to remove sequencing
adapters, low-quality sequences, and short reads from the dataset. The
quality-trimmed reads were mapped to the B. terrestris genome
assembly (Bter_1.0, GCA_000214255.1, downloaded from
metazoa.ensembl.org) using STAR (Dobin et al. 2013). A gene-specific
count table was created from the mapping files using HTseq (Anders et
al. 2015) and differentially expressed genes were identified by edgeR
(Robinson et al 2010), respectively. Functional annotation of theB. terrestris genome was based on annotation release v. 102
(available in gff3 format from metazoa.ensembl.org). Since 4,975 of the
12,008 genes did not contain any functional annotation, the functional
annotation workflow of the OmicsBox platform was used to update the
published annotation with additional information. In brief, the coding
sequences of unannotated genes were used to extract functional
annotations from refseq_protein database (www.ncbi.nlm.nih.gov) and
InterProScan (www.ebi.ac.uk). These we then
fed into the GO mapping and annotation tools of the pipeline and finally
merged to the existing functional annotations. Gene Set Enrichment
Analyses (GSEA; Subramanian et al. 2005) were performed using ranked
list of genes (rank = sign(logFC) * -log10(p-value); FC: fold-change)
and gene sets defined by Gene Ontology’s functional annotations. For the
functional network analysis of enriched GO terms we used ClueGo (v.
2.5.9; Bindea et al. 2009) and CluePedia (v. 1.5.9; Bindea et al. 2013)
plugins in Cytoscape (v. 3.9.1; Shannon et al. 2003). We used the
packages ggplot2 (Wickham 2016), ggpubr (Kassambra 2020)
and pheatmap (Kolde 2019) to plot transcriptome data in R 4.2.1
(R Core Team 2022).
Results
Effect of pollutants on the bumblebee gut microbiome
Amplicon sequencing of the bacterial 16S rDNA fragments yielded a total
of 2,425,928 raw reads. After quality filtering and removal of
unassigned sequences, we also removed samples with a sampling depth
below 3900 reads (n=7), all from DEP treatment, to ensure adequate
sampling depth (13 DEP replicate samples remained in the analysis). In
the remaining samples we obtained 1,856,025 16S rDNA gene sequences with
a mean of 25,425 reads per sample (n=73), corresponding to 468 amplicon
sequence variants (ASVs). Sample-based rarefaction curves suggest a
sufficient sequencing depth for a representative coverage of the
microbiome as most of the samples reach a plateau (Figure A1 ).
Taxa abundance
On the genus level, the most common bacterial taxa (> 1 %
in at least one treatment) were: Gilliamella, Snodgrassella,
Lactobacillus, Asaia, Bombiscardovia, Methylorubrum andBombilactobacillus. The relative abundance of the most common
genera for each sample shows a different microbial composition in the
DEP treatment compared to the other treatment groups (Figure
1 ).