No recombination suppression in asexually produced males of Daphnia pulex

Abstract

Obligate parthenogenesis (OP) is often thought to evolve by disruption of reductional meiosis and suppression of crossover recombination. In the crustacean Daphnia pulex, OP lineages, which have evolved from cyclical parthenogenetic (CP) ancestors, occasionally produce males that are capable of reductional meiosis. Here, by constructing high-density linkage maps, we find that these males show only slightly and nonsignificantly reduced recombination rates compared to CP males and females. Both meiosis disruption and recombination suppression are therefore sex-limited (or partly so), which speaks against the evolution of OP by disruption of a gene that is essential for meiosis or recombination in both sexes. The findings may be explained by female-limited action of genes that suppress recombination, but previously identified candidate genes are known to be expressed in both sexes. Alternatively, and equally consistent with the data, OP might have evolved through a reuse of the parthenogenesis pathways already present in CP and through their extension to all events of oogenesis. The causal mutations for the CP to OP transition may therefore include mutations in genes involved in oogenesis regulation and may not necessarily be restricted to genes of the “meiosis toolkit.” More generally, our study emphasizes that there are many ways to achieve asexuality, and elucidating the possible mechanisms is key to ultimately identify the genes and traits involved.

Introduction

The mechanisms of evolutionary transitions to obligate parthenogenesis (OP) remain poorly understood, though it has become clear that these transitions more often occur through modifications of meiosis rather than through replacing meiosis by mitosis (Lenormand et al., 2016; Lynch & Conery, 2000; Simon et al., 2003; Vanin, 1985). These meiosis modifications include various mechanisms that maintain ploidy, sometimes coupled with a suppression of crossover recombination, and allow egg development without fertilization (Archetti, 2010; Bell, 1982). One modification that is often observed is a disruption of segregation by a suppression (or abortion) of the first meiotic division (Bell, 1982; Schurko et al., 2009b; Simon et al., 2003). Together with a suppression of crossover recombination, this meiosis modification leads to clonal offspring that are genetically undistinguishable from offspring produced by mitosis (Archetti, 2010; Boyer et al., 2021; Engelstädter, 2017; Moritz & Haberl, 1994; Rey et al., 2011).

Given the central position of meiosis in sexual life cycles, it is often hypothesized that OP evolves through mutations occurring in one or several meiosis genes (Hofstatter et al., 2020; Schurko & Logsdon, 2008; Simon et al., 2003). This hypothesis has led to the definition of a “meiosis toolkit,” a set of meiotic genes that are conserved among eukaryotes and whose function is restricted to meiosis (Hofstatter et al., 2020; Schurko & Logsdon, 2008). However, many OP species have partially asexual ancestors, in which OP may evolve by reusing the pathways for parthenogenetic reproduction that are already present in the partially asexual ancestors and by extending them to the entire life cycle (Simon et al., 2003; van der Kooi & Schwander, 2014). Considering this alternative mechanism for the evolution of OP is important because transitions to OP are particularly common in species with partial asexuality (Hebert, 1981; Kramer & Templeton, 2001; Simon et al., 2002) and because these species are often used as models to study the evolution of OP, by comparing them with their partially asexual sister clades (Schwander & Crespi, 2009; Stelzer, 2011).

A prominent example in which OP has evolved from partially asexual ancestors is the freshwater crustacean Daphnia pulex (more precisely, diploid lineages of a specific clade of this species complex, see Materials and Methods). Daphnia pulex has both cyclical parthenogenetic (CP) and OP lineages, with CP being ancestral to OP (Hebert et al., 1988, 1989; Tucker et al., 2013; Xu et al., 2013). Both CP and OP share a phase of subitaneous egg production (direct development of eggs, without diapause), during which females produce live-born offspring by parthenogenesis. Parthenogenesis occurs through a modified meiosis where the first division is aborted after homologous pairing and proceeds directly to the second division (Hiruta et al., 2010). This modified meiosis is nonreductional and nonrecombining, and thus leads to clonal offspring (Hebert & Crease, 1980, 1983). Parthenogenetically produced descendants of a given female thus constitute a clonal lineage (except for some rare recombination or gene conversion events, leading to some loss of heterozygosity, Flynn et al., 2017; Omilian et al., 2006). In both CP and OP, offspring can be males or females (because sex is environmentally determined), except in some particular, female-only lineages (so-called NMP for “non-male producing” lineages, as opposed to MP for “male-producing” lineages). Male production rates in male-producing OP and CP lineages can be vastly different. Male-producing CP lineages produce many males before engaging in sexual reproduction (Innes et al., 2000; Wolinska & Lively, 2008). In contrast, OP lineages usually produce almost only females and few so-called “rare OP males” (Hebert & Crease, 1983; Innes et al., 2000; Lynch, 1984; Wolinska & Lively 2008). CP and OP differ also in the mode of diapause egg production. In CP females (independently of whether they are MP or NMP), diapause egg production is sexual, whereas it is parthenogenetic in OP (Hebert, 1978; Hebert & Crease, 1980, 1983). The switch to diapause egg production typically concerns only some females in a population at any given time. The switch is induced by high density and environmental factors correlated with the end of the season favorable for adult growth (Fitzsimmons & Innes, 2006; Innes & Dunbrack, 1993; Innes & Singleton, 2000; Kleiven et al. 1992).

In OP D. pulex, a candidate gene of the meiosis toolkit, Rec8-B, has been identified with an upstream insertion and a frameshift mutation. This modified Rec8-B has been hypothesized to be the causal factor converting CP into OP (Eads et al., 2012). Indeed, from fission yeast to mammals, Rec8 is essential for reductional chromosome segregation during meiosis I (Lee et al., 2003; Watanabe & Nurse, 1999). In addition, at least in yeast, Arabidopsis, Caenorhabditis and mice, Rec8 also influences crossover recombination, with mutants showing suppressed or reduced crossover recombination (Lambing et al., 2020; Pasierbek et al., 2001; Xu et al., 2005; Yoon et al., 2016). Also note that, in mice, Rec8 plays a similar role in male and female meioses (Xu et al., 2005).

