Factors that increase reproductive variance among individuals act to reduce effective population size (Ne), which accelerates loss of genetic diversity and decreases efficacy of purifying selection. These factors include sexual cannibalism, offspring investment, and mating system. Pre-copulatory sexual cannibalism where the female consumes the male prior to mating exacerbates this effect. We performed comparative transcriptomics in two spider species, the cannibalistic Trechaleoides biocellata and the non-cannibalistic T. keyserlingi, to generate genomic evidence to support these predictions. First, we estimated heterozygosity and found that genetic diversity is relatively lower in the cannibalistic species. Second, we calculated dN/dS ratios as a measure of purifying selection, higher dN/dS ratio indicated relaxed purifying selection in the cannibalistic species. These results are consistent with the hypothesis that sexual cannibalism impacts operational sex ratio and demographic processes, which interact with evolutionary forces to shape the genetic structure of populations. However, other factors such as the mating system and life-history traits contribute to shape Ne. Comparative analyses across multiple contrasting species-pairs would be required to disentangle these effects. Our study highlights that extreme behaviours such as pre-copulatory cannibalism may have profound eco-evolutionary effects. 

Spider collections and RNA extractions

We collected six adult males from each species (T. keyserlingi and T. biocellata), consisting of two males from three populations where the two species co-occur in Uruguay 1) Parque San Miguel (33°41′51″S 53°32′00″W) Rocha; 2) Quebrada de los Cuervos (32º55’39’’S 54º27’25’’W), Treinta y tres and 3) Paso Centurión (32°8′2″S 53°47′55″W) Cerro Largo. We additionally collected two males from Paratrechalea ornata, from Quebrada de los Cuervos, as an outgroup. We chose only males to avoid any biases related to possible sex differences in gene expression. Once in the laboratory, we flash froze them in liquid nitrogen and stored them in -80ºC until RNA extraction. We extracted RNA from half of the cephalothorax using RNAasy Kit–QIAGEN following the manufacturer’s recommendations. Then, we stabilized the extractions in an ethanol (EtOH) precipitation. Library preparation (TruSeq Stranded mRNA LT Sample Prep Kit) and sequencing (PE, 101bp, Illumina platform) were done by Macrogen Inc.

Quality control and de novo transcriptome assembly 

Prior to de novo transcriptome assembly, we quality trimmed reads with Trimmomatic v.0.36 (Bolger et al., 2014) using default parameters. We performed de novo species assemblies using Trinity v2.12.0 with default parameters (Grabherr et al., 2011). We used BUSCO v5.2.2 (Benchmarking Universal Single-Copy Orthologous) (Manni et al., 2021) to evaluate the proportion and quality (complete, fragmented and duplicated) of the spider orthologs genes present in each assembly (Parasteatoda tepidariorum reference). We then used CD-HIT (Li & Godzik, 2006; Fu et al., 2012) with EST mode with a threshold of 0.95 on each assembly to reduce redundancy. From now on, we will refer to these last CD-HIT outputs assemblies as just assembly. We show individual raw reads statistics in Table S1, and BUSCO completeness results for species assemblies in Figure S1.

Genetic diversity

We aimed to assess heterozygosity for each Trechaleoides species and P. ornata. First, we mapped all Trechaleoides and P. ornata individuals reads to their corresponding species assembly using bowtie2 (Langmead & Salzberg, 2012). Second, we performed the SNP calling and created the vcf files using Bcftools package (Danecek et al., 2021) with a combination of mpileup (-g) and call (multialellic-caller (-m)) methods. We filtered sites using (Bcftools package), keeping only those with cover higher than 20X and mapping quality (as output by bowtie) over 20. We calculated Hobs as the number of heterozygous sites within the number of the total sites for each individual. To test the difference in Hobs between species, we used R software (R Team Core, 2022) and performed a Generalized Linear Mixed Model with Binomial distribution (GLMM (B)) including species as a fixed effect and population as a random effect. 

Selection intensity

We aimed to estimate dN/dS ratios in species with and without pre-copulatory sexual cannibalism. First, we used Transdecoder v5.5.0 (Haas & Papanicolaou, 2016) on each species assembles to predict the longest isoforms and the nucleotide coding sequences (CDS), then we used these nucleotide CDS to identify putative 1-to-1 orthologs shared between species using Orthofinder v2.5.4 (Emms & Kelly, 2019). Afterwards, we translated the nucleotide sequences to their corresponding amino acid sequences and aligned them using MUSCLE v 3.8 (Edgar, 2004). We created nucleotide sequence alignments using the protein alignment as a reference in PAL2NAL v14 (Suyama et al., 2006) and performed a test of positive selection at the species level using the Codeml program implemented in the PAML package (Yang, 2007). We evaluated the evolutionary rate of each terminal branch in the species tree, by concatenating all individual genes alignments and estimating rates of nonsynonymous (dN) substitutions, synonymous (dS) substitutions, and their dN/dS ratio for each branch. With PAML, we bootstrapped codon columns from the concatenated alignment (n = 1000), and using the branch-free ratio model (model = 1) we estimated an independent dN, dS and dN/dS for each terminal branch in the phylogeny for each bootstrapped alignment (Yang & Nielsen, 1998). We, then, estimated 95% confidence limits of dN, dS, dN/dS for each terminal branch using the distribution of the 1000 bootstrapped estimates using R (R Team Core, 2022). For the model, we used the phylogenetic tree inferred by Orthofinder: ((T. biocellata, T. keyserlingi) P. ornata).