Changes in rearing conditions rapidly modify gut microbiota structure in Tenebrio molitor larvae

The gut microbiota of multicellular organisms has been shown to play a key role in their host biology. In mammals, it has an invariant component, responsible for establishing a mutualistic relationship with the host. It also contains a dynamic fraction which facilitates adaptation in response to changes in the environment. These features have been well described in mammals, but little is known about microbiota stability or plasticity in insects. We assessed changes in microbiota composition and structure in a reared insect after a change in rearing conditions. We reared Tenebrio molitor (Coleoptera, Tenebrioninae) larvae for five days in soil samples from two river banks and analyzed their gut microbial communities by a metabarcoding technique, using the V3-V4 region of the 16S rRNA gene and the housekeeping gene gyrB. We found that soil-reared insects had a significantly more diverse microbiota than the control insects and that insects reared in soil from different sites had significantly different microbiota. We confirmed this trend by absolute quantification of the two mains fluctuating taxonomic groups: the Enterobacteriaceae family and the Pseudomonas genus, dominant in the soil-reared insects and in the control insects, respectively. Our results suggest the existence of a resident microbiota in T. molitor gut, but indicate that rearing changes can induce rapid and profound changes in the relative abundance of some of the members of this resident microbiota.

gut microbiota is a key component, with major effects on host physiology. For example, the mammalian gut microbiota has been the object of many studies on digestive functions with health implications (Belizário and Napolitano 2015).
The composition of the mammalian gut microbiota displays both plasticity and invariant features.
The core microbiota, which consists of the microorganisms common to the majority of individuals within a population, is generally defined as the most prevalent of the microbial species detected (Shetty et al. 2017). This common fraction of the microbiota plays a fundamental role in supporting the mutualistic symbiotic relationship with the host (Candela et al. 2012). For example, changes in the human core microbiota are associated with physiological perturbations, such as obesity and Crohn's disease (Turnbaugh et al. 2009;Hedin et al. 2015). However, another key feature of the mammalian gut microbiota is its plasticity, i.e. its ability to change in composition and structure. In humans, dietary changes induce a remarkable degree of variation in gut microbiota in terms of both phylogenetic and functional composition (Candela et al. 2012). These changes depend on various factors including host age, sex, genetic make-up, immune and health status (Shetty et al. 2017), but also exposure to environmental bacteria, geographic origin and climate (Candela et al. 2012). It has been suggested that this plasticity of the human gut microbiota facilitates rapid responses to environmental change, resulting in rapid ecological adaptation (Alberdi et al. 2016).
Most studies on the gut microbiota concern mammals. However, the use of mammals, and more generally of vertebrates, in experimental approaches raises numerous practical, financial and ethical issues. Large-scale experiments require model organisms that are easy to manipulate and can be obtained in large numbers. Insects are interesting experimental models in this respect. Although 2016;Welte et al. 2016;Shao et al. 2017;Raymann and Moran 2018). The gut microbiota of nonsocial insects is principally acquired from the environment through feeding (Engel and Moran 2013). Its composition depends on environmental conditions and diet in both laboratory and wild individuals (Chandler et al. 2011;Montagna et al. 2015;Staudacher et al. 2016). For example, it has been shown for some coleopteran species that microbiota changes with geographical location, environmental condition, and life stage (Huang and Zhang 2013;Montagna et al. 2014).
One potential limit of these previous studies is that they used either insects from the wild, which cannot be controlled for many of their characteristics, or lab-reared insects, which are controlled but have a poorly diversified microbiota. Here we used laboratory-reared T. molitor larvae and mimicked a soil environment by rearing the larvae for five days in different soil samples. We assessed the changes in gut microbiota composition after acclimatization to soil samples and demonstrated a large shift in gut microbial structure. We showed in addition that different soil samples induced different modifications in insect microbiota, and that the observed plasticity was probably dependent on changes in the abundance of some of the resident OTUs.