Eads et al. (2012) hypothesized that the modified Rec8-B allele has a repressing effect on the translation of other Rec8 mRNAs. This idea was inspired by the fact that the upstream insertion probably originated by a transposition of a mobile element, which might lead to Rec8 translational silencing via piRNA pathways. This hypothesis may also explain its dominant effect on the reproductive phenotype (the allele is heterozygous in all known OP lineages and absent from CP lineages). Yet, rare males produced by OP females are able to undergo reductional meiosis at least in some lineages and to successfully fertilize CP females (Innes & Hebert, 1988; Xu et al., 2015a) despite carrying the modified Rec8-B allele (OP males are genetically identical to OP females, being environmentally determined). These observations are difficult to reconcile with the Rec8-B hypothesis of OP evolution. Most importantly, because all ancestral Rec8 genes are expressed in both males and females (Eads et al., 2012; Schurko et al., 2009a) and because the modified Rec8-B allele is expressed in OP females (albeit at a low rate, Eads et al., 2012), at least one additional mutation would be required to silence it in males in order to achieve the “female-limited meiosis suppression” hypothesized by several authors (Hebert et al., 1988, 1989; Lynch et al., 2008; Schurko et al., 2009b; Eads et al., 2012; Ye et al., 2021a).

An alternative possibility is that the Rec8-B mutation affects only part of the OP phenotype. Indeed, at least three mutations are required in Arabidopsis to induce obligate asexuality (including one in Rec8, D’Erfurth et al., 2009) and previous results suggest that several genomic regions are involved also in D. pulex (Lynch et al., 2008). In the evolution of parthenogenesis, suppression of recombination is often observed as a different step, likely achieved secondarily by mechanisms that differ from those leading to nonreductional meiosis (Dukić et al., 2019; Hiruta et al., 2010; Lenormand et al., 2016). Recombination suppression and reductional meiosis can be largely separately regulated, which can be seen from the frequent occurrence of achiasmate species (Satomura et al., 2019), which have a perfectly reductional spermatogenesis in the absence of crossover recombination. The Rec8-B mutation may thus suppress recombination, whereas other mutations with female-limited expression may govern nonreductional meiosis during diapause egg production (e.g., those that govern parthenogenesis also in CP). Under this scenario, rare OP males are expected to undergo reductional spermatogenesis but have suppressed crossover recombination. However, whether or not recombination occurs during these male meioses is unknown.

In OP females, we know that recombination is absent or very low during parthenogenetic oogenesis of diapause eggs (Hebert & Crease, 1980, 1983). Yet, recombination rates during sexual oogenesis of diapause eggs in CP females are unknown. Indeed, recombination in CP has so far only been studied through sex-average and male-specific linkage maps (Cristescu et al., 2006; Xu et al., 2015a), but never specifically in females. Direct evidence for a difference in recombination rate during oogenesis of diapause eggs between CP and OP females is therefore still lacking. Furthermore, we do not know whether OP males recombine at similar or reduced rates compared to CP males. Here, we study the recombination rates during spermatogenesis in rare OP males and in CP males, as well as during oogenesis of diapause eggs in CP females (“recombination” referring to crossover recombination throughout, unless otherwise stated). We use these estimates to elucidate how recombination rate variation in males and females was affected by the CP to OP transition and discuss how these results inform different scenarios for the genetic underpinnings of OP evolution, including the involvement of Rec8.

To estimate these recombination rates, we performed two crosses to produce linkage maps of each of the four parents (one OP male, one CP male, and two CP females). We deliberately chose an OP lineage known to produce males whose spermatogenesis is reductional to specifically determine if the presence of the Rec8-B mutation may have other effects on meiosis, notably on recombination. In order to maximize the number of offspring, we used a mass mating approach with female-only CP lineages. Each cross involved crossing numerous females from a CP NMP lineage (a different clone in each of the two crosses) with either males from an OP lineage (OP x CP cross) or males from another CP lineage. Using Restriction-site Associated DNA sequencing (RAD-seq) on parents and offspring of each cross, we constructed highly saturated linkage maps and investigated recombination rate during gamete production in each of the four parents.

Materials and methods

Material

The four parent clones used for this study originate from three populations (called LPB, STM, and TEX) in temperate North America (Great Lakes, Quebec, see Supplementary Table S1). They belong to a panel of clones used in earlier phylogenomic and cytological studies (Tucker et al., 2013; Xu et al., 2015a). This is important, as the taxonomy of the group is complex and these earlier studies clearly show that these are diploid clones of the “Panarctic” clade (sometimes called American D. pulex, Colbourne et al., 1998, or Daphnia cf pulex sensu Hebert, 1995; Mergeay et al., 2006; Ye et al., 2021b), which is a distinct species from European D. pulex s.s and not even its sister clade (Colbourne et al., 1998; Dufresne, 2011; Taylor & Hebert, 1993). The diploid OP clones of this clade have a history of contagious asexuality, where new OP lineages arise by transmission of asexuality genes via rare males (Crease et al., 1989; Hebert et al., 1989, 1993; Paland et al., 2005), though a hybridization event with the sister clade D. pulicaria occurred around the time of the initial evolution of asexuality (Innes & Hebert, 1988; Xu et al., 2013). Other OP lineages exist in this species complex (Beaton & Hebert, 1988; Dufresne & Hebert, 1995; Ma et al., 2019; Mergeay et al., 2008; Vergilino et al., 2009). These are often polyploid, likely produced by hybridization among various taxa, and are not studied here. To be consistent with other recent literature on the clade (e.g., Maruki et al., 2022; Ye et al., 2021), we use the name D. pulex throughout, but note that our study strictly relates only to the diploid lineages of said Panarctic clade.

