Transcriptomic dissection of termite gut microbiota following entomopathogenic fungal infection

Termites are social insects that live in the soil or in decaying wood, where exposure to pathogens should be common. However, these pathogens rarely cause mortality in established colonies. In addition to social immunity, the gut symbionts of termites are expected to assist in protecting their hosts, though the specific contributions are unclear. In this study, we examined this hypothesis in Odontotermes formosanus, a fungus-growing termite in the family Termitidae, by 1) disrupting its gut microbiota with the antibiotic kanamycin, 2) challenging O. formosanus with the entomopathogenic fungus Metarhizium robertsii, and finally 3) sequencing the resultant gut transcriptomes. As a result, 142531 transcripts and 73608 unigenes were obtained, and unigenes were annotated following NR, NT, KO, Swiss-Prot, PFAM, GO, and KOG databases. Among them, a total of 3,814 differentially expressed genes (DEGs) were identified between M. robertsii infected termites with or without antibiotics treatment. Given the lack of annotated genes in O. formosanus transcriptomes, we examined the expression profiles of the top 20 most significantly differentially expressed genes using qRT-PCR. Several of these genes, including APOA2, Calpain-5, and Hsp70, were downregulated in termites exposed to both antibiotics and pathogen but upregulated in those exposed only to the pathogen, suggesting that gut microbiota might buffer/facilitate their hosts against infection by finetuning physiological and biochemical processes, including innate immunity, protein folding, and ATP synthesis. Overall, our combined results imply that stabilization of gut microbiota can assist termites in maintaining physiological and biochemical homeostasis when foreign pathogenic fungi invade.


Introduction
Microbial symbionts inhabit the guts of many insects, including beetles, silkworms, fruit flies, and termites (Cheng et al., 2017;Garcia-Robles et al., 2020;Li et al., 2021;Vikram et al., 2021). Gut microbiota can be seen as a virtual organ whose existence is indispensable to the life of the host (Gravitz, 2012). Termites possess a wide variety of gut symbionts, including but not limited to bacteria, protists, and fungi, which play many important roles in physiological and biochemical processes such as social alarm, hygienic responses, metabolic capacity, and immunity (Xiong, 2022).
Termites are social insects that live in colonies consisting of three basic caste types: workers, soldiers, and reproductives (He et al., 2018). Because they often live in the soil or in decaying wood, pathogens found in these substrates pose a significant threat to termite colonies (Peterson and Scharf, 2016a). How termites contend with these threats, through mechanisms such as individual and social immunity, has gradually become a new hot spot for research (Avulova and Rosengaus, 2011;Cremer et al., 2018;Hong et al., 2018).
The defensive function of gut microbiota against invading pathogens has been studied in other insects (Zhou et al., 2020;Deng et al., 2022). Serratia and other symbionts render mosquitoes resistant to plasmodium infection through activation of the host immune system or by secreting antimalarial enzymes (Xi et al., 2008;Ramirez et al., 2012;Bai et al., 2019;Gao et al., 2021). In Drosophila, gut microbiota help flies fight pathogens by mediating the renewal of intestinal epithelial cells, and expression of host immune-related genes is altered after the removal of intestinal microorganisms (Buchon et al., 2009;Douglas, 2018).
The lower termites possess particularly robust defenses against fungal pathogens. Individual workers may ingest a large number of fungal spores while grooming nestmates, which become inactive in the gut and do not germinate (Chouvenc et al., 2009). An antimicrobial peptide with antifungal activity, Termicin, was discovered, in addition to synthesis by gut protozoa of the multifunctional β-1,3-glucanases, which break down fungal cell walls (Rosengaus et al., 2014). Similarly, intestinal symbionts protect lower termites from entomopathogenic fungi by producing antifungal compounds and collaborating with host endogenous compounds containing antifungal peptides or enzymes (Peterson and Scharf, 2016b).
In Odontotermes formosanus, a higher termite, only bacterial symbionts are found in the gut (Liu et al., 2013;Hu et al., 2019). Comparably less research has been done on the contributions of bacterial symbionts to termite immunity relative to protists (Peterson and Scharf, 2016a). Previous research on termite immunity has revealed differences in species diversity and relative abundance of the gut microbiota between healthy and pathogen-infected termites (Liu et al., 2015). In addition, the composition and diversity of intestinal symbionts influences gene expression in the termite gut (Xi et al., 2008;Gao et al., 2021). Therefore, it is important to investigate the impacts of bacterial symbionts on termite physiology and immunity.
In this study, we changed the abundance and load of gut microbiota in the hindgut of the higher termite, O. formosanus, which is widely distributed in China, with a broad-spectrum antibiotic, kanamycin, a bacterial antibiotic which binds to the bacterial ribosome 30S subunit to inhibit protein synthesis (Shu et al., 2011). In our previous 16s rDNA sequencing experiment, we found kanamycin significantly impacted intestinal composition, especially Firmicutes and Spirochete. Then, we infected termites with a pathogenic fungus, Metarhizium robertsii, and measured the expression changes of genes related to immunity and energy metabolism in intestinal tissue. Analyzing the effects of specific gene expression changes following immune challenge contributes to further understanding of the symbiotic relationship between gut microbiota and their termite hosts and can also provide reference and inspiration for related research in other insects.

