To test if independent adaptation to the same or similar hosts could lead to convergent gene expression patterns in the venom gland of parasitoid wasps from different lineages, we investigated the venom gland transcriptomes of 19 hymenopteran species spanning six superfamilies, including five Drosophila parasitoid wasps that represent four independent adaptations to the Drosophila hosts (Additional file 1: Table S1). These five Drosophila parasitoid wasps include three larval Drosophila parasitoids: Leptopilina boulardi (Figitidae), L. heterotoma (Figitidae) and Asobara japonica (Braconidae), and two pupal Drosophila parasitoids: Trichopria drosophilae (Diapriidae) and Pachycrepoideus vindemmiae (Pteromalidae). Notably, high-quality genomes are available for all species except T. drosophilae. We first performed orthology inference analysis using the protein sequences of the 19 species, which allowed us to identify 1479 single-copy orthologous genes shared by all species. Subsequently, we calculated the expression level (TPM) of each gene in each sample and finally generated a normalized expression matrix of the venom gland (TPM10K, see “Methods” for a detailed description) to ensure comparability across species. If independent adaptation to the same host can result in convergent gene expression patterns in the venom gland of parasitoid wasps, it is reasonable to expect that the transcriptome of the venom gland in Drosophila parasitoid wasps will exhibit similarities. Specifically, we anticipate that the gene expression tree of the venom gland will deviate from the species phylogenetic tree, leading to the clustering of Drosophila parasitoid wasps together [2, 21].

Principal component analysis fails to cluster Drosophila parasitoid wasps

To explore the overall gene expression patterns of venom glands across parasitoid wasps, we first performed principal component analysis (PCA) on both genome sequence data and expression data based on the 1479 one-to-one orthologues genes. We observed a different pattern between the two PCA plots. The 19 hymenopteran species were strictly separated according to taxon for the genome sequence data (Fig. 1A). However, the transcriptome data of venom glands was contrary to our initial expectation, the Drosophila parasitoids data failed to cluster together, nor exhibited clear separation from other species (Fig. 1B, Additional file 2: Fig. S1). We noticed that T. drosophilae was separated from other species, which could be an artifact of the lack of reference genome. We further examined the correlation between principal components (PCs) with trait data, including host species and taxonomy (Fig. 1C). Our analysis revealed that PC2 possessed a significant positive correlation with lineages, indicating the dominant role of lineage differences in driving venom gland gene expressions. However, we did not observe any PCs correlated with host species significantly (Fig. 1C). To understand if higher similarity exists within Drosophila parasitoids, compared to parasitoid wasps of other host species on the venom gland transcriptomes, we computed the Spearman rank correlation coefficient ρ. As a result, there was no more resemblance observed in venom glands between the Drosophila parasitoids and non-Drosophila parasitoids (Fig. 1D).

Fig. 1

Principal component analysis (PCA) of venom gland transcriptomes from 19 Hymenoptera species. A–B PCA on genomic sequences and on expression data based on the 1479 one-to-one single copy genes. See Additional file 3 for individual data values. C Heatmap showing the Pearson’s correlation coefficient of the PC1-5 eigenvectors from venom gland gene expression with the phylogeny and host species of parasitoid wasps. P values by two-sided Pearson’s correlation test (**P < 0.001). D Comparison of the similarity of venom glands between Drosophila parasitoids (n = 4) and non-Drosophila parasitoids (n = 14). The Spearman rank correlation coefficient ρ was used to scale the similarity of the venom gland transcriptomes, and no significant difference was observed between Drosophila parasitoids and non-Drosophila parasitoids (two-sided Wilcoxon rank sum test, P > 0.05). See Additional file 3 for individual data values

Gene expression phylogeny of venom glands is similar with the species phylogeny