Using the four parental clones, we performed two mapping crosses. The first cross, “CP x CP,” was carried out using males of the CP clone TEX-1 and females of the CP clone LPB-87, while the second cross, “OP x CP,” was carried out using rare males of the OP clone STM-2 and females of the CP clone TEX-114. Both crosses were thus interpopulation crosses, and the fact that males of TEX-1 were used in one cross and females of TEX-114 in the other, allowed comparing male and female maps between clones from the same population. Both clones used as females (LPB-87 and TEX-114) are female-only(NMP) clones, that is, they are unable to produce males and thus participate in sexual reproduction only as females (Galimov et al., 2011; Innes & Dunbrack, 1993; Tessier & Cáceres, 2004; Ye et al., 2019). The use of NMP clones meant that mass-mating, needed to induce the production of ephippia (capsules containing diapause eggs), could be performed without occurrence of within-clone mating (i.e., all offspring were produced by obligate outcrossing between the two clones): to initiate a given cross, we introduced males of the father clone into a mass culture of the mother clone. Specifically, we regularly (about once every two weeks) introduced a small number of males into two 10 L aquaria containing mass cultures of females (one for each of the two crosses), across a period of six (CP x CP) to eight (OP x CP) months. In total, 165 males were used for the CP × CP cross and 299 males for the OP × CP cross. Both crosses produced several thousands of ephippia, which were collected and stored at 4 °C in the dark for two months or longer (necessary to break the diapause). Differences in male numbers used and in the duration of ephippia production were explained by the fact that many ephippia from the OP × CP cross were empty (i.e., did not contain any viable embryos) and because we wanted to ensure that we would be able to obtain a sufficient number of hatchlings for linkage analysis in each cross.

Hatching was induced by bathing ephippia for two hours in demineralized water (tap water purified by reverse osmosis), eight minutes exposure to a 2.6% bleach solution, followed by abundant rinsing with demineralized water (Paes et al., 2016; Retnaningdyah & Ebert, 2012). The ephippia were then exposed to high light for 24 hr and subsequently placed to standard laboratory conditions, in artificial Daphnia medium (Klüttgen et al., 1994), kept at 19 °C under a 16:8 hr light:dark photoperiod. The hatching vials were carefully inspected every two days for hatched juveniles, and any juvenile present was isolated individually in a new vial to initiate a clonal culture. We obtained a total of 104 clonal cultures of F1 offspring from the CP × CP cross (i.e., hatchlings that survived to adulthood and established a clonal culture by parthenogenesis). However, due to low hatching success, only 44 clonal cultures of F1 offspring of the OP × CP cross were obtained. All parent and offspring clones were kept under standard conditions in the laboratory, fed with 0.25 ml of a solution of microalgae (5 g/L of lyophilized Tetraselmis chuii and 500,000 cells/ml of Shellfish Diet 1800).

DNA extraction and RAD-sequencing

One batch (offspring clones) to three batches (parent clones) of 15–20 individuals of the same clone were collected, frozen in liquid nitrogen, and stored at −80 °C. Total genomic DNA was extracted from each batch using the DNeasy Blood & Tissue kit (Qiagen). DNA concentration and quality were examined by electrophoresis on 1% agarose gels and with a Qubit 3.0 (high sensitivity) fluorometer. The replicate batches of the parent clones were extracted and sequenced separately to increase sequencing depth (reads from all replicates of a given parent were pooled prior to analysis). Library construction was carried out according to the RAD-sequencing protocol described by Svendsen et al. (2015). The libraries were sequenced on four Illumina HiSeq2500 lanes, using 100 bp single-end sequencing by the Montpellier GenomiX platform (MGX, Montpellier, France).

SNP calling and filtering

Raw sequencing data were demultiplexed with Stacks v.2.41 (Catchen et al., 2013) using process_radtags. Reads were aligned to the D. pulex reference genome V1.1 (Colbourne et al., 2011) using BWA (version: bwa-0.7.17-r1188), and reads with a mapping quality of 30 or less were removed using Samtools v1.7 (Li et al., 2009). This procedure resulted in 5,217 to 4,695,427 reads per F1 of both crosses (Table S2). Even though most F1 were well-covered (83% of F1 had more than one million reads, Table S2), also low-coverage F1 were kept because the downstream analyses in Lep-Map3 (Rastas, 2017), specifically take genotype likelihoods into account, and removal of low-coverage samples is not recommended for these analyses. Parents were all highly covered with 3,381,813–6,000,981 reads per parent clone (all three replicates per parent clone combined).

The Stacks module “gstacks” with default parameters was used (--model marukilow and --var-alpha .05) to call SNPs and to infer genotype likelihoods. SNP markers were named according to their location, that is, scaffold name and base pair position in the reference genome. SNP markers were filtered using the module “population,” with 0.25 as the maximum proportion of missing values allowed per SNP marker across all F1 of a given cross. After this filtering step, 40,975 SNP markers were retained in the CP × CP cross and 41,917 SNP markers in the OP × CP cross.

Linkage maps construction and analysis

Linkage maps

