Phenotypic plasticity rather than genotype drives reproductive choices in Hydra populations

Facultative clonality is associated with complex life cycles where sexual and asexual forms can be exposed to contrasting selection pressures. Facultatively clonal animals often have distinct developmental capabilities that depend on reproductive mode (e.g., negligible senescence and exceptional regeneration ability in asexual individuals, which are lacking in sexual individuals). Understanding how these differences in life history strategies evolved is hampered by limited knowledge of the population structure underlying sexual and asexual forms in nature. Here we studied genetic differentiation of coexisting sexual and asexual Hydra oligactis polyps, a freshwater cnidarian where reproductive mode‐dependent life history patterns are observed. We collected asexual and sexual polyps from 13 Central European water bodies and used restriction‐site associated DNA sequencing to infer population structure. We detected high relatedness among populations and signs that hydras might spread with resting eggs through zoochory. We found no genetic structure with respect to mode of reproduction (asexual vs. sexual). On the other hand, clear evidence was found for phenotypic plasticity in mode of reproduction, as polyps inferred to be clones differed in reproductive mode. Moreover, we detected two cases of apparent sex change (males and females found within the same clonal lineages) in this species with supposedly stable sexes. Our study describes population genetic structure in Hydra for the first time, highlights the role of phenotypic plasticity in generating patterns of life history variation, and contributes to understanding the evolution of reproductive mode‐dependent life history variation in coexisting asexual and sexual forms.