Experimental insects and setup
O. formosanus colonies were collected from Lion Mountain (Huazhong Agricultural University, Wuhan, Hubei, China) in July 2020. Termites were transferred from pieces of rotting wood to 9 cm diameter Petri dishes with filter paper and placed in darkness for 24 h to stabilize the community. Afterwards, termites were maintained in 9 cm diameter Petri dish in 24-h darkness at 25°C ± 1°C and were fed with filter paper soaked with sterile water that was changed once every 2 days. M. robertsii ARSEF#2575 was obtained from the Institute of Plant Physiology and Ecology, Chinese Academy of Sciences, Shanghai. M. robertsii spores were grown on Potato Dextrose Agar medium (PDA) (200 g potato, 20 g glucose, 15-20 g agar, 1000 mL distilled water) and washed with 0.1% Tween 80 solution. Spore concentration was measured with a hemocytometer, then spores were stored at 4°C as a spore suspension. The concentration of antibiotics used was determined by previous research (Peterson et al., 2015).
Before fungal infestation, we separated termite colonies into four experimental groups (control groups/CT: treated with sterilized water; Kana groups/AT: treated with 5% kanamycin; Kana_ infected groups/IAA: treated with 10 8 spores/mL suspension at 48 h time point after treatment with 5% kanamycin; Infected groups/MI: treated with 10 8 spores/mL suspension at 48 h time point), with every group consisting of three biological replications containing 30 termites each. 5% kanamycin means that the ratio of powder to distilled water is 1:20. Termites were reared on a filter paper soaked with 1 mL of the corresponding solution for their treatment group. Before dissecting, all termite samples were placed on ice and 75% ethanol was used to remove any surface microorganisms. We use anatomical forceps to tear open the abdominal epidermis of termites under a posture microscope, pull out the intestine, and ensure its integrity then transferred to an empty Eppendorf tube. The entire process is immersed in 1× sterile PBS (Phosphate Buffer Solution). All gut samples were stored at −80°C until RNA extraction.

RNA extraction
All samples were homogenized in Trizol (Ambion/Invitrogen, United States) following the manufacturer's protocol, purified by isopropanol precipitation, and dissolved in DEPC water. Each treatment contains 30 termite samples and sets up 3 biological replicates. RNA concentration was measured using 1% agarose gel electrophoresis and a Nano-Drop 2000 spectrophotometer. RNA Frontiers in Physiology frontiersin.org 02 completion was measured using an Agilent 2,100 bioanalyzer. Then, RNA products were stored in a −80°C cryogenic refrigerator.

