Dissecting the Pandora’s box: preliminary phylogenomic insights into the internal and external relationships of stink bugs (Hemiptera: Pentatomidae)

Abstract

Introduction

The stink bugs (Pentatomidae) compose the fourth most speciose family of true bugs. Their astonishing ecological and phenotypic diversity make the group ideal to investigate macroevolutionary questions related to phenotypic evolution and its underlying genetic bases (Rebagliati et al. 2005, Genevcius et al. 2017, 2022, Roell and Campos 2019). Nevertheless, our capability of approaching these questions using stink bugs as models has been severely hampered for decades due to the absence of a stable and well resolved phylogenetic hypothesis.

The first formal phylogenetic study to include a reasonable amount of stink bug species was Grazia et al. (2008). The authors found poor support for both the position of Pentatomidae as well as its internal relationships. Several studies in the following decade included stink bugs as terminals, but focused on relationships among the families (e.g., Wu et al. 2016; Liu et al. 2019). Analyses using Sanger sequencing found strikingly incongruous results, especially regarding the position of Pentatomidae. The first molecular phylogeny focusing on stink bug taxa was released in 2017, although with very strict taxonomic range (Bianchi et al. 2017). Genevcius et al. (2021) further inferred a phylogeny for the tribe Chlorocorini that included a reasonable range of tribes and subfamilies as outgroups. The authors found major conflicts with the groupings proposed in the traditional taxonomic classification in tribes and subfamilies (Rider et al. 2018).

A milestone for the stink bug systematics was the molecular study by Roca-Cusachs et al. (2022). They studied more than a hundred species encompassing several lineages with four molecular markers, representing the first attempt to investigate the phylogenetic history of stink bugs as a whole (Fig. 1A). In their analyses, well supported clades were those more terminal in the phylogeny, especially those among species of the same genus. Unfortunately, very few clades at higher levels (tribes, subfamilies) were well supported, even considering a moderate bootstrap threshold (Fig. 1B). A common reason for this lack of resolution is an insufficient amount of molecular data, which usually yields poorly resolved trees at deeper nodes in ancient groups such as the Pentatomidae. Therefore, a solid phylogeny-based classification for the stink bugs will most likely be possible only with the inclusion of genome-scale data. Here, I tackle this issue by conducting the first phylogenomic study focusing on the internal and external relationships of stink bugs. By using exclusively public data, my main goal was to provide a preliminary assessment of how genome-scale data may help overcome the above-mentioned issues. My specific aims are 3-fold: (i) to test if genome-scale data may recover well supported clades in different hierarchical levels (e.g., genera, tribes, etc.); (ii) to test how a phylogenomic tree compare to recent Sanger-derived phylogenies, especially for the poorly supported clades from these works; (iii) to attempt to resolve the controversial position of Pentatomidae among other families of Pentatomoidea.

Fig. 1.

Phylogenetic hypotheses for the stink bugs: Sanger-derived tree adapted from Roca-Cusachs et al. (2022) (A and B), and the phylogenomic trees resulted from the present work (C and D). A) a simplified version of the ML tree by Roca-Cusachs et al. (2022), where some species were pruned for visual purposes and to highlight relationships among major groups; monophyletic tribes and subfamilies are reduced to triangles. B) same tree as A), but clades with bootstrap supports below 70 (Roca-Cusachs et al. 2022) are collapsed. C) ML phylogenomic tree inferred from a supermatrix approach (i.e., genes concatenated into a single alignment) in IQ-TREE. D) coalescent phylogenomic tree inferred from individual gene trees in ASTRAL. For C) and D), only clade supports below 100/1.0 are exhibited on the right side of nodes; diamonds in terminal branches are species with reference genomes, while circles represent transcriptomic data.

Open in new tabDownload slide

Materials and Methods

