Six colonies of the dampwood termite H. sjostedti were collected in Mangshan National Forest Park, Hunan, China. A section of the wood log which was infested by H. sjostedi was collected. The logs were kept in a plastic container (51 × 35 × 30 cm) at 25°C in the laboratory, and water was sprayed onto the log regularly to maintain the moisture of the wood.

Wooden slices (15 × 10 × 0.5 cm) and sticks (15 × 1.5 × 0.5 cm) of pine wood were soaked in water overnight. Two wooden sticks were placed at the bottom of a plastic box (19 × 13 × 7 cm), about 6 cm apart from each other. Then, a wooden slice was placed on the wooden sticks. Afterwards, two wooden sticks were placed on the wooden slice, with the same spacing of 6 cm. Using this method, three pieces of wooden slices were piled together. Several holes were punched on the lid of the plastic box for aeration. A number of 200 workers (6th instar or 7th instar) and five soldiers were selected randomly and kept in the box. All the boxes were maintained at 25°C. The termites were reared for 6 months before sample collection. After 6 months of captivity, termites showed an obvious sign of uric acid accumulation (Xing, et al., 2014; Arango et al., 2017). In summary, there were two groups of termites: the freshly-collected termites (F) and the laboratory-maintained termites (L).

Termites were anesthetized on ice. An individual termite was sterilized in 70% alcohol and then washed in sterilized water. Its gut was pulled out with one tweezer holding the head and another pulling the gut out from the end of the abdomen. Then, the termite was dissected using tweezers. The fat body was harvested using tweezers and stored at −80°C individually. Fat body from the freshly-collected workers (FBF) was collected. Afterwards, the hindgut paunch was broken with tweezers. The gut fluid from a hindgut was sucked with a pipette, and immediately immersed into liquid nitrogen and stored at −80°C individually. Two groups of hindgut fluid were used in this experiment: one group was the hindgut fluid from the freshly-collected workers (HFF) and the other was that from the laboratory-maintained workers (HFL). There were six colonies for each group, with one colony as a replicate (Jiang et al., 2022).

The hindguts of workers were dissected. Three hindguts were pooled together in the experiment, and their DNA was extracted using the Omega E.Z.N.A.^® Stool DNA Kit (D4015-01) according to the manufacturer’s protocol. Two groups of DNA were extracted, namely, the DNA from the hindguts of the F group and the DNA from the L group. There were five colonies for each group, with one colony as a replicate.

The fat body and hindgut fluid were pooled together, respectively, and freeze-dried using a lyophilizer. The weights of the freeze-dried samples were weighed and methanol-acetonitrile-water (2:2:1, v/v/v; internal standard = 20 mg/L) was added in proportion to the weight in order to equalize concentrations of the samples. 2-chloro-L-phenylalanine was used as the internal standard. The solutions were homogenized by a grinder with metal beads, followed by treatment in ultrasound bath for 10 min and incubation for 1 h at −20°C. The samples were then centrifuged at 12,000 rpm for 15 min at 4°C. The resulting supernatants were dried in a centrifugal vacuum concentrator, and then the dried samples were dissolved in acetonitrile-water (1:1, v/v), vortexed for 30 min, and treated with ultrasonic bath for 10 min. Then, the supernatants were collected by centrifugation. The quality control (QC) sample was prepared by pooling 10 μL of the supernatants from all the samples. The QC sample was analyzed to provide robust quality assurance. The supernatants were analyzed with Acquity I-Class PLUS coupled with Xevo G2-XS QTof. The injection volume was 1 μL. Separation was performed on an Acquity UPLC HSS T3 column (1.8 μm, 2.1*100 mm; Waters) using water (0.1% formic acid; mobile phase A)-acetonitrile (0.1% formic acid; mobile phase B) at a flow rate of 400 μL/min. The elution gradient program was as follows: 0 min, 98% A, 2% B; 0.25 min, 98% A, 2% B; 10.0 min, 2% A, 98% B; 13.0 min, 2% A, 98% B; 13.1 min, 98% A, 2% B; 15.0 min, 98% A, 2% B. Under the supervision of the acquisition software, the Waters Xevo G2-XS QTOF high resolution mass spectrometer collected primary and secondary mass spectrometry data in MSe mode (MassLynx V4.2, Waters). Dual-channel data acquisition can be performed on both low collision energy and high collision energy at the same time throughout each data acquisition cycle. For a mass spectrum, the low collision energy was 2 V, the high collision energy range was 10∼40 V, and the scanning frequency was 0.2 s. The ESI ion source’s parameters were as follows: capillary voltage: 2000 V (positive ion mode) or 1500 V (negative ion mode); cone voltage: 30 V; ion source temperature: 150°C; desolvent gas temperature 500°C; backflush gas flow rate: 50 L/h; desolventizing gas flow rate: 800 L/h.

