Comparative genomics reveals evolutionary drivers of the dietary shift in Hemiptera

Abstract

Hemiptera insects exhibit a close relationship to plants and demonstrate a diverse range of dietary preferences, encompassing phytophagy as the predominant feeding habit while a minority engages in carnivorous or haematophagous behaviour. To counteract the challenges posed by phytophagous insects, plants have developed an array of toxic compounds, causing significant evolutionary selection pressure on these insects. In this study, we employed a comparative genomics approach to analyse the expansion and contraction of gene families specific to phytophagous insect lineages, along with their adaptive evolutionary traits, utilising representative species from the Hemiptera order. Our investigation revealed substantial expansions of gene families within the phytophagous lineages, especially in the Pentatomomorpha branch represented by Oncopeltus fasciatus and Riptortus pedestris. Notably, these expansions of gene families encoding enzymes are potentially involved in hemipteran-plant interactions. Moreover, the adaptive evolutionary analysis of these lineages revealed a higher prevalence of adaptively evolved genes in the Pentatomomorpha branch. The observed branch-specific gene expansions and adaptive evolution likely contribute significantly to the diversification of species within Hemiptera. These results help enhance our understanding of the genomic characteristics of the evolution of different feeding habits in hemipteran insects.

Introduction

Hemiptera encompasses several crucial agricultural pests such as aphids and whiteflies, as well as significant sanitary pests including Triatoma rubrofasciata, Rhodnius prolixus, and Cimex lectularius (Rosenfeld et al., 2016; Liu et al., 2019; Huang et al., 2021; Rijal et al., 2021). Therefore, the study of Hemiptera holds great significance. There are 97,000 to 103,590 known species representing over 10% of all known insect species (Johnson et al., 2018). The distinguishing characteristic that differentiates Hemiptera from other insects is the variation in mouthparts, which has evolved from chewing mouthparts in their ancestors to piercing-sucking mouthparts to enhance adaptability (Li et al., 2017). These specialised piercing-sucking mouthparts have undergone continuous modification and adaptation, enabling different Hemiptera species to feed on plant phloem, blood, or other fluids (Li et al., 2017).

Hemipterans exhibit a remarkable biodiversity, spanning a wide range of life histories and habitats, encompassing both terrestrial and aquatic ecosystems (Li et al., 2017). The diversity of environments drives their exploration of various food sources such as plants, arthropods, fungi, and vertebrate blood (Rider, 1996). Based on previous phylogenetic study of Hemipteroid, there is a temporal correlation between the evolution of true bugs in Hemipteroids and a shift in feeding habits, implying a possible transition from phytophagy to predation (Sweet, 1979; Johnson et al., 2018). The successful evolution of dietary shifts in Hemiptera insects can be attributed to the interaction between trophic niche and genome-level changes (Weirauch et al., 2019). The environment and prevailing conditions are primary factors that influence genetic variation within insect populations, which leads to the divergence of different species. Among various environmental factors, the trophic niche, which denotes the principal source of nourishment, plays a pivotal role in moulding phenotypic and genomic alterations in insects (Johnson et al., 2018).

Studies investigating dietary variations in Hemiptera revealed that within this insect order, Acyrthosiphon pisum, whitefly, Laodelphax striatellus, Oncopeltus fasciatus, and Riptortus pedestris are classified as phytophagous insects (Li et al., 2017). On the other hand, T. rubrofasciata, R. prolixus, C. lectularius, and Gerris buenoi are carnivorous insects (Rosenfeld et al., 2016; Armisén et al., 2018; Liu et al., 2019; Huang et al., 2021; Rijal et al., 2021). Except G. buenoi, the habitats of these insects are primarily semi-aquatic (Armisén et al., 2018). Plants employ diverse defence strategies against phytophagous insects, they have developed various mechanisms of insect resistance specifically targeting Hemiptera to protect themselves (Santamaria et al., 2013). Consequently, this dynamic interplay exerts selective pressure on phytophagous Hemiptera species.

Detoxification of plant secondary metabolites is a crucial adaptation strategy in the evolutionary dynamics of phytophagous insects, enabling them to neutralise or mitigate the effects of these compounds (Seppey et al., 2019). Several gene families play pivotal roles in this process, including cytochrome P450 monooxygenases, carboxylesterases (CEs), UDP-glycosyltransferases (UGTs), glutathione S-transferases (GSTs), and ATP-binding cassette (ABC) transporters (Voelckel 2014). These genes are instrumental in the structural adjustment, modification, and transport of toxic compounds. In addition to producing compounds that are resistant to insects, plants also produce protein inhibitors targeting specific mechanisms of insect resistance. Notable gene families associated with insect resistance mechanisms include endopeptidases, such as cysteine (CYSs) and serine (SERs) proteases, along with more specialised enzymes like glycoside hydrolases (GHs) (Pauchet et al., 2010; McKenna et al., 2016). Furthermore, chemosensory systems, including the olfactory and gustatory sensors, play a pivotal role in the adaptation of phytophagous insects to their host plants (Goldman-Huertas et al., 2015).