I conducted this phylogenetic study using exclusively public genome-scale data, including transcriptomes and high-quality reference genomes. The final dataset comprised amino acid data of benchmarking universal single-copy orthologs (BUSCO). For the species with RNA-seq data available, I used seven already assembled transcriptomes and assembled 19 transcriptomes de novo from raw reads. I conducted the assembly in rnaSPAdes v. 3.15.4 (Bushmanova et al. 2019) with default parameters (Supplementary Table S1). Resulting hard filtered transcripts were submitted to a filtering step, aiming at removing shorter isoforms and reducing transcript redundancy. I did the latter by applying a clustering algorithm in CD-HIT v. 4.8.1 (Fu et al. 2012) using a 90% identity cutoff. For the 10 remaining species, I used the released reference genome assemblies. The final dataset included a total of 40 terminal taxa (Supplementary Table S1), from which 34 were ingroups (Pentatomoidea) and six were outgroups (Lygaeoidea, Coreoidea, Miroidea, Aradoidea, and Miroidea). I rooted the tree in Miridae.

I used BUSCO v. 5.3.2 (Simão et al. 2015) to identify universal single-copy orthologs in the assemblies against the Hemiptera Odb10 database. In summary, the method uses tools of gene prediction and BLAST searches to identify orthologous sequences in the transcriptomes/genomes from a database of 2,510 hemipteran genes that are known to remain single-copy in evolutionary course. Each ortholog found in a transcriptome/genome is then classified as complete (single or duplicated) or fragmented. Genes that composed the final dataset (i.e., the supermatrix) were those found to be complete and single-copy.