Library preparation and sequencing
Library preparation and sequencing were performed by Novogene (Tianjin, China) on an Illumina NovaSeq 6,000 (Illumina Inc., San Diego, United States) platform. mRNA was isolated using Oligo (dT) magnetic beads (Thermo Fisher Scientific, Waltham, Mass., United States). Subsequently, the obtained mRNA was fragmented using NEB Fragmentation Buffer (New England Biolabs (Beijing) LTD., Beijing, China). Using fragmented mRNA as a template and random oligonucleotides as primers, the first strand of cDNA was synthesized using the M-MuLV reverse transcriptase system, and the RNA strand was degraded using RNase H (TAKARA Bio, Beijing, China). In the DNA polymerase I system, dNTPs (TAKARA Bio, Beijing, China) were used as primers to synthesize the second strand of cDNA. The purified doublestranded cDNA underwent end-repair, A-tailing, and ligation of sequencing adapters. AMPure XP beads (Beckman Coulter, Shanghai, China) were used to screen 250-300bp cDNAs, then PCR amplification was performed and AMPure XP beads were used again to purify the PCR products, which constituted the complete library. After the library was constructed, preliminary quantification was carried out using a Qubit2.0 Fluorometer (Thermo Fisher Scientific, Waltham, Mass., United States), then, library concentration was diluted to 1.5 ng/µL. Next, insert size was assessed with an Agilent 2,100 bioanalyzer (Agilent Technologies, Beijing, China). The raw data were submitted to NCBI (National Center for Biotechnology Information) databases (Accession numbers: PRJNA926100).

Differential gene expression analysis
The transcripts were assembled using Trinity as reference sequence and clean data were mapped to the reference sequence using RSEM and bowtie2 with default parameters (Li and Dewey, 2011). Then, read counts were obtained by calculating the bowtie2 result and the FPKM value of each gene was derived. Differences in gene expression were analyzed using the DESeq2 R package (1.6.3) with default parameters. Finally, | log2FoldChange| > 1 and P-adj < 0.05 were set as standards to identify differentially expressed genes (DEGs) (Wang et al., 2010;Love et al., 2014).

KEGG and GO enrichment of DEGs
Enriched GO terms among differentially expressed genes were identified using GOseq (Young et al., 2010) and KOBAS3.0 was used to analyze the statistical enrichment of DEGs among KEGG pathways (Mao et al., 2005).

Quantitative real-time PCR (qRT-PCR)
A total of 1.0 μg of RNA was reverse-transcribed in vitro using a PrimeScript ™ RT reagent kit, following the manufacturer's instructions. The cDNA products were used as a template for real-time qPCR. RT-qPCR was used to detect gene expression levels. The primers used in this study were designed by NCBI Primer-BLAST software (https://www.ncbi.nlm.nih.gov/tools/ primer-blast/) and are shown in Table 1. The reaction mixtures were prepared using a NovoStart ® SYBR qPCR SuperMix Plus kit (Novoprotein Technology Ltd., China), following the manufacturer's instructions. A 10 μL qPCR reaction system was used, including 5 μL of 2 × NovoStart ® SYBR qPCR SuperMix Plus, 0.5 μL of upstream and downstream primers, 1.5 μL of the template, and 2.5 μL of ddH2O. The reactions were performed on a LightCycler ® 96 System (Roche, Basel, Switzerland). The following qPCR protocol was used: 1 cycle at 95°C for 5 min, followed by 40 cycles at 95°C for 20 s, and 54°C for 40 s. The 2 −ΔΔCT method was adopted to calculate relative gene expression levels. Each group was repeated three times. β-actin was used as the internal reference gene (Supplementary Table S4).

Kanamycin accelerates termite death by compromising its gut microbiota
To verify which antibiotics can play a significant role in removing termite intestinal microorganisms, we fed termites with four different common antibiotics and measured the change in total bacteria abundance at multiple time points (0, 24, 48, and 72 h) using qRT-PCR (primers shown in Table 1, termite β-actin as reference gene). The total bacterial changes in the guts of termites exposed to four different antibiotic treatments are shown in Figure 1. The largest change was observed in the kanamycintreated individuals, particularly at the 24 h time point (p < 0.0001) ( Figure 1A a). Based on the results of this assay, kanamycin was ultimately selected as the antibiotic to be used in the following experiments.
Frontiers in Physiology frontiersin.org Moreover, there was no significant difference in the survival rate of kanamycin treated termites at 48 h compared to controls ( Figure 1B). At 72 h, mortality was significantly higher in the kanamycin-treated group (p < 0.0001). This result suggests that termites remain in a healthy state of life for at least 48 h after intestinal microorganism inhibition. Therefore, we chose to treat termites with antibiotics for 48 h before M.robertsii infection in subsequent experiments.

