M.D. Jérémie Bettoni
email: jeremiebettoni@gmail.com

Keywords: Osteoradionecrosis, Free flap, Revascularization surgery, Management, Complication, Morbidity

Reduction of morbidity of the revascularization surgery in the management of mandibular osteoradionecrosis by basilar edge preservation.

Introduction

The chronic complications of bone free flap revascularization surgery are mainly characterized by skin fistulization of the osteosynthesis material, which sometimes leads to reactivation of the osteoradionecrotic phenomenon. The objective of the study is to evaluate the benefit of mandibular basilar edge preservation in bone reconstructive surgery in irradiated areas performed for the treatment of advanced mandibular osteoradionecrosis.

Materials and methods

A retrospective monocentric study conducted between 2003 and 2018 including all patients undergoing revascularization surgery for the treatment of advanced osteoradionecrotic lesion with respect to the basilar margin was conducted.

Results

8 patients (7 males and 1 female, aged 50 to 63 years) who had a marginal mandibulectomy with reconstruction by bone free flap or composite free flap were included. The stability of the reconstruction (jonction native mandible/bone free flap) was achieved on average by using 1.75 [range 0-4] mini plates (Medartis® Modus 2.0, Medartis AG, Basel, Switzerland). During the follow-up (30 ± 13 months) no chronic complication related to a dissociation of the osteosynthesis material or a reactivation of the osteoradionecrotic phenomenon were identified.

Conclusions

Text

Introduction

First described by Régaud [START_REF] Régaud | Sur la sensibilité du tissu osseux normal vis-à-vis des radiations X et Y, et sur le mécanismes de l ost oradion crose[END_REF] in 1922 as "a vascular disorder following bone tissue irradiation", mandibular osteoradionecrosis (ORN) is a serious complication of head and neck oncological irradiation. Although the development of prophylactic measures and the advent of intensity modulated radiotherapy (IMRT) have reduced its incidence rate from 30% in 1968 to less than 10% today 2 , the management of osteoradionecrotic lesions remains controversial. In a review of literature published in 2015, Rice et al. [START_REF] Rice | The management of osteoradionecrosis of jaws A review[END_REF] propose to prioritize the therapeutic arsenal, ranging from long-term antibiotic therapy to radical surgery with free flap reconstruction. The results of decades of pathophysiological reflection, we gladly cite the theories of the "3H" of Marx [START_REF] Marx | Osteoradionecrosis : a nex concept of its pathophysiology[END_REF] ("Hypoxy, Hypocellularity, Hypovascularization"), the "2I" of Meyer [START_REF] Meyer | Infectious diseases of the jaws[END_REF] ("Irradiation, Infection"), the proposal of a vascular genesis by Bras et al. [START_REF] Bras | Osteoradionecrosis of the mandible : pathogenesis[END_REF] ("Radiation-induced obliteration of the inferior alveolar artery") or the radio-induced fibrosis of Delanian et al 7 ("Radiation-induced fibroatrophic process"). Rice et al. [START_REF] Rice | The management of osteoradionecrosis of jaws A review[END_REF] recommend that conservative management (hyperbaric oxygen therapy, long term antibiotic therapy, pentoclo protocol or sequestrectomy with direct mucosal closure) should be preferred and that radical management, consisting of resection of the pathological bone with reconstruction by microanastomosis flap, should be reserved to fracturial lesions or advanced subcanals refractory to medical treatment. D 'Hauthuille et al. 8 characterized this procedure as a "revascularization surgery", considering that the supply of healthy vascularized tissue constitutes a physiological response to the devascularization of osteoradionecrotic bones. Despite a high rate (86% and 100%) of effectiveness on the osteoradionecrotic phenomenon 9-15 , reconstructive surgery in irradiated areas has the reputation of complex surgery, supported by its high complication rate between 18% and 75%.

Among these complcations, 7% to 57.5% are secondary to disassociation of the osteosynthesis material 10,11,[START_REF] Baumann | Free flap reconstruction of osteoradionecrosis of mandible: a 10-year review and defect classification[END_REF][START_REF] Bettoni | Role of « revascularization surgery » in the management of mandibular osteoradionecrosis lesions : a report of 21 years of experience[END_REF] complicated by cutaneous or mucosal fistulization or reactivation of the osteoradionecrotic phenomenon.

Considering as a starting hypothesis that these chronic complications are related to a lack of primary stability of the reconstruction, leading to cortico-periosteal necrosis at the necks of the osteosynthesis screws resulting in a chronic inflammatory phenomenon on irradiated tissues. We decided to adapt our reconstruction technique.

The objective of this work is to evaluate the benefits of increasing the contact surfaces between the free-flap and the native bone, by preservation of the basilar edge during mandibular bone reconstruction in ORN context.

Material and method

From January 2003 to December 2018, we included all patients who received revascularization surgery for ORN lesions involving the mandibule below the inferior dental canal respecting the jaw basilar edge (with or without skin fistula) accessible to marginal mandibulectomy and requiring reconstruction by bone free flap to restore mandibular solidity. The feasibility of preserving jaw continuity was assessed preoparatively with radiographic exams (dental panoramic or X-Ray CT-scan); and intraoperatively by the absence of clinical osteolic involvement of the basilar edge.

Each osteoradionecrotic lesion was characterized by its topography (Urken classification [START_REF] Urken | Oromandibular reconstruction using microvascular composite free flaps. Report of 71 cases and new classification scheme for bony, soft-tissue, and neurologic defects[END_REF] ), the period between the end of irradiation and the revascularization surgery (months), as well as the modalities of reconstructive surgery (type of free flap, recipient vessels, modalities of microsurgical anastomoses). For each synthesis, the number and references of the material used were collected.

Revascularization surgery was considered effective in the absence of progressive continuation of the osteoradionecrotic phenomenon. Any local recurrence after a period of remission or osteoradionecrotic lesion at a distance from the first surgical site was considered as new episode of ORN.

During follow-up, acute complications occurring within the first 30 post-opretive days were differentiated from chronic complications occurring more than one month after surgery.

Ethical consideration

This study was performed in accordance with the rules of the local ethical committees.

Results

Between 2003 and 2018, 8 patients were treated with this modified revascularization surgery, consisting in marginal mandibulectomy with preservation of the basilar edge (Figure 1 and Supplementary material 1) and reconstruction by bone free flaps (7 men, 1 woman, aged between 50 and 63 years), for the management of advanced mandibular osteoradionecrotic lesions. The characteristics of this population are detailed in Table 1.

Concerning the topography of ORN, bone damage were located on the body of the mandible (8) and/or the mandibular angle (4) and/or the mandibular symphysis (2). In 3 patients, the resection of the pathological bone by piezosurgery® (Mectron S.p.a, Carasco, Italy) allowed a preservation of the lower alveolar nerve. Reconstruction for restoring the strength of the bone infrastructure involved 1 fibula free flap, 1 iliac crest free flap and 6 composite fibula free flaps with skin paddle (Table 2). Except in 1 patient (ID 3), where the strength of the construct was obtained solely by embedding the bone iliac crest free flap within the jaw resection site optimizes by suturing between the external oblique muscle and the native mucosa with absorbable suture (Ethicon Vycril® 3/0, Johnson & Johnson, New Brunswick, USA) (Figure 2), all osteosynthesis were performed on average by 2.6±1.9 miniplaques (Medartis® Modus 2.0, Medartis AG, Basel, Switzerland) (Figure 3). Excluding the miniplates necessary for osteosynthesis of segmental osteotomies of the bone free flap, 1.75 (range 0-4) miniplaques were used on average to stabilize the assembly between native mandibular bone and bone flap.

The mean operating time for revascularization surgery was 6h56±1h17 min.

Regarding the microsurgery, all the anatosmoses were performed in a terminoterminal fashion.

During the clinical follow-up (30 ± 13 months), we reported a post-operative cervical hematoma for 1 patient (ID 7) on post-operative day 1, requiring a revision of the anastomosis and allowing the flap to be saved, and one local infection secondary to disunion for another patient (ID 1). Nevertheless, there were no chronic complication related to the disassociation of osteosynthesis material (Table 3).

Discussion

In this study, we aimed to evaluate the benefits of preserving the basilar edge during mandibular bone reconstruction in an ORN context. We reported the management of height patient with an overall acute complication rate of 11.2%, and no chronic complication (dislocation of osteosynthesis material, or new episode of the osteoradionecrotic phenomenon observed).

Our study showed a lower complication rate than the others studies dealing with revascularization surgery in the treatment of mandibular ORN 9-15 . While the reduction in the rate of acute complication may be explained by early management of ORN lesions before a fracture stage, it seems that basilar edge concervation allows avoiding the occurrence of chronic complications related to osteosynthesis material (disassociation, fistulization or reactivation of the osteoradionecrotic process). Indeed, inspired by the tenon-mortise assembly, the significant increasing of the contact surface between the bone flap and the native bone seems to provide better primary stability for the assembly. Furthermore, the mechanical stress and tension activate osteogenesis in vitro and in vivo through inflammatory chemokines [START_REF] Burger | Mechanical stress and osteogenesis in vitro[END_REF][START_REF] Sato | Mechanical tension-stress induces expression of bone morphogenetic protein (BMP)-2 and BMP-4, but not BMP-6, BMP-7, and GDF-5 mRNA, during distraction osteogenesis[END_REF][START_REF] Song | Mechanical Stress Regulates Osteogenesis and Adipogenesis of rat Mesenchymal Stem Cells through PI3K/Akt/GSK-3 / -Catenin Signaling Pathway[END_REF] . We think that this procedure reduces the mechanical stresses, inflammatory processes at the neck of the osteosynthesis screws and the micromovements at the origin of corticostero-periosteal necrosis, leading to a disassociation of the material and chronic inflammation. Thus, skin fistulization and bone exposure, leading to the reactivation of the osteoradionecrotic processes, are avoided.

By reducing the number of miniplates required for robust assembly, this procedure reduces the use of inert materials in an irradiated area. In addition to the economic aspects, from a technical point of view the preservation of the basilar edge makes it possible to perform mandibulectomies with preservation of the alveolar nerve diverted from its canal and spatialized into the mandibular outer table. In addition, the preservation of the native bone rail makes it possible to guide spatialization, the modalities of segmental osteotomies and to prevent the risk of condyle malposition leading to mandibular laterodeviation or limitation of the buccal opening. Finally, from an aesthetic point of view, the preservation of the basilar edge respects the oval of the face and limits subsequent asymmetries.

