Evolution of reproductive strategies in the species-rich land snail subfamily Phaedusinae (Stylommatophora: Clausiliidae)

,


Introduction
The evolution of novel phenotypes that provide access to new ecological niches can increase the diversification rate of lineages (Yoder et al., 2010).One key innovation, which may stimulate reproductive isolation and rapid speciation, is viviparity, a reproduction strategy with extended investment per single progeny and release of live young (Zeh andZeh, 2000, Helmstetter et al., 2016).Transitions from oviparity to viviparity attracted much attention as important evolutionary changes in the phylogenies of several animal taxa, mainly vertebrates (e.g.Reynolds et al., 2002, Mank and Avise 2006, Warren et al., 2008, Pyron and Burbrink, 2014, Furness and Capellini, 2019, Horreo et al., 2019).Viviparity has originated around 150 times among vertebrates, while species with bimodal reproduction are extremely rare (Horreo et al., 2019).Yet, the selective factors that led to the evolution of viviparity are not well understood.In squamates, viviparity might have evolved as an adaptation to cold climate where viviparous females can use thermoregulatory behaviour to shorten embryonic developmental time and to reduce the exposure of embryos to stressful environmental conditions (Shine, 1983).This so-called cold-climate hypothesis, linking viviparity with lower temperatures characteristic for high latitudes and altitudes, acquired much but not unequivocal support (Hodges, 2004, Watson et al., 2014, Ma et al., 2018).This explanation is likely not applicable to most non-squamate animals, yet the herpetological perspective predominates in studies on the evolution of parity.It is thus necessary to employ research on further, phylogenetically distant animal groups, where transitions to viviparity are less known but common (e.g.Meier et al., 1999, Wong et al., 2013, Ostrovsky et al., 2016), to provide further insights into selection pressures driving switches in reproductive mode.In aquatic gastropods, for example, the evolution of viviparity may have had an effect on ecological and evolutionary success.The relatively common occurrence of this reproductive mode in freshwater snails as opposed to marine snails (Glaubrecht, 2006), can be seen as an adaptation to a harsher and less predictable habitat (Köhler et al., 2004).
While it is widely accepted that oviparity is the ancestral condition in Gastropoda and the most common mode of reproduction in land snails, the fertilized eggs of some taxa undergo a partial or complete embryonic development in the reproductive tract (Tompa, 1979ab, Heller, 2001, Ostrovsky et al., 2016).This strategy is known as embryo-retention but it corresponds to viviparous reproduction sensu Blackburn (1999) if live young are released; the latter has also been termed 'ovoviviparity' (Tompa, 1979ab, Nordsieck, 2002, Maltz and Sulikowska-Drozd, 2008).The lack of suitable oviposition sites in certain habitats, such as exposed rock-walls, was hypothesised as a factor stimulating the evolution of viviparity in land snails (Baur, 1994).The frequent occurrence of this reproductive mode in arboreal gastropods from tropical regions, e.g.Achatinellidae and Partulidae, may have the same cause (Cowie, 1992, Baur, 1994, Price et al., 2015).
Data on reproductive traits in land snails are fragmentary but suggest independent origins of embryo-retention and viviparity in approximately 30 families of Stylommatophora (Tompa, 1979a).On a finer scale, however, the evolution of parity has never been investigated in a phylogenetic framework.The species-rich stylommatophoran land snail family Clausiliidae (door snails) provides a suitable model for such research, as three modes of reproduction occur in this group: oviparity, embryo-retention (i.e.oviparity with extended embryonic development in the reproductive tract) and viviparity (e.g.Loosjes, 1953, Likharev, 1962, Maltz and Sulikowska-Drozd, 2008, Szybiak et al., 2015).Clausiliids have already been the focus of evolutionary studies on shell chirality and complex apertural barriers, including the unique structure that is the clausilium (e.g.van Moorsel et al., 2000;Uit de Weerd et al., 2004;Uit de Weerd et al., 2006;Kornilios et al., 2015).
The Phaedusinae, which comprise the three tribes Serrulinini, Phaedusini and Synprosphymini (Bouchet et al., 2017), and which include around 600 known species (Nordsieck, 2007;Hunyadi andSzekeres, 2016, Páll-Gergely andSzekeres, 2017), are the most speciose subfamily of clausiliids (Nordsieck, 2007).The range of this subfamily consists of two isolated parts: (1) western Eurasia (Pontic-Caucasian-Hyrcanian region) harbouring Serrulinini, and (2) East and Southeast Asia harbouring Phaedusini and Synprosphymini.Molecular studies have so far yielded only limited insight into the phylogenetic relationships within Phaedusinae and respective diversification ages.Uit de Weerd and Gittenberger (2013) showed that the East and Southeast Asian tribe Phaedusini constitutes a well-supported clade nested within the western Eurasian Serrulinini.According to their study, the lineage leading to the Phaedusini arose 22.7 Ma in western Eurasia and spread eastwards prior to the aridification of the interlaying region.Pre-Pleistocene fossils of Phaedusinae are known from western Eurasia (Europe north of the Tethys) only (Nordsieck, 2015).
The diversity of reproduction modes in Phaedusinae was reported by Loosjes (1953), Nishi (1980), Minato (1983Minato ( , 1994)), Nordsieck (2001, 2002), and Sulikowska-Drozd et al. (2020).Nordsieck (1998Nordsieck ( , 2001Nordsieck ( , 2002) ) regarded certain shell characters to be correlated with reproductive modes and distinguished in his taxonomical system at that time the oviparous Megalophaedusini and the 'ovoviviparous' Phaedusini tribes; the latter group was characterised by a broad clausilium plate.In his subsequent publications, this taxonomical division was not maintained (Nordsieck, 2007(Nordsieck, , 2011)).Meanwhile, evidence was provided that a wide patency of apertural barriers is a common adaptation to the delivery of live young in clausiliids of different genera (Sulikowska-Drozd et al., 2014).This is congruent with growing awareness of common homoplasy in these complex structures (van Moorsel et al., 2000, Uit de Weerd et al., 2004).Recent findings indicated a link between phylogeny and reproduction mode for Phaedusinae from Japan (Motochin et al., 2017).Among six main clades, two included only oviparous and two only 'ovoviviparous' taxa; in the remaining two clades both strategies occurred (Motochin et al., 2017; note that embryo-retaining species were not distinguished and possibly regarded as oviparous).
Given this background, the present study aims to reconstruct the evolution of reproductive strategies in the door snail subfamily Phaedusinae based on a combination of time-calibrated molecular phylogenetics, reproduction mode examinations, and ancestral state reconstruction.

Material studied
Material was collected from various parts of Eurasia (Fig. 1a).Our dataset includes representatives of all three Phaedusinae tribes: Serrulinini, Phaedusini and Synprosphymini.Locality and specimen details including specimen images are available in the Barcode of Life Data Systems (BOLD) dataset DS-PHAEDUSA [dx.doi.org/https://doi.org//10.5883/DS-PHAEDUSA].
Individual clausiliids used for DNA sequencing were either collected in the wild (F0 generation) or from the F1 generation of laboratory colonies (see DS-PHAEDUSA).For 12 Phaedusa angustocostata and Phaedusa ramelauensis specimens, DNA isolates were provided by F. Köhler (see Köhler and Mayer, 2016 for extraction methodology) and for 11 of them, COI sequences from Köhler and Mayer (2016) were used (see DS-PHAEDUSA).
The present study follows the suprageneric classification of Bouchet et al. (2017).At genus-level and below, we follow Motochin et al. (2017) for all taxa revised by them, and Nordsieck (2007) for the remaining taxa, with the single exception of the genus Renschiphaedusa, which is regarded as a synonym of Phaedusa following Köhler and Mayer (2016).Subgenus and subspecies names are only used for selected taxa.

DNA processing and analysis
DNA was isolated from a total of 141 ethanol-preserved individuals.We amplified seven markers: two mitochondrial markers, the barcoding fragment of the cytochrome c oxidase subunit I gene (COI, ca.650 bp) and a fragment of the 16S ribosomal RNA gene (ca.390 bp), and five nuclear markers, fragments of the 28S ribosomal RNA (ca.950 bp), histone H3 (326 bp) and histone H4 (266 bp) genes, and the complete internal transcribed spacers 1 (ITS1, ca.650 bp) and 2 (ITS2, ca.300 bp); see Appendix, Tables S1 and S2 for details.
Obtained sequence data were tested for contamination via BLASTN (Altschul et al., 1990) and verified as Clausiliidae.Heterozygous sites were only observed for ITS1 and ITS2; in such cases, only one version was used.GENEIOUS 6.1.4(https://www.geneious.com;Biomatters 2013) was used for sequence assembly, editing and quality trimming (the latter included the removal of the partial H3-H4 spacer region from the H3 and the H4 sequences).All generated and processed sequences were deposited in GenBank and are available in DS-PHAEDUSA including metadata.Altogether, 682 new DNA sequences were generated in this study (125 sequences of COI, 126 of 16S, 98 of 28S, 71 of ITS1, 52 of ITS2, 103 of H3 and 107 of H4).
Alignments of sequence data, including those for subsequent analyses, were performed using MAFFT with automatic settings selecting the best algorithm for each dataset (Katoh et al., 2005; all alignments are deposited in the Dryad repository).The full datasets of all single markers (see DS-PHAEDUSA) were tested for substitution saturation with DAMBE 5.3 (Xia, 2013) using the index proposed by Xia et al. (2003).The test revealed little saturation under the assumption of a symmetrical tree for all markers and for most markers also under the assumption of an asymmetrical tree (Table S3).Molecular distances were presented as trees using the Neighbor-Joining method (Saitou and Nei, 1987) with the Kimura 2-parameter model (Kimura, 1980) and 1000 bootstrap replicates (Felsenstein, 1985) in MEGA X (Kumar et al., 2018) for the full datasets of each single-marker (Fig. S1).In order to assess potential discrepancies in phylogenetic signal between different genes, individual phylogenies were reconstructed using the full datasets of each single-marker in a Maximum Likelihood (ML) approach through RAxML 8.2.8 (Stamatakis, 2014).The best-scoring ML trees were produced using the GTR + I + G substitution model.Bipartition information was drawn from the phylogenies obtained with the rapid hillclimbing tree search algorithm.Statistical supports were estimated with thorough bootstrap tests set to 1000 repetitions.These analyses did not reveal substantially conflicting patterns between individual markers, while they provided a generally poor resolution for deeper divergences (Fig. S2).

Reconstruction of time-calibrated phylogeny
For the reconstruction of a time-calibrated phylogeny, the full dataset of sequences (153 individuals; see DS-PHAEDUSA) was reduced.This reduced dataset (58 individuals; see individuals labelled in black in Fig. S3) comprised a single individual per taxon except for Megalophaedusa pinguis platyauchen, Phaedusa cf.lypra, Phaedusa stenothyra, Leptacme sykesi and Megalophaedusa (cf.) martensi, which showed relatively high levels of intraspecific variation, and Phaedusa ramelauensis, where euphallic and aphallic specimens were studied (see DS-PHAEDUSA).For these taxa, 2-3 individuals were included.
In order to apply the molecular clock calibration of Uit de Weerd and Gittenberger (2013) to our data, we combined the reduced dataset of 28S, H3 and H4 markers, with the dataset from Uit de Weerd and Gittenberger ( 2013) and carried out an initial phylogenetic analysis following Uit de Weerd and Gittenberger (2013) (see Appendix, Tables S4, S5 and Fig. S3).In addition, a ML analysis was performed in RAxML 8.2.8 using the same settings as for the single-markers phylogeny reconstructions (Fig. S5a).Based on the results from this initial analysis, a set of three secondary calibration points was selected for following main analysis (Table S4, Fig. S3).
For the main phylogenetic analysis, we combined the reduced dataset with five individuals from the subfamily Clausiliinae from Uit de Weerd and Gittenberger (2013) (see Appendix, Fig. S3).The phylogeny was reconstructed using Bayesian Inference in BEAST 2.5 (Bouckaert et al., 2019), using all 7 markers produced in the study, for details see Appendix.All single markers were treated as separate partitions following the partitioning scheme selected by PartitionFinder 2.1.1.(Lanfear et al., 2012), evolutionary models were selected using bModeltest (Bouckaert and Drummond, 2017; results in Table S5).The Effective Sampling Size of each parameter was verified to be above 200 in Tracer 1.7 (Rambaut et al., 2018).In order to assess bootstrap support values (BSVs), we used the ML approach in RAxML 8.2.8 with the same settings as for the single-markers phylogeny reconstructions.

Identification of reproductive strategies
Living individuals of 33 Phaedusinae taxa (see taxa with data source 'Laboratory culture' in Table S6) were held in laboratory culture at the Dept.Invertebrate Zoology and Hydrobiology, University of Lodz, Poland, in order to identify their reproductive strategy.Observations were conducted on individuals kept in climatic chambers (POL-EKO KK350 TOP) under controlled temperature (15-21 • C) and natural light conditions.Individuals from different populations were bred separately, summing up to 42 laboratory colonies.Snails were kept in pairs or groups of three in plastic containers (300 cm 3 ) with cellulose tissue and fragments of limestone as source of calcium, and were fed with lettuce, cucumber and zucchini.Humidity in the containers was almost constant and close to 100% due to regular spraying.During the first months of the study, the containers were inspected every 2-3 days for the presence of eggs and neonates, which were immediately transferred to separate containers with damp cellulose tissue.Freshly deposited eggs were checked every three days for hatching.Later, the reproduction was monitored weekly and progeny was preserved in ethanol.Additionally, mature individuals (3-10 per taxon) were preserved in ethanol and dissected to determine whether they retained developing embryos in the reproductive tract.Each taxon was observed for at least one year.Based on these long term observations and dissections, taxa were classified into three categories according to their reproductive mode (see Sulikowska-Drozd, 2009): viviparous taxa (litters of neonates observed regularly in culture; large shelled embryos found in dissected animals), oviparous taxa (clutches of eggs observed regularly in culture; egg incubation lasts 2 weeks or longer; shelled embryos never found in dissected animals) and embryo-retaining taxa (clutches of eggs observed regularly in culture; egg incubation lasts shorter than 2 weeks; eggs with shelled embryos found in dissected animals, see Fig. 1b for examples).
We furthermore attempted to classify taxa according to their reproductive mode by means of X-rays of individuals (see taxa with data source 'X-ray studies' in Table S6).The X-ray studies were carried out in the μ-CT laboratories of the Institute of Palaeontology, Polish Academy of Science, Warsaw, Poland, and of the University of Silesia, Chorzów, Poland.This method allowed to rule out oviparous reproduction whenever X-rayed individuals were gravid with shelled embryos (see Sulikowska-Drozd et al., 2020).For distinguishing embryo-retaining and viviparous taxa, embryo size and size distribution were taken into account; embryonic shells of more than 2 whorls were regarded as evidence for viviparity.All specimens that were used for X-ray study were also utilized for DNA sequencing.For Phaedusa angustocostata and P. ramelauensis, viviparity was assumed based on Köhler and Mayer (2016).

Reconstruction of ancestral reproduction strategies
While changes in reproductive strategy took place in various animal groups over evolutionary times, such transitions are presumably relatively rare events.Based on this assumption, a maximum parsimony approach was chosen to reconstruct the evolution of reproductive strategies among Serrulinini and Phaedusini.Synprosphymini were not included in this analysis given that no data on their reproduction were available.The analysis was run in the program Mesquite 3.6 (Maddison and Maddison, 2018).'Oviparous', 'embryo-retaining' and 'viviparous' were coded as ordered states in a multistate analysis, with embryoretaining treated as an intermediate state between oviparous and viviparous.The character states were mapped on a subset of 1000 trees from the main analysis using multiple maximum parsimony reconstructions (MPR option) per tree.The results were summarized on the consensus tree from the main analysis.In addition, we calculated the number of transitions between the alternative reproductive strategies for each of the 1000 trees, using the Pscores command in PAUP* 4.0 (Swofford, 2001).

Time calibrated reconstruction of phylogeny
Our phylogenetic analysis (Fig. 1c; see Fig. S4 for HPD and Fig. S5b for ML results) supports the monophyly of a group formed by the Serrulinini and the Phaedusini (clade SE + PH; BPP = 1.0,BSV = 96).The phylogenetic relationships of the Synprosphymini remain less clear; our analysis positions them as sister group of the Garnieriinae, however, without good support (BPP = 0.92, BSV = 57).
Within clade SE + PH, the studied representatives of the Serrulinini form two successively basal lineages.The more basal of them consists of Pontophaedusa funiculum and diverged from its sister group in the Late Eocene, 37.3 Ma (95% HPD: 48.9-25.8).The less basal lineage, which is composed of Caspiophaedusa perlucens and Pravispira semilamellata, diverged from its sister group in the Oligocene, 33.2 Ma (95% HPD: 43.6-23.3).The remaining lineages form the well-supported clade PH (BPP = 1.0,BSV = 93, 95% HPD: 36.8-20.0),which corresponds to the tribe Phaedusini.This group diversified into eight major lineages (PH1-PH8) in the Oligocene, with largely unresolved interrelationships.Five of those lineages (PH1, PH3, PH4, PH7 and PH8) lead to supraspecific clades, each well-supported in our analysis (BPP ≥ 0.96, BSV 66-100).Clade PH1, which is composed of species from Vietnam, China and East Timor, diversified mainly through the Miocene.In this clade, Euphaedusa porphyrea is the sister group of all other taxa, which form two subclades.Clade PH2 only consists of Acanthophaedusa ookuboi from China.Lineages PH3-PH6 are composed of species living in Vietnam.The most species-rich group of these in our analysis is PH4 with a diversification spanning through the Miocene.Lineage PH3 contains Hemiphaedusa polydona and Hemiphaedusa billeti, which split in the Miocene.Lineage PH5 only consists of P. stenothyra and lineage PH6 is represented by the single species Oospira antibouddah.Lineage PH7, which diversified during the Miocene, includes only species found in Taiwan.Clade PH8, which began to diversify 25.1 Ma (95% HPD: 17.41-32.46),comprises numerous species from Japan, Vietnam, and Taiwan and consists of four lineages (J1-J4) with poorly resolved interrelationships.Clade J1, which includes several species from all three areas, diversified in the Miocene according to our analysis.Within the clade, Tauphaedusa sheridani from Taiwan and Tauphaedusa tau from Japan form the sister group of all other taxa.Lineage J2 consists of species of the genus Reinia from Japan and Taiwan, which diverged in the Pleistocene, and lineage J3 consists only of Changphaedusa horikawai from Taiwan.Clade J4, which includes only species from Japan, diversified in the Miocene.This clade consists of two subclades that correspond to the genera Stereophaedusa and Megalophaedusa.

Reproduction strategies in studied Phaedusinae and reconstruction of ancestral reproduction strategies
Among the studied Phaedusinae taxa, we identified 16 as oviparous, 18 as viviparous and three as embryo-retaining; we left 12 taxa without conclusively determined reproductive strategy (Fig. 1a, c, Table S6).
Focussing on the well-supported clade of Phaedusini and Serrulinini (clade SE + PH), none of the reproduction modes was uniquely synapomorphic for a lineage within this clade.According to our ancestral state reconstruction (Fig. 1c, Table S7), clade SE + PH evolved from an oviparous most recent common ancestor.This state was found as most parsimonious for the respective node in each of the 1000 trees.Our data are less clear about the ancestral reproductive strategy of the Phaedusini (clade PH).The most often found ancestral state for this group across the 1000 trees examined was oviparous (50%), followed by viviparous (39%) and embryo-retaining (11%).Of the major subclades distinguished within this clade, clade PH1 evolved from a viviparous ancestor, according to our analyses, a condition that was reversed in the oviparous Oospira vanbuensis.Clade PH3 consists entirely of taxa without conclusively determined reproductive strategy.The ancestral states of clades PH4, PH7, and PH8 differ among trees.For clade PH8, the result of our analysis was 40% oviparous, 14% embryo-retaining and 43% viviparous (3% node absent).The total number of transitions in reproductive strategies within clade SE + PH was between 13 and 16.

Reproductive modes of recent Phaedusinae and inference of ancestral strategies
We showed that similar reproduction modes evolved independently among the studied Phaedusinae.Both oviparity and viviparity occur frequently in not closely related lineages of the East and Southeast Asian Phaedusini.According to our analysis multiple switches must have taken place during the group's diversification, which indicates a high evolutionary plasticity.
If viviparity would have evolved as an adaptation to cold climates, as it was suggested for squamates in the cold-climate hypothesis (Shine, 1983), one may expect viviparous snail species to be more common in colder regions.The majority of viviparous stylommatophoran species, however, occurs in the tropics and subtropics (Tompa, 1979b) and thermoregulatory behaviour has never been recorded in land snails.
Acknowledging that a test for potential correlations between temperature and reproduction mode is beyond the scope of the present paper, our study does not reveal an obvious latitudinal pattern in the distribution of different strategies among Phaedusini.This may suggest that no strategy is clearly superior in the macroclimatic context.However, oviparous and viviparous species may differ in their distribution at microhabitat level.Motochin et al. (2017) noticed a correlation between viviparous reproduction and an arboreal lifestyle among Phaedusini from Japan.As none of the oviparous species studied by Motochin et al. (2017) was regarded as a tree-dweller, it may be assumed that the snails require soil for egg-laying.This would be in line with the hypothesis of Baur (1994), that viviparity can allow reproducing successfully in habitats unavailable for oviparous land snail species.Outside Japan, however, this pattern is not universal given that the oviparous species Oospira vanbuensis is arboreal (personal observation).Besides, we assume that so-called arboreal clausiliids retreat to the soil during the dry or winter season, both in Japan and in Indochina.
While viviparous animal species may produce less progeny than their oviparous relatives, the survival rate of their offspring is presumably higher (see Wourms, 1981 for fishes).For clausiliids, however, such an evolutionary trade-off has not yet been found (Sulikowska-Drozd, 2009, Szybiak et al., 2015, Sulikowska-Drozd et al., 2018a,b); for species differing in reproduction mode, the available evidence does not suggest considerable fecundity differences (neither observed in laboratory conditions nor in population dynamics in the field).
Oviparity is likely the ancestral mode of reproduction in stylommatophoran land snails (Tompa, 1979ab, Heller, 2001) and it predominates in clausiliids (Maltz and Sulikowska-Drozd, 2008).All studied Serrulinini species are exclusively oviparous and, to our knowledge, no viviparous species have been identified in this tribe.Two oviparous lineages are basal to the Phaedusini in our phylogenetic analysis (Pontophaedusa and Caspiophaedusa/Pravispira) and, based on our reconstruction, this reproductive strategy is ancestral for the whole clade SE + PH.It should be noted, however, that the Serrulinini exhibit a significant variation in the structure of egg envelopes, which is unique among oviparous clausiliids.In contrast to the typical, partly calcified clausiliid egg, where discrete crystals of CaCO 3 are dispersed randomly in the jelly layer (Tompa, 1976, Maltz andSulikowska-Drozd, 2008), Pontophaedusa funiculum produces a hard egg envelope (Páll-Gergely and Németh, 2008).The spiral deposition of calcium crystals in the egg membranes of Caspiophaedusa and Pravispira also stands out from the usual pattern (Stolarski et al., 2019).Further insights in the evolution of reproduction strategies among the Phaedusinae might be gained from the Synprosphymini, whose reproduction biology remains still unknown.Future studies should thus focus particularly on this group.
The reproductive mode of the most recent common Phaedusini ancestor remains ambiguous in our analysis (only 50% of the trees support oviparity) due to the insufficient resolution of the relationships between lineages PH1-PH8.The proportion between species of different reproductive modes in the whole clade is yet to be surveyed.
For the Phaedusini subclade PH8, which contains species from the Japanese Archipelago and Taiwan, our dataset partly overlaps with the molecular phylogeny of Motochin et al. (2017).In the subset of Japanese clausiliids included in our analysis, oviparous species predominate, yet, all sublineages (J1-J4) include viviparous taxa.Our results confirm the relatively close relationship of the exclusively viviparous genus Tauphaedusa and the reproductively diverse genus Zaptyx.According to Motochin et al. (2017), the genus Zaptyx includes oviparous and 'ovoviviparous' taxa.The former reproductive mode predominates, which could suggest its ancestral origin.The distribution of viviparity among different subgenera of Zaptyx (Motochin et al., 2017) may indicate a high evolutionary plasticity of parity-modes in this group.During laboratory observations, we found the intermediate strategy of embryoretention for two Zaptyx species (Z.hachijoensis and Z. kikaiensis), which were classified as 'ovoviviparous ' and oviparous, respectively in Motochin et al. (2017).This is in line with our hypothesis on the significant role of embryo-retention in the repeated evolution of viviparity in Phaedusini (see below).Our ancestral state reconstruction showed that viviparity or embryo retention could be ancestral for the entire J1 lineage (59% and 39%, respectively).In contrast, a predominance of oviparity was found for the ancestor of lineage J4 (71%).From the two groups within this lineage, Megalophaedusa and Stereophaedusa, the former includes only oviparous species, while the latter includes a T. Mamos et al. similar number of 'ovoviviparous' and oviparous species according to Motochin et al. (2017).Based on our analyses, which include only four Stereophaedusa species, oviparity is the most probable ancestral mode in this group and at least two switches between reproductive strategies have to be assumed.
A remarkable finding of our study concerns clade PH1, which occurs in Southeast Asia.All except one of the species in this clade are viviparous, and our analysis assumes a viviparous most recent common ancestor of the clade in 100% of the analysed trees.The oviparous Oospira vanbuensis is deeply nested in clade PH1 with high support values, which indicates a reversal from viviparity to oviparity during the evolution of this species.To our knowledge, this is the first evidence for such a reversal of parity in land snails.
A reversal to a previous reproductive strategy seems to be very rare in animals.In vertebrates, viviparity has evolved repeatedly but it is highly disputed if oviparity can re-evolve (Blackburn, 2006, 2015, Pyron and Burbrink, 2014, King and Lee, 2015, Recknagel et al., 2018, Furness and Capellini, 2019).According to King and Lee (2015), a reversal may be possible only for a finite time window after the evolution of viviparity, as an accumulation of mutations irreversibly degrades the genetic pathways related to oviparity.In squamates, for example, the evolution of viviparity requires the loss of the eggshell, which protects the embryo from desiccation but restricts gas exchange if the embryo is kept in the reproductive tract (Griffith et al., 2015).In consequence, during reversal towards oviparity, the developing embryo may not be drought resistant enough to survive in ambient environments (Griffith et al., 2015, Shine, 2015).
In contrast, the egg envelope of land snails, which provides mechanical support and calcium necessary for the embryonic shell formation, does not hinder gas exchange and because of this, does not protect from desiccation (Tompa, 1979ab).As a consequence, the egg envelope may keep its basic structure during the shift towards viviparity.The embryos of some viviparous clausiliids (e.g.Alinda biplicata, Idyla bicristata and Stereophaedusa jacobiana) are initially enveloped in an egg membrane with calcium carbonate crystals, which resembles the structure of the partly calcified egg of oviparous taxa (Maltz and Sulikowska-Drozd, 2012, Sulikowska-Drozd et al., 2018a,b, 2019).The presence of these crystals suggests that the ability of egg shell production could be maintained in these species.In individuals of Reinia variegata, Tauphaedusa tau and Changphaedusa horikawai that were dissected in the present study, however, no calcium carbonate crystals were present in the egg membrane.While these data are still preliminary, such histological differences between viviparous taxa from different clades may be related to the multiple independent origins of their reproduction strategy.Possible differences in egg membrane structure and in the sources of embryo nutrition among viviparous taxa from different phylogenetic lineages, should be investigated in future studies.Furthermore, future research should focus on potential differences between the reproduction of Oospira vanbuensis and other oviparous Phaedusini in order to examine whether the oviparity of this species really results from a reversal of parity.
The frequent reproduction mode switches during Phaedusini evolution could have been promoted by a high fitness of species with an intermediate strategy.Embryo-retention, was recorded in two major Phaedusini lineages (clades PH7 and PH8).Selection towards it was likely involved in the repeated shifts of reproductive modes.Embryoretention could represent a stepping stone, from which selection can act in either direction, maintaining evolutionary plasticity.In embryoretaining species, the genetic heritage of their oviparous ancestors may still be preserved, while adaptations to a delayed delivery of embryos have evolved.These adaptations may include serous subepithelial cells in the free oviduct, which secrete products that provide oncotic pressure and create a suitable environment for embryos (Maltz et al., 2017), and a widening of the shell channel, which allows passing of shelled embryos (Sulikowska-Drozd et al., 2014).
The evolutionary plasticity regarding their reproductive strategy could have provided the Phaedusini access to new ecological niches and thus boosted their diversification into all main lineages in a relatively short time window.The exact factors, which promote shifts between strategies, however, are yet to be discovered.

Phaedusinae systematics and biogeography
The present study challenges the division of the Phaedusinae into the tribes Serrulinini, Phaedusini, and Synprosphymini (Bouchet et al., 2017).In our molecular phylogeny, the Serrulinini species Pontophaedusa funiculum forms the sister group to the Phaedusini and all other Serrulinini.Our analysis thus confirms the paraphyly of Serrulinini, suggested already by Nordsieck (2007) and shown by Uit de Weerd and Gittenberger (2013), and endorses the classification of Pontophaedusa in a separate tribe as suggested by Bouchet et al. (2017).Our phylogenetic reconstruction also showed that the Synprosphymini represent a phylogenetically very distinct lineage, which may be more closely related to the Garnieriinae than to other Phaedusinae.Given that the respective node was not well-supported and that the knowledge of deeper relationships in the family Clausiliidae remains limited (see Uit de Weerd and Gittenberger, 2013), however, future studies expanding the set of molecular markers are needed to confirm this result.Many Phaedusini genera that inhabit Southeast Asia (e.g.Oospira, Phaedusa) do not represent monophyletic groups in our analysis and require a taxonomic revision.The changes introduced to the Phaedusinae system by Motochin et al. (2017) are supported by our study.
Almost 200 Phaedusinae species and subspecies occur on the Japanese Islands (Minato, 1994, Motochin et al., 2017), of which 16 were included in our analysis.Our findings indicate that the most recent common ancestor of the clade that includes all Japanese taxa (clade PH8) occurred at the end of the Oligocene, ca. 25 Ma, while intense diversification of Japanese lineages started about 15 Ma, which correlates to the full isolation of the archipelago due to the opening of the Okhotsk Sea (ca.20 Ma) and the Sea of Japan (ca. 15 Ma) (see Tojo et al., 2017 and references therein).The Taiwanese species Tauphaedusa sheridani clusters among the Japanese taxa, which suggests dispersal from Japan to Taiwan.Potential migration routes between Japan and continental Asia via Taiwan developed about 5 Ma, and again during Pliocene and Pleistocene glaciation peaks (Shih et al., 2006).The relatively recent split between the Japanese and Taiwanese members of the genus Reinia (ca. 2 Ma) suggests that dispersal along an ice-age land bridge may have been involved.However, migration might have also occurred over sea probably via ocean currents, as it was suggested by Motochin et al. (2017) for clausiliid species from isolated islands.The same mechanism was also proposed for other terrestrial animal groups (Kobayashi et al., 2014, He et al., 2015).Notably in our phylogeny, Zaptyx hachijoensis, which is endemic to the oceanic Izu Islands southeast of Japan's main island Honshu, is relatively closely related to Messageriella gargomyi from southern Vietnam.This relationship suggests dispersal in the opposite direction to the major Kuroshio ocean current, probably via migrating birds, a mechanism that has already been suggested for other clausiliids (Gittenberger et al., 2006).

Conclusion
The present study represents the first reconstruction of the evolutionary history of parity in land snails with advanced analytical tools.Oviparity is likely the ancestral reproduction strategy of the clade consisting of Phaedusini and Serrulinini (clade SE + PH).The diversification of Phaedusini is associated with many switches in reproduction mode.The evolutionary success of this group might thus not only depend on the evolution of viviparity but also on the maintenance of various strategies.

Fig. 1 .
Fig. 1.Sampling localities and phylogeny of studied taxa.(a) Sampling localities of studied Phaedusinae specimens with information on reproductive strategies (legend in figure); coloured ellipses indicate geographical areas.(b) Examples of Phaedusinae species with different reproduction strategies including respective eggs or neonates; from the left: Oospira vanbuensis, O. formosensis, and Tauphaedusa sheridani.(c) Time-calibrated Bayesian maximum clade credibility tree of Phaedusinae, Garnieriinae and outgroup taxa (not shown) based on sequence data from 7 markers (main analysis).Node ages were inferred from common ancestral heights.Bayesian posterior probabilities (only values ≥ 0.50 shown) and bootstrap values (only values ≥ 50 shown) from Maximum Likelihood framework (see Fig. S5b) are provided at the respective nodes.For outgroup taxa, highest posterior density, and sample IDs see Fig. S4.Circles at branch tips indicate present reproduction modes (no circle: unknown), pie charts at nodes indicate ancestral reproduction strategies based on the maximum parsimony reconstruction (black: oviparous, grey: embryo-retaining, white: viviparous).Taxon names are highlighted in colour according to the geographical area of sampling localities, as shown in (a).