Linkage maps were constructed using Lep-MAP3 (Rastas, 2017). Relationships between parents and offspring in each family were confirmed through the IBD (identity by descent) module in Lep-MAP3. The module “ParentCall2” was used to recall missing or erroneous parental genotypes based on genotype likelihoods of the offspring, as well as to remove noninformative markers (i.e., markers that were homozygous in both parents). The module “Filtering2” was used to remove strongly distorted markers (p-value < .0001, as recommended by the software for single-family data). These filtering steps reduced the numbers of retained markers to 25,951 and 32,654 for the CP × CP and the OP × CP cross, respectively. The stronger reduction in the number of markers in the CP × CP cross is explained by a higher proportion of distorted markers (21%) compared to the OP × CP cross (9%). Marker distortion (i.e., deviations from expected Mendelian segregation ratios among offspring) may represent genotyping errors (e.g., due to repetitive regions) or result from selection. Besides the high likelihood that many of them are erroneous, the exclusion of highly distorted markers is also motivated by fact that they are difficult to map (Hussain et al., 2017; Rastas, 2017).

The initial assignment of markers to linkage groups (LGs) followed the previous linkage map of D. pulex (Xu et al., 2015a), which was based on the same reference genome. Specifically, all markers on scaffolds that were present on the previous map, were assigned to the corresponding LGs of these scaffolds in the previous map. Second, we used the module “JoinSingles2All” to add markers on unmapped scaffolds (lodLimit=18). After the subsequent ordering steps, the initial assignment of markers to LGs was re-evaluated and corrected (if needed) using Lep-Anchor (see below). To order markers within each LG and to estimate linkage map distances, we used the module “OrderMarkers2.” The analyses were conducted separately for each parent of the two crosses using a pseudo-test cross design (Grattapaglia & Sederoff, 1994): SNPs that were heterozygous in both parents of a given cross (“ab × ab” SNPs) were used for the maps of both parents, while “ab × aa” and “aa × ab” SNPs (heterozygous only in the mother or only in the father) were used only for the maternal or paternal maps, respectively.

Finally, we used Lep-Anchor (Rastas, 2020) to detect potential assembly errors (“chimeric scaffold”), split them, if needed, and rerun the Lep-MAP3 pipeline using the split scaffolds. We ran three rounds of Lep-Anchor + Lep-MAP3 on the maps, until no further chimeric scaffolds were detected. This procedure identified 19 cases of likely assembly errors (assignment of parts of scaffolds to two distinct LGs or to different parts of the same LG, separated by a gap of at least 20 cM, Supplementary Table S3). The final maps were based on 15,577 SNPs for LPB-87 (female of the CP x CP cross), 13,733 SNPs for TEX-1 (male of the CP x CP cross), 16,492 SNPs for TEX-114 (female of the OP × CP cross), and 21,405 SNPs for STM-2 (male of the OP × CP cross).

Physical distances between markers

To estimate physical distances between markers, we performed a final ordering and orientation of scaffolds, using two additional rounds of Lep-Anchor + Lep-MAP3 with SNP markers from all four linkage maps combined. This resulted in a single ordering of the scaffolds containing at least one informative marker in at least one of the four maps. Using this ordering, we estimated physical distances (in bp) between markers, using a custom R script, assuming no gaps between adjacent scaffolds and forward orientation of scaffolds whose orientation could not be determined based on the information of the linkage maps.

Integrated linkage map

Based on the single physical ordering of the scaffolds, we also produced a single linkage map (“integrated linkage map”) using information of both crosses. First, a sex-averaged linkage map (using the option “sexAveraged=1” in the module “OrderMarkers2”) was produced for each of the two crosses. Second, these two sex-averaged maps were combined by averaging. Specifically, for each physical position, we estimated the cM position by linear extrapolation of the nearest markers in each sex-averaged map using a custom R script and averaged these values to obtain the integrated map.

Recombination rate

Genome-wide recombination rate (cM/Mb) was estimated by summing cumulative genetic lengths of all LGs and dividing them by the total length of the D. pulex genome, 197.3 Mb (Colbourne et al., 2011) or, alternatively, by the sum of the physical lengths of all anchored scaffolds 148.3Mb. An average recombination rate for each LG was estimated using the total genetic length of a given LG, divided by the sum of the physical lengths of the scaffolds anchored on that LG.

We also estimated the within-LG recombination parameter, $r¯intra$ (Veller et al., 2019), which, in addition to the number of crossover events also takes into account their locations to estimate the average amount of shuffling of genes that occurs within a chromosome per meiosis [central and widely-spaced crossovers generate more shuffling than tightly-spaced or terminal crossovers, (Veller et al., 2019)]. To estimate $r¯intra$ we used the MATLAB script from Veller et al. (2019), considering, as measure of physical length, the total length of anchored scaffolds. Following Veller et al. (2019), we also estimated $r¯inter$⁠, which is the probability of allele shuffling due to random assortment (i.e., segregation).

Comparison of recombination rate between maps

To visualize variation in recombination rates within LGs and to compare this variation among the different parents, we used Marey maps, which plot cumulative genetic distances (cM with respect to the first marker) against cumulative physical distances (Mb with respect to the first marker) for each marker of a given LG. The Marey maps were constructed using all markers. To quantitatively compare recombination rates between the four parents, we then used a subset of the data only (the “reduced dataset”), with truncated LGs in order to ensure identical terminal positions for all four maps. Specifically, the Mb position of the most interior terminal markers on any of the four maps was used (one at each LG end), and, in maps where no marker was present at that specific physical position, the cM position was estimated by linear extrapolation of the cM positions of the two nearest markers. The cM position of all markers was subsequently adjusted (by subtracting the cM position of the first marker) to ensure that the corrected cM-position of the first marker was zero. The reduced dataset included 2.3% fewer base pairs compared to the integrated map. This means that potential differences in recombination rates in the very terminal parts of LGs (where recombination is least efficient in shuffling genes) cannot be assessed with these data. We tested for differences in genetic lengths of LGs among the four parents, using ANOVA (based on the genetic length estimate of each LG in each parent in the reduced dataset and after verifying that the assumptions of ANOVA were met). We used pair-wise post-hoc tests with the adjusted Holm method to investigate pairwise differences between any pair of parents.

