Apriona swainsoni is a pest widely distributed in southeast China, which is harmful to the economical and natural trees (Zhang S. K et al., 2022). Totally 150 fourth or fifth-instar larvae of A. swainsoni, which are in the middle stage of larva development with typical physiology phenotype and metabolism status, were collected in the C. decapetala wood branches. The insect species were recognized by the specified head and pronotum of the larva of A. swainsoni, and four gut segments of all the collected larvae were dissected for further examination. Each gut segment for multi-omics analyses contains four samples of biological repeats. Every gut segment sample consists of 8–12 dissected gut segment organs to meet the quantity requirement of multi-omics analyses. Then the C. decapetala branches were used for validating C/N ratio of larvae food. The gut segments with lumen contents and plant branches were stored at −80°C for subsequent analysis.

QIAamp^® Fast DNA Stool Mini Kit (Qiagen, Hilden, Germany) was used to isolate total DNA of samples following the manufacturer’s instructions. A NanoDrop2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States) and agarose gel electrophoresis were used to assess DNA concentration and integrity. The concentration of extracted DNA was 20–68 ng/μl. Agencourt AMPure XP beads (Beckman Coulter Co., United States) were used to clean up the DNA fragmented with S220 Focused-ultrasonicators (Covaris, United States). Afterward, libraries were constructed using the TruSeq Nano DNA LT Sample Preparation Kit (Illumina, San Diego, CA, United States). Metagenome sequencing and analysis were performed by OE Biotech Co., Ltd. (Shanghai, China).

Trimmomatic was used to trim and filter FastQ sequences (Bolger et al., 2014). Host pollution control was performed due to the DNA was extracted from a host-related sample. To align the post-filtered pair-end reads with the host genome, BWA was used (Li and Durbin, 2010) and the aligned reads were discarded. After valid reads were obtained, the genome was assembled using SOAPdenovo2 (Luo et al., 2012). Use gaps inside the scaffold as breakpoints to interrupt the scaffold into new contigs (Scaftig), and these new Scaftigs with lengths > 200 bp (or 500 bp) were retained. ORF prediction of assembled scaffolds using prodigal was performed and translated into amino acid sequences. The non-redundant gene sets were built for all predicted genes using CDHIT (Fu et al., 2012). The clustering parameters were 95% identity and 90% coverage. The longest gene was selected as the representative sequence of each gene set. Clean reads of each sample were aligned against the non-redundant gene set (95% identity) using Bowtie2 (Langmead and Salzberg, 2012), and the abundant information of the gene in the corresponding sample was counted. Then, the gene set representative sequence (amino acid sequence) was annotated with NR, KEGG, COG, SWISSPROT, and GO database with an e-value of 1e-5. Taxonomy for the species was obtained from the NR Library’s corresponding taxonomy database, and the abundance of the species was calculated using the corresponding abundance of the genes. At each taxonomy level of Domain, Kingdom, Phylum, Class, Order, Family, Genus, and Species., abundance statistics were performed to construct the abundance profile.

The gene sets were compared with the CAZy database (Drula et al., 2022) using the corresponding tool hmmscan to obtain information on the carbohydrate-active enzyme corresponding to the gene. R software (v3.2.0) was used for the PCA analysis and plotting of species or functional abundance spectra.

For bin analysis, Megahit (Li et al., 2015) was used to assemble clean sequencing data and remove hosts to obtain contigs. Bowtie2 (Langmead and Salzberg, 2012) and Samtools (Li et al., 2009) were respectively used for comparison and format conversion. After obtaining the depth data of Contigs, Contigs above 1,500 bp were used for bin classification and Metawrap (Uritskiy et al., 2018) software was used. Three built-in software algorithms, Maxbin2 (Wu et al., 2016), MetaBAT2 (Kang et al., 2019), and Concoct (Lu et al., 2017), are used to perform macro factor component bins, and the best results are obtained from highly similar bins, and non-redundant bins are output. The Lineage_WF function in ChecKM software (Parks et al., 2015) was applied to evaluate the quality of the box sorting results, and to obtain the completion degree, pollution degree, species classification, and other information of each box. Bins with a screening completion degree greater than 50% and a contamination degree less than 10% were used for subsequent analysis. Based on NCBI_tax and NCBI_nt databases, Taxator-TK (Droge et al., 2015) was used to annotate each contig species and then estimate the overall bin species.

