We collected seeds of Ca. sclerophylla and Ca. tibetana infested with weevils during October and November 2022 when the leaves were falling and after the weevil larvae had parasitized the seeds. The seeds were kept in plastic boxes in the laboratory where emerging larvae were collected. Healthy C. bimaculatus larvae were collected from Ca. sclerophylla and Ca. tibetana. Curculio davidi larvae were collected from Ca. tibetana. Individual larvae from different host species were collected and labeled; thus, C. bimaculatus feeding on Ca. sclerophylla were labeled as BS; C. bimaculatus feeding on Ca. tibetana were labeled as BT, and C. davidi feeding on Ca. tibetana were labeled as DT (Table 1). The mature larvae were collected and surface sterilized by washing with ddH₂O for 3 min to remove surface dust, then soaking in 75% alcohol for 3 min to disinfect the body surface, rinsing with 2% NaClO solution for 30 s, and rinsing three times in ddH₂O, each time lasting 1 min, to remove residual reagents. After the treatment, each larva was placed in a sterile 2.0 ml centrifuge tube containing two small steel beads (φ = 3 mm) and stored at −80°C. Individual larvae were ground with a tissue grinder (60 Hz, 55–60 s; Tissuelyser-48, Shanghai Jinxin Industrial Development Co., Ltd.) for DNA extraction. Total DNA was extracted from each larva using the TIANamp Genomic DNA Kit (TIANGEN Biotech Co., Ltd.) according to the manufacturer’s instructions. DNA extracts were checked on 1% agarose gels. Since weevil larvae are very similar in appearance, we amplified the DNA fragments of the coding regions of mitochondrial cytochrome oxidase subunits I (COI) to identify each larva to species. The sequences were amplified using the PCR primers LCO1490 (5′-GGTCAACAAATCATAAAGATATTGG-3′) and HCO2198 (5′-TAAACTTCAGGGTGACCAAAAAATCA-3′) (Folmer et al., 1994). Amplification of mtDNA was performed as follows: an initial 1-min denaturation (98°C), followed by 30 cycles of 30 s at 94°C, 30 s at 48°C, and 1 min at 72°C, and a final extension cycle of 5 min at 72°C. All PCRs were performed in a total volume of 25 μL using 2 × Taq PCR MasterMix (TIANGEN Biotech Co., Ltd.), according to the manufacturer’s instructions.

The DNA of three larvae with the same labels was combined as one sample. Amplification of the V3-V4 regions of the 16S rRNA gene was performed using the primers 338F (5′-ACTCCTACGGGAGGCAGCA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) (Claesson et al., 2010). The PCR comprised an initial denaturation at 95°C for 3 min, followed by 30 cycles of 95°C for 30 s, 55°C for 30 s, 72°C for 45 s, and a final extension step at 72°C for 10 min. All PCRs were performed in a 20 μL total volume using TransStart^® FastPfu DNA Polymerase according to the manufacturer’s instructions. PCRs were performed in triplicate. The PCR products were purified, pooled in equimolar concentrations, and sequenced on the Illumina MiSeq PE300 platform (Illumina, San Diego, CA, USA) at the Majorbio Bio-Pharm Technology Co., Ltd (Shanghai, China). Paired-end sequences were merged into a single sequence with a length of ∼300 bp using FLASH 1.2.11. The obtained sequences were imported into the QIIME2. The DADA2 plugin in was used to denoise sequencing reads, resulting in high resolution amplicon sequence variants (ASVs) for downstream analysis. The consensus sequences for the ASVs were classified with a classify-sklearn classifier trained against the SILVA 16S rRNA gene reference (release 138.2) database.

