Female P. antipodarum are ovoviviparous (Winterbourn, 1970). We determined whether female snails were reproductively active by either observing embryos through the shell (only possible in light-shelled lineages) or by dissecting the female and opening the brood pouch. Only the ovaries from females that contained developing offspring were extracted. We experienced major challenges in finding enough reproductively active sexual females (a common challenge in laboratory-raised P. antipodarum), which constrained our ability to replicate within and across sexual lineages. All dissections were performed in RNAlater (Thermo Fisher Scientific, Waltham, MA, United States). The dissected ovaries were submerged in 150 μL of RNAlater (Sigma-Aldrich) in a sterile microfuge tube, incubated at 4°C for 24 h to allow the RNAlater to saturate tissue, and then stored at −80°C until RNA extraction. P. antipodarum ovaries are so small (<0.01 mm³ in volume) that adequate RNA isolation for RNA sequencing at the time of data collection in 2014 required pooling of ovarian tissue from at least six females per sample. Accordingly, for each sample other than the Selfe sexuals, ovaries from six females from the same lineage were pooled into a sample for RNA sequencing, with two biological replicates/lineage (12 females total for each lineage). For the Selfe sexuals, the lack of reproductively active females meant that we were only able to sequence one pool of six females.

We homogenized ovary tissue in Trizol reagent (Invitrogen, Carlsbad, CA, United States) using a TissueLyser II (Qiagen, Valencia, CA, United States) with a 5-mm bead and 2 rounds of shaking at 20 cycles/sec. for 2 min. Following the addition of 0.2 volumes of chloroform, we used Phase Lock Gel (5PRIME, Gaithersburg, MD, United States) to achieve phase separation. Next, we isolated total RNA from the aqueous phase with RNA Clean & Concentrator-5 columns (Zymo Research, Irvine, CA, United States) and treated this RNA with Dnase I (Qiagen) to remove contaminating DNA. We then quantified total RNA by fluorometry with Ribogreen (Invitrogen) on a Turner TBS-380 (Promega, Madison, WI, United States), and we assayed RNA quality with an Experion (Bio-Rad, Hercules, CA, United States).

We then used the Illumina TruSeq RNA Library Prep Kit v2 (Illumina, San Diego, CA, United States) to generate each RNA sequencing library from 1 μg of total RNA and sequenced this RNA on an Illumina HiSeq 2000 (2 × 100 bp). This sequencing effort resulted in a mean (±SD) number of reads per run of 28,610,000 (±16,410,000).

We used the Fastx-toolkit (Gordon and Hannon, 2010; version 0.0.14) to filter out low-quality reads by retaining only those in which at least 75% of bases had PHRED scores ≥ 20 (mean ± SD of reads following filtering was 27,700,000 ± 15,340,000) and assessed the quality of the reads with FastQC (Andrew, 2010; version 0.10.1). Next, we used Trinity (Grabherr et al., 2011; Hass et al., 2013; version 20131110) for in silico normalization of each read pair (maximum read coverage of 30 and otherwise default settings) to reduce memory use and improve computation ease for de novo assembly. We then pooled the normalized reads from all of the lineages into one assembly using a group pair distance of 350 and otherwise default settings. Following de novo assembly, we used the Trinity plug-in tool TransDecoder (Hass et al., 2013) to identify the most likely coding sequences of the transcriptome and filter out misassembled transcripts. Next, we used CD-HIT-EST (Huang et al., 2010; version 4.6.4) to hierarchically cluster contigs based on 95% identity to reduce redundancy in our transcriptome. Finally, to identify possible contaminant sequences, we used the Blobology (Kumar et al., 2013) pipeline script ‘blast_taxonomy_report.pl’ to assign each transcript to a kingdom and domain based on BLASTX results and removed any non-metazoan transcripts. We annotated the final transcriptome using BLASTX via Blast+ (Camacho et al., 2009; version 2.6.0) against the NCBI non-redundant protein database with an E-value cutoff of 1e-5. Finally, we used Blast2GO (Conesa et al., 2005; version 3) to assign gene ontology (GO) terms, InterPro IDs, and KEGG pathways to each transcript.