With only 8 cases reported in 15 years, our experience highlights the problem of early diagnosis and management of osteoradionecrotic lesions. Although it is known that the more complex the reconstruction, the higher the complication rate [START_REF] Kroll | Choice of flap and incidence of free flap success[END_REF] , it is still difficult to propose free flap reconstruction in irradiated areas given the risks of acute and chronic complications (Table 3). It therefore seems essential to consider mandibular ORN as a real public health problem and to redefine an effective and documented therapeutic algorithm. Because recently the questioning of the effectiveness of the main conservative treatments by Annane et al. [START_REF] Annane | Hyperbaric oxygen therapy for radionecrosis of the jaw : a randomized, placebocontrolled, double-blind trial from the ORN96 study group[END_REF] , Madrid et al [START_REF] Madrid | Osteoradionecrosis : An update[END_REF] or Martos-Fernandes et al [START_REF] Martos-Fernandez | Pentoxifylline, topopherol, and clodronate for the treatment of mandibular osteoradionecrosis: a systematic review[END_REF] , requires a profound reflection on the therapeutic strategy to be adopted concerning mandibular osteoradionecrosis. It is therefore paradoxical that having proven its effectiveness (range 86% -100%) 9-15 , revascularization surgery is not only considered as a treatment of recourse. Some authors [START_REF] De Montreuil | Limited mandibular osteoradionecrosis. Early conservative surgical treatment[END_REF][START_REF] Hauthuille | Part III: free periosteal flaps as treatment for mandibular osteoradionecrosis[END_REF] propose to extend its indications immediately after the failure of conservative treatments. Considering the refractoriness of medical treatments as a gravity criterion in its own right, requires a modification of the therapeutic strategy. In 2009, D'Hauthuille et al. [START_REF] Hauthuille | Part III: free periosteal flaps as treatment for mandibular osteoradionecrosis[END_REF] proposed the use of the periosteal flap in the event of refractory superficial damage during sequestrectomy and direct closure. This attitude based on the osteoinductive [START_REF] Takato | Vascularized periosteal grafts: an experimental study using two different forms of tibial periosteum in rabbits[END_REF] and neoangiogenic 27 potential of periosteal flaps makes it possible to offer effective early treatment of osteoradionecrotic phenomena. But above all, this therapeutic sequence is part of a revolution that reserves effective therapy only in the event of an outdated or fractured lesion where the complexity of the reconstruction exposes the patient to a high complication rate [START_REF] Kroll | Choice of flap and incidence of free flap success[END_REF] . If the objective of mandibular ORN management is to stop the osteoradionecrotic phenomenon by controlling complications, it seems essential to identify ORNs refractory to concerted treatment and to propose early revascularization surgery in a multicentric study to provide a sufficient level of evidence to load the current paradigm.

The basilar margin preservation offers the prospect of a reduction in chronic complications and aesthetic sequelae associated with bone reconstruction.

Nevertheless, this attitude requires close collaboration between the various actors (radiotherapists, dentists, general practitioners and surgeons) allowing early diagnosis and effective and appropriate graduated management. Moreover, given its technical constraints, mandibular ORN surgery should be centralized in expert centers because any failure penalizes the success of subsequent interventions exposing to real therapeutic impasses.

Conclusion

Given its complexity, revascularization surgery must be conceptualized as a unique and curative procedure that restores or preserves the solidity of the mandibular infrastructure and bone sealing essential for oral health. Our study is in favour of preserving the basilar edge where it is possible. This procedure seems to offer a better primary stability during assembly, allowing a reduction in both stress on the osteosynthesis equipment and the resulting chronic complications (fistula and/or reactivation of the osteoradionecrotic phenomenon). But this preliminary results have to be confirmed with more cases or by studies with higher level of evidence.

7 Delanian S, Lefaix JL. The radiation-induced fibroatrophic process : therapeutic perspective via the antioxidant pathway. Supplementary material 1: X-ray CT scan 3D reconstruction at 1 year postoperative (ID 5).

Captions to illustrations

Figure 1 :Figure 2 :

 12 Figure 1: (a.)Peroperative photography of a marginal piezosurgery® mandibulectomy

Figure 3 :

 3 Figure 3: 3D reconstruction of a one week post-operative CT scan of a

Table 1 :

 1 Characteristics of the study population.

Table 2 :

 2 Modalities for reconstructing bone loss according to each patient included *Urken classification : S-Simphysis / B-Body / R-Angle **this difference is due to the neef of two type of osteosynthesis miniplates, some are used to fix the segmental osteotomies of the bone free flap whereas the others were used to stabilize the assembly between native mandibular bone and bone flap.

Table 3 :

 3 Literature review of acute and chronic complications in series of patients managed for mandibular osteoradionecrosis. (nk: not know)

27-Ortak T, Ozdemir R, Uysal A, Ulusoy MG, Sungur N, Sahin B et al. Osteogenic capacities of periost grafts, periost flaps and prefabricated periosteal flaps : experimental study. J Craniofac Surg 2005;16:594-600 (PMID: 16077303).

Osteoradionecrotic lesions involving at least 2 Urken anatomical units 7

Skin fistula 2


Romain Feron
email: romainferon91@gmail.com

Qiaowei Pan

Ming Wen

Boudjema Imarazene

Elodie Jouanno

Jennifer Anderson

Amaury Herpin

Laurent Journot

Hugues Parrinello

Christophe C Klopp

Verena A Kottler

Alvaro S Roco

Kang Du

Susanne Kneitz

Mateus Adolfi

Catherine A Wilson

Braedan Mccluskey

Angel Amores

Thomas Desvignes

Frederick W Goetz

Ato Takanashi

Mari Kawaguchi

H William

Detrich Iii

Marcos Oliveira

Rafael Nobrega

Takashi Sakamoto

Masatoshi Nakamoto

Anna Wargelius

Ørjan Karlsen

Zhongwei Wang

Matthias Stöck

Robert M Waterhouse

Ingo Braasch

John H Postlethwait

Manfred Schartl

Yann Guiguen
email: yann.guiguen@inrae.fr

Affiliations

RADSex: a computational workflow to study sex determination using Restriction Site-Associated DNA Sequencing data

Keywords: sex determination, RAD-Sequencing, fish, computational workflow, visualization

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.

INTRODUCTION

Sexual reproduction is widespread in animals [START_REF] Goodenough | Origins of eukaryotic sexual reproduction[END_REF] and, in gonochoristic species, involves the segregation of the male and female gonadal functions into separate individuals for their entire lives [START_REF] Bachtrog | Sex Determination: Why So Many Ways of Doing It[END_REF] . Males produce small, mobile gametes while females produce large, immobile gametes, and this difference has fueled the divergent evolution of morphological, physiological, and behavioral traits between the two sexes. Individuals acquire sex-specific traits during sexual development, which starts with sex determination, the process that controls whether an individual develops male or female reproductive organs. Sex determination can be triggered by genetic factors (Genetic Sex Determination, GSD), like in virtually all mammals and birds [START_REF] Bachtrog | Sex Determination: Why So Many Ways of Doing It[END_REF] , by environmental factors (Environmental Sex Determination, ESD), like temperature in many non-avian reptilia [START_REF] Pezaro | The ecology and evolution of temperature-dependent reaction norms for sex determination in reptiles: a mechanistic conceptual model[END_REF] , or a combination of the both, as in several fish species [START_REF] Piferrer | Genetic, endocrine, and environmental components of sex determination and differentiation in the European sea bass (Dicentrarchus labrax L.)[END_REF] . When genetic factors are involved, sex determination is entirely or partially controlled by one or multiple master sex determining (MSD) genes located on sex chromosomes. Genetic sex determination systems can involve male heterogamy (XX/XY and XX/X0), female heterogamy (ZZ/ZW and ZZ/Z0), or multiple loci on different chromosomes (polygenic systems) [START_REF] Bachtrog | Sex Determination: Why So Many Ways of Doing It[END_REF] . In male-and female-heterogametic systems, the sex chromosomes can be morphologically different, i.e. heteromorphic, or similar, i.e. homomorphic [START_REF] Bachtrog | Sex Determination: Why So Many Ways of Doing It[END_REF] . Initially, the overwhelming majority of knowledge on the structure and evolution of sex chromosomes came from studies in mammals [START_REF] Wallis | Sex determination in mammals--before and after the evolution of SRY[END_REF] and in Drosophila [START_REF] Salz | Sex determination in Drosophila: The view from the top[END_REF] , where suppression of recombination between the X and Y chromosomes led to degeneration of the Y.

These findings have spawned both theoretical and empirical interest in sex chromosomes as models to study the consequences of recombination suppression and associated processes [START_REF] Bachtrog | The temporal dynamics of processes underlying Y chromosome degeneration[END_REF][START_REF] Charlesworth | The degeneration of Y chromosomes[END_REF][START_REF] Corcoran | Introgression maintains the genetic integrity of the mating-type determining chromosome of the fungus Neurospora tetrasperma[END_REF][START_REF] Doorn | Turnover of sex chromosomes induced by sexual conflict[END_REF][START_REF] Gu | Conserved Patterns of Sex Chromosome Dosage Compensation in the Lepidoptera (WZ/ZZ): Insights from a Moth Neo-Z Chromosome[END_REF][START_REF] Huylmans | Global Dosage Compensation Is Ubiquitous in Lepidoptera, but Counteracted by the Masculinization of the Z Chromosome[END_REF][START_REF] Muyle | Rapid De Novo Evolution of X Chromosome Dosage Compensation in Silene latifolia, a Plant with Young Sex Chromosomes[END_REF][START_REF] Peichel | The master sex-determination locus in threespine sticklebacks is on a nascent Y chromosome[END_REF] . More recent studies investigating sex chromosomes in insects [START_REF] Blackmon | Sex Determination, Sex Chromosomes, and Karyotype Evolution in Insects[END_REF] , non-avian reptilia [START_REF] Modi | Sex chromosomes and sex determination in reptiles[END_REF] , amphibians [START_REF] Miura | Sex Determination and Sex Chromosomes in Amphibia[END_REF] , and fishes [START_REF] Kikuchi | Novel sex-determining genes in fish and sex chromosome evolution[END_REF] have found many homomorphic sex chromosomes displaying varying levels of differentiation, with the extreme case of a sex locus restricted to allelic variation of a single nucleotide as reported in the Japanese pufferfish [START_REF] Kamiya | A trans-species missense SNP in Amhr2 is associated with sex determination in the tiger pufferfish, Takifugu rubripes (fugu)[END_REF] . These results called into question the single unified concept of sex chromosome evolution and highlighted the importance of obtaining a broader understanding of sex determination and sex chromosomes in many species across the tree of life. The first step in this process is identifying sex-specific genomic sequences, i.e. sequences from the sex locus that are only found in one of the two sexes because of allelic divergence between the sex chromosomes or because of a large insertion in the hemizygous chromosome. These sex-specific sequences can then be aligned to a reference genome to locate the sex locus, identify candidate MSD genes and other genes involved in sex determination, and characterize patterns of differentiation between the sex chromosomes. In addition, facilitating the identification of such sequences has practical applications: sex is an important factor in ecological [START_REF] Benestan | Sex matters in massive parallel sequencing: Evidence for biases in genetic parameter estimation and investigation of sex determination systems[END_REF] and conservation studies [START_REF] Ancona | Estimating adult sex ratios in nature[END_REF] , as well as agriculture (Al-Ameri, Al-Qurainy, Gaafar, [START_REF] Al-Ameri | Molecular Identification of Sex in Phoenix dactylifera Using Inter Simple Sequence Repeat Markers[END_REF][START_REF] Liao | Development of male-specific markers and identification of sex reversal mutants in papaya[END_REF][START_REF] Spigler | Genetic mapping of sex determination in a wild strawberry, Fragaria virginiana, reveals earliest form of sex chromosome[END_REF] and animal production [START_REF] Dan | Genetic Differentiation and Efficient Sex-specific Marker Development of a Pair of Y-and X-linked Markers in Yellow Catfish[END_REF][START_REF] Yano | The sexually dimorphic on the Y-chromosome gene (sdY) is a conserved male-specific Y-chromosome sequence in many salmonids[END_REF]) , yet the sex of an individual cannot always be easily determined by non-invasive methods.