To investigate potential differences in the map length among the four parents at smaller scales, we divided each LG into three zones of equal length, each of them being composed of five windows (of equal length). Linkage map positions of the boundary positions of zones and windows were estimated for each map by linear extrapolation of the linkage map positions of the nearest markers. We first tested whether specific LGs showed differences in genetic length among the parents. Second, we investigated whether specific zones within LGs showed such differences. Finally, to test for differences in crossover occurrences, independently of the total map length of the LG, we normalized all four maps to the same genetic length. Using this normalized dataset, we again tested for differences among the four parents restricted to specific LGs or specific zones within LGs. Due to the frequent occurrence of windows without any crossover events, the ANOVA assumptions were not met. We thus analyzed the truncated data (both normalized and non-normalized) with pairwise, nonparametric tests: for each LG or each zone, we performed a Wilcoxon rank test (ZIW) modified for zero-inflated data (Wang et al., 2021). To test for effects of specific LGs, 72 (12 LGs*6 pairs) pairwise tests were performed. The effect of specific zones within LG was assessed through 216 (12 LGs*3 zones*6 pairs) pairwise tests, using windows as units of replication. The p-values were adjusted according to the Benjamini and Hochberg (1995) correction for multiple tests. Note that ten zones could not be tested because all windows within these zones showed zero recombination in the given pair.

Results

Comparison of recombination rate among maps

Overall genetic length

The linkage maps of all four parents, including the OP male, were highly similar (Figure 1, Table 1). This can be seen visually on the Marey maps (Figure 1), where flat, horizontal sections correspond to regions with low or no recombination and steep sections to regions with high recombination rates. The overall genetic length of a LG is represented by the cM position of the last marker (50 cM corresponding to one expected crossover), but quantitative comparisons of the overall genetic length between maps were based on the reduced data set (truncated to identical terminal positions). We found a slight but significant variation in genetic length among the four maps (ANOVA, F = 3.59, p = .02, Table 2), with only one of the pairwise post-hoc tests being significant (OP male vs. CP female LPB-87, p = .01, Table 2). Regarding sex-differences, the map length of the CP male (TEX-1) was slightly (6% on average) but nonsignificantly lower than the map lengths of the two CP females (Figure 2, Table 2, Supplementary Table S4). We thus find evidence for at most slight (and nonsignificant) heterochiasmy and recombination rates and patterns (Figure 1, Supplementary Figure S1) that are highly similar between males and females. Regarding the difference between CP and OP, the genetic length of the OP male was 11.9% lower than that of the CP male and 15.5% lower compared to the mean of the three CP parents, especially on the larger LGs (Supplementary Table S4, Figure 2). This subtle (though compared to the CP male nonsignificant) reduction in recombination in the OP male is not consistently found in all LGs (Supplementary Figure S1). Indeed, the ranking of the genetic lengths of LGs varied among parents (Figure 2), suggesting a possible interaction between LG length and parent. However, we were unable to test this interaction (because the length of each LG was only estimated once in each parent), and purely statistical effects (unplanned comparison or regression toward the mean) cannot be excluded.

Genetic length of specific LGs and zones within LGs

We tested whether the differences in genetic length among the four maps were driven by just some of the LGs or even more narrowly by just some zones within LGs. None of the LGs differed significantly among parents in any of the pairwise comparisons (after correcting for multiple testing; Supplementary Table S5). Only LG 9 showed a tendency for being shorter (in terms of genetic length) in the OP male, compared to each of the three CP parents (Figures 1 and 2, Supplementary Table S5). Two zones within LGs showed significantly different genetic lengths among maps (Supplementary Table S5): the middle zone of LG 7 was significantly longer (i.e., had more recombination, p_adj < .003) in the OP male compared to each of the three CP parents, and the middle zone of the LG9 showed significant differences (p_adj < .003) between most pairs, being shorter in the OP male than in the CP parents (Supplementary Table S5).

Normalized maps

We used the normalized dataset to test for differences in the localization of crossovers, independent of the length of the maps. Again, none of the LGs showed a significant difference in any of the pairwise comparisons and the only two zones that showed significant differences were the same ones already identified with the nonnormalized maps (Supplementary Table S5).

Integrated D. pulex map

Map characteristics

Our analysis resulted in an integrated D. pulex map (Supplementary Figure S2, Table 3) based on a single scaffold ordering and sex-average maps from both crosses (see Supplementary Table S6 for the corresponding detailed map characteristics of each individual map). The integrated map contains 345 of the 5,191 scaffolds of the Xu et al. (2015a) assembly (Supplementary Table S3). Note that the LG numbering is equivalent to the one in Xu et al. (2015a), but we added suffixes “_1”, “_2”, or “_3” for scaffolds that were split during the Lep-Anchor analysis (i.e., due to evidence that these likely are chimeric scaffolds). The total length of the 345 anchored scaffolds is 148.3 Mb (Table 3 and Supplementary Table S2), which represents 75.2% of the combined length of all scaffolds of the reference genome used here (Colbourne et al., 2011). The mapped physical lengths of LGs ranged from 5.7 Mb (LG 11) to 17.2 Mb (LG 10, Table 3). The four individual maps were on average 3.1% shorter than the integrated map, missing, on average, 50 (range 41 to 64), mostly smaller scaffolds (Supplementary Table S6, Supplementary File S1). Our integrated D. pulex linkage map contains a total of 48,615 SNP markers (Table 3), with an average intermarker distance of 0.02 cM (Table 3). The total map length is 1,174 cM with the different LGs spanning between 69.96 cM (LG 4) and 129.47 cM (LG 10, Table 3). The two sex-averaged Marey maps of each cross as well as the integrated Marey map are represented in Supplementary Figure S2.

