Host and venom evolution in parasitoid wasps: does independently adapting to the same host shape the evolution of the venom gland transcriptome?

Abstract

Background: Venoms have repeatedly evolved over 100 occasions throughout the animal tree of life, making them excellent systems for exploring convergent evolutionary novelty. Growing evidence supports that venom evolution is predominantly driven by prey or host-related selection pressures, and the expression patterns of venom glands reflect adaptive evolution. However, it remains elusive whether the evolution of expression patterns in venom glands is likewise a convergent evolution driven by their prey/host species.

Results: We utilized parasitoid wasps that had independently adapted to Drosophila hosts as models to investigate the convergent evolution of venom gland transcriptomes in 19 hymenopteran species spanning ~ 200 million years of evolution. Comparative transcriptome analysis reveals that the global expression patterns among the venom glands of Drosophila parasitoid wasps do not achieve higher similarity compared to non-Drosophila parasitoid wasps. Further evolutionary analyses of expression patterns at the single gene, orthogroup, and Gene Ontology (GO) term levels indicate that some orthogroups/GO terms show correlation with the Drosophila parasitoid wasps. However, these groups rarely include genes highly expressed in venom glands or putative venom genes in the Drosophila parasitoid wasps.

Conclusions: Our study suggests that convergent evolution may not play a predominant force shaping gene expression levels in the venom gland of the Drosophila parasitoid wasps, offering novel insights into the co-evolution between venom and prey/host.

Keywords: Comparative transcriptome; Convergent evolution; Gene expression; Parasitoid wasps; Venom.

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

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, d_T) 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 d_T 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

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

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