For each group, DNA extracted from three individuals was selected and pooled to create a composite biological sample, with three such pools (representing three biological replicates) being collected for each group. In total, nine samples were used for metagenomic analyses of the weevils from the three groups. The composite samples were randomized such that the same sample names appeared in the microbiome and metabolomics analyses, but they represented different mixtures of individuals with the same labels. Based on the methods already described (Yuan et al., 2021; Zhou et al., 2023), libraries were prepared by fragmenting the extracted DNA to ∼350 bp using Covaris ME220 (Covaris, USA) and sequencing the paired-end reads on the Illumina HiSeq X-Ten platform at the Majorbio Bio-Pharm Technology (Shanghai, China), and generated a sequencing depth of approximately 12G per sample. Low-quality raw reads, sequences with low-quality bases (length < 50 bp or with a quality value < 20), and adapters were removed using Fastp v0.20.0.¹ About 12 Gb of raw data were obtained from each sample. MEGAHIT v1.1.2² was used to assemble high-quality reads, and contigs shorter than 300 bp were discarded. After gene prediction with Prodigal v2.6.3,³ CD-Hit67⁴ was used to cluster the sequences to create a non-redundant gene catalog with 90% sequence identity and 90% coverage. High-quality reads were aligned to the non-redundant gene catalogs to calculate gene abundance with 95% identity using SOAPaligner⁵ (version 2.21). Taxonomic analyses were performed by comparisons (BLASTp) against the NCBI-NR database at the e-value cutoff of 1e^–5 using Diamond v2.0.13.⁶ Based on the results of taxonomic annotation, a non-redundant gene catalog of bacteria was selected for functional annotation by comparisons (BLASTp) against the KEGG, and COG databases at the e-value cutoff of 1e^–5 using Diamond v2.0.13.

Using assembled metagenomic contigs with a minimum length of 1,000 bp, binning was performed on single-sample data using the CONCOCT tool Version 0.5.0.⁷ All resulting bins were subsequently dereplicated using dRep Version 2.2.9.⁸ The average nucleotide identity (ANI) analysis was performed based on the reference genome using fastANI version 1.0. Primary clusters were established based on a Mash ANI threshold of ≥ 90%. Secondary clustering was then performed at an ANI threshold of ≥ 99%, requiring a genome overlap of ≥ 10% between genomes. Quality evaluation of the resulting metagenome-assembled genomes (MAGs) was conducted using CheckM Version 1.0.12,⁹ following the criteria of completeness ≥ 50% and contamination < 10%. The obtained medium-quality MAGs and genome of one representative symbiotic bacterium Sodalis were further analyzed using antiSMASH to predict biosynthetic gene clusters (BGCs).

All seeds without weevil holes of every group were selected from the collected seeds for determination of metabolite contents. Each group of seeds collected were ground into the powder. All content determinations of the seed samples were carried out in triplicate.

The fat content was determined by the Soxhlet extraction method (Latimer and Latimer, 2023). Approximately 5 g of finely ground seed samples were placed in a thimble and subjected to continuous extraction with petroleum ether for 8 h. After extraction, the solvent was evaporated, and the remaining lipid content was weighed to determine the fat percentage.

The protein content was determined by the Kjeldah method (Hagerman and Butler, 1978). 2 g of seed powder was digested in concentrated sulfuric acid with a catalyst, converting organic nitrogen into ammonium sulfate. The digest was then neutralized with sodium hydroxide, and the released ammonia was distilled into a boric acid solution and titrated with hydrochloric acid. The nitrogen content was converted to protein content using a factor of 5.3.

Soluble sugar and starch contents were determined using anthrone-sulfuric acid colorimetry (Latimer and Latimer, 2023). In this assay, 2 g of seed material was hydrolyzed into simple sugars, which reacted with anthrone reagent in the presence of concentrated sulfuric acid, producing a blue-green complex. The absorbance of the reaction mixture was measured spectrophotometrically at 620 nm, and the concentrations of soluble sugar and starch were determined using calibration curves prepared with standard reference (glucose for soluble sugar, starch for starch).

The flavonoid content of the seeds was determined spectrophotometrically using rutin as a standard (Zhishen et al., 1999). A weighed portion of 1 g of seed material was extracted using ethanol, and the total flavonoid content was estimated by its reaction with Aluminum Nitrate Solution and Potassium Acetate Solution, forming a colored complex. The absorbance was recorded at 420 nm, and the flavonoid concentration was calculated based on a standard curve constructed with rutin.