Until recently, such analyses were not readily feasible mainly due to technical barriers preventing the precise and high-throughput identification of sex-specific genomic sequences, especially in non-model species, but advances in sequencing technologies have enabled the study of genetic sex determination in a much broader spectrum of taxa by comparing genomes from phenotypically distinct males and females. In essence, this process is akin to comparing genomic differences between two populations or between a mutant and a wild type genotype, and therefore methods from molecular genetics, population genetics, and ecology can be applied. A popular representational approach to comparing the genetics of populations is Restriction Site-Associated DNA Sequencing, or RAD-Seq [START_REF] Andrews | Harnessing the power of RADseq for ecological and evolutionary genomics[END_REF] . RAD-Seq generates short sequences for a small but consistent fraction of the genome and thus allows the sequencing and comparison of multiple individuals from several populations at relatively low cost [START_REF] Davey | Genome-wide genetic marker discovery and genotyping using next-generation sequencing[END_REF] and without requiring additional genomic resources. RAD-Seq has been successfully used to identify sex-specific sequences in non-model species from diverse taxa, including fish [START_REF] Drinan | Identification of Genomic Regions Associated With Sex in Pacific Halibut[END_REF] , amphibians [START_REF] Bewick | A large pseudoautosomal region on the sex chromosomes of the frog Silurana tropicalis[END_REF] , non-avian reptilia [START_REF] Gamble | Using RAD-seq to recognize sex-specific markers and sex chromosome systems[END_REF][START_REF] Gamble | The Discovery of XY Sex Chromosomes in a Boa and Python[END_REF][START_REF] Gamble | Restriction Site-Associated DNA Sequencing (RAD-seq) Reveals an Extraordinary Number of Transitions among Gecko Sex-Determining Systems[END_REF][START_REF] Gamble | XX/XY Sex Chromosomes in the South American Dwarf Gecko (Gonatodes humeralis)[END_REF][START_REF] Gamble | Identification of sex-specific molecular markers using restriction site-associated DNA sequencing[END_REF][START_REF] Nielsen | Dynamic sex chromosomes in Old World chameleons (Squamata: Chamaeleonidae)[END_REF][START_REF] Nielsen | ZZ/ZW Sex Chromosomes in the Endemic Puerto Rican Leaf-Toed Gecko (Phyllodactylus wirshingi)[END_REF] , invertebrates [START_REF] Carmichael | Identification of a sex-linked SNP marker in the salmon louse (Lepeophtheirus salmonis) using RAD sequencing[END_REF][START_REF] Mathers | Transition in sexual system and sex chromosome evolution in the tadpole shrimp Triops cancriformis[END_REF][START_REF] Pratlong | Evidence for a genetic sex determination in Cnidaria, the Mediterranean red coral (Corallium rubrum)[END_REF] , and plants [START_REF] Kafkas | Identification of sex-linked SNP markers using RAD sequencing suggests ZW/ZZ sex determination in Pistacia vera L[END_REF] , and to identify the sex chromosomes and the sex locus in some of these species [START_REF] Wilson | Wild sex in zebrafish: loss of the natural sex determinant in domesticated strains[END_REF] . Although dedicated pipelines have been developed to analyze RAD-Seq data specifically for sex determination [START_REF] Gamble | Identification of sex-specific molecular markers using restriction site-associated DNA sequencing[END_REF] , most of the above mentioned studies have used STACKS (J. [START_REF] Catchen | Stacks: an analysis tool set for population genomics[END_REF][START_REF] Catchen | Stacks: building and genotyping Loci de novo from short-read sequences[END_REF][START_REF] Hohenlohe | Population genomic analysis of model and nonmodel organisms using sequenced RAD tags[END_REF][START_REF] Rochette | Deriving genotypes from RAD-seq short-read data using Stacks[END_REF] to cluster RAD-Seq reads into polymorphic markers subsequently filtered with custom in-house scripts. STACKS, however, requires substantial computational resources to run on large datasets and depends on multiple parameters that can greatly influence the outcome of the analysis [START_REF] Paris | Lost in parameter space: a road map for stacks[END_REF][START_REF] Rodríguez-Ezpeleta | Population structure of Atlantic mackerel inferred from RAD-seq-derived SNP markers: effects of sequence clustering parameters and hierarchical SNP selection[END_REF][START_REF] Shafer | Bioinformatic processing of RAD-seq data dramatically impacts downstream population genetic inference[END_REF] , for instance, in the context of sex determination, the classification of a marker as a sex-specific sequence ( e.g. from a sex-specific insertion) or as a sex-specific allele ( e.g. an allele specific to one sex in a polymorphic locus) [START_REF] Utsunomia | Evolution and conservation of Characidium sex chromosomes[END_REF] . In addition, the lack of visualization tools for RAD-Seq data makes the interpretation of results from some datasets difficult, for instance, datasets containing individuals with mis-assigned sex, strong population structure, or sex-bias in sequencing depth.

To overcome these limitations, we developed RADSex, a RAD-Seq data analysis computational workflow specifically designed to study sex determination. RADSex is simple to use, makes no assumptions whether sex-biased markers are sex-biased sequences or sex-biased alleles in a polymorphic locus, requires few resources, and offers visualization tools in the form of an R package. To demonstrate the relevance of RADSex, we analyzed a previously published dataset from the Japanese medaka, Oryzias latipes , for which we replicated previous findings, including the identification of two sex-reversed individuals, and we uncovered a previously overlooked Y-specific polymorphism in the sequenced population. We then used RADSex to analyze new datasets that we generated from 15 fish species spanning multiple systematic orders and compared the runtime and memory usage to that of STACKS. In six of these species, we identified the sex-determination system as well as multiple sex markers, and we identified the sex chromosome in one species for which a reference genome was available. Our results show that RADSex is well-suited to study genetic sex-determination in non-model species, outperforming STACKS in speed, resource usage, and ease of use. Furthermore, our multi-species analysis provides insights into the genetic mechanisms of sex-determination in fishes and highlights some of the limitations of RAD-Seq in detecting small sex loci.

MATERIAL AND METHODS

Overview of the RADSex analysis workflow

The underlying principle of RADSex is to group identical RAD-Seq reads from all individuals in a dataset into non-polymorphic markers, and then consider the presence or absence of each marker in each individual. Markers created by RADSex thus differ from that of other software such as Stacks (J. M. [START_REF] Catchen | Stacks: building and genotyping Loci de novo from short-read sequences[END_REF] which typically attempt to reconstruct genotypes from polymorphic markers. The RADSex workflow ( Fig. 1 The next step is to identify markers significantly associated with sex. The distrib command uses the table of marker depths generated with process to compute the distribution of markers between males and females. Phenotypic sex for each individual is given by a user-supplied tabulated file with the identifier of each individual in the first column and its sex in the second column. In all radsex commands, a marker is defined as present in a given individual if the depth of this marker in this individual is higher than a user-specified threshold given by the parameter --min-depth (or -d) ( Supp. Fig. 1 ). Following this definition, the distrib command enumerates marker presence in every possible combination of number of males and number of females ( e.g. for a population of one male and two females: 0♂1♀; 0♂2♀; 1♂0♀; 1♂1♀; 1♂2♀), and calculates the probability of association with sex for each combination using Pearson's chi-squared test of independence with Yates' correction for continuity. The associated p-value is obtained from the Cumulative Density Function of the chi-squared distribution with one degree of freedom, as implemented in samtools [START_REF] Li | The Sequence Alignment/Map format and SAMtools[END_REF] When a genome assembly is available, markers can be aligned to it using the radsex map command in order to locate genomic regions differentiated between sexes. Markers obtained from either the process , signif , or subset commands are aligned to the genome using the BWA mem library [START_REF] Li | Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM[END_REF] , and only markers uniquely aligned with a mapping quality higher than a user-specified threshold are retained. Two metrics are computed for each retained marker: 1) the probability of association with sex p association with sex , estimated from a chi-squared test with Bonferroni correction as described for the distrib command, and 2) the sex-bias ( S), defined as S = ♂ / ♂ total -♀ / ♀ total , where M and F are the number of males and females in which the marker is present, and M total and F total are the total number of males and females in the population, inferred from the user-supplied sex information file. Sex-bias thus ranges from -1 for a marker present in all females and absent from all males to +1 for a marker present in all males and absent from all females, and is zero for a marker present in the same proportion of males and females. Using the output of radsex map ( Supplementary table 3 ), the radsex_map_circos() function of sgtr generates a circular plot in which sex-bias, on the top track, and -log 10 (p association with sex ), on the bottom track, are plotted against genomic position using R's circlize package (Z. [START_REF] Gu | circlize implements and enhances circular visualization in R[END_REF] . Each sector of the plot represents a linkage group, with unplaced scaffolds, which are important to include due to frequent assembly problems with sex chromosomes, appearing in the last sector. Each metric can also be displayed in a Manhattan plot using the radsex_map_manhattan() function; in addition, a linear plot showing both metrics for a given genomic region can be generated using the radsex_map_region() function from sgtr .

Finally, radsex provides two commands to assess the distribution of marker depths in the entire dataset. First, the freq command computes the number of individuals in which each marker is present for a given minimum depth value and generates a count table of markers present in each number of individuals. Second, the depth command computes the minimum, maximum, median, and average marker depth for each individual using markers present in more than 75% of individuals (or a user-specified threshold). Results from radsex freq and radsex depth can be visualized with the radsex_freq() and radsex_depth() functions from sgtr .

Sample collection

General information on the different species, species collectors, and samples is given in Supp. Atlantic cod ( Gadus morhua ): Adult Atlantic cods were sampled from an aquaculture population raised for a nutrition trial [START_REF] Karlsen | Copepods enhance nutritional status, growth and development in Atlantic cod (Gadus morhua L.) larvae -can we identify the underlying factors[END_REF] . They were sampled at 2.5 years-old in January Walleye ( Sander vitreus ): Fin clips were obtained from individuals from the Spirit Lake strain collected by gillnet in April 2015 from Clear Lake, Iowa (USA). These fish were part of an annual broodstock collection to obtain gametes for juvenile production used in stocking inland lakes.

The sex of the fish was determined by a combination of the presence of expelled semen or ovulated eggs, dissection for fish that were dead, and size for a small number of fish. EDTA pH 8, 0.2% SDS, 0.2 M NaCl) containing 200 μg/ml proteinase K for three hours at 80°C, and 500 μL phenol was added to each sample. After mixing and incubating at room temperature for 10 min, 500 μ L chloroform/isoamyl alcohol (24:1) was added. Samples were incubated with occasional inversion at room temperature for 10 min and then centrifuged for 10 min at 5000 g (5°C). The upper layer was transferred to a new tube and one volume of chloroform/isoamyl alcohol (24:1) was added. The samples were kept at room temperature for 10 min with occasional inverting, and subsequently centrifuged for 10 min at 5000 g (5°C). The supernatant was transferred to a glass vial on ice and 2.5 times the sample volume of cold 100% ethanol was carefully added. The DNA was spooled with a glass rod and dissolved in TE buffer pH 8.

