Repeated Evolution of Asexuality Involves Convergent Gene Expression Changes

Abstract Asexual reproduction has evolved repeatedly from sexual ancestors across a wide range of taxa. Whereas the costs and benefits associated with asexuality have received considerable attention, the molecular changes underpinning the evolution of asexual reproduction remain relatively unexplored. In particular, it is completely unknown whether the repeated evolution of asexual phenotypes involves similar molecular changes, as previous studies have focused on changes occurring in single lineages. Here, we investigate the extent of convergent gene expression changes across five independent transitions to asexuality in stick insects. We compared gene expression of asexual females to females of close sexual relatives in whole-bodies, reproductive tracts, and legs. We identified a striking amount of convergent gene expression change (up to 8% of genes), greatly exceeding that expected by chance. Convergent changes were also tissue-specific, and most likely driven by selection for functional changes. Genes showing convergent changes in the reproductive tract were associated with meiotic spindle formation and centrosome organization. These genes are particularly interesting as they can influence the production of unreduced eggs, a key barrier to asexual reproduction. Changes in legs and whole-bodies were likely involved in female sexual trait decay, with enrichment in terms such as sperm-storage and pigmentation. By identifying changes occurring across multiple independent transitions to asexuality, our results provide a rare insight into the molecular basis of asexual phenotypes and suggest that the evolutionary path to asexuality is highly constrained, requiring repeated changes to the same key genes.


Introduction
Sexual reproduction is extremely costly. Sex is less efficient than asexuality for transmitting genes to future generations (Maynard Smith 1971) and in order to outcross, an individual has to find a partner, forgo foraging, and risk contracting sexually transmitted diseases and predation while mating (Bell 1982;Lehtonen et al. 2012). Yet, the overwhelming number of sexual, as compared with asexual, animal and plant species (Avise 2008;van der Kooi et al. 2017) indicates that sexual reproduction is highly advantageous. Identifying potential advantages conferred by sex has motivated decades of research and a rich body of work on the evolution and maintenance of sexual and asexual reproduction has been produced (reviewed in Bell 1982;Lewis 1987;West et al. 1999;Otto 2009;Neiman et al. 2017). In contrast, little is known about the molecular underpinnings required to evolve asexuality from sexual ancestors (Neiman et al. 2014). Yet, these molecular underpinnings have the potential to provide insights into the processes involved in the evolution of asexuality, and to help understand how sex is maintained. For example, sex is more easily maintained if asexuality evolves gradually in a sexual population than if it emerges suddenly via major effect mutations (Templeton 1982;Burt 2000;Schwander et al. 2010).
Some insight into the genetic basis of asexuality has been gained from studies of individual asexual lineages (Innes and Hebert 1988;Lynch et al. 2008;Eads et al. 2012;Jaqui ery et al. 2014), but a broad comparative framework for exploring common principles of the molecular basis of asexuality is lacking. For example, a major unresolved question is whether independent transitions to asexuality involve similar or different molecular changes. To address these shortcomings, we explored the molecular underpinnings of asexuality in stick insects of the genus Timema, a genus of wingless, herbivorous insects native to the West coast of North America and the mountains of the Desert Southwest. This group is uniquely suited for comparative studies of asexuality, as asexuality has evolved at least seven times independently (Schwander et al. 2011; fig. 1), allowing us to study convergence across replicate transitions from sexual to asexual reproduction. Furthermore, close sexual relatives are at hand for each asexual lineage for comparison. All asexual Timema species reproduce via obligate parthenogenesis (Schwander and Crespi 2009), meaning that they evolved the ability to produce unreduced eggs which develop without fertilization by sperm. Additional phenotypic changes evolved convergently as adaptations to asexuality, including a reduced sperm storage organ, and reduced sexual pheromone production . Thus, asexual Timema females are less attractive to sexual males , which use both airborne and contact signals to identify suitable mates (Nosil et al. 2007;Arbuthnott and Crespi 2009;Schwander, Arbuthnott, et al. 2013), and even when copulations between sexual asexual females and males from sister-species are forced under laboratory conditions, eggs are not fertilized .
To capture molecular changes associated with the evolution of asexuality we performed whole-body and tissuespecific transcriptome sequencing (RNA-seq) on females from five sexual and five asexual Timema species ( fig. 1). We chose two different tissues, the reproductive tract and legs, to identify the molecular mechanisms underlying the production of asexual offspring (reproductive tract), and adaptations to a celibate life (e.g., reduction of various different sexual traits in the reproductive tract and legs). Note that the reproductive tract and leg samples actually represent a collection of tissues, but we refer to them as tissues throughout for brevity. Whole-body samples were included as they allow us to identify important changes that may be missing in the tissue-specific transcriptomes. Using this approach, we identified convergent expression changes which were likely driven by selection. We also observed changes specific to each sexual-asexual species-pair which typically showed concerted changes across tissues, consistent with being a product of drift (Blekhman et al. 2008;Liang et al. 2018). Finally, to complement our expression analyses, we examined patterns of molecular evolution in genes showing convergent expression changes following a transition to asexuality.

Transcriptomes and Orthology
Reference transcriptome assemblies for each species were generated previously (Bast et al. 2018). Bast et al. (2018) also identified 3,010 one-to-one orthologs, which were used as our transcriptome reference. For each tissue, orthologs with low expression (counts per million <0.5 in two or more libraries per species) were filtered prior to expression analyses. Thus, the final number of orthologs kept for analyses of whole-body, reproductive tract, and leg samples was 2,984, 2,753, and 2,740, respectively.