We used the final version of the ovary transcriptome as a reference for differential gene expression analysis. We mapped each set of reads to this reference using Tophat2 (Kim et al., 2013; version 2.0.9) and compared gene expression with Cuffdiff 2 (Trapnell et al., 2010, 2012; Kim et al., 2013; version 2.1.1) via two approaches. The first gene expression analysis that we conducted, hereafter “sex-asex,” treated each of the five sexual and six asexual sequence samples as biological replicates for the conditions of “sexual” or “asexual,” with the goal of broadly comparing gene expression between sexual and asexual P. antipodarum. This approach was intended to have sufficient power for detecting changes in gene expression that discriminated between sexuals and asexuals, with the caveat that we might also inflate the likelihood of bias or error (e.g., detecting false positives and negatives) for differential expression by treating samples from the same lineage as independent data points and with uneven representation between original population sources (e.g., two sexual Alexandrina lineages and no sexual Sarah lineages). Furthermore, this analysis does not account for the possible effect lake of origin could have on gene expression differences. Despite its limitations, this analysis was the best we could do with available snails to take initial steps toward our focal question regarding expression differences between sexual and asexual P. antipodarum.

With this logic and goal in mind, we sought to narrow down our candidate gene field by using Cuffdiff 2 to identify the most consistently differentially expressed sequences between sexual and asexual P. antipodarum, here treating all six laboratory lineages (represented by two sequence samples each, except for the single sample from the sexual Selfe) as independent conditions and conducting 15 pairwise gene expression comparisons, again with a false-discovery rate applied via Benjamini-Hochberg correction for multiple tests (q < 0.05). This approach was intended to reduce the possible noise introduced in the sex-asex differential expression comparison and identify sets of genes broadly relevant to differences in reproductive mode and/or specific to distinct transitions to asexuality. We acknowledge that we lose power in this “pairwise” approach relative to the sex-asex comparison because of low or absent replication per lineage, which prevents us from more confidently characterizing population or lineage-specific changes in gene expression and the addressing the possibility of distinct transcriptional signatures of asexuality between different asexual lineages. We therefore relied on these pairwise results primarily for global expression profiles to evaluate the roles of population history and reproductive mode on gene expression and aid our effort to identify candidate genes by highlighting the genes with minimal expression variation outside of the sexual vs. asexual comparison. We similarly used CuffDiff 2 to compare gene expression between the two biological replicates for each of the five lineages for which we had replicated sequence data to assess within-lineage variability. The main caveat that applied to this analysis is that these particular CuffDiff 2 comparisons lack replication (e.g., one sequenced sample from the asexual Alexandrina lineage vs. another sample from the same lineage). Nevertheless, this comparison appears to demonstrate consistent expression between each set of replicates (Supplementary Figure 2) relative to much higher variability across lineages. This result indicates that the lack of replication for the sexual Selfe lineage may not have strongly affected our ability to make comparisons between this lineage and the five lineages for which we did have biological replicates.

To identify transcripts that are particularly likely to be directly associated with sexual vs. asexual reproductive mode, we searched GO annotations to identify 702 genes with putative roles in (sexual) reproduction, oogenesis, embryogenesis, meiosis, germline, ovaries, female-specific functions, mitosis, chromosomes, and cell cycle, hereafter referred to as “focal sequences” for subsequent analyses. This collection of sequences is far from comprehensive and was used as a means to test for potential expression differences that were not captured with the GO-enrichment analysis of the sequences differentially expressed between sexual and asexual snails. For example, a greater proportion of this “focal sequence” set may be differentially expressed between sexuals and asexual than the rest of the transcriptome. To address this possibility, we also collected 702 transcripts for which we used a random number generator to select from the transcriptome as a whole as a control set (both sets of transcripts are listed in Supplementary Tables 5, 6). We used a chi-square test to determine if the “focal” and “control” sets of genes had different proportions of differentially expressed genes (“significant” vs. “non-significant”) in the sex-asex expression analysis results.

We used F_ST analyses to determine relative levels of genetic differentiation among lineages and to examine within vs. among-lineage differentiation. We also used outlier analyses to identify genes with consistently high differentiation in sexual vs. asexual comparisons. We first filtered the ovary transcriptome to the 11,428 transcripts that were expressed (FPKM > 0) by all six lineages. We then again normalized read counts of each read pair with Trinity’s in silico read normalization to a maximum read coverage of 60 to ease memory demands for downstream analyses and use comparable coverage depths for SNP analysis; we otherwise used default settings. Next, we used Tophat2 (with the same parameters as in the gene expression analyses) to map reads from each lineage to this reduced reference assembly. We then used Picard Tools¹ (version 1.99) and elements of the GATK (McKenna et al., 2010; version 3.4) pipeline according to recommendations for using RNA-seq in population genomics (De Wit et al., 2015) and following Bankers et al. (2017) to prepare the BAM files produced by Tophat2 for variant discovery. In short, we used Picard Tools AddOrReplaceReadGroups to tag reads with sample identifiers and MarkDuplicates to identify and remove duplicated reads. We then used the GATK SplitNCigarReads to reassign mapping qualities and then realigned indels with RealignerTargetCreator and IndelRealigner. Next, we used SAMtools (Li et al., 2009; version 0.1.19) to create mpileup files from the processed bam files. From these bam files, we calculated F_ST for all SNPs using Popoolation2 (Kofler et al., 2011; version 1201) for all pairwise comparisons between lineages and between biological replicates from the same lineage.

