Sexually dimorphic venom proteins in long-jawed orb-weaving spiders (Tetragnatha) comprise novel gene families

Collection

Individuals of T. versicolor were collected by hand from three populations: along the southern fork of Strawberry Creek on the UC Berkeley Campus (UCB) (37.872°N, 122.262°W), Little Sugar Creek, Binford Farm (BF), Crawfordsville, Indiana (40.061°N, 86.853°W), and Greenville, North Carolina (ECU) (35.626°N, 77.409°W).

Venom extraction

Live specimens were transported to Lewis & Clark College where venom was extracted by electrostimulation (Binford, 2001). We obtained venom samples from 19 females and 16 males from UCB, nine females and two males from BF, and nine females and five males from ECU. To compile sufficient protein amounts for proteomics analyses, venom was pooled within sexes for each of these three populations.

RNA Isolation

Trancsriptomic analyses were performed using RNA isolated only from Indiana (BF) specimens. To capture some breadth of transcriptional timing after emptying venom glands, surviving spiders (10 females and two males) were divided into two groups within each sex, and venom glands were isolated from five females and one male two and three days after venom extraction. The glands extracted two and three days after milking were pooled within sexes and processed and analyzed separately between sexes for all subsequent analyses. To extract glands, spiders were anesthetized with CO₂ and venom glands were removed by dissection and flash-frozen immediately in liquid nitrogen. Total RNA was isolated by grinding tissues in TRIzol® reagent (Life Technologies, Carlsbad, CA, USA), adding chloroform (200 μL per mL of TRIzol®), mixing by inversion, and incubating for 20 min at 4 °C. The tube was centrifuged at 4 °C for 15 min at 14,000 rpm. An equal volume of cold 100% ethanol was added to the RNA-containing upper aqueous phase. The solution was then passed through an RNeasy® Mini Spin Column (Qiagen, Chatsworth, CA, USA) and purified, according to RNeasy® protocols.

Illumina RNA sequencing and quality control