The library construction and sequencing were as follows. To identify the amplicons from different samples, the barcodes were added to the primers 27F and 1492R (27F: 5′-AGRGTTTGATYNTGGCTCAG-3′; 1492R: 5′-TASGGHTACCTTGTTASGACTT-3′). The primers were used to amplify the full length 16S rRNA genes. The thermal cycling conditions were 95°C for 2 min, followed by 25 cycles of 98°C for 10 s, 55°C for 30 s, and 72°C for 1 min 30 s, and then followed by a final extension of 72°C for 2 min. The amplification products were purified, quantified and homogenized to generate the SMRTbell libraries. After library quality control, the qualified libraries were sequenced on PacBio Sequel Ⅱ platform. The bam format files generated by the Pacbio Sequel Ⅱ were converted into CCS (Circular Consensus Sequencing) files by smrtlink. CCS produces a highly accurate consensus sequence from multiple reading of a single SMRTbell molecule (Wenger et al., 2019). The sequences from each sample were identified by the barcode sequences and converted into the fastq sequences.

The original data collected by MassLynx were imported into Progenesis QI for retention time calibration, peak identification, peak extraction and peak integration, and peak alignment, etc. The metabolites were identified based on accurate mass, retention time, and tandem MS data using the METLIN database and a self-established database (Biomarker Technologies Corporation, Beijing). For theoretical fragment identification, the parent ion mass deviation was set within 100 ppm and fragment ions mass deviation was within 50 ppm. The quantification of the metabolites was processed using MetaboAnalystR (Chong and Xia, 2018). The metabolites were then annotated using the Kyoto Encyclopedia of Genes and Genomes database (KEGG) (Kanehisa and Goto, 2000), Human Metabolome Database (HMDB) (Wishart et al., 2018), and Lipid Metabolites and Pathways Strategy (LIPID MAPS) (Fahy et al., 2007). Afterwards, the fold change of each metabolite between two groups was calculated.

To respectively explore the metabolic differences between FBF and HFF, and those between HFF and HFL, the orthogonal partial least squares discriminant analysis (OPLS-DA) was performed (Thévenot et al., 2015). The OPLS-DA is a supervised multivariate data analysis method. The orthogonal variates unrelated to categorical variates were removed from the data set. To validate the robustness of the OPLS-DA model, a random permutation test was performed. To screen the differentially abundant metabolites, a combination of fold change, P-value in the student’s t test, and the variable importance in projection (VIP) value of the OPLS-DA was used. The criteria were FC > 1, P < 0.05, and VIP > 1. The differentially abundant metabolites were summarized in a volcano plot. The differentially abundant metabolites were analyzed with hierarchical clustering analysis. To identify the function of the differentially abundant metabolites, the KEGG pathway enrichment analysis was performed (Yu et al., 2012). Fisher’s exact test was used to perform the KEGG pathway enrichment analysis.

The original bam format data were converted into CCS files by smrtlink through CCS identification, CCS filtration, and Chimera removal as follows. Raw CCS sequences were generated by identification barcode sequences using lima v1.7.0. Then, cutadapt 1.9.1 was used to remove the primer sequences and to conduct length filtration to obtain the clean CCS sequences (Martin, 2011). UCHIME v4.2 was used to remove the chimera to generate the effective CCS sequences (Edgar et al., 2011). The sequences for different samples were identified by the barcode sequences and transformed into fastq format data. The quality of the sequences was evaluated in terms of the number of the sequences in each processing stage.

To explore the taxonomic composition of the hindgut bacteria, the 16S amplicon sequences were first categorized into operational taxonomic unit (OTU) based on sequence similarity. The OTUs were obtained based on 97% similarity with Usearch (Edgar, 2013). One OTU corresponded to one representative sequence. Venn diagram was used to display the relationship among replicates in a group (Chen and Boutros, 2011). To identify the taxonomy of each OTU, the taxonomic annotation and analysis were performed. Naive Bayes classifier was used to annotate the characteristic sequence using the SILVA database (release 132) (Quast et al., 2012). Then, the composition for each sample was calculated at the phylum, class, order, family, genus, and species level. The abundance at each taxonomic level was calculated using QIIME (Bolyen et al., 2019). The bacterial community composition histogram was drawn using R.

Alpha diversity reflects the richness and diversity of each species in a given community. The evaluation indices include Chao1, ACE, Shannon, Simpson, and PD whole tree. Chao1 and ACE were used to evaluate the species richness. Shannon index, Simpson index, and PD whole tree were used to evaluate the species richness and evenness. The alpha diversity indices were calculated using QIIME2 (version 2020.6). The differences of the alpha diversity indices between F and L group were tested by Mann-Whitney U test. To evaluate the similarity of species diversity between the treatments, the beta diversity analysis was performed using Bray Curtis distance metric. Beta diversity is used to compare the similarity of the species diversity between samples. Principal coordinates analysis (PCoA) was used to visualize the beta diversity analysis using R (Dubois et al., 2010).

