Sequences and samples quality control
For the 286 stink bug PCR samples of the bacterial 16S rRNA V4 sequence, the Illumina MiSeq sequencing produced approximately 7.5 million reads in one lane, and the mean number of raw reads ranges from about 7,000 to 53,000 per sample for different locations. Reads filtering procedures including merging reads and trimming low-quality bases removed large percentages of reads from most locations, with 33%-73% reads remaining, but the stringency of filtering left high-quality sequence data for subsequent analyses (Supplemental Table 2). After OTU identification, low-quality samples were filtered out for downstream analysis, leaving a total of 192 samples out of 286, with 75 samples from E. heros collected from 10 sites in Brazil and 117 samples from P. guildinii collected from three sites in Brazil and four sites in the US (Supplemental Table 2).
The sufficiency of sequencing coverage was evaluated by alpha diversity and rarefaction analyses. The depth of sequencing was sufficient as shown by P. guildinii rarefaction plot that began to plateau by 1000–1500 reads per sample (Supplemental Fig. 1B). Similarly, the E. heros rarefaction plot also began to plateau by 1000–1500 reads per sample for most locations, although the Teresina and Uberlandia locations might have benefitted from additional sequencing as those plotlines still had a slightly positive slope at 1500.
Microbiome characterization of Euschistus heros
Comparing the different collection sites of E. heros using the non-metric multidimensional scaling (NMDS) analysis, the bacterial communities present in the gut of these stink bugs clustered according to the following ecoregions of Brazil: tropical dry forests, tropical moist forests, and tropical savanna (Fig. 1A and B, PERMANOVA P = 0.0001). Therefore, subsequent analyses were based on these ecoregions.
Three metrics were used to estimate the alpha diversity of the microbiome among all the stink bug samples, i.e. the number of observed OTUs, Chao1, and Shannon. The species richness, based on the number of observed OTUs, is different among the ecoregions (p-value = 0.015). E. heros collected in tropical dry forests (p-value = 0.0059) and tropical savanna (p-value = 0.034) had significantly higher observed OTUs than those collected in tropical moist forests. This result was confirmed by the richness estimator Chao1 (p-value = 0.0004). E. heros collected in sites from tropical dry forests (p-value = 3.8e-05) and Tropical Savanna (p-value = 0.012) had more richness of species than E. heros collected in tropical moist forests. However, Shannon’s diversity, an index that takes into account both abundance and evenness of the species present, did not show significant difference among ecoregions. A noticeable observation would be that the E. heros microbiome from tropical dry forest showed the highest richness of species (OTUs and Chao1) but the lowest Shannon diversity index, indicative of low evenness/equitability of species within this bacterial community (Fig. 1C).
The most abundant bacterial families in the intestinal microbiome of E. heros collected in tropical dry forests were Bradyrhizobiaceae, Enterobacteriaceae, and Sphingomonadaceae. Additionally, E. heros from this ecoregion contained the lowest number of bacterial families. The most abundant bacterial families identified from insects collected in the tropical moist forests region were: Acetobacteraceae, Bacillaceae, Bradyrhizobiaceae, Enterobacteriaceae, Enterococcaceae, Hyphomicrobiaceae, Methylobacteriaceae, Polyangiaceae, and Sphingomonadaceae. Among these families, the Acetobacteraceae and Bacillaceae were unique to insects sampled from the tropical moist forests. From the tropical savanna insects, we identified the top bacteria families being the Bradyrhizobiaceae, Enterobacteriaceae, Enterococcaceae, Hyphomicrobiaceae, Methylobacteriaceae, Moraxellaceae, Polyangiaceae, and Sphingomonadaceae (Fig. 1D).
Comparing the bacterial families across ecoregions, the Bradyrhizobiaceae, Enterobacteriaceae, and Sphingomonadaceae were found in the intestinal microbiome of insects from all three regions. However, striking differences in abundances could be observed for the Bradyrhizobiaceae and the Enterobacteriaceae. More specifically, stink bugs from the tropical moist forests had a significantly higher mean proportion of the Bradyrhizobiaceae than dry forests (p-value = 0.012) and savanna (p-value = 0.016) (Fig. 1E). Although Enterobacteriaceae was the most abundant family in stink bugs from all three ecoregions (Fig. 1D), it was more abundant in the tropical dry forests, followed by the tropical savanna and moist forests (Fig. 1E). The Moraxellaceae only existed abundantly in stink bugs from tropical savanna (~ 20%) whereas in the other two ecoregions they were not detected, or present at a very low level.
Recent discoveries indicated that E. heros from the north and the south of Brazil were composed of two distinguished lineages. With the expansion of soybean in the Brazilian Cerrado, these two lineages could meet and produce hybrids more adapted to the climates and even to new hosts, such as cotton [39, 40]. Thus, we performed analyses to investigate the relationships between the stink bug intestinal microbiomes and the proposed lineages of E. heros: north (Teresina/PI, Palmeirante/TO), south (Piracicaba/SP, Anhembi/SP and Ponta Grossa/PR), and hybrids (the other localities). The NMDS analysis showed that the bacterial species present in our samples clustered significantly by stink bug lineages (PERMANOVA p-value = 0.0001) (Supplemental Fig. 2A). The species richness, based on Observed OTU and Chao1 index, is higher in samples from the north lineage than the south lineage, and the hybrids have an intermediate species richness. The species diversity, based on the Shannon index, is bigger on the north and hybrid lineages than the south lineage of E. heros (Supplemental Fig. 2B).
Analyzing the relative abundance of the 10 most abundant bacterial families in each lineage, we can see that the hybrid lineages had more bacterial families than the north and south lineages and the proportion of these families present in the gut of the hybrids lineages are more evenly distributed than in pure lineages (Supplemental Fig. 2C). Furthermore, the family Bacillaceae is unique in the hybrids lineages and the Acetobacteraceae family is unique in the south lineage.
Comparing the bacteria families among lineages, the family Bradyrhizobiaceae was the most abundant in E. heros from hybrids lineages than either of the pure lineages. The family Hyphomicrobiacea was the most abundant in E. heros from hybrid lineages than north lineages. The families Moraxellaceae and Sphingomonadaceas were the most abundant in hybrid lineages than south lineages and the Enterobacteriaceae family was the most abundant in both pure lineage than hybrids (Supplemental Fig. 3).
Microbiome characterization of Piezodorus guildinii
The comparison of the different collection sites of P. guildinii by NMDS analysis (Fig. 2B) revealed that the bacterial species present in the gut of these stink bugs also clustered (p-value of 0.0001) according to the ecoregions of Brazil and the United States: temperate broadleaf forests, temperate conifer forests, temperate grasslands, and tropical savanna (Fig. 2A). Therefore, subsequent analyses were based on these ecoregions.
The bacterial species diversity across ecoregions were significantly different based on the alpha diversity index, the number of OTUs, and Chao1 index (p-value = 0.0072 and p-value = 0.0035, respectively) (Fig. 2C). On average, P. guildinii collected in temperate grasslands had more bacterial species than P. guildinii collected in the other ecoregions. Considering the richness estimator Chao1, the values were also significantly different (p-value = 0.0035). P. guildinii from tropical savanna were richer in species than insects from temperate broadleaf forests (p-value = 0.0025). The Shannon’s diversity showed a significant difference among ecoregions (p-value = 0.018). P. guildinii from temperate conifer forests (p-value = 0.015) and tropical savanna (p-value = 0.044) had higher equitability than temperate broadleaf forests, confirming the lowest species diversity from insects of temperate broadleaf forests at Shannon index (Fig. 2C).
When focusing on the microbial abundance of P. guildinni at the family level, P. guildinii from temperate broadleaf forests had mostly gut microbiota bacteria of the families Acetobacteraceae, Aurantimonadaceae, Bradyrhizobiaceae, Enterobacteriaceae, Enterococcaceae, Moraxellaceae, Rhodocyclaceae, Sinobacteraceae, and Sphingomonadaceae. Among these families, the Sinobacteraceae family was found only in insects from temperate broadleaf forests and it existed at a low proportion (Fig. 2D). Six major families Acetobacteraceae, Aurantimonadaceae, Bradyrhizobiaceae, Enterobacteriaceae, Moraxellaceae, and Sphingomonadaceae were found in the insect guts from temperate conifer forests and one major family, Enterobacteriaceae, was found in the intestinal microbiome of insects from temperate grasslands forests. The insects collected in tropical savanna had in the intestinal microbiome the families Acetobacteraceae, Aurantimonadaceae, Bradyrhizobiaceae, Deinococcaceae, Enterobacteriaceae, Enterococcaceae, Moraxellaceae, Rhodocyclaceae, and Sphingomonadaceae. The Rhodocyclaceae family was found only in insects from this ecoregion and temperate broadleaf forests in low proportion (Fig. 2D).
Comparing the bacterial families present in the intestine of P. guildinii insects shared by more than one ecoregions, the Acetobacteraceae family was significantly more abundant in insect guts from temperate broadleaf forests than insect guts from tropical savanna with a p-value at 0.044. The Aurantimonadaceae was more abundant in P. guildinii from temperate conifer forests than temperate broadleaf forest and tropical savanna (p-value 0.015 and 0.029). The Bradyrhidobiaceae family was more abundant in P. guildinii from temperate broadleaf forests than tropical savanna. Lastly, the Enterobacteriaceae family, as a major component in all four ecoregions, exhibited a lower proportion in temperate conifer forests than the other three ecoregions, but the difference did not reach statistical significance either, due to large variation (Fig. 2E).
Since the species P. guildinii co-occurred in both Brazil and the US, we further explored the association between its intestinal microbiome and the country of origin. The NMDS analysis (Fig. 3A) showed distinct groups of bacterial species in insects from Brazil and the US with a p-value of 0.0107. According to the number of observed OTU (p-value = 0.015) and the richness estimator Chao1 (p-value = 0.0087), P. guildinii sampled from Brazil were richer in intestinal bacterial species than insects from the US (Fig. 3B).
The most abundant bacterial families in P. guildinii collected in Brazil were Acetobacteraceae, Aurantimonadaceae, Bradyrhizobiaceae, Deinococcaceae, Enterobacteriaceae, Enterococcaceae, Moraxellaceae, Rhodocyclaceae, and Sphingomonadaceae. In comparison, the most abundant bacterial families in P. guildinii collected in the United States were Acetobacteraceae, Aurantimonadaceae, Bradyrhizobiaceae, Enterobacteriaceae, Enterococcaceae, Moraxellaceae, Rhodocyclaceae, Sinobacteraceae, and Sphingomonadaceae (Fig. 3C). However, the Deinococcaceae family occurred only in stink bugs collected from Brazil, and the Sinobacteraceae occurred only in insects collected from the US (Fig. 3C). The families Acetobacteraceae, Bradyrhizobiaceae, and Moraxellaceae were significantly more abundant in P. guildinii from the United States than Brazil, and the family Rhodocyclaceae was more abundant from Brazilian insects (Supplemental Fig. 4). Enterobacteriaceae, the most abundant bacterial family in both countries, constituted about 80% of the gut microbiome of both countries and showed no significant difference (Fig. 3C; Supplemental Fig. 4).
Analysis of E. heros and P. guildinii combined
We further combined the microbiome data of E. heros and P. guildinii and clustered samples by ecoregions. With PERMANOVA, the bacterial species present in stink bugs displayed distinct compositions by different ecoregions with a significant p-value at 0.0001(Fig. 4B). As seen in the NMDS analysis (Fig. 4B), samples from tropical dry forests clearly grouped together and were separated from samples from tropical moist forests and tropical savanna. Therefore, subsequent analyzes were based on ecoregions.
According to the number of observed OTUs, we found that stink bugs from tropical dry forests contained more abundant bacterial species than stink bugs from the tropical moist forests (p-value = 0.0059) and tropical savanna (p-value = 0.041). The same occurs when the richness estimator Chao1 is compared (p-value = 3.8e-05 and p-value = 0.0042, respectively). Lastly, Shannon’s index showed that the species diversity of all three ecoregions did not differ significantly from each (Fig. 4C).
The most abundant bacterial families present in the digestive tracts of E. heros and P. guildinii in tropical dry forests were Bradyrhizobiaceae, Enterobacteriaceae, and Sphingomonadaceae, with Enterobacteriaceae comprising more than 80% of the bacterial community. In turn, the abundant families in insects from tropical moist forests were Acetobacteraceae, Bacillaceae, Bradyrhizobiaceae, Enterobacteriaceae Enterococcaceae, Hyphomicrobiaceae, Rhodocyclaceae, and Sphingomonadaceae (Fig. 4D). Insects from tropical savanna had gut bacteria of the families Acetobacteraceae, Bacillaceae, Bradyrhizobiaceae, Deinococcaceae, Enterobacteriaceae Enterococcaceae, Hyphomicrobiaceae, Moraxellaceae, Rhodocyclaceae, and Sphingomonadaceae. The families Deinococcaceae and Moraxellaceae only existed in the tropical savanna group (Fig. 4D). Comparing the bacterial families present in two or more ecoregions, the Bradyrhizobiaceae family was significantly more abundant in insects from tropical moist forests than tropical dry forests and tropical savanna. Enterobacteriaceae family was most abundant in insects from tropical dry forests, followed by insects from tropical savanna (Supplemental Fig. 5).
Comparing between the species of stink bugs E. heros and P. guildinii, the NMDS analysis showed that these pentatomids had different bacterial species in their intestines, which can be separated into two groups with a p-value equal to 0.0001 (Fig. 5A). Analyzing the diversity indexes, the number of observed OTUs showed that E. heros had more bacterial species than P. guildinii. Based on Shannon’s index, the bacterial community in the gut of E. heros had significantly higher diversity than P. guildinii, which also suggests a relatively even distribution of species in the E. heros intestines (Fig. 5B).
The composition of bacterial communities in the two stink bugs was similar at the family level (Fig. 5C). Enterobacteriaceae was the major family and constitute more than 50% of the bacterial species in both stink bugs. The two other major bacterial families shared by the two stink bug species were Bradyrhizobiaceae and Moraxellaceae. However, differences can be observed as well. Hyphomicrobiaceae occurred only in E. heros and Deinococcaceae occurred only in P. guildinii (Supplemental Fig. 6).
The phylogenetic relationship of the OTUs from the top 10 most abundant bacterial families was identified to understand their distribution (Fig. 6). Analyzing the OTUs among P. guildinii from Brazil and the United States, 13 unique OTUs belonged only to Brazil and 19 unique OTUs belonged only to the United States. Most OTUs were not related among countries. The families Enterococcaceae (2 OTUs), Acetobacteraceae (1 OTU), Enterobacteriaceae (15 OTUs), Aurantimonadaceae (1 OTU), and Bradyrhizobiaceae (1 OTU) belonged to P. guildinii from both countries. This difference may be related to bacterial diversity among environments. On the other hand, comparing the OTU from the ten most abundant bacterial families among E. heros and P. guildinii from Brazil, 16 unique OTUs belong only to E. heros and 7 unique OTUs belong only to P. guildinii. The families Enterococcaceae (2 OTUs), Bacillaceae (2 OTUs), Acetobacteraceae (1 OTU), Rhodocyclaceae (2 OTUs), Moraxellaceae (3 OTUs), Enterobacteriaceae (12 OTUs), Sphyngomonadaceae (7 OTUs), and Bradyrhizobiaceae (3 OTUs) belonged to both stink bug species. This difference is most interesting, E. heros had twice the number of unique bacteria than P. guildinii, which may be favoring the adaptation of this species over the other since it is the most difficult to control today (Fig. 6).