The contents of the intestinal cavity of larvae from the four main regions of the digestive tract were homogenized and centrifuged at 21,000 g for 10 min at 4°C. Then the supernatant was transferred to methanol/chloroform aqueous solution, and was remove the upper water-soluble phase which was put into a labeled glass bottle with completely dry for metabolomics analysis. Then, the settled protein was added to 1 ml of ice-cold methanol. The sample was centrifuged again at 10,000 g for 10 min to precipitate the protein group, the methanol was removed, and the sample was placed in the fume hood to dry. Sample lysate was added to the precipitation, and the protease inhibitor PMSF was added to make the final concentration of 1 mM. The solution was centrifuged twice at 12,000 g at 4°C for 10 min to obtain the supernatant. The supernatant is the total protein solution, which shall be determined for protein concentration and stored at −80°C for later use. The concentration of extract protein was 1.2–8.2 μg/μl.

A Q-Exactive HF mass spectrometer (Thermo, United States) equipped with a Nanospray Flex source (Thermo, United States) was used to perform the metaproteome identification. Samples were loaded and separated by a C18 column (15 cm × 75 µm) on an EASY-nLCTM 1200 system (Thermo, United States). The flow rate was 300 nl/min and the linear gradient was 75 min (0–63 min, 5%–45% B; 63–65 min, 45%–90% B; 65–75 min, 90% B; mobile phase A = .1% FA in water and B = .1% FA in ACN). The AGC target value was set at 3e6, and full MS scans were acquired between 350 and 1,500 m/z with a mass resolution of 60,000. With a collision energy of 32, the 20 most intense peaks in MS were fragmented by higher-energy collisional dissociation (HCD). MS/MS spectra were obtained with a resolution of 45,000 with an AGC target of 2e5 and a max injection time of 80 ms. The Q Exactive HF dynamic exclusion was set for 30.0 s and run under positive mode. The raw data was thoroughly searched against the sample protein database using ProteomeDiscoverer (v.2.4). The database was searched with the specificity of Trypsin digestion. In the database search, alkylation on cysteine was considered a fixed modification. For the protein quantification method, the TMT labeling method was selected. The false discovery rate (FDR) was set at .01 globally, and at least two peptides were required for each protein group to be quantified.

An accurately weighed sample of 60 mg was transferred to a 1.5 ml Eppendorf tube. Two small steel balls were added to the tube. For each sample, an internal standard (2-chloro-l-phenylalanine in methanol, .3 mg/ml) and an extraction solvent (4/1, v/v) were added. Samples were stored at −80°C for 2 min and then grinded at 60 Hz for 2 min, ultrasonicated at ambient temperature (25°C–28°C) for 10 min, and stored at −20°C for 30 min. The extract was centrifuged at 15,000 g, 4°C for 15 min. With the aid of a freeze concentration centrifugal dryer, 300 μl of supernatant was dried in a brown and glass vial. 400 μl mixture of methanol and water (1/4, vol/vol) was added to each sample. Then the samples were vortexed for 30 s, then placed at 4°C for 2 min. Samples were centrifuged at 15,000 g, 4°C for 5 min. The supernatants from each tube were collected using crystal syringes, filtered through .22 μm microfilters, and transferred to LC vials. The vials were stored at −80°C until LC-MS analysis. QC samples were prepared by mixing aliquots of all samples to be a pooled sample. The analytical instrument for this experiment is a liquid mass coupled system composed of Dionex U3000 UHPLC ultra-performance liquid tandem QE Plus high-resolution mass spectrometer.

Progenesis QI software (Waters Corporation, Milford, United States) was used to analyze the acquired LC-MS raw data with the following parameters. The precursor tolerance was set to 5 ppm, fragment tolerance was set to 10 ppm, and retention time (RT) tolerance was set to .02 min. Internal standard detection parameters were deselected for peak RT alignment, isotopic peaks were excluded for analysis, the noise elimination level was set at 10.00, and the minimum intensity was set to 15% of base peak intensity. The Excel file was obtained with three-dimension data sets including m/z, peak RT and peak intensities, and RT–m/z pairs were used as the identifier for each ion. The resulting matrix was further reduced by removing any peaks with a missing value (ion intensity = 0) in more than 50% of samples. Data QC (reproducibility) was performed by the internal standard. Based on public databases such as http://www.hmdb.ca/; http://www.lipidmaps.org/ and self-built databases progenesis QI software was used to identify the metabolites. The differential metabolites among groups were determined based on the combination of a statistically significant threshold of p values from a two-tailed Student’s t-test on the normalized peak areas and variable influence on projection (VIP) values, where metabolites with and p values less than .05 and VIP values larger than 1.0 were considered as differential metabolites.