Recombination rates

The estimated genome-wide recombination rate of the integrated map is 7.92 cM/Mb or 5.95 cM/Mb (ranging from 5.26 to 6.29 cM/Mb among the four linkage maps), depending on whether the total linkage map length was divided by the total length of anchored scaffolds or by the estimated total genome size (197.3 Mb) of D. pulex (Table 3). The genome-wide intrachromosomal recombination parameter $r¯intra$ across all LGs is .0164, while the interchromosomal recombination parameter $r¯inter$ is .45. Recombination rates of individual LGs varied between 6.2 cM/Mb (LG 3) and 12.8 cM/Mb (LG 6, Table 3, Supplementary Figure S2), and the intrachromosomal recombination parameter $r¯intra$ ranged between 3 × 10⁻⁴ (LG 11) and 3 × 10^–3 (LG 10, Table 3). The $r¯intra$ was positively correlated with the total genetic length (in cM) across LGs (Pearson r = .83, d.f. = 10, p = .0007) but negatively correlated with the recombination rate (in cM/Mb) (Spearman ρ = −.68, d.f. = 10, p = .01). As evident from the individual Marey maps (Figure 1), recombination rate varied extensively within LGs. In most LGs, we detected a large region with zero or almost zero recombination, putatively the peri-centromeric regions (Svendsen et al., 2015), although centromere locations are unknown in D. pulex. In contrast, recombination rates were high towards the ends of the LGs (Figure 1 and Fig. S2).

Discussion

No recombination suppression in OP males

The main goal of this study was to examine how recombination rate in both sexes was affected by the CP to OP transition. Our results demonstrate that recombination is not absent in OP males. Rather, the OP male showed similar levels of recombination compared to both the CP male and the CP females. While recombination rate was slightly lower than in CP, this effect was mainly local, being largely explained by LG 9, which corresponds to one of the regions that affect asexuality itself (Eads et al., 2012; Lynch et al., 2008; Tucker et al., 2013; Xu et al., 2013). The asexuality-determining regions are highly heterozygous in OP, because they introgressed from a related species (Heier & Dudycha, 2009; Tucker et al., 2013; Xu et al., 2013, 2015b). High heterozygosity (i.e., high divergence between homologs) may cause a local reduction in recombination, as demonstrated in other species (Lukacsovich & Waldman, 1999). Overall, our results indicate that OP males can be fully functional in the sense of producing sperm by canonical meiosis with at most subtly reduced recombination compared to CP males (see below for a discussion of reduced fertility). This contrasts with OP females, in which diapause egg production is clonal (or nearly so), based on the nonsegregation of allozymes and absence of loss of heterozygosity (Hebert & Crease, 1980; Hebert et al., 1989; Innes & Hebert, 1988). Similarly, recombination and segregation are absent (or present only at extremely low levels, 0.5% loss of heterozygosity per generation, including gene conversion events, Omilian et al., 2006) during the subitaneous phase of CP and OP females (Flynn et al., 2017; Keith et al., 2016; Omilian et al., 2006; Xu et al., 2011). Note that OP females cannot be crossed and were therefore not included in the present study. Overall, the presence of recombination in OP males and during oogenesis of diapause eggs in CP females but not in OP females shows that recombination suppression in the CP to OP transition is sex-specific. It also indicates that the only major change in recombination rate that occurred during the transition concerns the oogenesis of diapause eggs, which evolved suppressed recombination alongside nonreductional meiosis in OP females.

Because we concentrated on a single, diploid OP lineage, we cannot exclude that males of other OP lineages do show suppressed recombination (especially those producing unreduced sperm, Xu et al., 2015b). Yet, just one OP lineage is sufficient to demonstrate that male recombination suppression is not shared and therefore is unlikely to be ancestral among all OP lineages of this clade. The OP males used here were previously known to exclusively produce haploid sperm (Xu et al., 2015b). It is therefore unlikely that the higher number of empty ephippia and lower hatching rate in the OP x CP cross compared to the CP x CP cross are explained by the presence of unreduced sperm causing low fertility or by triploid or aneuploid offspring with low viability. Nonetheless, in our design, we can only estimate recombination in those offspring that hatched and we cannot formally exclude that nonviable offspring were fathered by sperm with suppressed recombination. However, low fertilization and hatching success have a number of plausible alternative explanations, including a generally decreased performance of rare OP males (independent of meiosis and recombination), hybrid incompatibilities, and suboptimal environmental conditions to break diapause. The two latter hypotheses may explain why hatching rate was low (albeit somewhat higher) also in offspring of the CP x CP cross, which was also an interpopulation cross.

Scenarios for the evolution of OP in D. pulex

Our finding that male recombination suppression is not shared among all OP lineages together with earlier findings showing the presence of reductional meiosis in males of the same OP lineage (Hebert & Crease, 1983; Innes & Hebert, 1988; Xu et al., 2013, 2015b) suggest that the modified Rec8-B allele does neither lead to suppressed recombination nor to nonreductional meiosis in males. The results therefore indicate that the causal genetic differences between CP and OP have female-limited effects. As mentioned above, this indicates that OP evolution in D. pulex is likely explained by a more complicated, multigenic scenario rather than by the Rec8-B insertion alone.