To evaluate the effect of lineage and lake of origin on population differentiation we compared mean pairwise F_ST values from within lineages vs. between lineages and within lakes vs. between lakes with Mann–Whitney U (MWU) tests. We also identified F_ST outlier SNPs as those SNPs with F_ST values above the upper 95% confidence interval of the mean F_ST per SNP for all pairwise comparisons between lineages using GraphPad Prism (version 7.0). We then compared sets of transcripts that contained outlier SNPs across the various pairwise comparisons to identify genes that appear to be consistently differentiated between sexual vs. asexual lineages. We conducted GO enrichment analysis with Blast2GO on sets of transcripts with outlier SNPs in from all pairwise comparisons and identified GO term enrichment for transcripts with outlier F_ST SNPs in both Alexandrina sexual vs. asexual lineage comparisons and the Selfe sexual vs. asexual lineage comparisons to identify the types of genes, if any, differentiating sexuals and asexuals across multiple populations.

We also used a phylogenetic approach to demonstrate the evolutionary relationships among the six lineages used in this study. First, we assembled transcriptomes for each lineage separately following the same methods described earlier. In brief, we used Trinity to assemble transcripts followed by transdecoder to retain only the coding sequence and CD-HIT-EST to cluster highly similar sequences. Next, we performed reciprocal blastn searches between pairs of transcriptomes with “reciprocalblast_allsteps.py” (Warren et al., 2014) using an E-value cutoff of 1e-20 to identify sequence pairs that were one another’s best hit for each blast search. We then identified sequences from these pairwise comparisons that were reciprocal best blast hits in all 15 comparisons, i.e., 1:1:1:1:1:1 best hits. We considered these 6,757 sequences orthologous and retained only these sequences for further analysis.

We then used MUSCLE (Edgar, 2004; version 3.5) to produce individual alignments for each of the 6,757 orthologous sequence sets. Because we had processed these sequences in TransDecoder, each sequence set was already in frame prior to alignment. Consequently, we kept only alignments in which no gaps were introduced to any sequence during the alignment process for subsequent analysis. Next, we concatenated the 1,225 single-sequence alignments into one file containing 1,135,416 aligned positions for phylogenetic analysis. We then used the 3rd position sites from this alignment within a maximum-likelihood framework and Kimura 2-parameter computed distances to generate an unrooted phylogeny in MEGA (Kumar et al., 2016; version 7) and used 1,000 bootstrap replicates to evaluate phylogeny support.

In addition to gene expression analyses, we used these RNA-seq data to evaluate the (1) population genetic differentiation on a per-SNP basis across the transcriptome to determine if lineages from the same lake are more similar genetically than between lakes and (2) phylogenetic history of the sequenced lineages, to address the hypothesis that population history has a strong influence on gene expression, regardless of reproductive mode.

These analyses revealed substantially lower F_ST values for comparisons within vs. between lineages, with the lowest mean F_ST/SNP value from an inbred sexual line (0.002723 ± 0.005571; Figure 1A), and the highest mean F_ST/SNP (∼35× greater than within the inbred sexual line) from the comparison of the asexual Alexandrina lineage and sexual Selfe lineage (0.09754 ± 0.162; Figure 1A). When all within-and between-lineage mean F_ST/SNP values are rank ordered, a clear pattern of distinctly lower values for within vs. between-lineage differentiation emerges (Supplementary Figure 2; p < 0.025 Bonferroni MWU). We also found that lineages from the same lake tend to be more similar to one another than to lineages from different lakes (Supplementary Figure 2; p < 0.025 Bonferroni MWU). The Alexandrina lineages in particular reveal a clear signal of increasing genetic differentiation across a spectrum of comparisons, with minimum genetic differentiation within lineages, intermediate values of genetic differentiation for Alexandrina sexual lineages, higher (but still intermediate) genetic differentiation for sexual lineages vs. the asexual lineage from Alexandrina, and the highest genetic differentiation for any Alexandria lineage vs. a lineage from a different lake (Figure 1A). We do not see a systematic grouping of the asexual vs. sexual comparisons for mean F_ST, suggesting that the triploidy of the asexual lineages does not systematically affect genetic differentiation. A good example of the absence of triploid/asexual grouping is provided by the Alexandrina sexual vs. asexual comparisons, which have lower mean F_ST than Alexandrina sexual vs. Selfe sexual comparisons.