Soil samples
We sampled soil from riverside land around Montpellier in the South of France ( Figure 1A): on the banks of the Hérault river near Causse-De-La-Selle (N43°49.884' E003°41.222'; CDS sample) and on those of the Lez river near Montferrier-sur-Lez (N43°40.801' E003°51.835'; MTF sample). Both soils had a sand-silt-clay composition typical of riversides on chalky substrata. The sand:silt ratio was higher for MFT than for CDS. We collected three soil subsamples from each plot. These subsamples were taken at a depth of 20 cm and were separated by a distance of 10 m. They were named CDS1, CDS2, CDS3 and MTF1, MTF2, MTF3 ( Figure 1B). The use of these six soil subsamples made it possible to compare the variability in microbiota composition both between and within plots. Each soil subsample was split into four portions, each of which was placed in a 1 L plastic box ( Figure 1C), in which it was mixed with heat-sterilized (20 min at 121 °C) wheat bran (1:3 (v/v) ratio, as previously described (Jung et al. 2014).

Insects
Larvae were provided by Micronutris (St-Orens, France) and fed with heat-sterilized bran before the experiment. As it was not possible to determine their precise developmental stage, but we used only larvae weighing between 20 and 50 mg, which should correspond to 13th or 14th instar individuals (Huang et al. 2011).

Rearing of T. molitor larvae in soil samples
We maintained laboratory-reared T. molitor larvae for five days in sterilized wheat bran mixed with soil samples. During this period, the larvae were incubated at 15 °C in the same humidity conditions. They were then starved for 24 hours ( Figure 1D) to exclude individuals that were infected with pathogens (which would have died within this 24 hours period) and to limit the risk that the DNA we extract comes from the larval alimentary bolus. This starvation period potentially induces a stress on insect larvae, which might in turn impact their microbiota. We imposed it on all insects, so that the potential bias it creates is identical in all treatments.
Control insects were reared in the same conditions than other insects except that they were incubated in sterile wheat bran, with no soil mixed. Control insects microbiota should thus be close to what it was for all insects before the experiment. per site) and 5 control insects. However, we failed to amplify 16S rRNA during PCR step for 2 samples, ending with 24 samples for CDS, 22 samples for MTF and 5 controls. Insect larvae were sterilized in 70% ethanol, rinsed in water and then killed. The guts of the larvae were dissected in sterile Ringer solution (Jung et al. 2014). Dissection tools were sterilized with 70% ethanol between insects. Dissected guts were placed in an Eppendorf tube with 100 µL of lysis solution and 1 µL lyzozyme (Quick Extract, Bacterial DNA extraction TEBU-BIO) and ground with 3 mm steel beads for 30 seconds at 20 Hz with a TissueLyzer (Qiagen). The resulting homogenates were incubated at room temperature for two days, then frozen in liquid nitrogen and heated at 95 °C to ensure that all the cells were lysed. DNA was prepared by the phenol-chloroform-alcohol and chloroform extraction method. The DNA was resuspended in sterile water and quantified with a NanoDrop spectrometer (Thermo Fisher Scientific). We performed extraction blanck controls using DNA-free water.