Marbled notothen (

All gDNA concentrations were quantified with a NanoDrop ND2000 spectrophotometer (Thermo scientific, Wilmington, Delaware) and/or Qubit (ThermoFisher, France) before processing gDNA for RAD libraries construction.

RAD-Sequencing and demultiplexing of RAD-Seq reads RAD libraries were constructed for each species from individual fish gDNA using a single restriction enzyme, SbfI, following the standard protocol as previously described [START_REF] Amores | Genome evolution and meiotic maps by massively parallel DNA sequencing: spotted gar, an outgroup for the teleost genome duplication[END_REF] . Libraries were sequenced as single end 100 bp reads on one lane of Illumina HiSeq 2500. Quality control and demultiplexing of the reads was performed with the process_radtags.pl wrapper script from Stacks version 1.44 (J. M. [START_REF] Catchen | Stacks: building and genotyping Loci de novo from short-read sequences[END_REF] using the following options: discard reads with low quality scores (-q), remove any reads with an uncalled base (-c), and rescue barcodes and RAD-seq reads (-r).

Demultiplexing results for each species are summarized in Supplementary Table 7 .

Performance measurements

Runtime and peak memory usage of radsex (version 1.0.0) and Stacks (version 2.2)

were measured with "/usr/bin/time -v <command>" on the Genotoul computational platform (Toulouse, France) using 16 threads. For radsex , performance was measured for the process , distrib , and signif commands run sequentially with a minimum depth of 1, as this value would require the longest runtime and memory. For Stacks, performance was measured for the denovo_map.pl pipeline with the following settings: "-M 3 -n 2 -X populations:--fstats", which represent a standard usage to identify sex-specific markers.

https://github.com/SexGenomicsToolkit/sgtr for sgtr . Both radsex and sgtr are also available on the conda channel Bioconda. A complete documentation for RADSex is available at https://sexgenomicstoolkit.github.io/html/radsex/introduction.html .

RESULTS

Performance comparison of RADSex and Stacks

The runtime and peak memory usage of radsex process, distrib, and signif were measured on datasets from 15 fish species that we generated for this study and compared to that of STACKS denovo_map.pl (version 2.2, (J. [START_REF] Catchen | Stacks: an analysis tool set for population genomics[END_REF] ). Using our configuration,

RADSex validation on a published dataset

To assess whether RADSex can recover previous results on sex determination, we analyzed a dataset which has been used to characterize a major sex-determination locus with a STACKS-based approach [START_REF] Wilson | Wild sex in zebrafish: loss of the natural sex determinant in domesticated strains[END_REF] . This dataset consists of RAD-Seq reads from 31 males and 30 females from the carbio strain of Oryzias latipes , the Japanese medaka, which has a known XX/XY sex determination system with chromosome 1 as the sex chromosome [START_REF] Matsuda | DMY is a Y-specific DM-domain gene required for male development in the medaka fish[END_REF][START_REF] Nanda | A duplicated copy of DMRT1 in the sex-determining region of the Y chromosome of the medaka, Oryzias latipes[END_REF] . We analyzed this O. latipes dataset with RADSex using a minimum depth of 10 (-d 10) to define the presence of a marker in an individual. This value of minimum depth was chosen based on the median sequencing depth in each individual computed with radsex depth : the average median sequencing depth was 34, and thus a stringent minimum depth of 10 can be used to discard markers from potential sequencing errors but still retain real markers. In total, we found 121,492 markers present in at least one individual with a minimum depth of 10, among which 194 markers were significantly associated with male phenotype ( radsex signif , p < 0.05, chi-squared test with Bonferroni correction, highlighted tiles in Fig 4 .A ). Among these 194 markers, 165 were present in at least 20 of the 31 males. In addition, a single marker present in 26 of the 30 females and 5 males was significantly associated with female phenotype ( Fig. 4.A ). However, we did not find any marker simultaneously present in more than 29 of the 31 phenotypic males and absent from all phenotypic females. To understand why no marker was found exclusively in all males given that medaka has an XY/XX sex determination system, we extracted individual depth for all 165 markers present in at least 20 males and none of the females using radsex subset . We showing zero or very low depth for the extracted markers. These two outlier phenotypic males had already been found to be sex-reversed genetic XX females by retrospective genotyping in [START_REF] Wilson | Wild sex in zebrafish: loss of the natural sex determinant in domesticated strains[END_REF] .

After correcting the designated phenotypic sex of these two individuals, radsex signif yielded 232 markers significantly associated with male phenotype and no marker significantly associated with female phenotype (p < 0.05, chi-squared test with Bonferroni correction), as expected for an XY/XX system. In addition, six males forming a sub-cluster within the male cluster from sgtr radsex_markers_depth() had low or null depth for almost half of the extracted markers, suggesting that several Y chromosome haplotypes could be present in the sequenced population.

To identify the sex chromosome and to delimit the sex-differentiated region, we aligned the RADSex markers to the assembly used by Overall, these results are consistent with Wilson et al. [START_REF] Wilson | Wild sex in zebrafish: loss of the natural sex determinant in domesticated strains[END_REF] , who found 248 RAD-seq reads associated with male phenotype mapping to a region ranging from 14.3 to 32.5 Mb on Ola1. However, our analysis revealed a previously unidentified Y-specific polymorphism in the sequenced population. This example highlights the effectiveness and versatility of RADSex and its visualization tools to identify, explore, and explain patterns in datasets that can be overlooked with traditional approaches.

Analysis of 15 new ray-finned fish RAD-sequencing datasets RAD-Seq datasets were generated for 15 ray-finned species in which no markers for genetic sex determination were available when we initiated this work. For each species, we used RADSex to search for sex-biased markers with a minimum marker depth (-d) of 1, 2, 5, and 10. We report results for a minimum marker depth of 10 ( Table 1 ) as it is the most stringent, except in three species with the lowest sequencing depths, for which we report results for a minimum marker depth of 2 ( Supp Fig. 2 ). We identified an XX/XY sex determining system in six species ( Cyprinus carpio , Gymnotus carapo , Plecoglossus atlivelis , Tinca tinca , Gadus morhua , and Poecilia sphenops ) with a few male and/or female outliers in four of these species ( Fig. 5 ). Markers significantly associated with phenotypic sex for each species are provided in Supplementary table 6 . An XX/XY SD system is in agreement with previous reports for C. carpio [START_REF] Gomelsky | Chromosome set manipulation and sex control in common carp: a review[END_REF] , G. carapo (da Silva, Matoso, Artoni, & Feldberg, 2014) , P. altivelis [START_REF] Watanabe | Detection of ayu sex-linked DNA markers using homologous clones[END_REF] , and G. morhua [START_REF] Haugen | Sex differentiation in Atlantic cod (Gadus morhua L.): morphological and gene expression studies[END_REF][START_REF] Whitehead | Ovarian development and sex ratio of gynogenetic Atlantic cod (Gadus morhua)[END_REF]) , but conflicts with the identification of P.

sphenops being female heterogametic (I. [START_REF] Nanda | Primitive sex chromosomes in poeciliid fishes harbor simple repetitive DNA sequences[END_REF] .

However, opposite population-specific sex determination systems have also been described in P. sphenops [START_REF] Volff | Variability of genetic sex determination in poeciliid fishes[END_REF] , and because previous studies were performed on ornamental fish and laboratory strains, previously derived conclusions may not apply to wild populations like the one we analyzed. Lastly, to the best of our knowledge, sex determination was not previously characterized in the tench, Tinca tinca .

In G. morhua , three markers were significantly associated with male phenotype, but they were also found in two females and none of them was present in more than 27 of the 34 males (Fig. 5.B) . Furthermore, multiple markers were found predominantly in males but not significantly associated with sex, with eight markers consistently present with depth higher than five in 20 males and absent from all females (these markers are highlighted with arrows in Fig.

6.A

). The three markers significantly associated with sex (red-outlined boxes) and most other markers strongly but not significantly associated with sex aligned to cod linkage group 11 (Gmo11) ( Fig. 6.B and 6.C ). Markers strongly associated with sex aligned to a region spanning from 0 to 20 Mb on Gmo11, and the three markers significantly associated with sex aligned between 10 to 15 Mb on this chromosome, in a region that was previously characterized as the region containing the sex locus in this species [START_REF] Kirubakaran | Characterization of a male specific region containing a candidate sex determining gene in Atlantic cod[END_REF][START_REF] Star | Genomic characterization of the Atlantic cod sex-locus[END_REF] .

These results confirm the previously identified XX/XY sex determining system and Gmo11 as the sex chromosome in Atlantic cod [START_REF] Haugen | Sex differentiation in Atlantic cod (Gadus morhua L.): morphological and gene expression studies[END_REF][START_REF] Whitehead | Ovarian development and sex ratio of gynogenetic Atlantic cod (Gadus morhua)[END_REF] , but our findings suggest an additional sex determination complexity in the studied Atlantic cod aquaculture population which had not been observed before on wild populations of Atlantic cod [START_REF] Kirubakaran | Characterization of a male specific region containing a candidate sex determining gene in Atlantic cod[END_REF][START_REF] Star | Genomic characterization of the Atlantic cod sex-locus[END_REF] . The incomplete sex-linkage of sex-biased markers may be explained by a Y chromosome polymorphism and/or the existence of sex-reversed genetic XX females. Although Y chromosome population differences have rarely been explored, at least in fish species, the results from our re-analysis of a medaka population suggest that such polymorphism may be more frequent than expected. In addition, female-to-male sex-reversal could be a consequence of intensive rearing aquaculture conditions leading to stress induced sex-reversal [START_REF] Geffroy | The Adaptive Sex in Stressful Environments[END_REF] , which has been observed in other fish species reared in laboratory or aquaculture facilities [START_REF] Pan | Identification of the master sex determining gene in Northern pike (Esox lucius) reveals restricted sex chromosome differentiation[END_REF] .

In the nine other species ( Table 1 and Supp. Fig. 3 ), we did not find any marker significantly associated with phenotypic sex, which indicates that these species either 1) have a small undifferentiated sex locus, 2) lack a genetic sex determination mechanism, or 3) rely on a complex combination of environmental and genetic sex determination factors that prevented the detection of a clear sex-biased signal in the dataset.

