Riptortus pedestris colony (mtCOI GenBank accession no. OR654302) was collected at the Institute of Wheat Research, Shanxi Agricultural University, China, in 2022 and 2023. The colony was maintained on soybean plants (Glycine max L. cv. Jindou 36). The R. pedestris colony harbors a few gut symbiotic bacteria, such as Burkholderia sp. (Takeshita et al., 2014), Caballeronia insecticola (Jouan et al., 2024), and Serratia marcescens (Lee and Lee, 2022).

The soybean plants were grown in a potting mix. Two or three soybean plants were grown in 1.5 L pots to the six-to-seven true-leaf stage to raise R. pedestris. All R. pedestris colonies and soybean plants were maintained in climate-controlled chambers at 25 ± 2°C with a 16 h light: 8 h dark regime and 50–60% relative humidity (RH). The LED fluorescent lights were used, and the light intensity in the walk-in chamber was approximately 400 μmol/m² s.

Newly hatched 2nd, 3rd, 4th, and 5th nymphs and female adults were collected and divided into G1, G2, G3, G4, and G5 (Figure 1A; Supplementary Table S1). These samples’ surface were sterilized in 70% ethanol, and then the midguts of each sample were dissected in phosphate-buffered saline (1 × PBS) (Figures 1B–F), immediately stored in liquid nitrogen for ten minutes and then stored at −80°C for further analysis and metagenomic sequencing. Three biological replicates were performed for each group.

Five groups of samples were sent to Novogene Bioinformatics Technology Co. Ltd. (Beijing, China) for shotgun metagenome sequencing. Every group of samples was homogenized, and DNA was extracted according to the manufacturer’s recommended protocol for the QIAamp^® DNA Microbiome Kit (Cat. No. 51704, QIAGEN, Germany). This kit removes the host DNA to avoid affecting the sequencing result. The extracted gut genomic DNA was randomly sheared into 350 bp fragments, and sequencing libraries were generated. The obtained fragments were end-repaired, A-tailed and further ligated with an Illumina adapter. The fragments with adapters were PCR amplified, size selected, and purified. The library quality was checked with a Qubit 3.0 fluorometer (Life Technologies; Grand Island, NY, United States), real-time PCR was used for quantification, and a bioanalyzer was used for size distribution detection. The quantified libraries were pooled and sequenced on the Illumina NovaSeq X Plus platform (Illumina; San Diego, CA, United States) and the Agilent 2,100 (Agilent, Santa Clara, CA) system, according to the effective library concentration and data amount needed.

Metagenomic sequencing analysis was performed as described by Zhang Z. et al. (2022). Briefly, 20 metagenomic DNA libraries were constructed via paired-end 150 bp (PE150) sequencing. Considering the possibility of host contamination in samples, clean data needs to be mapped to the host database to filter out reads that may come from host origin via Bowtie2 software (Karlsson et al., 2013).

The richness and evenness of the gut bacteria were estimated via alpha diversity, including the Ace, Chao1, Shannon and Simpson indices, via QIIME 2 (Fung et al., 2021). Differences in the bacterial composition among the different developmental stages were calculated via principal coordinate analysis (PCoA) via nonparametric testing. Linear discriminant analysis effect size (LEfSe) (LDA score = 4, p < 0.05) was used to identify biomarkers with significantly different species between groups (Segata et al., 2011). Bar plots and heatmaps were generated for bacterial diversity and gene number differences in the different developmental stages via the R (version 4.0.3) package ggplot and pheatmap. To evaluate statistically significant differences between the groups, one-way ANOVA was performed with GraphPad Prism 7 at the 0.05 level.

The potential metabolic functions of the gut bacterial communities were predicted via DIAMOND (Feng et al., 2015). The commonly used functional databases currently include the Kyoto Encyclopedia of Genes and Genomes (KEGG), evolutionary genealogy of genes: Nonsupervised Orthologous Groups (eggNOG), and Carbohydrate-Active enzymes Database (CAZy). To predict the functions of the gut microbiota communities in each developmental stage, the Kyoto Encyclopedia of Genes and Genomes (KEGG) and carbohydrate-active enzymes (CAZy) databases were used to predict microbial functions and functional pathways.