To further dissect the trajectory of venom gland transcriptome evolution, we reconstructed the expression tree using the neighbor-joining (NJ) method based on the Spearman coefficient distance matrix (Additional file 2: Fig. S2). First, we compared the topology of the expression tree with the species phylogenetic tree constructed by the maximum likelihood method based on the protein sequences of the 1479 single-copy orthologous genes (Fig. 2A). We found consistency between the expression tree and the species tree (Fig. 2B). The expression tree correctly clustered the parasitoid wasps from the same superfamily, except for a few outlier species with special traits. Notably, it accurately resolved most of the relationships in Chalcidoidea except for the fig wasp C. solmsi, which shifts to adapting on plant hosts. In particular, the five Drosophila parasitoids were placed near their close relatives instead of being clustered together themselves. This suggests that the alternations in venom gland gene expression accumulate over evolutionary time, and closely related species tend to exhibit similar venom gland expression atlases.

Fig. 2

Comparative transcriptome analysis of venom glands from 19 Hymenoptera species. A Species phylogeny of 19 Hymenoptera species. The maximum-likelihood tree was reconstructed by IQ-TREE based on the multisequence alignments of 1479 one-to-one orthologous genes, and branches in red represent four independent origins of adapting to Drosophila hosts. All nodes received 100% bootstrap support. B Topology comparison between the species tree and the venom gland expression tree. Left is the species phylogeny obtained from A, and right is the expression tree reconstructed by the NJ method based on the Spearman distance matrix of the venom gland. C Frequency distributions of topology distance (PH85 distance, dT) between the genome tree and the random tree, transcriptome tree, transcriptome tree with the 50% most highly expressed genes, transcriptome tree based on OG-level expression matrix, and transcriptome tree based on GO-level expression matrix. Each distribution was obtained based on 10,000 random trees or bootstrapped trees (n = 10,000). Dotted line represented the observed dT between genome tree and other trees based on the original data. See Additional file 3 for individual data values. D–E Pairwise comparison of the distance matrix between the species tree and the expression tree based on gene-level expression matrix (n = 171) (D) and based on the OG-level expression matrix (n = 171) (E). The two-sided Mantel test was used to test for significant correlation between the distance matrices. See Additional file 3 for individual data values

To quantify the differences in topology between the expression tree and species phylogenetic tree, we employed the PH85 distance (dT), which is twice the number of internal branches defining different bipartitions of the tips. The calculated distance between the expression tree topology and the species tree topology was determined to be 20 (Fig. 2C). For comparison, we utilized the bootstrap method with 10,000 replications to compare dT between the transcriptome tree and the species tree with that between the transcriptome tree and the randomly generated tree. We found the dT between the expression tree and the species tree was significantly smaller than that between the random tree and the species tree (average dT = 30, P < 0.001), and no significant change when using the 50% most highly expressed genes (Additional file 2: Fig. S2). Additionally, none of the obtained topologies of random tree exhibited a dT value below 20 (Fig. 2C). These results indicated that the gene expression phylogeny of venom glands closely resembled the species phylogeny and this observation was unlikely to be attributed to chance.

In addition to comparing the topologies, we assessed the correlation between the branch lengths of the expression tree and the species tree. Our results revealed a significant positive correlation between the expression tree of venom glands and the corresponding species tree (R = 0.364, P = 0.013, Mantel test; Fig. 2D). This finding further indicates that the overall divergence of gene expression in the venom glands of parasitoid wasps increases with evolutionary time.

Considering that we might miss information from multicopy gene families due to gene loss or gene duplication (i.e., one-to-many and many-to-many orthologs), we also examined another dataset containing 3614 orthogroups (OGs) of orthologous genes shared by all species. Given that OG contains a group of genes derived from a common origin, we used a sum of TPM values in each OG to represent expression abundances at the OG level and generated a normalized (TPM10K) expression matrix to repeat the expression evolutionary analysis. We observed a reduction in the dT value from the species tree to 18 when using the OG-level expression tree (Fig. 2C, Additional file 2: Fig. S2), and there was a stronger positive correlation of the branch lengths between the species trees and the expression tree (Fig. 2D). This is unsurprising because the expression level of orthogroups is composed of many different genes with variable expression profiles and hence distinguish changes in gene expression patterns less sensitively.