Convergent Gene Expression Changes
We identified convergent gene expression changes between sexual and asexual species by modelling gene expression as a function of species-pair (see fig. 1), reproductive mode (sexual or asexual), and their interaction in edgeR (Robinson et al. 2010). In such a model, convergence is indicated by an overall effect of reproductive mode (FDR < 0.05), but no interaction (FDR > 0.05; supplementary table 1, Supplementary Material online). Approximately four times as many genes changed convergently in the reproductive tract (7%; 203/2,754) and legs (8%; 206/2,737) as compared with the whole-body (2%; 57/2,985), perhaps reflecting the relative difficulty in identifying expression changes in complex tissue assemblies such as whole-bodies (Johnson et al. 2013 Liang et al. 2018). If such changes are nonneutral they are likely to be deleterious due to pleiotropy. As a consequence, gene expression changes that occur in parallel across different tissues are more likely a result of drift rather than selection (Blekhman et al. 2008;Liang et al. 2018). We found that convergent changes between sexual and asexual species were highly tissue-specific. Only 22 of the convergent genes in the reproductive tract (203) and legs (206) overlapped between the two tissues, a value not significantly greater than expected by chance (table 1). There was also little overlap between convergent genes in the two tissues and whole-bodies (table 1, supplementary fig. 2, Supplementary Material online). This pattern is consistent with the idea that convergent expression changes in asexuals are driven by selection rather than by drift. This interpretation of selection also predicts that convergent genes will be involved in divergent functions in each tissue which, as we show below, is indeed the case.
To directly test whether the convergent changes are driven by selection, we modeled convergent change as an Ornstein-Uhlenbeck (OU) process. Selection-driven expression changes are identified in this framework by comparing the likelihood of two models: a drift-model where expression changes are modelled as simple Brownian motion, and an adaptiveoptima model where expression changes are modelled as the result of both Brownian motion and selection toward specified adaptive-optima (see Butler and King 2004;Cressler et al. 2015). In this case, the second model specified that asexual species had a different adaptive-optimum than sexual species. As each asexual species is phylogenetically independent, this tests for convergence of expression. Using this approach, we found that the adaptive-optima model had a significantly better fit than the drift model for 70-80% of convergently expressed genes, depending on tissue (supplementary tables 2 and 3, Supplementary Material online, FDR < 0.05). This strengthens our interpretation that the majority of the convergent changes between sexual and asexual species is driven by selection.
Finally, we compared the two-state adaptive-optima model, where all asexual species share the same optimum, to a model where the optimum can vary between asexual

Functional Processes of Convergently Expressed Genes
To detect convergence at the process level, we performed gene set enrichment analyses (GSEA) for each tissue separately. Briefly, we scored Gene Ontology (GO) terms according to the rank of convergent expression change of genes annotated to the terms; GO terms were then called significant if they had a better average rank than expected by chance (see Materials and Methods). For each tissue, >100 GO terms are enriched (P < 0.05), providing strong support for convergence of biological processes between asexual species (supplementary tables 4-6, Supplementary Material online). This signal is not dependent on any threshold at the gene level, and thus provides information on convergence at the process level due to small but consistent contributions from many genes. Consistent with the gene expression results, enriched GO terms were generally tissue specific; we found no significant overlap between GO-terms enriched in the legs and reproductive tract (11 shared terms, FDR ¼ 0.123), between the legs and whole-body (4 shared terms, FDR ¼ 0.799), or between whole-body and reproductive tract samples (ten shared terms, FDR ¼ 0.064).
To reduce the number of enriched GO terms to examine we semantically clustered enriched GO terms using ReviGO (Supek et al. 2011; supplementary tables 7-9, Supplementary Material online). The annotations of convergent changes in the reproductive tract reflect the convergent evolution of parthenogenesis in asexual Timema, as they were linked to meiosis (meiotic spindle organization, meiosis II, centrosome duplication, meiosis I cytokinesis, meiosis II cytokinesis), and reproduction (growth of a germarium-derived egg chamber, sperm individualization, gamete generation). However, convergent changes were also linked to neuron development (neurogenesis, neuron development, neuron recognition), as well as several GO terms involved in development and metabolic processes for which the link to asexuality is less clear. In legs we identified GO terms involved in immune defense (response to fungus, regulation of production of molecular mediator of immune response, regulation of antimicrobial peptide production, regulation of humoral immune response), which may be because asexual females are no longer susceptible to the costs associated with diseases transmitted from sexual interactions (Knell and Webberley 2004). Convergent changes were also linked to sex determination (primary sex determination; soma, primary response to X:A ratio), which may control changes in the expression of sexual traits, and several metabolic processes. In whole-body samples we find some reproduction associated terms (courtship behavior, male mating behavior, male courtship behavior, sperm storage, regulation of ovulation) as in the reproductive tract, and behavioral, and immune related terms (immune response-regulating cell surface receptor signaling pathway) as in legs, but also some unique terms relating to the cuticle (ecdysone, pupal chitin-based cuticle development).