I conducted the phylogenetic analyses using two approaches: a maximum likelihood (ML) with a concatenated dataset of all genes (supermatrix) and a multi-species coalescent analysis based on individual gene trees. The orthology prediction analyses were conducted using the pipeline BUSCO_phylogenomics.py (https://github.com/jamiemcg/BUSCO_phylogenomics). For the supermatrix analysis, the pipeline selects all single-copy complete genes that are found in a given percentage of species, which I chose to be 80% for this study. Then, translated amino acids for each BUSCO protein are aligned separately using Muscle v. 5.2 (Edgar 2004), followed by a curation step with trimAl v. 1.4 (Capella-Gutiérrez et al. 2009) to correct for alignment issues. Finally, all proteins are concatenated into a single supermatrix. Following my criterion of gene sampling, the final dataset comprised 1,915 single-copy orthologs with 915,964 amino acids. Individual coverage levels for each taxon are expressed in Supplementary Table S1.

I chose the best substitution model for the entire protein alignment with ModelFinder (Kalyaanamoorthy et al. 2017) using the AICc criterion. The chosen model was a JTT amino acid model with rate heterogeneity among sites. I inferred the phylogeny with IQ-TREE v. 2.0.7 (Nguyen et al. 2015), which uses a ML approach. Clade support was calculated with 1,000 replicates of ultrafast bootstrap.

The multi-species coalescent analysis was run in ASTRAL-II with default parameters, and clade support levels are expressed as local posterior probabilities. Scripts, matrices and trees resulting for this study are deposited in public repositories (figshare.com/projects/Pentatomidae_Phylogenomics/215380 and github.com/bgenevcius/pentatomidae_phylogenomics).

Results

The phylogenetic analyses of genome-scale data resulted in fully resolved topologies with overall high support (Fig. 1C and D). For both the ML and species-tree coalescent (STC) trees, only 3 clades showed bootstrap below 80 or posterior probabilities below 0.8 (Fig. 1C and D). In both analyses, all families that had more than one representative were recovered as monophyletic, except for Pentatomidae and Tessaratomidae. The first was not monophyletic due to the position of Chalcocoris rutilans, which appeared as the sister group of Plataspididae in the ML analysis and sister to Tessaratomidae + Dinidoridae in the STC analysis. Tessaratomidae was paraphyletic due to the position of Aspongopus chinensis (family Dinidoridae), which comprised a lineage within the Tessaratomidae (Fig. 1C and D). Most other family-level relationships were different between the 2 analyses. For example, the Cydnidae was closely related to Acanthosomatidae in the ML analysis (Fig. 1C), while it was the sister group of Thyreocoridae according to the STC analysis.

The sister group of Pentatomidae, according to the ML analysis, was a clade composed by all remaining families (Fig. 1C). This result was similar in the STC analysis (STC), however with the exclusion of Plataspididae, which appeared as the sister group of the Pentatomidae and all other species of Pentatomoidea. Within Pentatomidae, both analyses were largely congruent. Three major groups (tribes or subfamilies) with more than one representative were monophyletic: Asopinae, Nezarini, and Strachiini (Fig. 1C and D). On the other hand, the tribe Carpocorini was split into a lineage with the neotropical species (Euschistus and Diceraeus) and another with the paleartic genus Dolycoris, which was closely related to Podopini and Aelini. The Neotropical Carpocorini were also recovered as the sister group of the remaining pentatomids. The most notable difference between the 2 methods with respect to Pentatomidae was the position of Erthesina (Halyini). The other differences are related to species or genera, for example, the relationships within Chinavia.

Notably, the relationships among major groups of Pentatomidae were highly supported (Fig. 1C and D). For example, part of the subfamily Pentatominae (the clade Pentatomini + Piezodorini + Cappaeini) was the sister group of the Asopinae; another part of Pentatominae, the Aelini plus Dolycoris, appeared as closely related to the genus Scotinophara, which belongs to the subfamily Podopinae. These results strongly support the subfamily Pentatominae as polyphyletic.

Discussion

Several studies have sought to infer relationships among families of Pentatomoidea in the last two decades, as reviewed by Bianchi et al. (2021). Each study employed different combinations of morphology, nuclear genes, and even complete mitochondrial genomes, and none of them fully agreed with each other. The position of Pentatomidae is illustrative of such incongruity, having been placed as the sister group of Plataspidae (Liu et al. 2019), Lestoniidae + Acanthosomatidae (Wu et al. 2017), Cydnidae + Tessaratomidae + Dinidoridae (Yuan et al. 2015), among others (e.g., Zhao et al. 2018; Roca-Cusachs et al. 2022). In cases as such, the most likely underlying reason is an insufficient amount of data or even the choice of markers per se, which may be a reflection of gene-tree and species-tree discordance (Forthman et al. 2019). My results point to novel possibilities, although there was no congruence between the ML and the STC approaches. Despite some discordance, both analyses indicate that there is no single family or a group of a few families as sister to Pentatomidae. Rather, both suggest a large group of several families as sisters to the stink bugs. The main discordance between both methods is whether the Plataspididae are included in this large group or not (Fig. 1C and D). Perhaps the most similar result to mine was that of Yuan et al. (2015), who recovered the same relationships among families as my ML tree. Interestingly, it seems that the position of Plataspididae has been the most labile among the above-mentioned studies. Considering that I included large amounts of molecular data and still did not resolve the position of this family, it is likely that only the addition of more taxa of Plataspididae and from the other families will allow for a full resolution of the family-level relationships.

While the position of Pentatomidae has been historically controversial, other relationships within the superfamily have shown stronger stability (Bianchi et al. 2021). A notable example is the proximity between Dinidoridae and Tessaratomidae (Yuan et al. 2015; Liu et al. 2019). My genome-scale analyses confirmed this hypothesis, also suggesting that Tessaratomidae is paraphyletic in relation to Dinidoridae. Further studies, including more terminals of these families, are encouraged to check the stability of these relationships. It is also noteworthy that the relationships within Pentatomidae were much more congruent between the ML and STC approaches compared to the relationships among families. This is likely due to the extreme antiquity of these family-level relationships (Johnson et al., 2018). To further address this issue, I suggest designing and implementing ultra-conserved elements approaches, which have shown promise in resolving deep phylogenetic relationships (Faircloth et al. 2015).

The long-lasting absence of a phylogeny for the stink bugs has led authors to consider it a Pandora’s Box (Roca-Cusachs et al. 2022). Recent studies have successfully demonstrated the stink bugs as a non-monophyletic group and that its internal relationships do not match the current classification (Genevcius et al. 2021, Roca-Cusachs et al. 2022). Nevertheless, no study to date has yet provided a satisfactory resolution at deeper nodes of the Pentatomidae tree, failing to confidently indicate which tribes and subfamilies are related to each other. By using a massive dataset of over 1,900 genes, I provide a phylogenomic tree that achieves good resolution and support across virtually all nodes (Fig. 1C). Some of my findings align with the recent Sanger study from Roca-Cusachs et al. (2022), for example, the monophyly of Asopinae, the relationship between Nezarini and Antestiini, the proximity of Rhaphigaster (Pentatomini) and Piezodorus (Piezodorini), and a few others (Fig. 1A–D). In my view, this indicates that part of the clades recovered by Roca-Cusachs et al. (2022) can indeed be trusted, despite their lower support. As expected, some relationships were strikingly different in my study. Notable distinctions were the position of Scotinophara (Podopinae: Podopini) within the Aelini and Carpocorini and the relative position of Asopinae among other clades that are present in both analyses. The strongly supported polyphyly of Pentatominae, the biggest subfamily of Pentatomidae, confirms the need for a revision in the classification of Rider et al. (2018).

Despite the high support levels of my phylogenomic tree, it is noteworthy that a few branches were overly short, especially in the ML analysis (Fig. 1C). I highlight the branch that splits the family nearly in half (one clade being the Asopinae, Halyini, Strachiini, Pentatomini, Cappaeini, and Piezodorini). This indicates that even a massive dataset of almost 2,000 genes may carry modest variation for particular clades. Further, it has been argued that traditional support metrics such as bootstrap and posterior probabilities may be outdated and unrealistic for phylogenomics (Thomson and Brown 2022). As such, I also caution on the use of such metrics as the only evidence of clade reliability. While my study was able to tackle the issue of limited molecular sampling from previous works, the addition of more taxa will be the next challenge for phylogenomic studies of stink bugs, as limited taxon sampling may also affect phylogenetic inference (Sanderson et al. 2010, but see Reddy et al. 2017). Future analyses expanding species sampling should have in mind that high-quality genome-scale data will be necessary, since abundant missing data may perturb phylogenetic inference (Roure et al. 2012).

In conclusion, I demonstrated that incorporating extensive molecular data into phylogenetic matrices may help achieve better supported trees, including the deeper nodes, for stink bugs and their allies. However, taxon sampling remains a limitation, as my study relied solely on publicly available data. To advance and refine these results, future research should focus on increasing the number of species with genome-scale data. In particular, I encourage efforts to sample key taxa with completely unknown positions and unique morphologies, such as Cyrtocorinae and Stirotarsinae (within Pentatomidae) and Serbaninae and Megarididae (within Pentatomoidea). With the acquisition of genomic data getting more accessible and cheaper (Bonetta 2010), another future avenue relates to methodological issues like orthology inference, gene-tree heterogeneity and model selection (Behura 2015). Regardless, I emphasize that my goal here was to build the groundwork for the future systematics of Pentatomoidea, rather than to comprehensively solve lineage relationships and classification issues. As such, I consider my 3 main initial goals to be successfully achieved. First, my results indicate that well supported clades, including deep nodes, are possible with the employment of -omics data. Second, my findings illuminate large incongruities between the phylogenomics and Sanger-based approaches, especially at the levels of families, subfamilies, and tribes. And third, by sampling a substantial portion of families as terminal taxa, I provided a fresh perspective on the position of Pentatomidae within the superfamily.

Disclaimer: Nagoya Protocol: No specimens were included in this project that required permits of any kind.

Acknowledgments

I thank Prof. Fernando P.L. Marques for providing access to computer resources and Dr. Renan Carrenho for comments on an early draft of the manuscript.

Author contributions

Bruno Genevcius (Conceptualization [lead], Data curation [lead], Formal analysis [lead], Investigation [lead], Methodology [lead], Writing—original draft [lead], Writing—review & editing [lead])

Funding

I thank Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) for the post-doctoral fellowship (Grant No. 2018/18184-4).

References