Regulation of host phenotypic plasticity by gut symbiont communities in the eastern subterranean termite (Reticulitermes flavipes)

Termites are eusocial insects living in colonies composed of three main caste phenotypes (workers, soldiers and reproductives), each with their own well-defined socio-behavioral roles (Lainé and Wright, 2003; Roisin, 2000; Wilson, 1971). Termite workers further have the unique characteristic capacity to molt into other caste phenotypes through the process known as caste differentiation, which is a form of phenotypic plasticity. Caste differentiation and development in Reticulitermes flavipes can follow either apterous/neuter or imaginal routes (Scharf et al., 2003) and can have three different developmental potentials that include status quo molts into next-instar workers, differentiation into a presoldier that later molts into a soldier (Zhou et al., 2006a), or differentiation into a supplementary reproductive individual. Research on the worker–soldier molt process is of recent interest as a non-destructive form (soldier) develops from the most economically destructive life stage (worker). Several endogenous and exogenous factors have been reported to regulate the worker–soldier molt. At the endocrine level, the sesquiterpenoid juvenile hormone (JH) acts as a hormone to differentiate workers into soldiers via an intermediate presoldier phenotype (Tarver et al., 2009; Zhou et al., 2007). In R. flavipes and most other insects, JH III is produced by a neuroendocrine gland known as the corpora allata (CA) closely associated with the brain (Yagi et al., 2005). The entire soldier differentiation process occurs in response to the elevated endogenous JH titers (Chan et al., 2011; Elliott et al., 2009; Elliott and Stay, 2008; Mao et al., 2005; Park and Raina, 2004). However, in other insects, elevated JH retains the immature phenotype through development resulting in extra-nymphal or larval instars (Nijhout et al., 2014; Smykal et al., 2014). In addition to regulating termite polyphenism, JH is also involved in vitellogenesis through inducing increased vitellogenin protein, insect diapause and other processes (Engelmann, 1984; Scharf et al., 2005; Yagi, 1976).

Besides JH, termite caste differentiation could be associated with their gut microorganisms. The termite gut is a habitat for numerous microorganisms from various taxa, including protists, Bacteria and Archaea, living in a consortium and contributing to termite physiology (Brune and Ohkuma, 2010; Ohkuma, 2008; Scharf and Tartar, 2008). Termite symbionts are crucial in cellulose/hemicellulose degradation, nitrogen fixation, acetogenesis and anti-fungal defense (Doolittle et al., 2008; Inoue et al., 2000, 1997; Peterson and Scharf, 2016a; Peterson and Scharf, 2016b). The literature is flooded with reports on the impact of gut symbionts on endocrine regulation in humans, mice and insects. The human gut bacterium Lactobacillus reuteri was reported to upregulate the neuropeptide hormone oxytocin to enhance wound-healing properties through nerve-mediated pathways (Poutahidis et al., 2013). Studies on human and mouse models have supported the conclusion that gut symbionts can affect immune and endocrine modulation in a variety of signaling pathways (e.g. Chen et al., 2015). Similarly, in some insects, gut symbionts were observed to influence the pheromones released, social behavior, host development, egg production, expression of hemolymph proteins, and hormones produced (Engl and Kaltenpoth, 2018; Lee et al., 2017). In honeybees, gut symbionts were also reported to increase body mass by mediating changes in insulin and vitellogenin signaling pathways (Zheng et al., 2017).

Several aspects of bacteria–termite symbiosis are yet to be understood. No reports are yet available on the interaction of termite gut symbionts with presoldier differentiation, which is regulated by JH produced by the endocrine system. In female bean bugs (Riptortus clavatus), the gut symbiont Burkholderia was reported to modulate the JH titer (Lee et al., 2019). Association with commensal bacteria such as Lactobacillus plantarum in Drosophila melanogaster has also been correlated with hormone signaling, providing strong evidence for the influence of gut bacteria on insect endocrine functioning (Storelli et al., 2011), while in termites, JH was reported to alter the expression of hundreds of symbiont genes in the gut (Sen et al., 2013). Another study has demonstrated that presoldier phenotypes have biased expression of genes coding for vitellogenin, a protein thought to bind JH (Scharf et al., 2003), while in other insects, gut symbionts were reported to influence the expression of vitellogenin (Zheng et al., 2017). Several studies on other organisms have further demonstrated that gut symbionts can impact hormone signaling pathways.