We next ask whether the expression pattern of the functional group in the venom glands of Drosophila parasitoid wasps shows convergent shifts. We assigned each gene to specific gene ontology (GO) terms and identified 16,862 GO groups shared by all 19 species. To quantify the expression level of each GO term, we summed the TPM10K values of all genes assigned to that particular GO term. This summation allowed us to generate a GO-level expression matrix. Our subsequent analysis using the GO-level expression data yielded a calculated dT from the species tree of 24, which is larger compared to both the gene-level and OG-level expression data analyses (Fig. 2C, Additional file 2: Fig. S2). Unexpectedly, we observed that the three larval stage Drosophila parasitoid wasps (L. boulardi, L. heterotoma, and A. japonica) formed a monophyletic group in the GO-level expression tree (Additional file 2: Fig. S2). This grouping was supported by a high bootstrap value of 92, indicating strong statistical support for the clustering of these species.

In summary, our expression evolutionary analysis revealed that the expression tree of venom gland was generally consistent with the species phylogeny. This suggests that the convergent adaptation to Drosophila may not contribute to the evolution of venom gland expression landscape of parasitoid wasps.

Phylogenetic analysis based on the expression of putative venom genes does not support clustering in Drosophila parasitoid wasps

Parasitoid wasps typically possess a relatively modest array of venom proteins, with reported numbers ranging from 59 to 210 across ten species analyzed in this study [17, 19, 20, 22,23,24,25,26,27,28]. Our above analysis, which included 1479 single-copy genes, 3614 orthogroups, and 16,862 GO terms likely exceeds the number of venom genes. This potential overrepresentation could obscure the detection of convergent signals in the venom gland transcriptomes of Drosophila parasitoid wasps. To address this issue, we first identified the highly expressed genes in the venom glands of each species. Orthogroups/GO terms containing these genes were assigned as venom highly expressed gene-containing groups (VHEOGs/VHEGOs). Among the 1479 single-copy orthogroups, 3614 orthogroups, and 16,862 GO terms shared by the 19 species, 560 single-copy orthogroups, 1202 shared orthogroups, and 13,880 shared GO terms were classified as VHEOGs/VHEGOs, respectively. We then performed phylogenetic analyses of the expression patterns focusing solely on these VHEOGs/VHEGOs. The results were consistent with previous findings. The five Drosophila parasitoid wasps did not cluster in the expression trees based on either dataset. The three larval Drosophila parasitoid wasps still formed a monophyletic group, but the bootstrap value was reduced to 48 when using the expression matrix based on 13,880 VHEGOs (Additional file 2: Fig. S3).

We further identified putative venom genes in each species and assigned relevant orthogroups and GO terms as venom gene-containing groups (VCOGs/VCGOs). Specifically, there were 106 single-copy orthogroups, 284 shared orthogroups, and 7563 shared GO terms classified as VCOGs/VCGOs, respectively. Similar phylogenetic analyses were carried out using the expression matrices based on these VCOGs/VCGOs. Again, the expression trees failed to cluster the five Drosophila parasitoid wasps from either dataset. Notably, we did not observe the clustering of the three larval Drosophila parasitoids when using the expression matrix based on the 7563 VCGOs (Additional file 2: Fig. S4). We then filtered our dataset to include only highly expressed genes in the venom gland and putative venom genes of the five Drosophila parasitoid wasps. This refinement resulted in 134 single-copy VHEOGs, 326 shared VHEOGs, 8235 shared VHEGOs, 26 single-copy VCOGs, 74 shared VCOGs, and 3914 shared VCGOs, respectively Subsequent analysis of the expression trees based on these datasets yielded similar results. We did not observe the clustering of the five Drosophila parasitoid wasps or the three larval Drosophila parasitoid wasps in either dataset (Additional file 2: Fig. S5, Additional file 2: Fig. S6).

