Research articlesEvidence for cryptic sex in parthenogenetic stick insects of the genus Timema
Abstract
Obligately parthenogenetic species are expected to be short lived since the lack of sex and recombination should translate into a slower adaptation rate and increased accumulation of deleterious alleles. Some, however, are thought to have been reproducing without males for millions of years. It is not clear how these old parthenogens can escape the predicted long-term costs of parthenogenesis, but an obvious explanation is cryptic sex. In this study, we screen for signatures of cryptic sex in eight populations of four parthenogenetic species of Timema stick insects, some estimated to be older than 1 Myr. Low genotype diversity, homozygosity of individuals and high linkage disequilibrium (LD) unaffected by marker distances support exclusively parthenogenetic reproduction in six populations. However, in two populations (namely, of the species Timema douglasi and T. monikensis) we find strong evidence for cryptic sex, most likely mediated by rare males. These populations had comparatively high genotype diversities, lower LD, and a clear LD decay with genetic distance. Rare sex in species that are otherwise largely parthenogenetic could help explain the unusual success of parthenogenesis in the Timema genus and raises the question whether episodes of rare sex are in fact the simplest explanation for the persistence of many old parthenogens in nature.
1. Introduction
Although sex is the most prevalent mode of reproduction in animals, many species are able to reproduce via female-producing parthenogenesis. These species can reproduce without males, either exclusively (obligate parthenogens) or facultatively (facultative and cyclical parthenogens) [1]. Parthenogenesis is associated with considerable short-term advantages relative to sexual reproduction, the two main ones being reproductive assurance (since parthenogenetic organisms only require one individual for reproduction) [2] and demographic advantage (since parthenogenic reproduction produces only females) [3,4]. However, these short-term advantages are predicted to be coupled with long-term costs, including the accumulation of deleterious alleles [5,6], and a slower rate of adaptive evolution [7,8]. These long-term costs are believed to result in the eventual extinction of obligately parthenogenetic species [4,6].
Despite the expected long-term costs of parthenogenesis, several parthenogenetic species are thought to have been reproducing solely parthenogenetically for millions of years [9,10]. This raises the question of how these old parthenogens have managed to escape the predicted long-term consequences of parthenogenesis thus far. One possible explanation would be that besides reproducing parthenogenetically, they occasionally reproduce sexually (cryptic sex). Indeed, bouts of rare sex in mainly parthenogenetically reproducing populations would be sufficient to overcome many of the long-term costs of parthenogenesis and could prevent the mutational deterioration of otherwise clonal lineages [11–13]. Consistent with this idea, several species originally thought to be obligately parthenogenetic have been recently shown to be able to exchange genetic material between individuals (cryptic gene flow), including the brine shrimp Artemia parthenogenetica [14], and the bdelloid rotifers Adineta vaga [15] and Macrotrachella quadricornifera [16,17], questioning the existence of ancient asexuality (see also [18]).
In bdelloid rotifers, males have never been observed [19], however patterns of allele sharing between individuals are compatible with some form of rare or non-canonical sex [15–17]. On the other hand, in Artemia brine shrimps it has been demonstrated that cryptic gene flow can be mediated by rare males [14]. Similar to Artemia, rare males are documented in many obligately parthenogenetic species [20–23], but they typically do not appear to reproduce sexually with parthenogenetic females (reviewed in [24]). In some species, males produced by parthenogenetic females can mate with females of related sexual strains, but while such matings can result in new parthenogenetic lineages via contagious parthenogenesis [25,26], they do not mediate gene flow within parthenogenetic lineages. Nevertheless, formal tests of cryptic gene flow within parthenogenetic lineages remain scarce, and rare males may therefore mediate cryptic gene flow in a broader panel of ‘obligately' parthenogenetic species than generally assumed.
In this study, we assess whether cryptic gene flow occurs in parthenogenetic species of the stick insect genus Timema (Phasmatodea). This genus comprises five described parthenogenetic species [27–29] which represent at least seven independently derived parthenogenetic lineages [30]. These lineages vary in age, and the oldest parthenogenetic Timema is thought to have been reproducing via parthenogenesis for over 1 Myr [30], making them an ‘old' asexual species. Parthenogenetic Timema are very homozygous as they most likely reproduce via a form of automictic parthenogenesis where heterozygosity is completely lost each generation [31,32]. Rare males have been found in four of the five described parthenogenetic Timema species [28,29,33], and these rare ‘parthenogenetic' males present normal male reproductive organs, and can mate and produce offspring with sexual sister-species females, although with reduced fertility [33]. It is not clear how such males are produced; however, they are thought to originate from the accidental loss of an X chromosome (aneuploidy) during parthenogenesis [33]. Given the XX/X0 mechanism of sex determination in Timema [34], if an XX parthenogenetic egg loses one X it could potentially develop into a male, a mechanism documented in other parthenogenetic insect species [35,36].
To test for cryptic sex in Timema, we use RADseq genotyping to estimate linkage disequilibrium (LD) in eight populations of parthenogenetic species. A population specific approach is appropriate for studies with these species, as they live in patchy habitats and have very limited dispersal [37]. LD can serve as a measure of genetic exchange as with strict, long-term parthenogenesis alleles will be transmitted in a single block. Note that many forms of automictic parthenogenesis involve some level of recombination and are associated with a rapid loss of heterozygosity after the inception of parthenogenesis in a sexual ancestor [38]. Under strict, long-term parthenogens, such recombination does however not affect population LD as heterozygosity can no longer decrease, either because it has eroded completely or because remaining heterozygosity tracts are mechanistically or selectively constrained [39,40]. In such cases, LD will only be broken by mutation or allele conversion, and is thus expected to be high. Sexual events on the other hand, will decrease LD, even when they are very rare. While we find no evidence for cryptic sex in six populations, we show that in two of the eight studied populations of parthenogenetic Timema, LD within and between chromosomes is very low, and there is a consistent LD decay along the chromosomes. In one of these populations, we also find several individuals with unusually high heterozygosity, indicative of recent sexual events. The results presented here are the first evidence for cryptic sex, most likely mediated by rare males, in ‘obligately' parthenogenetic Timema and contribute to the growing evidence that rare or non-canonical sex might have to be considered when studying the long-term persistence of (mostly) parthenogenetically reproducing species.
2. Material and methods
(a) Samples and sequencing
We looked for evidence of gene flow in two populations from each of four parthenogenetic species of Timema (T. genevievae, T. shepardi, T. monikensis, T. douglasi), in a total of eight populations collected between 2010 and 2018. A population of the sexual species T. cristinae was also included to establish a sexual reference for LD and heterozygosity patterns based on RADseq genotypes. We used a population specific approach because Timema are wingless, and it has been estimated that they can travel a maximum of 128 m per generation [37], which is much less than the distance between the two studied populations in each parthenogenetic species (6, 31, 43 and 219 km apart on a straight line for T. monikensis, T. shepardi, T. douglasi and T. genevievae, respectively).
We genotyped 24 females for each of the eight parthenogenetic populations (192 females in total) and 24 individuals (12 females and 12 males) for the sexual population (see also electronic supplementary material, table S1). Additionally, because we found evidence for recent gene flow in one population of T. monikensis (see below) we included additional samples of T. monikensis from different collection years for both populations under study, which included four males from the FS population (electronic supplementary material, table S1).
DNA was extracted from the head or legs (electronic supplementary material, table S1) using the Qiagen DNeasy Blood and Tissue Kit, following the manufacturer's instructions. ddRAD libraries were prepared following the protocol from [41], with the enzymes EcoRI and MseI, and a 200–450 bp size selection after addition of Illumina TruSeq indexes. For logistical reasons, a few samples underwent a different size selection (300–500 bp). Ninety-six samples were pooled in each library thanks to barcoded adapters (see protocol in [41]). Libraries were multiplexed by pairs on two Illumina lanes using Illumina TruSeq indexes iA06 or iA12 and 150 bp single-end sequencing was performed using Illumina Hiseq 2500 at the Lausanne Genomic Technologies Facility. Reads are available under BioProject accession number PRJNA798556.
(b) Reads processing
All scripts used in this pipeline are available online (https://github.com/SusanaNFreitas/cryptic_gene_flow). We used the STACKS software (version 2.5) [42] to de-multiplex the reads with the ‘process radtags' module from the STACKS suite. Reads were trimmed from adaptors and filtered for minimum length (80 bp) with Cutadapt 2.3 [43] and aligned against the corresponding reference genomes [31] with BWA. Sam files were sorted and converted to bam using samtools 1.4 excluding multiple alignments and PHRED quality score of 30 or less. FreeBayes (v. 1.3.3) [44] was used to call SNPs, which were then filtered by read depth (with GATK, v. 4.2.0.0 [45], keeping only positions with DP > 8 and DP < 200), genotype quality (using a custom python script, keeping SNP genotypes with QUAL > 30) and missing data (with vcftools [46], only positions with 75% data kept). The vcfallelicprimitives function of vcflib [47] was used to break MNPs in the vcf file into SNPs, and only biallelic SNPs were kept (eliminating indels and multi-allelic positions with custom unix oneliners). Low coverage samples were excluded from further analyses (electronic supplementary material, table S1).
(c) Population genetics
From the RADseq reads mapped to each species' genome, we recovered 25 777, 17 997, 14 326 and 10 617 SNPs for T. genevievae, T. monikensis, T. douglasi and T. shepardi, respectively. These were then used to infer the number of clones in our populations by estimating pairwise genetic distances between individuals, visualized as neighbour-joining networks (figure 2). To explicitly look for cryptic sex, we then measured LD and LD decay in each of the eight parthenogenetic populations (figures 3–5; and electronic supplementary material, figures S2–S9) and estimated relative heterozygosity for each individual (figure 6).
The total number of SNPs per species was estimated with the R package ADEGENET 2.1.3 [48,49]. Genetic distances were estimated using the Euclidean distance with the dist() function. LD was estimated with plink 1.9 [50], the r2 value was estimated within and between linkage groups with options --r2 and --inter-chr, and LD decay within linkage group with options --r2, --ld-window 999999 and --ld-window-kb 8000. Only SNPs with minor allele frequency >0.2 were used to estimate r2. For LD decay, we only used SNPs with available chromosomal coordinates [31,51], which were based on chromosomal information from the species T. cristinae (assembly v1.3 [52]). The best-fit line estimated on the LD decay plots was generated with the geom_smooth() function using the gam method, or the LOESS method whenever fewer than 1000 data points were available (chromosomal interactions). The best-fit line was only estimated for chromosomes with more than 100 observations points. Because the X chromosome was misassembled in the v1.3 T. cristinae assembly, we used the positional information for SNPs on all chromosomes except the X [31,51]. r2 values per population were plotted with the package ggplot2 [53] in R. In the plots of the chromosomal LD, the black dot represents the mean r2 per chromosome. Ninety-five per cent confidence intervals were approximated using the standard error (s.e. = (standard deviation)/(square root of the sample size)), adding (for the upper limit) or subtracting (for the lower limit) to the mean twice the respective s.e. Relative heterozygosity was estimated for each individual by dividing the number of heterozygous positions by the total number of called polymorphic positions.
3. Results
(a) Genetic diversity patterns: sexual versus parthenogenetic
Sexual Timema populations, as exemplified by the T. cristinae population included here, present high relative heterozygosity and genotype diversity (electronic supplementary material, figure S1, see also [31]). As expected for a sexual population, where recombination and segregation of chromosomes break LD, r2 values both within and between chromosomes were very low (figure 1a), even between SNPs separated by small distances (figure 1b).
Figure 1. LD in the sexual species T. cristinae. (a) Violin plots of the LD values within (purple) and between (red) chromosomes (linkage groups). Black dot represents the mean LD (r2) and the black bar the approximate 95% confidence interval. (b) LD decay for each linkage group. LD decay was estimated for all chromosomes excluding the X, for which we did not have positional information (see Material and methods). Lines of best fit (produced using the gam method) are shown in blue.
In contrast to the sexual population with high genotype diversity, populations reproducing solely by parthenogenesis are prone to undergo recurrent sweeps. This should result in the fixation of one or few clones (genotypes), with minor genetic differences between individuals within clonal lineages. In Timema, strict parthenogenesis would further be associated with very low heterozygosity for all individuals in a population, since their parthenogenesis mechanism leads to the complete or largely complete loss of heterozygosity every generation [31,32]. On the other hand, if parthenogenetic Timema populations are undergoing rare sex we would expect to see an increase in the genotypic diversity, and perhaps some individuals with elevated heterozygosity, which would indicate that they were produced via sex.
We found two distinct patterns in the genetic diversity estimates (pairwise distances, LD and heterozygosity) among the populations of the parthenogenetic species. In six of the eight populations, we found no evidence for cryptic sex, as the observed patterns clearly matched the expectations for a purely parthenogenetic population: genotype diversities were very low (one or two clones per population; figure 2), and LD was high (figure 3a; electronic supplementary material, figures S2–3, S4-A and S5-A), and did not decay over increasing genetic distances (electronic supplementary material, figures S6–S7, S8-A and S9-A). We found clear evidence for gene flow in the remaining two populations, with high diversity of genotypes (figure 2), low LD (figure 4) and evidence of LD decay (figure 5).
Figure 2. Neighbour-joining networks of pairwise distances between individuals within parthenogenetic Timema species: (a) T. genevievae, (b) T. monikensis, (c) T. douglasi and (d) T. shepardi. Individuals of the two studied populations in each species are distinguished by brown and blue labels. Figure 3. LD in the parthenogenetic species T. genevievae (population HW20). (a) Violin plots of the LD values within (purple) and between (red) chromosomes (linkage groups). Black dot represents the mean LD (r2) and the black line the approximate 95% confidence interval. (b) LD decay for each linkage group. LD decay was estimated for all chromosomes excluding the X, for which we did not have positional information (see Material and methods). Lines of best fit are shown in blue made using the gam method (chromosomes 1–8 and 10) and LOESS (chromosomes 9, 11). The line of best fit was only estimated when at least 100 datapoints were available. LD values for T. genevievae—population HW20 are presented here as an illustrative example but see electronic supplementary material, figures S2–S9 for all parthenogenetic populations. Figure 4. LD values in parthenogenetic populations showing evidence for gene flow. LD values within (purple) and between (red) chromosomes (linkage groups) represented in violin plots, with the black dot representing the mean LD (r2) and the black bar the approximate 95% confidence interval. (a) T. douglasi, population Manchester. (b) T. monikensis, population FS. Figure 5. LD decay per chromosome in parthenogenetic populations showing evidence for gene flow. LD decay was estimated for all chromosomes excluding the X, for which we did not have positional information (see Material and methods). Lines of best fit are shown in blue. (a) T. douglasi, population Manchester. Lines of best fit were estimated with the gam method for all chromosomes except for chromosome 12 for which the LOESS method was used. (b) T. monikensis, population FS. Lines of best fit were estimated with the gam method for all chromosomes except for chromosomes 1, 5 and 12 for which the LOESS method was used.
(b) Populations with no evidence for sex
The six populations that presented no evidence for cryptic sex (hereafter obligately parthenogenetic populations), had either one clonal lineage, namely T. genevievae—Antonio and T. shepardi—Elk Mountain, or two clonal lineages, T. genevievae—HW20, T. shepardi—Spring Road, T. douglasi—Black Bart Rd and T. monikensis—Sycamore (figure 2). Even though we estimated LD values for all cases (electronic supplementary material, figures S2–S9), estimates for populations with only a single clonal lineage are not informative given the extremely low polymorphism they present. Nevertheless, the presence of a single clone in those populations is in itself consistent with the expectation for the (at least short-term) absence of sex. For the four populations with two clonal lineages, LD values within and between chromosomes were consistently very high (figure 3a; electronic supplementary material, figures S2–S5), indicative of a tight linkage between polymorphic positions. There was also no evidence of LD decay (figure 3b; electronic supplementary material, figures S6–S9), consistent with linkage between polymorphic positions and the absence of recombination. Finally, individuals from all six obligately parthenogenetic populations (including the two for which LD estimates are not informative due to low polymorphism) presented consistently low heterozygosity (figure 6).
Figure 6. Relative heterozygosity per individual in each of the eight studied populations. (a) T. genevievae, (b) T. monikensis, (a) T. douglasi, and (d) T. shepardi. On the x-axis are the populations for the four parthenogenetic species studied, plus the males in the T. monikensis—FS population. On the y-axis is the relative heterozygosity estimated as the median between all chromosomes.
(c) Populations with evidence for sex
The patterns in the remaining two populations (T. douglasi—Manchester, and T. monikensis—FS) are, however, much different and best explained by (rare) sexual reproduction between parthenogenetic individuals. Both populations present a high diversity of genotypes, with eight different clonal lineages in T. douglasi—Manchester, but no clonal lineages and only independent genotypes (21) in T. monikensis—FS (figure 2). Regarding LD within and between chromosomes, both populations presented much lower values of r2 when compared with the strictly parthenogenetic populations (figure 3; electronic supplementary material, figures S2–S5), yet values were still elevated relative to the sexual T. cristinae population (figure 1). LD decay is also apparent in these two populations, with estimates consistently decreasing with increasing distance between markers (figure 5). Furthermore, for T. monikensis—FS, r2 values within linkage groups were higher than r2 values between different chromosomes (figure 4b), suggesting that segregation contributes more strongly to the high genotype diversity in this population than recombination.
(d) Sexually produced individuals
Given the evidence for cryptic gene flow in the T. douglasi—Manchester and T. monikensis—FS populations, we investigated whether we could identify sexually produced individuals. Parthenogenetically produced Timema individuals are largely or completely homozygous [31,32]. By contrast, sex between different clones would manifest as elevated heterozygosity in the resulting offspring. In other words, the presence of individuals with elevated heterozygosity would indicate events of sex in the previous generation.
As expected, the relative heterozygosity values for the individuals in the six fully parthenogenetic populations were all very low, varying between 0 and 0.05. By contrast, four out of the 21 T. monikensis females from the population FS showed elevated relative heterozygosity (figure 6), corroborating our previous conclusions of cryptic sex in this population. The heterozygosity levels of these four individuals are as expected if these individuals were produced from crosses between individuals with low heterozygosity from the same population (electronic supplementary material, figure S10). On the other hand, none of the 20 T. douglasi females from the Manchester population had elevated heterozygosity, indicating that none of these females were produced sexually.
Finally, to investigate whether gene flow in T. monikensis—FS is mediated by rare males, and to assess whether these males are produced sexually or stem from accidental X chromosome losses, we also estimated relative heterozygosity for the four available males from that population (collected in a different year). If produced sexually, males would present elevated heterozygosity, similar to the four heterozygous females, but if resultant from X loss, males would present low heterozygosity levels. All four males were heterozygotic outliers, similar to the four highly heterozygous females (figure 6) indicating they were likely produced via sex.
(e) Additional analyses
To corroborate our conclusions of cryptic gene flow in the T. douglasi—Manchester and the T. monikensis—FS populations we performed three additional complementary analyses. First, we found no evidence that highly heterozygous individuals were triploid, as would be expected from the rare fertilization of parthenogenetically produced diploid eggs (electronic supplementary material, figure S11). Second, we tested if our LD findings in the T. monikensis—FS population were driven by the heterozygous outliers. To do this, we repeated our LD analyses with heterozygous outlier individuals removed and found that LD and LD decay remained very similar (electronic supplementary material, figure S12). Finally, to test whether LD and heterozygosity in T. monikensis populations (both Sycamore and FS) were consistent between different years, we analysed additional individuals from previous collection years (namely, 2013 and 2015). We found similar results to our main analyses: no evidence for sex in Sycamore, and strong evidence for cryptic sex in the FS population from low LD (electronic supplementary material, figure S13), evident LD decay (electronic supplementary material, figure S14) and heterozygotic outlier individuals (electronic supplementary material, figure S15).
4. Discussion
While obligately parthenogenetic lineages are expected to be short lived, many are thought to have been able to persist in the long term, sometimes even for millions of years ([9,10], but see [18]). It is still unknown how parthenogenetic lineages manage to persist for so long and evade the costs of lacking sex. However, rare events of sexual reproduction have been suggested as an explanation for the persistence of long lived parthenogens, and evidence for rare gene flow in ‘obligately' parthenogenetic species is accumulating [14–17].
Rare events of sex in obligately parthenogenetic species can be difficult to demonstrate, since unless we find direct evidence (such as direct observations of sexual offspring in broods of parthenogenetic females), we need to rely on indirect evidence. Indirect evidence may include the presence of functional males, unexpectedly high genetic diversity in parthenogenetic populations or signatures of recombination (e.g. recombining haplotypes, low LD). Here, we searched for signatures of recombination and sex in parthenogenetic species of the genus Timema by estimating genotype diversity within populations, LD, LD decay with increasing marker distance and heterozygosity. Some forms of parthenogenesis can involve recombination and segregation [38,54], but in homozygous parthenogens such as Timema [31,32], these two mechanisms have little to no effect on genotype diversity. In such homozygous parthenogens, signatures of recombination can therefore provide very strong evidence for sexual reproduction.
While we found no evidence for cryptic sex in six populations, we detected signatures of recombination and cryptic sex in two populations of parthenogenetic Timema, the population Manchester of T. douglasi and the population FS of T. monikensis. Cryptic sex, at least in T. monikensis, is almost certainly mediated by rare males, known to occur in Timema parthenogens. Three lines of evidence indicate that cryptic sex within parthenogenetic populations is a much more likely explanation for our findings than ‘contagious parthenogenesis', where males produced by parthenogenetic lineages mate with sexual females and thereby generate new parthenogenetic strains via introgression. First, sexual and parthenogenetic sister species do not overlap geographically in Timema, with the closest known populations being separated by hundreds of kilometres [55]. Second, sexual and parthenogenetic populations are genetically diverged [31], such that introgression would generate genetic distances beyond the intra-population distances we describe (figure 2). Finally, the heterozygosity of the putatively sexually produced individuals fit the expected heterozygosity modelled for sex within populations (electronic supplementary material, figure S10).
The uncovered patterns in the populations with signatures for cryptic sex, Manchester and FS populations, reflected different scenarios. Regarding the Manchester population (T. douglasi), LD within and between chromosomes was generally low (and similar), and a clear pattern of LD decay was evident across all linkage groups. Genotype diversity was relatively high (i.e. we found eight different ‘clones' among the 20 genotyped individuals), and clones were often represented by multiple individuals (one to five per clone). Finally, the heterozygosity of all females was as low as in the strictly parthenogenetic populations, indicating they were produced via parthenogenesis. These findings are best explained by rare sexual events that occurred in the past, leading to a high diversity of genotypes, but also allowing for the spread of these genotypes into clonal lineages.
The second population with evidence for gene flow, the FS population of T. monikensis, likely features more, and more recent sexual events than the Manchester population of T. douglasi. More sexual events could result via a more frequent production of males, higher propensity of females to reproduce sexually when mated, and/or for longer periods of time, for instance. The 21 genotyped females represented 21 distinct genotypes, as would be expected for a sexual but not for a largely parthenogenetic population. Consistent with the high genotype diversity, LD decay was evident for all chromosomes and LD values were generally very low. Finally, four out of the 21 females (approx. 19%) were characterized by elevated heterozygosity, and were most likely produced from a sexual cross between different genotypes (electronic supplementary material, figure S6). Timema parthenogenetic females were previously screened for their ability to fertilize eggs by using crosses with males from related sexual species (T. cristinae males for T. monikensis females of the FS population), but none of the 322 genotyped hatchlings were sexually produced [33]. Even though the absence of sexually produced offspring in those crosses could partly stem from species divergence and associated incompatibilities, the reciprocal crosses (T. cristinae females mated to T. monikensis males) result in large numbers of hybrid offspring [33]. Given our new results indicating T. monikensis are capable of reproducing sexually, this suggests that sexually produced offspring were too rare to be detected in the analysed sample of hatchlings, yet we here report that as many as 19% of the adult females of the FS population were sexually produced. One possible explanation for such a frequency difference could be that sexually produced females have a survival advantage over parthenogenetically produced ones. If that is the case, then the maintenance of mostly parthenogenetic reproduction in the FS population in spite of this apparent selection for sex will be an interesting focus for future research.
Even though we do find evidence for genetic exchange in two parthenogenetic populations, most of the populations analysed here are de facto parthenogenetic, and do not present any evidence of gene flow at present or in the recent past. Two scenarios could potentially explain this patchy occurrence of sex. In the first scenario, a parthenogenetic all-female population can produce males spontaneously (via e.g. X aneuploidies), and parthenogenetic females are then able to mate and produce at least some sexual offspring. Such sexual offspring would comprise males, which would then allow for the propagation of males in the population. In the second scenario, males and sexual reproduction were originally present in all populations since the split with the closest sexual ancestor and were gradually lost or even extinct in the majority of recent parthenogenetic populations. In this scenario, populations that have males at present have somehow maintained rare sex and sexually produced offspring (and thus males) perhaps through selection. In both scenarios, de facto parthenogenetic populations have been reproducing asexually for the last generations, but we cannot rule out these populations have not undergone episodes of (rare) sexual reproduction in the distant past.
Why sex occurs in some populations but not others is unclear and may also be a dynamic process. For example, if sex can emerge spontaneously in a population it may persist for multiple generations before going extinct due to the short-term advantages of parthenogenesis. In species that inhabit fire-prone habitats, such as Timema, an ecologically important short-term advantage of parthenogenesis is the ability to rapidly recolonize areas following fire with a single (female) individual. As such, the recurrent selection for parthenogenesis following fire, combined with the inherent stochasticity of extinction–recolonization dynamics of fire-prone habitats, may drive a patchy distribution of sex in Timema parthenogenetic populations.
As mentioned above, males are known in four out of the five described parthenogenetic species of Timema, including T. monikensis [27–29,33]. Our findings of cryptic sex in one population of T. monikensis now raise the question of whether males in other Timema species could also mediate cryptic sex in some populations, albeit much more rarely than in the FS population of T. monikensis. In combination with other recent reports of gene flow in parthenogenetic animal species [14–17], our findings suggest that episodes of rare sex in largely parthenogenetic species may occur more frequently than generally appreciated. Thus, it would be useful to screen additional parthenogenetic species for signatures of sex, especially species with only suggestive evidence that sexual reproduction may still be ongoing, such as parthenogenetic species with functional males, or unexpectedly high levels of genetic diversity. Males have been found in several parthenogenetic species (e.g. [20–23]), and many examples exist for high genetic diversity in parthenogenetically reproducing populations. Both these findings have been used to suggest cryptic gene flow in parthenogenetic species [56,57], but neither are reliable evidence for sex. The presence of males itself may not necessarily mean that those males can mate with parthenogenetic females, since traits required for sexual reproduction are typically vestigialized in parthenogenetic females [24]. Furthermore, males can be dysfunctional, especially in older parthenogenetic lineages where mutations may have accumulated in pathways for exclusively male functions [24]. Males produced by parthenogenetic females could also mate with sexual females, such as in the cases of contagious parthenogenesis [22,58]. Accordingly, high genetic diversity in parthenogenetic species could be explained by multiple transitions to parthenogenesis [59–61] and not indicate cryptic sex. However, future studies of such species could reveal cryptic sex in at least some of them. Finally, despite the evidence for cryptic sex in some parthenogenetic species, evidence for long-term exclusive parthenogenesis in eukaryotes has also been reported for protozoans [62] and oribatid mites [63]. Since both these groups are characterised by extremely large population sizes, these findings could be an indication of support for theory showing diminishing benefits of sex with large population size [6,64].
In conclusion, independently of how many additional parthenogenetic species will be discovered to reproduce sexually (on rare occasions), the fact that at least some presumed obligate parthenogens can reproduce sexually questions our understanding about the plasticity of sexual and parthenogenetic reproduction modes in animals. When classifying a species as sexual, we expect that all offspring are produced by sex, yet many species are capable of spontaneous parthenogenesis [54]. Accordingly, when classifying a species as obligately parthenogenetic we rarely consider occasional and/or non-canonical sex. While most reproductive systems in animals fit into a bimodal pattern of ‘sexual' or ‘parthenogenetic', we also concur with [14], that reproductive systems may contain more variation than previously assumed, and believe this could in many ways explain how old parthenogens escaped the long-term costs of the ‘lack of sex'.
Ethics
This work did not require ethical approval from a human subject or animal welfare committee.
Data accessibility
Raw sequence reads have been deposited in NCBI's sequence read archive under the bioproject: PRJNA798556. Scripts for the analyses in this paper are available at: https://doi.org/10.5281/zenodo.8014163 [65].
Supplementary material is available online [66].
Declaration of AI use
We have not used AI-assisted technologies in creating this article.
Authors' contributions
S.F.: conceptualization, data curation, formal analysis, investigation, methodology, project administration, visualization, writing—original draft, writing—review and editing; D.J.P.: data curation, formal analysis, investigation, methodology, writing—review and editing; M.L.: investigation, methodology; Z.D.: investigation, methodology; T.S.: conceptualization, funding acquisition, project administration, supervision, writing—review and editing.
All authors gave final approval for publication and agreed to be held accountable for the work performed therein.
Conflict of interest declaration
We declare we have no competing interests.
Funding
This study was supported by the European Research Council (Consolidator Grant No Sex No Conflict) and Swiss FNS grant no. 31003A_182495 to T.S.
Acknowledgements
The authors thank Roger K. Butlin for his constructive comments to the methods used here, and Chloe Larose and Bart Zjilstra for support during fieldwork. Most analyses were run using the computing infrastructure at the University of Lausanne maintained by the former Vital-IT platform of the SIB Swiss Institute of Bioinformatics (SIB) and the DCSR. We would like to acknowledge funding from the European Research Council (Consolidator Grant No Sex No Conflict) and Swiss FNS grant 31003A_182495 to T.S.
References
- 1.
Bell G . 1982 The masterpiece of nature: the evolution and genetics of sexuality. Berkeley, CA: University of California Press. Google Scholar - 2.
Gerritsen J . 1980 Sex and parthenogenesis in sparse populations. Am. Nat. 115, 718-742. (doi:10.1086/283594) Crossref, Web of Science, Google Scholar - 3.
- 4.
- 5.
Muller HJ . 1964 The relation of recombination to mutational advance. Mutat. Res. 1, 2-9. (doi:10.1016/0027-5107(64)90047-8) Crossref, Web of Science, Google Scholar - 6.
Gabriel W, Lynch M, Bürger R . 1993 Muller's ratchet and mutational meltdowns. Evolution 47, 1744-1757. (doi:10.2307/2410218) Crossref, PubMed, Web of Science, Google Scholar - 7.
Hill WG, Robertson A . 1966 The effect of linkage on limits to artificial selection. Genet. Res. 8, 269-294. (doi:10.1017/S0016672300010156) Crossref, PubMed, Web of Science, Google Scholar - 8.
Barton NH . 1995 Linkage and the limits to natural selection. Genetics 140, 821-841. (doi:10.1093/genetics/140.2.821) Crossref, PubMed, Web of Science, Google Scholar - 9.
Judson OP, Normark BB . 1996 Ancient asexual scandals. Trends Ecol. Evol. 11, 41-46. (doi:10.1016/0169-5347(96)81040-8) Crossref, PubMed, Web of Science, Google Scholar - 10.
Schurko AM, Neiman M, Logsdon JM . 2009 Signs of sex: what we know and how we know it. Trends Ecol. Evol. 24, 208-217. (doi:10.1016/j.tree.2008.11.010) Crossref, PubMed, Web of Science, Google Scholar - 11.
Green RF, Noakes DLG . 1995 Is a little bit of sex as good as a lot? J. Theor. Biol. 174, 87-96. (doi:10.1006/jtbi.1995.0081) Crossref, Web of Science, Google Scholar - 12.
Bengtsson BO . 2009 Asex and evolution: a very large-scale overview. In Lost sex (edsSchön I, Martens K, van Dijk P ), pp. 1-19. Dordrecht, The Netherlands: Springer. Crossref, Google Scholar - 13.
D'Souza TG, Michiels NK . 2010 The costs and benefits of occasional sex: theoretical predictions and a case study. J. Hered. 101, S34-S41. (doi:10.1093/jhered/esq005) Crossref, PubMed, Web of Science, Google Scholar - 14.
Boyer L, Jabbour-Zahab R, Mosna M, Haag CR, Lenormand T . 2021 Not so clonal asexuals: unraveling the secret sex life of Artemia parthenogenetica. Evol. Lett. 5, 164-174. (doi:10.1002/evl3.216) Crossref, PubMed, Web of Science, Google Scholar - 15.
Vakhrusheva OA 2020 Genomic signatures of recombination in a natural population of the bdelloid rotifer Adineta vaga. Nat. Commun. 11, 6421. (doi:10.1038/s41467-020-19614-y) Crossref, PubMed, Web of Science, Google Scholar - 16.
Signorovitch A, Hur J, Gladyshev E, Meselson M . 2015 Allele sharing and evidence for sexuality in a mitochondrial clade of bdelloid rotifers. Genetics 200, 581-590. (doi:10.1534/genetics.115.176719) Crossref, PubMed, Web of Science, Google Scholar - 17.
Laine VN, Sackton TB, Meselson M . 2022 Genomic signature of sexual reproduction in the bdelloid rotifer Macrotrachella quadricornifera. Genetics 220, iyab221. (doi:10.1093/genetics/iyab221) Crossref, PubMed, Web of Science, Google Scholar - 18.
Schwander T . 2016 Evolution: the end of an ancient asexual scandal. Curr. Biol. 26, R233-R235. (doi:10.1016/j.cub.2016.01.034) Crossref, PubMed, Web of Science, Google Scholar - 19.
Birky CW . 2010 Positively negative evidence for asexuality. J. Hered. 101, S42-S45. (doi:10.1093/jhered/esq014) Crossref, PubMed, Web of Science, Google Scholar - 20.
Simon JC, Blackman RL, Gallic JFL . 1991 Local variability in the life cycle of the bird cherry-oat aphid, Rhopalosiphum padi (Homoptera: Aphididae) in western France. Bull. Entomol. Res. 81, 315-322. (doi:10.1017/S0007485300033599) Crossref, Web of Science, Google Scholar - 21.
Butlin R, Schön I, Martens K . 1998 Asexual reproduction in nonmarine ostracods. Heredity 81, 473-480. (https://www.nature.com/articles/6884540#:~:text=Asexual%20reproduction%20has%20evolved%20repeatedly,Clonal%20diversity%20is%20highly%20variable) Crossref, Web of Science, Google Scholar - 22.
Delmotte F, Leterme N, Bonhomme J, Rispe C, Simon J-C . 2001 Multiple routes to asexuality in an aphid species. Proc. R. Soc. Lond. B 268, 2291-2299. (doi:10.1098/rspb.2001.1778) Link, Web of Science, Google Scholar - 23.
Tvedte ES, Ward AC, Trendle B, Forbes AA, Logsdon JM . 2020 Genome evolution in a putatively asexual wasp. bioRxiv, 2020.12.23.424202. (doi:10.1101/2020.12.23.424202) Google Scholar - 24.
van der Kooi CJ, Schwander T . 2014 On the fate of sexual traits under asexuality. Biol. Rev. 89, 805-819. (doi:10.1111/brv.12078) Crossref, PubMed, Web of Science, Google Scholar - 25.
Innes DJ, Hebert PDN . 1988 The origin and genetic basis of obligate parthenogenesis in Daphnia pulex. Evolution 42, 1024-1035. (doi:10.2307/2408918) Crossref, PubMed, Web of Science, Google Scholar - 26.
Maccari M, Amat F, Hontoria F, Gómez A . 2014 Laboratory generation of new parthenogenetic lineages supports contagious parthenogenesis in Artemia. PeerJ 2, e439. (doi:10.7717/peerj.439) Crossref, PubMed, Web of Science, Google Scholar - 27.
Sandoval CP, Vickery VR . 1996 Timema douglasi (Phasmatoptera: Timematodea), a new parthenogenetic species from SouthWestern Oregon and Northern California, with notes on other species. Can. Entomol. 128, 79-84. (doi:10.4039/Ent12879-1) Crossref, Web of Science, Google Scholar - 28.
Vickery VR, Sandoval CP . 2001 Descriptions of three new species of Timema (Phasmatoptera: Timematodea: Timematidae) and notes on three other species. J. Orthoptera Res. 10, 53-61. (doi:10.1665/1082-6467(2001)010[0053:DOTNSO]2.0.CO;2) Crossref, Google Scholar - 29.
Vickery VR, Sandoval CP . 1999 Two new species of Timema (Phasmatoptera: Timematodea: Timematidae), one parthenogenetic, in California. J. Orthoptera Res. 8, 45-47. Crossref, Google Scholar - 30.
Schwander T, Henry L, Crespi BJ . 2011 Molecular evidence for ancient asexuality in Timema stick insects. Curr. Biol. 21, 1129-1134. (doi:10.1016/j.cub.2011.05.026) Crossref, PubMed, Web of Science, Google Scholar - 31.
Jaron KS 2022 Convergent consequences of parthenogenesis on stick insect genomes. Sci. Adv. 8, eabg3842. (doi:10.1126/sciadv.abg3842) Crossref, PubMed, Web of Science, Google Scholar - 32.
Larose C, Lavanchy G, Freitas S, Parker DJ, Schwander T . 2023 Facultative parthenogenesis: a transient state in transitions between sex and obligate asexuality in stick insects? Peer Community J. 3, e60. (doi:10.24072/pcjournal.283) Crossref, Google Scholar - 33.
Schwander T, Crespi BJ, Gries R, Gries G . 2013 Neutral and selection-driven decay of sexual traits in asexual stick insects. Proc. R. Soc. B 280, 20130823. Link, Web of Science, Google Scholar - 34.
Schwander T, Crespi BJ . 2009 Multiple direct transitions from sexual reproduction to apomictic parthenogenesis in Timema stick insects. Evolution 63, 84-103. (doi:10.1111/j.1558-5646.2008.00524.x) Crossref, PubMed, Web of Science, Google Scholar - 35.
Pijnacker LP, Ferwerda MA . 1980 Sex chromosomes and origin of males and sex mosaics of the parthenogenetic stick insect Carausius morosus Br. Chromosoma 79, 105-114. (doi:10.1007/BF00328476) Crossref, Web of Science, Google Scholar - 36.
Blackman RL, Hales DF . 1986 Behaviour of the X chromosomes during growth and maturation of parthenogenetic eggs of Amphorophora tuberculata (Homoptera, Aphididae), in relation to sex determination. Chromosoma 94, 59-64. (doi:10.1007/BF00293530) Crossref, Web of Science, Google Scholar - 37.
Sandoval C . 2000 Persistence of a walking-stick population (Phasmatoptera: Timematodea) after a wildfire. Southwest. Nat. 45, 123. (doi:10.2307/3672452) Crossref, Web of Science, Google Scholar - 38.
Engelstädter J . 2017 Asexual but not clonal: evolutionary processes in automictic populations. Genetics 206, 993-1009. (doi:10.1534/genetics.116.196873) Crossref, PubMed, Web of Science, Google Scholar - 39.
Smith NMA 2019 Strikingly high levels of heterozygosity despite 20 years of inbreeding in a clonal honey bee. J. Evol. Biol. 32, 144-152. (doi:10.1111/jeb.13397) Crossref, PubMed, Web of Science, Google Scholar - 40.
Jaron KS, Bast J, Nowell RW, Ranallo-Benavidez TR, Robinson-Rechavi M, Schwander T . 2021 Genomic features of parthenogenetic animals. J. Hered. 112, 19-33. (doi:10.1093/jhered/esaa031) Crossref, PubMed, Web of Science, Google Scholar - 41.
Brelsford A, Dufresnes C, Perrin N . 2016 High-density sex-specific linkage maps of a European tree frog (Hyla arborea) identify the sex chromosome without information on offspring sex. Heredity 116, 177-181. (doi:10.1038/hdy.2015.83) Crossref, PubMed, Web of Science, Google Scholar - 42.
Catchen J, Hohenlohe PA, Bassham S, Amores A, Cresko WA . 2013 Stacks: an analysis tool set for population genomics. Mol. Ecol. 22, 3124-3140. (doi:10.1111/mec.12354) Crossref, PubMed, Web of Science, Google Scholar - 43.
Martin M . 2011 Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal 17, 10-12. (doi:10.14806/ej.17.1.200) Crossref, Google Scholar - 44.
Garrison E, Marth G . 2012 Haplotype-based variant detection from short-read sequencing. arXiv, 1207.3907. (doi:10.48550/arXiv.1207.3907) Google Scholar - 45.
Rimmer A, Phan H, Mathieson I, Iqbal Z, Twigg SR , WGS500 Consortium.Wilkie AO, McVean G, Lunter G . 2014. Integrating mapping-, assembly- and haplotype-based approaches for calling variants in clinical sequencing applications. Nat. Genet. 46, 912-918. (doi:10.1038/ng.3036) Crossref, PubMed, Web of Science, Google Scholar - 46.
Danecek P 2011 The variant call format and VCFtools. Bioinformatics 27, 2156-2158. (doi:10.1093/bioinformatics/btr330) Crossref, PubMed, Web of Science, Google Scholar - 47.
Garrison E, Kronenberg ZN, Dawson ET, Pedersen BS, Prins P . 2022 A spectrum of free software tools for processing the VCF variant call format: vcflib, bio-vcf, cyvcf2, hts-nim and slivar. PLoS Comput. Biol. 18, e1009123. (doi:10.1371/journal.pcbi.1009123) Crossref, PubMed, Web of Science, Google Scholar - 48.
Jombart T . 2008 adegenet: An R package for the multivariate analysis of genetic markers. Bioinformatics 24, 1403-1405. (doi:10.1093/bioinformatics/btn129) Crossref, PubMed, Web of Science, Google Scholar - 49.
Jombart T, Ahmed I . 2011 adegenet 1.3-1: New tools for the analysis of genome-wide SNP data. Bioinformatics 27, 3070-3071. (doi:10.1093/bioinformatics/btr521) Crossref, PubMed, Web of Science, Google Scholar - 50.
Purcell S 2007 PLINK: a tool set for whole-genome association and population-based linkage analyses. Am. J. Hum. Genet. 81, 559-575. (doi:10.1086/519795) Crossref, PubMed, Web of Science, Google Scholar - 51.
Parker DJ, Jaron KS, Dumas Z, Robinson-Rechavi M, Schwander T . 2022 X chromosomes show relaxed selection and complete somatic dosage compensation across Timema stick insect species. J. Evol. Biol. 35, 1734-1750. (doi:10.1111/jeb.14075) Crossref, PubMed, Google Scholar - 52.
Nosil P, Villoutreix R, Carvalho CF, Farkas TE, Soria-Carrasco V, Feder JL, Crespi BJ, Gompert Z . 2018. Natural selection and the predictability of evolution in Timema stick insects. Science 359, 765-770. (doi:10.1126/science.aap9125) Crossref, PubMed, Web of Science, Google Scholar - 53.
Wickham H . 2016 Ggplot2: elegant graphics for data analysis. New York, NY: SpringerVerlag. Crossref, Google Scholar - 54.
White MJD . 1973 Animal cytology and evolution, 3rd ed. Cambridge, UK: Cambridge University Press. Google Scholar - 55.
Law JH, Crespi BJ . 2002 The evolution of geographic parthenogenesis in Timema walking-sticks. Mol. Ecol. 11, 1471-1489. (doi:10.1046/j.1365-294X.2002.01547.x) Crossref, PubMed, Web of Science, Google Scholar - 56.
Belshaw R, Quicke DLJ, Völkl W, Godfray HCJ . 1999 Molecular markers indicate rare sex in a predominantly asexual parasitoid wasp. Evolution 53, 1189-1199. (doi:10.2307/2640822) Crossref, PubMed, Web of Science, Google Scholar - 57.
Fontcuberta A, Dumas Z, Schwander T . 2016 Extreme genetic diversity in asexual grass thrips populations. J. Evol. Biol. 29, 887-899. (doi:10.1111/jeb.12843) Crossref, PubMed, Web of Science, Google Scholar - 58.
Simon J-C, Leterme N, Latorre A . 1999 Molecular markers linked to breeding system differences in segregating and natural populations of the cereal aphid Rhopalosiphum padi L. Mol. Ecol. 8, 965-973. (doi:10.1046/j.1365-294x.1999.00648.x) Crossref, PubMed, Web of Science, Google Scholar - 59.
Fontaneto D, Boschetti C, Ricci C . 2008 Cryptic diversification in ancient asexuals: evidence from the bdelloid rotifer Philodina flaviceps. J. Evol. Biol. 21, 580-587. (doi:10.1111/j.1420-9101.2007.01472.x) Crossref, PubMed, Web of Science, Google Scholar - 60.
Janko K, Kotusz J, Gelas KD, Šlechtová V, Opoldusová Z, Drozd P, Choleva L, Popiołek M, Baláž M . 2012 Dynamic formation of asexual diploid and polyploid lineages: multilocus analysis of Cobitis reveals the mechanisms maintaining the diversity of clones. PLoS One 7, e45384. (doi:10.1371/journal.pone.0045384) Crossref, PubMed, Web of Science, Google Scholar - 61.
Capel KCC, Toonen RJ, Rachid CTCC, Creed JC, Kitahara MV, Forsman Z, Zilberberg C . 2017 Clone wars: asexual reproduction dominates in the invasive range of Tubastraea spp. (Anthozoa: Scleractinia) in the South-Atlantic Ocean. PeerJ 5, e3873. (doi:10.7717/peerj.3873) Crossref, PubMed, Web of Science, Google Scholar - 62.
Weir W 2016 Population genomics reveals the origin and asexual evolution of human infective trypanosomes. eLife 5, e11473. (doi:10.7554/eLife.11473) Crossref, PubMed, Web of Science, Google Scholar - 63.
Brandt A 2021 Haplotype divergence supports long-term asexuality in the oribatid mite Oppiella nova. Proc. Natl Acad. Sci. USA 118, e2101485118. (doi:10.1073/pnas.2101485118) Crossref, PubMed, Web of Science, Google Scholar - 64.
Gordo I, Charlesworth B . 2000 The degeneration of asexual haploid populations and the speed of Muller's ratchet. Genetics 154, 1379-1387. (doi:10.1093/genetics/154.3.1379) Crossref, PubMed, Web of Science, Google Scholar - 65.
Freitas S, Parker DJ . 2023SusanaFreitas/cryptic_geneflow_scripts: v.1.0.0. Zenodo . (doi:10.5281/zenodo.8014163) Google Scholar - 66.
Freitas S, Parker DJ, Labédan M, Dumas Z, Schwander T . 2023Evidence for cryptic sex in parthenogenetic stick insects of the genus Timema . Figshare. (doi:10.6084/m9.figshare.c.6824105) Google Scholar