The goal of this study was to investigate whether termite gut bacteria can influence JH-dependent presoldier differentiation in R. flavipes, with the overarching hypothesis being that specific gut symbionts are connected to the presoldier differentiation process. Our specific objectives were to (i) study the effects of JH and antibiotic treatments on presoldier emergence and/or survivorship and protist numbers, and (ii) seek to connect specific gut bacteria with the above conditions using 16S rRNA gene sequence surveys of the gut microenvironment. Our findings provide a novel example of the influences of gut microbial symbionts on host phenotypic plasticity.

This study is based on a preliminary experiment where JH III and kanamycin were used to induce presoldier differentiation and results were similar to those obtained in the present study. We opted to use methoprene in place of JH III because it is consistently available with high purity, whereas JH III is not. Termites [Reticulitermes flavipes (Kollar 1837)] from colonies reared in the laboratory were used in bioassays. There were two phases for this experiment. In the first phase, paper towel sandwiches were pretreated with either kanamycin (50% w/v in water) or water alone, and were then fed to replicate groups of termites for 2 days. For the second phase, termites were then moved to paper towel sandwiches that were treated with either acetone or methoprene (30 µg per dish) in acetone (Chem Service Inc., West Chester, PA, USA) and dried in a fume hood prior to use. Methoprene and kanamycin were used to induce soldier caste differentiation and manipulate gut microbial load, respectively. Based on previous studies, kanamycin significantly reduced and altered the gut microbial load with minimal impacts to the host (Peterson et al., 2015), while exogenous application of methoprene has been reported to induce the genes in JH signaling pathways and similarly induce presoldier differentiation (Du et al., 2020; Tarver et al., 2012). Also, preliminary experiments with JH III demonstrated similar presoldier differentiation as seen in the present study with methoprene. Phase 1 pre-treatments lasted 2 days, which was identical for all the experiments; however, the duration of the final treatment was dependent on the type of experiment, as detailed in the following sections. Treatment abbreviations are as follows: Kan (kanamycin alone), Meth (methoprene alone), KanMeth (kanamycin in phase 1 and methoprene in phase 2) and UnCtl (untreated control), which received water in phase 1 and acetone (evaporated) in phase 2.

The assessment of presoldier emergence and survivorship was done to observe the differences among treatments at the whole-insect ‘holobiont’ level. There were four biological and three technical replicates each for this assessment (12 replicates total). For this bioassay, 25 worker termites were provided with kanamycin- or water-treated paper towel sandwiches in Petri dishes for 2 days (phase 1). After 2 days, 20 termites were transferred for final treatments as detailed above (phase 2) and observed for 15 days. Termites emerging from molting with transparent elongated mandibles were scored as presoldiers. Termites that were immobile and could not move their appendages when disturbed by soft forceps were considered dead. Survivorship in untreated controls, even in smaller groups, is quite high and invariable. We have tried working with larger groups in the past, and treatment-induced mortality is still a problem, coupled with an inability to score assays quickly and without causing additional, confounding assay stress. Therefore, we used small groups of workers for this experiment.

Each bioassay treatment described above was replicated 6 times and 35 termites were placed in each Petri dish for pre-treatment (phase 1). For final treatments (phase 2), 30 termites were transferred to Petri dishes with methoprene- or acetone-treated paper towel sandwiches for 9 days, which was the day of emergence of at least one presoldier in the methoprene treatment. The hind guts from five live termites each from assay days 0, 2, 6 and 11 were dissected. The guts were squeezed gently into 500 μl phosphate buffered saline (PBS, pH 7.4). A 10 μl suspension was used to count the protist numbers using a phase contrast microscope at 20× magnification and a Sedgwick Rafter counting cell (SPI Supplies, West Chester, PA, USA).