Similar results emerged from the phylogenetic approach. We aligned orthologous sequences from lineage-specific transcriptome assemblies and concatenated the gapless alignments to generate a phylogeny for these lineages. We identified 6,757 orthologous sequences that were reciprocal best-blast hits across all six lineages. We aligned these sequences individually and then concatenated the 1,225 gapless sequence alignments into one alignment with 378,472 third codon sites. The resulting phylogeny is very similar to the biogeographic history, with the six lineages grouped into an Alexandrina clade and a Selfe and Sarah clade (Figure 1B). This result is also consistent with our F_ST analysis: the asexual Selfe and Sarah lineages had the lowest mean F_ST of any cross-lake comparison (0.0728 vs. 0.0847–0.0975 for others) and was lower than the mean F_ST from the Selfe sexual vs. the asexual Selfe snails (0.0827). We speculate that the close relationship of the asexual Sarah and Selfe lineages might reflect a recent invasion of Lake Sarah by an asexual lineage from Selfe; Neiman and Lively (2004) found that all 35 sampled snails from Sarah harbored the same recently derived mitochondrial haplotype, which was also found in 19/23 snails collected from Selfe. The very close proximity of these two lakes (22 km as the bird flies) is consistent with this scenario.

The pairwise analysis of gene expression allowed us to perform a variety of comparisons of patterns of gene expression (e.g., sex vs. asex, sex vs. sex, etc.) and identify transcripts consistently differentially expressed between sexual and asexual lineages (FPKM and statistical results from pairwise analysis listed in Supplementary Table 8). Mean proportions of transcripts with differential expression did not differ across sexual vs. sexual or asexual vs. asexual lineage comparisons, sexual vs. asexual lineages, or comparisons of any two lineages (Supplementary Figure 3). With the caveat that these pairwise comparisons necessarily involved some pseudoreplication and will not be able to identify lineage-specific drivers of transitions to asexuality, these results indicate that expression profiles between sexuals and asexuals do not seem to be any more different, in terms of proportion of transcripts differentially expressed, than are any two sexual or asexual lineages. These comparisons also revealed that within-lineage comparisons (i.e., differential expression analysis between biological replicates) were much less likely to reflect differential expression (>14× smaller proportion of transcripts) than any across-lineage comparison (Supplementary Figure 3). While more extensive analyses, including RNA-seq of many individuals from the same lineage raised separately in order to account for culture effects, will be needed to further test the influence of genetic history vs. maternal effects and other shared environmental factors, our results support a scenario where snails from the same lineage tend to express genes at very similar levels to one another relative to snails from other lineages.

A dendrogram-based visualization of overall gene expression profiles revealed that lineages originating from the same lake populations tend to cluster relative to lineages from other lakes (Figure 2A), regardless of reproductive mode. This result suggests that lake of origin is a more important predictor of expression profiles than reproductive mode. We obtained similar results with a multidimensional scaling (MDS) plot (Figure 2B and Supplementary Figure 4), in which the M1 axis groups the three Alexandrina lineages together and splits the sexuals and asexuals on the M2 axis. The M2 axis grouped the gene expression profiles of the Selfe sexual lineage with the Sarah asexual lineage, and the M1 axis separated the Selfe and Sarah lineages from the Alexandrina lineages. Principal component analyses resulted in similar groupings, where PC2 separates these lineages into an Alexandrina group and a Selfe and Sarah group (Figure 2C and Supplementary Figure 5). As noted above, lakes Sarah and Selfe are physically much closer to each other than either lake is to Lake Alexandrina. Altogether, the relationship of these expression profiles is consistent with earlier studies in P. antipodarum suggesting that geographic distance and biogeography are associated with population-level differences in traits (e.g., Verhaegen et al., 2018) including gene expression (e.g., Bankers et al., 2017). Notably, the earlier evidence for population-level clustering of gene expression in P. antipodarum were from field-collected snails (Bankers et al., 2017), while the snails analyzed in this study had been raised in common garden conditions for several generations, suggesting that these field effects persist even under laboratory culture and thus likely represent inherited expression traits. Our results therefore add another layer of understanding of gene regulation in this system.