RNA extractions from the BF population of T. versicolor males (two) and females (10) were shipped to the Genomic Services Lab at HudsonAlpha (Huntsville, AL, USA), where cDNA libraries were prepared from total RNA (poly-A isolation, Illumina TruSeq RNA Library Prep Kit v2), and 50 bp paired-end Illumina HiSeq 2500 RNA-seq was used to generate sequence reads. All QC trimming and assemblies were done using a pipeline provided by the UC Berkeley Museum of Vertebrate Zoology (https://github.com/MVZSEQ). Quality and GC content of the resulting paired-end reads was assessed using the FastQC v0.10.0 program (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). TRIMMOMATIC v0.36 (Bolger, Lohse & Usadel, 2014) and CUTADAPT v1.7.1 (Martin, 2011) were used to clean up the sequence data: nucleotides below a quality threshold of 20 were trimmed from the ends of sequences, and sequences shorter than 36 nucleotides (after trimming) were discarded. The reads were aligned to a custom library of bacterial sequences to remove prokaryotic contamination using BOWTIE2 v2.1.0 (Langmead & Salzberg, 2012). Individual paired-end files were resynchronized, removing any paired-end sequences only present in one of the two files. Left and right reads that overlap were combined into a single longer read to aid in downstream assembly using FLASH v1.2.7 (Magoč & Salzberg, 2011).

Transcriptome assembly and ORF prediction

The resulting male and female files were assembled separately, by sex, and combined for a general T. versicolor venom transcriptome. The processed reads for each sex and the combined read files were assembled using the TRINITY pipeline v2.0.6 (http://trinityrnaseq.sourceforge.net/) with default parameters except the following, group_pairs_distance = 999 and min_kmer_cov = 2. High-confidence open reading frames (ORFs) (i.e., likely coding sequences) were obtained for each gene in the transcriptome using TRANSDECODER r20140704 (Haas et al., 2013). A minimum protein length of 30 amino acids was used in ORF predictions. The completeness of each assembly was assessed via BUSCO v1.1 (Simão et al., 2015).

Read mapping to identify transcriptome dimorphisms

The processed reads from each sex were mapped against the combined assembly to identify genes that are unique to either sex, with particular emphasis placed on male-only transcripts. The mapping was performed using STAR v2.4.2a (Dobin et al., 2012) and default parameters. A custom python script was used to generate a BED file from the combined transcriptome and BEDTOOLS v2.18.1 “multicov” (Quinlan & Hall, 2010) was used to generate counts of reads from each sex mapping to combined assembly transcripts. An additional round of mapping with BOWTIE v1.1.1 (Langmead et al., 2009) followed by GFOLD v1.1.4 (Feng et al., 2012) analysis identified differentially expressed transcripts between males and females (GFOLD cutoff = 2, which approximates a log₂ fold change of two).

Functional annotation and GO enrichment analyses

The male, female, and combined assemblies were annotated via the TRINOTATE pipeline (Haas et al., 2013). This approach comprises the following steps. All contigs (BLASTX) and predicted amino acid sequences (BLASTP) were searched against the Swissprot database (downloaded 23-iv-2015). Protein domains were identified by running a HMMER v3.1b2 (Finn, Clements & Eddy, 2011) search against the PFAM (Bateman et al., 2004; downloaded 23-iv-2015) database, and signal peptides, indicating secreted proteins, were discovered with SignalP v4.1 (http://www.cbs.dtu.dk/services/SignalP/). TMHMM v2.0c (Krogh et al., 2001) was used to annotate transmembrane domains. Finally, RNAMMER v1.2 (Lagesen et al., 2007) identified rRNA transcripts. Results from database searches were loaded into a sqlite database, and GO terms were applied and used in downstream analyses.

GO enrichment analyses were performed on the combined assembly using the subset of genes endemic to male spiders, as verified by the proteomic analyses outlined below. Two data sets were analyzed, (1) all male specific proteins and (2) high molecular weight (>43 kDa) male proteins. This was done using scripts provided with the TRINITY and TRINOTATE software and the R Bioconductor package “GOseq” v1.18.0 (Young et al., 2010).

Gene family reconstruction via Markov clustering

Predicted ORFs from the male, female, and combined assemblies were combined into a single FASTA file, a BLAST database was created, and the sequences were searched (BLASTP; e-value cutoff = e10⁻⁵) against themselves (i.e., an ALLvsALL BLAST). The results were clustered into putative gene families using the Markov Clustering Algorithm (MCL v14.137; Enright, Van Dongen & Ouzounis, 2002) with default parameters and an inflation value of 2.0. The resulting e-values were negative log transformed, and the results were separated using a heuristically chosen cutoff (−1.91) and visualized in CYTOSCAPE v3.0.1 (Shannon et al., 2003). Clusters representing putative gene families were used in subsequent analyses. Initial clustering identified a single cluster of particular interest, cluster six, comprising high molecular weight proteins that was subsequently subclustered (negative log transformed and heuristically chosen edge weight cutoff = − 0.93).

Proteomic analyses

Crude venoms were dissolved in a standard buffer (5 mM CaCl₂/50 mM Tris, pH 8), pooled by sex, and shipped to the Arizona Proteomics Consortium. Venom-expressed proteins were separated by size using SDS-PAGE (12%, Fig. 1C). To increase resolution, each lane was divided into three sections and digested with trypsin followed by a clean-up step using C18 ZipTips (Millipore, Burlington, MA, USA). Tryptic peptides were analyzed using an LTQ Orbitrap Velos mass spectrometer (Thermo Fisher Scientific, Waltham, MA, USA), and the resulting MS/MS data were searched using SEAQUEST on DISCOVERER (V 1.3.0.339; Thermo Fisher Scientific, Waltham, MA, USA) against masses of theoretical fragments from a database that included our translated transcriptome sequences and all chelicerate sequences in NCBI (downloaded 4/23/2015), totaling 171,068 sequences. Matches required a fragment ion mass tolerance of 0.80 Da and a parent ion tolerance of 10.0 ppm; oxidation of methionine and carbamidomethyl of cysteine were specified in Seaquest as variable modifications. The Seaquest output was organized in Scaffold (V 4.4.3; Proteome Software Inc., Portland, OR, USA). Peptides were identified with 90% minimum threshold and 0% false discovery rate (FDR), and proteins were identified with 100% minimum threshold and 0% FDR.

Data archiving

All raw read data are available through the NCBI short read archive (SRA accession number SRP118124). Results and data files associated with proteomics (https://dx.doi.org/10.6084/m9.figshare.5378308; https://dx.doi.org/10.6084/m9.figshare.5378299), gene family clustering (https://dx.doi.org/10.6084/m9.figshare.5420986), gene ontology/annotations (https://dx.doi.org/10.6084/m9.figshare.5421655), and differential expression (https://dx.doi.org/10.6084/m9.figshare.5421517) are available via Figshare.

Transcriptome

The final trimmed read files comprised 43,550,457 ♀ (20,370,960 left; 20,370,960 right; 2,808,555 merged), 48,457,662 ♂ (22,783,475 left; 22,783,475 right; 2,890,675 merged), and 92,008,100 combined (male and female reads added together) reads. The Trinity assemblies produced 16,799 ♀ (N50 = 745; 9,904,215 total bases), 24,351 ♂ (N50 = 464; 10,664,050 total bases), and 38,021 combined (N50 = 661; 20,752,384 total bases) contigs. BUSCO estimations of completeness show the combined assembly out performs the sex specific assemblies in capturing core single copy orthologs (% missing: 87 ♀; 91 ♂; 79 combined). These completeness results likely have such high percentages of missing core orthologs due to sequencing material from such highly specific venom gland tissues. Unless otherwise noted, the combined transcriptome was used in the remainder of analyses. ORF predictions produced 23,624 putative peptides, and functional annotations were obtained for 8,075 out of 28,241 Trinity identified “genes” in the combined assembly, not to be confused with contigs that comprise isoforms and alleles of individual genes.

A flow chart summary of steps we used to identify the set of transcripts (“entities”) that are venom-expressed (detected in the proteome) and dimorphic is presented in Fig. 1B. We use “entities” to refer to the unique sets of venom components and define entities as the total population of clusters (homologous groups) and singleton sequences. MCL analysis identified 10,581 entities in the combined transcripts. Of these 87% are dimorphic—53% are found only in males and 47% in females. GFOLD analysis identified 3,800 (out of 23,634—16.07%) differentially expressed transcripts, indicating sexual dimorphisms in a wide range of genes including, but not limited to, venom cocktail peptides/proteins (1,430 female upregulated; 2,370 male upregulated).

Proteome

Crude venom separations of male and female venoms show that the profiles of expressed proteins are quite different (Fig. 1C), evidence that is supported by proteomics. LC-MS/MS produced 3,205 spectra that corresponded to 62 distinct proteins in 31 clusters, which correspond to 31 entities (Fig. 1B), only 0.29 % of total entities. Of these, nine are male-specific, non-metabolic proteins; at least eight of which have no significant homology (e ≤ 10⁻⁵) with any sequences in databases searched including NCBI and Arachnoserver (Fig. 1C; Table 1A; Herzig et al., 2010; https:// www.ncbi.nlm.nih.gov/).

The results below detail the entities that were dimorphic and confirmed as venom components through detection in proteomes.

Cluster 6

The most abundant, dimorphic cluster, the sixth most highly represented in the combined transcriptome (66 homologous proteins, 38,922 mapped reads), included the proteins with the highest number of distinct proteins (12) and represent 30.7 % of MS spectra in the male proteome (Table 1). Transcripts in this cluster appear to code for proteins of sizes consistent with the large proteins unique to male venoms (Fig. 1C). While none of the assembled transcripts are full length (initiating methionine through stop codon), individual transcripts in this cluster translate into proteins ranging from 41 to 1,093 aa. Multiple sequence alignment of these proteins generates a consensus sequence of 1,158 aa in length, and the longest single transcript in the alignment (1,093 aa) has a predicted MW of 128.95 kDa. The homologs in this cluster are grouped by MCL into eight sub-clusters and five singletons (Fig. 2).

Proteins in cluster 6 contain multiple repeating units, each of which has ∼150 aa and a conserved pattern of 10 cysteine residues (Fig. 3). Submission of a multiple alignment of these repeats to the Fold and Function Assignment (FFAS) server (http://ffas.sanfordburnham.org/ffas-cgi/cgi/ffas.pl; PMID 15980471) returns strong evidence of distant homology to Argos, a 223-residue antagonist of epidermal growth factor receptor signaling with a known structure (PDB ID 3c9a). Argos contains three small β-sheet rich domains, the first two of which correspond to one of the sequence repeats present in the cluster 6 proteins, with a similar pattern of 10 cysteine residues making five disulfide bonds (Fig. 3B; Klein et al., 2008). The third Argos domain corresponds to an extra half-repeat with six of the 10 cysteines and three disulfide bonds (Fig. 3B). Argos uses contacts from multiple domains to bind and sequester small protein ligands of EGF receptors. Based on the FFAS score (−15), it is quite likely that the cluster 6 proteins have a similar fold, disulfide pattern, and domain organization to Argos; however, the sequence homology to Argos is so distant (<20% sequence identity between Argos domains 1 and 2 and any sequence repeat in the Tetragnatha proteins as shown in Fig. 3A) that a functional similarity is much less certain.

Other male-specific venom proteins

Seven additional male-specific clusters include transcripts that code for proteins that span 20–37 kDa. These are sizes that correspond to predicted molecular weights of full-length proteins within each cluster (Fig. 1C). Cluster numbers are labeled next to individual bands in Fig. 1C and correspond to the relative rankings based on representation in the combined transcriptome: cl-2448 (>36 kDa), cl-217 (37 kDa), cl-218 (35 kDa), cl-1113 (26 kDa), cl-9248 (>26 kDa), cl-7374 (>24 kDa), and cl-2111 (21 kDa). All were detected in the proteome, though not as abundantly as cluster 6. Proteins encoded by genes in all of these clusters correspond to novel gene families without homologous sequences present in existing databases, including GenBank and Arachnoserver. Proteins in each of these clusters have multiple cysteine residues (8–21), and with the exception of clusters 217 and 218, they have at least two C × C motifs.

In addition to identifying novel families of proteins in these venoms, we also identified several sequences with evidence of homology to known enzymes. Venom proteins in males hit three large peptidases: M13 metalloendopeptidase (specifically neprilysin, ∼90 kDa), M14 carboxypeptidase M (sequence not full length, but estimated to be >55 kDa), and S8 peptidase (specifically neuroendocrine convertase, ∼72 kDa). Only one large peptidase—M13 metalloendopeptidase—was identified in the venom proteome of females; however, homologous sequences to each were found in both transcriptomes. Both S8 and M14 peptidase family members are known to be involved with activation and processing/regulation of hormones, respectively, whereas M13 peptidases are comprised of GluZincins, a superfamily of peptidases that act on molecules <∼40 aa (MEROPS peptidase database, http://merops.sanger.ac.uk/index.shtml) and have also been reported in venom of a trapdoor spider (Undheim et al., 2013). The degradative enzymes hyaluronidase and phospholipase A2 (PLA2) were identified in the venom of females, constituting ∼1% of the proteome. Cluster 19 contains peptides corresponding to various hemocyanin subunits and represents the most abundant set of proteins detected in the female venom proteome (Table 1B). Each subunit varies in size, and the most prominent was subunit G, which is predicted to be ∼72 kDa (Fig. 1C). The female venom is also rich in small cysteine-rich peptides, corresponding to clusters 2335, 846, and 8,293 identified in the proteome (Table 1A).

Within gene families observed only in male proteomes, we recovered interesting patterns of sexually dimorphic expression and potential functions. Despite some mRNA from females mapping to these transcripts, their peptides were not detected in the female proteome. This could be due to a lack of translation following transcription or perhaps the proteins are not present in the venom. GO enrichment analyses performed on all male-only peptides, based on the proteomics analysis, indicated potential non-feeding functions of these proteins. Four GO terms were significantly enriched, two relating to hormone functions—GO:0016486 (BP peptide hormone processing; FDR p-value = 0.0143), GO:0008237 (MF metallopeptidase activity; FDR p-value = 0.0349), GO:0042445 (BP hormone metabolic process; FDR p-value = 0.0349), and GO:0006518 (BP peptide metabolic process; FDR p-value = 0.0497). Half of the enriched GO terms were specifically associated with hormone functions, suggesting this venom-based sexual dimorphism could be involved in sexual communication.