Resistance genes are commonly present in the gut microbiota. Therefore, research on resistance genes has received widespread attention (Martínez et al., 2015). To further investigate the hazard mechanism of soybean stinkbugs, we conducted a statistical analysis of resistance genes at different developmental stages of R. pedestris via antibiotic resistance ontology (ARO) from the comprehensive antibiotic resistance database (CARD) (Jia et al., 2017).

The 2nd-, 3rd-, 4th-, and 5th- instar nymphs and female adults were collected and then surface sterilized in 70% ethanol. Individual insects were subsequently dissected in 1× PBS, and a midgut region with crypts was carefully cut out (Figures 2A,B). The individuals or groups of insects were homogenized via a JXFSTPRP-CL-48 homogenizer (Shanghai Jingxin Industrial Development Co., Ltd.). DNA was extracted as described previously (Kikuchi et al., 2005) with some modifications. The DNA was then quantified via a Micro Drop spectrophotometer (Bio-DL, Shanghai, China). The abundance of the gut symbionts Burkholderia, C. insecticola (thereafter, Caballeronia), Enterococcus sp. (thereafter, Enterococcus) and S. marcescens (thereafter, Serratia) were quantified via qPCR using a LineGene 9,600 Fluorescent Quantitative PCR Detection System (Bioer, Hangzhou, China) with 2 × SYBR Green master mix (Selleck, United States) and their primers were shown in Supplementary Table S1. The qPCR primers were designed with Primer Premier 5.0 software (Premier Biosoft International, Palo Alto, CA, United States). The thermal cycling protocol included initial denaturation for 30 s at 95°C, followed by 40 cycles of 5 s at 95°C and 30 s at 60°C, with melting curve analysis confirming that only the specific products were amplified. Two technical replicates were performed for each of the 6 biological replicates for R. pedestris.

qRT–PCR was used to quantify the expression of 3 genes involved in the ecdysone-related pathway. Total RNA was extracted from the detected gut samples via the TRIzol reagent method (Mei5bio, Beijing). The RNA was then quantified via a Micro Drop spectrophotometer. The RNA was first-strand reverse transcribed via a Vazyme kit (HiScript III RT SuperMix for qPCR (+gDNA wiper), Nanjing, China). The qRT–PCR primers were designed with Primer Premier 5.0 software. The primers used for the ecdysone-related genes were shown in Supplementary Table S1. qRT–PCR was performed as follows: denaturation at 95°C for 5 min, followed by 40 cycles of 15 s at 95°C, 30 s at 60°C and 30 s at 72°C. Two technical replicates were performed for each of the three biological replicates.

Standard curves were constructed, and the amplification efficiency was estimated as described previously (Luan et al., 2015). The elongation factor 1 α (EF1α) gene of R. pedestris served as an endogenous control gene for normalization via qPCR and qRT–PCR (Kim et al., 2015b). The Primer was shown in Supplementary Table S1. The relative symbiont density and gene expression were calculated via the 2^−ΔCt method (Schmittgen and Livak, 2008).

To survey the influence of gut bacteria on host development, 2nd-instar nymphs were treated with tetracycline hydrochloride (Biotopped, Beijing, China) to eliminate the gut symbiont Burkholderia. The 2nd- instar nymphs were collected from the hatching cage, and released into each feeding chamber, and then fed distilled water containing 0.05% ascorbic acid solution (w/v) supplemented with tetracycline hydrochloride at 1, 0.1 or 0.01 mg/mL for 6 days (Treatment). The control insects were administered 0.05% ascorbic acid solution not supplemented with antibiotics (CK). The artificial diets were renewed every day. Following the diet regimens, these treated nymphs were transferred to soybean plants. The R. pedestris strains with reduced Burkholderia titres (-BRp) were obtained via antibiotic treatment, and the control R. pedestris strains (+BRp) were obtained via the addition of an ascorbic acid solution that was not supplemented with antibiotics. The molting and developmental periods were statistically analysed. To evaluate the significant differences between the different treatments, t-tests were performed with GraphPad Prism 7 at the 0.05 level.

