Genomic and transcriptomic analyses support a silk gland origin of spider venom glands

Abstract

Background: Spiders comprise a hyperdiverse lineage of predators with venom systems, yet the origin of functionally novel spider venom glands remains unclear. Previous studies have hypothesized that spider venom glands originated from salivary glands or evolved from silk-producing glands present in early chelicerates. However, there is insufficient molecular evidence to indicate similarity among them. Here, we provide comparative analyses of genome and transcriptome data from various lineages of spiders and other arthropods to advance our understanding of spider venom gland evolution.

Results: We generated a chromosome-level genome assembly of a model spider species, the common house spider (Parasteatoda tepidariorum). Module preservation, GO semantic similarity, and differentially upregulated gene similarity analyses demonstrated a lower similarity in gene expressions between the venom glands and salivary glands compared to the silk glands, which questions the validity of the salivary gland origin hypothesis but unexpectedly prefers to support the ancestral silk gland origin hypothesis. The conserved core network in the venom and silk glands was mainly correlated with transcription regulation, protein modification, transport, and signal transduction pathways. At the genetic level, we found that many genes in the venom gland-specific transcription modules show positive selection and upregulated expressions, suggesting that genetic variation plays an important role in the evolution of venom glands.

Conclusions: This research implies the unique origin and evolutionary path of spider venom glands and provides a basis for understanding the diverse molecular characteristics of venom systems.

Keywords: Adaptive traits; Comparative transcriptomics; Gene co-expression networks; Gene selection pressure; Genomics.

Fig. 1

Core network module of the common house spider venom glands. a Cytoscape plot of the core network module. We termed the largest module that contained 2765 genes in the venom glands the core network, which has great functional relevance to the venom glands. Of these, 17 hub genes with the highest degrees of connection and their connections were visualized. Line thicknesses indicate the interaction strengths, and circle sizes represent the connection degrees. b GO enrichment of genes in the core network module. The complete enrichment results are shown in Additional file 1: Table S6

Fig. 2

The expression patterns between the spider venom glands and other tissues across multiple species. a PCA using the expression levels of 1983 orthologs from ten species. The shapes and colors of the points represent species and tissues, respectively. The red dotted circles represent the clustering of glandular tissues. b PCA of the 2388 ortholog expression matrix across six species in the dotted box. c PCA of the 3952 ortholog expression matrix across four species, whose branches are shown in light blue. d Module preservation between the common house spider venom glands and other tissues. Z_summary > 10 implies strong preservation; Z_summary values between 2 and 10 indicate weak to moderate evidence of preservation; if Z_summary < 2, there is no evidence that the module is preserved. The red dot (module 1) represents the core network. e The observed pairwise semantic similarity (SS) scores and permutated ones between the common house spider venom glands and other tissues. Of these, the fifth density plot represents the high similarity between the spider ovary and fruit fly ovary (this value was used as a control in our analysis). The vertical lines show the observed pairwise SS values. The shades show 1000 permutated SS values with 95th and 90th percentiles labeled. f Workflow for similarity index (SI) calculations among the DUGs of glandular tissues. See the “Methods” section for the meaning of the equation. SI, similarity index; DUGs, differentially upregulated genes. g SI comparisons among the DUGs of the common house spider venom glands and glandular tissues. The comparison result between the spider ovary and fruit fly ovary is labeled in blue, showing high similarity (this value was also introduced as a control). The dotted line in the violin plot represents the median value. **FDR < 0.01 (Mann‒Whitney U test). See Additional file 2: Fig. S7–S10 for other comparison results

Fig. 3