Transcriptome sequencing, assembly, and annotation
In this study, a total of 108.6 Gb of raw data were obtained from Illumina sequencing. After quality control, a total of 105.8 Gb of clean data was used for further analysis (Table 1). There were 142531 transcripts and 73608 Unigenes obtained by Trinity; their N50 length was 3,655 and 2746bp, respectively (Supplementary Tables S1,S2).
Annotated genes were placed into one of three groups based on GO classification: BP (biological process), CC (cellular component), or MF (molecular function) (Supplementary Figure S1). Among BP genes, cellular process and metabolic process were the most common classifications. Among CC genes, cellular anatomical entity, intracellular, and protein-containing complex were the most common classifications. Among MF genes, binding and catalytic activity were the most common classifications.

Identification of DEGs following antibiotic treatment and/or pathogen exposure
To identify the immune effect of the inhibition of gut microorganisms on termites exposed to entomopathogenic fungi, we used DESeq2 software to analyze differences in gut gene expression between termites from four different treatment groups. A corrected p-value of 0.05 and log2fold-change of ±1 was set as thresholds for significant differential expression. Overall, 43 unigenes were upregulated and 57 were downregulated in MI compared to CT termites; 2,204 were upregulated, and 1610 were downregulated in IAA compared to CT termites; 2,711 were upregulated and 2,244 were downregulated in IAA compared to MI termites; 1160 were upregulated and 1133 were downregulated in IAA compared to AT termites; and 125 were upregulated and 204 were downregulated in AT compared to MI termites (Supplementary Figure S3). When termites are infected by M. robertsii, their gut microbiota can adapt to the external environment to a certain extent. Similarly, following antibiotics treatment, the impact on uninfected termites is smaller. A total

FIGURE 1
The Effect of Antibiotics on Termites (A) The effects of antibiotic treatment on bacterial symbiont abundance in the termite gut. The qPCR results for total bacteria in the guts of termites treated with 5% Kanamycin (a), 5% Ampicillin (b), 5% Streptomycin (c), and 2.5% Tetracycline (d). β-actin was used to normalize the data, which are shown as mean ± standard error; the mean was obtained from three independent replications. The 2 −ΔΔCT method was used to calculate the relative expression level. Statistical differences between time points, determined with an unpaired t-test, are shown via GraphPad Prism8. Different letters indicate significant differences (B means p < 0.0001, C means p < 0.05). (B) Number of dead termite individuals fed with 5% kanamycin.AT termites were fed with 5% kanamycin, while CT termites were fed with sterile water. The experiments were performed in three biological replicates. Kanamycin significantly increased mortality compared with the control group after 72 h (p < 0.05). The experiments were performed in three biological replicates (n = 30 per replicate). Statistical differences between treatment at every time points, determined with an paired t-test, are shown via GraphPad Prism8.
Frontiers in Physiology frontiersin.org of 22 DEGs were identified between the MI: CT group, 115 between the AT: CT group, 555 between the IAA:AT group, and 3,070 between the IAA:MI group (Figure 3).

Enrichment analysis of DEGs reveals the existence of gut symbionts that maintain termite intestinal stability during fungal invasion
After identifying the numbers of DEGs among different pairs of the four treatment groups, we also examined the placement of these DEGs among KEGG pathways (Figure 4).
When we examined DEGs identified in the AT:CT group, it was found that most metabolic pathways did not change in kanamycin treating 48 h, indicating that antibiotic treatment had no negative effect on gene expression or the physiological status of the termite intestinal tract itself. Examination of DEGs identified in the MI:CT group also yielded a low number of metabolic pathways showing change despite infection by M. robertsii. DEGs identified in the IAA:AT and IAA:MI  Frontiers in Physiology frontiersin.org groups suggest significant changes in several metabolic pathways, including the cAMP and calcium signaling pathways. In addition, MAPK signaling pathway were enriched in IAA: CT groups by KEGG enrichment (Supplementary Figure S4A). These results demonstrate that the presence of intestinal microorganisms can alleviate the metabolic effect of fungal invasion to a certain extent. The GO classification of DEGs identified in each comparison group was also examined ( Figure 5). DEGs identified in the AT:CT and MI:CT comparisons showed minimal changes within the CC and MF categories. DEGs placed in the BP category mostly involved membrane-related reactions. DEGs identified in the IAA:AT and IAA:MI comparison groups were spread between all three GO categories.
The results of KEGG and GO enrichment indicate that a change in intestinal symbionts alters host gene expression in the gut and causes the transformation of many metabolic pathways. However, unexpectedly, the main enriched pathways were not directly related to immunity, but indirect pathways such as signal transduction, indicating that gut microorganisms may affect host immunity through other forms of regulation.