During lineage divergence, insects undergo genomic changes that contribute to adaptive development. Identifying these genes and their associations with phenotypic variations in different insect species is crucial for understanding the process of speciation (Hurst, 2009). Genomic changes such as gene duplication events lead to the expansion of gene families, resulting in the emergence of new gene copies with similar or identical functions (Kondrashov, 2012). Additionally, other genomic changes, such as point mutations in genes, alter their functions. Branch-specific gene expansion refers to the expansion of gene families in specific lineages and encompasses various mechanisms, including but not limited to adaptive evolution (Innan and Kondrashov, 2010). Despite these mechanisms are not exclusively related to adaptation, the production of different gene copies through duplication provides opportunities for natural selection to operate.

Materials and methods

Data sources and quality assessment

This study encompassed the genomes of 13 Hemiptera species carefully selected to represent five suborders of Hemiptera, ensuring a well-balanced sampling across the suborders of Cimicomorpha, Pentatomomorpha, Gerromorpha, Fulgoromorpha, and Sternorrhyncha. Sequence data of R. prolixus, C. lectularius, and A. pisum (CDS and protein sequences) was obtained from Ensembl Metazoa. Additionally, sequence data of O. fasciatus and R. pedestris (CDS and pep sequences) were acquired from InsectBase. Sequence of G. buenoi was obtained from the I5 K-pilot project, while T. rubrofasciata, L. striatellus, Sogatella furcifera, and Trialeurodes vaporariorum were obtained from the I5 K-pilot and GigaDB databases. Sequence data for the whitefly Mediterranean cryptic species (MED) and whitefly Middle East-Asia Minor 1 cryptic species (MEAM1) was obtained from the Whitefly Genome Database (http://www.whiteflygenomics.org/) and sequence data for Diaphorina citric was obtained from NCBI.

For phylogenetic and divergence time analyses, we introduced the fruit fly Drosophila melanogaster as an outgroup. The genomic data of D.melanogaster were obtained from Ensembl Metazoa. Detailed information about the genomic data sources pertinent to this study can be found in the list provided in supplementary Table S1. The quality assessment of all employed genomic datasets was conducted using BUSCO (version 5.2.2) (Waterhouse et al., 2018) with the arthropoda_odb10.2020-09-10 dataset and the model tran.

Phylogenetic analysis of species

We employed the genomes of the 13 Hemiptera species and one outgroup mentioned earlier, resulting in a total of 14 species, to identify direct homologues. To limit redundancy arising from selective splice variants, we retained gene models encoding the longest protein sequence for each locus. Direct homologous gene analysis was conducted using DIAMOND in conjunction with OrthoFinder (version 2.5.4) (Emms and Kelly, 2019). Single-copy genes were identified based on the OrthoFinder results and subsequently utilised for downstream analyses, including the construction of phylogenetic trees, estimation of divergence times, and investigation of gene family contraction and expansion.

For the phylogenetic tree construction, a total of 22,723 protein sequences from the single-copy gene families were utilised. MUSCLE (version 3.8.95) (Edgar, 2004) was used to generate multiple sequence alignments with default parameters for the protein sequences within each single-copy family. The resulting alignments were combined into a super alignment matrix, which was then employed for phylogenetic tree reconstruction using the PROTGAMMAJTT model in the RAxML software (Stamatakis, 2014).

Functional annotation

The protein Pfam family was identified by InterProScan (Jones et al., 2014). The uniref50 database was annotated using Blastp with an e-value threshold set at 1e-20. GO term annotation was performed using eggnog (Suzek et al., 2015). Orthologous groups (OGs) were preliminarily selected based on the annotation in the uniref50 database. Subsequently, candidate OGs were identified by matching the uniref50 annotation and the Pfam and GO term annotation results.

Analysis of contraction and expansion of gene families

Based on the orthogroups (gene families) information derived from the identified orthologous genes using OrthoFinder (Emms and Kelly, 2019), the expansion and contraction of orthologous gene families were assessed using CAFE v4.2.1 (De Bie et al., 2006). This software utilises birth and death processes to model the gene gain and loss events throughout the phylogenetic history process.

Adaptive evolutionary analysis

Furthermore, we conducted adaptive evolutionary selection analyses on the 13 Hemiptera species selected in this study. These analyses aimed to examine the differential occurrence of positive Darwinian selection on specific branches of the phylogeny for individual homologous genes. For this purpose, single-copy orthologous genes generated by OrthoFinder were employed for selective pressure analysis. Multiple sequence comparisons were performed using MACSE (Ranwez et al., 2018), and subsequent analyses utilised the branching site model of the codeml program within the PAMLpackage (https://github.com/Hua-CM/BatchPAML) (Yang, 2007).

To identify positively selected genes, we employed Bonferroni correction and Benjamini-Hochberg false discovery rate (FDR) control to select genes with FDR-corrected p-values below 0.05. We performed comparative sequence recombination checks using the Pairwise Homoplasy Index (PHI) (Bruen et al., 2006), Neighbour Similarity Score (NSS), and Maximum Chi-Square tests implemented in the PhiPack program to exclude false-positive selection events resulting from sequence recombination. Recombination events were considered present when the p-value (q-value) for PHI was less than 0.05 and when at least one other test supported the occurrence of recombination. In this study, we systematically examined the branches of Pentatomomorpha, Fulgoromorpha, Sternorrhyncha, and Cimicomorpha for adaptive evolution using the aforementioned process. Enrichment analysis was conducted to identify genes associated with adaptive evolution, following the procedure described earlier. We utilised the web-based tool agriGO (systemsbiology.cau.edu.cn/agriGOv2) (Tian et al., 2017) and the FlyBase database to assess the enrichment of genes in specific Gene Ontology (GO) terms (Thurmond et al., 2019). Furthermore, the KOBAS (Xie et al., 2011) and BlastKOALA software (Kanehisa et al., 2016) were employed to statistically evaluate the enrichment of genes in the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (Kanehisa and Goto, 2000).

Results

Representative sampling of Hemiptera insects

To investigate the adaptation of Hemiptera to host plants, we conducted a sampling of 13 Hemiptera species, representing various suborders. Specifically, our sampling included five species from the Sternorrhyncha suborder, three species from the Cimicomorpha suborder, two species from the Pentatomomorpha suborder, two species from the Fulgoromorpha suborder, and one species from the Gerromorpha suborder. These species inhabited different habitats, with G. buenoi being semi-aquatic and the remaining species terrestrial (Table S1). The feeding habits of these species varied and consisted of three categories: phytophagous, predaceous, and haematophagous. Notably, within the Cimicomorpha suborder alone, all three feeding habits were represented, among which T. rubrofasciata and R. prolixus were predaceous, C. lectularius being haematophagous. O. fasciatus and R. pedestris were phytophagous (Table S1). Therefore, our sample encompassed different Hemiptera suborders, diverse habitats, and a range of feeding habits, including a substantial number of phytophagous species (Table S1). We assessed the integrity of the genome sequences obtained from the sampled species using BUSCO analysis, which revealed sequence integrity ranging from 74.8% to 99.9% across the sampled Hemiptera species and outgroups. With these sequences, we identified 22,723 directly homologous genes and constructed a time-calibrated species phylogeny using 221 single-copy gene sequences. The sequences from these orthologous groups (OGs) were further subjected to functional annotation, particularly focusing on gene families involved in phytophagous insect-plant interactions. Among these candidate OGs, we identified 72 candidates from eight gene families, including 11 UDP-glycosyl transferases (UGTs), 61 cytochrome P450s (P450s), 13 carboxylesterases (CEs), 1 glutathione S-transferase (GST), 4 serine proteases (SERs), 2 cysteineproteases (CYSs), 21 ATP-binding cassette transporters (ABCs), and 3 glycoside hydrolases (GHs) (Table 1).

Table 1. Candidate gene categories and key terms/identifiers are used to select from the complete annotation sequences using InterProScan

InterProScan Scan-derived categories for Inclusion in OGs are considered for inclusion in candidate (OGs), OGs are required to have at least one sequence matching both UniRef and InterProScan entries. Additionally, supplementary Gene Ontology terms are taken into account.

Phytophagous insects show more frequent gene expansion

According to the CAFE analysis, the predicted λ (gain) value for gene expansions was 0.0020 gain/gene/million years, and the μ (loss) value for gene contractions was 0.0036 loss/gen e/million years (De Bie et al., 2006). Out of the total 22,723 orthologous groups (OGs) examined, the root of the phytophagous insect suborder exhibited a higher frequency of genes. In the root of the Sternorrhyncha suborder, 21 genes underwent expansion, while 5 genes were contracted. Similarly, within the Fulgoromorpha suborder, 64 genes experienced expansion, while 18 genes underwent contraction. In the Cimicomorpha suborder, which includes two predaceous and one haematophagous insect species, all three species are non-phytophagous, and the root exhibited expansion of 3 genes and contraction of 39 genes. Conversely, in the Pentatomomorpha suborder, which comprises two phytophagous insect species, the root showed an expansion of 82 genes and a contraction of 9 genes. The number of gene contractions and expansions in each branch is depicted in fig. 1. Significance analysis using a p-value threshold of <0.01 revealed that 11.11% of the 2525 OGs identified by CAFE exhibited significant contractions or expansions.

Figure 1. The research quantitatively analyses phylogenetics and gene family expansions using genomic datasets. The species tree is based on 221 single-copy orthologous genes, with strong bootstrap support (99% or 100%). Data are from sequenced genomes. (OGs) expand (+) and contract (-) at each subsequence root. The study included herbivorous, semi-aquatic, and blood-feeding insects.

Adaptive evolution is more widespread in the phytophagous branch Pentatomomorpha

We specifically examined the phytophagous insect branches of Pentatomomorpha, Fulgoromorpha, and Sternorrhyncha for adaptive evolution, as well as the non-phytophagous insect branch of Cimicomorpha. The results respectively a total of 18, 3, 1, and 3 genes undergoing adaptive evolution in the branches of Pentatomomorpha, Fulgoromorpha, Sternorrhyncha, and Cimicomorpha (FDR < 0.05) (Table 2). Among these genes, one member of the ABC gene family, ACYPI010103, exhibited adaptive evolution in the Pentatomomorpha branch with FDR value 7.6046E-08 (Table 2). The GO functional enrichment analysis of these adaptively evolved genes revealed that these genes were enriched in 28 different GO terms such as skeletal muscle myosin thick filament assembly, myosin filament assembly, and striated muscle myosin thick filament assembly (Table S2). In terms of KEGG functional enrichment, the adaptively evolved genes were found to be enriched in GTP-binding proteins (Table S3).

Table 2. Result of evolutionary selection pressure analysis on Pentatomomorpha, Fulgoromorpha, Sternorrhyncha, and Cimicomorpha branches

In contrast, no significantly enriched GO or KEGG terms were observed for the adaptively evolved genes in the Fulgoromorpha and Sternorrhyncha branches (P < 0.05). For the non-phytophagous Cimicomorpha branch, the three adaptively evolved genes were found to be enriched in a single GO term: brown fat cell differentiation (Table S2). However, no significantly enriched pathways were detected when analysing KEGG enrichment (Table S3).

Discussion

In this study, we employed comparative genomics to investigate the genomic characteristics, including gene family contraction, expansion and adaptive evolution, across different feeding branches of the Hemiptera suborder. We analysed whitefly MED cryptic species, whitefly MEAM cryptic species, and T. vaporariorum from the polyphagous Sternorrhyncha suborder, as well as A. pisum from the monophagous suborder. In addition, O. fasciatus and R. pedestris, belonging to the Pentatomomorpha suborder, were also monophagous. R. pedestris, commonly known as the bean bug, primarily feeds on leguminous plants, particularly G. max which is commonly referred to as soybean (Huang et al., 2021). On the other hand, O. fasciatus, also known as the large milkweed bug, is a specialised seed feeder and has been observed to in 13 Hemiptera species (Panfilio et al., 2019). Among these, three are predaceous: T. rubrofasciata and R. prolixus from the Cimicomorpha suborder, and G. buenoi from the Gerromorpha suborder (Armisén et al., 2018; Liu et al., 2019; Huang et al., 2021). There is also a haematophagous species, C. lectularius, which belongs to the Cimicomorpha suborder. The remaining nine species were phytophagous, consisting of L. striatellus and S. furcifera from the monophagous Fulgoromorpha suborder, which generally undergo nymphal development on a limited number of host species within the Asclepias genus. To test the hypothesis that gene expansion and adaptive evolution contribute to dietary changes in Hemiptera, we analysed gene family expansions and adaptive evolutionary events in the 13 Hemiptera species. It is widely recognised that phytophagous insects require specific genes for the detoxification of plant secondary metabolites and digestion of plant tissues on their host plants, to achieve the purpose of adapting to different environments comparing to predaceous and haematophagous insect (Heidel-Fischer and Vogel, 2015; Simon et al., 2015). Both gene family expansion and adaptive evolution are important molecular mechanisms underlying organismal adaptation to varying environments (Innan and Kondrashov, 2010; Kondrashov, 2012). Previous studies on beetles have demonstrated the significance of gene family expansion and adaptive evolution in the dietary shift to phytophagy (Seppey et al., 2019). In our study, we observed gene family expansion in the phytophagous insect branch. However, regarding adaptive evolution, a consistent association between adaptive evolution and gene family expansion was only detected in the Pentatomomorpha suborder. Therefore, adaptive evolution may only play a minor role in the intrinsic mechanism of gene family expansion in Hemiptera. Our analysis identified the phytophagous trophic niche as a driving force behind gene family expansion in Hemiptera. However, the expansion of these genes did not lead to functional divergence (Innan and Kondrashov, 2010; Kondrashov, 2012).

There is a shared anatomical feature of a specialised piercing and sucking mouthpart within the Hemiptera insects. However, these insects have undergone diversification to exploit a wide range of food sources, including seeds and plant tissues (phytophagy), and even vertebrate blood (haematophagy) (Li et al., 2017). As a result of this diversification, numerous hemipterans have become significant agricultural pests or vectors of human diseases. Consequently, extensive efforts have been devoted to genome sequencing in these species, aiming to unravel their genomic composition and elucidate their complex biological characteristics (Li et al., 2017; Johnson et al., 2018).

When conducting comparative analyses of datasets from different species, it is crucial to ensure that all analysed species have comparable gene contents. To mitigate the risk of false-positive gene identification, we implemented an additional filtering step based on the method proposed by Seppey et al. (2019). This process combined InterProScan results using Pfam or InterPro identifiers, along with Gene Ontology information. Additionally, we integrated the relevant keywords from the carefully selected UniRef dataset for further filtration (Table 1). Although this stringent filtering approach might have potentially excluded some candidate orthogroups (OGs) from our analysis, it was crucial in ensuring the accuracy and validity of the OGs used in this study.

In addition to the gene expansion in Pentatomomorpha, the positive selection analyses also identified a higher number of genes in the branch that underwent positive selection for adaptive evolution. A noteworthy gene in this context is the ABC transporter, which is known to be involved in the detoxification of plant secondary metabolites (Yazaki, 2006). This finding suggests that selective pressure has played a role in the expansione expansion of the gene family associated with detoxification enzymes (Yazaki, 2006; Calla, 2021).

Emerging research indicates that ABC transporters are also involved in insect detoxification of pesticides, Bt and plant secondary metabolites (Yazaki, 2006; Heckel, 2012; Merzendorfer, 2014). In addition to the ABC transporter, several other adaptively evolved genes in the Pentatomomorpha branch are associated with muscle development, including skeletal muscle myosin thick filament assembly, myosin filament assembly, and muscle cell development, and others (Table 3). The adaptive evolution of these genes provides a mechanistic foundation for the enhanced locomotor capacity of Pentatomomorpha. Several mechanisms contribute to gene family expansion, with adaptive evolution being one of the driving factors, but it may not be the predominant factor (Innan and Kondrashov, 2010; Seppey et al., 2019). Although a significant number of gene family expansions and adaptive evolutions were observed in the Pentatomomorpha branch, we did not find the same level of cooperativity of gene family expansions and adaptive evolutions in Hemiptera as observed in Coleoptera (Seppey et al., 2019). Therefore, for gene family expansion, neutral or purifying selection forces may play a more prominent role.

Table 3. Candidate OGs with CAFE p-value <0.01 were identified

OGs were categorised as Cytochrome P450 (P450), Carboxylesterase (CE), Glutathione S-Transferase (GST), and Cysteine Protease (CY5) based on OrthoDB v8 functional annotations.

Neutral or purifying selection is likely the primary mode of occurrence for gene family expansions in Hemiptera. It can be thought that genes associated with nutrition and detoxification might result from gene duplication rather than the emergence of novel functions (Panfilio et al., 2019; Volonté et al., 2022).

Conclusions

In this study, we analysed branch-specific gene family expansions and adaptive evolutions by comparing genomic data from representative species of Hemiptera. Our study revealed a higher occurrence rate of gene family expansions in the phytophagous branch, particularly within the Pentatomomorpha branch. Furthermore, through further analysis of adaptive evolutions within each branch, we identified a greater abundance of genes that underwent adaptive evolutions in the Pentatomomorpha branch.