Frass samples obtained from five larvae and wood samples were loaded in 8 mm × 5 mm tin capsules using 10–20 mg for nitrogen quantification and 2 mg for carbon quantification. A zero blank autosampler connected to an ECS 4010 Element Analyzer (Costech Analytical) coupled to a Delta V Plus isotope ratio mass spectrometer (Thermo Fisher Scientific) was used to identify the carbon and nitrogen concentration of the samples. The analytical precisions on concentration yields were 70.14 ± 1.29 wt% (1 s; n = 7) for carbon and 4.83 ± .09 wt% (1 s; n = 10) for nitrogen based on the laboratory atropine standard sample.

The four gut region samples of larvae were isolated and dissected in a sterile environment, and the guts were immediately flushed with 4% paraformaldehyde. After flushing, the samples were fixed in 5 times the volume of 4% paraformaldehyde buffer for at least 72 h. After dehydration and paraffin embedding, 5 mm gut tissue sections were taken for Hematoxylin-eosin (H & E) staining. The morphology of gut tissue was observed under an upright microscope (OLYMPUS BX51, Tokyo, Japan).

The gut tissues from each gut region were isolated in a sterile environment. Total RNA was extracted from gut tissue using a Trizol reagent (Ambion, United States). The concentration of extracted total RNA was 35–80 ng/μl. FastQuant RT Kit reverse transcription kit (TIANGEN, Beijing, China) was applied to reverse RNA into complementary DNA (cDNA). NifH and ureC genes were amplified by RT-PCR using primers as follows: nifH forward primer: AAAGGYGGWATCGGYAARTCCACCAC, nifH Reverse primer: TTGTTSGCSGCRTACATSGCCATCAT; ureC forward primer: ATHGGYAARGCNGGNAAYCC, ureC Reverse primer: GTBSHNCCCCARTCYTCRTG. The primers were designed by Primer Premier 5.0 and synthesized by Sangon Biotech (Shanghai, China).

The data were processed and analyzed by SPSS 23.0 and GraphPad Prism 6.02 statistical software. Based on the normality test and homogeneity of variance, the Student-t-test was used when two groups of data were compared (Welch’s t-test was used when the variance was not homogeneous). One-way ANOVA analysis of variance with the Tukey-Kramer method was used when multiple groups were compared.

Through anatomical analysis, we found that the larva of A. swainsoni is about 3 cm long, and its intestinal length is about 11–13 cm long (Figure 1A, B). According to morphological characteristics, the gut can be divided into four segments: foregut (FG), midgut (MG), anterior hindgut (AHG), and posterior hindgut (PHG) (Figures 1C, D). According to HE staining results, the diameter of MG is small, and the intestinal wall of MG is dense and thick. The diameter of PHG is large, and the intestinal wall of PHG is loose and thin. These features of FG and AHG are fallen in between MG and PHG. Thus, the presented results suggested that the morphology of four gut segments is also different, and may function diversely in A. swainsoni (Figure 1E).

To understand the effects of different gut segments on larval habitat adaptation, we sequenced and identified the metagenome, metaproteome, and metabolome of four gut segments of the A. swainsoni larva. The results of metagenomic principal component analysis (PCA) at the species level suggested that there were no significant differences between samples from different gut segments (Figure 2A), but at the gene number level, there were significant differences between gut segments (p < .05), except for foregut and midgut, and foregut and hindgut pairs (Figure 2B). According to taxonomy barplot comparisons, samples from different segments of the larvae differ at both phylum and species levels (Figures 2C, D). Top 15 phylum comparison results showed that there are differences in four gut segments since the eighth rich abundance phylum, i.e., Basidiomycota, Actinobacteria, Microsporidia, Planctomycetes, Cossaviricota, Candidatus, and Uroviricota. Top 15 species comparison results showed that there are differences in four gut segments since the fifth rich abundance species especially Gibbsiella quercinecans, Enterobacter cloacae, Flavobacterium sp., and Enterobacter sp. Notably, the results of PCA showed that the metaproteome and metabolome of four segments of A. swainsoni larvae were aggregated together and significantly different from other segments, suggesting that the proteome and metabolome of four segments of A. swainsoni larvae were significantly different, which may be related to the functional differences of different segments.