qRT-qPCR analysis
To further evaluate the DEGs identified in the transcriptome and better explain the role of intestinal microorganisms in termite  Supplementary Table S4 and Blast results are shown in  Supplementary Table S5. Gene expression patterns were divided into five groups. Eight genes showed lower expression in the AT and IAA treatments than in the MI treatment, while expression in the CT and MI treatments was relatively equal ( Figure 6A). These eight genes are mainly involved in the process of reverse transcription and cell division. The second group of genes, four in total, showed upregulation in the MI treatment, downregulation in the IAA treatment, and no change in the AT treatment ( Figure 6B). These genes included ApoA2, Cal-5like, and Hap. The third group of genes, three in total, were downregulated in all non-CT treatments and included genes involved in sodium and potassium ion transport ( Figure 6C). The change in the internal environment of termites treated with antibiotics and/or pathogens or the influence of external stimuli may inhibit this process. Next, the gene chitinase5 mainly functions as a chitinase and participates in the decomposition of β-1,4-glucan ( Figure 6D), which is the main component of plant cellulose, indicating that the antibiotic treatment does not affect the ability of termites to degrade cellulose at this time point. In addition to the genes that cannot be annotated, in these three genes, hsp70 showed that all the treatment groups were upregulated. But ATP-syn expressed in MI was higher than CT or AT and low expressed in the IAA group.

Discussion
Various symbionts inhabit the intestinal tracts of termites, and these symbionts play an indispensable role in the lives of their hosts, aiding in nutrient supply, digestion, nitrogen metabolism, and more (Liu et al., 2013;Peterson and Scharf, 2016a). However, the role of intestinal microorganisms in the immunity of their termite hosts, such as when encountering a fungal pathogen, is unclear (Peterson and Scharf, 2016b). We exposed O. formosanus termites to antibiotics, pathogenic M. robertsii, or both in sequence, and then measured gut gene expression. Termites exposed to antibiotics followed by a fungal pathogen showed significantly larger changes in gene expression than those exposed only to one or the other (Figure 3), suggesting that gut symbionts can buffer their hosts against foreign pathogens.  (D) no changed in all groups and (E) the remained genes. β-actin was used to normalize the data, which are shown as mean ± standard error; the mean was obtained from three independent repeats. The 2 −ΔΔCT method was adopted to calculate the relative expression level. The statistical difference relative to the control (CT) group, determined with an unpaired t-test, is shown via GraphPad Prism8. Asterisks indicates statistically significant differences (* = p < 0.05; ** = p < 0.01; *** = p < 0.001). (CT: treated with sterilized water as control treatment; AT: treated with 5% kanamycin antibiotic; MI: treated with 108/mL M. robertsii conidia suspension at 48-h; IAA: treated with 5% kanamycin, then treated with 108/mL M. robertsii conidia suspension at 48 h).