The guts from different developmental stages of R. pedestris (Figures 1A–F) were used for metagenomic sequencing. Twenty metagenomic DNA libraries constructed from the G1 to G5 groups were sequenced on the Illumina NovaSeq X Plus platform. As shown in Supplementary Table S1, there were 6,049.31–6,906.23 Mb of raw data from the G1 to G5 groups. After these data were filtered, 6,030.28–6,799.76 Mb of clean data were obtained from the G1 to G5 groups (Supplementary Table S1). The Q20s (%) of the G1 to G5 groups ranged from 98.10 to 99.03% (Supplementary Table S2). Furthermore, the GC contents (%) of the G1 to G5 groups ranged from 40.89 to 58.82% (Supplementary Table S2). From the G1 to G5 groups, the effective percentage (%) was greater than 98% in each group (Supplementary Table S2). A PCoA plot based on a phylum-level relative abundance profile revealed that axis 1 (PCoA1) explained 82.14% of the variability and that axis 2 (PCoA2) explained 12.78% of the variability. The PCoA plot revealed that the samples from different developmental stages were almost completely separated (Figure 1G). The above sufficient sequencing data indicated that these data were sufficient to estimate the gut bacterial diversity of R. pedestris fully.

We subsequently compared the number of genes associated with the gut microbiota at different developmental stages of R. pedestris. According to Spearman correlation coefficient analysis, gene abundance among samples was positively correlated, and there was a high correlation between different groups of samples (Supplementary Figure S1). In the five stages, the minimum number of genes was detected in the 4 N group compared with the other groups (Figure 1H, Wilcoxon rank-sum test, p < 0.001). However, there is no direct proportional relationship between the number of genes in a species and its richness. To screen for species with significant differences between groups, we conducted a phylogenetic analysis of the gut species between different groups via LEfSe at LDA > 4 and p < 0.05 (Figure 1I). In the 2 N and 3 N, the species richness was relatively limited. At the class level, Bacilli was more abundant in 2 N and at the order level, Lactobacillales was more ample in 3 N. In the 4 N and 5 N guts, there were individually 10 microbial communities at different taxonomic levels. At the order level, Hyphomicrobiales and Enterobacterales were more aboundant in G4 gut. At the family level, Burkholderiaceae was more plentiful in the 5 N midgut. Interestingly, the microbial communities in the FA intestine were concentrated at the phylum level (Figure 1I).

To further clarify the relative abundance of species in each sample, the abundances of gut microbial gene families and species diversity were quantified. The numbers of genes annotated from 2 N to FA were 26,860, 5,851, 505, 935, and 38,018, respectively, and 22,893 genes were collectively annotated (Figure 3A, Supplementary Table S4). Among all the annotated species, 3 genera had the highest abundance among the five groups, including Serratia, Enterococcus and Lactococcus. Caballeronia and Providencia presented high proportions of 2 N, 3 N and 5 N (Figure 3B).

The alpha diversity indices mainly include the Chao1 index, ACE index, Shannon index, and Simpson index. The ACE and Chao1 indices are commonly used to estimate unobserved species richness, whereas the Shannon and Simpson indices are used to measure species diversity. In the FAs, the number of unknown microorganisms was significantly greater than that in the G1 to G4 groups according to the ACE and Chao1 analyses (Figures 3C,D, p < 0.001; Supplementary Table S3), so the microbial community was relatively smaller than that in the other groups. According to Krona visualization analysis, up to 36% of the microorganisms in FA were unknown (Supplementary Figure S2E). However, according to the results of the Shannon and Simpson analyses, the microbial diversity of the 2 N–4 N groups were obviously greater than that of the other two groups (Figures 3E,F, p < 0.05; Supplementary Table S3). These results were consistent with the Krona visualization results (Supplementary Figures 2A–C).