One possibility is that the modified Rec8-B allele is the key OP mutation, but another mutation restricts its expression to females. It is unclear why such silencing of the modified Rec8-B allele would evolve and immediately be beneficial in OP, but it might have been favored secondarily to promote contagious asexuality (by allowing OP males to produce reduced sperm and to transmit asexuality). This intriguing possibility might explain the variation in male meiosis phenotypes among OP lineages (Lynch et al., 2008; Xu et al., 2015b), with meiotic defects potentially restricted to lineages that do not carry the secondary mutation silencing the expression of the modified Rec8-B. The drawback of this scenario is that it relies on relatively weak selection as events of contagious asexuality are likely to be rare. The complexity of the scenario is, however, consistent with the identification of multiple chromosome regions potentially being implicated in OP determination (Lynch et al., 2008; Ye et al., 2021a). Furthermore, additional potential candidate genes have recently been identified, including genes with allele-specific expression (Ye et al., 2021a) and with general expression differences compared to the two parental species (Xu et al., 2022). A difficulty in investigating OP evolution is that causal genetic changes are not easily distinguishable from secondary ones. Indeed, a secondary loss of male functions in OP may be expected due to a relaxation of selection pressure (Innes et al., 2000; van der Kooi & Schwander, 2014; Wolinska & Lively, 2008). It may therefore be interesting to investigate whether OP lineages with meiotic defects are derived or ancestral to lineages that produce haploid sperm (Xu et al., 2015b). Loss of function mutations may occur secondarily also in other genes that are no longer under strong selection pressure (Normark et al., 2003), and allele-specific expression is predicted to increase secondarily in clonal OP according to recent population genetic theory (Fyon & Lenormand, 2018).

Another possibility is that the modified Rec8-B allele is noncausally related to the OP phenotype. OP might have evolved by reusing the subitaneous parthenogenesis pathway already present in the ancestral CP and extending it to oogenesis of diapause eggs. Interestingly, during subitaneous oogenesis, all Rec8 paralogs are transcribed at high levels (Eads et al., 2012; Schurko et al., 2009a), suggesting that it is possible to achieve a nonreductional meiosis without altering them (of course it is possible that they are all repressed at the translational level by an unknown repressor). As parthenogenesis in CP is specific to subitaneous oogenesis, there is no a priori reason to believe that spermatogenesis should be affected by the extension of parthenogenesis to oogenesis of diapause eggs in OP. This scenario requires very few mutational steps, since a single regulatory switch to the subitaneous parthenogenetic pathway may lead to asexual production of diapause eggs, and therefore OP.

Because this modification is likely to be minor compared to evolving an entirely new parthenogenesis pathway, it may be a common route to evolve OP in Daphnia and other CP-OP systems. In aphids, OP has evolved through a genetic change that prevents individuals from entering the diapause egg production in temperate regions with mild winters (Dedryver et al., 2013; Simon et al., 2002, 2010). The identified candidate region in the pea aphid (Acyrthosiphon pisum) contains genes involved in photoperiod sensitivity (Jaquiéry et al., 2014). Similarly, in rotifers, OP evolution is thought to be caused by a genetic change that prevents individuals from responding to chemical signals that induce sexual reproduction in CP (Stelzer, 2008; Stelzer et al., 2010). In contrast to OP aphids and rotifers, OP Daphnia still produce diapause eggs, and the main difference between CP and OP is in how diapause eggs are produced. It therefore seems unlikely that the mechanism involves an altered sensitivity to environmental signals. Nonetheless, the general principle may be similar. Rather than evolving a distinct parthenogenetic pathway restricted to oogenesis of diapause eggs, the transition to OP may be accomplished either by secondarily evolving a suppression of the part of the life cycle involving sex (such as in aphids and rotifers) or by adopting the existing parthenogenesis pathways to other parts of the life cycle. In line with this idea, gene expression is very similar between subitaneous oogenesis in CP and oogenesis of diapause eggs in OP D. pulex (Schurko et al., 2009a). Importantly, such transition may occur due to genes entirely unrelated to recombination or meiosis, suggesting that the meiosis toolkit may not be the only potential candidates.

In this second scenario, the Rec8-B mutation would just be a secondary, inconsequential and fortuitous loss of function which is difficult to reconcile with the facts that (a) the mutation occurs in all extant OP lineages (i.e., it must already have been present in their most recent common ancestor), that (b) Rec8 repression is known to affect major processes involved in parthenogenetic meiosis modifications (homologous recombination, sister chromatid orientation, and cohesion at meiosis I, Lee et al., 2003; Pasierbek et al., 2001; Watanabe & Nurse, 1999; Xu et al., 2005), and that (c) a plausible scenario exists for dominant Rec8 repression by the modified Rec8-B allele. Given the difficulties with each of the two scenarios, perhaps the most likely explanation is that OP evolution in D. pulex is explained by a combination of both scenarios, such as initial re-use of the subitaneous parthenogenesis pathway for the oogenesis of diapause eggs, followed by the Rec8-B mutation, which may have been positively selected to finetune the OP phenotype in some unknown way. A direct involvement of the Rec8-B mutation in the hypothesized pathway-change for diapause egg production seems unlikely, given the molecular functions of Rec8 and the fact that these functions are similar in male and female meiosis (Eads et al., 2012; Schurko et al., 2009a). Obviously, these points remain highly speculative and will have to be elucidated by future studies.

No striking heterochiasmy in CP D. pulex

