Genetics, age, and diet influence gut bacterial communities and performance of black soldier fly larvae (Hermetia illucens)

Abstract

Background

The gut microbiota of black soldier fly larvae (BSFL, Hermetia illucens) play a crucial role in recycling various organic waste streams. This capability is linked to the presence of a potential common core microbiota in BSFL. However, subjective thresholds for defining core taxa and the difficulty of separating genetic and environmental influences have prevented a clear consensus in the literature. We analysed the gut bacterial communities of two genetically distinct BSF lines (wild type (WT) and lab-adapted line (LD)) raised on ten different diets based on common agricultural by-products and food waste in Southeast Asia.

Results

High-throughput 16S rRNA gene sequencing revealed that gut bacterial communities were significantly influenced by genetics (p = 0.001), diet (plant/meat-dominated; p = 0.001), larval age (p = 0.001), and the interactions between all three (p = 0.002). This led us to investigate both common core taxa and lineage-specific core taxa. At a strict > 97% prevalence threshold, four core taxa were identified: Providencia_A_732258, an unclassified genus within the family Enterococcaceae, Morganella, and Enterococcus_H_360604. A relaxed threshold (> 80% prevalence) extended the core to include other potential common core taxa such as Klebsiella, Proteus, and Scrofimicrobium. Our data suggest that Proteus, Scrofimicrobium, Corynebacterium, Vagococcus_B, Lysinibacillus_304693 (all LD), and Paenibacillus_J_366884 (WT) are lineage-specific rather than members of a common core (> 90% prevalence in either LD or WT, with prevalence significantly different between lines (p ≤ 0.05)). Positive correlations were observed between several core genera and larval performance in LD, typical of a highly optimized lab-adapted line. Interestingly, only members of the genus Providencia appeared to play a crucial role in most aspects of larval performance in both genetic lineages.

Conclusion

Our study demonstrates that the gut microbiota of BSFL is influenced by genetic factors, diet composition, larval age, and their interactions. We identified a distinct lineage-specific core microbiota, emphasizing genetic background’s role. Future studies should apply a standardized high prevalence threshold of at least > 90% unless there is a valid reason for relaxation or sample exclusion. The consistent association of Providencia spp. with larval performance across both genetic lines highlights their crucial role in the BSFL gut ecosystem.

Introduction

Hermetia illucens, commonly known as the black soldier fly (BSF), has garnered attention for its exceptionally high nutritional content, sustainable and safe mode of rearing, and valuable biomaterials, which make it a promising candidate for versatile applications in the food and feed industry. It has been shown that feeding substrates, particularly their macronutrient composition (i.e., protein, fat, and carbohydrate content) influence the development and nutritional value of BSF larvae (BSFL; e.g., protein and lipid content) [5, 65]. In other insects, it has been well established that the gut microbiota, a dynamic and diverse community of microorganisms residing in the gastrointestinal tract, profoundly influences diverse aspects of host physiology such as nutritional uptake, immunity, development, and ecological interactions [45, 60]. In BSFL, the gut microbiota has been shown to have promotive effects on growth and development [35]. Therefore, in recent years, several studies have explored the influence of biotic (e.g., initial substrate bioburden) and abiotic (e.g., density, temperature, water content, pH) factors in shaping the gut microbial community of BSFL as well [15, 44]. The presence of a potential “core” gut microbiota together with the microbial community acquired from the diet are thought to contribute to the substrate versatility of BSFL. However, while some studies suggest that there is a core microbiota essential to the larva’s survival [31], others suggest that there is no core microbiota shared amongst all BSFL [15] or that the taxonomic composition of the core microbiota differs with experimental variables [25]. Adding to the difficulty, different studies used different criteria to define a core taxon, and much debate remains around the question as to whether a genus can be considered part of a “core” if it is not 100% prevalent in all samples [15, 21].

A possible missing link to this controversy may be underlying genetic factors of the BSF populations used across different studies. So far only a handful of studies have investigated genotype (G), and genotype-by-environment (G × E) interactions on BSFL performance (e.g., biomass gain) and body composition (e.g., protein and lipid content) [51, 67]. However, genetic factors appear to contribute significantly to the variation observed in the gut microbiota of BSFL. In 2020, a study investigated the genetic variation and diversity of gut microbes among geographically different BSF populations [27]. Greenwood et al. [21] revealed that both genetics and diet influenced the larval gut microbiota of two geographically different BSF populations with varying duration in captivity [21]. This finding aligns with similar studies conducted on other wild vs. lab-adapted animal species such as Drosophila melanogaster and Caenorhabditis elegans [50], whereby the gut microbiota underlies both host genetic pressure as well as environmental pressures (e.g., diet) [11, 32]. However, the limited diversity of diets tested in the context of BSFs, predominantly comprising generic meat or plant-dominated options but not heterogeneous in the Asian context like Singapore’s food wastes rich in oils, spices, and carbohydrates, as well as unspecified sampling timepoints or larval age hamper a more comprehensive understanding. Additionally, assumptions about genetic diversity in some studies further impede the exploration of the dynamic interaction between host genetics and the gut microbiota of BSFL. Understanding the factors that shape the gut microbiota, particularly the interplay between taxon variations and host genetics, is essential for harnessing the potential of BSFs in applications such as waste management and sustainable resource utilization.

In a recent study, we used BSFL from two genetically variant populations fed on more than ten different diets under standardised conditions and showed that genetic variation had a significant effect on the plasticity in larval development and bioconversion ability as well as on larval nutritional profiles [65]. The cause of the observed variation is uncertain, leading us to the investigation of the gut microbial aspect. Therefore, the primary objective of this study was to (a) examine the gut microbiota composition of the BSFL of the two different lineages on different diets over time, and (b) obtain insight into the contribution of different gut bacterial taxa to larval performance. The secondary aim was to explore the presence of core bacterial taxa shared across or within BSFL populations that may facilitate such adaptive mechanisms previously reported. We hypothesize that if a common core gut microbiota exists in BSFL, this core will be consistently discovered across genetic lines, diets, and larval age.