DISCUSSION

We demonstrate that RADSex is an efficient and, thanks to its simple usage and visualization options, user-friendly tool to identify sex-determination systems and sex-biased markers from RAD-Seq data. Unlike other popular RAD-Seq analysis software, e.g. STACKS (J.

M. [START_REF] Catchen | Stacks: building and genotyping Loci de novo from short-read sequences[END_REF] and PyRad [START_REF] Eaton | PyRAD: assembly of de novo RADseq loci for phylogenetic analyses[END_REF] , which catalog polymorphic loci, RADSex creates monomorphic markers. RADSex therefore counts each allele at a locus defined as polymorphic by other software as a separate marker. Working with monomorphic markers allows statistical comparison of populations using straightforward presence/absence tests, yielding results that are comparable between datasets. In contrast, the outcome of grouping sequences into polymorphic markers can be sensitive to parameter values [START_REF] Paris | Lost in parameter space: a road map for stacks[END_REF][START_REF] Rodríguez-Ezpeleta | Population structure of Atlantic mackerel inferred from RAD-seq-derived SNP markers: effects of sequence clustering parameters and hierarchical SNP selection[END_REF][START_REF] Shafer | Bioinformatic processing of RAD-seq data dramatically impacts downstream population genetic inference[END_REF] . Leveraging this property of monomorphic markers, RADSex offers a simple and reproducible way to identify de novo genetic differences between females and males and to locate these differences on a genome when a reference assembly is available. In addition to this simplicity, RADSex requires few resources and provides comprehensive visualization tools to assist the interpretation of its output.

In particular, the tile plot visualization provides a comprehensive overview of the RAD-Seq results across all individuals and allows the user to quickly identify straightforward cases of genetic sex determination systems. For instance, XX/XY or ZZ/ZW monofactorial systems can be easily detected in absence of any major environmental effects on sex, as demonstrated by our analyses on the ayu ( P. altivelis ) and the tench ( T. tinca ) datasets: in both species, the tile plot revealed a complete sex-linkage of several markers significantly associated with male phenotype, indicating a male heterogametic sex determination system. Besides these straightforward cases, the tile plot visualization is most useful to identify the sex-determination system in more complex datasets, for instance in the presence of outliers and sequencing biases. These complex datasets can then be further explored with the marker depths heatmap visualization, for instance to confirm whether sex-biased markers are absent or present in specific individuals. Using this heatmap, we were able to quickly identify sex-reversed individuals responsible for the non-complete sex-linkage of sex-biased markers in a publicly available medaka dataset [START_REF] Wilson | Wild sex in zebrafish: loss of the natural sex determinant in domesticated strains[END_REF] , as well as outliers in our own common carp ( C. carpio ), Atlantic cod ( G. morhua ), banded knifefish ( G. carapo ), and common molly ( P. sphenops ) datasets. These outliers could arise from human errors in assessing phenotypic sex; in particular, in the common molly, males are externally indistinguishable from females until puberty, which occurs late in some males, and immature testes resemble non-reproductively active ovaries. Because in this species sex was determined from secondary sex characters and macroscopy of the gonads, but not histology, it is possible that the relatively high number of females in which we found male-biased markers were immature males that were mis-sexed as females. However, the presence of outliers could also hint at a more complex sex determination system, for instance polygenic sex-determination [START_REF] Moore | Polygenic sex determination[END_REF] , autosomal modifiers that are polymorphic in the population [START_REF] Wu | The fate of autosomal modifiers of the sex-ratio trait in drosophila and other sex-linked meiotic drive systems[END_REF] , or simple monofactorial genetic sex determination systems with some sex-reversed individuals triggered by environmental factors (Baroiller & D'Cotta, 2016;[START_REF] Dupoué | Some like it dry: Water restriction overrides heterogametic sex determination in two reptiles[END_REF][START_REF] Wessels | ddRADseq reveals determinants for temperature-dependent sex reversal in Nile tilapia on LG23[END_REF] . With a high occurrence of outliers, sex-biased markers may no longer be significantly associated with phenotypic sex, yet these markers can still be biologically relevant.

Because RADSex creates non-polymorphic markers, the number and distribution of markers between individuals is only affected by a single parameter controlling the minimum depth to consider a marker present in an individual. An overly high minimum depth value can lead to false negatives, i.e. markers considered absent from an individual because of insufficient sequencing depth, and a very low value can create false positives, i.e. markers considered present in an individual because of sequencing errors. Moreover, because the probability of association with phenotypic sex for a marker is adjusted using Bonferroni correction, the stringency of the test for significance increases with the total number of markers identified in the dataset, which is affected by the minimum depth parameter. In the 15 datasets that we generated, median sequencing coverage in an individual was between 30x and 100x, except for three species for which it was lower than 30x ( Supp. Fig. 2 ). For all species except the three species with a low sequencing coverage, we found that minimum depth values between 5 and 10 provided a good balance between stringency and minimizing false negatives. However, some datasets may require using low minimum depth values, for instance the three datasets for which median individual sequencing depth was low, or species with large genome sizes like some sharks [START_REF] Hara | Shark genomes provide insights into elasmobranch evolution and the origin of vertebrates[END_REF] or some amphibians [START_REF] Sclavi | Genome size variation and species diversity in salamanders[END_REF] for which sequencing all individuals of a population with sufficient depth can be costly. In such cases, the memory efficiency of RADSex enables the analysis of large datasets with minimum depth as low as one, which can require a lot of memory with other software.

The numbers of markers significantly associated with phenotypic sex can provide an estimate of the size of the sex determining region, with a low or high number of sex-specific markers reflecting a small or a large non-recombining region, respectively. It is worth noting, however, that other factors can affect this number, for instance the level of differentiation between males and females within the non-recombining region, or bias in GC content within this region [START_REF] Sigeman | Insights into Avian Incomplete Dosage Compensation: Sex-Biased Gene Expression Coevolves with Sex Chromosome Degeneration in the Common Whitethroat[END_REF][START_REF] Smeds | Evolutionary analysis of the female-specific avian W chromosome[END_REF] potentially leading to an over-or under-representation of RAD-Seq markers depending on the restriction enzyme used for library construction. An extreme case is a complete absence of markers associated with phenotypic sex for any value of minimum depth, which we observed in nine of our datasets. These null results may be linked to two key features of sex determination in teleosts. First, early cytological studies have shown that the majority of teleost fish carry homomorphic sex chromosome pairs [START_REF] Devlin | Sex determination and sex differentiation in fish: an overview of genetic, physiological, and environmental influences[END_REF] , and more recent genomic analyses have revealed that small sex loci are common in teleosts, with the extreme case of a single SNP being the sole difference between the sexes in Fugu [START_REF] Kamiya | A trans-species missense SNP in Amhr2 is associated with sex determination in the tiger pufferfish, Takifugu rubripes (fugu)[END_REF] . In such cases when there is very little differentiation between the sex chromosomes, the fragmented resolution of RAD-Seq means that the the likelihood to obtain RAD sequences from the sex locus can be very low. Although this problem might be overcome by increasing the number of RAD-seq reads, for example by using a restriction enzyme that cuts more frequently, genome-wide approaches like pooled sequencing of males and females [START_REF] Feron | Characterization of a Y-specific duplication/insertion of the anti-Mullerian hormone type II receptor gene based on a chromosome-scale genome assembly of yellow perch, Perca flavescens[END_REF][START_REF] Gammerdinger | Comparative analysis of a sex chromosome from the blackchin tilapia, Sarotherodon melanotheron[END_REF][START_REF] Wen | Sex chromosome and sex locus characterization in the goldfish, Carassius auratus[END_REF] or individual whole genome sequencing [START_REF] Star | Genomic characterization of the Atlantic cod sex-locus[END_REF] may be necessary to identify the sex-determining region. Second, a complete absence of sex-biased markers can reflect a complex sex determination system, for instance a polygenic system or a system involving a strong environmental effect that would weaken the association between genetic markers and phenotypic sex. In fish, many environmental components can interfere with genetic sex determination mechanisms (Kikuchi andHamakuchi 2013, Heule et al., 2014). Although the visualization tools included in RADSex facilitate the identification of sex reversed individuals, the association between phenotypic sex and genomic regions may still be difficult to detect when sex-reversal is frequent, especially when associated with phenotyping uncertainty. For such systems, pooled sequencing may not be well-suited because male and female pools will each contain heterogeneous genotypes, and individual sequencing remains expensive despite decreasing sequencing costs. A hybrid strategy using both RAD-Seq and pooled sequencing may prove the most efficient for these complex sex-determination systems, as demonstrated in the case of the goldfish [START_REF] Wen | Sex chromosome and sex locus characterization in the goldfish, Carassius auratus[END_REF] , a species with a strong thermal effect on sex determination [START_REF] Goto-Kazeto | Temperature-dependent sex differentiation in goldfish: Establishing the temperature-sensitive period and effect of constant and fluctuating water temperatures[END_REF] . Finally, in some species sex is determined entirely by environmental factors and therefore does not involve any genetic mechanism [START_REF] Martínez-Juárez | Mechanisms related to sexual determination by temperature in reptiles[END_REF] .

One of the strengths of RAD-Seq is the ability to compare populations without a reference genome, which may not always be available for non-model species. However, when a genome is available, locating sex-biased markers provides information on the size and level of differentiation of the sex-determining region and can help to identify candidate genes involved in sex determination. The manhattan plot, circos plot, and contig plot visualizations included in RADSex quickly display the location of the sex-biased markers over the whole genome or over selected regions to identify the sex-determining region. Furthermore, the sex-bias metric and probability of association with sex computed with RADSex can be more effective than F ST computed on polymorphic loci to detect variability in the sex-determining region. This was illustrated by our reanalysis of a publicly available Japanese medaka dataset [START_REF] Wilson | Wild sex in zebrafish: loss of the natural sex determinant in domesticated strains[END_REF] , in which we detected multiple Y chromosome haplotypes that were not found in the original study. Intra-specific sex chromosome polymorphisms have been reported in different taxa, including mammals in which Y chromosome polymorphisms are widely used for human evolution studies [START_REF] Jobling | The human Y chromosome: an evolutionary marker comes of age[END_REF] and some fruit flies in which the Y chromosome polymorphism is thought to control male fitness [START_REF] Chippindale | Y chromosome polymorphism is a strong determinant of male fitness in Drosophila melanogaster[END_REF] . In fish, sex chromosome polymorphisms have been reported in different populations of the guppy, Poecilia reticulata , and in other closely related species of guppies [START_REF] Nanda | Sex chromosome polymorphism in guppies[END_REF] , and different sex chromosomes have even been found within the same species, for instance in tilapias [START_REF] Cnaani | Genetics of sex determination in tilapiine species[END_REF] . However, to the best of our knowledge, there are no reports of intra-specific population differences of a single sex chromosome in fish, and our results provide the first example of such a Y chromosome polymorphism in this clade.

As part of this study, we used RADSex to investigate the sex determination systems of 15 species sampled broadly across the ray-finned fish phylogeny (Sup. figure X). We identified XX/XY SD systems in six of these 15 species, and we did not find evidence of a ZZ/ZW SD system in any of the other nine species. Although our sampling is too limited to make general inferences about teleost sex determination, the predominance of XX/XY over ZZ/ZW SD systems in our datasets is in agreement with previous findings that in teleosts, transitions from ZZ/ZW to XX/XY SD systems are more frequent than the reverse [START_REF] Pennell | Transitions in sex determination and sex chromosomes across vertebrate species[END_REF] , which would result in XY SD systems being prevalent in this clade.

All our datasets were generated using single-digest RAD-Seq, but the workflow accepts double-digest RAD-Seq data as well. Furthermore, although all species included in this study are fish, we did not use any fish-specific assumptions in RADSex's implementation, and therefore the computational workflow can be applied to RAD-Seq data from any species. Finally, we specifically developed RADSex to study sex determination, and this design decision is reflected in the wording of this manuscript. However, the computational workflow was designed to be generic, and therefore both the radsex software and the sgtr R package could be used for other large and major quantitative trait loci with contrasting binary phenotypes. A table containing the sequencing depth of each marker in each individual is computed using the process command from radsex . This table is then used as input for the distrib command to compute the distribution of markers between males and females, and this distribution is visualized with the radsex_distrib() function of sgtr . Markers significantly associated with sex are extracted with signif and their depth in each individual is visualized with the radsex_markers_depth() function of sgtr. When a genome assembly is available, markers can be aligned with radsex map and results are visualized with the radsex_map_circos() function of sgtr . Each plot shows the distribution of RADSex markers between phenotypic males (horizontal axis) and phenotypic females (vertical axis), with color intensity for a tile indicating the number of markers present in the corresponding number of males and females. Tiles for which association with phenotypic sex is significant (Chi-squared test, p<0.05 after Bonferroni correction) are highlighted with a red border. A) Tile plot for an XX/XY sex determination (SD) system showing markers present in all males and absent from all females in the bottom right corner. B) Tile plot for a ZZ/ZW SD system showing markers present in all females and absent from all males in the top left corner. C) Tile plot for an XX/XY SD system with five male outliers showing markers present in all but five males and absent from all females at position (15, 0). D) Tile plot for a ZZ/ZW SD system with two male and three female outliers showing markers present in all but three females and absent from all but two males at position (2, 17). Each row corresponds to a marker, and each column corresponds to an individual. Males (blue labels) and females (red labels) clustered in two separate groups, except for two male outliers, identified by the green arrows, which clustered with the females. In addition, six males forming a sub-cluster within the male cluster had low or null depth for half of the markers. C : Manhattan plot showing -log10(P) of a chi-squared test on number of males and number of females, after Bonferroni correction, for all RADSex markers aligned to the genome of O. latipes , with odd linkage groups (LGs) in yellow with light background and even LGs in blue with shaded background. Unplaced scaffolds were combined in a super scaffold (U) by order of decreasing size. The dashed line indicates a p-value of 0.05 after Bonferroni correction. Most markers significantly associated with sex were aligned to LG01 (Ola1) or to unplaced scaffolds. D : The -log10(P) of a chi-squared test on the number of males and number of females is plotted against genomic position on LG01 of O. latipes (Ola1) for all RADSex markers aligned to this LG. The dashed line indicates a p-value of 0.05 after Bonferroni correction. Markers significantly associated with sex were aligned to a continuous region between 14.7 and 31.9 Mb. Markers present in all males apart from the two outliers and absent in all females were aligned in a continuous region spanning from 17 to 23 Mb (green box), while male-biased markers absent from the sub-cluster of six males were aligned in two regions, from 14 to 17Mb and from 24 to 31.9 Mb (purple boxes).

Marker

  2014 after sedation with an overdose of MS222 (200 mg/l) and final euthanasia by applying a sharp blow to the head. Phenotypic sex was determined based on the gonad macroscopy and fin clips were sampled and stored in ethanol until gDNA extraction. Sampling was carried out within the Norwegian animal welfare act guidelines, in accordance with the Animal Welfare Act of 20th December 1974, amended 19th June, 2009, at a facility with permission to conduct experiments on fish (code 93) provided by the Norwegian Animal Research Authority (FDU, www.fdu.no). As this experiment did not involve any treatment apart from euthanasia, no specific permit was required under the Norwegian guidelines. Ayu ( Plecoglossus altivelis ): Samples were collected in October 2013 in the Jinzu river at Toyama prefecture, Japan. Sexing was performed based on secondary sex characters (male dorsal and anal fins). Sample processing was performed with the approval of the Institutional Animal Care and Use Committee of the Tokyo University of Marine Science and Technology.

  Notothenia rossii ): Reproductively active individuals were collected in April-May 2013 and 2014 by bottom trawling and baited traps deployed from the ARSV Laurence M. Gould in the vicinity of Low Island and Dallmann Bay on the West Antarctic Peninsula. Live fish were transferred immediately from the trawl net or trap to the ship aquaria and transported to the aquatic facilities at Palmer Station, Antarctica, where they were maintained in flow-through seawater tanks at ambient temperatures of -1 to +1 °C. Following euthanasia using MS-222 (Syndel, Ferndale, WA, USA), individuals were dissected to identify the sex by gonad examination and a fragment of the gonads was preserved in Bouin's fixative for histological verification of the sex at the University of Oregon. A fin clip was also sampled and stored in 80% ethanol for gDNA extraction. All procedures were performed according to protocols approved by the Institutional Animal Care and Use Committees (IACUC) of the University of Oregon (#13-27RRAA). Tench ( Tinca tinca ) and common carp ( Cyprinus carpio ): Tenchs ( Tinca tinca ) were sampled in January 2017 from an aquaculture stock from the fish production unit of the fishing federation of Ille-et-Vilaine (Pisciculture du Boulet, Feins, France). A local aquaculture strain of common carp, Cyprinus carpio (U3E-INRA experimental aquaculture strain) reared in outdoor ponds was sampled in November and December 2013. For both species, phenotypic sex was recorded on each animal either by the presence of expelled semen or eggs or by macroscopic examination of the gonad after euthanasia with a lethal dose of MS-222 and dissection of the fish. Fin clips were taken on each animal and stored in ethanol 90% until gDNA extraction. For both common carp and tench, sampling and euthanasia procedures conformed to the principles for the use and care of laboratory animals, in compliance with French ("National Council for Animal Experimentation" of the French Ministry of Higher Education and Research and the Ministry of Food, Agriculture, and Forest) and European (European Communities Council Directive 2010/63/UE) guidelines on animal welfare. No specific permit was required under French guidelines for a simple euthanasia and tissue sampling. Genomic DNA (gDNA) extraction Genomic DNA (gDNA) of Lepisosteus oculatus , Plecoglossus altivelis , Gadus morhua , Notothenia rossii , Sander vitreus , Tinca tinca , Carassius auratus , and all Danios were extracted from fin clips stored in ethanol using NucleoSpin Kits for Tissue (Macherey-Nagel, Duren, Germany) following the producer's protocol. Genomic DNA of Gymnotus carapo , Gymnocorymbus ternetzi , Hippocampus abdominalis , and Poecilia sphenops were extracted from individual fin clips stored in ethanol or from pooled organs of individual fish for P. sphenops with a phenol/chloroform protocol. Fin clips or tissues were lysed in 1 ml extraction buffer (0.1 M

  RADSex's runtime was linearly and positively correlated to the number of RAD-seq reads, ranging from 139 seconds for Tinca tinca , whose dataset contained ~130 M. RAD-seq reads , to 1106 seconds (18 min and 26 s) for Sander vitreus , whose dataset contained ~262 M. RAD-seq reads ( Fig. 3.A ). RADSex was on average 42 times faster than STACKS on the same dataset, ranging from 26 times faster for Danio choprae to 88 times faster for Danio albolineatus ( Fig 3.B ). RADSex's peak memory usage was also linearly and positively correlated to the number of RAD-seq reads, ranging from 2.4 Gb for Poecilia sphenops to 17 Gb for S. vitreus ( Fig 3.C ) . RADSex used on average 1.9 times less memory than STACKS to analyze the same dataset, ranging from 1.2 times less memory for D. choprae to 3.4 times less memory for D. albolineatus ( Fig 3.D ). It is worth noting that for radsex, the peak memory usage and most of the runtime were in the process command, which needs to be run only once for each dataset because filtering parameters ( i.e.minimum depth for marker presence) can be applied in the later steps of the RADSex pipeline.

  clustered both markers and individuals separately based on these depth values and generated a heatmap of the results using the radsex_markers_depth() function of sgtr ( Fig. 4.B ). As expected, phenotypic males and females clustered in two separate groups, except for two phenotypic males (green arrows in Fig. 4.B ) that clustered with the phenotypic females,

  54 (23%) aligned to unplaced scaffolds, and 44 (19%) were not uniquely aligned with mapping quality higher than 20 ( Fig.4.C ). On Ola1, markers significantly associated with male phenotype aligned to a continuous region spanning from 14.7 Mb to 31.9 Mb ( Fig.4.D ). Within this region, markers found in all males and no females aligned between 17 and 23 Mb (green box in Fig.4.D) and markers present in all but the six males identified in the clustering step and absent from all females aligned between 14.7 Mb and 31.9 Mb (purple boxes in Fig. 4.D ). This result supports the hypothesis that multiple Y chromosome haplotypes are present in this medaka population, with one region located between 17 and 23 Mb on Ola1 showing strong X/Y differentiation in all Y haplotypes a wider region located between 14.7 and 31.9 Mb on Ola1 showing strong X/Y differentiation in only one Y haplotype.

FIGURESFigure 1 :

 1 FIGURES

Figure 2 :

 2 Figure 2: Tile plots showing the theoretical result of distrib in four scenarios.Each plot shows the distribution of RADSex markers between phenotypic males (horizontal axis) and phenotypic females (vertical axis), with color intensity for a tile indicating the number of markers present in the corresponding number of males and females. Tiles for which association with phenotypic sex is significant (Chi-squared test, p<0.05 after Bonferroni correction) are highlighted with a red border. A) Tile plot for an XX/XY sex determination (SD) system showing markers present in all males and absent from all females in the bottom right corner. B) Tile plot for a ZZ/ZW SD system showing markers present in all females and absent from all males in the top left corner. C) Tile plot for an XX/XY SD system with five male outliers showing markers present in all but five males and absent from all females at position (15, 0). D) Tile plot for a ZZ/ZW SD system with two male and three female outliers showing markers present in all but three females and absent from all but two males at position (2, 17).