| INTRODUC TI ON
Sexual reproduction is the most widespread mode of reproduction in multicellular eukaryotic species as more than 99% of these organisms reproduce sexually (Bell, 1982). Sexual populations are expected to have an advantage in heterogeneous environments due to more frequent genetic reshuffling (Peck et al., 1999) and can take advantage of recombination to enhance the spread of favourable mutations, thus allowing them to survive in a broader range of ecological conditions (Pound et al., 2004). Asexual (clonal) reproduction, where offspring are derived from a single parent, is by comparison less common, but has been described in most major taxonomic groups, except mammals and birds (Avise et al., 1992).
Asexual populations enjoy the benefit of a higher multiplication rate because they do not pay the twofold cost of sex (Maynard-Smith, 1978), although they are thought to be less resilient towards changes in environmental conditions (Haldane, 1932). Some organisms are thought to have reproduced for millions of years without sex (Welch & Meselson, 2000), although more and more cases of rare sexual reproduction are observed in species that were previously considered completely asexual (Maynard-Smith et al., 1993;Tibayrenc & Ayala, 2002). Consequently, most species where clonal reproduction is present today can be considered facultatively sexual (or facultatively asexual).
The presence of both sexual and asexual forms of reproduction (i.e., facultative asexuality) within the life cycle has profound effects on the demography, genetics and life history of these organisms (Halkett et al., 2005). Asexual and sexual reproduction are considered to be distinct strategies that, on average, will have different fitness outcomes depending on environmental conditions, and hence selection is expected to shape them differently (e.g., Gardner & Mangel, 1997). Accordingly, life history strategies differ markedly depending on the mode of reproduction in many facultatively clonal species. Asexual organisms ranging from plants to sponges, cnidarians, flatworms, polychaetes and ascidians often have higher levels of somatic maintenance, lower senescence and considerably extended lifespans (summarized in Nilsson Sköld & Obst, 2011, but see Martínez and Levinton, 1992), compared to sexual individuals belonging to the same species or closely related taxa. This is perhaps best seen in cnidarians (Yoshida et al., 2006), planarians (Baguñà, 1998) and annelids (Zattara & Bely, 2016), where asexual individuals are capable of extensive regeneration of lost body parts and seem to avoid age-related declines thanks to their pluripotent adult stem cells (Aboobaker, 2011;Elliott & Sánchez Alvarado, 2013;Newmark & Alvarado, 2002;Rink, 2013;Tan et al., 2012;Valenzano et al., 2017). Conversely, individuals/species that follow a sexual strategy or switch from clonal to sexual reproduction often lose regeneration ability and have increased rates of ageing (Baguñà et al., 1999;Kobayashi & Hoshi, 2002;Krois et al., 2013;Tan et al., 2012). For instance, in the freshwater cnidarian Hydra oligactis, asexual individuals can be maintained in the laboratory for years without agespecific increases in mortality or declines in fecundity (Brien, 1953;Martínez et al., 2010;Tomczyk et al., 2015). However, switching from asexual to sexual reproduction results in a senescence-like process that consists of depletion of adult stem cells, loss of regeneration ability and neurogenesis, decrease in body size, and ultimately an increased rate of mortality (Brien, 1953;Sebestyén et al., 2018;Tomczyk et al., 2017Tomczyk et al., , 2019Yoshida et al., 2006).
Explaining how reproductive mode-dependent life history strategies evolved is relatively easy when sexual and asexual strategies occur in allopatry and inhabit different environments. For instance, in small freshwater animals capable of clonal reproduction, asexual and sexual strategies are often genetically fixed (Ament-Velásquez et al., 2016;Pongratz et al., 1998Pongratz et al., , 2003Simon et al., 2003), with asexual lineages often restricted to marginal (Peck et al., 1998) or ephemeral habitats (Dudycha & Hassel, 2013), and occupying distribution areas at higher latitudes than their sexual relatives (Hörandl, 2009).
The evolution of distinct life histories is more puzzling when alternative strategies coexist in the same environment at the same time.
When genetically distinct sexual and asexual lineages co-occur, theory predicts that asexual lineages should be eliminated because they are slow at integrating favourable mutations but fast at accumulating deleterious ones (Hadany & Beker, 2003;Kondrashov, 1988;Muller, 1964). However, in some cases sexual and asexual lineages occurring in sympatry do not compete with each other for the same habitats or resources but rather specialize on distinct ecological niches (Barraclough et al., 2003;Maynard-Smith, 1978;Vrijenhoek, 1984), resulting in the stable coexistence of the two forms. Such specialization can lead to speciation and generation of cryptic speciesmorphologically similar but genetically more or less distinct entities (Birky & Barraclough, 2009;Mayr, 1948;Peccoud et al., 2009)-with reproductive isolation further reinforced by differences in reproductive strategies among sexual and asexual lineages. Overall, this will result in genetically distinct sexual and asexual lineages specializing on different niches, having distinct life histories and possibly on the path of speciation (an example of this occurs in planarian species; Leria et al., 2020).
Facultative clonality usually occurs in highly variable environments, where asexual and sexual modes of reproduction are observed depending on environmental conditions. In such facultative clonal organisms, the presence of both asexual and sexual strategies within the same population is the result of phenotypic plasticity in mode of reproduction. Phenotypic plasticity can be defined as the ability of a genotype to express different phenotypes when exposed to different environmental conditions (Pigliucci et al., 2006), but the internal state of an individual (e.g., condition, the state of immune system, age) may also play an important role in the formation of phenotype (Hadany & Otto, 2007). Phenotypic plasticity in mode of reproduction can occur either when (i) the same individuals switch between sexual/asexual strategies during their lifetime (i.e., individual plasticity, such as seen in cyclical parthenogens), (ii) a clonal lineage expresses sexual/asexual strategies in different individuals of the clone (i.e., clonal plasticity) or (iii) a combination of these two, when individuals within a clone switch between reproductive modes but vary in their propensity to do so. We can find good examples of cyclical parthenogenesis (i) in aphids (Moran, 1992), in naidine annelids (Learner et al., 1978) and in cladocerans (Decaestecker et al., 2009). Clonal plasticity (iii) occurs in planarians (Pongratz et al., 2003) and cladocerans (Dudycha & Hassel, 2013). The third category (iii) may include some cladocerans (Decaestecker et al., 2009), but there might be many organisms in this group (possibly including hydras as well; . Sexual reproduction in these species is often associated with unfavourable and stressful conditions. Indeed, theory suggests that in facultatively clonal organisms the mode of reproduction should depend both on ecological conditions (Gardner & Mangel, 1997;Sakai, 1995) and intrinsic physiological state (Hadany & Otto, 2007), such that high-condition individuals facing favourable conditions (e.g., high food availability, low competition, few stressors) should reproduce clonally, while otherwise sex should prevail (Gardner & Mangel, 1997;Hadany & Otto, 2007;Sakai, 1995). Hence, reproductive mode and the associated life history might represent a plastic response to favourable or stressful environments, as these distinct conditions could occur simultaneously within the same environments in the form of microhabitats.
In this study, we aimed to understand population genetic structure underlying variation in reproductive mode within Central European H. oligactis populations. H. oligactis displays marked reproductive mode-dependent variation in life history. This species reproduces asexually under favourable conditions (warm conditions) and does not show signs of senescence, but in unfavourable conditions (cooling), polyps start sexual reproduction (Reisa, 1973), then undergo post-reproductive senescence and many of them die within a few months after initiating sexual reproduction (Yoshida et al., 2006). The level of sexuality could be genetically determined in this species as well, because, when kept under standard conditions in the laboratory, H. oligactis strains express differences in the probability of initiating sexual reproduction and in post-sexual survival rates Tomczyk et al., 2015). However, actual data on population structure in this species (or any other Hydra species), to our knowledge, does not exist and this lack of data precludes any inference on the interrelationship of different reproductive mode strategies. While several phylogenetic studies in Hydra species have been carried out over the last two decades, these studies employed markers that do not provide sufficient resolution to explore genetic patterns within Hydra populations (Kawaida et al., 2010;Martínez et al., 2010;Schwentner & Bosch, 2015). Notably, however, one of these studies (Schwentner & Bosch, 2015) suggested that several cryptic species might exist within the genus Hydra, raising the possibility that asexual and sexual forms of H. oligactis coexisting in the same population might be distinct cryptic species.
To describe population genetic structure in H. oligactis we employed restriction-site associated DNA sequencing (RAD-Seq, Andrews et al., 2016;Baird et al., 2008;Davey & Blaxter, 2010). The use of RAD-Seq in population genomics is becoming increasingly popular and it has been successfully applied in cnidarians to infer phylogeographical patterns in sea anemones (Reitzel et al., 2013), chimerism among colonial hydrozoans (Chang et al., 2018), and other phenomena, such as identification of introgressive hybridization (Combosch & Vollmer, 2015) and sex determination (Pratlong et al., 2017) in corals. Specifically, we aimed to obtain detailed genetic data on population structure of this model system and thus to obtain more insight on the evolution of distinct reproductive strategies.
We asked whether differences in reproductive mode of H. oligactis polyps inhabiting the same population are caused by genetic differentiation or phenotypic plasticity. This could work in the following possible ways. (i) If there is notable genetic differentiation among reproductive mode categories (asexual and sexual), we expect that sexually reproducing individuals would cluster into different groups from asexual ones. (ii) If phenotypic plasticity plays a major role in shaping the mode of reproduction, individuals of the same genotype may follow different reproduction modes. (iii) Alternatively, both genetic differentiation and phenotypic plasticity might be involved; in this case we expect variable proportions of sexual/asexual individuals within distinct genotypes.