We used SNPs with outlier F_ST values to identify transcripts and gene categories that are significantly differentiated between sexual and asexual P. antipodarum. For this analysis, we used Blast2GO to determine functional enrichment for F_ST outlier-containing transcripts for all lineage pairwise comparisons. To identify gene classes that might differ according to reproductive mode and to minimize the contribution of lineage effects to this classification effort, we identified GO terms that were overrepresented in sets of transcripts with outlier F_ST SNPs harbored by both Alexandrina sexual lineages compared to the Alexandrina asexual lineage, while excluding overrepresented terms from the comparison of the two Alexandrina sexual lineages. This collection of 78 GO terms included nine categories related to mitosis, cell cycle, and DNA damage repair and 17 categories related to respiration and mitochondrial functions. The latter group is of interest in light of prior evidence for a higher load of deleterious mutations in the mitochondrial genomes of asexual vs. sexual P. antipodarum (Neiman et al., 2010; Sharbrough et al., 2018). We also conducted GO function enrichment analysis for transcripts with F_ST outliers for the Selfe sexual vs. asexual lineages and the Selfe sexual lineage vs. the Sarah asexual lineage. This analysis revealed that Selfe sexual vs. asexual transcripts with F_ST outliers were enriched for 277 GO terms. Similar to the Alexandrina sexual/asexual analysis, while many (e.g., 27/87 cellular component categories) of these overrepresented categories were also associated with mitochondrial function, GO functions related to “single fertilization,” “fertilization,” “zona pellucida receptor complex,” “binding of sperm to zona pellucida,” and “sperm-egg recognition” were also overrepresented for the Selfe sexual/asexual comparison. Finally, 20 GO terms relating to mitosis, cell cycle, and DNA damage repair were overrepresented in the Selfe sexual vs. Sarah asexual F_ST comparison.

The only GO term categories enriched for transcripts differentially expressed in sexual vs. asexual lineages and transcripts with F_ST outlier SNPs for sexual vs. asexual comparisons in both groups of lineages (i.e., lineages from Alexandria and ones from Selfe and Sarah) were scavenger receptor activity, cargo receptor activity, and carbohydrate binding. These receptor activities are notable for being part of immune responses and therefore could represent genes that are rapidly evolving in response to pressure from parasites (e.g., Ameline et al., 2021), perhaps not unexpected in this textbook model system for host-parasite coevolution (e.g., Emlen and Zimmer, 2020). Also notable is that P. antipodarum experiences extensive across-population variation in the frequency of parasite infection (e.g., Lively, 1987) as well as population-specific transcriptomic responses to trematode infections (Bankers et al., 2017), meaning that different populations likely experience different selection for immune function.

Chitin-related enzymes were both overrepresented as differentially expressed between sexual and asexual snails in the broad sex-asex comparison and had outlier F_ST SNPs in comparisons of sexual vs. asexual lineages in lakes Alexandrina and Selfe. “Chitin catabolic process” and “chitinase activity” were also overrepresented GO terms for transcripts with outlier SNPs for the Selfe sexual vs. Sarah asexual F_ST comparison. Various chitinases are active throughout mollusk development and might play a role in processes from shell development to cell proliferation (Badariotti et al., 2006; Weiss et al., 2006; Schönitzer and Weiss, 2007; Yonezawa et al., 2016), meaning that the differential expression of these genes in sexual and asexual P. antipodarum (and, often, overexpression in the latter) is of particular interest in light of the higher growth rate of asexual snails (Larkin et al., 2016). Both sexual and asexual P. antipodarum are ovoviviparous, featuring internal fertilization followed by internal embryo brooding until birth, so by sampling ovarian tissues of reproductively active snails, we might have also captured expression of early embryos.

The GO terms “chitin binding” and “chitin metabolic process” were also overrepresented in genes differentially expressed between cyclical and obligate parthenogenic monogonont rotifers (Hanson et al., 2013), though this difference between sexual and asexual forms of these rotifers might instead be related to the major differences in egg casing between the sexually produced resting eggs and the eggs produced by asexual females (Hanson et al., 2013). Differences in egg casings in sexual vs. asexual settings may contribute to why chitin binding genes are upregulated in sexual relative to asexual Daphnia similoides (Zhang et al., 2016). It is also possible that the apparently adaptive nature of at least some of the morphological plasticity in shell shape/size in asexual P. antipodarum in the United States and Europe (e.g., Kistner and Dybdahl, 2013; Verhaegen et al., 2018) might also be linked to the expression of chitin-related genes.