Figure 3 :

 3 Figure 3: Performance of the RADSex workflow ( process, signif, and distrib ) and comparison with Stacks . A : Total runtime of RADSex (in minutes) plotted against the number of reads for each dataset. B : Peak memory of RADSex in gigabytes plotted against the number of RAD-seq reads for each dataset C : Runtime ratio between Stacks (version 2.2) and RADSex showing that RADSex was 25 to 80 times faster than Stacks . Datasets are sorted by number of RAD-seq reads. D : Peak memory ratio between Stacks and RADSex showing that RADSex used 1.2 to 3.4 times less memory than Stacks . Datasets are sorted by number of RAD-seq reads.

Figure 4 :

 4 Figure 4: Validation of the RADSex approach on the Oryzias latipes dataset . A : Tile plot showing the distribution of RADSex markers between males (horizontal axis) and females (vertical axis) using the original sex assignment and a minimum depth of m = 10. Color intensity indicates the number of markers present in the corresponding number of males and females, and a red border indicates significant association with sex (Chi-squared test, p<0.05 after Bonferroni correction). B : Heatmap showing individual depth for all markers present in at least 20 males and absent (-d 10) from all females.Each row corresponds to a marker, and each column corresponds to an individual. Males (blue labels) and females (red labels) clustered in two separate groups, except for two male outliers, identified by the green arrows, which clustered with the females. In addition, six males forming a sub-cluster within the male cluster had low or null depth for half of the markers. C : Manhattan plot showing -log10(P) of a chi-squared test on number of males and number of females, after Bonferroni correction, for all RADSex markers aligned to the genome of O. latipes , with odd linkage groups (LGs) in yellow with light background and even LGs in blue with shaded background. Unplaced scaffolds were combined in a super scaffold (U) by order of decreasing size. The dashed line indicates a p-value of 0.05 after Bonferroni correction. Most markers significantly associated with sex were aligned to LG01 (Ola1) or to unplaced scaffolds. D : The -log10(P) of a chi-squared test on the number of males and number of females is plotted against genomic position on LG01 of O. latipes (Ola1) for all RADSex markers aligned to this LG. The dashed line indicates a p-value of 0.05 after Bonferroni correction. Markers significantly associated with sex were aligned to a continuous region between 14.7 and 31.9 Mb. Markers present in all males apart from the two outliers and absent in all females were aligned in a continuous region spanning from 17 to 23 Mb (green box), while male-biased markers absent from the sub-cluster of six males were aligned in two regions, from 14 to 17Mb and from 24 to 31.9 Mb (purple boxes).