Convergently Expressed Genes in Whole-bodies Show Evidence for Sexual Trait Decay
Several of the enriched functional processes described above are suggestive of sexual trait decay. Under this scenario we expect a reduction of purifying selection on genes underlying sexually dimorphic traits in asexual species, indicated by an increased accumulation of nonsynonymous changes.
The power to detect differences in pN/pS or dN/dS between gene sets in asexuals is low, as genes are inherited as a single linkage group. Nevertheless, we found that genes showing convergent changes in expression in whole-bodies showed elevated pN/pS and dN/dS when compared with the genomic background (permuted t-test P-value for: pN/pS < 0.0001, dN/dS ¼ 0.0084, fig. 3, supplementary figs. 4 and 5, Supplementary Material online), consistent with the idea of sexual trait decay. Sexual trait decay is further We conducted several additional analyses to check the robustness of our results and corroborate our interpretations. Firstly, we examined in detail the associated functional annotations of candidate gene sets for which there was very strong evidence for convergent changes, and secondly, we used cross-species mapping to examine expression changes occurring across the whole transcriptome, rather than only in the subset of genes we identified as single copy orthologs between the ten species. Both approaches support the results from our original analyses and are described below.

Strongly Convergent Candidate Genes and Their Function
Although all the convergent genes we identified showed an overall shift in expression across the five species-pairs, often expression change in one or two of the pairs was small (<1.2fold change). We defined top candidate genes as convergent genes for which the absolute log 2 fold change in expression was >0.25 ($1.2-fold change) for all species-pairs. Most of these top genes showed convergent shifts in the reproductive tract (36 genes, relative to 4 and 15 genes for whole-body and legs, respectively; fig. 4   MBE Material online). The functions of these candidate genes largely reflected the functional processes identified for the full set of convergently expressed genes, and highlight a number of key genes potentially involved in producing asexual phenotypes. For the reproductive tract four genes are involved in meiotic spindle formation and centrosome organization (OG-513, OG-1448, OG-1488, OG-314). In particular we find two genes (OG-1448, OG-1488) belonging to a family of Elovl proteins that mediate elongation of very-long-chain fatty acids, including an ortholog to Drosophila melanogaster gene bond, which effects spindle formation and has been shown to be important for meiotic, but not mitotic, cytokinesis (Szafer-Glusman et al. 2008). In particular, D. melanogaster males defective for bond commonly display two to four nuclei in spermatids causing sterility. Female bond mutants are also infertile (Szafer-Glusman et al. 2008), although the mechanism is unknown. The other two genes (OG-513, OG-314) have roles in centrosome function, including an ortholog to poc1 which is involved in centrosome formation (Blachon et al. 2009). Six genes (OG-758, OG-2002, OG-1478, OG-1993, OG-2686 were annotated with reproduction associated terms which may be responsible for the convergent reproductive changes we observe between asexual and sexual females. Interestingly, one gene, OG-511, is an ortholog to glucose dehydrogenase which is important for sperm storage in female D. melanogaster (Bloch Qazi et al. 2003). Finally, we find that 11 genes (OG-1195, OG-1478, OG-1841, OG-2197, OG-2808, OG-366, OG-445, OG-511, OG-705, OG-712, OG-758, OG-810) have annotations to the nervous system. The majority of these appear to be sensory in nature, and in particular seven are annotated with the Convergent Gene Expression Changes . doi:10.1093/molbev/msy217 MBE GO term "sensory perception of pain." Changes in these genes may represent changes associated with female receptivity and postmating behavior in asexual females, which are targets of substances in the male ejaculate (Heifetz and Wolfner 2004;Sakai et al. 2009;Heifetz et al. 2014). For leg samples three genes (OG-1651, OG-2048, and OG-1081) are involved in immune defense. In particular orthologs of both genes (Trx-2 and MP1) are involved in the activation of melanization in response to fungal and bacterial infection in D. melanogaster (Tang et al. 2006;Jin et al. 2008). Three genes are involved in cuticle development (OG-2221, OG-2738, and OG-2995. Orthologs of two other genes (OG-1371 and OG-2031) are involved in male specific behaviors (male courtship behavior and intermale aggressive behavior) in D. melanogaster (CaMKII and Fkbp14; Mehren and Griffith 2006;Edwards et al. 2009). Since these genes are also expressed in females, changes to their expression may have resulted from the release of intralocus sexual conflict.
Whole-body samples had only four strong candidate genes, and all either have no annotation or have only broad GO-terms annotated. One potentially interesting gene, OG-2188, has an ortholog (CG12237) that has been associated with female sterility in D. melanogaster (Sopko et al. 2014). Finally, the remaining candidate genes across all tissues were either unannotated (12 genes) or only have very broad GOterms annotated (ten genes).

Cross-species Mapping
Using only the 3,010 genes with 1-to-1 orthologs across all species could impact our ability to detect convergent changes since we only use a relatively small fraction of the total number of transcripts in each assembly (23,435-37,847; supplementary table 10, Supplementary Material online). To investigate more genes, we mapped reads from all samples to genes from each species which had a reciprocal-best-blasthit between sexual-asexual sister species (which includes the 1-to-1 orthologs analyzed above). This approach generated ten different data sets (one for each species assembly), with between 15,500 and 17,583 genes. After filtering out genes with low expression (using cpm, see Materials and Methods) in each data set, this approach allowed us to examine between 2.43 and 3.12 (dependent on species and tissue) times more genes than using the 1-to-1 orthologs (supplementary table 10, Supplementary Material online). Results from this approach qualitatively confirmed the results found using only the 1-to-1 orthologs: the percentage of genes showing a convergent expression ranged from 4% to 5% for whole-body samples and 6-8% for the reproductive tract and leg samples, dependent on which of the species transcriptome was used (supplementary table 10, Supplementary Material online), and GSEA produced similar enriched GO terms (supplementary tables 11-13, Supplementary Material online).

Species-pair Specific Changes
The approach taken above allowed us to identify genes which showed convergent changes in expression across independent transitions to asexuality. This approach will not identify expression changes confined to a single or a few species-pairs. Changes occurring in only a minority of species-pairs are clearly not convergent at the gene expression level; however, these changes could be convergent at the functional process level, whereby species-pair specific changes in gene expression are involved in common functional processes between species-pairs (Rittschof et al. 2014;Berens et al. 2015). To test this, we compared each asexual species to its closest sexual relative and called differentially expressed (DE) genes from each pairwise comparison.
The number of significantly DE genes between each pair varied greatly depending on species-pair and tissue, with a generally greater number of genes DE in leg tissue (59-626, supplementary fig. 6, Supplementary Material online). This greater number in leg tissue is likely due to the smaller variation between replicates (common biological coefficient of variation was lowest for legs: whole-body ¼ 0.314, reproductive tract ¼ 0.340, legs ¼ 0.238), as tissue differences disappeared when a fold-change threshold was applied (supplementary fig. 7, Supplementary Material online). There were no genes that showed overlap between all sexual-asexual species-pairs in any tissue ( fig. 5A). Examination of overlaps between pairs of sexual-asexual species-pairs found some overlapping genes, but these were close to the expectation by chance ( fig. 5B, for all levels see supplementary table 14, Supplementary Material online). The majority of the DE genes also showed a significant interaction between speciespair and reproductive mode in the model used to identify convergently changing genes (whole-body ¼ 69%, reproductive tract ¼ 66%, and legs ¼ 81%), corroborating the finding that the vast majority of the DE genes is species-pair specific. Note the species-pair by reproductive mode interactions do not appear to be generated by one specific species-pair as generally genes DE between one species-pair were not DE between the other four species-pairs. The DE genes of the different species pairs are not involved in convergent functional processes. Species-pair-specific genes were enriched for a number of GO terms (presented in supplementary table 15, Supplementary Material online); however, no GO terms were found to overlap between all pairs in any tissue (supplementary fig. 8, Supplementary Material online). Examination of overlaps between pairs of sexual-asexual species-pairs found some overlapping GO terms, but these were close to the expectation by chance (supplementary table 16, Supplementary Material online). This pattern remained even when a more liberal approach, whereby related GO-terms were considered as a unit, was applied (supplementary fig. 9A, Supplementary Material online). This overall lack of overlap suggests that species-pairspecific genes are not involved in producing convergent phenotypes. Instead, species-pair-specific genes are the product either of lineage-specific selection (for instance, if different species use different mechanisms to achieve parthenogenesis) or of drift. These two processes are difficult to disentangle, but our results are more consistent with drift rather than with lineage-specific selection. Indeed, species-pair specific genes showed similar changes in gene expression across tissues, in contrast to the mainly tissue-specific changes uncovered for convergently changing genes. The overlap of species-pair Parker et al. . doi:10.1093/molbev/msy217 MBE specific genes between tissues was significantly greater than expected by chance (table 1).
Finally, these results were reproduced when examining a much larger set of genes (genes with reciprocal-best-blast-hits between species-pairs, see above) as both genes DE between each species-pair, and their enriched GO terms, showed little overlap (supplementary figs. 10 and 11, and supplementary tables 17 and 18, Supplementary Material online).

Discussion
Asexuality has convergently evolved numerous times across the tree of life, and a large body of research focuses on the reasons why sexual reproduction persists in the face of competition from asexual lineages. In contrast, the molecular underpinnings of transitions from sexual to asexual reproduction remain largely unknown (Neiman et al. 2014). In this study, we examined gene expression changes associated with transitions to asexuality across five independently evolved asexual lineages, in whole-bodies, reproductive tracts and legs. The changes we observe provide, for the first time, insights into the convergent evolution of asexuality at the molecular level.
We found evidence for convergent changes in gene expression in all three tissues. Four lines of evidence suggest that these changes are a product of selection. Firstly, parallel  Convergent Gene Expression Changes . doi:10.1093/molbev/msy217 MBE changes across multiple independent transitions represent strong evidence of selection and thus are unlikely to be due to drift (Losos 2011;Natarajan et al. 2016). Secondly, convergent changes were primarily tissue-specific. This finding is consistent with selection, because expression changes due to drift are likely to be correlated across tissues (Blekhman et al. 2008;Liang et al. 2018). Indeed, the different functional roles of reproductive tracts and legs make it unlikely that selection would drive changes in the same genes in all tissues. Thirdly, expression changes of convergent genes inferred over the Timema phylogeny fit better with a model of adaptive change than with models of simple drift or of independent change in asexuals. Finally, the functional processes of convergently expressed genes mirror the changes observed at the phenotypic level, supporting the interpretation that these genes contribute to the convergent phenotypic changes observed between sexual and asexual females.
The overall amount of convergence is striking, particularly in the reproductive tract and legs with $8% of genes showing a convergent shift in expression. Such a large amount of convergence suggests that the path from sexual reproduction to asexuality is strongly constrained, requiring changes to the same genes and biological processes in order to produce asexual phenotypes.

Convergent Changes in Gene Expression Reveal the Mechanisms Underlying the Production of Asexual Offspring
Asexuality is a complex adaptation that includes two major components: the ability to produce viable asexual offspring, and secondary adaptive changes that would not have been selected for in sexual species (e.g., the reduction of costly sexual traits). A key change necessary for the production of asexual offspring is the ability to produce unreduced eggs (Engelst€ adter 2008). Convergently expressed genes in the reproductive tract were enriched for changes in meiosis, and in particular meiotic spindles, which are key for the proper division of cells during meiosis. Mutations in meiotic spindles have been shown to result in unreduced meiotic products in D. melanogaster, and specifically in two genes (bond, Szafer-Glusman et al. 2008 andpelo, Eberhart andWasserman 1995) which show convergent changes in expression in asexual Timema. As such we suggest that these changes may underlie the nonreduction of eggs in asexual Timema. An alternative hypothesis is that since Timema reproduce parthenogenetically (and thus likely no longer recombine) changes in meiotic genes represent trait decay. Although possible, previous work has shown that, in fact, meiotic genes are not only retained in asexual lineages without damaging mutations, but often appear to be subject to selection for changes in expression, via duplication or differential upregulation of promoters (Srinivasan et al. 2010;Hanson et al. 2013;Raborn et al. 2016;Warren et al. 2018). Taken together with our results, we suggest that modifications to meiotic genes, specifically those that that disrupt meiotic cell division, are key in overcoming a major barrier to the evolution of asexuality: the production of unreduced eggs.
The production of unreduced eggs is not the only barrier to producing offspring asexually. In most species, sperm transfer essential components for the formation of a functioning centrosome (Schatten 1994;Manandhar et al. 2005). This paternal contribution represents a second key barrier in the evolution of parthenogenesis in many systems (Engelst€ adter 2008). However, in phasmids the centrosome is assembled without any contribution from sperm in both sexual and asexual species (Marescalchi et al. 2002). This may act as pre-adaptation for asexuality in stick insects and account, in part, for the large number of asexual stick insect species.
A final barrier to asexual offspring production in many systems is egg activation. In many species mature oocytes are arrested at a specific stage (e.g., at metaphase II in mammals, and metaphase I in most insects), and must be activated by sperm to re-enter the cell-cycle (Stricker 1999;Engelst€ adter 2008;Nishiyama et al. 2010). In insects however, egg activation does not require sperm as activation is induced by the transit through the reproductive tract (Sartain and Wolfner 2013). Despite this, ovulation and egg-laying rates are strongly tied to mating (Eberhard 1996;Gillott 2003) meaning this signal must be modified in order for asexual insects to have normal levels of fecundity. In insects, the signal to a female that she has successfully mated is likely detected by sensory neurons in her reproductive tract (Yapici et al. 2008). Consistent with this, we find changes in gene expression linked to sensory neurons in the reproductive tract of asexual females, which may act to cue high levels of ovulation without mating. Alternatively, these changes may represent the decay of these neurons since they are no longer needed to detect mating events, or these changes may result from cessation of sexual conflict. Sensory neurons in the reproductive tract are known targets of substances in the male ejaculate (Heifetz and Wolfner 2004;Sakai et al. 2009;Heifetz et al. 2014) to induce the release of eggs and to reduce female receptivity (Gillott 2003). This manipulation is countered by female resistance adaptations which are likely costly, meaning that, following a transition to asexuality, there will be selection against them.

Convergent Changes in Gene Expression Show Evidence for the Decay of Female Sexual Traits
Sexual traits in asexual females are often observed to be reduced or lost (van der Kooi and Schwander 2014). For instance, in insects, females typically produce pheromones as a sexual cue to attract males (Greenfield 2002), and this cue has been repeatedly reduced or lost in several asexual species (see van der Kooi and Schwander 2014), including Timema . Such trait decay can be the result of either reduced purifying selection acting on traits that are now selectively neutral, or selection to reduce the cost of producing sexual traits. In asexual Timema reproductive decay has been primarily attributed to selection rather than reduced purifying selection, as reproductive trait decay in very young asexual lineages is as extensive as in old ones .
Convergent gene expression changes underlying the decay of reproductive traits are mostly observed in Timema Parker et al. . doi:10.1093/molbev/msy217 MBE whole-bodies. In particular, we find enrichment of terms associated with sperm storage and sexual behavior. Changes in the legs were less obviously associated with reproductive trait decay, however we do find changes in genes involved in cuticle development, pigment biosynthesis, sensory perception of touch, and changes in sexual behavior. These changes could represent the reproductive decay of both sexual cues (e.g., cuticular hydrocarbons and pigmentation which are both important for mate choice in insects [reviewed in Hunt and Sakaluk 2014]), and their detection (via sensory receptors on the leg [reviewed in Dahanukar et al. 2005]). In addition, we also find changes in genes associated with sex determination in the soma, including sexlethal, a master-feminizing switch in Drosophila (Cline et al. 2010) which may have a major influence on the development of many sexual traits in the legs.
Although we focus on expression, it is possible that the decay of sexual traits is also evident at the sequence level. By examining the coding regions of genes, we found evidence for reduced purifying selection acting on the sequence of genes showing convergent expression changes in the whole-body. This suggests, that in some cases, the reduction of sexual traits may be accomplished by both expression and sequence changes, which potentially act interactively to produce a phenotypic change.
Unexpectedly, we also find changes to immune function in the legs and whole-body, the majority of which show downregulation in asexual females. A possible explanation for this is that asexual females are likely to face a reduced number of immune challenges compared with sexual females due to the elimination of sexually transmitted diseases, the costs of which can be considerable, even shaping the evolution of many aspects of an organism's life history, such as mate choice, mating rate, and sexual signal investment (Kokko et al. 2002;Knell and Webberley 2004). As such we suggest asexual females may be reducing the allocation of resources to immune function due to the absence sexually transmitted diseases. This effect may be particularly strong in solitary species such as Timema, where the majority of socially transmitted diseases comes from sexual interactions.

Species-pair Specific Changes
In addition to convergent changes, we also identified many species-pair specific gene expression changes. In contrast to convergent genes, species-pair specific genes showed common shifts in expression across tissues, and inconsistent associations with functional processes between speciespairs, that were largely unrelated to asexual phenotypes. Taken together, these results suggest that the majority of changes we observe from a single sex-asex species-pair comparison is due to drift rather than selection. Our findings thus highlight the problem of drawing inferences on the causes or consequences of asexuality from the examination of only a single transition to asexuality, whereas examining several transitions allows us to disentangle adaptive changes and those due to drift.
Overall, we find evidence for a striking number of convergent changes across five transitions to asexuality. Previous studies that have examined convergent expression changes across the genome have found little evidence of convergence at the gene level (<1%; Rittschof et al. 2014;Berens et al. 2015), underlining the surprisingly large amount of convergent gene expression (up to 8%) we find here. The amount of molecular convergence to expect, however, is dependent on several factors including the complexity of the phenotype, and the size of the mutational target (Martin and Orgogozo 2013). For example, we find that a key change required for asexual reproduction, the production of unreduced eggs, likely requires changes to meiotic spindle regulation. The pathways that govern meiotic spindle regulation are relatively small in number (Bennabi et al. 2016), meaning that only a small minority of genes are likely able to confer the relevant changes, making the chance of molecular convergence for this trait relatively high. In contrast, the observed reduction of sexual traits could be produced by changes to numerous genes and pathways (i.e., there is a large mutational target) making convergent molecular changes for these traits less likely. Despite this, our and previous studies examining trait loss have also demonstrated a high amount of convergence (Whittall et al. 2006;Gross et al. 2009;Partha et al. 2017;Sackton et al. 2018), implying that certain genes have a disproportionate role in not only the convergent evolution of novel phenotypes, but also in their convergent loss (Martin and Orgogozo 2013;Stern 2013).

Samples
Females for whole-body samples were collected from the field as juveniles in spring 2013. All individuals were then raised in common garden conditions (23 C, 12h:12 h, 60% humidity, fed with Ceanothus cuttings) until eight days following their final molt. Prior to RNA extraction, individuals were fed with artificial medium for two days to avoid RNA contamination with gut content and then frozen at À80 C. Individuals used for tissue-specific samples were collected in spring 2014 as juveniles and raised in the same common-garden conditions as whole-body samples. For leg samples three legs were used from each individual (one foreleg, one midleg, and one hindleg). Reproductive tracts were dissected to consist of ovaries, oviducts and spermatheca. Note the same individuals were used for leg and reproductive tract samples. Collection locations for all samples are given in supplementary table 19, Supplementary Material online.

RNA Extraction and Sequencing
The three biological replicates per species and tissue consisted of 1-9 individuals per replicate, which were combined prior to RNA extraction (207 individuals in 90 replicates in total; see supplementary table 19, Supplementary Material online). RNA extraction was performed by freezing samples in liquid nitrogen followed by addition of Trizol (Life Technologies) before being homogenized using mechanical beads (Sigmund Lindner). Chloroform and ethanol were then added to the samples and the aqueous layer transferred to RNeasy MinElute Columns (Qiagen). RNA extraction was then Convergent Gene Expression Changes . doi:10.1093/molbev/msy217 MBE completed using a RNeasy Mini Kit following the manufacturer's instructions. RNA quantity and quality was measured using NanoDrop (Thermo Scientific) and Bioanalyzer (Agilent). Strand-specific library preparation (one library per replicate) and single-end sequencing (100 bp, HiSeq2000) were performed at the Lausanne Genomic Technologies Facility.
The 90 libraries produced a total of just over 3 billion singleend reads. Four whole-body and six tissue-specific libraries produced significantly more reads than the average for the other samples. To reduce any influence of this on downstream analyses, these libraries were sampled down to approximately the average number of reads for whole-body or tissue-specific libraries respectively using seqtk (https://github.com/lh3/seqtk Version: 1.2-r94; last accessed December 9, 2018).

Transcriptome References
De novo reference transcriptome assemblies for each species were generated previously (Bast et al. 2018). For our analyses, we used the 3,010 one-to-one orthologs present in all ten Timema species as identified by Bast et al. (2018). Identified ortholog sequences varied in length among different species. Since length variation might influence estimates of gene expression, we aligned orthologous sequences using PRANK (v.100802, default options;Löytynoja and Goldman 2005) and trimmed them using alignment_trimmer.py (Parker 2016) to remove overhanging gaps at the ends of the alignments. If the alignment contained a gap of >3 bases then sequence preceding or following the alignment gap (whichever was shortest) was discarded. Three genes were discarded at this stage as the trimmed length of sequence was <300 bp. These trimmed sequences were then used as reference transcriptomes for read mapping. Note that genes with significant BLAST hits to rRNA sequences were removed from the transcriptome references prior to mapping.

Read Trimming and Mapping
Raw reads were trimmed before mapping. Firstly CutAdapt (Martin 2011) was used to trim adapter sequences from the reads. Reads were then quality trimmed using Trimmomatic v 0.36 (Bolger et al. 2014): first clipping leading or trailing bases with a phred score of <10 from the read, before using a sliding window from the 5 0 end to clip the read if four consecutive bases had an average phred score of <20. Following quality trimming any reads <80 bp in length were discarded. Quality-trimmed reads from each library were then mapped separately to the reference transcriptome using Kallisto (v. 0.43.1;Bray et al. 2016) with the following options -l 210 -s 25 -bias -rf-stranded for whole-body samples and -l 370 -s 25bias -rf-stranded for tissue specific samples (the -l option was different for whole-body and tissue specific samples as the fragment length for these libraries was different).

Differential Expression Analysis
Expression analyses were performed using the Bioconductor package EdgeR (v. 3.18.1; Robinson et al. 2010) in R (v. 3.4.1; R Core Team 2017). Analyses were done separately for each tissue. Genes with counts per million <0.5 in 2 or more libraries per species were excluded from expression analyses. Normalization factors for each library were computed using the TMM method in EdgeR. To estimate dispersion, we then fit a generalized linear model (GLM) with negative binomial distribution with the terms species-pair, reproductive mode and their interaction. We used a GLM likelihood ratio test to determine significance of model terms for each gene by comparing appropriate model contrasts. P-values were corrected for multiple tests using Benjamini and Hochberg's algorithm (Benjamini and Hochberg 1995), with statistical significance set to 5%. Using this approach, we classified genes as convergently DE when there was a significant effect of reproductive mode (FDR < 0.05) but no interaction effect of speciespair by reproductive mode (FDR > 0.05). DE genes within each species-pair were identified using pairwise contrasts between each sexual and asexual pair.
To determine if genes DE within each species-pair and tissue show greater than expected number of overlapping genes we used the SuperExactTest package (v. 0.99.4; Wang et al. 2015) in R which calculates the probability of multi-set intersections. When examining multiple intersections Pvalues were multiple test corrected using Benjamini and Hochberg's algorithm implemented in R.
To test if the observed number of convergent genes was significantly greater than expected by chance we performed a permutation test where, for the read counts of each gene, we randomly switched the assignment of reproductive mode (sexual or asexual) within a species-pair. Note that all biological replicates from a particular group were always assigned to the same reproductive mode (i.e., in the event of a switch, all sexual replicates were assigned as asexual, and vice versa). This process was repeated to produce 10,000 permuted data sets, which were then run through the gene expression pipeline described above to generate a distribution of the number of convergent genes we expect to find by chance.

OU Based Models
OU based models were fit using the R package ouch (v. 2.11-1; King and Butler 2009). We fit three models to each convergently expressed gene: a drift-model where changes are modelled as Brownian motion, a two-state adaptive-optima model where asexual species had a different adaptiveoptima than sexual species, and a multi-asexual optima model where sexual species and each asexual species had different optima (supplementary fig. 12, Supplementary Material online). Initial values of sqrt.alpha and sigma for the adaptive-optima models were set to 1. Log-likelihood ratio tests were used to determine if the two-state adaptiveoptima model was a significantly better fit to observed expression values than the drift-model, and if the multi-asexual optima model was a significantly better fit to observed expression values than the two-state adaptive-optima model. Pvalues were corrected for multiple tests using Benjamini and Hochberg's algorithm (Benjamini and Hochberg 1995), with statistical significance set to 5%. Mean observed expression values (log 2 counts per million) for each species and gene were calculated using EdgeR (see above). We used RAxML (v. 8.2.8, options: -p 12345 -m GTRCAT -T 6; Stamatakis 2014) Parker et al. . doi:10.1093/molbev/msy217 MBE to produce a maximum likelihood tree for use in ouch using the concatenated one-to-one ortholog alignments produced in Bast et al. (2018).

GO Term Analysis
Genes were functionally annotated using Blast2GO (version 4.1.9; Götz et al. 2008) as follows: sequences from each sexual species were compared with BlastX to either NCBI's nrarthropod or D. melanogaster (drosoph) databases, keeping the top 20 hits with e-values <1 Â 10 À3 . Interproscan (default settings within Blast2GO) was then run for each sequence, and the results merged with the BLAST results to obtain GO terms. This produced two sets of functional annotations, one derived from all arthropods and one specifically from D. melanogaster. The D. melanogaster GO term annotation generated around four times more annotations per sequence than NCBI's nr-arthropod database. We therefore conducted all subsequent analyses using the GO terms derived from D. melanogaster but note that results using the annotations from all arthropods were qualitatively the same (see supplementary fig. 9B, Supplementary Material online).
We conducted GSEA using the R package TopGO (v. 2.28.0; Alexa and Rahnenfuhrer 2016) using the elim algorithm to account for the GO topology. GSEA identify enriched GO terms in a threshold-free way, by finding GOterms that are overrepresented at the top of a ranked list of genes. For comparisons within a species-pair, genes were ranked by FDR; to identify enrichment of convergent genes, genes were ranked by FDR value for reproductive mode, with the FDR value for genes that showed a significant lineage by reproductive mode set to 1. GO terms were considered to be significantly enriched when P < 0.05. Enriched GO terms were then semantically clustered using ReviGO (Supek et al. 2011) to aid interpretation.
The significance of overlapping GO terms was determined using SuperExactTest as described above. The hierarchical nature of GO terms generates a bias towards finding a significant amount of overlap, since enrichment terms are nonindependent. It is however possible that the complexity of the GO term hierarchy could lead to convergent functional processes being overlooked. For instance, if a GO term is enriched in one comparison, but its parent term is enriched in another comparison, then there would be no apparent overlap. To address this, we also looked at the amount of 'linked overlap' of GO terms, whereby significant GO terms were first clustered together based on parent or child terms.
For the GO term enrichment analyses of convergently DE genes we used only the annotation from Timema bartmani as it had the most number of sequences annotated. Annotations to each of the other species were very similar to those from T. bartmani, with 80% of annotations being identical across all five species annotations. The remaining 20% of sequences were typically characterized by an additional term in one or more of the species. For comparisons within a lineage we used the annotation of the sexual species in that lineage. Although the annotations are very similar across all ten species the small differences in annotation could create differences in the amount of overlap observed between contrasts (e.g., if a term is annotated to an ortholog in one annotation but not another). To examine this, we repeated the analysis using only annotations from T. bartmani. This produced a virtually identical result (supplementary fig. 9C, Supplementary Material online) as when using the speciespair specific annotations.

Polymorphism and Divergence
To test for differences in the rate of evolutionary divergence between gene categories, we used dN/dS ratios for each of the one-to-one orthologs from Bast et al. (2018). Briefly, Bast et al. (2018) aligned each of the one-to-one orthologs with M-Coffee (Wallace et al. 2006) and trimmed them with Gblocks (Talavera and Castresana 2007). These alignments were then used as input for codeml of the PAML package (Yang 2007) to generate maximum likelihood estimates of dN/dS for each terminal branch in the phylogeny (using the "free model"). To obtain an estimate for pN/pS for each ortholog, reads from the whole-body libraries (based on a pool of three individuals) for each asexual species were mapped to the reference using RSEM/bowtie2 with default parameters and fragment length mean ¼ 200 fragment length sd ¼ 100 (Li and Dewey 2011;Langmead and Salzberg 2012). Samtools v1.2 was then used to create an mpileup file, which was filtered with VarScan v2.3.2 (minimum coverage ¼ 20, minor allele frequency ¼ 10%, and minimum average phred quality ¼ 20) to obtain SNPs. To identify nonsynonymous and synonymous segregating polymorphisms we identified the n-fold degenerate positions following Li et al. (1985) from which pN, pS, and (pN/pS) could be calculated per gene. Comparison of mean pN/pS and dN/dS between convergent and nonconvergent (background) genes was conducted using a permutation t-test (number of permutations ¼ 10,000) in R.

Cross-species Mapping
All of the above analyses used only the one-to-one orthologs. To examine a larger fraction of the transcriptome we produced species-pair references by using a reciprocal BLAST between the assemblies of sexual-asexual sister species (T. bartmani -Timema tahoe, Timema cristinae -Timema monikensis, Timema poppensis -Timema douglasi, Timema californicum -Timema Shepardi, and Timema podura -Timema genevievae) (BlastN, minimum e-val ¼ 0.00001, minimum query coverage ¼ 30%). Prior to this step potential contaminants were filtered from these by blasting transcripts to local versions of the nt (using BlastN, default options except task ¼ BlastN, max_target_seqs ¼ 10) and nr (using BlastX, default options except, max_target_seqs ¼ 10) databases (downloaded: 07/08/2016) using NCBI's BLAST client (v. 2.2.30þ). BLAST hits with an e-value > 0.0000001 were discarded. The remaining BLAST hits were used to assign a phylum to sequences if !50% of BLAST hits came from one phylum (in the event of a tie, the taxa with the highest e-value was used as a tiebreaker). Transcripts that were assigned to a nonarthropoda phylum were discarded (note that transcripts with no BLAST hits or that blasted to mixed phyla were retained). This filtering removed between 4% and 8% of Convergent Gene Expression Changes . doi:10.1093/molbev/msy217 MBE transcripts (see supplementary table 10, Supplementary Material online). Reads of each species were then mapped to each species-pair reference in the same way as for the 1-to-1 orthologs. Differential expression analyses and GO-term enrichment analyses were then repeated as described above.

Supplementary Material
Supplementary data are available at Molecular Biology and Evolution online.

Data
Raw reads have been deposited in the SRA. Accession codes are given in supplementary table 19, Supplementary Material online. Scripts for the analyses in this paper are available at: http://dx.doi.org/10.5281/zenodo.2025853.