Limited evidence for convergent evolution in the expression profiles of individual orthogroups/GO terms among Drosophila parasitoid wasps

The results that lack a strong convergent signal of global expression patterns in the venom gland of Drosophila parasitoid wasps caused us to examine if individual genes, orthogroups, and functional groups showed similar patterns which may correspond to the convergent adaptation to Drosophila hosts. To this end, we first constructed the expression tree of each 1479 one-to-one single copy orthologues genes by the NJ method based on the Euclidian distance matrix, respectively. We next screened each of the topologies and tended to identify expression trees in which the five Drosophila parasitoid wasps formed a monophyletic group, which was particularly effective in detecting adaptive evolution in gene expression data [2, 21]. We did not observe any genes that meet this criterion. Only one gene (Frizzled-2) where the expression tree showed a monophyletic group of four Drosophila parasitoid wasps (L. boulardi, L. heterotoma, A. japonica, and T. drosophilae) and three genes (ADP-ribose pyrophosphatase, nuclear pore membrane glycoprotein, and aconitate hydratase) where the expression tree showed a monophyletic group of three Drosophila parasitoid wasps (L. boulardi, L. heterotoma, and A. japonica/P. vindemmiae) (Fig. 3A, Additional file 1: Table S2). This observation falls significantly below the expected values (P = 0.00016, binomial test), indicating that there is a lack of genes displaying a similar expression pattern in the venom gland of Drosophila parasitoid wasps. We next performed the same analysis on each of the 3124 orthogroups shared by the 19 species. We identified one additional orthogroup (OG0004615: serine/threonine protein kinase) where the expression tree clustered the four Drosophila parasitoid wasps (L. boulardi, L. heterotoma, A. japonica, and T. drosophilae) and nine orthogroups where the expression trees exhibited a monophyletic grouping of three Drosophila parasitoid wasps (Additional file 1: Table S3), significantly below the expectation values (P = 6.5e-07, binomial test). Collectively, these findings suggest that, consistent with the overall expression patterns in the venom gland of Drosophila parasitoid wasps, convergent signal also could not be detected at both single-gene and orthogroup levels.

Fig. 3

Representative cases of expression trees showing monophyletic grouping in Drosophila parasitoid wasps. A Genes exhibiting a monophyletic grouping of four or three Drosophila parasitoid wasps in the expression tree. B GO terms exhibiting a monophyletic grouping of four Drosophila parasitoid wasps in the expression tree. The red branches in the tree indicate clustering of Drosophila parasitoid wasps. The bar plot illustrates the scaled expression levels of each gene or GO term across various species. These genes or GO terms demonstrate differential expression, either higher or lower, within the Drosophila monophyletic group compared to other species (Red). See Additional file 3 for individual data values