For the purpose of analyzing the hindgut bacteria difference between the F group and the L group, Metastats analysis (White et al., 2009) and LEfSe (Line Discriminant Analysis (LDA) Effect Size) (Segata et al., 2011) were performed. The FDR adjusted P-value (Q value) was used to determine significance in the Metastats analysis (Q < 0.05). LEfSe is an analysis method to detect and interpret the high-dimensional biomarker, such as genes, pathways, and clades. It is used to detect the features (or biomarkers) which are most likely to interpret the between-group differences, and the extent to which these biomarkers contribute to the difference. The data were subjected to LEfSe analysis to find out the discriminative features in the samples. To study the correlation between different bacteria, Spearman rank correlation was performed. The correlation results whose ρ > 0.1 and P < 0.05 were used to construct the correlation network (Deng et al., 2012). The correlation network at the genus level was drawn using Python.

A total of 2,424 metabolites were identified in the positive mode, and 1,542 metabolites in the negative mode. The metabolites in both modes were combined together, after which the KEGG, HMDB, and LIPID MAPS databases were used to annotate all the metabolites. The quantification and annotation information of the metabolites is summarized in Supplementary Table S1. After checking the list of all the metabolites, the downstream metabolites of the degradation of uric acid were not detected.

FBF and HFF were separated clearly from each other on the OPLS-DA score plot (Figure 1A). Both the R2Y and Q2Y values exceeded 0.9 (R2Y = 1, Q2Y = 0.966), indicating the dependability of the OPLS-DA model. In addition, the result that R2Y value equaled to 1 represented the good accuracy of the OPLS-DA model. Therefore, the differentially abundant metabolites can be screened using VIP. According to the permutation test, the y-intercept of the Q2Y regression line was negative, which indicated that there was no overfitting in the OPLS-DA model (Figure 1B). To summarize, the original OPLS-DA model was quite robust, and was able to explain the metabolic difference between FBF and HFF. The fact that all the red dots were beneath the blue dots suggested that the training dataset and test dataset were relatively independent.

There was also a clear separation between HFF and HFL in the OPLS-DA plot (Figure 1C). The OPLS-DA model had good reliability and accuracy (R2Y = 0.992, Q2Y = 0.717). Therefore, the VIP can be used to screen the differentially abundant metabolites between HFF and HFL. No overfitting of the OPLS-DA model was observed according to the permutation test (Figure 1D). The HFF vs. HFL OPLS-DA also had considerable robustness. In addition, the training dataset was relatively independent of the test dataset, as all the red dots were below the blue dots.

According to the cluster analysis of the differentially abundant metabolites, the six replicates from FBF were more likely to cluster together, while the six replicates from HFF formed their own cluster (Figure 2A). There were 1,443 metabolites with increased abundance, 621 with reduced abundance, and 1,902 with unaltered abundance in FBF, with HFF as a reference (Figure 2B). The differentially abundant metabolites between FBF and HFF are summarized in Supplementary Table S2. The concentrations of the amino acids, including L-valine, L-tyrosine, and L-arginine, were higher in HFF compared with that in FBF. D-fructose 1-phosphate and shikimic acid, the precursors for phenylalanine, tyrosine and tryptophan biosynthesis, were also more concentrated in HFF than in FBF. In addition, the higher concentrations of metabolites in HFF were also observed for the vitamins, cofactors, and their precursors, including vitamin K1, vitamin A, NADP, coenzyme B, (1R,6R)-6-hydroxy-2-succinylcyclohexa-2,4-diene-1-carboxylate, γ-tocotrienol, and 7,8-dihydrobiopterin. Some metabolites associated with fatty acid biosynthesis also had higher concentrations in HFF, including α-linolenic acid, oleic acid, decanoic acid, and stearic acid. In addition, the concentrations of two antibiotics (cephamycin C and geneticin) were higher in HFF compared with that in FBF. Finally, certain metabolites associated with the ovarian steroidogenesis and steroid biosynthesis were more concentrated in HFF than in FBF, including estrone glucuronide, 21-hydroxypregnenolone, pregnenolone, 15 (S)-HpETE, and testosterone.