We observed SNP F_ST outliers in two meiosis- related genes, rad51 and msh2, among the three Alexandrina lineages. In particular, we found that for comparisons involving any two Alexandrina lineages, each of these two genes had at least one outlier SNP. For msh2, the two sexuals shared the same outlier SNP when compared with the asexual lineage. This SNP did not appear as an outlier when the two Alexandrina sexual lineages were compared to one another. Neither gene had outlier SNPs for the Selfe sexual vs. asexual lineages or the Selfe sexual vs. Sarah asexual lineage comparisons, and msh2 was not differentially expressed among any of the Selfe lineages. There was no clear pattern or association with reproductive mode for rad51. It is important to note that msh2 and rad51, along with other meiosis-related genes, tended to be expressed at relatively low levels. This low expression likely limits our statistical power to detect significant expression differences in the absence of much deeper coverage than available in our study. While we cannot exclude the possibility that we did not capture enough tissue and/or the correct developmental stage for a rigorous analysis of meiotic gene expression in P. antipodarum, it is relevant to point out that absence of evidence for expression differences in meiosis genes between sexual and asexual lineages have also been reported for obligate vs. cyclical parthenogenetic strains of the monogonont rotifer Brachionus calyciflorus (Hanson et al., 2013) and during meiosis and parthenogenesis in Daphnia pulex (Schurko et al., 2009).

The sex-asex differential gene expression analysis identified 1,103 transcripts with significantly different (FDR q < 0.05) expression levels between the two reproductive modes. Of these transcripts, 660 (3.39% of the total transcripts tested) were expressed at significantly higher levels by asexual P. antipodarum, and 443 (2.28% of the total transcripts tested) were expressed at significantly higher levels by sexual P. antipodarum (Figure 3A and Supplementary Figure 6). This general pattern of more transcripts with significantly higher expression in asexuals relative to sexuals extended to our pairwise comparisons between sexual and asexual lineages. Here, at lineage-level differential expression comparisons, asexuals overexpressed, on average (±SD), 2293 (±336.5) transcripts and underexpressed, on average, 1897 (±66.4) transcripts relative to sexuals (Figure 3B).

We also used the pairwise analysis to assess whether and which transcripts were consistently differentially expressed between sexual and asexual lineages to do the best we could to account for the limitation that our sex-asex comparison could not control for lake or lineage. We narrowed down the pairwise expression results to transcripts only differentially expressed in sexual vs. asexual lineage comparisons, excluding transcripts differentially expressed in any sexual vs. sexual or asexual vs. asexual comparison. This exclusion could inadvertently eliminate results relevant to population or lineage-specific transitions to asexuality, but because we did not generate data sufficient for such characterization, we chose to focus on identifying genes with the strongest evidence for reproductive-mode driven expression differences. This approach enabled us to in part account for the caveats of the sex-asex differential expression analysis (e.g., not structured to consider population history or lineage, which could lead to the appearance of sexual vs. asexual differential expression being driven by lineage-specific overexpression) to determine if common expression differences exist between sexual and asexual P. antipodarum regardless of genetic background. In light of the design limitations with respect to sequencing replication per lineage and the uneven representation of lake populations, we caution that these results are not definitive but instead represent the likeliest set of genes related to reproductive differences that we could identify with the data we generated.

From the pairwise comparisons, we identified 347 sequences only differentially expressed between reproductive modes, and not within reproductive modes. 86 of these sequences (24.8%) were also identified in the sex-asex analysis. Of those 347 sequences, 89 sequences were upregulated by at least two sexual lineages relative to asexuals or at least two asexual lineages relative to sexuals; 60 of those 89 sequences (67.4%) were identified in the sex-asex approach. Of the 23 sequences that were differentially expressed in all sexual vs. asexual pairwise comparisons 21 were identified with the sex-asex approach (91.3%). The increasing degree of overlap between the sex-asex results and specificity of the pairwise results indicate that the pairwise approach was useful in identifying consistently differentially expressed genes in sexual vs. asexual snails. Furthermore, the pairwise results demonstrate that a more thorough consideration for genetic background will need to be taken into account when designing future RNA-seq experiments aimed at uncovering regulatory differences in sexual vs. asexual (and different ploidy levels) gene expression in P. antipodarum.