For the GO-level analysis, we extended the same topology examination to all the 16,862 GO terms shared by the 19 species. The results were consistent with the individual gene level analysis. We did not find any GO terms where the expression tree showed a monophyletic group of all five Drosophila parasitoid wasps, and only four GO terms (pH elevation, intracellular pH elevation, positive regulation of amyloid-beta formation, and positive regulation of amyloid precursor protein catabolic process) where the expression tree showed a monophyletic group of the remaining four Drosophila parasitoid wasps, excluding T. drosophilae (Fig. 3B, Additional file 1: Table S4). Interestingly, we found 74 GO terms where the expression tree showed a monophyletic group of three of the five Drosophila parasitoid wasps, in which the clustering of three larval stage Drosophila parasitoid wasps (L. boulardi, L. heterotoma, and A. japonica) was significantly higher than the expectation of 7.4 (16 out of 74, P = 0.002, binomial test, Additional file 1: Table S4). To determine if the observation of 16 GO terms which clustering the three larval stage Drosophila parasitoid wasps was caused by chance, we enumerated all the 969 possible combinations of selecting three parasitoid wasps from the 19 hymenopteran species and counted the number of GO terms that supported each combination. We ranked each combinations based on the total number of supported GO terms and found that 698 out of 969 combinations (72%) had a supported GO number ≦ 16 (Additional file 2: Fig. S7). This suggested that the three larval stage Drosophila parasitoid wasps exhibited a higher degree of clustering in the GO-level expression tree compared to other species. We next meticulously examined the expression patterns of each of the 16 GO terms and found that the three larval stage Drosophila parasitoid wasps showed either higher (14 GO terms) or lower (2 GO terms) expressions than other species (Fig. 3B). These GO categories primarily encompassed transporter, central nervous system myelination, and peptide hormone processing (Additional file 2: Fig. S8). Although we found a small number of GO terms where the expression tree showed a monophyletic grouping of the three larval stage parasitoid wasps, the majority of GO terms did not exhibit convergent signals in the Drosophila parasitoid wasps. Further analysis revealed that 4 of the 16 GO terms (GO:0016486: peptide hormone processing, GO:0022010: central nervous system myelination, GO:0032286: central nervous system myelin maintenance, and GO:0080132: fatty acid alpha-hydroxylase activity) contained genes previously identified as highly expressed in the venom gland of the three larval Drosophila parasitoid wasps. Only one GO term (GO:0016486: peptide hormone processing) contained a putative venom gene (neprilysin) identified in A. japonica (Additional file 1: Table S5). These findings suggest that the observed convergence in the expression of these 16 GO terms in the three larval Drosophila parasitoid wasps may be attributed to other conserved genes rather than actual venom compositions.

In addition to directly comparing the topology of the expression tree to the species tree for each orthogroup/GO term, we used Bayesian phylogenetic generalized linear mixed models (MCMCglmm models) to identify orthogroups/GO terms whose expression levels correlate with host species. Phylogenetic relationships and different genes within the same orthogroups or GO terms were treated as random effects in each model. Among the 1479 single copy genes, we identified 19 genes whose expression levels showed significant correlation with Drosophila parasitoid wasps (Additional file 1: Table S6). One of these genes (Frizzled-2) was previous identified as forming a monophyletic grouping of four Drosophila parasitoid wasps. However, only two genes (Pvin00505 and Pvin06853, see details in Additional file 1: Table S6) were found to be highly expressed in the venom gland of P. vindemmiae, and none of these genes were identified as putative venom genes across the five Drosophila parasitoid wasps. For the 3124 shared orthogroups, we identified an additional 52 orthogroups where the expression patterns significantly correlated with Drosophila parasitoid wasps (Additional file 1: Table S7). Among these, 14 orthogroups contained genes previous identified as highly expressed genes in the venom gland of Drosophila parasitoid wasps. Five orthogroups contained genes previously defined as putative venom genes in the Drosophila parasitoid wasps. However, none of these orthogroups contained genes highly expressed in the venom gland or putative venom genes across the five Drosophila parasitoid wasps. In the GO-level analysis, we identified 65 GO terms with expression levels significantly correlated with Drosophila parasitoid wasps (Additional file 1: Table S8). Among these, 23 GO terms contained genes previously identified as highly expressed in the venom gland of Drosophila parasitoid wasps, with 12 GO terms harboring putative venom genes in Drosophila parasitoid wasps. Nevertheless, only two GO terms (GO:0005471: ATP: ADP antiporter activity; GO:1990737: response to manganese-induced endoplasmic reticulum stress) contained highly expressed genes in the venom gland across all five Drosophila parasitoid wasps, and none of these GO terms contained putative venom genes across all five Drosophila parasitoid wasps.

