# A rapid rate of sex-chromosome turnover and non-random transitions in true frogs

Daniel L Jeffries, Guillaume Lavanchy, Roberto Sermier, Michael J Sredl, Ikuo Miura, Amaël Borzée, Lisa N Barrow, Daniele Canestrelli, Pierre-André Crochet, Christophe Dufresnes, et al.

► To cite this version:

Daniel L Jeffries, Guillaume Lavanchy, Roberto Sermier, Michael J Sredl, Ikuo Miura, et al. A rapid rate of sex-chromosome turnover and non-random transitions in true frogs. Nature Communications, 2018, 9 (1), 10.1038/s41467-018-06517-2. [DOI: 10.1038/s41467-018-06517-2](https://doi.org/10.1038/s41467-018-06517-2)

HAL Id: hal-04886557

https://hal.science/hal-04886557v1

Submitted on 14 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

ARTICLE

DOI: 10.1038/s41467-018-06517-2

OPEN

# A rapid rate of sex-chromosome turnover and non-random transitions in true frogs

Daniel L. Jeffries ID <fnref n="1" />, Guillaume Lavanchy<fnref n="1" />, Roberto Sermier<fnref n="1" />, Michael J. Sredl², Ikuo Miura³, Amaël Borzée ID ⁴, Lisa N. Barrow⁵, Daniele Canestrelli ID ⁶, Pierre-André Crochet⁷, Christophe Dufresnes<fnref n="1" />, Jinzhong Fu⁸, Wen-Juan Ma<fnref n="1" />, Constantino Macías Garcia ID ⁹, Karim Ghali<fnref n="1" />, Alfredo G. Nicieza¹⁰,¹¹, Ryan P. O'Donnell², Nicolas Rodrigues ID <fnref n="1" />, Antonio Romano¹²,¹³, Íñigo Martínez-Solano¹⁴, Ilona Stepanyan<fnref n="15" />, Silvia Zumbach¹⁶, Alan Brelsford¹⁷ & Nicolas Perrin ID <fnref n="1" />

The canonical model of sex-chromosome evolution predicts that, as recombination is suppressed along sex chromosomes, gametologs will progressively differentiate, eventually becoming heteromorphic. However, there are numerous examples of homomorphic sex chromosomes across the tree of life. This homomorphy has been suggested to result from frequent sex-chromosome turnovers, yet we know little about which forces drive them. Here, we describe an extremely fast rate of turnover among 28 species of Ranidae. Transitions are not random, but converge on several chromosomes, potentially due to genes they harbour. Transitions also preserve the ancestral pattern of male heterogamety, in line with the 'hot-potato' model of sex-chromosome transitions, suggesting a key role for mutation-load accumulation in non-recombining genomic regions. The importance of mutation-load selection in frogs might result from the extreme heterochiasmy they exhibit, making frog sex chromosomes differentiate immediately from emergence and across their entire length.

<footnote n="1">¹Department of Ecology and Evolution, University of Lausanne, CH-1015, Lausanne, Switzerland. ²Arizona Game and Fish Department, 5000 W. Carefree Highway, Phoenix, AZ 85086, USA. ³Amphibian Research Center, Hiroshima University, Higashi-Hiroshima 739-8526, Japan. ⁴Division of EcoScience and Department of Life Sciences, Ewha Womans University, Seoul 03760, Republic of Korea. ⁵Museum of Southwestern Biology, MSC03 2020, 1 University of New Mexico, Albuquerque, NM 87131, USA. ⁶Department of Ecological and Biological Science, University of Tuscia, O1100, Viterbo, Italy. ⁷CEFE, CNRS, University Montpellier, University Paul Valéry Montpellier 3, EPHE. IRD, Route de Mende, Montpellier 34293, France. ⁸Department of Integrative Biology, University of Guelph, Guelph, ON N1G 2W1, Canada. ⁹Instituto de Ecología, Universidad Nacional Autónoma de México, Ciudad Universitaria, Mexico City, 04500, Mexico City, Mexico. ¹⁰Research Unit of Biodiversity, UO-CSIC-PA, 33006, Mieres, Spain. ¹¹Department of Biology of Organisms and Systems, Universidad de Oviedo, 33600, Mieres, Spain. ¹²Consiglio Nazionale delle Ricerche, Istituto per i sistemi agricoli e forestali del mediterraneo, Via Patacca 84, I-80056, Ercolano, NA, Italy. ¹³MUSE - Museo delle Scienze, Sezione di Zoologia dei Vertebrati, corso del Lavoro e della Scienza 3, 38122, Trento, Italy. ¹⁴Department of Biodiversity and Evolutionary Biology, Museo Nacional de Ciencias Naturales, CSIC, José Gutiérrez Abascal, 2, Madrid 28006, Spain.</footnote>

<footnote n="15">¹⁵Scientific Center of Zoology and Hydroecology, National Academy of Science, Republic of Armenia, P. Sevak str. 7, 0014, Yerevan, Armenia. ¹⁶Info Fauna - karch, UniMail, Bellevaux 51, 2000, Neuchâtel, Switzerland. ¹⁷Evolution, Ecology, and Organismal Biology Department, University of California Riverside, Riverside, CA 92521, USA. These authors contributed equally: Daniel L. Jeffries, Guillaume Lavanchy. Correspondence and requests for materials should be addressed to D.L.J. (email: dljeffries86@gmail.com)</footnote>
Sex chromosomes experience very different evolutionary forces compared to the rest of the genome, owing to their differential occurrence in males and females. A male-beneficial mutation occurring close to a newly evolved male-determining locus, for instance, is expected to spread, even if highly detrimental to females, because linkage with the sex-locus makes it more likely to be transmitted to sons than daughters. The canonical model of sex-chromosome evolution predicts that a region of suppressed recombination will then expand outwards from the sex-determining locus, possibly via inversions, thus increasing linkage with sexually antagonistic genes¹,². Over time, however, this loss of recombination will lead to the accumulation of potentially deleterious mutations, repetitive elements and gene loss in sex-limited chromosomes (Y or W) due to Hill-Robertson interactions³ and Muller's Ratchet⁴. This progressive build-up of differences between gametologs (i.e. between X and Y, or between Z and W) is thought to have given rise to the highly heteromorphic sex chromosomes prevalent in many clades, including mammals, birds, flies and some snakes⁵⁻⁷.

However, more and more examples of homomorphic sex chromosomes are now being discovered, notably in the majority of poikilothermic vertebrates, i.e. in most fish, lizards and amphibians⁸,⁹. This homomorphy implies that gametologs have recombined relatively recently in these taxa, and one main hypothesis for how this can occur is through recurrent sex-chromosome turnovers, i.e. a switch in the chromosome pair used for sex determination¹⁰,¹¹. In systems where a turnover has recently occurred, sex chromosomes will be relatively young; thus gametologs will not have had time to substantially differentiate from one another. Indeed, frequent turnovers have been documented in several lineages of fish¹²,¹³, amphibians¹⁴⁻¹⁷ and reptiles¹⁸, and are thought to have contributed to the homomorphy of sex chromosomes in these species. However, it should be noted that turnover rates have sometimes been inferred from changes in the patterns of heterogamety¹⁸,¹⁹, which can grossly underestimate true values. Indeed, several groups seemingly show constant heterogamety despite high turnover rates¹³,¹⁴,²⁰,²¹.

Despite the importance of turnovers for sex-chromosome evolution, one important question remains unanswered: why do some taxa fall into the evolutionary trap of extensive sex-chromosome differentiation and degradation, while others can avoid this by continuously swapping them? Answering this question requires knowledge of the evolutionary forces driving turnovers. A turnover can occur if the existing sex-locus is translocated to an autosome²²,²³, or if a new gene acquires the sex-determining role via mutation²⁴,²⁵. The subsequent fixation of the new sex chromosome is thought to be mediated by one (or a combination) of four main evolutionary forces²⁶, namely (i) genetic drift²⁷⁻²⁹; (ii) sex-ratio selection, induced by sex biases arising, e.g., from meiotic-drive elements or endoparasites³⁰⁻³²; (iii) sexually antagonistic selection on a gene linked to the new sex determinant²⁴,³³; and (iv) deleterious mutation load accumulating on the non-recombining Y or W chromosome³⁴,³⁵.

The predicted patterns of turnovers under each of these forces are distinct from each other. Turnovers induced by mutation-load selection are predicted to strictly preserve the heterogametic sex, because a switch (e.g. from an XY to a ZW system) would require fixation of the heavily loaded sex-limited chromosome (Y or W) as an autosome³⁵. Under drift-induced turnovers, patterns of heterogamety are generally expected to be maintained 2-4 times more often, as X or Z chromosomes are more likely to be fixed than Y or W chromosomes due to their higher initial frequency (75% vs. 25%). This ratio, however, decreases along with effective population size, and can even reverse in case of extreme polygyny, where XY becomes more likely to transition towards ZW than towards XY²⁹.

We also expect to see changes in the heterogametic sex via turnovers induced by either sex-ratio selection or sex-antagonistic selection: the new system of heterogamety should, in both cases, be independent of that preceding the transition, and will specifically depend on whether the initial sex bias was towards males or females³¹,³² or the new autosomal sex-antagonistic mutation is male- or female-beneficial²⁴,³³. The relative importance of these driving forces for turnovers is hotly debated, but by testing for the predicted patterns of heterogamety outlined above, we can gain insights into which forces are the dominant drivers of sex-chromosome turnovers in a given system.

One particularly intriguing observation from studies documenting sex-chromosome turnovers is that some specific genome regions are used for sex determination more often than others³⁶. For instance, two species of turtles have independently recruited the same autosome for sex determination³⁷; two lineages of medaka fishes have independently co-opted *Sox3* as their sex determiner³⁸; and in sticklebacks, LG12 has been linked to sex at least twice independently³⁹. Several other examples exist which, together, imply that certain regions of the genome are predisposed to recruitment for sex determination, because they harbour important genes with the potential to determine sex³⁸,⁴⁰, or to fix sexually antagonistic alleles³⁹. However, to our knowledge, no statistical support has yet been provided to support this claim, likely due to the low number of species available to test it.

In this study, we examine the frequency and patterns of sex-chromosome turnover among 28 species of true frogs (Ranidae) to test whether heterogamety is more likely to be preserved or switched during turnovers and to test whether some chromosomes are indeed more likely to be recruited for sex determination than others. This family is thought to have a very high rate of sex-chromosome turnover, and anecdotal evidence points to some chromosomes being more likely to be recruited for sex determination than others¹⁴, making it an ideal system in which to study these questions. However, as previous work was phylogenetically blind, it has not been possible to distinguish the patterns of sex-chromosome recruitment from phylogenetic correlation. Here, we use Restriction site Associated DNA sequencing (RADseq) to search for sex-determination system and the identity of the sex chromosomes in 19 frog species. We couple this new data with that from ground-truthed literature to assess patterns of turnovers across a total of 28 species in a phylogenetic framework. Our results show firstly that Ranidae exhibit an extremely fast rate of sex-chromosome turnover, and that some chromosomes are indeed significantly more likely to be recruited for sex determination than others. Secondly, we find that in all but one case, male heterogamety is conserved during turnovers, suggesting that mutation-load selection is an important driving force in this system. Finally, we show that, in stark contrast with the conventional model of sex-chromosome evolution, frog sex chromosomes rapidly diverge along almost their entire length due to the highly reduced rate of recombination in males, which further supports a role for the accumulating mutation load and might explain why Ranidae undergo such frequent sex-chromosome turnovers relative to other taxa.

## Results

### Identifying sex-linked markers and sex-determination systems.

In total, we produced RADseq libraries for 19 species of Ranidae (Supplementary Table 1) and also re-analyzed the previously published RADseq dataset from *Rana arvalis*⁶⁸. The final number of assembled and retained RADtags per species dataset ranged from 5251 to 146,854 (Mean = 61,152, SD = 27,485), which contained between 2351 and 71,462 Single Nucleotide Polymorphisms (SNPs) (Mean = 29,410, SD = 21,947) and the mean
locus read-depth per dataset was 13.9–36.8 (Mean 23.8, SD = 6.7).

Our screens for sex linkage yielded markers that fit the expectations for either an XY or a ZW system in all datasets, pointing to the widespread occurrence of false positives. We therefore used an in silico approach to test whether observed numbers of putative XY or ZW markers differed from random expectations, by permuting sex across samples 1000 times and running tests for sex-linked markers on each permutation (Supplementary Note 1). The resulting null distribution was then compared with the number of sex-linked markers identified under real sex assignments. Of the 20 species for which we analyzed RADseq data, 12 had sex-linked marker sets that passed these sex-permutation validation tests for at least one method (Supplementary Fig. 1, Supplementary Table 1). In several cases, however (see *R. montezumae*, *R. arvalis*, *R. italica* and *R. kukunoris*, Supplementary Fig. 1), sets of sex-linked markers supporting both XY and ZW systems passed validation in the same species. In *R. montezumae*, *R. italica*, and *R. kukunoris*, there were enough of such markers that their position could be found in the reference genome (see below). In all three cases, both XY- and ZW-like markers aligned to the same chromosome (Supplementary Fig. 2a). As it is biologically implausible that both XY and ZW sex-determination systems simultaneously exist on the same chromosome, we hypothesized that the ZW signal actually resulted from markers that were either X-limited (resulting in male hemizygosity) or had an allele specific to the paternal X chromosome. Such loci can produce almost exactly the same signal as ZW markers. For these four species, where sample availability allowed, we tested this hypothesis with several supplementary analyses. The results of these tests alongside independent lines of evidence, which are detailed in Supplementary Note 2 and Supplementary Fig. 2, allowed us to be confident that these species do indeed exhibit XY heterogamety.

The final number of sex-linked markers confidently identified per dataset varied greatly, ranging from 34 to 1925 (Supplementary Table 1). Only these markers were carried forward in the analyses and are shown in bold in Supplementary Table 1. Failure to find validated sex-linked marker sets in the eight remaining species (five of which, from our literature search, have been shown to be male heterogametic) could have been due to a number of factors, including small sample sizes (Supplementary Note 3, Supplementary Fig. 3c), the presence of several Y haplotypes in the male data (Supplementary Note 3, Subsample 2 in Supplementary Fig. 3b), extremely undifferentiated sex chromosomes (which can only be identified using families with sexed offspring) or non-genetic sex determination⁴¹.

**Identifying the sex chromosome.** In order to identify the chromosome to which sex-linked RADtags belonged, we first anchored scaffolds from an existing fragmented *R. temporaria* genome assembly to the high-quality assembly of *Xenopus tropicalis* via intermediate alignment to high-density linkage maps from males of six *R. temporaria* families (see details in Supplementary Note 4 and Supplementary Fig. 4). Linkage maps contained a total of 15,313 SNPs and allowed us to assign 408.14 Mb (9.07%) of the *R. temporaria* assembly to *X. tropicalis* chromosomes (Supplementary Table 2). Aligning RADtags from our study species to these chromosome-anchored *R. temporaria* scaffolds allowed us to identify the chromosome pair to which they belonged. To control for mis-assembly, mis-alignment or any biases towards particular chromosomes in the RADseq data, we also aligned 1000 randomly chosen subsets of markers (each of the same size as the set of sex-linked markers) to the genome for each dataset. We only accepted a sex chromosome

identification if the number of sex-linked markers which aligned to that chromosome was significantly higher than the 99th percentile of the distribution of the randomly chosen marker subsets. Of the 12 species for which sex-linked markers were confidently identified, seven species passed these alignment validation tests (Supplementary Fig. 5, Supplementary Table 1). In these species, the number of sex-linked markers that aligned to the sex chromosome ranged from 16 (*R. dalmatina*) to 100 (*R. italica*), and in all cases, only one of the 13 chromosomes was identified as being sex-linked. We note that if a set of sex-linked markers met our validation conditions, this further increased our confidence in the sex linkage of these markers, as the likelihood of false positives clustering on a single chromosome (only when using correct sex assignments) is infinitesimal.

The reasons for the unsuccessful alignments in the five remaining species for which sex-linked markers were available were either a high divergence time between that species and *R. temporaria* (e.g. *Pelophylax perezi*, approx. 55 M years diverged) and subsequently a low overall alignment rate, a low number of sex-linked markers identified (e.g. *R. latastei*, 94 RADtags) or a combination of the two (e.g. *R. pipiens*, *R. tarahumarae*) (Supplementary Fig. 6). However, in one species, *R. kukunoris* (Muyu population), the reason was likely a high number of false-positive sex-linked markers in the dataset (Supplementary Figs. 1, 6).

Across the seven species for which a sex-linked chromosome could be reliably identified, three chromosomes out of the possible 13 were used as the sex chromosome, Chr01 (*R. arvalis*, *R. italica*, *R. japonica*, *R. montezumae*), Chr03 (*R. iberica*) and Chr05 (*R. dalmatina*, *R. kukunoris*) (Fig. 1, Supplementary Fig. 5). The results from the literature search, confirmed by mapping of enzyme or microsatellite markers to the *X. tropicalis* genome, provided sex-chromosome identities for 12 additional species and, importantly, also corroborated the finding that the same three chromosomes are more commonly used for sex determination than any others (Fig. 1, Supplementary Table 3). Of the four species that were present in both our study and those from the literature, we were only able to identify the system of heterogamety in two species (*R. japonica* (East) and *R. pipiens* (West)) and the sex chromosome in one species (*R. japonica* (East)). In all cases they agreed. This literature search also revealed several cases of intraspecific polymorphism. *R. temporaria* mostly uses Chr01 throughout its distribution range, but Chr02 has also been shown to segregate with sex (along with Chr01) in one Swedish population⁴². In *R. pipiens*, both Chr02 and Chr05 have been described as sex-linked in different populations, representing an intraspecific turnover between lineages of Eastern and Western USA⁴³. Other examples of intraspecific turnover have been described in *Pelophylax porosus*, which uses Chr05 in the Okayama race (subsp. *porosus*) and Chr03 in the Nagoya race (subsp. *brevipodus*), and *R. japonica*, which uses Chr01 in Western Japan but Chr03 in Eastern Japan. Interestingly, no linkage was found between sex and any of these two chromosomes in one *R. japonica* population from Akita (Northern Honshu), suggesting another, as yet undetected intraspecific turnover²⁰. Finally, Chr08 is potentially involved in three turnover events in *Glandirana rugosa*, the first being a turnover to Chr08 from the ancestral state, and following that, two independent homologous turnovers from an XY to ZW sex-determination system¹⁴.

**Rate and patterns of sex-chromosome turnovers in Ranidae.**

We placed the sex-chromosome identities for each species onto a phylogenetic tree produced by combining data from Yuan et al.⁴⁴ and Pyron and Wiens⁴⁵. The topology and node dates of our tree
<figure><img src="image_3.png" /></figure>
**Fig. 1** Sex chromosome turnovers across 28 true frog species. Sex-determination system and sex-chromosome identities come from both RADseq (Supplementary Table 1) and literature (Supplementary Table 2) data. Karyotype (top left, chromosomes not to scale) shows the number of species using each chromosome for sex determination and colours correspond to arrows, node pie charts and tips. Coloured arrows show the branch on which inferred turnovers occur based on the stochastic mapping analyses (Supplementary Fig. 8) and the pie charts at nodes represent the proportion of simulated trees in stochastic mapping with each of the states at that node. Tips with two colours represent intraspecific turnover events, with the transition described by the coloured arrow after the species name. Grey represents unknown sex chromosome identities in both tips and turnover arrows 1 and 2 and question marks at tips represent an unknown system of heterogamety. As there was not enough high-quality sequence information for *P. porosus*, it could not be included in the phylogenetic reconstruction or stochastic mapping. Its position here is inferred from that in ref. ⁴⁵

conformed very closely to those of both previous studies (Supplementary Fig. 7). We estimated ancestral states throughout the tree using stochastic mapping analyses, which allowed us to infer a minimum of eight sex-chromosome turnover events within the tree (Supplementary Fig. 8). When combined with the three intraspecific sex-chromosome turnover events in *R. japonica*, *R. pipiens*, and *P. porosus* and the two independent XY to ZW turnovers in *G. rugosa*, this equates to 13 sex-chromosome turnovers within approximately 55 M years (Fig. 1). In Supplementary Note 5, we outline an approach to calculate an approximate rate of transitions throughout our tree, which we envisage can be used to compare rates across different studies. Using this approach, we estimate a rate of approximately 0.02 turnovers per million years. This translates into 1 turnover for every 50 M years of independent evolutionary time (i.e. we would expect to see one turnover between two species which were 25 M years diverged from one another).

Although our ancestral states reconstruction was not able to confidently infer the state for the root of our phylogeny, it points strongly to Chr01 being the ancestral state for the *Rana* genus, which diversified approximately 40 M years ago. From this state, there have been two independent turnovers to Chr03 and four independent turnovers to Chr05. Ancestral state reconstruction in *Pelophylax* was based only on two species (i.e. too low to confidently infer ancestral states), but the sex-linkage of Chr03 in this clade indicates at least one more independent turnover event to Chr03 (plus one additional transition to Chr05 within *P. porosus*). Due to the ambiguous state at the root of the tree, the number and placement of these turnovers is unclear. Six hypotheses exist for the turnovers at this point in the tree (between nodes a, b and c, Supplementary Fig. 9); the most parsimonious (requiring the fewest turnover events) is represented on Fig. 1 (turnovers 1 and 2).

The number of times each chromosome was recruited for sex determination was highly variable. Chr01 was used most often (eight species), but in all cases was inferred to be the ancestral sex chromosome. Chr03 and Chr05 were both used by five species; however, Chr05 showed five independent recruitments, whereas only two or three (see turnover 1) were inferred for Chr03. Chr02 was used by two species, which represent two independent recruitments, and Chr08 was used by one species but, as mentioned above, is involved in three independent turnovers within the *G. rugosa* lineage.

As chromosome size varies greatly in Anuran genomes, we tested whether this pattern of sex-chromosome recruitment was consistent with expectations of a random model, taking into account the number of genes on each chromosome. Regression analysis showed no significant relationship between the number of times chromosomes had been recruited and the number of genes they contained ($R^2 = 0.028, p = 0.64$). By taking 1000 random samples of 13 genes (representing the 13 independent recruitment events) from the genome in silico, we show that the observed pattern of recruitment represents a significant departure from random expectation for Chr05, which has been independently recruited many more times than expected ($p < 0.001$;

Fig. 2). Furthermore, despite being the fourth largest chromosome in the genome, Chr04 was never found to be sex-linked in the species examined here. Taken together, these results clearly show a non-random pattern of recruitment of chromosomes for sex determination. Heterogamety shifts during transitions were not random either—eleven out of 13 transitions kept the ancestral system of male heterogamety. This differs significantly from binomial expectations when assuming an equal probability for both types of transitions ($p < 0.001$), but not when assuming a 3/4 probability of maintaining heterogamety, as expected under drift-induced transitions ($p \sim 0.1$).

**Patterns of sex-chromosome differentiation in Ranidae.** As our RADseq data provided multiple sets of sex-linked RADtags, we took the opportunity to examine their distribution along the sex chromosome, in an effort to elucidate the drivers of sex-chromosome turnovers among these species. As the previously anchored scaffolds were aligned to male linkage maps during the anchoring procedure, order within linkage groups was not reliable (because recombination is rare in male meiosis). We therefore created another linkage map combining data from female recombination patterns only across our six *R. temporaria* families (Supplementary Note 4, Supplementary Fig. 10). This linkage map contained 10,853 RADtag markers, confidently identified 13 linkage groups as expected (Supplementary Fig. 10) and showed the typical drastic difference between male and female recombination patterns previously observed for this species⁴¹,⁴⁶. Expanding upon this previous work, we show explicitly that almost the entirety of the female linkage group falls into the non-recombining section of the homologous male linkage group (Supplementary Fig. 10), implying that almost the entire length of male chromosomes did not recombine in these six families (but recombined more than in females at chromosome tips). This map allowed us to order 328 Mb (7.31%) of the *R. temporaria* genome. We then aligned sex-linked RADtags from the seven species for which we could confidently identify a sex chromosome to these ordered scaffolds and plotted their relative positions on each linkage group (Fig. 3). In all cases, sex-linked RADtags mapped across almost the entirety of the linkage groups, with no clear islands of divergence that could be associated to a sex-determining region.

## Discussion

We have shown here that the rate of turnovers in true frogs is extremely fast, with at least 13 turnover events observed among 28 species within approximately 55 M years. In some cases, this has resulted in sex chromosomes that are likely less than 3.6 M years old (between *R. blairi* and *R. berlandieri*), or even younger in the case of intraspecific turnovers (e.g. *R. japonica*, *R. pipiens*, *G. rugosa* or *P. porosus*). This is, to our knowledge, the fastest rate of sex chromosome turnover observed to date. Several other systems have been identified as exhibiting frequent sex-chromosome turnovers, for example the genus *Oryzias* (seven different sex chromosomes identified among 14 species, which
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Non-random recruitment of chromosomes for sex determination. Grey bars represent the number of genes per *X. tropicalis* chromosome. Violin plots show the expected number of times each chromosome would be recruited if genes were recruited for sex determination at random (1000 replicates). Black dots represent the number of times each chromosome was observed as recruited from our data</figcaption></figure>

diverged approximately 65 M years ago<sup>12,25,47)</sup> and Salmonids (six different sex chromosomes (often fused with autosomes) among eight species of salmon, trout, charr and whitefish<sup>48,49</sup>). However, these systems lack ancestral state reconstruction analyses and thus, the number and phylogenetic position of turnovers is not currently known. This unfortunately precludes a formal comparison between the rate observed in Ranidae to that of the other systems. However, we have devised an approach (Supplementary Note 5) that we envisage can be used to compare average rates of turnover among studies. Therefore, we hope that future characterization of transitions using ancestral state reconstruction in other systems will allow for this comparison, which would be highly informative for identifying the forces responsible for turnovers.

The rapid rate of turnover in Ranidae observed here would undoubtedly help to maintain sex chromosome homomorphy. This is, however, only a minimum estimate of turnover rate; it is likely that better knowledge of ancestral states would reveal even more turnovers than we are able to characterize here. Furthermore, our approach in principle cannot detect situations where the new sex determiner arises on the already sex-linked linkage group ('homologous' turnovers) unless it was accompanied by a change in heterogamety (XY to ZW), as was the case in *G. rugosa*. Also, we cannot rule out that several turnovers have occurred on the same branch in our phylogeny. Both phenomena would act to increase the observed rate of turnovers.

In order to infer the major forces driving the high rate of turnovers in this system, we asked whether patterns of heterogamety were conserved during turnovers. Of the 13 turnover events observed here, 11 of them were XY to XY, with the two XY to ZW turnovers both being in *G. rugosa*. Even if one or more of the four species for which no sex-determination system could be found with either RADseq or in existing literature (*R. macrocnemis*, *R. chensinensis*, *R. ornativentris* and *R. uenoi*) were in fact ZW, this would still leave a clear bias towards XY heterogamety in this family. This bias runs against expectations from models where either sex-ratio selection or sex-antagonistic selection drive turnovers<sup>31–33,50</sup>. Two counterarguments could be made here. First, mechanisms are not mutually incompatible; if sexual antagonism was driving turnovers but mutationally loaded Y chromosomes cannot be fixed, then heterogamety would be preserved as well, but the causal force would not be the deleterious load. Second, sex-antagonistic selection could also bias turnovers towards XY systems in Ranidae: given the much reduced recombination in males, a male-determining mutation occurring on an autosome is more likely to benefit from linkage with a sex-antagonistic gene than is a female-determining

mutation<sup>50,51</sup>. Although both arguments are correct in principle, the point must also be made that sex-antagonistic selection is not expected to generate the continuous cycles of turnovers as documented here. An autosomal male-beneficial mutation might indeed trigger an initial XY-to-XY transition but should strongly oppose any further transition once sex-linked<sup>34,35</sup>. Furthermore, genomic investigations do no support a role for sex-antagonistic genes in sex-chromosome evolution in *R. temporaria*<sup>52</sup>. Overall, data in hand seem more consistent with a role for mutation-load selection or possibly drift (or a combination of the two). It is worth noting that the only known exceptions to the general pattern (both in *G. rugosa*) have explicitly been assigned to sex-ratio selection: the ZW races likely result from two hybridization events between highly diverged XY lineages of *G. rugosa*, and lab crosses between these parental lineages indeed show a male bias in the progeny, which should suffice to select for a dominant feminizing mutation<sup>14</sup>. If this interpretation is correct, then transitions left to be explained all keep the original pattern of male heterogamety, which further increases the likelihood of the mutation-load model over the genetic-drift model.

If mutation load is indeed the driving force for the preservation of XY heterogamety in Ranidae, then the sex-specific patterns of recombination might help explain why the rate of sex-chromosome turnover is so high in frogs relative to many other taxa. Frogs exhibit extreme patterns of heterochiasmy (different recombination rates between males and females): the ratio of female-to-male map length has been estimated at ~70 in *R. temporaria* based on microsatellite markers<sup>53</sup>, a remarkable result given that values normally range 0.5–2<sup>54–56</sup>, with a few outliers in fishes (up to 8.26 in the Atlantic salmon *Salmo salar*<sup>57</sup> and down to 0.135 in the Japanese flounder *Paralichthys olivaceus*<sup>58</sup>). From our present RADseq data, this strong sex difference stems from the fact that all chiasmata in male meiosis form at chromosome extremities (Supplementary Fig. 10). A similar pattern was recently described in stickleback, although the clustering of chiasmata towards the telomeres was not as extreme<sup>51</sup>. As meiotic recombination depends on phenotypic sex (not on genotypic sex), X and Y chromosomes may still recombine occasionally in the rare sex-reversed XY females that occur when alleles at the sex locus show incomplete penetrance<sup>59</sup>, preventing X-Y differentiation. This possibly accounts for our inability to identify sex-linked markers in some of our samples. However, recombination will immediately cease along almost the entire length of the X and Y chromosomes as soon as full-penetrance alleles are fixed, which is clearly the case in some *R. temporaria* populations<sup>42,60,61</sup> as well as in populations of other Ranid species that display extensive X-Y differentiation (Fig. 3). The sudden linkage of a very
<figure><img src="image_5.png" /><figcaption>**Fig. 3** Chromosome-wide distribution of sex-linked RAD markers in seven true frog species. Density plots to the left of linkage maps show the marker density along the linkage group. Coloured bars are histograms of the number of sex-linked markers along the linkage group, with bins of 2 cM. Picture credits: Jan Jezek (R. arvalis), Astolinto (R. italica), Yasunori Koide (R. japonica), Jim Rorabaugh (R. montezumae), DLJ (R. iberica), NR (R. dalmatina), AB (R. kukunoris)</figcaption></figure>

high number of genes is expected to generate particularly strong Hill-Robertson interferences. As a result, deleterious mutations will immediately begin to accumulate, and will do so along the length of the entire chromosome (except at the tips). This process is in stark contrast with the progressive expansion of the non-recombining region outwards from the sex determiner posited by the canonical model of sex-chromosome evolution, and does not require the appearance of sexually antagonistic genes to favour recombination arrest along the sex chromosome. In the absence of sexually antagonistic alleles, the sex-limited chromosome is expected to very quickly become less fit than any potential proto-sex-chromosome that would arise. It then follows that mutation-load selection would quickly favour a turnover as soon as a new option for sex determination is available (a particularly efficient version of the 'hot-potato' model of sex-chromosome turnover³⁵). It would be interesting to quantify deleterious mutations on differentiated Y chromosomes and examine sex-chromosome differentiation and turnover in other systems with similar patterns of heterochiasmy, to see if this hypothesis holds.

A central question of this study was: are some chromosomes more likely to be recruited for sex determination than others? We show here that recruitment departs significantly from random expectations in true frogs. Of the 13 chromosome pairs, only five are used for sex determination, with three of them being recruited multiple times independently—most notably Chr05 which has been recruited an astonishing five times. Chr01, used as sex chromosome in the majority of species investigated here, has also been co-opted independently in species from two highly diverged clades of frogs (Hylidae and Bufonidae)⁶². This raises the question: why are some chromosomes, for instance Chr01, Chr03 and Chr05, more likely to be recruited for sex determination than others? Is it that some chromosomes are better at sex determination than others⁶³? One explanation for this could be the presence of particular genes on these chromosomes that play

important roles in the sex-determination cascade. Indeed, several genes (or their paralogs) have been recruited independently more than once as master sex-determiners, suggesting that they are particularly good at fulfilling this role. These so called 'usual suspects'⁶⁴ include *Sox3*, which maps to Chr08 (involved in three transitions within *G. rugosa*). *Sox3* determines sex in several medaka (*Oryzias*) species³⁸ and is a strong candidate in *G. rugosa*⁶⁵, while its paralog *Sry* is the master sex-determining gene in therians⁶⁶. The transcription factor *Dmrt1*, which plays a crucial role in the sex-determination cascade throughout all animals (and maps to Chr01), has been independently co-opted for sex-determination in birds⁶⁷ and some fishes⁶⁸ and it is a strong candidate in *R. temporaria* and in some tree frogs (*Hyla*) ⁶⁹,⁷⁰. Furthermore, its paralogs *DmY* and *DmW* determine sex in some medaka fishes (*Oryzias*) and clawed frogs (*Xenopus*), respectively³⁸,⁷¹. These paralogs have functions related to those of their parent genes (e.g. the female determiner *DmW* is thought to block the binding of *Dmrt1* to its target sequence, thus suppressing the male determining cascade⁷²) lending even more credence to the idea that the specific functions of these genes make them particularly proficient at the role of master sex-determiner. Interestingly, Chr05 in frogs harbours *Foxl2* (forkhead box L2), also one of the most important genes in the vertebrate sex-determination cascade⁷³. It codes for a transcription factor essential for ovarian development and has also been implicated in the suppression of testis formation. *Foxl2* interacts directly with *Dmrt1*; the male-determining *Dmrt1* allele blocks the expression of *Foxl2*, and in turn the development of ovaries, thus producing males. Indeed, studies in both tilapia⁷⁴ and zebrafish⁷⁵ have shown that the knockdown of *Foxl2* results in female to male sex reversals. Thus, a plausible hypothesis for the recurrent recruitment of Chr05 for sex-determination is that a loss-of-function mutation in *Foxl2* uncouples its activity from that of *Dmrt1*, allowing it to take over as the primary male determining gene.
Accordingly, Chr05 might not just be good at the role of sex chromosome, but more specifically it might be good at taking that role from Chr01.

Still, we cannot rule out a translocation of the same sex
determiner to new chromosomes. In the case of the jumping sex-
determining gene in salmonids, Lubieniecki et al.⁴⁸ found that
three chromosomes, all of which are sex-linked in some popu-
lations, harboured sequences that were homologous to the regions
flanking the sex determiner. These regions facilitated the trans-
location of the sex determiner to these locations, producing a
recruitment bias towards some chromosomes (as well as a pattern
of strongly preserved male heterogamety) similar to the result
seen here. We note, however, that the two transcription factors
identified as serious candidates in our dataset (respectively Sox3
in G. rugosa and *Dmrt1* in R. temporaria) did not jump in this
specific instance (*Dmrt1* is autosomal in G. rugosa, and *Sox3* is
autosomal in R. temporaria). We also note that convergent
recruitment of the same genes should not be considered as an
alternative to the hot-potato model to account for the main-
tenance of heterogamety. We rather see the accumulation of
deleterious mutations as the ultimate cause (selective force), and
the convergent recruitment of genes from the sex-determination
cascade as the proximate cause. Convergence per se does not
impose maintenance of heterogamety, as mutation of any gene in
the male cascade can be either masculinizing or feminizing (a
mutation upregulating *Dmrt1* expression will be masculinizing,
while a mutation downregulating its expression will be feminiz-
ing). Future work should concentrate on validating candidate sex-
determining genes in this family.

In conclusion, we show here through the combination of novel
and existing data that sex-chromosome turnover among frogs
occurs at a very high rate. This may be a result of the extreme
heterochiasmy that characterizes frogs, which should induce a
rapid accumulation of deleterious-mutation load, and thereby
quickly select for a turnover (a frog version of the ‘hot-potato’
model³⁵). This would account for the consistent maintenance of
male heterogamety despite frequent transitions. Intriguingly,
these turnovers converge on only a few chromosomes, implying
that they are particularly adept at the role of sex determination,
perhaps due to specific genes that they harbour.

**Methods**

**Sample collection and RADseq library preparation.** In total, we generated

RADseq data for 19 species of true frogs (Family: Ranidae). The taxonomy is still debated within this family⁷⁶, so here we opt to use the simpler early nomenclature also recently used by Yuan et al.⁴⁴. Under this nomenclature, these species fall into two genera, *Rana* s.l. (sometimes split into *Rana* and *Lithobates*) and *Pelophylax*. We obtained male and female samples of each species (mean *N*<sub>FEMALES</sub> = 22.2, *N*<sub>MALES</sub> = 21.7), totalling 736 individuals (Supplementary Table 1).

Samples collected for this study were done so under the following permits: *R. berlandieri*, *R. chiricahuensis*, *R. tarahumarae* and *R. yavapaiensis* samples were collected by M.J.S., under the authorities in Arizona Revised Statutes, Title 17, and under Section 6 authorities granted by the Endangered Species Act. *R. sphenocephala* samples were collected by L.N.B. under permits from the FL Fish and Wildlife Conservation Commission (LSSC-09-322). *R. montezumae* samples were collected by C.M.G. under permits from Secretaría del Medio Ambiente y Recursos Naturales (SEMARNAT) (SGPA/DGVS/02919/15). *R. latastei* samples were collected by G.L., W.J.M. and K.G. under permits from Ufficio della natura e del paesagio, Ticino. *R. italica* was collected by A.R. under permits issued by the Italian Ministry of Environment (0003951/PNM of 03/03/2015). *R. dalmatina* was sampled by D.C. with permits from the Ministry of Environment, Italy (Prot.0007727). *R. iberica* was sampled by A.G.N. under permits from the Picos de Europa National Park and the Principality of Asturias (Spain). *P. perezi* was sampled by I.M.S. under permit 10/016997.9/13 from Consejería de Medio Ambiente, Comunidad de Madrid (Spain). *P. saharicus* samples were collected by PAC under permits from Haut Commissariat aux Eaux et Forêts et à la Lutte Contre la Désertification. No permits were needed for sampling of *R. uenoi* (according to Korean law) or for *R. japonica*, *R. ornativentris*, *R. chensinensis* and *R. kukunoris* (according to Chinese law).

On collection, individuals were also phenotypically sexed, either
morphologically, in many cases using the presence or absence of nuptial pads or
vocal sacs, or by dissection and examination of gonads⁷⁷.

DNA was sampled from collected frogs by fixing tissue (tail clips, toe clips, leg muscle, whole tadpoles or whole frogs) in 70-90% ethanol or tissue buffer, or by taking buccal cell swabs. We then extracted DNA using the Qiagen® DNeasy® Blood & Tissue Kit and produced double-digest RADseq libraries for each species according to the protocol described in Brelsford et al.⁷⁸. In brief, we digested genomic DNA using restriction enzymes SbfI-HF and MseI and then ligated adapters containing unique barcodes of 4-8 bases separately to each sample before amplifying the libraries using PCRs of 20 cycles. We then pooled PCR products and performed a gel-based size selection, isolating fragments of approximately 400-500 bp. The finished RADseq libraries were single-end sequenced on an Illumina HiSeq 2500.

**Raw data processing and SNP calling.** Illumina raw reads were quality checked using FASTQC v0.10.1⁷⁹ and were demultiplexed by individual barcode using the *process_radtags* module of STACKS v1.48⁸⁰. STACKS modules *Ustacks*, *Cstacks*, *Sstacks* and *Populations* were then run separately. For all datasets, we ran initial tests using multiple values for the core STACKS parameters (*Ustacks*: -M, -m and *Cstacks*: -n), but in all cases, the default values provided a good balance between data quantity and quality, thus these defaults were used for all species. In the final analyses, SNP markers were retained if they were present in at least 75% of both males and females. Additionally, we specified a minimum minor allele frequency of 0.05 and a maximum observed heterozygosity threshold of 0.75. The latter acts to remove over-merged stacks resulting from repetitive elements, as each paralogous copy of the repeat will look like a separate allele in such stacks, these false-loci are likely to have a much higher number heterozygous calls across individuals compared to either Hardy-Weinberg expectations (Population data) or Mendelian segregation expectations (Family data, except for loci homozygous for different alleles between the two parents, which are thus not sex-linked). Note that there is the possibility of removing sex-linked markers using such a filter if the number of the heterogametic sex makes up >75% of the samples. However, this was not the case for any of our sample sets. The number of high-quality SNPs retained per dataset ranged from 2351 to 71,462 (Mean = 31,515, SD = 22,021).

**Identifying sex-linked markers.** Final SNP datasets for each species were bioinformatically screened for signatures of sex-linkage, for both male (XY) and female (ZW) heterogametic systems, using the three approaches previously described in Brelsford et al.⁸¹. These approaches are described briefly below. For the sake of simplicity, we refer to only XY systems, but in all cases the reverse is true for ZW systems.

Approach one screens for loci with a Y-specific SNP. In an XY system, these SNPs will be present at a frequency of 0.5 in males and 0 in females. We therefore consider a locus to be sex-linked if the X-allele frequency is ≥0.95 in females and ≥0.4 and ≤0.6 in males, allowing a margin of error for a few incorrect genotype calls or phenotypic sex assignments.

Approach two also screens for Y-specific alleles, but instead uses the differences in patterns of heterozygosity between males and females. Under this approach, a locus is considered sex-linked if it is homozygous in all females and heterozygous in at least half of the male samples. This method is likely to return many of the same sex-linked markers as the first approach but is more stringent in its conditions. For example, an allele frequency of 0.5 in males may be possible by chance if the sample set is small, or if there is population structure within males leading to consistently fixed loci among them. Thus, it is expected that approach two will be less sensitive to false positives in such situations.

Approach three screens for sex-limited RADtags that are specific to the Y chromosome. In this approach, RADtags are assessed for their presence or absence in one sex and are considered sex-linked if they are completely absent in the homogametic sex and present in at least half of the heterogametic sex. Several mechanisms can result in such sex-specific loci: a deletion on the X chromosome, a mutation in the restriction site on the X chromosome (leading to a sex-specific null allele), a mutation on the Y creating a novel Y-specific restriction site, or under-merging of a RADtag in Stacks due to high divergence between X and Y alleles. In the latter case, only the Y allele will be sex specific and thus, can still be used as a sex-linked marker.

In several species, samples were collected from multiple populations or lineages. It was therefore necessary to split these data into subsets and screen each one for sex-linked markers separately (see Supplementary Table 1), in order to control for the potential noise produced by population structure.

**Identifying the sex chromosome.** In order to identify the sex chromosome in each species, it was necessary to locate the RADtags identified as sex-linked on a reference genome. The most reliable chromosome-level assembly in frogs is that of *X. tropicalis*⁸² and, being an outgroup, it provides a good reference with which to unambiguously name chromosomes. However, due to the amount of divergence between Ranidae and *X. tropicalis* (210 M years⁸³), mapping rates of 92 bp RADtags to this reference genome are almost zero. We therefore anchored as many scaffolds as possible from an existing fragmented *R. temporaria* genome assembly⁴¹ to the *X. tropicalis* assembly. We created linkage maps using RADseq data from six
*R. temporaria* families from Rodrigues et al.⁶¹ (60–90 offspring per family) and used these, alongside *Nanorana parkeri* and *X. tropicalis* genomes, to assign scaffolds to their homologous *X. tropicalis* chromosome (henceforth referred to as chromosome-assigned scaffolds). Note that synteny is known to be highly conserved between *N. parkeri*, *X. tropicalis* and *Ranidae*⁴⁶,⁶²,⁸⁴,⁸⁵. Linkage mapping and assignment of scaffolds to chromosomes is detailed in the Supplementary Note 4 and Supplementary Fig. 3.

RADtags identified as sex-linked were aligned to the chromosome-assigned scaffolds in order to identify the chromosome to which they belonged. All alignments mentioned above were performed using blastn v2.3.0⁸⁶, with hits only retained if their e-value was below $1 \times 10^{-20}$ and, in the case of multiple matches, at least five orders of magnitude lower than that of the next best hit (to account for the repetitive nature of Anuran genomes⁸⁴).

We combined our new data with already-published sex-chromosome identities for other species from the literature (Supplementary Table 3). However, as many studies use cytogenetic methods, which are not completely reliable in the case of such homomorphic sex chromosomes, and to avoid confusion from ambiguous chromosome naming, we included only studies that used sex-linked allozymes to identify the sex chromosome. The location of the enzymes used could then be reliably confirmed by searching them against the well annotated *X. tropicalis* genome (Xenbase.org⁸²). Of the 19 species for which new data were generated in this study, five were present in the retained studies from our literature search (*R. arvalis*, *R. japonica*, *R. berlandieri*, *R. sphenocephala*, *R. pipiens*).

**Validating sex-linked marker and sex-chromosome identification.** When screening large SNP datasets for sex-linked loci, it is imperative to account for the possibility of false positives, i.e. loci which fit the expected pattern of sex-linkage due to artifacts or simply by chance, without being truly sex-linked. Here, we used a permutation approach to increase the confidence around sex-linked marker and sex-chromosome identification, while accounting for the sample size, sex-skew and biological characteristics of each dataset. For a given dataset, we first used true phenotypic sex assignments, and the three approaches detailed above, to find putative sex-linked markers. We then calculated the empirical null distribution for the rate of false positive discovery by permuting male and female assignments randomly across the sample set 1000 times, each time screening for sex-linked markers with all three approaches. We retained only sex-linked marker sets where the number found using the real sex assignments fell outside of the 99th percentile of the null distribution, i.e. had a probability of <0.01 of occurring by chance.

We also used this permutation approach to identify which dataset properties were the most important predictors of the number of sex-linked markers or false positives found. To do this, we extracted various subsamples from the already published RADseq dataset from *R. arvalis*⁸¹. We chose this dataset for several reasons: first, the sex-linked markers in this dataset (identified using the same three approaches as used here) have been validated via the use of PCR. Second, the original sample size (29 males and 19 females) was large enough to allow for several rounds of downsampling. And third, exploratory PCA clustering of individuals using only the sex-linked markers shows evidence for two major Y haplotypes present among males (referred to here as Hap 1 and Hap 2). This dataset therefore allowed us to test four potential predictors for the success of sex-linkage screens: (1) the effect of overall sample size (via a series of subsamples from 19 males (M), 19 females (F) down to 5 M, 5 F), (2) a skew in sample size towards males (18 M, 9 F), (3) a skew towards females (9 M, 18 F), and (4) the effect of having multiple sex-chromosome haplotypes in the same dataset (6M<sup>Hap1</sup>, 6M<sup>Hap2</sup>, 12 F). The latter factor has not previously been considered when assessing the stringency of screens for sex-linked markers.

Putative sex-linked marker sets that fell outside of the 99th percentile of the empirical null distribution were then used in the next step of the analyses, finding the location on the reference genome using the chromosome-assigned scaffolds. However, as the presence of false positives in sex-linked marker sets could not be ruled out, and because the amount of divergence between the target species and *R. temporaria* chromosome-assigned scaffolds was often high (5–41 M years), we expected some noise when aligning markers stemming from either false mappings, mapping of false positive sex-linked loci which are actually autosomal, or true translocations that have occurred since the divergence of the target species and the *R. temporaria* reference genome. Such errors could lead to false signals of sex linkage for chromosomes other than the true sex chromosome. Thus, it was also necessary to calculate the expected false-positive mapping rate, per dataset and per *R. temporaria* chromosome. To do this, we first mapped the putative sex-linked markers for a given dataset to the ordered *R. temporaria* scaffolds. We then took 1000 random subsamples of RADtags from the Stacks catalogue of loci for that dataset, with each subsample being of the same size as the set of sex-linked markers. Each subsample was mapped to the chromosome-assigned *R. temporaria* scaffolds to give a distribution of the expected mapping rate for each chromosome.

For final sex-chromosome identification, we required at least ten sex-linked markers to align to a single chromosome and the number of sex-linked marker alignments to fall outside of the 99th percentile of the random marker alignment distribution for that chromosome.

**Assessing patterns of sex-chromosome turnover.** In order to test whether sex-chromosome recruitment was random, or if some chromosomes are more likely to

be recruited, it was necessary to place all sex-chromosome identities onto a phylogenetic tree. A phylogeny was therefore produced using a combination of data from two previous phylogenetic studies⁴⁴,⁴⁵. We used all ten loci in ref.⁴⁴, which included four mitochondrial genes (12S, 16S, cytb and nd2) and six nuclear genes (rag1, rag2, bdnf, slc83a3, tyr and pomc). Of the 28 species in our study, 22 were included in ref.⁴⁴. For the remaining six (*R. italica*, *P. lessonae*, *P. perezi*, *P. porosus*, *G. rugosa* and *P. saharicus*), we used data for 12S, 16S and cytb from ref.⁴⁵. Gen-bank accession numbers for all sequences used can be found in Supplementary Table 4.

Sequences were first aligned separately for each gene, using MUSCLE⁸⁷ as implemented in Geneious v11.1.5⁸⁸. Alignments were then manually trimmed, excluding regions with high amounts of missing data or ambiguous alignment. The ten trimmed gene alignments were used to make a supermatrix, which was partitioned according to gene and codon position. Best-fit substitution models were then found for each partition using PartitionFinder2 (Supplementary Table 5)⁸⁹–⁹¹. Yuan et al.⁴⁴ showed that maximum parsimony, maximum likelihood and Bayesian inference all produced very similar trees, thus here we used only a Bayesian approach as implemented in BEAST ²⁹². We also incorporated the four fossil-based node calibrations used in Yuan et al.⁴⁴ in order to date the tree, using a relaxed uncorrelated clock⁹³ and a Yule birth–death model. Trees were calculated using MCMC chain lengths of 100 M, with sampling every 10,000 iterations, ensuring that the effective sample size (ESS) of each parameter in the model was above 200.

Using this phylogenetic framework, we inferred ancestral states using a stochastic mapping approach, implemented in the R package Phytools v0.6-44⁹⁴. The best model for transition rates between these states was identified by comparing the likelihood scores (using Akaike Information Criterion (AIC)) for three different transition rate models, equal rates (ER), symmetrical (SYM), all rates different (ARD). To reduce noise and increase the power of this analysis, we assumed that only the five chromosomes that had been confidently identified as a sex chromosome in at least one of our species could determine sex; thus we considered five possible states in the reconstruction. Because of sex-chromosome polymorphism within some species, we coded the states at the tips as a probability matrix, giving the different states equal probabilities in that species. For species where we could not detect the sex-determination system, all five possible states (sex-chromosomes) were given equal probabilities. Using these inputs, we reconstructed ancestral states along each branch for 1000 simulated trees. At any given point in the tree, the most likely state was the one with the highest number of simulated trees supporting it at that position, thus, turnovers could be placed on branches where there was a switch in the state with the most support. We also repeated this analysis with all ten possible chromosome states included, which confirmed that the use of only five states did not bias ancestral states towards any particular chromosomes.

Finally, using the number of turnovers and the identities of the ancestral/ derived chromosomes for each one, we tested whether or not the observed pattern might be expected under a random recruitment model. We first used an ordinary least-squares regression to test for a correlation between the number of times a chromosome had been recruited and the number of genes on that chromosome (taken from the homologous chromosome of *X. tropicalis*). We then simulated the random expectation for the number of recruitments per chromosome by repeatedly (100 replicates) sampling a subset of *N* genes (*N* = number of sex-chromosome recruitments observed) from the genome and comparing these distributions to the observed number of recruitments for each chromosome.

**Assessing patterns of sex-chromosome divergence.** As our RADseq data provided multiple sets of sex-linked RADtags, we took the opportunity to examine their distribution along the sex chromosome, in the hope of elucidating the drivers of sex-chromosome turnovers among these species. Recombination occurs much less frequently in male frog meiosis, and because Y chromosomes are almost always found only in males, we hypothesized that, instead of the Y chromosome diverging from the X only at a small non-recombining region around a sex-determining, or sex-antagonistic genes, the frog Y chromosomes would diverge from the X across almost their entire length, except for the tips of chromosomes where recombination is known to occur⁴¹. To test for this expected pattern of divergence between gametologs, we first produced an additional linkage map in which we combined female-informative markers from all six *Rana temporaria* families into a single map using Lep-Map³⁹⁵ (see 'Linkage mapping' in Supplementary Note 4). While this approach resulted in a map with fewer markers than the six separate male maps produced with MSTmap, using only female linkage information provided more accurate relative marker positions within linkage groups. We used this map to order *R. temporaria* scaffolds along each *X. tropicalis* chromosome, and hereafter refer to these scaffolds as 'ordered scaffolds'. These could then be used to find the relative location of sex-linked markers on the sex chromosomes identified across different species by simply aligning sex-linked markers to the ordered *R. temporaria* scaffolds using blastn.

**Code availability.** All custom scripts and functions used can be found in the github repository https://github.com/DanJeffries/My_misc scripts. Python notebooks detailing the use of these scripts in data exploration and sex-linked marker screens
for all species can also be found on Figshare (https://doi.org/10.6084/m9.fishshare.6949331), please refer first to the README.txt distributed with these notebooks.

## Data availability

All sequencing data used in this study can be found (demultiplexed by sample) on the NCBI Sequence Read Archive under the Bioproject accession PRJNA478189. Stacks intermediate files for each species, linkage maps and alignment and tree files for the phylogeny can be found on Figshare (http://dx.doi.org/10.6084/m9.fishshare.6949331) along with fasta files containing putatively sex-linked markers where identified.

Received: 14 April 2018 Accepted: 29 August 2018

Published online: 05 October 2018

## References

1. Charlesworth, D., Charlesworth, B. & Marais, G. Steps in the evolution of heteromorphic sex chromosomes. *Heredity* **95**, 118–128 (2005).

2. Kirkpatrick, M. How and why chromosome inversions evolve. *PLoS Biol.* **8**, e1000501 (2010).

3. Hill, W. G. & Robertson, A. The effect of linkage on limits to artificial selection. *Genet. Res.* **8**, 269–294 (1966).

4. Muller, H. J. The relation of recombination to mutational advance. *Mutat. Res.* **106**, 2–9 (1964).

5. Bergero, R. & Charlesworth, D. The evolution of restricted recombination in sex chromosomes. *Trends Ecol. Evol.* **24**, 94–102 (2009).

6. Vicoso, B. & Bachtrog, D. Numerous transitions of sex chromosomes in Diptera. *PLoS Biol.* **13**, e1002078 (2015).

7. Vicoso, B., Emerson, J. J., Zektser, Y., Mahajan, S. & Bachtrog, D. Comparative sex chromosome genomics in snakes: differentiation, evolutionary strata, and lack of global dosage compensation. *PLoS Biol.* **11**, e1001643 (2013).

8. Schartl, M. Sex chromosome evolution in non-mammalian vertebrates. *Curr. Opin. Genet. Dev.* **14**, 634–641 (2004).

9. Eggert, C. Sex determination: the amphibian models. *Reprod. Nutr. Dev.* **44**, 539–549 (2004).

10. Bachtrog, D. et al. Sex determination: why so many ways of doing it? *PLoS Biol.* **12**, e1001899 (2014).

11. Beukeboom, L. W. & Perrin, N. *The Evolution of Sex Determination* (Oxford University Press, Oxford, UK, 2014).

12. Kikuchi, K. & Hamaguchi, S. Novel sex-determining genes in fish and sex chromosome evolution. *Dev. Dyn.* **242**, 339–353 (2013).

13. Phillips, R. B. Evolution of the sex chromosomes in salmonid fishes. *Cytogenet. Genome Res.* **141**, 177–185 (2013).

14. Miura, I. An evolutionary witness: the frog *Rana rugosa* underwent change of heterogametic sex from XY male to ZW female. *Sex. Dev.* **1**, 323–331 (2007).

15. Furman, B. L. S. & Evans, B. J. Sequential turnovers of sex chromosomes in African clawed frogs (*Xenopus*) suggest some genomic regions are good at sex determination. *G3 (Bethesda)* **6**, 3625–3633 (2016).

16. Dufresnes, C. et al. Sex-chromosome homomorphy in Palearctic tree frogs results from both turnovers and X-Y recombination. *Mol. Biol. Evol.* **32**, 2328–2337 (2015).

17. Miura, I. Sex determination and sex chromosomes in amphibia. *Sex. Dev.* **11**, 298–306 (2017).

18. Gamble, T. et al. Restriction site-associated DNA sequencing (RAD-seq) reveals an extraordinary number of transitions among gecko sex-determining systems. *Mol. Biol. Evol.* **32**, 1296–1309 (2015).

19. Hillis, D. M. & Green, D. M. Evolutionary changes of heterogametic sex in the phylogenetic history of amphibians. *J. Evol. Biol.* **3**, 49–64 (1990).

20. Sumida, M. & Nishioka, M. Geographic variability of sex-linked loci in the Japanese brown frog, *Rana japonica*. *Sci. Rep. Lab. Amphib. Biol.* **13**, 173–195 (1994).

21. Ezaz, T., Sarre, S. D., O’Meally, D., Graves, J. A. M. & Georges, A. Sex chromosome evolution in lizards: independent origins and rapid transitions. *Cytogenet. Genome Res.* **127**, 249–260 (2009).

22. Woram, R. A. et al. Comparative genome analysis of the primary sex-determining locus in salmonid fishes. *Genome Res.* **13**, 272–280 (2003).

23. Traut, W. & Willhoefft, U. A jumping sex determining factor in the fly *Megaselia scalaris*. *Chromosoma* **99**, 407–412 (1990).

24. Roberts, R. B., Ser, J. R. & Kocher, T. D. Sexual conflict resolved by invasion of a novel sex determiner in Lake Malawi cichlid fishes. *Science* **326**, 998–1001 (2009).

25. Myosho, T. et al. Tracing the emergence of a novel sex-determining gene in medaka, *Oryzias luzonensis*. *Genetics* **191**, 163–170 (2012).

26. van Doorn, G. S. Patterns and mechanisms of evolutionary transitions between genetic sex-determining systems. *Cold Spring Harb. Perspect. Biol.* **6**, a017681 (2014).

27. Bull, J. J. & Charnov, E. L. Changes in the heterogametic mechanism of sex determination. *Heredity* **39**, 1–14 (1977).

28. Veller, C., Muralidhar, P., Constable, G. W. A. & Nowak, M. A. Drift-induced selection between male and female heterogamety. *Genetics* **207**, 711–727 (2017).

29. Saunders, P. A., Neuenschwander, S. & Perrin, N. Sex chromosome turnovers and genetic drift: a simulation study. *J. Evol. Biol.* **31**, 1413–1419 (2018).

30. Cordaux, R., Bouchon, D. & Grève, P. The impact of endosymbionts on the evolution of host sex-determination mechanisms. *Trends Genet.* **27**, 332–341 (2011).

31. Kozielska, M., Weissing, F. J., Beukeboom, L. W. & Pen, I. Segregation distortion and the evolution of sex-determining mechanisms. *Heredity* **104**, 100–112 (2010).

32. Kozielska, M., Pen, I., Beukeboom, L. W. & Weissing, F. J. Sex ratio selection and multi-factorial sex determination in the housefly: a dynamic model. *J. Evol. Biol.* **19**, 879–888 (2006).

33. van Doorn, G. S. & Kirkpatrick, M. Turnover of sex chromosomes induced by sexual conflict. *Nature* **449**, 909–912 (2007).

34. Blaser, O., Grossen, C., Neuenschwander, S. & Perrin, N. Sex-chromosome turnovers induced by deleterious mutation load. *Evolution* **67**, 635–645 (2013).

35. Blaser, O., Neuenschwander, S. & Perrin, N. Sex-chromosome turnovers: the hot-potato model. *Am. Nat.* **183**, 140–146 (2014).

36. Graves, J. & Peichel, C. L. Are homologies in vertebrate sex determination due to shared ancestry or to limited options. *Genome Biol.* **11**, 205 (2010).

37. Montiel, E. E., Badenhorst, D., Tamplin, J., Burke, R. L. & Valenzuela, N. Discovery of the youngest sex chromosomes reveals first case of convergent co-option of ancestral autosomes in turtles. *Chromosoma* **126**, 105–113 (2017).

38. Myosho, T., Takehana, Y., Hamaguchi, S. & Sakaizumi, M. Turnover of sex chromosomes in celebensis group medaka fishes. *G3 (Bethesda)* **5**, 2685–2691 (2015).

39. Ross, J. A., Urton, J. R., Boland, J., Shapiro, M. D. & Peichel, C. L. Turnover of sex chromosomes in the stickleback fishes (Gasterosteidae). *PLoS Genet.* **5**, e1000391 (2009).

40. Ezaz, T., Srikulnath, K. & Graves, J. A. M. Origin of amniote sex chromosomes: an ancestral super-sex chromosome, or common requirements? *J. Hered.* **108**, 94–105 (2016).

41. Brelsford, A., Rodrigues, N. & Perrin, N. High-density linkage maps fail to detect any genetic component to sex determination in a *Rana temporaria* family. *J. Evol. Biol.* **29**, 220–225 (2016).

42. Rodrigues, N., Vuille, Y., Brelsford, A., Merilä, J. & Perrin, N. The genetic contribution to sex determination and number of sex chromosomes vary among populations of common frogs (*Rana temporaria*). *Heredity* **117**, 25–32 (2016).

43. Wright, D. A. & Richards, C. M. Two sex-linked loci in the leopard frog, *Rana pipiens*. *Genetics* **103**, 249–261 (1983).

44. Yuan, Z.-Y. et al. Spatiotemporal diversification of the true frogs (genus *Rana*): a historical framework for a widely studied group of model organisms. *Syst. Biol.* **65**, 824–842 (2016).

45. Pyron, R. A. & Wiens, J. J. A large-scale phylogeny of Amphibia including over 2800 species, and a revised classification of extant frogs, salamanders, and caecilians. *Mol. Phylogenet. Evol.* **61**, 543–583 (2011).

46. Palomar, G. et al. Comparative high-density linkage mapping reveals conserved genome structure but variation in levels of heterochiasmry and location of recombination cold spots in the common frog. *G3 (Bethesda)* **7**, 637–645 (2017).

47. Setiamarga, D. H. E. et al. Divergence time of the two regional medaka populations in Japan as a new time scale for comparative genomics of vertebrates. *Biol. Lett.* **5**, 812–816 (2009).

48. Lubieniecki, K. P. et al. Genomic instability of the sex-determining locus in Atlantic salmon (*Salmo salar*). *G3* **5**, 2513–2522 (2015).

49. Sutherland, B. J. G., Rico, C., Audet, C. & Bernatchez, L. Sex chromosome evolution, heterochiasmry and physiological qtl in the salmonid brook charr *Salvelinus fontinalis*. *G3 (Bethesda)* **7**, 2749–2762 (2017).

50. van Doorn, G. S. & Kirkpatrick, M. Transitions between male and female heterogamety caused by sex-antagonist selection. *Genetics* **186**, 629–645 (2010).

51. Sardell, J. M. et al. Sex differences in recombination in sticklebacks. *G3 (Bethesda)* **8**, 1971–1983 (2018).

52. Ma, W.-J. et al. Tissue specificity and dynamics of sex-biased gene expression in a common frog population with differentiated, yet homomorphic, sex chromosomes. *Genes* **9**, 294 (2018).

53. Rodrigues, N., Betto-Colliard, C., Jourdan-Pineau, H. & Perrin, N. Within-population polymorphism of sex-determination systems in the common frog (*Rana temporaria*). *J. Evol. Biol.* **26**, 1569–1577 (2013).

54. Burt, A., Bell, G. & Harvey, P. H. Sex differences in recombination. *J. Evol. Biol.* **4**, 259–277 (1991).
55. Lenormand, T. & Dutheil, J. Recombination difference between sexes: a role for haploid selection. *PLoS Biol* **3**, e63 (2005).

56. Lorch, P. D. Sex differences in recombination and mapping adaptations. *Genetica* **123**, 39–47 (2005).

57. Moen, T., Hoyheim, B., Munck, H. & Gomez-Raya, L. A linkage map of Atlantic salmon (*Salmo salar*) reveals an uncommonly large difference in recombination rate between the sexes. *Anim. Genet.* **35**, 81–92 (2004).

58. Coimbra, M. R. M. et al. A genetic linkage map of the Japanese flounder, *Paralichthys olivaceus*. *Aquaculture* **220**, 203–218 (2003).

59. Rodrigues, N., Studer, T., Dufresnes, C. & Perrin, N. Sex-chromosome recombination in common frogs brings water to the fountain-of-youth. *Mol. Biol. Evol.* **35**, 942–948 (2018).

60. Rodrigues, N., Merilá, J., Patrelle, C. & Perrin, N. Geographic variation in sex-chromosome differentiation in the common frog (*Rana temporaria*). *Mol. Ecol.* **23**, 3409–3418 (2014).

61. Rodrigues, N., Vuille, Y., Loman, J. & Perrin, N. Sex-chromosome differentiation and 'sex races' in the common frog (*Rana temporaria*). *Proc. Biol. Sci.* **282**, 20142726 (2015).

62. Brelsford, A. et al. Homologous sex chromosomes in three deeply divergent anuran species. *Evolution* **67**, 2434–2440 (2013).

63. O'Meally, D., Ezaz, T., Georges, A., Sarre, S. D. & Graves, J. A. M. Are some chromosomes particularly good at sex? Insights from amniotes. *Chromosome Res.* **20**, 7–19 (2012).

64. Herpin, A. & Schartl, M. Plasticity of gene-regulatory networks controlling sex determination: of masters, slaves, usual suspects, newcomers, and usurpators. *EMBO Rep.* **16**, 1260–1274 (2015).

65. Oshima, Y., Naruse, K., Nakamura, Y. & Nakamura, M. Sox3: a transcription factor for *Cyp19* expression in the frog *Rana rugosa*. *Gene* **445**, 38–48 (2009).

66. Foster, J. W. et al. Evolution of sex determination and the Y chromosome: *SRY*-related sequences in marsupials. *Nature* **359**, 531–533 (1992).

67. Smith, C. A. et al. The avian Z-linked gene *DMRT1* is required for male sex determination in the chicken. *Nature* **461**, 267–271 (2009).

68. Chen, S. et al. Whole-genome sequence of a flatfish provides insights into ZW sex chromosome evolution and adaptation to a benthic lifestyle. *Nat. Genet.* **46**, 253–260 (2014).

69. Ma, W.-J., Rodrigues, N., Sermier, R., Brelsford, A. & Perrin, N. *Dmrt1* polymorphism covaries with sex-determination patterns in *Rana temporaria*. *Ecol. Evol.* **6**, 5107–5117 (2016).

70. Brelsford, A., Dufresnes, C. & Perrin, N. Trans-species variation in *Dmrt1* is associated with sex determination in four European tree-frog species. *Evolution* **70**, 840–847 (2016).

71. Bewick, A. J., Anderson, D. W. & Evans, B. J. Evolution of the closely related, sex-related genes DM-W and DMRT1 in African clawed frogs (*Xenopus*). *Evolution* **65**, 698–712 (2011).

72. Yoshimoto, S. & Ito, M. A ZZ/WZ-type sex determination in *Xenopus laevis*. *FEBS J.* **278**, 1020–1026 (2011).

73. Capel, B. Vertebrate sex determination: evolutionary plasticity of a fundamental switch. *Nat. Rev. Genet.* **18**, 675–689 (2017).

74. Li, M.-H. et al. Antagonistic roles of *Dmrt1* and *Foxl2* in sex differentiation via estrogen production in tilapia as demonstrated by TALENs. *Endocrinology* **154**, 4814–4825 (2013).

75. Yang, Y.-J., Wang, Y., Li, Z., Zhou, L. & Gui, J.-F. Sequential, divergent, and cooperative requirements of *Foxl2a* and *Foxl2b* in ovary development and maintenance of zebrafish. *Genetics* **205**, 1551–1572 (2017).

76. Pauly, G. B., Hillis, D. M. & Cannatella, D. C. Taxonomic freedom and the role of official lists of species names. *Herpetologica* **65**, 115–128 (2009).

77. Haczkiewicz, K. & Ogielska, M. Gonadal sex differentiation in frogs: how testes become shorter than ovaries. *Zoolog. Sci.* **30**, 125–134 (2013).

78. Brelsford, A., Dufresnes, C. & Perrin, N. High-density sex-specific linkage maps of a European tree frog (*Hyla arborea*) identify the sex chromosome without information on offspring sex. *Heredity* **116**, 177–181 (2015).

79. Andrews, S. *FastQC: A Quality Control Tool For High Throughput Sequence Data*. Babraham Bioinformatics. http://www.bioinformatics.babraham.ac.uk/projects/fastqc (2010).

80. Catchen, J., Hohenlohe, P. A., Bassham, S., Amores, A. & Cresko, W. A. Stacks: an analysis tool set for population genomics. *Mol. Ecol.* **22**, 3124–3140 (2013).

81. Brelsford, A., Lavanchy, G., Sermier, R., Rausch, A. & Perrin, N. Identifying homomorphic sex chromosomes from wild-caught adults with limited genomic resources. *Mol. Ecol. Resour.* **17**, 752–759 (2016).

82. Karpinka, J. B. et al. Xenbase, the *Xenopus* model organism database; new virtualized system, data types and genomes. *Nucleic Acids Res.* **43**, D756–D763 (2015).

83. Pyron, R. A. Biogeographic analysis reveals ancient continental vicariance and recent oceanic dispersal in amphibians. *Syst. Biol.* **63**, 779–797 (2014).

84. Sun, Y.-B. et al. Whole-genome sequence of the Tibetan frog *Nanorana parkeri* and the comparative evolution of tetrapod genomes. *Proc. Natl Acad. Sci. USA* **112**, E1257–E1262 (2015).

85. Hotz, H., Uzzell, T. & Berger, L. Linkage groups of protein-coding genes in western palaearctic water frogs reveal extensive evolutionary conservation. *Genetics* **147**, 255–270 (1997).

86. Altschul, S. F., Gish, W., Miller, W., Myers, E. W. & Lipman, D. J. Basic local alignment search tool. *J. Mol. Biol.* **215**, 403–410 (1990).

87. Edgar, R. C. MUSCLE: multiple sequence alignment with high accuracy and high throughput. *Nucleic Acids Res.* **32**, 1792–1797 (2004).

88. Kearse, M. et al. Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. *Bioinformatics* **28**, 1647–1649 (2012).

89. Lanfear, R., Frandsen, P. B., Wright, A. M., Senfeld, T. & Calcott, B. PartitionFinder 2: New methods for selecting partitioned models of evolution for molecular and morphological phylogenetic analyses. *Mol. Biol. Evol.* **34**, 772–773 (2017).

90. Lanfear, R., Calcott, B., Ho, S. Y. W. & Guindon, S. Partitionfinder: combined selection of partitioning schemes and substitution models for phylogenetic analyses. *Mol. Biol. Evol.* **29**, 1695–1701 (2012).

91. Guindon, S. et al. New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. *Syst. Biol.* **59**, 307–321 (2010).

92. Bouckaert, R. et al. BEAST 2: a software platform for Bayesian evolutionary analysis. *PLoS Comput. Biol.* **10**, e1003537 (2014).

93. Drummond, A. J., Ho, S. Y. W., Phillips, M. J. & Rambaut, A. Relaxed phylogenetics and dating with confidence. *PLoS Biol.* **4**, e88 (2006).

94. Revell, L. J. phytools: an R package for phylogenetic comparative biology (and other things). *Methods Ecol. Evol.* **3**, 217–223 (2012).

95. Rastas, P. Lep-MAP3: robust linkage mapping even for low-coverage whole genome sequencing data. *Bioinformatics* **33**, 3726–3732 (2017).

## Acknowledgements

We thank E.A. Luna, J. Leuenberger, A. Athanasiades, J. Buser, A.D. King, A.H. McCall, S.J. Sferra, M. Biaggini, G. Sánchez-Montes and Dario Dominighetti for their help with collection and processing of samples. We also thank Philippe Geniez and the collection B. E.V. of the C.N.R.S. for the *P. saharicus* samples, Nicolas Salamin for advice on phylogenetic analyses and Barret Phillips for comments on the manuscript. This work was primarily funded by two grants (31003A_166323 & CRSII3_147625) awarded by the Swiss National Science Foundation to N.P. All computations were performed at the Vital-IT (http://www.vital-it.ch) Center for high-performance computing of the SIB Swiss Institute of Bioinformatics

## Author contributions

N.P. and A.B. conceived the project and were involved with study design alongside D.L.J. and G.L. D.L.J., G.L., W.J.M. and R.S. carried out the DNA extractions and RADseq library preparations. D.L.J. and G.L. performed the analyses. D.L.J., G.L. and N.P. wrote the manuscript. I.M., I.S., W.J.M., M.J.S., J.F., N.R., C.D., K.G., C.M.G., I.M.S., A.B., D.C., A.R., R.P.O., L.N.B., A.G.N., P.A.C. and S.Z. provided the samples and commented on the manuscript.

## Additional information

Supplementary Information accompanies this paper at https://doi.org/10.1038/s41467-018-06517-2.

**Competing interests:** The authors declare no competing interests.

**Reprints and permission** information is available online at http://npg.nature.com/ reprintsandpermissions/

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_6.png" /></figure>

**Open Access** This article is licensed under a Creative Commons

Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2018