Similar to Bankers et al. (2017), where almost no transcripts (3/62,862) were consistently differentially expressed between infected and uninfected P. antipodarum across source populations, we see relatively few transcripts consistently differentially expressed in all sexual vs. asexual lineage comparisons (Figure 4). Of those 23 transcripts, 10 were upregulated in sexuals relative to asexuals, and 13 transcripts were upregulated in asexuals relative to sexuals in every comparison (Figure 4, listed in Supplementary Tables 9–14). We applied the SuperExactTest (Wang et al., 2015), which demonstrated that these overlaps in sequences consistently expressed higher in sexuals vs. asexuals and vice versa were larger than expected by chance (10 vs. 0.00180219 in sexuals and 13 vs. 0.0007407894 in asexual; p = 2.031013e-35 and p = 6.984346e-53). Only four of the 10 consistently upregulated transcripts in sexuals were annotated. One of these transcripts was annotated as kinetochore protein spc24, an essential component of the NDC80 complex, which is necessary for chromosome segregation (McCleland et al., 2003). In yeast, the regulation of NDC80 subunits is integral for meiosis I to occur (reviewed in (Chen and Ünal, 2021). Of potential relevance to the P. antipodarum system is that a lineage of baker’s yeast with a unique non-synonymous mutation in NDC80 was also found to have a particularly high incidence of polyploidy (83% vs. 34% of all 144 strained analyzed) (Zhu et al., 2016). Also noteworthy is that one of the 13 transcripts consistently upregulated in asexual relative to sexual lineages is annotated as DNA primase large subunit, which is involved with DNA synthesis and replication. This last result is of interest in light of the polyploid status of asexual P. antipodarum, which have more DNA to replicate during cell division than diploid counterparts. While these examples of genes that are differentially expressed across reproduced modes are of potential relevance to our overarching goals regarding the identification of candidate loci involved in the transition to asexuality, they are not informative enough to construct a model of how asexuality is achieved in P. antipodarum. Characterizing the molecular evolution of these genes across a diverse and representative sample of P. antipodarum sexual and asexual lineages would be a worthwhile next step in evaluating the potential that these candidate loci play a role in the transitions to asexuality.

These results are consistent with the apparent lack of systematic differences in gene expression between sexual and asexual reproductively active female P. antipodarum, including our “focal” gene set with GO terms associated with reproduction. We used a chi-square test to evaluate the proportion of differentially expressed genes in the “focal” vs. “control” set of genes from the sex-asex expression analysis. We found that the “focal” genes were less often differentially expressed than the random set of genes (χ² = 13.3601, p < 0.05). Of the “focal” genes, 0.92% were differentially expressed, compared to the 5.29% in the random set. The random set of genes had a very similar amount of differentially expressed genes detected transcriptome-wide (5.67%). This sampling of genes likely involved with reproduction and cell cycle processes is not exhaustive, but this evidence for relatively few differentially expressed genes from this large set of genes suggests that sexual and asexual P. antipodarum reproduction are not globally different at the level of gene expression. More rigorous evaluation of sexual and asexual gene regulation during reproduction would benefit in identifying genes in P. antipodarum that are upregulated during reproduction or exclusively expressed in ovaries.

Similar results were observed in comparisons of gene expression between sexual and asexual species of Artemia brine shrimp, where only 60 genes were differentially expressed in female gonads (Huylmans et al., 2021). Like brine shrimp, P. antipodarum may experience little reprogramming of reproductive pathways following transitions to asexuality, possibly because asexuality is recently derived in these species (Neiman et al., 2005; Rode et al., 2021) compared to much longer divergence times between sexual and asexual Timema species (Schwander et al., 2011), which have a high degree (∼8% of all genes) demonstrating convergent expression differences between reproductive modes (Parker et al., 2019). Extrapolation from asexual brine shrimp must be done with caution, however, because unlike (to the extent that we can discern) asexual P. antipodarum, asexual brine shrimp reproduce via automixis, in which meiosis occurs. Recent developments regarding the extent of cryptic sex occurring in asexual Artemia also raises the question of whether these lineages can truly be considered to be asexual (Boyer et al., 2021). The likelihood that apparently asexual lineages engage in rare sex is increasing in light of a growing body of evidence for such events, even extending to the bdelloid rotifers (Vakhrusheva et al., 2020; Laine et al., 2022; Paul et al., 2022), (in)famous as the most “ancient” of all asexuals. If long-term asexual persistence relies on occasional sex, then successful asexual lineages could plausibly represent taxa that can still engage in occasional recombination (or recombination-like) processes. To date, there is no evidence for recombinational processes in asexual P. antipodarum, though the occasional production of male offspring by asexual female P. antipodarum (Neiman et al., 2012) hints that this might be possible. On the flip side, that these males produce sperm that have strikingly aberrant morphology compared to sexual counterparts (Jalinsky et al., 2020) suggests that such recombination within asexual lineages might be extremely rare. Similarly, taxa thought to be obligate sexuals that give rise to asexuals may have cryptic asexual reproductive capabilities (i.e., facultative asexuality rather than obligate sexuality) and gradually transition to obligate asexuality, as had been hypothesized in Timema stick insects (Schwander and Crespi, 2009) and more recently supported with genetic evidence (Larose et al., 2022). There is nevertheless no evidence to date for facultative asexuality in P. antipodarum.

Using Blast2GO and Fisher’s Exact test (FDR q < 0.05), we identified 12 GO terms that were significantly overrepresented in the set of transcripts differentially expressed between sexual and asexual lineages (Figure 3C). We did not detect any statistical overrepresentation of GO terms related to reproduction, sex, meiosis/mitosis, chromosome movement, or any specific function we might expect to differ between diploids and triploids in sets of differentially expressed transcripts in sexual vs. asexual comparisons but we did identify 12 GO-terms enriched in differentially expressed genes. Several of these genes were annotated for functions related to metabolism, including the functions of “chitin-metabolic activity,” “glucosamine-containing compound metabolic process,” “amino sugar metabolic process,” and “aminoglycan metabolic process.” Chitin and aminoglycan metabolism were enriched GO terms for DE genes in cyclical vs. obligate parthenogenic rotifers (Hanson et al., 2013). Other RNA-seq studies in asexual vs. sexual animals have also reported metabolic GO terms enriched in sets of differentially expressed genes (Gallot et al., 2012; Liu et al., 2014; Parker et al., 2019). RNA sequencing has revealed evidence for distinct metabolic states, particularly with respect to amino acid metabolism, in sexual vs. asexual stages of flatworms (Sekii et al., 2019). Carter et al. (2012) found evidence for different energy investments between obligate and cyclical asexual aphids that indicate greater investment toward reproduction in obligate asexuals. Reminiscent of these findings, asexual female P. antipodarum grow more rapidly and reach reproductive maturity earlier than sexual females (Larkin et al., 2016). Our results may therefore be indicative of regulatory differences in metabolic processes in sexual vs. asexual P. antipodarum.

Polyploidy could also contribute to the differential expression observed between sexual and asexual P. antipodarum and enrichment in metabolic processes. Polyploidy and large genome size are also expected to influence standard metabolic rate via increasing cell size (Cavalier-Smith, 1978), as has been observed in polyploid fish (Maciak et al., 2011). It is unclear whether higher-ploidy P. antipodarum actually have larger somatic cells (Neiman et al., 2009), with the exception of sperm cells (Jalinsky et al., 2020). The overwhelming majority of literature on gene expression changes following polyploidization refer to allopolyploids, whereas whether autopolyploidy influences gene expression is not as well studied but is generally not associated with major shifts in gene expression (reviewed in Parisod et al., 2010). Long-term genome reorganization during diploidization following autopolyploidy, however, may be accompanied by significant changes in gene regulation while duplicated sequences diverge (e.g., Gundappa et al., 2022). Despite wide DNA content variation in P. antipodarum asexual polyploids (Neiman et al., 2011; Million et al., 2021), a trend toward DNA reduction in asexuals does not emerge (Million et al., 2021), indicating that polyploidy is maintained. The relatively recent derivation of most asexual lineages of P. antipodarum (Neiman et al., 2005), however, indicates that we might not yet be able to observe major global differences in gene expression between diploids and polyploids. Altogether, whether reproductive mode or ploidy variation is more important to gene regulation-related metabolic and physiological differences in P. antipodarum will need to be carefully evaluated in downstream studies comparing diploid, triploid, and tetraploid snails.

Along these lines, we cannot exclude the possibility that there exist differences in gene expression between sexual and asexual P. antipodarum that went undetected as a function of, for example, relatively small magnitude of differences or as a consequence of our focus on reproductively active females. Other studies have used RNA-seq to compare gene expression in reproductive tissue and other tissues/body regions in sexual and asexual lineages to distinguish the genes that may be most directly related to reproduction (Parker et al., 2019; Huylmans et al., 2021). This approach could also be useful to distinguish between ploidy and reproductive mode driving differences in gene expression. In future efforts to characterize the genetic causes of asexuality in P. antipodarum, improvements in microdissection of ovaries would be particularly useful to more likely capture tissue enriched with active oogenesis, as would characterization of critical time points in reproduction, and particularly, determining where meiosis predominantly occurs (i.e., certain regions of the ovaries).

The discovery of distinguishing features for egg production and development between sexual and asexual P. antipodarum might also provide a powerful means of establishing how to proceed. At this time, we are unaware of any differences in early development or egg formation between sexual and asexual P. antipodarum; profiling gene expression across development may be a useful avenue to explore. While we focused on likely protein-coding genes in this study, it is also possible that non-coding RNAs could play an important role in the differences in sexual and asexual reproduction. Indeed, a central role for non-coding RNAs has been demonstrated in Boechera (Mau et al., 2013), for which the non-coding RNA UPGRADE2 which appears to be responsible for unreduced pollen formation. Finally, we also cannot rule out the possibility that posttranscriptional regulation might be involved with sexual vs. asexual reproduction in P. antipodarum.