Methods

Genetic lineages of BSF

The experiment to assess the influence of lineage on gut microbiota employed two genetically distinct lines: a lab-adapted line (LD), established in 2018 by crossing Southeast Asian and North American cultures and selectively bred on National University of Singapore dining hall surplus food (NUS) [17], and a wild-caught line (WT) from City Sprouts, Singapore (1° 17′ 04.7″ N; 103° 49′ 12.4″ E) in September 2020, which was maintained on commercial chicken feed (CF; PK Agro-industrial Products M Sdn Bhd; Johor, Malaysia) for four months and on NUS for two generations prior to the commencement of the experiment (for further information, refer to Zhang et al., [65].Genetic variability of both populations was confirmed by restriction-site associated DNA sequencing (RAD-seq; Figure S1). Both populations have been maintained in natural light under ambient levels of relative humidity (64 – 96%) and temperature (~ 27 °C). The BSF nurseries for hatching have been maintained on CF.

Diet preparation and cage set-up

Diet preparation and cage set up were as previously described [65]. Briefly, four replicates were set up for each diet, and a total of 220 larvae were used for each replicate. A total of ten different diets, namely CF (chicken feed, control diet), OKA (okara), PKM (palm kernel meal), SBM (soybean meal), RIB (rice bran meal), P (coconut (25%) + okara (75%)), MLP (fishcake (5%) + duck offals (20%) + coconut (5%) + okara (70%)), M (fishcake), MHP (fishcake (15%) + duck offals (60%) + coconut (5%) + okara (20%)), and NUS (university dining hall surplus food) were used (Table 1). The selected substrates were predominantly plant-based, with just a few meat-based ones, due to their abundant availability from an agricultural company that produces cooking oil and soy products. Additionally, some of these substrates are nutritionally poor for BSFL development. By using a range of carbohydrate to protein ratios, we observed developmental effects on BSFL (pure duck offal was excluded due to the high mortality rate of the larvae on the diet). As discussed in Zhang et al. [65], genetic variation and diet significantly impacted larval development plasticity, bioconversion ability, and larval nutritional profiles. A single feeding regime was selected in order to ensure a practical approach for potential future industrial applications.

Table 1 Summary of the composition of diets tested to observe the influence of various diets on the gut microbiota of two genetic populations of BSFs

Sample collection and gut isolation

To set the baseline, three LD and three WT larvae were collected at the experimental start date (day 0) and snap frozen on dry ice. Three larvae were collected from each of three replicates from each test diet on day 5, day 10, day 15 and day 20 (experimental endpoint: larvae entering prepupal phase) of the experiment and snap frozen as well to analyse the changes in the gut microbiota over time. All samples were stored at -80 °C until further processing.

The dissection procedure followed the technique outlined by Klammsteiner et al. [31]. Three whole larval guts per replicate were excised and pooled together to ensure sufficient material for subsequent DNA extraction. Before and after each excision, 70% ethanol was used to sterilize all surfaces and instruments. Isolated gut samples were stored at -80 °C until further processing.

Microbial community profiling

DNA was extracted from isolated gut samples using the QIAamp^® PowerFecal^® Pro DNA Kit (QIAGEN) following the manufacturer’s handbook, employing the vortex adapter method for the homogenization of tissue samples. Before the final centrifugation step, 50 µl of C6 solution were added onto the MB Spin Column filter membrane and the column left to incubate at room temperature for 5 min to ensure sufficient time was provided for the DNA to completely dislodge.

Polymerase chain reaction (PCR) amplification was done targeting the V3-V4 region using primer pair 338 F (5′-ACT CCT ACG GGA GGC AGC AG-3′) and 806R (5′-GGA CTA CHV GGG TWT CTA AT-3′) [66]. Unique 12-bp barcodes were added to the forward primer for multiplexing (Table S1). For each PCR reaction, 10 µl of Q5^® High-Fidelity DNA Polymerase (New England Biolabs, Singapore), 0.5 µl of forward and reverse primer each (10 µM), 50 to 100 ng DNA template, and water to reach a final volume of 25 µl were used. Each reaction was set up in triplicates to minimise amplification bias and no template controls were included to ensure absence of contamination. The PCR cycling conditions were as follows: denaturation for 3 min at 95 °C, 33 cycles of 30 s at 95 °C, 30 s at 55 °C, 45 s at 72 °C, followed by final extension at 72 °C for 10 min (Bio-Rad T100 Thermal Cycler). Verification of successful amplification was conducted through gel electrophoresis. PCR amplicons (~ 500 bp) were then purified using the QIAquick^® PCR Purification Kit (Qiagen, cat# 28106), quantified using the QuantiFluor^® ONE dsDNA System (Promega) and Quantus™ Fluorometer (Promega), and pooled together in equimolar concentrations for library preparation (four separate libraries were generated, each representing amplicons obtained from 32 DNA samples). The pools were then sent to Macrogen (Singapore) for TruSeq Nano DNA Library Prep followed by Illumina MiSeq 300PE amplicon sequencing. Out of the 126 samples, one sample (WTD5CFR1) failed sequencing QC and could not be included in downstream analyses.

Bioinformatics analysis

Demultiplexed paired-end reads were imported into QIIME2 (version 2021.11) (Bolyen et al., 2019). Denoising was done using the DADA2 pipeline via the q2-dada2 plugin (Callahan et al., 2016). Sequences were carefully trimmed at the 5’-end by exactly the length of the barcode plus linker plus forward-primer and at the 3’-end by the length of the reverse-primer and truncated towards the end of the read according to quality scores prior to merging and chimera checking. A classifier pre-trained on “Greengenes2 2022.10” was used for taxonomic assignment [39].

Statistical analysis

Artifacts generated from QIIME2 were imported into RStudio (version 4.1.2) for statistical analysis using unrarefied data as rarefying artificially limits the sample size thus limiting the power of statistical hypothesis testing [24, 40]. Good’s coverage (rarefaction curve analysis) was performed using OTU output at the genus level using the formula 1 - (F1/N) where F1 is the number of singleton OTUs and N is the sum of counts for all OTUs.

Microbial compositions were analysed using packages phyloseq (version 1.38.0) and qiime2R (version 0.99.6, github). Shannon α-diversity was calculated for a quantitative view of microbial community richness followed by ANOVA to test for statistical significance of variables such as age, genetic lineage and diet. Statistical differences in microbial richness between diets were validated using Wilcoxon rank-sum test corrected by Benjamini-Hochberg (FDR) method (reported as q-value).

For evaluation of β-diversity, the Bray-Curtis dissimilarity metric was used. Student’s t-test was performed to test for a significant difference between all within-lineage and between-lineage dissimilarities. Principal Coordinates Analysis to determine correlation of variables to gut bacterial composition was done based on the Bray-Curtis dissimilarity matrix for individual diets sub-categorised by genetic line. Taxonomic abundance bar plots were generated using Microsoft excel (version 2406). Redundancy analysis (RDA) was performed, and significance of models was validated using microviz (version 0.11.0) [4]. Spearman’s Rank correlation coefficients were derived using Hmisc (version 5.1.1) to explore the correlations of (a) dietary macronutrients, and (b) larval performance indicators with gut bacterial taxa (raw data obtained from Zhang et al., [65] corrected for tied ranks using “rcorr” function [3, 56]. ANOVA was applied to calculate the contribution of each dietary macronutrient to larval gut microbiota composition. SIMPER (Similarity Percentage) analysis, a multivariate statistical approach, was used to detect the dissimilarities in bacterial taxa observed between genetic populations, and between D5 and D10 larvae followed ANOVA for statistical significance. Core microbiota analysis was performed at the genus level, applying several thresholds (prevalence = 50–99.9% of samples (in 10% increments)) and a low abundance limit (detection = > 0) to prevent the exclusion of rare taxa.

Results

Sequencing coverage and bacterial diversity in the gut of BSFL

The primary objective of this study was to examine the composition of the gut microbiota in relation to genetics, diet, and larval age, and to decipher the involvement of specific taxa in the use of dietary macronutrients and larval performance. Approximately 30,000 raw reads per sample were collected from Illumina MiSeq. After thorough quality assessments and filtering, the average sequence yield per sample, suitable for further analysis, was 22,777 ± 10,053 reads (standard deviation (SD)). Good’s coverage was determined to be 99.99% ± 0% as calculated from the QIIME sample-frequency output at the genus level, and 99.98% ± 0% at the OTU level. This metric provides insight into the completeness of sampling in capturing the bacterial diversity present in the dataset. Shannon α-diversity was evaluated based on diet, age of larvae, and genetic line (Fig. 1).

Fig. 1

Bacterial community diversity of BSFL as assessed using the Shannon α-diversity metric before (D0) and after the introduction of test diets (D5, top; D10, bottom) for Line D (LD, black bars) and the wild-type line (WT, white bars). The diets are arranged from plant-dominated diet types (CF to P) on the left to meat-dominated diet types (M to NUS) on the right. Diet abbreviations are provided in Table 1

While overall, a slightly higher diversity was observed for LD as compared to WT (p = 0.018; Student’s t-test), the differences decreased with larval age (D5: p = 0.08) and were insignificant at D10 (D10: p = 0.21). Similarly, no significant differences in bacterial richness were observed between experimental days for each genetic line (p > 0.05). However, differences in bacterial richness were observed between genetic lines fed on the same diet as well as within lines between certain diets (p < 0.05). Higher bacterial richness was observed on plant-dominated (CF, OKA, PKM, RIB, SBM, P and MLP) as compared to meat-dominated diets (M, MHP and NUS; p < 0.001). This was observed independently for both genetic lines. To further investigate the differences in α-diversity between individual diets and genetic lines, pairwise Wilcoxon rank-sum tests were performed, and results corrected for FDR. In both WT and LD, significant differences were observed between certain plant-dominated diets and meat-dominated diets; for example, MLP, OKA, P, and PKM (both lines), as well as RIB (only LD) were significantly different to M (all q < 0.05), and MLP, and P (both lines), as well as OKA (only LD) were significantly different to MHP (all q < 0.05). In WT however, no significant differences were observed between CF and other diets (all q > 0.05). Furthermore, differences were observed even between individual plant-dominated diets in WT whereby MLP induced a higher bacterial richness than RIB and SBM (both q < 0.05). A similar trend was observed in LD as well, however, with interesting additional differences observed between diets. Within plant-dominated diets, there was a difference in richness induced by CF as opposed to MLP and P, and P as opposed to RIB and SBM (all q < 0.05). Notably, the bacterial richness induced by P (heterogeneous pure plant diet) was significantly higher in LD as compared to WT (q < 0.05). Amongst all the diets tested, NUS (the most heterogenous diet) induced the least bacterial richness consistently across both genetic lineages.

Influence of genetics and larval age on gut bacterial community composition

We further aimed to investigate how genetic and environmental factors influenced bacterial community composition in the larval gut. Overall, Firmicutes D and Proteobacteria were the most abundant phyla with relative abundances of 42.3% ± 18.8% and 34.4% ± 19.5% in LD and 44.4% ± 18.2% and 29.3% ± 18.3% in WT, respectively, followed by Firmicutes A (LD: 7.5% ± 9.9%, WT: 11.3% ± 14.4%), Actinobacteroidota (LD: 7.7% ± 10.4%, WT: 6.6% ±11.8%), and Bacteroidota being the least abundant phylum (LD: 7.3% ± 10.1%, WT: 6.6% ± 7.8%) (Fig. 2). Firmicutes D and Proteobacteria were predominantly present in D0 samples, together accounting for 77.6% ± 4.7% and 80.5% ± 11.4% in LD and WT bacterial communities, respectively. It was obvious that, despite having been fed the same diet, the bacterial communities of each of the genetic lineages were enriched in different phyla. For example, when fed on meat-dominated diets, Proteobacteria were more abundant in LD as compared to WT (LD: 55.9% ± 15.1%, WT: 44% ± 13%, p = 0.016). Inversely, Firmicutes D were significantly more abundant in WT larvae when fed on meat-dominated diets (LD: 36.3% ± 13.9%, WT: 45.3% ± 11.2%, p = 0.041). Campylobacterota appeared to be a rare taxon in general (maximum relative abundance of only 6.7% across all samples). However, interestingly, this phylum was observed in most samples collected from LD regardless of the diet being offered (except for CF and NUS), while absent from WT.

Fig. 2

Overall relative abundance of taxa at the phylum level across both genetic lineages (LD: black bars; WT: white bars), aggregated irrespective of age or diet

Multivariate analysis was done using bacterial profiling data at genus-level resolution to determine if genetic lineage, larval age, or diet alone explained gut bacterial differences, or if they exerted combined effects in shaping the BSFL bacterial communities. Interestingly, the combined effect of genetic line and larval age showed the strongest effect on shaping the gut bacterial community (p = 0.001, residual degrees of freedom = 80). Based on the Bray-Curtis dissimilarity metric, average within-line dissimilarity was higher for WT than LD (0.78 ± 0.15 vs. 0.73 ± 0.14; p < 0.001 (Student’s t-test)). PCoA was used to further visualise genetic differences and temporal changes inbacterial community structure across different diets (p = 0.001, Table S2; Fig. 3).

Fig. 3

Principal Coordinates Analysis of the BSFL bacterial gut microbiota in two different BSF populations feeding on ten different diets and sampled at three different time points. The Bray-Curtis dissimilarity metric was used. Refer to Table 1 for the details of the diets. The filled circles depict LD and open circles depict WT

D0 samples from the two genetic lineages showed distinct clustering. D5 and D10 samples from the two genetic lineages fed on P, MLP, and M showed similar bacterial community compositions respective to their genetic lineages, suggesting a divergence of bacterial communities in the two lines resulting from differential adaptation to these substrates. However, D5 and D10 samples from larvae fed on MHP clustered together irrespective of genetic lineage, except for D5 WT larvae, which clustered away from the other samples. For two of the diets, PKM and RIB, D5 samples displayed proximity to D0 samples, rather than D10 samples as observed for the remainder of diets. This suggested a slower adaptation of bacterial communities to these fibre-rich plant substrates. Overall, the diets showed to be an important driving factor resulting in some larvae displaying stronger genetic effects and others displaying stronger age effects in the shaping of their gut microbiota.

SIMPER (Similarity Percentage) analysis, a multivariate statistical approach, revealed specific bacterial genera underlying the observed differences between the genetic lineages (Table 2; Table S3).

Table 2 SIMPER (similarity percentage) analysis was used to identify genus-level taxa that showed significantly different relative abundances between the two genetic lineages, larval age (D5 and D10), and between both genetic lineages across the tested diets. The table displays the average relative abundances ± SD (%) of significant bacterial genera. Only those diets on which certain bacterial genera were differentially enriched between the genetic lines are shown. The significance codes (Signif. codes) are provided for the corresponding p-values. The complete list of diets is presented in Table S3

Levilactobacillus and Morganella were present in higher abundances in LD compared to WT (all p < 0.05), while Paenibacillus_J and Lacrimispora were more abundant in WT than LD (all p < 0.05). Further investigation of the most influential genera on genetic lineage based on diets showed Morganella to be significantly higher in MHP fed LD larvae (LD: 21% ± 15.9%, WT: 5.4% ± 5.1%; p < 0.05; Table 2). Interestingly, within the genus Paenibacillus, Paenibacillus_C and Paenibacillus_J appeared to be differentially abundant in P and MLP fed larvae. In contrast to Paenibacillus_J, Paenibacillus_C was significantly more abundant in LD when fed on P or MLP (both p < 0.001). However, no significant differences were observed between genetic lineages when larvae were fed on CF, RIB, PKM, NUS, and M.

Following the same statistical approach, we next identified genera that were significantly associated with the age of the larvae (i.e., D5 and D10), regardless of lineage and diet (Table 2). Lactiplantibacillus, Klebsiella, Levilactobacillus, and Lysinibacillus (all p < 0.05) exhibited significantly higher abundance in D5 samples as compared to D10 samples. Conversely, Dysgonomonas, Ignatzschineria, Scrofimicrobium, Sphingobacterium, and Morganella (all p < 0.05) were significantly more abundant in D10 samples.

Influence of diet on gut bacterial community composition and correlation of bacterial taxa with dietary macronutrients

Multivariate analysis also revealed a significant effect of diet on the composition of bacterial communities in the gut of individual BSFL populations (p = 0.001, Table S2). Members of Firmicutes A, for example, were most highly abundant in MHP, MLP, P and PKM diets (> 15% in relative abundance whilst their abundance in other diets was below 5%; Fig. 4A).

Fig. 4

Relative abundances of the most abundant bacterial taxa found in the gut of BSFL fed on the ten different diets. A) Relative abundances of most dominant phyla across all tested diets. B) Relative abundances of top 15 most abundant bacterial genera across all tested diets. “Other < 1%” represents all remaining taxa, including those with < 1% relative abundance in the samples

In addition, Proteobacteria were higher in abundance in meat-dominated diets, and lower in plant-dominated diets (meat-dominated: 50% ± 15.1%, plant-dominated: 25% ± 15.5%, F_{(1, 117)} = 66.2, p < 0.001), whereas Actinobacteroidota and Bacteroidota were observed in plant-dominated diets, but were almost absent from meat-dominated diets (plant-dominated: 10.1% ± 12.3% and 10% ± 9.6%, respectively; meat-dominated: 0.4% ± 0.3% and 0.9% ± 1.9%, respectively, F_{(1, 117)} = 24.5 p < 0.001, F_{(1, 117)} = 31.5 p < 0.001).

Not only did we detect an effect of diet category (meat-dominated vs. plant-dominated) on shaping the gut bacterial communities (Fig. 4B) but dietary macronutrients also appeared to be significantly correlated with certain bacterial genera (Table S4, Figure S2) resulting in different gut bacterial communities between diets (RDA analysis, Fig. 5). For example, both types of plant-dominated diets, MLP and P, appeared to harbour members of Dysgonomonas (MLP: 11.9 ± 9.0% and P: 11.8% ± 8.6%), Lacrimispora (MLP: 11.9% ± 7.6%, P: 8.4% ± 6.1%), Paenibacillus_J (MLP: 7.5% ± 7.7%, P: 14.3% ± 7.2%) and Paenibacillus_C (MLP: 3.1% ± 3.2%, P: 6% ± 7.1%), whereas these genera were present at considerably lower relative abundances in samples from all other diets (< 2.5%). This was in congruence with the RDA results, whereby an association of these four distinct genera was observed with fat and protein content of the diet (Fig. 5, Figure S2).

Fig. 5

Redundancy analysis (RDA) biplot showing the correlation between bacterial taxa and dietary macronutrients. The taxonomic data were transformed using the Hellinger transformation with a scaling factor of 2. Dietary macronutrients as well as the top 10 bacterial taxa are illustrated as arrows, denoting both their direction and the degree of association with the samples. Sample positions relative to these arrows indicate their similarity or dissimilarity based on the variables. Arrows pointing towards sample clusters indicate strong influence of a variable on those samples. The angles between arrows reflect correlations among variables; smaller angles denote stronger correlations, while larger angles signify weaker ones. Taxa situated in the direction of the macronutrients may be significantly impacted by them

Genera such as Scrofimicrobium and Sphingobacterium seemed to be more prominent in larvae fed on plant- as compared to those fed on meat-dominated diets (meat-dominated: 0.3% ± 0.3% and 0.1% ± 0.4%, plant-dominated: 8.9% ± 12.2% and 2% ± 4.8% respectively, F_{(1, 117)} = 17, p < 0.001 and F_{(1, 117)} = 5.8, p < 0.05; Table 1; Fig. 4B). Sphingobacterium also displayed a strong positive association with carbohydrate content in the homogeneous plant-dominated diets. Inversely, Providencia_A was generally higher in abundance in meat-dominated diets (meat-dominated: 32.3% ± 20.6%, plant-dominated: 8% ± 8%, F_{(1, 117)} = 85.3, p < 0.001). Similarly, Morganella and Proteus showed significantly higher relative abundances on meat-dominated diets compared to plant-dominated diets (p < 0.05 and p < 0.001, respectively). Members of the genus Lactiplantibacillus and Levilactobacillus were highly abundant in CF and in other plant-dominated diets but interestingly in M (pure meat) as well. Klebsiella, one of the most prominent genera across all samples, did not show any significant association with plant or meat-dominated diets (p > 0.05). Similarly, members of the family Enterococcaceae (classified and unclassified) were present in larvae on all diets.

Bacterial community structure and larval performance

Spearman’s Rank correlation test was used to assess for any significant relationships between larval performance and gut bacterial taxa across the different diets as observed by Zhang et al. [65] (Table 3). For this, bacterial community composition at D10 was used as this closely corresponded to the larvae’s harvesting age used for performance analyses [65]. Due to the identified genetic influence on bacterial community composition, the two different genetic lines were analysed separately. Out of several parameters analysed, for the current study, biomass conversion ratio (BCR) of the larvae, prepupal weight (PW), protein and fat conversion rate, prepupal protein content and prepupal fat content were selected for the correlation analysis [65] due to their relevance in industrial applications.

Table 3 Spearman’s Rank coefficient was calculated to explore the correlation between larval performance (biomass conversion ratio (BCR), prepupal weight (PW), protein and fat conversion rates, prepupal fat, and protein content) and gut bacterial taxa at the genus level on experimental day 10. A positive value indicates a positive correlation, while a negative value indicates a negative correlation. Significance of correlating parameters was determined using Student’s t-test followed by calculation of p-values (levels of significance are indicated by asterisk(s)). Only significant taxa (p < 0.05) are shown, with their corresponding Spearman’s Rank coefficients (R-values). Taxa that were present at less than 1% relative abundance in all samples were summarized and named “<1%”

For both lines, several taxa were identified that correlated either positively or negatively with the various parameters of larval performance. In LD, the genus that showed the most and some of the strongest positive correlations with larval performance was Providencia. Providencia was positively correlated with BCR, PW, protein conversion, fat conversion, as well as prepupal fat content (Table 3). Members of the genus Proteus also showed significant positive correlations with BCR, PW, and fat conversion rate. The genus Bacillus was positively correlated with PW and BCR, while Klebsiella demonstrated a positive correlation with protein conversion, and Morganella was positively correlated with prepupal protein content. In contrast, Scrofimicrobium, Dysgonomonas, and Sphingobacterium exhibited negative correlations with most larval performance metrics, except for protein conversion rate whereas Corynebacterium exhibited negative correlations with protein conversion rate. Paenibacillus and Lacrimispora correlated negatively with fat conversion rate.

In WT, fewer significant correlations were observed compared to LD. The only two genera that showed positive correlations with larval performance were Providencia as well as Peptostreptococcus (Table 3). Providencia displayed positive correlations with all larval performance metrics, except for protein conversion rate and prepupal protein content. While members of the genus Peptostreptococcus showed positive correlations with BCR, PW, and protein conversion rate. On the contrary, Scrofimicrobium exhibited negative correlations with PW and fat conversion rate, while Dysgonomonas was negatively correlated with only fat conversion rate and Sphingobacterium was negatively correlated with only PW. Morganella displayed negative correlations with several performance parameters namely protein conversion rate, prepupal fat and protein content while Lactiplantibacillus was negatively correlated with BCR.

Prevalence threshold determines the range of core bacterial taxa

Thus far, several studies have identified a core microbiota in BSFL. However, the methods used differed between studies. Authors opted to include either (1) prevalence data only (detection of taxa in a defined percentage threshold of samples), (2) abundance data only (relative abundance of taxa is set as the threshold across the samples thus uncommon and rare taxa that fall below this threshold will be omitted), or (3) a combination of both [42, 49]. To investigate the effect of methodology on the ability to detect an overall core across all samples in this data set as well as to investigate the possibility of lineage-specific core taxa notwithstanding age, and diet, the percent prevalence was calculated for each genus-level taxon across all samples (Fig. 6A), and within each fly lineage (Fig. 6B).

Fig. 6

Analysis of shared genera across all samples (A), and across the two different genetic lines used in this study (B) notwithstanding experimental day (larval age) or diet. Dashed lines indicate the prevalence ranges. Only taxa that showed a prevalence of > 50% across all samples or in at least one line are shown. Taxa that showed a significantly different prevalence between the two lineages, but with one lineage fulfilling the prevalence criterion of > 90% are marked as significant with asterisks as follows: *: p ≤ 0.05, **: p ≤ 0.01, ***: p ≤ 0.001

The threshold plays a crucial role in data interpretation, as stricter thresholds naturally detect a lower number of shared “core” genera. While as many as 30 genera were shared across > 50% of all samples, only two genera, namely Providencia_A_732258 and an unclassified genus within the family Enterococcaceae, were shared across 100% of all samples regardless of genetic lineage, diet, or age of the larvae and represented the core in this study. Moreover, Morganella and Enterococcus_H_360604, occurred in 100% of LD samples and > 95% of WT samples. These may be argued to belong to the core microbiota as well. Relaxation of the prevalence threshold to > 80% extended the core by a further nine genus-level taxa (Fig. 6A). Furthermore, it was evident that certain bacterial genera showed significantly different prevalences depending on the genetic lineage of their host, suggesting the existence of lineage-specific core microbiota. At a prevalence threshold of > 90%, Proteus, Scrofimicrobium, Corynebacterium, Vagococcus_B, and Lysinibacillus_304693 were significantly more prevalent in LD, while Paenibacillus_J_366884 was significantly more prevalent in WT (p ≤ 0.05; Student’s t-test; Fig. 6B).

Discussion

BSF gained popularity in recent years for various potential applications [53]. Despite several studies that have investigated the effects of genetics on performance of BSFL [27, 67], only one study thus far analysed gut bacterial communities of BSFL in relation to their genetics and performance [21]. It has not been explored yet to what extent genotype, larval age, diet, and the interaction of these factors shape gut bacterial communities of BSFL. This study provides to fill this gap by examining the influence of genetics, diet, and larval age in shaping gut bacterial communities, thereby potentially altering gut metabolic capabilities, and allowing BSFL to adapt to various plant- and meat-dominated substrates, which ultimately determines performance.

Genetics and environment influence the gut microbiota of BSFL

Studies on other insect species have shown that the gut microbiota are influenced by both genetics and environment (G × E interaction) [7]. This appears to be the case for BSFs as well. Differences between lines, diets, as well as larval age were found to significantly influence gut bacterial diversity and/or community structure. While domesticated insects are generally associated with decreased gut microbial diversity, potentially limiting their adaptability to changing environments compared to wild populations with ongoing horizontal microbiota transfer, wild populations have often been reported to possess a higher microbial diversity [33]. Here, the less-adapted WT line overall showed a lower α-diversity, however, this was largely driven by the D5 data, whereas no significant differences in α-diversity were observed between lines at D10. In fact, overall β-diversity was higher for WT, as demonstrated by its significantly higher average within-line Bray-Curtis dissimilarity compared to LD. In view of the high performance of NUS-fed LD larvae, their low bacterial richness may hint to a specialized adaptation of LD to this specific environment by retaining only those bacteria beneficial for its fitness [22, 33]. So far, only one study compared a plant-dominated substrate (Housefly Gainesville diet) to a meat-dominated substrate (hatchery waste) using a single line of BSFs and found a higher richness in larval gut bacteria on the meat-dominated substrate [2]. This contrasts with previous studies on mammalian and human hosts, in which intestinal bacterial diversity was generally reported to be higher in herbivores and hosts feeding on plant-dominated diets [34, 63]. In line with these works, our study seems to suggest that plant-dominated substrates foster a higher bacterial richness as compared to meat-dominated diets (based on data from two distinct BSF lines reared on ten different feed substrates available in Southeast Asia).

Members of several genera, namely Levilactobacillus and Morganella (LD) as well as Paenibacillus and Lacrimispora (WT) showed significant differences between lines irrespective of age or diet suggesting the potential impact of a genetic component on their prevalence and abundance in BSFs. While species of Levilactobacillus and Lacrimispora have not yet been specifically highlighted as potential closely-host associated symbionts, species of Paenibacillus form close associations with insect pathogenic nematodes of the genus Heterorhabditis [16]. While insect pathogenicity of some Morganella species cannot be ruled out, mutualistic relationships of Morganella with different insect hosts have also been observed [59]. Further investigation is needed to validate the mechanisms through which these taxa are influenced.

Expectedly, several taxa were found to be significantly correlated with larval age. In line with previous dietary studies on BSFL, a succession of bacterial communities over time was observed as a result of dietary adaptation [2, 46]. This needs to be considered when analysing and interpreting the potential roles of different bacterial taxa in larval performance.

In terms of diet composition, an optimal protein to carbohydrate ratio is reported to be most crucial in the development of BSFL [5]. Even though identical diets (with identical protein and carbohydrate content) were administered to both lineages, the divergence of the bacterial communities between LD and WT observed on D5 and D10 confirmed that G × E interaction modulates both lineages in their dietary adaptability, thereby forming specific bacterial communities that determine larval growth and development [35]. From the current data, we can conclude that different lines show different performance on the same substrate due to differential underlying gut bacterial communities. Therefore, new feedstocks should be tested with several and matched with the best-performing fly line. Further adaptation of gut-associated bacteria to the new feedstock can then lead to optimized growth and performance as previously shown for the adaptation of whiteflies to new host plant species [52]. Future studies employing metagenomics may reveal the biochemical pathways utilized by the gut microbiota of genetically distinct lineages that underlie performance differences when reared on the same substrate.

Correlation of bacterial taxa with larval performance

As established above, a multitude of factors including genetic variation, diet, and age influence the composition of gut bacterial communities in BSF larvae, and these in turn determine larval performance. We used the data of Zhang et al. [65] to identify bacterial taxa that were positively or negatively correlated with larval performance across the different plant- and meat-dominated substrates.

Although the contributions of Morganella to gut functioning in BSFL remains largely unexplored, certain members of both Morganella and Providencia are known to participate in urea hydrolysis wherein urea is broken down into ammonia which can be utilised by cells [19]. Moreover, incorporating urea or ammonia into the diets of BSFL has been observed to enhance larval biomass [48]. In addition, a prior investigation revealed a positive association between dietary fat and Providencia, implying that members of Providencia could play a significant role in lipid and protein conversion within the gut of BSFL [1]. Providencia has also been shown to be a dominant player throughout BSFL instars promoting growth and development (i.e., larval weight gain, development of prepupae, and eclosion) [35]. In this study, although Morganella and Providencia did not appear to associate with any of the dietary macronutrients, in LD, Providencia correlated positively with BCR, PW, protein, fat conversion, and prepupal fat content while Morganella with protein content. Higher relative abundances of Providencia were observed specifically in certain plant-dominated diets (as well as NUS). This was corroborated by previous findings, whereby Providencia has been shown to possess enzymes that enable the degradation of xylan present in plant cell walls [54]. In WT, only Providencia correlated positively with performance. Highest relative abundances of Providencia were observed in NUS fed larvae for both genetic lines suggesting a potentially higher release of ammonia from food waste could have assisted in their weight gain [28]. Confoundingly, Morganella was the only genus that positively correlated with prepupal protein content in LD while negatively correlating with WT. Studies on BSF have reported Morganella to induce immune responses [37, 43].While specific studies directly linking Morganella to altered protein storage in insects are limited, infections from pathogenic bacteria in general often induce stress responses and immune reactions that alter the host’s normal metabolic processes, including protein synthesis and storage [35]. Furthermore, based on the 16S rRNA gene data alone, we could not determine which Morganella species was or were responsible for the observed correlations. It is acknowledged that different species or strains within this genus may exhibit differential protein-accumulating capabilities, which warrants further investigation.

In addition to Providencia, Peptostreptococcus showed a positive correlation with several larval performance metrics in WT. Although specific studies on insects are lacking, human research indicates that Peptostreptococcus produces enzymes that break down complex proteins into simpler amino acids and metabolites for absorption [14]. Additionally, Peptostreptococcus species are involved in fermentation processes that produce short-chain fatty acids (SCFAs) like acetate, propionate, and butyrate, which are known to be crucial for energy, gut health, and metabolic processes [9]. While studies on BSF have shown the presence of Peptostreptococcus in rearing substrates [8, 61], its association with BSFL gut has not been explored. Given its positive correlation with BCR, PW, and protein conversion rate in WT larvae, this association may be a promising area for further exploration.

Sphingobacterium showed high relative abundances in larvae fed with PKM and RIB and was associated with dietary carbohydrate. It appears that Sphingobacterium thrives on these lignin-rich diets, possibly due to its ability to produce enzymes aiding in the breakdown of lignin [47]. Previously, a psychrotolerant strain of Sphingobacterium demonstrated lipase activity [29]. Thus, in the context of BSFs, Sphingobacterium has previously also been hypothesized to contribute to the decomposition of lipid-rich substrates such as plant and animal fats [30]. In the current study, Sphingobacterium was negatively correlated with prepupal fat content, with lowest relative abundances in larvae fed lipid-rich diets (NUS, P, and MLP; [65]).

As one of the suggested “core” genera, Proteus has previously been emphasised to play an important role in BSF oviposition [19], and has been shown to enhance BSFL biomass production [38]. On the contrary, Li et al. [35], showed that Proteus inoculation into the gut of BSFL led to a reduction in biomass gain. These seemingly opposite findings may be explainable by possible experimental or methodological differences such as the use of different BSF lineages, or different species or strains of Proteus. The positive correlation found between Proteus relative abundance and PW, BCR, as well as fat conversion rate in LD (but not WT) aligns with the observations of Mazza et al. [38]. Interestingly, a study conducted by Cifuentes et al. [12] demonstrated the inhibition of the growth of Proteus in the gut of BSFL by Bacillus. In our study, inhibition was not observed in LD; in fact, both Proteus and Bacillus positively correlated with PW and BCR in LD larvae. However, in WT, high abundances of Bacillus generally co-occurred with low abundances of Proteus, apart from a single sample in which Proteus showed a high relative abundance despite the presence of Bacillus. It may be conceivable that bacterial populations in highly-adapted laboratory fly lines housed under standardized conditions are similarly optimized as their hosts and adapted to co-exist, partition resources, or even co-operate, while antagonistic behaviour between populations may be more pronounced in wild BSF lines frequently subjected to a changing environment, overall resulting in a higher resilience of these communities [13, 18]. Moreover, in WT larvae, Proteus did not show any correlation with performance. This suggests that Proteus plays multifaceted roles in the physiology of its host, and harnessing its full potential for larval performance may depend on a variety of factors. Further research is needed to decipher the specific functions of this bacterial genus especially those pertaining to the performance of BSFL.

Scrofimicrobium is a relatively recently recognised genus [62]. However, since its discovery, it has been found in several studies on BSFL gut microbiota [1, 25] and has been suggested to function similarly to its close relative Actinomyces in its antibacterial properties and degradation of plant material through enzymatic activities [64]. However, in this study, a negative correlation between Scrofimicrobium and larval performance was found for both genetic populations. The frequent occurrence of Scrofimicrobium in BSFL studies warrants further investigation. It may suggest that BSFL vary in their microbial associations for nutrient acquisition depending on the administered diets [15].

Klebsiella is ubiquitous in the environment and commonly known as a food spoilage microbe [10]. It is also one of the most common genera found in the gut of BSFL [15]. It remains to be elucidated whether the strains found in the gut of BSFL are the same or different to those potentially present in certain diets. Amongst the species of the genus Klebsiella, Klebsiella oxytoca stands out for its pectinolytic activity [58]. Thus, previous studies on BSFL, which reported a notable prevalence of Klebsiella in larvae reared on substrates rich in cellulose and hemi-cellulose suggested a potential role in breaking down these complex dietary sugars [20]. Contrastingly, human and animal studies have shown that a high intake of complex dietary carbohydrate and fibre inhibited colonisation of the colon by Klebsiella pneumoniae [23]. In our study, Klebsiella did not associate significantly with either plant or meat-dominated diets; highest relative abundances were found in NUS (both lines), PKM (WT), and RIB (LD). While the genus exhibited no correlations in WT, it correlated positively with protein conversion rate in LD. Interestingly, in larvae fed NUS, elevated levels of Klebsiella co-occurred with optimal larval performance across various parameters in both LD and WT populations. In future, the use of sequencing technologies allowing for species or even strain-level resolution may help to disentangle these confounding results and to establish exactly which members of this genus are contributing positively or negatively to larval performance.

Common and lineage-specific core bacterial taxa

Core microbial taxa have been assumed to play important roles in BSF physiology. For example, the BSFL’s core gut bacteria together with the environmentally modulated microbial community may contribute to the larvae’s ability to adapt to nutritionally diverse diets. To date, there is no clear consensus as to the existence of a universal “core” due to the use of non-standardized methods to define “core” microbial taxa in different studies. In a recent meta-analysis, five core genera were reported despite none of them showing a prevalence of 100% [25]. This suggests that BSFs do not harbour a core microbiota in the strict sense of heritable symbionts that engage in a stable mutualistic relationship with their host, as observed for some insect taxa [15, 41]. Generally, we agree with using a strict prevalence threshold of 100% to define a core taxon, meaning that a taxon should occur in all samples to be recognised as a core taxon. However, there may be justifiable explanations for the apparent absence of a taxon from a sample, in which case a relaxation of the prevalence threshold or exclusion of the sample may be advisable. Valid reasons to relax the prevalence threshold or exclude samples from core analyses may include experimental (e.g., antibiotics treatment, presence of toxic compounds), methodological (e.g., DNA extraction, PCR, or sequencing bias against a certain taxon or under-sampling due to low sequencing depth, as evident from the coverage index), or biological biases (e.g., illness of individual larvae). The five “core” genera, as reported at the 80% prevalence threshold by IJdema et al. [25], were also observed at the > 80% prevalence threshold in the current study. Out of these, Enterococcus, Morganella, and Providencia were shared across > 97% of all samples and across 100% of samples for LD and > 95% for WT when distinguishing between lineages. These findings corroborate the important roles of these taxa as potential core microbiota in BSFL. For example, Providencia has been reported to be involved in the conversion of lipids and proteins in BSFL [1] and is positively correlated to oviposition in adults [55]. This has been further confirmed in our study where Providencia appeared to participate in several aspects if not all for larval performance for both genetic lines emphasising its importance in larval growth. Providencia and other potential core genera identified previously, although their respective functionalities are not yet fully understood, seem to be interconnected [55]. Lenient thresholds of as low as 50% prevalence as used in some studies can be problematic when aiming at delineating a truly common core microbiota. Instead, separating the data by meta data categories, such as genetic lineage, may reveal additional category-specific core taxa. Here, for the first time, we show that in a controlled rearing experiment, the genetic background of BSF influences the gut microbiota. The concept where host-associated microbial communities evolve in a pattern that mirrors host evolutionary relationships, or phylosymbiosis, describes the existence of lineage-specific core microbiota in various species such as Drosophila fruit flies, Aedes, and Culex mosquitoes, whereby the intraspecific transfer of microbiota within a lineage can either lead to an increase or decrease in host fitness and performance [6, 33]. Unfortunately, much is unknown regarding the evolutionary trajectory of BSFs and the functionality of potentially coevolved host-associated gut bacteria. The new findings presented in the current study on the diversity and structure of gut bacterial communities in BSFL as well as on the potential impacts of individual bacterial (core) taxa on larval performance, underscore the necessity for further research in this domain, in particular, uncovering the routes of transmission of the highly prevalent taxa. This future investigation is imperative to enhance our comprehension of the functional and potentially partially co-evolutionary relationship between BSFL and their gut microbiota.

Conclusion

We have exemplified that two genetically different populations of BSFL not only differed in their ability to adapt to various diets but appeared to perform differentially in plant/meat-dominated diets. Besides harbouring different genera between populations, the shifts of these genera in plant- or meat-dominated diets (G × E interaction) have not been previously reported. Apart from detecting several genus-level bacterial taxa that were shared at > 97% prevalence across two genetically distinct lineages, we also observed for the first time that additional lineage-specific core taxa existed, which showed high prevalence in one line but not in the other. We recommend that for any new study careful consideration should be given to potential biases within the data set prior to the delineation of core taxa. Unless there is valid reason for the relaxation of the prevalence threshold or exclusion of samples, a standardised high prevalence threshold of at least > 90% if not higher should be applied in future studies. Further exploration of the genetic factors influencing the gut microbiota and the potential role of certain genera for larval performance is crucial for understanding the mechanisms through which the host regulates its microbial community. This knowledge can pave the way for selective breeding strategies aimed at enhancing desirable traits related to waste decomposition efficiency, nutrient cycling, and resistance to pathogens.