Figure 5 :

 5 Figure 5: Distribution of markers between males and females for six species (A-F) with markers significantly associated with sex . The number of markers, represented by tile color intensity, is plotted for every combination of number of males (horizontal axis) and number of females (vertical axis) for each species. Combinations of number of males and number of females for which association with sex is significant (p<0.05, Chi-squared test with Bonferroni correction) are highlighted in red. Distributions were computed with a minimum depth of 10 (-d 10).

Figure 6 :

 6 Figure 6: Results of RADSex analysis on Gadus morhua . A : Tile plot showing the presence of markers in each individual for all markers present with depth higher than five in at least 20 males and in fewer than two females. Each row corresponds to a marker and each column corresponds to an individual, with females in red and males in blue. B: Manhattan plot showing -log10(P) of a chi-squared test on the number of males and number of females for all RADSex markers aligned to the genome of G. morhua , with odd linkage groups (LGs) in yellow with light background and even LGs in blue with shaded background. Unplaced scaffolds were combined in a super scaffold (U) by decreasing order of size. The dashed line indicates a p-value of 0.05 after Bonferroni correction. C : -log10(P) of a chi-squared test on the number of males and the number of females for all RADSex markers aligned to LG11 of G. morhua . The dashed line indicates a p-value of 0.05 after Bonferroni correction. Markers significantly associated with sex (p < 0.05) and other markers strongly associated with sex were aligned to the first 20 Mb of the chromosome. The position of the candidate master sex determining gene is indicating by a red arrow.

see Fig. 4.B for an example from a real dataset). Markers obtained with the signif command can also be used as input for the radsex_markers_depth() function from sgtr if they were exported

  as a table of marker depths.

	Tiles for which association with sex is significant are highlighted with a red border. Simplified
	examples of tile plots for several sex-determining systems are presented in Fig. 2 .
	Next, RADSex can filter and cluster markers based on depth across individuals to
	identify consistent absence or presence of markers in some individuals. Markers can be filtered
	based on presence in males and females using the radsex subset command and exported in
	the same format as the table of marker depths ( Supp. Fig. 1 ). After extracting a subset of
	markers, a heatmap showing the depth of each marker in each individual can be generated
	using the radsex_markers_depth() function from sgtr . This function optionally allows the user to
	cluster both individuals and markers based on depth values and display the resulting
	cladograms along the heatmap. Distances are computed using the base R dist function and
	clustering is performed with the base R hclust	function; both distance calculation ( e.g.
	euclidean, maximum, binary ...) and clustering methods ( e.g.	complete, average, centroid ...)
	can be specified by the user. An application of the subset command is the detection of outliers
	in the sequenced population, e.g.	sex-reversed individuals or individuals with misassigned sex.
	In such cases, some markers may have biased presence in individuals from one sex without
	being statistically associated with phenotypic sex. When clustering individuals based on depth
	for these sex-biased markers, real outliers will be grouped with individuals from the other sex
	(		. Bonferroni correction is applied by
	multiplying the p-value by the total number of markers present in at least one individual. Markers
	for which association with sex is significant, i.e. p	< 0.05 (or a user-specified threshold) after
	Bonferroni correction, can be obtained with the signif command and exported either in the same
	format as the marker depths table or in fasta format. The results of distrib	( Supplementary
	table 2 ) can be visualized using the radsex_distrib() function from sgtr , which generates a tile
	plot with the number of males on the horizontal axis, the number of females on the vertical axis,
	and the number of markers represented by a tile's color for each number of males and females.

Table 5 .

 5

Danios: Several species of the genus Danio (Danio aesculapii, Danio albolineatus, Danio choprae, Danio kyathit): were obtained from Eugene Research Aquatics, LLC, Eugene, Oregon, USA. Individuals were euthanized in MS-222 for 10 minutes prior to dissection for sex identification and fin clipping for DNA extraction . Individuals were characterized as male or female based on gonad morphology.

  Spotted gar ( Lepisosteus oculatus ): Adult spotted gars (total length > 40 cm) were caught by electrofishing in 2014 and 2015 from natural populations around Thibodeaux, Louisiana (Bayou Chevreuil and Atchafalaya basin). Individuals were dissected to confirm sex based on gonadal morphology and fin clips were sampled and stored in 96% ethanol until gDNA extraction. Spotted gar sampling was approved by the University of Oregon Institutional Animal Care and Use Committee (Animal Welfare Assurance Number A-3009-01, IACUC protocol 12-02RA).

	euthanized with benzocaine hydrochloride solution of > 250 mg/L and preserved into glasses
	containing 100% ethanol for later gDNA extraction and histological analysis. Phenotypic sex
	was determined from histological processing of the gonads following standard protocols: the
	material was dehydrated in a series of alcohol, embedded in Leica historesin (methacrylate
	glycol) (Leica Microsystems, Wetzlar, Germany), sectioned with 3μm thickness, and stained
	with Hematoxylin/Eosin for examination under microscope. All procedures were consistent with
	Brazilian legislation regulated by the National Council for the Control of Animal Experimental
	(CONCEA) and Ethical Principles in Animal Research (Protocol n. 671-CEUA).

Pot-bellied seahorse ( Hippocampus abdominalis ): Adult pot-bellied seahorses were obtained from Seahorse Ways Co. Ltd., Kagoshima, Japan. Individuals were first euthanized in ethyl 3-aminobenzoate methanesulfonate (MS-222) for 5 min before tissue collection. Sex phenotypes were characterized based on gonadal morphology and fin clips were sampled and stored in 90% ethanol until gDNA extraction Common molly ( Poecilia sphenops ): Individuals were sampled from a laboratory strain (WLC No: 2418) kept by closed colony breeding, established from wild animals collected in the Rio Jamapa. Fish were kept and sampled in accordance with the applicable European Union and National German legislation governing animal experimentation. Experimental protocols were approved through an authorization (568/300-1870/13) of the Veterinary Office. Phenotypic sex was determined by macroscopy of gonads and secondary sex characters, including male pigmentation (black spotting, androgen dependent pattern, not expressed in females), male gonopodium (sex-specific anal fin morphological differentiation), and female gravidity spot.

Black widow tetra ( Gymnocorymbus ternetzi ): Individuals were collected in the Guaporé River (Pontes e Lacerda, Mato Grosso, Brazil) in December 2015. Fish were caught using trawl nets (3 mm mesh size, 25 x 2.5 m) which were placed at two points of the river. Samples were then Banded knifefish ( Gymnotus carapo ): Samples were collected along the Piracicaba River (Piracicaba, São Paulo, Brazil) in December 2014. The specimens were collected using sieves with square (2 x 1 m) and round (1.2 m diameter) shape meshes. Fish were then transported to the Department of Structural and Functional Biology, Institute of Biosciences, Botucatu (São Paulo, Brazil). Samples were euthanized with benzocaine hydrochloride solution of > 250 mg/L and preserved into glasses containing 100% ethanol for later gDNA extraction and histological analysis. Phenotypic sex was determined from histological processing of the gonads following standard protocols: the material was dehydrated in a series of alcohol, embedded in Leica historesin (methacrylate glycol) (Leica Microsystems, Wetzlar, Germany), sectioned with 3 μm thickness, and stained with Hematoxylin/Eosin for examination under microscope. All procedures were consistent with Brazilian legislation regulated by the National Council for the Control of Animal Experimental (CONCEA) and Ethical Principles in Animal Research (Protocol n. 671-CEUA).

Table 3 :

 3 example of map output for four markers. In order: marker ID, contig and position of the contig where the marker was aligned, sex bias, p-value of association with sex and significativity of this p-value for each marker.

	Contig	Position	Bias	P	Signif
	0	Chr1		241356	0.0689655 0.422281	False
	1	Chr5		17595862	-0.0851293 0.604401	False
	2	Chr7		136598	0.931035 1.69831e-12 True
	3	Scaffold1		365231	0.0150862 1	False
	4	Chr13		5397682	0.0592672 0.671211	False
	Supplementary Species	radsex		radsex peak	Stacks	Stacks peak
		runtime (s)	memory (Mb)	runtime (s)	memory (Mb)
	Cyprinus carpio	205		3,412		8,408	7,912
	Danio aesculapii	615		10,639		27,753	18,074
	Danio albolineatus	699		11,803		61,282	39,542
	Danio choprae	559		9,130		14,453	10,801
	Danio kyathit	650		8,866		24,851	16,460
	Gadus morhua	310		5,138		8,890	7,912
	Gymnocorymbus	648		15,566		55,996	20,091
	ternetzi
	Gymnotus carapo	499		7,944		13,152	11,845
	Hippocampus	198		3,173		6,101	8,716
	abdominalis
	Lepisosteus oculatus	383		4,956		13,584	6,293
	Notothenia rossii	258		5,026		8,692	8,087
	Plecoglossus altivelis	215		3,602		8,085	5,651
	Poecilia sphenops	158		2,612		7,431	6,667
	Sander vitreus	1,106		17,026		38,384	29,822
	Tinca tinca	141		3,381		4,831	8,062

Scientific name

ACKNOWLEDGEMENTS

We thank all colleagues who helped us improve the beta version of RADSex by testing it and providing feedback. This project was supported by funds from the "Agence Nationale de la Recherche'', the "Deutsche Forschungsgemeinschaft'' (ANR/DFG, PhyloSex project, 2014-2016, SCHA 408/10-1, MS ), the National Institutes of Health (USA) grant R01GM085318, JHP), and the National Science Foundation (USA) Office of Polar Programs (grants PLR-1247510 and PLR-1444167 to HWD and grant OPP-1543383 to JHP, TD, and HWD). The MGX core sequencing facility was supported by France Genomique National infrastructure, funded as part of "Investissement d'avenir" program managed by Agence Nationale pour la Recherche (contract ANR-10-INBS-09). RF was partially supported by Swiss National Science Foundation grant PP00P3_170664 to RMW. Common carp were provided by the PEARL INRA 1036 U3E experimental facilities that are supported by the ANAEE-France National Infrastructure. We thank Allyse Ferrara, Quenton Fontenot, and the Bayousphere Lab at Nicholls State University (Thibodeaux, LA) for support with spotted gar sample collection and