The tannin content of the seeds was determined spectrophotometrically using gallic acid as a standard (Latimer and Latimer, 2023). 2 g of seed powder was extracted using acetone-water, and tannins were quantified based on their reaction with Sodium Tungstate-Sodium Molybdate Mixed Solution and Sodium Carbonate Solution, producing a color change proportional to tannin concentration. The absorbance was measured at 765 nm, and tannin content was determined using a gallic acid standard curve.

The protein content of Ca. tibetana was significantly higher than that of Ca. sclerophylla. The contents of soluble sugars, flavonoids, and tannins of Ca. sclerophylla were significantly higher than those of Ca. tibetana. There were no significant differences in fat or starch contents between the two tree species (Table 2). Data were analyzed by t-test using SPSS software (SPSS Inc., USA). Statistical significance was set to P < 0.05 (Rakić et al., 2006).

After subsampling each sample to an equal sequencing depth, 42,883 sequences were obtained per sample. After annotation, 599 ASVs were obtained belonging to 17 phyla, 26 classes, 67 orders, 110 families, and 204 genera.

An association network was used to visualize the relationships between ASVs and the different host insect or host plant combinations. Of these ASVs, only 2.84% were shared among all groups; 88.6% were associated with only one group. Approximately 1.5% of the ASVs were shared by only two groups from the same host plant Ca. tibetana, and 5.01% were shared by two groups of the same host insect C. bimaculatus. There were 12.69% and 19.03% unique ASVs detected in the groups from C. davidi and C. bimaculatus, respectively, feeding on Ca. tibetana, and 56.93% unique ASVs were found in C. bimaculatus feeding on Ca. sclerophylla (Figure 1A). The α-diversity of the samples at the ASV level was estimated by the Chao, Ace, Simpson, and Shannon indices (Supplementary Figure 1A and Supplementary Table 1). The Chao or Ace indices of C. bimaculatus feeding on Ca. sclerophylla (BS) were higher than the other two groups (Kruskal–Wallis H test, P < 0.05). The Simpson indices of C. davidi and C. bimaculatus feeding on Ca. tibetana (DT and BT) were both significantly higher than that of C. bimaculatus feeding on Ca. sclerophylla (BS; Kruskal–Wallis H test, P < 0.001), and the Simpson indices of C. davidi feeding on Ca. tibetana and C. bimaculatus feeding on Ca. tibetana were significantly lower than that of C. bimaculatus feeding on Ca. sclerophylla (Kruskal–Wallis H test, P < 0.01). There were no significant differences between the Ace or Chao indices of weevils feeding on the same plants Ca. tibetana, but the Shannon and Simpson indices were significantly different between the two groups (Kruskal–Wallis H test, P < 0.05). The α-diversity suggested that the bacterial community richness levels and the bacterial community diversity in the weevils feeding on Ca. sclerophylla was greater than that in the weevils feeding on Ca. tibetana.

The association network of genera showed that the Proteobacteria, Firmicutes, and Actinobacteriota were the three most abundant phyla in the bacterial communities in weevils feeding on Fagaceae seeds (Figure 1B). Of these genera, 32 were shared between all groups; three genera were only shared by groups from the same host plant Ca. tibetana; 29 genera were only shared by groups from the same host insect C. bimaculatus. The C. bimaculatus feeding on Ca. sclerophylla had the largest number of unique genera of 74. There were 11 main genera with a total abundance > 1% in nine samples (Supplementary Figure 1B and Supplementary Table 2). The prevalence of Sodalis, Candidatus Curculioniphilus, Methylorubrum, Rickettsia, Paenibacillus and Sphingomonas was 100% (Supplementary Table 2). However, the genera were distributed unevenly among the different groups. Candidatus Curculioniphilus, the primary endosymbiont of Curculio, had an average abundance > 2.6% in all groups. The Ca. sclerophylla feeders showed a greater intrasample bacterial diversity of genera than the Ca. tibetana feeders. However, the bacterial composition of C. bimaculatus differed from that of C. davidi. The β-diversity of samples at the genus level showed that two groups of C. bimaculatus were clustered in composition and community structure, while C. davidi was an independent cluster (Supplementary Figure 1C). The main genera of C. davidi were different from those of C. bimaculatus.