Heterochiasmy (higher recombination rate in one sex than in the other, suggesting at least partially different recombination pathways between sexes) is commonly observed in many species (Burt et al., 1991; Lenormand, 2003; Lenormand & Dutheil, 2005). The second main result, the comparison of sex-specific linkage maps of CP D. pulex, revealed slightly reduced recombination in the CP male compared to both CP females. While this may indicate slight, female-biased heterochiasmy, the effect was nonsignificant and nonconsistent across LGs. Similar results have recently been obtained in D. pulicaria, the sister clade of the lineages used here: in D. pulicaria, slight male-biased heterochiasmy was found overall, though also with variable effects among different LGs (Wersebe et al., 2022). The most informative contrast for our heterochiasmy analysis is the comparison of the CP male (TEX-1) with the CP female from the same population (TEX-114). Indeed, LPB-87 has a nonrecombining region on LG6 that is not shared by any of the other parents and points to the possibility of population differences in recombination. The use of NMP females may also have affected local recombination patterns, as we cannot exclude that incompatibilities occurred between genes of the fathers and genes in the NMP-determining region on LG 1 (homologous to the NMP-determining region of the map used by Ye et al., 2019). Outside this region, MP and NMP lineages do not differ genetically (Ye et al., 2019). Thus, despite the fact that the male–female contrast was confounded in our design with the contrast between MP and NMP (both female clones were NMP), the absence of any striking heterochiasmy suggests that little if any influence of the use of NMP clones on recombination patterns. Overall, our results together with those on D. pulicaria (Wersebe et al., 2022) suggest that recombination does not differ strongly between males and females of this species complex.

A new reference map for D. pulex

Our sex-specific and integrated maps constitute important additions to existing genomic resources for D. pulex. Previously, two D. pulex linkage maps had been published, a sex-average map based on microsatellite data (Cristescu et al., 2006) and a male-specific map, based on single sperm genotyping (Xu et al., 2015a). A third map, which was published as an annex of a new genome paper (Ye et al., 2017), is likely erroneous, as it predicts, on average, over eight crossovers per LG and meiosis, as opposed to a bit over two in our map and that of Xu et al. (2015a). We therefore compare our results, mainly to those of Xu et al. (2015a), which were based on the same genome assembly as ours (Colbourne et al., 2011). The Xu et al. (2015a) map anchored 187 scaffolds (131.9 Mb) with an average intermarker distance of 0.87 cM, while our integrated map anchors 345 scaffolds (148.3 Mb) with 0.02 cM between markers, on average. The main improvement thus comes from the mapping of many additional, mostly smaller scaffolds. Furthermore, while there was a high degree of collinearity between the maps, we corrected 19 likely assembly errors (chimeric scaffolds) and placed the part-scaffolds back to the linkage map. Still, about one fourth of the total assembly (197 Mb) remains unmapped, either due to smaller scaffolds containing no SNPs or filtered scaffolds due to a low mapping score, and perhaps also due to the presence of contaminant scaffolds in the reference genome.

Regarding the genome-wide recombination rate, the estimates from our study and that of Xu et al. (2015a) are very similar (7.9 cM/Mb and 7.3 cM/Mb, respectively). These estimates are also similar to those from other Daphnia species (D. pulicaria, 7.4 cM/Mb, Wersebe et al. (2022) and D. magna, 6.8 cM/Mb, Dukić et al. (2016)), suggesting conservation of recombination rates in the genus.

Regarding the individual maps, the main differences among parents (a part from the few, region-specific ones discussed further up) appear to be due to a small group of terminal markers per LG. These differences may, however, be artefactual because the estimation of recombination rate is less reliable for terminal markers. Indeed, to counter the well-known fact that erroneous genotype information artificially increases recombination rate, Lep-MAP3 (Rastas, 2017) uses information on several flanking markers to smoothen spikes in apparent recombination rates due to unreliable markers. As a consequence, the estimated total genetic length of LGs may be rather sensitive to these errors in terminal parts (where they are less easily smoothened) as well as to the inclusion or not of an additional terminal marker (because most LGs exhibited higher recombination rates in peripheral parts than in central ones).

The high prevalence of peripheral crossovers likely has also contributed to the observed low $r¯intra$ (within-LG recombination parameter) because terminal recombination contributes only little to effective gene shuffling. The excess of recombination in peripheral parts was mainly noted in (physically) larger LGs, a pattern also observed in many other animal and plant species (Haenel et al., 2018). This pattern might amplify the well-known negative relationship between the recombination rate (cM/Mb) and the physical size of LGs, caused by the constraint of at least one crossover per LG and meiosis (Hunter, 2007; Mather, 1938). It might thus also contribute to the observed positive and negative correlations of $r¯intra$ with cM length and cM/Mb recombination rate across LGs, which are likely explained by the same factors.

Conclusion

We found that the CP to OP transition in D. pulex, which involves an almost complete suppression of recombination in females, did not lead to a substantial reduction in male recombination, and that recombination is very similar between male and female CP. These findings may be explained by the evolution of OP either through a restriction of effects of the Rec8-B mutation to females, through an extension of the subitaneous parthenogenesis pathway to oogenesis of diapause eggs, or through a combination of the two. A re-use of pre-existing parthenogenesis pathways may be a common way to evolve OP in species with mixed sex–asex reproductive systems, suggesting that genes of the meiosis toolkit may not necessarily be the only candidates for the underlying changes.

Data availability

Linkage maps including all anchored markers and their positions are given in File S1. Scripts are available at the URL https://github.com/ces940208/Scripts. RAD-seq data on parents and offspring of the two crosses are available on NCBI SRA database under the BioProject accession number PRJNA967230.

Author contributions

C.H., C.M., and T.L. conceived the study. C.M. conducted the analysis.

Conflict of interest: The authors declare no conflict of interest.

Acknowledgments

The authors are grateful to Michael Lynch for kindly providing the D. pulex strains used. We thank Roula Zahab, Marie-Pierre Dubois and the GEMEX platform for help with the construction of the libraries. We thank the MGX-Montpellier GenomiX platform for the sequencing. We thank Pasi Rastas and Enrique Ortega-Abboud for helpful comments during the maps’ construction. This work was funded by grant ANR-17-CE02-0016-01, GENASEX, from the Agence Nationale de la Recherche.

References