Frontiers in Physiology frontiersin.org
This relationship between gut microbiota and immunity has been demonstrated previously in the lower termite Reticulitermes flavipes, which shows greatly increased sensitivity to the fungal pathogen Beauveria bassiana following antibiotic treatment (Peterson and Scharf, 2016b). A similar relationship is observed in mosquitoes, which possess midgut microorganisms that mediate host resistance to foreign pathogens, including plasmodium (Dong et al., 2009;Aida et al., 2013;Gao et al., 2021), viruses (Xi et al., 2008), and pathogenic fungi (Wei et al., 2017). Most of these symbionts directly affect the gene expression of the Toll or Imd immune pathway of the host.
The results of our GO and KEGG enrichment analyses indicated that identified DEGs were not involved in direct immune pathways, but rather in signaling pathways, including the calcium signaling pathway (Figures 4, 5). Antibiotics alone are not considered immune activation conditions for termites, so simple antibiotic treatment will not cause immune response in termites. Calcium (Ca 2+ ) is a second messenger that participates in the modulation of various biological processes, including cell survival, cell proliferation, apoptosis, and the immune response (Kong et al., 2021). There is increasing evidence that Ca 2+ signaling is beneficial for JAK-STAT activation induced by different stimuli and that it plays a vital role in regulating the activation of NF-κB to facilitate gene transcription (Wang et al., 2008;Liu et al., 2016). In addition, we found that the expression of genes related to melanogenesis, which is involved in the innate immunity of insects (Koike and Yamasaki, 2020), changed dramatically in the IAA treatment group (Figure 4). Our GO enrichment analysis also showed a lack of DEGs involved in direct immune pathways ( Figure 5). Rather than directly impacting immune pathway gene expression, gut bacteria may contribute to host immunity through other means, such as by affecting host signal transduction and melanogenesis or by changing intestinal oxygen levels at the time of pathogen infestation (Nappi and Vass, 1993;Nappi et al., 2009).
Using qRT-PCR analysis, we measured the gene expression changes of the top 20 most significant DEGs identified from our transcriptomic study ( Figure 6). We found that antibiotic treatment did not affect the expression of chitinase in the host ( Figure 6D). At the same time, expression of the reverse transcription-related genes RT-like 1, RT-like 2 and RT-like 3 in the intestinal tract were decreased after antibiotic treatment ( Figures 6A,E). Studies have shown that intestinal microbial imbalance can cause cell cancerization and changes in DNA levels, especially in adverse environments (Sobhani et al., 2019). It can be speculated that the presence of pathogenic bacteria does not affect the reverse transcription process of termites, but that the imbalance of intestinal microorganisms leads to this change. Expression of the apolipoprotein gene APOA2 was upregulated in MI termites but downregulated in IAA termites ( Figure 6B). Apolipoproteins participate in the innate immune regulation of insects, demonstrating antibacterial activity and cooperation with antimicrobial peptides (Rahman et al., 2006;Hanada et al., 2011;Feingold and Grunfeld, 2012;Kamareddine et al., 2016). The pattern of expression that was observed in MI and IAA termites suggests that the balance of intestinal microorganisms is necessary for termites to resist pathogen invasion. Interestingly, genes associated with calcium signal pathway were highly enriched in our KEGG analysis, and calpain ( Figure 6B) is related to apoptosis and cell division (Shi et al., 2022). Finally, haptoglobin ( Figure 6B) contains trypsin-like serine protease, which is also related to immunity (Pászti-Gere et al., 2021;Yang et al., 2021). Combined with the above results, we speculate that antibiotic treatment slows down the apoptosis of host cells and increases the expression of apolipoprotein on the surface of cell membranes.
We also measured HSP70 and ATP-syn gene expression in all four treatment groups ( Figure 6E). These genes are involved in physiological stress processes, aiding primarily in the folding and stretching of newly synthesized proteins and in the repair of misfolded proteins (Fernández-Fernández and Valpuesta, 2018). Under the stimulation of fungal infection, the expressions of these genes increased. It is speculated that when insects encounter pathogens, HSP70 is upregulated to protect against cellular injury (Tang et al., 2012). Expression of ATP-syn was increased in MI termites, while decreased in IAA termites. We speculate that M. robertsii infection led to an increase of protein synthesis in the host, resulting in a large amount of ATP production and consumption. Downregulation of ATP-syn in the IAA group suggests that intestinal symbionts play a role in host ATP synthesis during pathogen invasion (Peterson and Scharf, 2016a).

Conclusion
In this study, termites exposed to antibiotics followed by a fungal pathogen showed significant changes in gene expression in comparison to those challenged exclusively by one factor. Notably, these changes affected genes associated with calcium signaling pathway and melanization process. These results suggest that intestinal symbionts may play a buffering role when the insect host encounters foreign pathogens, stimulating the expression of immune related genes, such as apolipoprotein. Symbionts may also affect the protein synthesis process and energy supply of the host through HSP70 and ATP synthesis. The gut microbiota provides essential health benefits to its host, particularly by maintaining homeostasis. Disruption of such homeostatic balance can lead to significant adverse impacts, including diseases and lethality. Given that gut microbiota and host immune system have coevolved to maintain homeostasis, genetic underpinnings govern the dialogues between the two is one of the emerging questions warrant immediate attention.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, PRJNA926100.