Furthermore, we ran an additional MCMCglmm model for each GO term to detect convergent signals among the three larval Drosophila parasitoid wasps. We identified 37 GO terms whose expression levels were significantly associated with the three larval Drosophila parasitoid wasps, with six of these GO terms overlapping with the 16 GO terms identified using the topology comparison method (Additional file 1: Table S9, Additional file 1: Table S10). We observed that nine GO terms contained genes previously identified as highly expressed in the venom gland across the three larval Drosophila parasitoid wasps. Four GO terms contained putative venom genes specific to one or two larval Drosophila parasitoid wasps, but not across all three larval Drosophila parasitoid wasps.

Altogether, our analyses of expression levels, employing both topology comparison and MCMCglmm models across individual genes, orthogroups, and GO terms, reveal that a small number of orthogroups or GO terms exhibit similar expression patterns in the five Drosophila parasitoid wasps or the three larval Drosophila parasitoid wasps. However, most of these identified orthogroups or GO terms do not contain highly expressed genes or putative venom genes across the five Drosophila parasitoid wasps or the three larval Drosophila parasitoid wasps. This suggests that convergent evolution is not a major force shaping gene expression levels in the venom gland of Drosophila parasitoid wasps.

Stabilizing selection is the primary evolutionary force driving gene expression divergence in the venom gland of parasitoid wasps

It is widely accepted that venom is an adaptive trait, and the diversity in venom composition and function is subject to strong natural selection as a result of adaptation towards specific diets [29, 30]. However, our transcriptome evolutionary analyses did not detect convergent signals in the evolution of gene expression in the venom gland of parasitoid wasps that have independently adapting to Drosophila hosts. In contrast, our findings indicated a phylogenetic signal in the evolution of global gene expression in the venom gland of parasitoid wasps. This has prompted us to investigate the underlying evolutionary forces that driving gene expression variation in the venom gland of parasitoid wasps. More specifically, we aim to identify if expression variation in individual genes in the venom gland of parasitoid wasps is evolving through genetic drift, stabilizing selection or directional selection. To achieve this, we tested, for each 1479 one-to-one single copy gene separately, the model fit to three common models of trait evolution along the time-calibrated species tree including Brownian motion (BM), a single-optimum Ornstein–Uhlenbeck (OU) process or an OU model with branch shift (OUM). The BM model assumes that the expression variation between species accumulated over time due to divergence and evolutionary rate, often likened to genetic drift in evolutionary modeling. OU model introduces a selective regime that pulls expression values towards an optimum, reflecting stabilizing selection. The OU model can be expanded to incorporate branch-specific events (OUM), resembling directional selection in specific lineages [31]. We found that 1115 out of the 1479 genes (75.3%) exhibited a better fit to the OU model compared to the other two models (Fig. 4A). This suggests that the majority of genes have evolved in accordance with the OU model, providing evidence that stabilizing selection has played a primary role in shaping the gene expression patterns in the venom gland of parasitoid wasps. However, we found that only 65 genes were identified to be putative venom protein-encoding genes, indicating the observation is likely driven by most genes not encoding venoms. Additionally, we conducted further analysis of our expression data by comparing the expected gene expression divergence under an OU model to the observed expression data. We employed a simulation approach where we generated expression levels for 10,000 genes along the known time-calibrated species phylogeny with a range of different α values under an OU model, respectively. Subsequently, we computed the mean squared expression distance for each pairwise species based on the simulated expression levels and compared them to the observed data. We found that our observed gene expression data had a close fit to the simulated data (Fig. 4B), further suggesting that the stabilizing selection is the primary evolutionary force driving gene expression divergence in the venom gland of parasitoid wasps.

Fig. 4

Gene expression dynamics in the venom gland across hymenopteran species. A The number of per-gene expression patterns that have evolved under a BM, OU, or OUM model of trait evolution. See Additional file 3 for individual data values. B Pairwise mean squared distance between P. puparum and other species across evolutionary time. Red: observed gene expression divergences. Blue: expression divergence from 10,000 simulated trajectories under OU model with different values of alpha (n = 10,000). See Additional file 3 for individual data values