Expression differentiation of gene families in the common house spider venom glands and silk glands. Our results clarified that expression divergences of gene duplication events (including TFs) are common in venom and silk glands. a Venn diagram of specifically expressed orthogroups in spider venom and silk glands. The orthogroup in which a specifically expressed gene is located is considered to be specifically expressed. A total of 358 and 862 gene orthogroups were specifically expressed in the venom and silk glands, respectively, with 79 orthogroups shared between both tissue types. b TF orthogroups shared between the venom and silk glands showing expansions in spiders. P-values indicate significant differences in gene numbers (Mann‒Whitney U test). c Relative expression abundances of specifically expressed toxin paralogs in the venom and silk glands. d Nine TFs showing high module membership with the toxin gene CRISP-3. These genes were specifically expressed in spider venom glands. Line thicknesses indicate the interaction strengths. Six genes labeled red represent the upregulation in the venom glands. e Six upregulated TFs associated with the toxin gene CRISP-3. Boxes show the range, upper and lower quartiles, and median. HOX domain, TFs containing the HOX domain; ASH1, achaete-scute homolog 1; Zfp65, zinc finger 65; AITX, kappaPI-actitoxin-Avd3c; PLA2, phospholipase A2; SjAPI, venom peptide SjAPI; CRISPs, cysteine-rich secretory proteins; CRISP-3, cysteine-rich secretory protein 3; Lbx1, transcription factor LBX1; NKX2-6, homeobox protein Nkx-2.6; PHOX2, paired mesoderm homeobox protein 2; pok, ets DNA-binding protein pokkuri; SOBP, sine oculis-binding protein homolog; Sox4, SRY-Box transcription factor 4; Thrb, thyroid hormone receptor beta-A; XlCGF17.1, gastrula zinc finger protein XlCGF17.1; Zasp, PDZ and LIM domain protein Zasp

Fig. 4

Gene variations in specific transcription modules of spider venom glands. a Heatmap of Spearman correlation coefficients between specific modules of glandular tissues. Venom gland-specific modules contain 166 orthologs found through the isa2 method. venom, venom gland; silk, silk gland; saliva, salivary gland; others, other glands, including the prothoracic, hypopharyngeal, and lymph glands of insects. b REVIGO clusters of the significantly enriched GO terms for venom gland-specific modules. Bubble sizes indicate the number of GO terms in each cluster, and the colors represent the corrected enrichment P-value on a log₁₀ scale. Similar clusters plot closer to each other. c Differential expression analysis shows up- and downregulated genes in spider venom glands. The venom glands were compared with spider silk glands (green box and circle) and all salivary glands (orange box and circle). False discovery rates (FDR) ≤ 0.01 are indicated by blue dots, while FDRs > 0.01 are indicated by black dots. Red dots signify the positively selected genes at the ancestor branch of spiders, as well as the genes in specific transcription modules. Digits in the Venn diagrams represent differentially expressed gene numbers in each group. Ptp36E, protein tyrosine phosphatase 36E; msps, mini spindles; DUGs, differentially upregulated genes; DDGs, differentially downregulated genes. d Time tree of ten species across arthropods. Yellow branches indicate the spiders used in our study; the red dot signifies the spider ancestor. MYA, million years ago. e GO enrichment of fourteen positively selected genes at the ancestral branch of spiders. These genes were contained in the venom gland-specific transcription modules and were differentially upregulated in the venom glands. Digits in the circles indicate the gene numbers enriched in the terms. f Specific mutation of the gene Ptp36E. Two sites (orange) of this gene were positively selected at the spider ancestor branch. This gene was under relaxed selection at three spider branches

Fig. 5

Hypotheses regarding the evolutionary derivation of spider venom glands. a Salivary gland origin hypothesis. We challenge this hypothesis based on the comparison results of module preservation, GO semantic similarity, and DUG similarity analyses (see Fig. 2). b Ancestral silk gland origin hypothesis. Our analyses prefer to support the concept that venom glands are likely derived from silk-producing glands present in early chelicerates. Previous assumptions were that spider silk glands evolved from accessory glands or were derived from the coxal glands [37, 38]. Modern spider silk glands may generate functional convergence with ancestral silk-producing glands, in turn resulting in high transcriptional similarities between spider venom glands and silk glands