16S and gyrB DNA amplification
We targeted the V3-V4 region of the 16S rRNA gene, which is classically used for bacterial identification in microbial ecology studies, as clean and complete reference databases are available for this region. We also used the bacterial housekeeping gene gyrB, to support the data for the 16S rRNA (Barret et al. 2015). The V3-V4 region of the 16S rRNA gene was amplified with the PCR1F_460 (5'-ACGGRAGGCAGCAG-3') / PCR1R_460 (5'-TACCAGGGTATCTAATCCT-3') primers (modified versions of the primers used in a previous study Klindworth et al. (2012)).
Amplification was performed with the MTP Taq polymerase (Sigma,, according to the manufacturer's protocol, with 1 µL of 1/10 diluted DNA extract for each sample. The PCR protocol used for these primers was 60 s at 94 °C, followed by 30 cycles of 60 s at 94 °C, 60 s at 65°C , 60 s at 72 °C, and then 10 min at 72 °C. The gyrB gene was amplified with primers described elsewhere: gyrB_aF64 5'-MGNCCNGSNATGTAYATHGG-3' and gyrB_aR353 5'-ACNCCRTGNARDCCDCCNGA-3' (Barret et al. 2015

Metabarcoding data treatment
Sequence data for both markers were analyzed with OBITools (Boyer et al. 2015). Raw paired-end reads were aligned and merged, taking into account the phred quality score of each base to compute an alignment score. Reads with a low alignment score (>50), containing unknown bases or with an unexpected size (outside 400 bp and 470 bp, and 230 bp and 260 bp after primer trimming for the 16S rRNA gene and gyrB respectively) were removed from the dataset. After primer trimming, singletons (i.e. sequences only found once in the dataset) were removed (Auer et al. 2017 Barret et al. (2015) for gyrB gene (OTU assignments are available in Additional File 5).

Quantitative PCR analysis
To check for changes in OTU abundances, we performed quantitative PCR (qPCR) on two randomly picked insects per soil subsample among those used in the metabarcoding analysis. The sampling probability for each sample was adjusted for the total number of 16S rRNA reads for the sample. The five DNA samples corresponding to control insects were all analyzed.
All qPCRs were performed in triplicate, with 3 µL of reaction mixture, on a LightCycler480 bp mol -1 ) and the mean weight of a double-stranded base pair M b p (660 g mol -1 bp -1 ) as follows:

GC N t a r g e t = N A × M D N A L D N A × M b p ×n t a r g et
Using the parameters of the curves linking GC N t a r g e t and C t in standard samples, we then estimated the GCN of target genes in our gut samples. This estimation was possible because PCR efficiency (PE) was very close to that for standard samples (Additional File 6).

Community analysis
All analyses were performed with R version 3.3.3 (R Core Team 2015) (see Additional File 7 and 8 for the overall workflow). We did not rarefy data (McMurdie and Holmes  We calculated the beta diversity distance matrix from the Jaccard and Bray-Curtis distances for presence/absence and relative abundance data, respectively, using the vegan package. We also computed Unifrac and wUnifrac distances for presence/absence and relative abundance data,

Incubation of T. molitor larvae with soil increases the richness and diversity of their gut microbiota
After cleaning, the total dataset of the metabarcoding experiment contained 792,395 sequences insects (MTF and CDS) was a 48 ± 13 OTUs whereas that of control insects (BRAN) was 25 ± 9 OTUs ( Figure 2B). The OTU richness of the gut microbiota therefore increased significantly after the incubation of the insects with soil samples (Chao1 index, soil vs. control: Wilcoxon rank sum test, W=221, p-value = 1e-3). A similar conclusion was drawn for analyses based on the Shannon index, which reflects relative OTU abundance within samples ( Figure 2B, soil vs. control: Wilcoxon rank sum test, W=216, p-value = 1e-3). Moreover, control insects harbored bacterial communities dominated by a very small number of dominant OTUs (low Shannon index ≃ 0.2 and low Pielou's eveness ≃ 0.02). OTU assignment identified these dominant OTUs as belonging to the Pseudomonadaceae family ( Figure 2C). By contrast, soil-reared insects harbored bacterial communities with more balanced relative OTU abundances (Shannon index ≃ 1.7).
The gut microbiota of these insects was dominated by Enterobacteriaceae, together with Pseudomonadaceae and other less frequent families, such as Moraxellaceae and Aeromonadaceae ( Figure 2C). This was confirmed by the analysis of Pielou's eveness index which was significantly lower in control insects than in soil-reared insects (Wilcoxon rank sum test, W=0, p-value = 7.6e-7).
Thus, five days in soil significantly increased the richness of the microbiota in the gut of T. molitor larvae, and modified the balance of OTUs present.
We also investigated the effect of soil treatments according to soil origin, by comparing the alpha diversity of CDS and MFT samples. The Chao1 and Shannon indices were significantly lower in MTF than in CDS samples ( Figure 2B; Chao1

Soil treatment induces a change in microbiota composition that is variable between soil sampling sites
We investigated the effect of soil treatment on insect microbiota, by calculating the beta-diversity between insect gut microbiota with various distance indices (Figure 3). We first calculated a distance based on pairwise Jaccard and Bray-Curtis distances. These two indices are complementary, because Jaccard distance depends purely on the presence/absence of OTUs, whereas Bray-Curtis distance also takes into account the number of reads for each OTU as a proxy for their relative abundance. We performed PCoA analysis on distance matrices ( Figure 3A) where control insects tended to cluster together. PERMANOVA analysis confirmed that community composition differed between soil-reared insects and control insects (13 to 19% of the variance explained by soil treatment, Table 1A).
The microbiota profiles of insects placed in soils from the same site (i.e. CDS or MTF) or in the same soil subsample (e.g. CDS1, CDS2 or CDS3) did not cluster together perfectly. However, a second PERMANOVA model for these samples identified two explanatory factors, soil sampling site (i.e. CDS or MTF) and subsample identity (e.g. CDS1, CDS2 or CDS3), as having a significant impact on gut community composition (Table 1B). Indeed, sample site explained 14 and 8% of the variance and soil subsample explained 17 and 18% of the variance, for the Jaccard and Bray-Curtis indices, respectively.
As reported above, the soil-reared insects had a microbiota dominated by Enterobacteriaceae ( Figure 2C). We thus estimated Unifrac distances, which take into account the phylogenetic distances between OTUs, and wUnifrac distances, which also take relative OTU abundance into account. With these corrections, the differences between control insects and soil-reared insects were significant only when relative OTU abundance was taken into account ( Figure 3; Table 1A). Subtle but significant effects of sample site and soil subsample on community composition were also observed with the Unifrac and wUnifrac indices (Figure 3; Table 1B).
Overall, our results show that soil treatment changes the community composition of the gut microbiota and that this change is detectable despite inter-individual variability. The bacterial communities present in the gut differ both between sample sites and between soil subsamples. Most of the changes in the microbiota concern the relative abundances

of OTUs
We then pooled all individuals of a given treatment to determine which OTUs are found in at least one individual for each treatment. The 47 OTUs found in control insects were also present in the insects of the soil treatment groups ( Figure 4A). The 44 OTUs common to all three conditions matched 97% of the reads for soil-reared insects (gray area in Figure 4B and Figure 4C). However, after soil treatment, Pseudomonas, the dominant OTU in control insects (98% of the reads) accounted for only 27 and 23% of the reads in CDS and MTF samples, respectively ( Figure 4C).
Conversely, Serratia species, together with the Enterobacter group, which accounted for less than 1% of sequence reads in controls, accounted for 35% and 43% of the reads for CDS and MTF, respectively.
For confirmation of our initial metabarcoding results, we performed a second metabarcoding analysis with another marker, a 300 bp region of the gyrB housekeeping gene (see Additional File 1). This single-copy marker has been shown to provide assignments to more precise taxonomic levels than the 16S rRNA gene (Barret et al. 2015). In accordance with the results obtained with the 16S rRNA gene marker, Pseudomonas was the dominant OTU in control insects (more than 99 % of the reads) and its relative abundance was lower in soil-reared insects (CDS: 14 % MTF: 17 % of the reads). The genus Serratia and the Enterobacter group accounted for less than 0.06 % of the reads in control insects and a large proportion of those for the insects in the two soil treatment groups (CDS: 57 % MTF:70 % of the reads).
Finally, we also identified with 16S rRNA 59 OTUs that were not detectable in control insects but were present at low abundance (3% of the reads) in at least one soil-reared insect (red dashed area in Figure 4B and Figure 4D). These OTUs may correspond to taxa that were absent from the insects before soil treatment, and that colonized the insect gut during incubation in soil. Alternatively, they may have been present in the control insects at densities below the PCR detection threshold. Their abundance would then have increased above this threshold during incubation, just like the abundances of Serratia or Enterobacter. Overall, our data strongly suggest that the main effect of soil treatment is a change in the relative abundances of OTUs, although low levels of bacterial colonization from soil cannot be ruled out.

The balance between members of the resident OTUs contributes to the variation of abundances after soil treatment
We assessed the variation of OTU balance after soil treatment further, by quantifying the bacterial taxonomic groups present in all treatments but with different relative abundances between the two contrasting sets of conditions studied (control versus soil-reared). We first characterized the gut resident microbiota in our larvae, as the OTUs present in at least 95% of our samples (following Wallis rank sum test, chi squared = 2.66, df = 2, p-value = 0.26), which suggests that the total number of bacteria was similar in our 17 samples. We then targeted a region of the 16S rRNA gene specific to the Pseudomonas genus, (Pse -16S: 251 nucleotides of the V3-V4 hypervariable region, with 4 to 7 copies per genome Bodilis et al. 2012), and a region of the rplP gene, region specific to the Enterobacteriaceae family (Entero-rplP : 185 nucleotides of the rplP gene, one copy by genome). The Pse -16S GCN in soil-reared insects was one tenth that in control insects ( Figure 5A).
Conversely, the Entero-rplP GCN was 100 times higher in soil-reared insects ( Figure 5B). Soil acclimation therefore seems to induce a decrease in Pseudomonas abundance and an increase in Enterobacteriaceae abundance. Our data suggest that the main effect of soil treatment is to modify the relative abundances of the resident bacterial communities of the gut microbiota.

Discussion
Rearing larvae in soil rather than in bran caused major changes in gut microbiota structure. Soilreared larvae have a richer and more diverse gut microbiota than control larvae. Despite considerable inter-individual variability, we found that the changes in community composition depended on both the site from which the soil was obtained, and the precise soil subsample used.
An analysis of the OTUs found in the different samples suggested that the main effect of the soil treatment was a change in the relative abundance of OTUs. We confirmed this trend by qPCR for the two main taxonomic groups displaying changes in abundance: the Enterobacteriaceae family and the genus Pseudomonas, which predominated in soil-reared insects and in the control, respectively.
Our rearing conditions (laboratory versus soil acclimatization) were associated with two types of gut microbial patterns, consistent with previous findings for laboratory-reared and wild insects. On the one hand, gut microbiota communities of laboratory-reared insects, which are usually intestinal tract showed differences in pH, redox potential and hydrogen contents, and were associated to different bacterial communities (Bauer et al. 2015). The ingestion of soil particles probably modifies some of these properties of the gut. The fact that the soil characteristics differed between the two sampling sites (low sand/silt ratio for Causse-De-La-Selle (CDS), and higher sand/silt ratio for Montferrier (MTF)) could thus explain in part their different impacts on T. molitor gut microbiota.
The changes in the gut bacterial population may depend not only on treatment, but also on the bacterial community initially present in the gut. Previous studies (Jung et al. 2014;Osimani et al. 2018) showed that a Spiroplasma species predominated in the gut microbiota of the larval lineage, even after and environmental change. Spiroplasma has been shown to be a heritable endosymbiont in Drosophila (Mateos et al. 2006). Similar effects were observed for other endosymbionts, such as

Wolbachia, Cardinium, Blattabacterium-like and putative Bartonella-like symbionts in mites
Tyrophagus putrescentiae following dietary changes (Erban et al. 2017). In all these case, endosymbiont seem to impede major shifts in the gut microbiota or conceal changes in frequencies that may occur in low-abundance OTUs. This effect is absent in our experiment, probably because the insects we used are associated to Spiroplasma or any other endosymbiotic bacteria.
Our results also provide interesting insight into the spatial variation of the gut bacterial community in insect populations. The differences observed after incubation in soil from different plots were consistent with the findings of other studies on coleopterans, in which the dissimilarity of the gut bacterial community between individuals is correlated with the distance between sampling sites (Adams et al. 2010). However, we also observed a difference in the gut microbiota between insects incubated with soils collected a few meters apart, at the same sampling site, and this difference was detectable despite high levels of inter-individual variation. Minor environmental differences thus have a detectable impact on the gut microbiota and structure this microbiota within insect populations over very small geographic scales.     Venn diagram of OTUs found in at least one insect from each treatment. B. Bar plot of the relative abundance of the 44 OTUs common to the three treatments (in gray) and the 59 OTUs found only in soil treatments (CDS and MTF) (red stripes). The taxonomic assignment of these OTUs is detailed in C. and D.. Insects from the various treatment were pooled for these bar plots: 5 insects for BRAN, 24 insects for CDS and 22 insects for MTF. The relative abundance of OTUs was calculated from the total number of reads for each insect pool. We show here taxonomic assignments to genus level or to the lowest taxonomic level, for which the bootstrap score was < 80%. Some OTUs differ in sequence, but were assigned to the same taxonomic group. These sequences are differentiated by a number. On each graph, the 15 OTUs with the largest relative abundance are shown in color and the others are grouped together in the "Others" category. OTU names followed by a star (*) belong to the Enterobacteriaceae family. Insects from the various treatments were pooled for these bar plots: 5 insects for BRAN, 24 insects for CDS and 22 insects for MTF. The relative abundance of OTUs was calculated from the total number of reads for each insect pool. We show here taxonomic assignments to genus level or to the lowest taxonomic level for which the bootstrap score was > 80%. Some OTUs differ in sequence but were assigned to the same taxonomic group. These sequences are differentiated by a number.

Tables
On each graph, the 15 OTUs with the largest relative abundances are shown in color and the others are grouped together in the "Others" category. OTU names followed by a star (*) belong to the Enterobacteriaceae family.

Additional file 2: Example of a microbiota pattern in PCR replicates
We checked the reproducibility of PCR, by performing three technical PCR replicates (the three bars of the chart) on a sample chosen at random, with the whole metabarcoding procedure performed separately for each replicate. We show here the results for the CDS1D3 sample.

Additional file 3: OBITools workflow for 16S rRNA analysis
RMD_OBITools_workflow_V3V4.pdf and RMD_OBITools_workflow_gyrB.pdf contain OBITools, bash and R scripts used to obtain the OTU abundance table from raw sequencing data for both the 16S rRNA and gyrB genes.

Additional file 7: Statistical analysis workflow
RMD_R_workflow.pdf contains R scripts used to perform statistical analysis and to produce the figures presented in this study.