Bioassays with 25 worker termites were set up for DNA extraction. There were four replications each for the treatments. Similar to the phenological study, 20 termites were taken for the second phase of the experiment. Live termite workers were sampled and dissected after at least one presoldier was sighted (on day 11) in methoprene bioassay treatments. We picked day 11 for sampling DNA because the survival of termites was high and only live termite workers were dissected to acquire sufficient high-quality DNA from gut tissues. Also, a histological study of fat body (Cornette et al., 2007) revealed that the production of fat body and protein granules starts within the first 3 days after hormonal treatment. These granules disappear within a few days of presoldier development, indicating that termite physiological changes can be documented up to the point of presoldier molting. Four replications of 10 guts per treatment were sampled for DNA extraction. The 10 guts from live termite workers were dissected, placed in 200 μl PBS and processed for DNA extraction using the DNeasy^® Blood and Tissue Kit (Qiagen, Valencia, CA, USA). Along with the experimental samples, water was also run through the DNA extraction protocol as an extraction control. The quality of DNA in the aliquots was assessed through gel electrophoresis and was quantified using a Nanodrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and later verified using fluorometric detection by the University of Minnesota Genomics Core. Thirty microliters (41.3–125.9 ng µl⁻¹) of extracted DNA was sent for library preparation and sequencing, and the remaining DNA was used for 16S rRNA gene quantitative PCR (qPCR) analysis.

The V3–V4 region of 16S rRNA gene was sequenced using bacterial universal primers 341F-CCTACGGGAGGCAGCAG and 806R-GGACTACHVGGGTWTCTAAT (Caporaso et al., 2011; Peterson et al., 2015). A two-step PCR was used for library preparation. The indexed PCR products were normalized to 10 ng µl⁻¹, and 10 µl of each sample was pooled and concentrated to ∼100 µl and cleaned using solid reverse phase immobilization (SPRI) purification. The cleaned library was quantified by Qubit and fragment analyzed by Agilent Tapestation, and the library was diluted to 2 nmol l⁻¹. The pooled samples were denatured with NaOH, diluted to 8 pmol l⁻¹ in Illumina HT1 buffer, spiked with 15% PhiX and heat denatured at 96°C for 2 min immediately before loading. A MiSeq 600 (1/8th lane Stowaway) cycle v3 kit was used to sequence the sample by the University of Minnesota Genomics Center, MN, USA.

The paired end (301 bp) Illumina reads were processed and analyzed using the QIIME2 pipeline (Bolyen et al., 2019). Divisive Amplicon Denoising Algorithm 2 (DADA2) was used to obtain high-quality reads following the filtering, trimming, denoising, dereplicating, merging and removing chimera steps (Callahan et al., 2016). During this process, a total of 1,559,624 demultiplexed sequences were processed to obtain 1,234,554 non-chimeric high-quality sequences in QIIME2. The moving picture tutorial in QIIME2 was followed to calculate the diversity matrices. To calculate the diversity matrices (alpha and beta), the sequences were subsampled from the feature table from phylogenetic diversity analysis at the sampling depth of 2000 sequences that removed one replicate from KanMeth and four replicates of extraction control, respectively. The alpha diversity measures used were observed operational taxonomic units (OTUs), Pielou’s evenness and the Shannon index. The beta diversities, phylogenetic Unifrac distances [weighted and unweighted (Lozupone et al., 2011)] and non-phylogenetic distances [Jaccard and Bray–Curtis] were calculated and visualized using R and QIIME2. The reference database SILVA (version 132) was customized following the instructions on training the classifier for the V3/V4 region of 16S rRNA gene. Taxonomy was assigned to the data using naïve Bayes trained Silva 13,299% OTU classifier bounded by the 341F/806R primers set as the reference reads (Quast et al., 2013). Bar graphs representing the relative abundance of amplicon sequence variant (ASV) at the phylum and family levels were generated using the R package phyloseq from QIIME2-generated phylotype tables and taxonomy files (McMurdie and Holmes, 2013). Finally, analysis of composition of microbiomes (ANCOM) was used to study differential microbiota composition at phylum and genus level in different treatments via the QIIME2 plugin (Mandal et al., 2015).