The cluster analysis of the differentially abundant metabolites revealed that the six replicates of HFF were more likely to cluster together (Figure 2C). By comparing metabolites of HFF and HFL, 691 upregulated metabolites, 267 downregulated metabolites, and 3,008 unchanged metabolites were observed in HFF, with HFL as a reference (Figure 2D). The differentially abundant metabolites between HFF and HFL are summarized in Supplementary Table S3. The concentrations of metabolites, which were higher in HFF than in FBF, increased in HFL after captivity, including but not limited to vitamin A and metabolites in the steroid hormone biosynthesis and ovarian steroidogenesis [estrone glucuronide, 21-hydroxypregnenolone, and 15 (S)-HpETE]. However, for some other metabolites whose concentrations were higher in HFF than in FBF, their concentrations decreased after captivity. These metabolites included the following metabolites: L-tyrosine, cephamycin C, two metabolites in the amino sugar and nucleotide sugar metabolism [UDP-N-acetyl-α-D-glucosamine (UDP-GlcNAc) and guanosine diphosphate mannose], metabolites associated with the ubiquinone and other terpenoid-quinone biosynthesis [γ-tocotrienol and (1R,6R)-6-hydroxy-2-succinylcyclohexa-2,4-diene-1-carboxylate], and metabolites associated with the carbon metabolism (coenzyme B and gluconolactone). Meanwhile, the metabolites associated with the steroid hormone biosynthesis, including cortexolone, 7a-hydroxyandrost-4-ene-3,17-dione, dehydroepiandrosterone, and trans-dehydroandrosterone, increased in concentration in HFL compared with that in HFF. However, no significant difference was observed for the above four metabolites between FBF and HFF. An antibiotic, neamine, also became more concentrated in HFL compared with that in HFF.

The KEGG pathway enrichment analysis was performed on the differentially abundant metabolites between FBF and HFF. According to the analysis, there were two enriched pathways, including amino sugar and nucleotide sugar metabolism (P = 0.025), and nicotinate and nicotinamide metabolism (P = 0.034) (Figure 3A). The differentially abundant metabolites between HFF and HFL were also subjected to KEGG pathway enrichment analysis. There were five pathways enriched according to the analysis, including steroid hormone biosynthesis (P = 0.002), linoleic acid metabolism (P = 0.024), caffeine metabolism (P = 0.037), porphyrin and chlorophyll metabolism (P = 0.037), and PPAR signaling pathway (P = 0.037) (Figure 3B).

A total of 116,869 raw CCS reads was generated through barcode sequences identification (Supplementary Table S4). The average number of the raw CCS reads in the F group (11,689 ± 414) was similar with that in the L group (11,685 ± 847). Then, the primer sequences and chimeric reads were removed. After length filtration, 99.2% ± 0.2% of effective CCS reads was obtained in F group, 98.7% ± 0.4% in L group. The average length of the clean CCS reads was 1,466 bp in F group, and 1,464 bp in L group.

A total number of 166 bacterial OTUs were identified by the 16S amplicon sequencing. The number of OTUs in each sample is shown in in Supplementary Figure S1. The highest was found in a F group (n = 120), while the lowest number of OTUs was found in a L group (n = 79). The F group shared 47 OTUs in common, while 39 OTUs were shared by the L group (Figures 4A, B).

Distribution of the genera for the F group and the L group is summarized in Figure 4C. For clear presentation of the data, only the top 10 genera were shown in the graph, and the rest were combined as others. Treponema was the most abundant genus in both the F group and the L group. The second and third largest genera in both groups were Endomicrobium and Desulfovibrio, respectively.

The alpha diversity indices are summarized in Table 1. No significant differences between F and L group were observed for all the five alpha diversity indices (P > 0.05). The beta diversity of between groups was assessed using PCoA analysis. A clear separation between F and L group was observed. According to the PCoA analysis, PC1, PC2, and PC3 accounted for the 52.40%, 18.23%, and 11.46% of the variation, respectively. The samples from F group clustered more closely than those from L group (Figure 4D).

According to the Metastats analysis, the abundance of five species was significantly different between F group and L group: uncultured Termite Group 1 (TG1) bacterium (Q = 0.03), Candidatus Symbiothrix dinenymphae (Q = 0.03), unclassified Desulfovibrio (Q = 0.03), uncultured actinomycete (Q = 0.03), and Eubacterium nodatum (Q = 0.04) (Figure 5). There were six discriminative features in F group, and four in L group (Figure 6A). The cladogram based on LEfSe analysis is shown in Figure 6B. According to the LEfSe analysis, the species which served a significant role in the L group, including Treponema endosymbiont of Euconympha sp. and uncultured Bacteroidetes bacterium. In the F group, three species were recognized as significant biomarkers, including uncultured TG1 bacterium, Ca. S. dinenymphae, and uncultured Spirochaetes bacterium. According to correlation analysis, the most abundant genus Treponema was negatively correlated with Candidatus Azobacteroides. The second most abundant genus Endomicrobium was negatively correlated with Colidextribacter. None of the two genera showed positive correlation with other genera (Figure 7).