Software used in the analyses

All analyses were performed using version 1.0.0 of radsex and figures were generated using version 1.0.1 of the sgtr R package in R version 3.5.2. Both radsex and sgtr are released under GPLv3 license; releases and source code are available at https://github.com/SexGenomicsToolkit/radsex for radsex and DATA ACCESSIBILITY All RAD-Sequencing experiments have been submitted to Genbank under the bioproject PRJNA548074. A computational workflow implementing all the analyses performed in this study, including generating figures, is available at https://github.com/RomainFeron/paper-sexdetermination-radsex . RADSex is released under GPLv3 license; source code, installation instructions, and documentation are available at https://github.com/SexGenomicsToolkit/radsexex . sgtr is released under GPLv3 license and available at https://github.com/SexGenomicsToolkit/sgtr .

AUTHOR CONTRIBUTIONS

RF designed and implemented RADSex and sgtr, with feedback from QP, MW, and YG. YG, JHP, and MS designed the PhyloSex project and advised on results interpretation. QP, MW, BI, JA, AH, KK, ASR, KD, SK, CK, JHP, MS and YG participated in the analysis of the results. HP and LJ prepared libraries and performed the sequencing. RF, QP, JHP, MS and YG drafted the manuscript. RF, QP, RMW, JHP, MS, and YG revised the manuscript. YG, MS, JHP, EJ, SK, MW, MA, CW, BM, AA, TD, FWG, MK, HWD, MO, RN, TS, MN, AW, ØK and IB collected, sexed and/or extracted and prepared gDNA samples. All authors approved the final manuscript.

TABLES

Species

SD system identified with radsex

Number of sex markers (d=)

Cyprinus carpio XX/XY with male outliers 7 (d = 10)

Danio aesculapii

No SD identified 0 (d = 10)

Danio albolineatus

No SD identified 0 (d = 10)

Danio choprae

No SD identified 0 (d = 10)

Danio kyathit

No SD identified 0 (d = 10)

Gadus morhua XX/XY with outliers 3 (d = 10)

Gymnocorymbus ternetzi

No SD identified 0 (d = 2)

Gymnotus carapo XX/XY with female outliers 8 (d = 10)

Hippocampus abdominalis

No SD identified 0 (d = 10)

Lepisosteus oculatus

No SD identified 0 (d = 2)

Notothenia rossii

No SD identified 0 (d = 10)

Plecoglossus altivelis XX/XY 47 (d = 10)

Poecilia sphenops XX/XY with female outliers 7 (d = 10)

Sander vitreus

No SD identified 0 (d = 2)

Table 1: Summary of RADSex results for the 15 datasets analyzed. Species where a sex-determination system was identified with radsex are highlighted with a grey background.

The number of markers significantly associated with phenotypic sex is also given for each species, with the value of ' d ' indicating to the minimal number of reads for a marker to be considered present in an individual in the analysis.

SUPPLEMENTARY MATERIAL

Supplementary Figure 1: example illustrating the computation of a markers depth table with radsex process and the distribution of markers between sexes with radsex distrib . RAD-Seq read depth in each individual is summarized in a table of marker depths with the process command from radsex . This table is used by the distrib command to compute the distribution of markers between males and females for a given minimum depth value to consider a marker present in an individual (-d). Supplementary Table 1: example of process output for three individuals and three markers. In order: marker ID, marker sequence, and depth of this marker in each individual. For instance, the marker with ID 1 corresponds to the sequence TGCAGG..CT (the sequence was shortened), which is present 17 times in the sequencing results of Individual_1, 13 times in individual_2, and 15 times in Individual_3.

Supplementary

Males

Supplementary Table 6:

Excel file containing all significant sex-biased marker sequences and the distribution of these markers in all the animals analyzed in Cyprinus carpio (7 markers, m = 10), Gadus morhua (3 markers, m = 10), Gymnotus carapo (8 markers, m = 10), Plecoglossus altivelis (4 markers, m = 10), Poecilia sphenops (7 markers, m = 2) and Tinca tinca (6 markers, m = 10).

Supplementary Table 7: number of phenotypic males and females, minimum depth to consider a marker present in radsex , and number of markers for each dataset. Markers retained correspond to the number of markers present with depth higher than Min. depth in at least one individual.


MD Matthieu Olivetto
email: matthieuolivetto@live.fr

MD Jebrane Bouaoud

MD Jérémie Bettoni

An asymptomatic palatal tumor

Description

A 62-year-old male consults for an asymptomatic palatal swelling progressing for about 2 months.

He is a retired informatician and has no medical history. Clinical examination reveals a large, fleshcoloured, slightly ulcerated and painless palatal mucosal nodular mass measuring 3.8 centimetres in diameter (Fig. 1). There are no adenopathies in the cervical lymph node areas and the other aspects of the examination along with nasofibroscopy do not reveal any abnormalities.

An injected head and neck computed tomography does not show any associated osteolysis

The anatomopathological examination of the biopsy reveals a squamous mucosa whose chorion is massively infiltrated by epithelioid or elongated tumor cells, with eosinophilic cytoplasm sometimes pigmented (Fig. 2). The nuclei are oval with a prominent nucleolus. The stroma is thin and fibrous. The surface coating is focally ulcerated. Tumor cells are marked with anti-PS100, anti-MelanA, anti-HMB45 and anti-SOX10 antibodies. They are not marked with anti-CKAE1/AE3 antibody.

What is your diagnosis?

Response

It is a mucosal melanoma. These aggressive malignant tumours derived from melanocytes are rare (0.8 to 3.7% of melanomas) [START_REF] Yde | Melanoma: a Literature Review[END_REF] and frequently discovered at a late stage due to their hidden locations. In more than 50% of cases they concern the cervico-facial area (66% in sinuso-nasal and 25% in the oral cavity) [START_REF] Yde | Melanoma: a Literature Review[END_REF].

The median age of discovery of 70 years hides a very wide distribution. Their incidence is stable, without ethnic predominance, identified risk factors or a clinical specific presentations. If melanic lesions are evoking the diagnosis, in 40% of cases the lesion is amelanotic, flesh-coloured or greyish [START_REF] Paulo | Primary amelanotic mucosal melanoma of the oronasal region: report of two new cases and literature review[END_REF]. Symptomatology may include obstructive nasal mass or epistaxis events [START_REF] Dréno | Sinonasal mucosal melanoma: A 44-case study and literature analysis[END_REF]. Oral locations are more frequently palatal or gingival with lymphatic lymph node metastases frequent at diagnosis (33% of cases) [START_REF] Yde | Melanoma: a Literature Review[END_REF].

The anatomopathological analysis of the diagnostic biopsy does not systematically find tumour pigmentation. Immunohistochemistry allows diagnosis to be made using a series of markers, some of which are very sensitive (S100 protein) and others specific (MelanA and HMB45) [START_REF] Yde | Melanoma: a Literature Review[END_REF].

Therapeutic management requires a surgical excision with safety margins of 2 centimetres. Cervical neck dissection is recommended in patients with oral cavity disorders, but the lesser lymph node invasion of nasosinus melanomas does not justify lymph node dissection in the absence of clinical or radiological cervical adenopathy [START_REF] Yde | Melanoma: a Literature Review[END_REF].

Systematic adjuvant irradiation appears to improve local control [START_REF] Yde | Melanoma: a Literature Review[END_REF][START_REF] Dréno | Sinonasal mucosal melanoma: A 44-case study and literature analysis[END_REF]. The 5-year survival rates range from 0 to 31% for sinusonasal locations and from 15 to 45% for the oral cavity [START_REF] Yde | Melanoma: a Literature Review[END_REF]. No adjuvant systemic therapy has been shown to be effective. In the metastatic stage, it is the advent of targeted therapies that could improve prognosis in the coming years [START_REF] Yde | Melanoma: a Literature Review[END_REF]. The classical histological characteristics of prognostic relevance in cutaneous melanomas are not sufficient to classify mucosal melanomas [START_REF] Yde | Melanoma: a Literature Review[END_REF][START_REF] Dréno | Sinonasal mucosal melanoma: A 44-case study and literature analysis[END_REF]. Their molecular alterations and therapeutic targets are distinct [START_REF] Williams | Update from the 4th Edition of the World Health Organization Classification of Head and Neck Tumours: Mucosal Melanomas[END_REF].

Figure 1

 1 Figure 1 Click here to download high resolution image

Figure 2

 2 Figure 2Click here to download high resolution image

Conflict of Interest

The authors declare that they have no interest links.