To investigate the composition of the R. pedestris bacterial community at different groups, the top 10 most abundant genera and species were analysed. At the genus level, the advantageous genera were Caballeronia and Enterococcus, followed by Serratia and Lactococcus at 2 N, 3 N, and 5 N. In contrast, Caballeronia and Enterococcus were significantly enriched in the 5 N and groups, respectively (Figure 4A; Supplementary Table S5). There were five genera whose relative abundances were dominant in different groups, namely, Caballeronia, Enterococcus, Serratia, Bartonella and Lactococcus (Figure 4B; Supplementary Table S6). At the species level, Enterococcus faecalis was significantly greater in the FA group than in the other groups (Figure 4C; Supplementary Tables S7, S8). The genus Enterococcus and the species Enterococcus faecalis were quite abundant in FA, and their proportions reached 34 and 14%, respectively (Supplementary Figure S2E). In addition, at both the genus and species levels, the number of bacteria in the gut microbiota was significantly lower in the FA stage than in the other developmental stages, but many unassigned bacteria were found (Figures 4A,C; Supplementary Figure S2E). In the other groups, Enterococcus sp., Serratia sp. and Caballeronia sp. were present in a significant proportion (Figure 4D), so these species could be biomarkers in the stinkbug gut. These results revealed that the microbial community members changed with the developmental stage of R. pedestris.

To clarify the predominant species at different developmental stages, linear discriminant analysis (LDA) effect size analysis (LEfSe) was performed to detect notable differences in the gut bacteria across the life stages of R. pedetris. The results of the rank sum test (LDA > 4, p < 0.05) revealed that the gut microbiota was the most abundant in the 5 N (Supplementary Figures 3A,B). On the basis of the heatmap analysis, Burkholderia-Paraburkholderia-Caballeronia was the predominant genus in the 5 N (Supplementary Figure 3B). Interestingly, compared with that in the other stages, Burkholderia density was greater in the 5 N stage (Supplementary Figure S4). Therefore, bacteria of the Burkholderia genus have great potential for application in the prevention and control of stinkbugs.

The abundances of microbial metabolic pathways were predicted via KEGG, and the functional differences in the gut microbiome were compared. A PCA plot based on gene metabolic function revealed that the samples enriched with 2 N to FA were almost completely spatially separated (Figure 5A). The results of the KEGG statistical analysis revealed that the most abundant genes in the gut microbiota of all the samples were associated with metabolism, particularly carbohydrate metabolism, amino acid metabolism, metabolism of cofactors and vitamins, energy metabolism and nucleotide metabolism (Figure 5B; Supplementary Table S9). These results revealed that the midgut-colonizing microbiota could provide essential amino acids and vitamins for the host, which are lacking in leguminous crops. In addition, genes involved in environmental information processing, including membrane transport and signal transduction, were largely enriched in the gut microbiome of all the samples (Figure 5B). This result further explained that the gut microbiota of the stinkbugs were obtained mainly from the surrounding environment. From 2 N to FA, membrane transport and signal transduction, carbohydrate metabolism, and amino acid metabolism play important roles in the gut microbiota (Figure 5C). A comparison of the two groups revealed that, in the FA group, most of the microbiota functions were related to carbohydrate metabolism, membrane transport, nucleotide metabolism, signal transduction, replication and repair, and translation, which suggested that FA harbored the most abundant microbial functions (Figure 5D).

Owing to the highest abundance of genes associated with carbohydrate metabolism in the guts of the five groups (Figure 5B), we predicted carbohydrate-active enzymes (CAZy) of the gut microbiota on the basis of the CAZy database to further understand their activity in the gut microbiota. The CAZy (version 2023.03) database is a professional-level database for studying carbohydrate enzymes, covering six main functional categories: glycoside hydrolases (GHs), glycosyl transferases (GTs), polysaccharide lyases (PLs), carbohydrate esterases (CEs), auxiliary activities (AAs), and carbohydrate-binding modules (CBMs). A PCA plot based on gene enzyme function revealed that the samples enriched from 2 N to FA were almost completely spatially separated (Figure 6A). At the first CAZy classification levle, ove 66% genes of GHs (69.6%) and GTs (66.7%) were matched, indicating that these genes play a prominent role in the gut of stinkbugs (Figure 6B). Besides, the genes of CBMs (13.2%) are also essential enzyme families in the gut microbiome of stinkbugs (Figure 6B; Supplementary Table S10). At the second CAZy classification levle, GT4 and GT2 were significantly enriched in the GT family from 2 N to FA. GH13 and CBM50 were also relatively abundant in the gut microbiome of stinkbugs (Figure 6C; Supplementary Table S11). At the third CAZy classification levle, lipopolysaccharide N-acetylglucosaminyltransferase (EC 2.4.1.56) was the most abundant enzyme in five groups (Supplementary Table S12). Notably, the enzyme genes encoding the GH and GT families were significantly more abundant in the gut microbiota of different developmental stages of R. pedestris than other CAZy functional classes (Figure 6D). Therefore, GHs and GTs play important roles in the continuous development of stinkbugs.