To determine crucial microbial organisms in nitrogen fixation and lignocellulose deconstruction across the four gut segments, we performed genome bins by metagenomic data. A total of 22 genome bins were obtained, with 10 bins probably composed of a single bacterial strain with genome contamination ratio lower than 10% and genome completeness larger than 50% (Supplementary Table S1). The taxonomy and the microbial composition of the whole metagenome of these 10 bins were determined through analysis of their bacterial markers (Supplementary Table S1). Notably, all 10 genomes are determined as bacteria. These 10 genome bins were further annotated to verify the existence of functional genes coding for lignocellulose deconstruction enzymes, and nitrogen fixation-related enzymes (Supplementary Table S2). These 10 genome bins were also annotated by the CAZY database to identify the genes coding for glycoside hydrolases enzyme (Supplementary Table S2).

We first analyzed the macro metabolome composition of different gut segments and found that a total of 29,130 metabolites were identified in the samples from four gut segments, among which 1,027 metabolites showed significant differences in content in different gut segments (Figure 3A). Further, we found that cellulose content in different gut segments differed significantly (p < .01), ranking 31st as the metabolite with significant difference (Supplementary Table S3), and its content was lower in midgut and hindgut, but higher in foregut and foregut (Figure 3B). Interestingly, we also found that hemicellulose content showed a significant gradient of decrease in the four gut segments (p < .001), suggesting that metabolism is continuously decomposed with peristalsis of food in the intestinal tract (Figure 3C). Due to the limitations of our metabolomics approach, we were unable to identify differences in lignin, a complex woody component, in different gut segments. Thus, our results indicated that the lignocellulose composition of the A. swainsoni larva is distinct in the four gut segments.

According to the results of metabolomics analysis, there were significant differences in lignocellulosic content in four gut segments, so we analyzed the differences between lignocellulosic degradation-related genes and proteases in four gut segments. Metagenomic results (Figure 4A) showed that genes related to lignocellulose metabolism exist in different gut segments. In general, the abundance of related genes was high in the hindgut. According to CAZY annotation analysis (Figure 4B), we found that glucoside hydrolase was also widely distributed in different gut segments. Overall, the abundance of GH genes in the PHG segments was greater than in others, but the abundance of GH1 and GH5 was similar in different gut segments. Interestingly, according to functional Metaproteomic analysis, the expression of cellulase responsible for cellulose metabolism was higher in the midgut segment than in other segments, which may lead to the phenomenon of low cellulose content in the midgut segment found in the metabolomics results. Furthermore, analysis of the expression levels of enzymes related to hemicellulose and lignin metabolism showed that these enzymes showed high expression at the poles of different gut segments and low expression in the middle, which may lead to the gradient decrease of hemicellulose content found in the metabolomics results. It is worth noting that our analysis found that the expression level of GH enzymes identified in the Metaproteomic was significantly higher in the midgut than in other gut segments, similar to the distribution characteristics of cellulase, suggesting that the midgut is the main place for the metabolism of cellulose and related metabolites and sugars. Interestingly, according to the annotation of bin genomes, the genes coding for cellulase are prevalent in bins of four gut regions but are the most abundant in PHG which might explain the lowest abundance of cellulose in PHG (Figure 2C). Therefore, these multi-omics findings suggest that the gut microbiome of A. swainsoni in four gut segments may cooperate to attribute to lignocellulose deconstruction.

Here we aim to determine whether the different gut segments of A. swainsoni larvae are involved in nitrogen fixation. We first found that the C/N ratio of woody food of A. swainsoni larvae was significantly lower than that of frass. The C composition results were consistent with C/N ratio results (Figure 5A, B). Through analyzing the abundance of nitrogen-fixing genes in metagenomes of four gut segments (Figure 5C), we found that in general, nitrogen-fixing genes showed high abundance at the poles (FG and PHG) and low abundance in the middle of four gut segments (MG and AHG), similar to the distribution of hemicellulose degradation genes. Consistently, although nitrogen fixing genes are founded in all bins of four gut regions, the bins in PHG carry the most amount of nitrogen fixing genes. Unfortunately, due to the limitations of the metaproteomic technique, we were unable to identify nitrogenase expression. As a substitute, we analyzed the expression of the nitrogen-fixing core gene nifH (Zehr et al., 2003) in four gut segments (Figure 5D), and surprisingly found that nifH only expresses in AHG and PHG segments. In conclusion, metagenomic and gene expression results showed that nitrogen fixation genes mainly existed and expressed in two hindgut segments, suggesting a featured division of nitrogen fixation in gut segments of A. swainsoni larvae.