Spearman’s correlations between bacterial genera with a total abundance > 1% are shown in Figure 2A. The prevalence of Sodalis, Candidatus Curculioniphilus, Methylorubrum, Rickettsia, and Paenibacillus was 100%. The abundance of Sodalis was significantly positively correlated with Rickettsia. Curculio davidi had higher abundances of Sodalis and Rickettsia than C. bimaculatus (Supplementary Figure 2). The abundance of Candidatus Curculioniphilus was significantly positively correlated with Stenotrophomonas. The abundance of Methylorubrum was significantly positively correlated with Pantoea, Wolbachia, Serratia, Paenibacillus, and unclassified Enterobacteriaceae. The Ca. sclerophylla feeders showed a high abundance of Methylorubrum, Pantoea, Wolbachia, and Serratia. The abundance of Pantoea was significantly negatively correlated with most of the main genera. The Ca. sclerophylla feeders had the largest number of genera.

A correlation analysis between the main genera and metabolite contents was visualized using a heatmap (Spearman’s correlations; Figure 2B). Among the four genera with 100% prevalence, the abundances of Paenibacillus and Candidatus Curculioniphilus had no significant association with the metabolite contents. The abundances of Methylorubrum and Pantoea in Ca. sclerophylla feeders showed a significant negative association with the protein content, and Methylorubrum was positively associated with soluble sugars and flavonoids (P < 0.05). Pantoea was positively associated with tannins (P < 0.05). The abundance of Rickettsia in C. davidi was significantly negatively associated with soluble sugars and flavonoids (P < 0.05). There were no genera significantly associated with fat or starch contents.

A total of 6,210,941 non-redundant genes were identified in the three groups of samples. After searching for non-redundant gene sequences against the NR database, 99,898 genes were annotated for bacteria. Bacteria comprised the largest percentage of all samples (Supplementary Figure 3A). We selected genes of bacteria based on the results of taxonomic annotation for functional annotation.

In the COG functional annotations of bacteria, 2,210 COGs belonging to 24 functions were predicted for all samples (Supplementary Table 3). The function with the largest number of COG annotations was [L]Replication, recombination, and repair, followed by [J]Translation, ribosomal structure and biogenesis, [M]Cell wall/membrane/envelope biogenesis, [G]Carbohydrate transport and metabolism, and [E]Amino acid transport and metabolism (Figure 3A). The category with the highest abundance in the COG annotations was cellular processes and signaling, followed by metabolism (Figure 3B). The most abundant function, [X]Mobilome: prophages, transposons, was the same for each group.

In the C. bimaculatus groups, the next highest abundances were [L]Replication, recombination, and repair and [J]Translation, ribosomal structure, and biogenesis. However, in the C. davidi group, the next highest abundances were [G]Carbohydrate transport and metabolism and [E]Amino acid transport and metabolism, suggesting that the metabolic activity of bacteria in C. davidi was stronger (Figure 3C).

For the KEGG functional annotations of bacteria, 2,744 KOs (KEGG Orthology) were predicted for all samples (Supplementary Table 4). The KEGG pathway with the largest number of reads was global and overview maps, followed by amino acid metabolism and carbohydrate metabolism (Figure 4A). The proportions of KEGG level 1 pathways were similar across the groups. The most abundant level 1 pathway of each group was the metabolism pathway (Figure 4B). Figure 4C shows the composition of 15 main level 2 pathways in each group. The most abundant level 2 pathway of each group was global and overview maps. The C. bimaculatus groups had higher abundances of the amino acid metabolism and nucleotide metabolism pathways. However, C. davidi had a higher abundance of the pathways of carbohydrate metabolism and membrane transport. The similar annotation results of COG and KEGG suggest that the C. davidi bacterial community has stronger carbohydrate metabolic activity. Overall, the composition of functional genes of C. bimaculatus feeding on Ca. sclerophylla and C. davidi feeding on Ca. tibetana was relatively similar in structure (Supplementary Figure 3B).