| Study system
Hydra oligactis (Pallas, 1766) is a small freshwater invertebrate of the phylum Cnidaria (Class Hydrozoa). They are sedentary predators of small zooplankton in the Northern Temperate zone. H. oligactis is widespread in both Eurasia and North America (Martínez et al., 2010) and is commonly found in various aquatic habitats (Holstein, 1995;Schuchert, 2010), but prefers colder and larger/deeper water bodies, due to its attenuated heat shock response (Brennecke et al.,). The life cycle of H. oligactis includes asexual (budding) and sexual reproduction (Reisa, 1973). The asexual phase generally occurs in spring and early summer, while sexual polyps can be found in late summer and autumn (Reisa, 1973), because strong temperature decrease in autumn induces sex (Lenhoff, 1983;Littlefield et al., 1991;Reisa, 1973). Sexual reproduction consists of the development of gonads (testes in males, eggs in females; i.e., this species has separate sexes).
Following fertilization, the zygote detaches from the parent animal and after hatching the asexual life phase starts again.

| Study populations and sampling of Hydra polyps
Sampling was focused on our long-term study site near Tiszadorogma  (Table 1). From each population, we attempted to collect polyps from multiple locations at the same sampling site (with a distance between locations of at least 2 m) to increase the chance of finding genetically different individuals within the populations (because hydras can reproduce asexually and polyps close to each other could be more likely to belong to the same genetic line). We also recorded the GPS coordinates of each collection point (Table S1). Hydra polyps were collected from free-floating and submerged macrophytes (most often based on morphology-tentacle length/body length, the presence of stalk, tentacle formation in buds (Schuchert, 2010). Next, the sexual state of animals was determined-female (polyp with differentiated eggs), male (polyp with differentiated testes) or immature (polyps that have clearly visible gonads, but these gonads are in early stages of development and hence sex could not be determined; Figure 1).
We selected 120 individuals from all samples of collected Eppendorf tubes to cover asexual, male, female and immature individuals simultaneously in all populations, by picking one of these types from each collection point where they were present. In populations where only asexual individuals were found, we randomly selected 3-5 asexual individuals. More than one third of the samples (38%) were selected from our focal population (M28).

| Drying and DNA extraction
Hydra oligactis polyps were dried using silica gel and stored at room temperature to preserve DNA quality. This low-cost and effective drying agent was successfully used by previous authors for plants (Chase & Hills, 1991), but also for insect (Hackett et al., 2000;Post et al., 1993) and mammal samples (Cserkész et al., 2016(Cserkész et al., , 2017. In practice, this was done by putting the hydra polyps one by one inside sterile pipette tips, and tips with polyps were individually placed in plastic bags containing about 6 g of dry silica gel. Genomic DNA from collected H. oligactis individuals was isolated using a standard mammalian nucleic acid extraction protocol   (Cserkész et al., 2016), details of which can be found in the Supporting Information (Methods S1).

| RAD-Seq library preparation and sequencing
All samples were prepared separately in three RAD libraries. Details of the library preparation protocol can be found in Methods S1. The quality and quantity of the library were checked with a Bioanalyzer (

| Sequence processing
Raw Illumina reads were processed using the stacks process_radtags pipeline (Catchen et al., 2013). First, data were demultiplexed and reads with adapter contamination were removed (allowing an adapter mismatch =2). We used process_radtags to remove reads with uncalled bases (-c) and low-quality scores (-q). Reads with barcodes having a single mismatch were rescued (-r). We also trimmed reads to a uniform length of 140 nt. For the forward reads, this was done by removing six nucleotides (the inline barcodes) from the start and four low-quality nucleotides from the end of the reads. For the reverse reads, we removed the first nucleotide from the beginning of the reads (this nucleotide is invariable and was added during library preparation) with trimmomatic (version 0.36; Bolger et al., 2014) and nine low-quality nucleotides from the end of the reads using the built-in trimming functions of process_radtags.
We also performed in silico species identification (mapping de-

| stacks parameter choice and estimation of error rates
The choice of parameters used for assembly of RAD loci can greatly affect the outcome of the analysis by affecting how loci and single nucleotide polymorphisms (SNPs) are identified within and across samples. Furthermore, assembly options can influence the rate of genotyping errors. Such genotyping errors arise due to unequal shearing, PCR errors and allele dropout (Mastretta-Yanes et al., 2015). Given these sources of error, a useful strategy is to sequence technical replicates from the same sample and to look across the assembly parameter space to find parameter combinations that maximize the number of polymorphic loci while also minimizing the ratio of SNPs that are erroneously identified in technical replicates. To do so, we randomly selected N = 16 individuals for resequencing (one from each population, with the exception of population M28, where four samples were resequenced).
Using this data set, we ran the stacks de novo assembly pipeline, varying the values of the minimum depth of coverage required to create a stack (-m; from 3 to 5), the number of mismatches allowed between stacks within individuals (-M; from 1 to 9) and the number of mismatches allowed between stacks between individuals (-n; from 1 to 9 at the same value as -M; Paris et al., 2017).
After running the pipeline with these values, we looked at the number of polymorphic loci shared by at least 80% of the sample set and estimated the homozygote and heterozygote allele error rate, using the software tiger (Bresadola et al., 2020). The final stacks de novo parameters were selected by aiming to maximize the number of widely shared loci while also minimizing the genotyping error rate. We also generated neighbour-joining tree of the replicates with R's ape (version 5.0) package (Paradis et al., 2004;R Core Team, 2020) to check that replicates correctly cluster together. To obtain the neighbour-joining tree we used the genetic distance matrix calculated using the dist.gene function in ape, with pairwise deletion of loci containing missing data. Finally, we also calculated the genetic distance of replicate pairs and plotted them on the histogram showing the distribution of between-individual genetic distances (all distances calculated using the dist.gene function; Figure S1). This allowed us to see the overall magnitude of sequencing/genotyping errors: if error rate is low, replicate pairs should have a low genetic distance relative to that of distinct genotypes (also see below section on the spectrum of genetic diversity).
After identifying the assembly parameters that maximize the number of polymorphic loci and minimize the error rate, we ran the de novo pipeline using the final parameters for the whole data set.
The resulting catalogue of loci was filtered using vcftools (Danecek et al., 2011), as follows: we required a minor allele count (--mac) of 3, a minimum genotype quality (--minGQ) of 20, a minimum depth for an allele count (--minDP) of 3 and a minor allele frequency (--maf) of 0.05. We also excluded loci with a mean depth greater than the 97.5 percentile in an attempt to remove potentially paralogous loci. The resulting set of loci was further restricted based on presence across the samples (--max-missing) by selecting the set of loci that was represented in at least 80% of the samples. To reduce bias due to linkage disequilibrium (LD) that arises when closely situated SNPs are analysed, we selected the first SNP from each RAD locus (using the -write-single-snp option in stacks).

| Clone detection and sibship reconstruction
We used two methods to infer clones and assign individuals within our sample to multilineage genotypes (MLGs). The first method relies on genetic similarity between individuals and uses the frequency distribution of pairwise genetic distances (the spectrum of genetic diversity; Rozenfeld et al., 2007) to identify clones. Theoretically, individuals belonging to the same MLGs should have zero genetic distance because they are derived from a single individual through asexual reproduction. In practice, however, genotyping errors and somatic mutations cause variation within MLGs that generates a distribution of genetic distances greater than zero (Kamvar et al., 2015). We plotted the spectrum of genetic diversity and assumed that the first peak of small genetic distances represents pairs of individuals belonging to the same MLGs, while the second peak of larger genetic distances represents variation between MLGs. We used cutoff_predictor from R's poppr package (Kamvar et al., 2014) to identify the cutoff value between the two peaks. The genetic distance matrix was obtained using the dist.gene function of R's ape (version 5.0) package (Paradis et al., 2004;R Core Team, 2020), with pairwise deletion of loci containing missing data.
While identifying clones by looking at the spectrum of genetic diversity can accurately assign most individuals to MLGs, it relies on an arbitrarily chosen threshold whose optimal value might depend on a number of factors, such as the number of markers used, or the mistyping rate (Wang, 2016). As a more objective way to infer clones, the threshold can be optimalized to take into account features of the data set at hand. We used the software colony (version 2.0.6.6; Jones & Wang, 2010) to infer clones using an optimalized threshold that takes into account mistyping rates, missing data and the number and allele frequencies of markers (Wang, 2016). The method implemented in colony uses a likelihood framework to assign individuals to candidate relationships of clone mates and close competitive relationships (e.g., full sibships) and has been shown to accurately identify individuals belonging to the same MLGs through simulations (Wang, 2016).
The likelihood-based approach implemented in colony version 2.0.6.6 (Jones & Wang, 2010) also identifies full-sib dyads, half-sib dyads or unrelated pairs by performing a joint inference of parentage and sibship (including clones) from multilocus genotypes, taking into account sequencing errors and uncertainties in SNP calling. All individuals were included as potential offspring in the analysis, because the presence of clonality in Hydra implies that generations can be overlapping and there is no unequivocal way to assign candidate parents.
These potential offspring were then assigned into clonal lineages and family clusters (groups of individuals linked through sibship). For the colony analysis, we used a full-likelihood-pair-likelihood score combined (FPLS) method, assumed a polygamous mating system for both parents and kept all other parameters at their default values.

| Genetic structure
To obtain a more fine-scale view of genetic structure within our sample that reflects genetic structure among unrelated individuals, in addition to the patterns detected by looking at closely related individuals, we employed two methods. First, we calculated pairwise kinship coefficients between samples, using the method described in Loiselle et al. (1995). were also calculated on this reduced data set (extracted from stacks output).

| Genetic structure vs. reproductive strategies
To test for the association between genetic structure and sexual propensity we counted the number of sexual (immatures, mature females and mature males) and nonsexual (asexuals and nonreproductives) individuals in the family clusters inferred in the colony analysis and tested whether these are identical in the different clusters by means of Fisher's exact test. Furthermore, we performed analysis of molecular variance (AMOVA) between sexual and asexual individuals. AMOVA was implemented in R's pegas (version 0.10) package (Paradis, 2010) using 1,000 permutations.

| Read statistics, species identification and decontamination
There were altogether 625.7 million raw paired-end reads (all descriptive statistics include replicates as well). From these raw reads, 94.2% were retained after filtering for low-quality reads, adapter contamination, ambiguous barcodes and ambiguous RAD-tags.
There were an average of 4.3 million reads per sample (range 0.8-

million).
Mapping reads from individual samples to the Hydra 2.0 genome showed a divergent distribution ( Figure S2), with N = 12 samples (10% of the total) showing a mapping rate >60%, with the remainder showing a mapping rate of <40%. We conclude that samples with >60% mapping rate are Hydra vulgaris individuals mistakenly identified as H. oligactis and therefore we excluded them from all subsequent analyses.
Running the stacks de novo pipeline with default settings identified 1.37 million RAD loci. In total, 57% of these loci showed no hit in the nt database, a further 19% mapped to cnidarian sequences.
The remainder mapped to other taxonomic groups and were filtered out to form a contaminants database. The top contaminants were Pseudomonadales and Burkholderiales ( Figure S3), two bacterial orders that are commonly found within the Hydra microbiome (Fraune et al., 2015). After removing presumed contaminant loci, the secondary GC peak diminished substantially ( Figure S3).
The proportion of read pairs identified to be PCR duplicates after removal of contaminants was 29.4%. Mean coverage after removing contaminants and PCR duplicates was 13.0× per sample (range: 4.0-25.4×).

| stacks parameter selection
The number of polymorphic loci shared by at least 80% of the technical replicates included in the analysis decreased as the minimum number of reads required to create a stack (-m) increased. While Allele error rates were ~0.01 for homozygote loci and ~0.04 for heterozygotes ( Figure S4). The higher heterozygote allele error rates probably reflect the fact that allelic dropout has a high contribution to the overall genotyping error in this data set (Bresadola et al., 2020;Wang, 2004 Figure S4).
Taking into account these considerations, we chose m = 3 for the minimum number of reads required to create a stack, because at higher values of m, the number of widely shared polymorphic loci was lower and heterozygote allele error rate was higher. Furthermore, we chose M/n = 3, because the number of widely shared polymorphic loci did not increase, while error rates did not decrease considerably by setting these parameters to a higher value.
Running stacks de novo with these final settings, we obtained 9,844 RAD loci after filtering. Individual missingness was 13.7% (range: 2.4%-88.7%). Two individuals had individual missingness >60% (one from population M44 and one from M38), while the rest had missingness <50%. These two individuals could be either mistaken H. circumcincta individuals or samples with low DNA extraction efficiency. In support of the latter hypothesis these two samples had the lowest coverage in our samples set (4.0 and 4.5, respectively). We removed these two individuals and repeated the selection of loci without them to obtain the final SNP data set. This resulted in a total of 11,319 RAD loci, with an average individual missingness of 13.5% (range 2.7%-43.3%).

| Clone detection
The spectrum of genetic diversity showed a discontinuous distribution of genetic distances and cutoff_predictor estimated a genetic distance of 0.13 below which individuals can be considered to belong to the same multilocus genotypes. Based on this, we identified N = 61 multilocus genotypes in the set of N = 106 polyps included in the analysis (Table S2). The same set of MLGs was identified in colony. The number of polyps within MLGs ranged from one to eight and with one exception (population M67), all populations contained more than one clone. Clone mates were always located in the same population (i.e., none of the MLGs was observed in more than one population).
Technical replicate pairs had a mean genetic distance of 0.07 and with one exception (one out of 16, 6.25%), genetic distance was less than the MLG threshold value of 0.13 identified by cutoff_predictor ( Figure S1). The single exception was sample R12_1_2, where the genetic distance between the two replicates was 0.25. This outlier value is probably explained by the fact that R12_1_2 was the sample with the lowest coverage (6.2) in the final sample set of N = 106 individuals. This suggests that for individuals with low coverage assigning clones might be more difficult and our estimates of the number of MLGs might be slightly overestimated.

| Kinship structure
The colony analysis grouped the 61 MLGs into 13 family clusters ( Figure 2; Table S2). Nine of these clusters were population-specific The DAPC analysis likewise indicated substantial overlap between populations (Figure 4). The most distinct individuals were found in populations M28 and M44 (Figure 4).

| Genetic structure and reproductive mode
The MSN showed that populations were clustered into one smaller central and four well-separated terminal groups ( Figure 5) but populations did not differ significantly from each other within the cluster.
For the most part, members of a population appeared in the same group, with a few exceptions where some individuals appeared in different branches (M25, M26, M52, M90; Figure 5).
The distribution of reproductive mode in the four family clusters identified by the colony analysis was not uniform. Four clusters contained only asexual individuals. In the remaining clusters, both sexual and asexual individuals were observed, with the proportion of sexual individuals varying from 8% to 60%. However, based on AMOVA, sexual and asexual individuals were genetically not differentiated (p = .189).
We found 18 MLGs in which more than one individual was sampled; eight of them contained only asexual polyps, while the rest contained both asexual and sexual individuals (Table S2). Of these,

| DISCUSS ION
The primary purpose of our study was to examine the genetic background and differentiation of different co-existing reproductive strategies (asexual and sexual) in Hydra oligactis natural populations.
We found that: (i) there is clear evidence for phenotypic plasticity in mode of reproduction: polyps inferred to be clones had contrasting modes of reproduction; and (ii) there was no genetic structuring with respect to the mode of reproduction: sexual polyps showed no difference from the asexual polyps. We also found high genetic relatedness between populations (second-order kinship detected between individuals belonging to distant populations) and we observed two cases of apparent sex change (males and females found within the same clonal lineages). In addition, we found no evidence for the existence of cryptic species in H. oligactis populations in the Carpathian Basin. We discuss the implications of these findings below.  In facultative clonal organisms, both sexual and asexual reproduction can occur within the same species. However, the distribution of reproductive modes (asexual vs. sexual) between individuals can be highly variable. In European populations of the planarian Schmidtea polychroa, for instance, there is clear genetic differentiation between different reproductive modes that occur in allopatry throughout most of their range, although hybridization between sexual and parthenogenetic haplotypes can be observed (Pongratz et al., 2003). In other facultatively clonal organisms, alternations between asexual and sexual phases occur predictably during the lifetime of an individual in response to specific environmental conditions (e.g., low temperatures and stressful conditions during autumn in cyclically parthenogenetic aphids and cladocerans : Decaestecker et al., 2009;Nespolo et al., 2009). H. oligactis resembles cyclical parthenogenetic cladocerans in that it switches to a sexual mode of reproduction under conditions favouring dormancy (Tessier & Caceres, 2004). Although previous studies have suggested that there may be some genetic differences between different genotypes in this species (Tökölyi, Kozma, et al., 2017;Tomczyk et al., 2015), such a difference was not observed in this study. Hydra individuals engaged in each reproductive strategy are phenotypically very different, but we did not find significant genetic differentiation between sexual and asexual polyps.
However, experiencing similar conditions in the same populations, Hydra polyps differed in their mode of reproduction. We found clear evidence for phenotypic plasticity in this variation, because in several cases genetically identical clones differed in reproductive mode (some were sexuals, others asexuals). This suggests that there are some internal factors or environmental cues other than temperature that might induce sexual reproduction in this species. First, variation in individual condition due to differences in the availability and predictability of food, population density, age, size or previous exposure to stressors and parasites are factors with a known effect on life history traits in this species (Tökölyi et al., 2016;Tökölyi, Kozma, et al., 2017;Tökölyi, Kozma, et al., 2017). Second, particular microhabitats might occur in the examined water bodies, which could provide favourable conditions for Hydra polyps with either sexual or asexual reproduction. Third, it is possible that the plasticity in the mode of reproduction within a clonal lineage is simply part of the life history of this species and could occur in all genets, regardless of any environmental impact, because the individual decision to reproduce sexually or asexually is random (i.e., a form of bet-hedging in the face of unpredictability; Simons, 2009;Slatkin, 1974).
The presence of asexual individuals in autumn suggests that clonality might have an adaptive advantage, which may be part of the reproduction strategy of this species. The advantage could stem from the fact that, contrary to hypotheses in the literature, Hydra polyps might survive the winter. This is reinforced by the fact that we also found four populations in which only asexually reproducing individuals were present, although this may be due to the random effect of our sampling. In cladocerans with a similar reproductive system it has already been observed that a strategy with a higher investment in asexual reproduction may be adaptive in a climate with relatively mild winters, where there is a low risk of freezing during the winter and selection for dormancy may be weak (Tessier & Caceres, 2004). Unfortunately, few data exist about reproductive modes in our studied species from other natural population, but field observations from two North American populations seem to support this hypothesis. In one of these populations, Bryden (1952) performed extensive observations over 3 years in Kirkpatricks Lake, Tennessee (36.3°N), without finding sexual individuals. By contrast, sexual individuals were frequently observed during the autumn at a more northerly population (Douglas Lake, Michigan; 45. 6°N;Welch & Loomis, 1924). However, it is also worth mentioning the possibility that other factors may play a role in the absence of sexually reproducing individuals in these populations, such as different food availability, parasites, different composition of Hydra microbiomes or competition with other Hydra species.
We found a high level of kinship between individuals collected from populations located tens of kilometres away (e.g., half-sibships in multiple cases; Figure 3). This suggests that geographical barriers (i.e., a lack of direct water connections separating these H. oligactis populations does not significantly obstruct the spread of hydras within such short geographical distances). Usually, such distances seem insuperable for a creature the size of a Hydra polyp (10-30 mm) without any other external dispersal facilitation mechanism (e.g., zoochory or exploiting abiotic environmental propagation factors such as wind or water). However, all our populations (except two) are located within the drainage basin of the Tisza river, and could be connected to the river during floods that occur with high regularity.
Hence, it is possible that H. oligactis polyps or resting eggs from the headwater region could travel through the Tisza river and generate a constant source of gene flow into the lowland populations studied by us. The possibility of this mode of dispersion is not confirmed with confidence by our data, because the pattern of genetic differentiation of populations does not match the pattern of geographical location of these populations. Hence, we consider that the most likely dispersal mode is that hydras (either live polyps or resting eggs) are transmitted from one water body to another by sticking to waterfowl (zoochory), as the biological disseminating role of waterfowl has been revealed in more and more species, especially plants and small invertebrates (van Leeuwen et al., 2012(van Leeuwen et al., , 2017. In addition, the genetic structure of populations also suggests that gene flow between populations is independent of their geographical location, which also supports a possibility of spread by birds, although empir-  et al., 2001;Freeland et al., 2000), possibly due to relatively low rates of sex and small propagule banks. This might be the case in Hydra as well, because in all cases described so far only a subset of individuals within the population seems to reproduce sexually and the number of eggs produced is relatively low (a few tens of eggs per female at most, Schuchert, 2010). The high incidence of clonality relative to sexual reproduction is also inferred from the population genetic statistics obtained for our populations (specifically, the high heterozygosity). Previous theoretical and practical studies (Balloux et al., 2003;Halkett et al., 2005;Meloni et al., 2013) have suggested that a high rate of clonal propagation increases the effective number of alleles and heterozygosity in a population, while an opposite effect is expected on genetic differentiation among populations and on genotypic diversity. Based on the above, it is likely that this species may follow a propagation strategy with a high rate of clonality. Individuals that reproduce asexually could reach large numbers, especially in years with mild winters, thereby reaching significantly higher numbers in the population in the spring or having an increased chance of reaching new habitats. In addition, some studies have already revealed that clonal spread appears to be a more successful reproductive strategy in a small, isolated population exposed to environmental stress (Meloni et al., 2013 (Bosch & David, 1986;Nishimiya-Fujisawa & Kobayashi, 2018;Siebert & Juliano, 2017). In H. oligactis, sex change has been observed in the laboratory (Littlefield, 1986), but is thought to be rare (Bosch & David, 1986). The samples included in this study contained five multilocus genotypes (sets of polyps inferred to be clonally derived from a common parent) with more than one polyp of known sex, and two of these multilocus genotypes showed evidence for sex change (i.e., they contained both male and female polyps). Although the presence of sequencing errors and somatic mutations implies that we cannot be fully certain that these individuals are clonally descended from a single asexual parent, they were inferred to be clones by two distinct methods: (i) the spectrum of genetic diversity method that assigns multilocus genotypes based on the distribution of genetic distances among pairs of individuals in the sample set (Rozenfeld et al., 2007) and (ii) colony, which uses a likelihood framework to separate clones from candidate relationships (first-order kinship; Jones & Wang, 2010). These observations suggest that sex change might be much more common under natural conditions: within the 61 identified clone lines, we were able to identify sex change in two cases (3.3%), but only five of these multilocus genotypes had more than one sexual polyp (which is a necessary requirement for sex change to be detected in this data set), and hence the frequency of sex-changed multilocus genotypes might be as high as 40% (two out of five), or even higher if we assume that not all sex change cases could be detected. According to the observations of Littlefield (1986), male to female transformation occurs at a low frequency when male strains are cultured at 22°C, but is not observed when they are cultured at 18°C. Hence, sex change might be induced by elevated temperatures. We have recorded water temperatures >25°C during the summer at our long-term study site near Tiszadorogma (population M28), where one of the sex change cases was inferred. Sex change in this species could be advantageous by allowing adjustment of the sex ratio to environmental or social conditions, but at present too little is known about this phenomenon in Hydra to infer more about its adaptive value. However, the possibility of sex reversal needs to be taken into account in studies of reproductive systems of natural Hydra populations and in the design of laboratory experiments.
Finally, summarizing the above, we have carried out the first population genomic study in the genus Hydra, which gives insight into patterns of genetic diversity in natural Hydra populations. We de-

DATA AVA I L A B I L I T Y S TAT E M E N T
Raw DNA sequence reads have been deposited in the U.S. National Center for Biotechnology Information (NCBI) Sequence Read Archive (Accession no. PRJNA542036). The R-script of our analysis was uploaded to Github (https://github.com/jtoko lyi/Hydra -RADSeq).