Annotate the gut microbiome based on the CARD database, we further explored the relationship between antibiotic genotypes and phenotypes in different age groups of the stinkbug. The PCA plot based on ARG abundance revealed that almost all of the samples from G1 to G5 were spatially separated (Figure 7A). In different developmental stages of R. pedestris, ARGs were extremely prevalent and exceeded 16,000 (Figure 7B). A total of 47 resistance genes were predicted from 2 N to 5 N (Supplementary Table S13). The tet59 gene encoding tetracycline and the vanY gene encoding D, D-carboxypeptidase antibiotic resistance were most prevalent in the gut microbiota of 2 N (Figure 7C). The vanH gene, encoding the lactate dehydrogenase antibiotic resistance gene, was most abundant in the gut microbiome of 3 N. The types of ARGs including vanG, CRP, rsmA, FosA8, KpnH, msbA, emrR, SRT-3, KpnF, tet41, ImrD, and vanY were the most diverse in the gut microbiota of 4 N. The ARGs including vanH, adeF, and qacG were more abundant in the gut microbiota of 5 N. The ARGs including emeA, IsaA, dfrE, AAC6-Ic, vanT, and efrA were more abundant in the gut microbiota of FA (Figure 7C). Our results revealed a close relationship between the gut microbiota and ARGs.

The mechanisms of antibiotic resistance in bacteria include antibiotic efflux, antibiotic inactivation, antibiotic target alteration, antibiotic target protection, antibiotic target replacement and reduced permeability to antibiotics. We conducted a correlation analysis between the richness of bacteria and ARGs. Among the R. pedestris gut microbiota, Pseudomonadota and Bacillota were the principal phyla related to ARGs (Figure 7D). The bacteria from Pseudomonadota were associated with approximately 32 ARGs, particularly antibiotic efflux genes. Bacteria from Bacillota are associated with approximately 15 ARGs, especially antibiotic efflux genes and antibiotic inactivation genes (Supplementary Table S14).

To identify the effects of gut bacteria on the development of R. pedestris, the relative abundance of Burkholderia, Caballeronia, Enterococcus and Serratia, the number of molted insects and the life span of the insects were detected and analysed. The low Burkholderia titre line (–BRp) was generated by treatment with tetracycline hydrochloride at 1 mg/mL (Figure 2C, p < 0.01), whereas the other two treatments did not significantly affect the Burkholderia titre. The Burkholderia titre was reduced by 75% in the –BRp stinkbug group compared with the +BRp stinkbug group (Figure 2C). Likewise, Compared with the treatment without tetracycline hydrochloride (CK), titre of Caballeronia was significantly diminished by 74% by treatment with tetracycline hydrochloride at 1 mg/mL (Treatment) (Figure 2D, p < 0.05), however the titre of Enterococcus and Serratia were not significant changes between CK and treatment. We found that the molt number of the treated nymphs was significantly lower in the –BRp group than in the +BRp stinkbug group (Figure 2E, p < 0.001). Moreover, under the –BRp, 74EF, E75, and Er treatments, three ecdysteroid-related genes were downregulated. Among them, the expression of 74EF and E75 was significantly downregulated (Figure 2F). The comparison of the developmental periods of the –BRp and + BRp stinkbugs revealed that the –BRp stinkbug remained the 2nd nymph until it died on the 12th day, whereas the +BRp stinkbug successfully developed into an adult (Figure 2G). Taken together, our data demonstrated that low titres of Burkholderia strongly affected the development of R. pedestris, indicating their importance in the host gut.