A final trophallaxis experiment was conducted to understand the effect of different treatments in a eusocial context on donor and recipient termites. Prior to the experiment, ‘donor’ termites were treated with Nile Blue dye 0.5% w/v (Sigma-Aldrich, St Louis, MO, USA) for 2 days. After 2 days, the donor termites were transferred into Petri dishes containing filter paper sandwiches treated with methoprene or kanamycin, with an exposure duration of 2 days. The methoprene and kanamycin treatments and the untreated controls were set on the same day. Similarly, some untreated termites were exposed to kanamycin in different Petri dishes; these termites later served as ‘recipients’ for methoprene-treated donors. The treated termite donors were then placed in new Petri dishes with untreated or kanamycin-treated recipients in the ratio of 1:1 with a total of 20 individuals. The donor–recipient combinations during the experiment were: kanamycin-treated donors+untreated recipients (Kan-D+Untreated-R), methoprene-treated donors+kanamycin-treated recipients (Meth-D+Kan-R) and methoprene-treated donors+untreated recipients (Meth-D+Untreated-R). Controls for this experiment included kanamycin (Kan control), methoprene (Meth control) and untreated (Untreated control). Owing to the total mortality of termites in some treatments, the experiment was ended after 11 days. The data on termite mortality and presoldier initiation were taken every day until 11 days, with the cumulative data from day 11 of the treatment being presented. Each treatment combination was replicated three times.