The abundance of KEGG level 2 pathways of metabolism was compared between the species to investigate the relative abundances of bacterial genera involved in metabolism-related symbiotic functions in the weevils. There were significant between-group differences in the abundances of five pathways (Figure 5A). The bacteria of Ca. sclerophylla feeders harbored a significantly higher abundance of global and overview map genes (Kruskal–Wallis H test, P < 0.05). The bacteria of C. bimaculatus harbored a significantly higher abundance of amino acid metabolism and nucleotide metabolism genes (P < 0.05). The bacteria of C. davidi harbored a significantly higher abundance of carbohydrate metabolism genes (P < 0.05). The group of C. bimaculatus feeding on Ca. tibetana had a higher relative abundance than the groups of C. davidi feeding on Ca. tibetana and C. bimaculatus feeding on Ca. sclerophylla in the biosynthesis of other secondary metabolites pathway (P < 0.05).

There was no significant variation in the high abundance of genera involved in five metabolic pathways (Figure 5B). Sodalis, Wolbachia, Klebsiella, and Rickettsia in the phylum Proteobacteria were the main contributors to these pathways. The contribution of Klebsiella was reduced in nucleotide metabolism, carbohydrate metabolism, and biosynthesis of other secondary metabolites pathways. Klebsiella contributed to pathways of global and overview maps and amino acid metabolism of C. bimaculatus but had lesser contributions in C. davidi. In C. davidi, Sodalis with the highest abundance (Supplementary Figure 2) was the main contributor to these pathways. Wolbachia contributed low levels to the amino acid metabolism and carbohydrate metabolism pathways of Ca. tibetana feeders while being the leading contributor to nucleotide metabolism and biosynthesis of other secondary metabolites pathways. The contributions of Wolbachia in the two groups of C. bimaculatus were similar.

The metagenomic binning yielded four medium-quality non-redundant MAGs, three of which were associated with BGCs. No MAGs were identified in the C. bimaculatus feeding on Ca. tibetana (Figure 5C). The MAG unrelated to secondary metabolite biosynthetic gene clusters was classified as genus Wolbachia and was exclusively found in the C. bimaculatus feeding on Ca. sclerophylla. The three BGC-associated MAGs were classified as genus Sodalis. AntiSMASH analysis identified three types of BGCs within these MAGs: hserlactone, phosphonate, and T3PKS (Supplementary Table 5). Compared to the representative symbiotic bacterium Sodalis pierantonius of rice weevil Sitophilus oryzae (SOPE), the three Sodalis MAGs identified in this study lacked thiopeptide-BGCs, but contained an additional T3PKS-BGCs.

The patterns of variation in the abundances of 50 key enzymes were identified to understand the functions of enzyme responses in weevil metabolism. Almost all key enzymes were abundant in the groups of C. davidi feeding on Ca. tibetana and C. bimaculatus feeding on Ca. sclerophylla, while the group of C. bimaculatus feeding on Ca. tibetana exhibited low abundance of the metabolic enzymes (Figure 6A). Terpenoids, glucosinolates, and alkaloids are among the best-studied classes of plant secondary compounds. Thus, two key pathways involved in the biosynthesis of plant secondary compounds were also compared. The abundance of ko00960 (tropane, piperidine, and pyridine alkaloid biosynthesis) was significantly higher in C. davidi (Kruskal–Wallis H test, P < 0.05). In contrast, the abundance of ko00232 (Caffeine metabolism) in the group of C. bimaculatus feeding on Ca. tibetana was significantly higher than that in the groups of C. davidi feeding on Ca. tibetana and C. bimaculatus feeding on Ca. sclerophylla (Figure 6B). Tyrosine and phenylalanine are key aromatic amino acids involved in the exoskeleton development of beetles. The abundance of ko00350 (Tyrosine metabolism) and ko00360 (Phenylalanine metabolism) pathways was compared. C. davidi had a significantly higher abundance of these two key pathways. For C. bimaculatus, Ca. sclerophylla feeders had significantly higher abundances of the two pathways (Figure 6C).