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 ).