Statistical analyses were performed in the QIIME2 interface and R (Bolyen et al., 2019; https://www.r-project.org/). The Kruskal–Wallis test was used to compare presoldier emergence, total mortality, average number of protists and alpha diversities among the treatment groups and trophallaxis experiment. Pairwise comparisons were made using the Mann–Whitney U-test and were corrected for Benjamini–Hochberg false discovery rate. Beta diversity matrices were subjected to permutational multivariant analysis of variance (PERMANOVA) and the PERMDISP test of homogeneity of dispersions to study statistical significance among groups. The graphs for 16S rRNA gene sequencing data were generated and visualized using R and Microsoft Excel. Canonical correspondence analysis (CCA) (ter Braak, 1986) was conducted to study the relationship between the bacterial genera and other variables – protist number, presoldier initiation and total mortality. CCA was implemented with the R package ‘vegan’. For the CCA analysis, four replicates from the presoldier initiation, total mortality and protist numbers were used to study their association with the taxa abundance.

The effect of treatments on cumulative percentage of presoldiers initiated (H=37.625, d.f.=3, P<0.001), cumulative percentage mortality (H=27.027, d.f.=3, P<0.001) on day 17 and average protist numbers per gut (H=19.048, d.f.=3, P<0.001) on day 11 were significant based on Kruskal–Wallis tests. Based on pairwise Mann–Whitney comparisons (P<0.05), the total percentage of presoldiers emerged was significantly higher in the Meth than in the KanMeth treatment, while no presoldiers emerged in the UnCtl and Kan treatments (Fig. 1A). From the same bioassays, survivorship of workers+presoldiers was also recorded. Total mortality in the KanMeth treatment was significantly higher than in any other treatments, and as expected, the UnCtl treatment had the lowest mortality. Mortality levels in the individual Kan and Meth treatments were not significantly different (Fig. 1B). Some variability in the mortality of termites from different replicates in the Kan, KanMeth and Meth treatments was observed, but this variability was not sufficient to overcome the statistically significant trends as reported. The variability could be due to the variation in group response to the treatments. Similar group variation in the mortality of termites against JH analog was observed by Hrdý et al. (2006) while treating different termite colonies with the same compound. The average number of protists per termite gut on day 11 was significantly higher in the UnCtl treatment in comparison with other treatments, while the Kan and KanMeth treatments had similar protist numbers (Fig. 1C). The average number of protists in the Meth treatment was significantly higher than in the Kan and KanMeth treatments, but still was greatly reduced in comparison with the UnCtl treatment. The average mortality of termites in all treatments on the day of protist and DNA sampling (day 11) was less than 50% (Fig. S1).

The two independent primer sets, 338F–518R and 515F–806R, used for qRT-PCR showed highly similar results, indicating that the Kan treatment had an expected result of reducing bacterial loads (∼2.5× decrease compared with UnCtl), but also unexpectedly that the Meth (∼1.2× decrease compared with UnCtl) and KanMeth treatments (∼2× decrease compared with UnCtl) had reduced bacterial loads as well.

Rarified 2000 sequence reads obtained from 16S rRNA gene sequencing were used to calculate alpha and beta diversity measures. Alpha diversity was estimated using observed OTUs for species richness (Fig. 2A), Pielou’s evenness (Fig. 2B) and the Shannon diversity index (Fig. 2C). The effect of different treatments on observed OTUs (H=12.829, d.f.=3, P=0.005), Pielou’s evenness (H=12.567, d.f.=3, P=0.005) and the Shannon index (H=12.567, d.f.=3, P=0.005) were significant based on Kruskal–Wallis tests. Based on pairwise Mann–Whitney comparisons (P<0.05), the UnCtl treatment demonstrated significantly higher species richness and evenness, followed by the Meth treatment, whereas the Kan and KanMeth treatments were not significantly different from each other. Similar results were obtained for average numbers of protists estimated per gut (Fig. 1C), suggesting a close connection between some bacteria and protist symbiont communities.

PERMANOVA statistics show that the weighted UniFrac analysis was statistically significant (pseudo-F=36.101, P=0.001), which means the effect of treatments was significant in creating the dissimilarities among the different treatment groups (Fig. S2). Also based on the PERMDISP test, the bacterial community structures among the replicates were not determined to be significantly different (PERMDISP test, F=0.7168, P=0.224). Pairwise PERMANOVA results further indicate all treatments differed in terms of community structure for weighted Unifrac. Similar results were observed for beta diversity matrices (Bray–Curtis and Jaccard). However, for the unweighted Unifrac analysis, the Kan and KanMeth treatments were not significantly different, whereas other pairwise comparisons were significant.

The taxonomic composition and relative abundance of major bacterial phyla at abundances greater than 0.1% varied greatly among treatments but was consistent among the replicates within treatments (Fig. 3). At the phylum level, in the UnCtl treatment, Spirochaetes, Firmicutes and Bacteroidetes were the most abundant (∼70%). In the Meth treatment, the most abundant phyla were Proteobacteria and Bacteroidetes (∼55%). In the KanMeth treatment, the most abundant phyla were Proteobacteria and Actinobacteria (∼85%). Finally, the Kan treatment was dominated by the phylum Actinobacteria (∼75%). Similar to phylum level results, family level relative abundances also varied widely among the treatments.

Based on the ANCOM test at the phylum level, the bacterial phylum Elusimicrobia was the most significantly abundant among treatments. In contrast, at the genus level, the following top 13 bacteria were significantly differentially abundant among treatments (centered log ratio, clr>200): uncultured delta proteobacterium (Proteobacteria), uncultured Synergistetes (Synergistetes), Lactococcus (Firmicutes), Propionivibrio (Proteobacteria), uncultured Epsilonproteobacteria bacterium (Epsilonbacteraeota), Tyzzerella 3 (Firmicutes), Treponema (Spirochaetes), Treponema 2 (Spirochaetes), Candidatus Endomicrobium (Elusimicrobia), unidentified Firmicutes (Firmicutes), Stenoxybacter (Proteobacteria), uncultured Bacteroidetes bacterium (Bacteroidetes) and Cerasicoccus (Verrucomicrobia). The observed log abundances tied to different genera are shown in Fig. S3, and the list of significant genera based on ANCOM hypothesis test with their respective clr and W scores is shown in Table S1.

The CCA model indicated that the variation in bacterial species abundance across treatments was significantly correlated with presoldier initiation, total mortality and protist numbers (chi-square, F_3,11=1.7603, P=0.006). The CCA plot shows that the number of protists and presoldiers initiated was strongly associated with the UnCtl and Meth treatments, respectively, whereas total mortality was strongly associated with Kan and KanMeth treatments (Fig. 4). Axes 1, 2 and 3 explain 38.6%, 12.03%, and 7.8% of the total constrained variation, respectively.

The effect of trophallaxis ‘rescue’ treatments on cumulative percentage of presoldiers initiated (H=10.692, d.f.=3, P=0.013) and cumulative percentage mortality (H=16.111, d.f.=5, P=0.006) were significant based on Kruskal–Wallis tests. No presoldiers emerged in the Meth-D+Kan-R and untreated control treatments. All the other treatment combinations that did not result in presoldier emergence were removed from the presoldier initiation data analysis. Based on pairwise Mann–Whitney comparisons (P<0.05), the total percentage of presoldiers emerged in meth controls was significantly higher than the treatment combination Meth-D+Untreated-R termites (Fig. 5). From the same bioassays, survivorship of termites was also recorded. Total mortality in Kan control and Meth-D+Kan-R were significantly higher than in other treatment combinations. The mortalities in Kan-D+Untreated-R and Meth-D+Untreated-R were not statistically different from one another (Fig. S4). Collectively, these results show that presoldier emergence is higher in methoprene-treated individuals held with control individuals containing intact microbiomes than in methoprene-treated individuals held with kanamycin-treated individuals.

This study experimentally quantified the effect of gut symbiont communities on the phenomenon of JH-mediated presoldier differentiation in termites. The study was conducted using methoprene, whose chemical structure closely approximates that of JH (Charles et al., 2011). Methoprene is also consistently available in a highly pure form, while JH III is not, which makes independent verification of findings for JH III difficult. These findings provide a novel example of gut microbial community influence on host phenotype. Four different treatments were used to investigate possible roles of bacterial and protist gut symbionts in the caste differentiation process. These treatments included kanamycin alone (Kan), kanamycin+the JH analog methoprene (KanMeth), methoprene alone (Meth) and an untreated control (UnCtl). The effect of these treatments on presoldier formation, survivorship and protist abundance were observed through phenotypic bioassays and protist quantification. Treatment effects on bacterial abundance and composition were studied through the use of quantitative PCR and16S rRNA gene sequence analysis, which together were used to infer alpha and beta diversities and differential abundance through data analysis procedures.

The JH analog methoprene demonstrates a similar effect as endogenous JH would on a termite worker, resulting in significant presoldier initiation (Zhou et al., 2006b). Kanamycin was previously shown to reduce both protist and bacterial symbiont load in the termite gut (Peterson et al., 2015), and a similar effect was observed in this experiment as well. Therefore, the higher percentage of presoldiers emerged in the Meth compared with the KanMeth treatment suggests that gut symbionts along with JH play a crucial role in presoldier initiation and their survival. In the survivorship assessment among different treatments, the higher percent mortality of termites in the KanMeth treatment could be linked to the double stress from antibiotic and methoprene/JH. That means exposure to kanamycin prior to methoprene reduced the gut microbial load (possibly the gut microbes needed to survive stress during the molting process), resulting in the significantly high mortality observed in phenological bioassays. Additionally, the mortality observed with the Meth-D+Kan-R treatment in our rescue experiment further supports that presoldier initiation within colonies is lethal to termites lacking intact gut microbiomes. The trophallaxis experiment further support that the microbiome is relevant for surviving the presoldier molt in termites.

The mortality observed in the Meth treatment was also higher than the mortality observed with JH III by Hrdý et al. (2006). The possible explanation of this higher mortality could be due to the strong affinity of methoprene at the JH receptor or the stress caused during and after presoldier molting that is normally alleviated by the presence of specific gut symbionts and trophallaxis exchange with other colony members. Another cause of bioassay mortality could be the excessive proportion of presoldiers initiated in our study, which are unable to eat or be fed by workers exposed to methoprene in Petri dishes. Lastly, similar mortality percentages were observed by Howard and Haverty (1978) while using a similar concentration of methoprene as ours. Protist numbers were significantly reduced in termites prior to molting in Meth treatments (Nalepa, 2017; Sen et al., 2013), while some protists were still intact. Termite gut undergoes defaunation during the molting process causing a reduction in the numbers of gut symbionts. This defaunation or reduced protist number could be another cause for the significant mortality we observed in Meth treated individuals. Similar to previous findings (Peterson et al., 2015), protist numbers were reduced in the Kan treatment, which corresponds well with the observed reduction of several bacterial taxa (Spirochaetes, Elusimicrobia, Firmicutes, Bacteroidetes) that are known to maintain symbiotic associations with protists.

Alpha diversity matrices from the 16S rRNA gene sequencing demonstrated that the UnCtl treatment had the highest species evenness and richness, which was expected. Other treatment groups, however, had reduced alpha diversity in ways that were not entirely expected. The reductions could be due to ‘SPRY’ genes, which have been documented in association with reduced gut symbiont numbers in termites when challenged with hormones and antibiotics (Scharf et al., 2017). The significantly higher species evenness and richness in the Meth compared with the Kan and KanMeth treatments strengthens our above conclusion, i.e. that certain gut bacteria are required for successfully completing the presoldier transition, which is well supported by both the presoldier differentiation and trophallaxis experiments. Even though the KanMeth treatment resulted in some presoldier emergence, the total presoldier emergence was reduced significantly in comparison with the Meth treatment. A likely explanation for this difference is that gut symbionts were selected at two levels by kanamycin and JH in the KanMeth treatment. Bacterial diversity was further investigated among treatments using beta diversity matrices. Weighted UniFrac analysis demonstrated that different treatment groups had microbial communities different from each other. Similar results were observed for Bray–Curtis and Jaccard analyses, indicating that each treatment resulted in significantly different bacterial community structures (Lozupone et al., 2007). Overall, from the alpha and beta diversity matrices, and supporting statistical analyses, it is clear that different hormonal and antimicrobial treatments selected for different gut symbiotic load and content. This determination provided rationale to look further into the types of bacteria selected by the different treatments as detailed below.

Several protists, Archaea and Bacteria inhabit the small micro-environment inside a termite gut and provide synergistic physiological benefits to the host (Ohkuma, 2008; Scharf and Peterson, 2021). In this experiment, the major bacterial phyla observed were Spirochaetes, Firmicutes and Bacteroidetes, comprising more than 70% of abundance similar to the previous findings from 16S V5–V6 sampling of the R. flavipes gut lumen (Boucias et al., 2013). Shifts in both bacterial composition and abundance were observed among different treatments, following a hologenome theory of evolution, i.e. the artificially introduced stress was able to change the gut microbial community (Zilber-Rosenberg and Rosenberg, 2008). Most notably, Kan and KanMeth treatments lacked Spirochaetes and Elusimicrobia and had a reduction in Firmicutes. In wood-feeding termites, Spirochaetes are associated with acetogenesis and cellulose-degrading protists, whose absence could be detrimental to the survival of termite. Reduction in xylose released from lignocellulosic diet materials was also observed when the abundances of Elusimicrobia and Firmicutes were reduced (Leadbetter et al., 1999; Lucey and Leadbetter, 2014; Otani et al., 2014; Peterson et al., 2015). This evidence from prior studies highlights the involvement of these bacterial phyla in lignocellulose digestion and ultimately in the survival of termites, which is in agreement with the current findings. Spirochaete bacteria are ectosymbionts to flagellate protists (Iida et al., 2000), whereas Bacteroidetes, Elusimicrobia and Firmicutes are endosymbionts to protists (Fröhlich and König, 1999; Noda et al., 2007, 2005; Ohkuma et al., 2007; Stingl et al., 2005). Because the Kanamycin treatment eliminated or reduced the endosymbionts and ectosymbionts, the disruption of the bacteria–protist symbiosis likely resulted in significant reductions in protist counts and, ultimately, decreased termite survivorship (Peterson et al., 2015; Stephens and Gage, 2020; Treitli et al., 2019).

Based on ANCOM analysis, the genus Cerasicoccus (Verrucomicrobia) differed significantly by being completely absent from the Meth, KanMeth and Kan treatments. The taxa Propionivibrio (Proteobacteria), uncultured Epsilonproteobacteria bacterium (Epsilonbacteraeota), Tyzzerella 3 (Firmicutes), Treponema (Spirochaetes), Treponema 2 (Spirochaetes), Candidatus Endomicrobium (Elusimicrobia) and unidentified Firmicutes (Firmicutes) demonstrated the highest abundance in the UnCtl treatment, while these taxa were completely absent or highly reduced in the Kan/KanMeth treatments and reduced in the Meth treatment, suggesting their importance in the survival and successful completion of the presoldier differentiation process. In contrast, the taxa uncultured delta proteobacterium (Proteobacteria), uncultured Synergistetes (Synergistetes) and Lactococcus (Firmicutes) were equally abundant in the Meth and UnCtl treatments but were completely absent in the Kan and KanMeth treatments. And although the taxa Stenoxybacter (Proteobacteria) and uncultured Bacteroidetes bacterium (Bacteroidetes) were slightly reduced in the UnCtl compared with the Meth treatment, they were completely absent in the Kan and KanMeth treatments, suggesting their importance in the successful completion of the presoldier differentiation process. The significantly high mortality observed in the Kan and KanMeth treatments can be associated with gut stress owing to the disappearance or reduction of these taxa. Another possible explanation based on the trophallaxis experiment could be that when termites were treated with kanamycin and methoprene, they were not able to provide the symbionts to each other needed for survival, resulting in high mortality; this possibility also explains the lower or lack of presoldier differentiation in these treatments and provides further eusocial context to our findings.

The CCA analysis also confirms the strong association of total mortality with the Kan and KanMeth treatments. The reduction of specific protist flagellate-associated bacterial groups (Elusimicrobia, Firmicutes, Spirochaetes) with the Meth treatment is in agreement with the fact that termite guts are voided prior to molting (Nalepa, 2017; Sen et al., 2013) and explains the mortality that might be present in the colony during the differentiation, which could be one of the reasons why there are typically under 3% soldiers in R. flavipes colonies. Based on the presoldier emergence results from the trophallaxis experiment, it is also possible that normally faunated termites have inhibitory effects on presoldier differentiation, which may also contribute to the low numbers of soldiers present in R. flavipes colonies. In humans, bacterial symbionts are documented to impact host lipid metabolism, the immune system and neuroendocrine pathways, which in turn can affect the nervous system and brain function through microbiota–gut–brain axis pathways (Forsythe and Kunze, 2013; Krishnan et al., 2015). Another study in early life stages of mice demonstrated that the acquisition of gut bacteria enhanced stress responses (Sudo et al., 2004). In the present study, the combined results from survivorship, protist numbers, 16S rRNA gene sequence data and the trophallaxis experiment indicated the reduced stress response capabilities of termites resulting from reduced gut symbiont abundance and diversity. In particular, the stress owing to kanamycin and methoprene was reduced by holding treated individuals with faunated termites in the trophallaxis experiment, as would occur under natural colony conditions. Based on the tightly linked associations among different experimental outcomes, it is clear that gut symbionts and their relative proportions play important roles during the soldier caste differentiation process. Prior experimental evidence reveals how JH is able to select for a differential abundance of fat-body-related CYP4 and hexamerin genes and influence their expression during presoldier initiation (Zhou et al., 2006a,b, 2007). Subsequent studies have shown how JH dramatically reduces protist flagellate gene expression and populations specifically in the gut (Scharf et al., 2017; Sen et al., 2013). The current study builds on these prior findings by showing that JH selects for gut symbionts at certain abundance levels, which in turn contributes to successful completion of the presoldier differentiation process.

Lastly, the differential bacterial compositions we identified here between the Kan and KanMeth treatments are one of the advancements provided by this study. Specifically, we found that JH decreases protists and associated bacteria, while the remaining bacteria help emerging presoldiers survive the transition. However, if the gut microbial balance is disrupted before JH application, termites cannot withstand the physiological stress caused by the molting process. These impacts of JH on the gut microbial system are now disentangled by our findings. As the current study demonstrates, gut symbionts under the influence of JH, besides having importance in termite physiological functions such as digestion, reproduction, immunity, nitrogen fixation and acetogenesis (Cleveland, 1923; Rosengaus et al., 2011; Scharf et al., 2011), are directly linked to presoldier differentiation and emergence and are important in the survival of termites under colony conditions. Moreover, the composition and abundance of these gut microbes may predict the phenotype of the termite host. However, the specific mechanisms involved await determination. Unraveling these mechanisms and the roles of specific bacteria is likely to provide important new insights into termite caste differentiation, and more broadly into the processes related to gut microbial control over host phenotypic plasticity.

We thank the University of Minnesota Genomics Center for library preparation and sequencing of the 16S rRNA gene; Hannah Stewart, Ethan Hillman and Brittany Peterson for assay development in earlier phases of this research; Dr Timothy Johnson and Manoj Pandey for assistance with the data analysis; and Zachery Wolfe for laboratory assistance.