Chromosome-Wide Evolution and Sex Determination in the Three-Sexed Nematode Auanema rhodensis

Trioecy, a mating system in which males, females and hermaphrodites co-exist, is a useful system to investigate the origin and maintenance of alternative mating strategies. In the trioecious nematode Auanema rhodensis, males have one X chromosome (XO), whereas females and hermaphrodites have two (XX). The female vs. hermaphrodite sex determination mechanisms have remained elusive. In this study, RNA-seq analyses show a 20% difference between the L2 hermaphrodite and female gene expression profiles. RNAi experiments targeting the DM (doublesex/mab-3) domain transcription factor dmd-10/11 suggest that the hermaphrodite sexual fate requires the upregulation of this gene. The genetic linkage map (GLM) shows that there is chromosome-wide heterozygosity for the X chromosome in F2 hermaphrodite-derived lines originated from crosses between two parental inbred strains. These results confirm the lack of recombination of the X chromosome in hermaphrodites, as previously reported. We also describe conserved chromosome elements (Nigon elements), which have been mostly maintained throughout the evolution of Rhabditina nematodes. The seven-chromosome karyotype of A. rhodensis, instead of the typical six found in other rhabditine species, derives from fusion/rearrangements events involving three Nigon elements. The A. rhodensis X chromosome is the smallest and most polymorphic with the least proportion of conserved genes. This may reflect its atypical mode of father-to-son transmission and its lack of recombination in hermaphrodites and males. In conclusion, this study provides a framework for studying the evolution of chromosomes in rhabditine nematodes, as well as possible mechanisms for the sex determination in a three-sexed species.

It is hypothesized that mixed systems such as androdioecy, gynodioecy (females/hermaphrodites) and trioecy are intermediate steps between dioecy and hermaphroditism in some systems (Weeks et al. 2006a). Although no gynodioecious nematodes are known, androdioecy evolved from dioecy multiple times during nematode evolution (Kiontke et al. 2004;Cho et al. 2004;Herrmann et al. 2006b;Pires-daSilva 2007). The transition from dioecy to hermaphroditism can be achieved in few steps. For instance, the down-regulation of only two genes (a sex-determining gene and a sperm-activation gene) in females of the dioecious Caenorhabditis remanei is sufficient to induce the development of selfing hermaphrodites (Baldi et al. 2009). Androdioecy has been described in a number of species, especially free-living, terrestrial nematodes (Maupas 1900;Sudhaus 1976;Herrmann et al. 2006a;Pires-daSilva 2007). However, it is unclear which evolutionary steps are necessary and whether this is an evolutionary stable mating strategy (Loewe and Cutter 2008;Chasnov and Chow 2002;Stewart and Phillips 2002;Cutter 2005).
By genetically manipulating C. elegans sex determination it is possible to model various mating systems (Hodgkin 2002). Mutant alleles of sex determination genes have been combined to create dioecious and trioecious strains (Cutter 2005;Stewart and Phillips 2002). Dioecious strains can be indefinitely maintained by introducing a mutation (e.g., fog-2 or spe-27) that prevents spermatogenesis in hermaphrodites, turning them into functional females (Schedl and Kimble 1988;Minniti et al. 1996). Synthetic trioecious populations, consisting of C. elegans males and females (feminised hermaphrodite mutants) mixed with wild type hermaphrodites are short-lived: males and females are rapidly outcompeted by selfing hermaphrodites (Stewart and Phillips 2002). This happens even under high mutational conditions, which are predicted to lead to a short-term advantage of obligate outcrossing over selfing (Cutter 2005). However, these trioecious populations were artificially created and it is possible that the males of C. elegans have lost some of their reproductive faculties. Naturally occurring trioecious nematodes, such as free-living nematodes of the genus Auanema (Maupas 1900;Félix 2004;Kanzaki et al. 2017) and entomopathogenic nematodes of the genus Heterorhabditis (Ciche 2007;Zioni Cohen-Nissan et al. 1992) are interesting systems to study the stability of trioecy and the mechanisms controlling the development of different sexual morphs within a same population. They also offer the opportunity to assess the costs and advantages of selfing vs. outcrossing.
To investigate the possibility of trioecy being stable in specific circumstances, we have been studying the sex determination system of Auanema nematodes (Shakes et al. 2011;Chaudhuri et al. 2015;Chaudhuri et al. 2011). A. rhodensis produces males, females and hermaphrodites, both by selfing and crossing (Félix 2004). A. rhodensis males are XO, whereas hermaphrodites and females are XX (Shakes et al. 2011;Chaudhuri et al. 2015;Chaudhuri et al. 2011). Due to a modified meiosis and spermatogenesis, A. rhodensis males produce functional spermatids with one X chromosome (haplo-X sperm), whereas sperm without the X chromosome (nullo-X sperm) are discarded (Shakes et al. 2011).
The meiosis program governing the X chromosome is also atypical in A. rhodensis XX hermaphrodites. We have previously shown that the X chromosome does not seem to recombine during hermaphrodite meiosis, leading to the production of mostly nullo-X oocytes and sperm with two X chromosomes (diplo-X sperm) (Tandonnet et al. 2018;Shen and Ellis 2018). Consequently, self-fertilization results mostly in XX progeny (hermaphrodites or females) ( Figure 1A). Males, having only haplo-X sperm, produce exclusively male cross-progeny when fertilizing the nullo-X oocytes of hermaphrodites (Tandonnet et al. 2018) ( Figure 1B). In females, the X chromosome recombines normally and thus most oocytes are haplo-X (nullo-X oocytes are also produced, although relatively rare). Cross-progeny of females with males are mostly hermaphrodites (Chaudhuri et al. 2015) ( Figure 1C).
It is not clear how the hermaphrodite vs. female sexual fate is determined. Selfing hermaphrodite mothers produce progressively fewer female progeny as they become older, whereas females produce almost exclusively hermaphrodite progeny throughout their lives (Chaudhuri et al. 2015). Hermaphrodites and females follow different developmental paths: hermaphrodites always pass through the dauer stage, whereas females (and males) do not (Félix 2004;Chaudhuri et al. 2011). Femalefated larvae that are forced to pass through the dauer stage (by the use of dauer-inducing chemicals) become adult hermaphrodites. On the other hand, blocking dauer formation of hermaphrodite-fated larvae results in the development of adult females . These results indicate that dauer formation is necessary and sufficient for hermaphrodite development in A. rhodensis .
Here we sequenced the genome of A. rhodensis with the motivation to uncover molecular mechanisms involved in the hermaphrodite vs. female fate, as well as to determine the consequences of the lack of recombination of the X chromosome in hermaphrodites and males. Previously, we determined that A. rhodensis has an unusual karyotype compared to most other Rhabditina nematodes, with six autosomes in addition to the X chromosome (Tandonnet et al. 2018). We show that the X chromosome is the smallest chromosome in A. rhodensis, is more polymorphic than autosomes and that the extra chromosome evolved from the fusion of parts of different ancestral chromosomes. Furthermore, we found that the product of the gene coding for the transcription factor dmd-10/11 seems to be required for the hermaphrodite fate.

Strains
We used A. rhodensis inbred strains APS4 and APS6 (Kanzaki et al. 2017) to produce F2-derived lines (F2Ls) and generate a genetic linkage map. Strains were maintained at 20°according to standard conditions as for C. elegans (Stiernagle 2006) on Nematode Growth Medium (3 g/L sodium chloride, 2.5 g/L bacto peptone, 17 g/L agar, 1 mM magnesium sulfate, 5 mg/L cholesterol, 1 mM calcium chloride, 25 mM potassium phosphate) (Brenner 1974). Plates were seeded with the Escherichia coli streptomycin-resistant strain OP50-1. Microbial contamination was prevented by adding 50 mg/mL of streptomycin and 10 mg/mL of nystatin to the NGM. The inbred strains were obtained by letting a population expand from a single individual hermaphrodite (picked at the dauer stage) every few generations. The strains APS4 and APS6 underwent 50 and 11 of such rounds of bottlenecking (expansion from a single hermaphrodite), respectively.
DNA extraction, sequencing and data pre-processing To extract nematode DNA with minimal bacterial contamination, we used the split plate method (Pires-daSilva 2013). Nematodes were cultured on one compartment of a 10 cm, two-compartment plate. The compartment with nematodes contained NGM seeded with E. coli OP50-1, and the second compartment contained M9 buffer. As the compartment with nematodes became crowded, dauer larvae migrated to the compartment with M9. The dauers were collected from 10 plates and washed twice with M9 buffer. The nematode pellet was stored at -80°. DNA was extracted and treated with RNAse using the Gentra Puregene Core A Kit (Qiagen) following the manufacturer's instructions. The DNA was dissolved in nuclease-free water for library preparation and sequencing.
The APS4 strain was chosen for genome assembly. Three independent Illumina paired-end (PE) libraries with insert sizes of 250 bp were sequenced at UT Southwestern (Dallas, Texas, USA) on an Illumina HiSeq 2500 (Table 1). Two Illumina mate-pair (MP) libraries with virtual insert sizes of $3 kb and $6 kb were constructed and sequenced on Illumina HiSeq 2500 at Edinburgh Genomics (University of Edinburgh, UK). The raw reads were processed to remove low-quality bases using Skewer (version 0.2.1, parameter settings "-Q 20 -l 51 -t 32") (Jiang et al. 2014). Error correction was performed using Fiona (version 0.2.1, "-nt 48 -g 60000000") (Schulz et al. 2014). We used Blobtools (Kumar et al. 2013) to remove microbial contamination.
To call variants in strain APS6, an Illumina paired-end library with insert size 450 bp was constructed and sequenced on Illumina HiSeq2500 at UT Southwestern (Table 1). Raw reads were preprocessed using Skewer as for the APS4 libraries (Jiang et al. 2014).

RAD-seq and Genetic Map Construction
A genetic map was constructed using markers obtained from restriction site-associated DNA sequencing (RAD-seq) of 95 F2-derived lines (F2Ls) originating from A. rhodensis inbred strains APS4 and APS6 (Kanzaki et al. 2017). To generate the F2Ls, crosses between APS4 females and APS6 males were performed to generate several F1 hermaphrodites, which were allowed to reproduce by selfing. Each progeny F2L was established from single F2 hermaphrodite progenitors, which were left to expand. For some lines, bottleneck events may have occurred after F2. We used the split plate method as described above to isolate DNA from the lines. This method relies on the isolation of dauers. Since dauers of A. rhodensis always develop into hermaphrodites Félix 2004), the DNA isolation was derived only from this sex. Paired-end RAD-seq using PstI restriction digestion was carried out for each of the parental strains and the 95 F2Ls (Baird et al. 2008). The raw RAD-seq reads were demultiplexed and lowquality regions were removed using process_RAD_tags from the Stacks package (version 1.35) (Catchen et al. 2013). We then used the denovo_map. pl Stacks pipeline to determine the genotype of each locus (region sequenced adjacent to the PstI cut site) for each progeny sample.
The genetic map was constructed using the R packages OneMap (version 2.0-4) (Margarido et al. 2007) and r/qtl (version 1.38-4) n  (Broman et al. 2003). A LOD (logarithm of odds) score of 20 and a recombination fraction of 0.5 were used as parameters to arrange the loci into linkage groups. The initial genetic map was refined by removing duplicated markers (markers with exactly the same genotype across all samples) and those with missing genotypes in 50% or more of the samples (function 'drop.markers'). Large gaps (loose markers) in the genetic map were fixed by dropping 3 markers. The Kosambi mapping function was used to determine the genetic distances between markers. However, the genetic distances could not be estimated precisely, as the level of recombination in the F2Ls is unknown.
Synteny analysis and identification of the X chromosome The software Chromonomer (version 1.07) (Amores et al. 2014) was used to anchor the genomic scaffolds to the genetic map, yielding a chromosomal assembly with scaffolds ordered, where possible, in each linkage group. The resulting chromosomal blocks were aligned to the C. elegans and Pristionchus pacificus genomes using PROmer (version 3.07) (Kurtz et al. 2004) with default parameters. Macro-synteny was visualized using Circos (version 0.69) (Krzywinski et al. 2009). One linkage group (LG5) aligned almost exclusively to the C. elegans X chromosome. We genotyped 5 polymorphic markers from this linkage group in F1 hybrid males (from an APS4 x APS6 cross) confirmed it to be the X chromosome (Tandonnet et al. 2018).
Intra-and inter-strain variant density was plotted along each chromosome using KaryoploteR (Gel and Serra 2017). The gene and variant densities of unanchored scaffolds (i.e., those not mapped to a linkage group) were not examined.
Gene ontology (GO) enrichment analyses were performed to examine possible GO terms found over-or under-represented in the X chromosome gene set vs. the autosomal one. For each enrichment analysis, we used a two-tailed Fisher's exact test (FDR , 0.05) implemented in the program Blast2GO (version 4.1.9) (Götz et al. 2008). The list of GO terms found enriched or depleted in the test set was then reduced to the most specific terms.
RNA extractions and transcriptomic analyses L2 females, converted females and hermaphrodites. Female-and hermaphrodite-fated L2 larvae were isolated by using synchronized progeny populations generated by hermaphrodite mothers. Briefly, dauers (fated to become hermaphrodites) were isolated and allowed to develop into adult hermaphrodites. After $12 h of egg laying the mothers were removed and the early eggs laid were left to hatch and grow until the L2 stage. During the L2 stage, females and hermaphrodites are distinguishable by their size and coloring: hermaphrodites are smaller, develop slower, are thinner and darker than females. Additionally, the female and hermaphrodite gonads are different in size during the mid-L1 stage, the female gonad being larger than that of the hermaphrodite. To convert hermaphrodite fated larvae into females, mid-L1 larvae with smaller gonads were isolated and grown with OP50-1 containing 200 nM dafachronic acid (DA) on NGM in individual wells of a 12 well culture plate. Once the larvae reached the L2 stage, the ones that had a female morphology were collected and used for RNA extraction. About 200 L2s of each sex were picked and transferred to an Eppendorf tube containing 200 mL of M9 buffer (Stiernagle 2006). The nematodes were washed 2-3 times in M9 buffer, allowing them to sink to the bottom of the tube under gravity between each wash. After the final wash, the maximum amount of M9 was removed and 200 mL of Trizol was added to the tube. The tube was placed at -80°immediately. For RNA extraction, nematodes were first freeze-cracked in liquid nitrogen (2-3 times). Trizol was added to make up the volume to 500 mL and nematodes were shaken with a few sterile 0.5 mm glass beads on a BeadBeater homogenizer (20 s, 3 times with 30 s intervals). Subsequently, RNA was extracted using a standard chloroform approach and the pellet was dissolved in DEPC treated water and stored at -80°until further use.
Males and mixed stages samples. The same protocol was used to extract RNA from adult males and animals from various stages (mixed stages). To obtain RNA from males, about 500 young adult males were picked in a 1.5 mL centrifuge tube containing 0.5 mL of M9 and washed twice with the buffer. M9 was then replaced with 0.5 mL Trizol and the tube frozen at -80°. For the mixed staged nematodes, M9 was gently added to 5 culture plates (6 cm) containing a healthy population of nematodes. We avoided disturbing the bacterial lawn and naturally let the nematodes start to swim in the buffer. This method reduced bacterial contamination in the samples during harvesting. These nematodes were then collected into a centrifuge tube washed twice with M9, and then frozen using liquid nitrogen. Tissues were homogenized for 1 min using a probe homogenizer. After chloroform extraction, the RNA was dissolved in DEPC treated water and stored at -80°.
We generated three biological replicates for each RNA-seq condition (L2 females, L2 converted females, L2 hermaphrodites, males and mixed stages). RNA-seq was performed on the Illumina HiSeq2500 platform, generating a mean of 19.7 million 100 base read pairs per replicate. General assessment of the RNA-seq libraries was performed using FASTQC (Andrews 2010). The raw reads from each library were preprocessed using Trimmomatic (version 0.36, "HEADCROP:15 SLIDINGWINDOW:5:20 MINLEN:20") (Bolger et al. 2014). The RNA-seq aligner STAR (version 2.4.2, "-sjdbOverhang 84") (Dobin et al. 2013) was used to align the processed reads of each library to the primary scaffolded genome assembly. Transcript abundances were obtained using FeatureCounts (from the SubRead Package, version 1.5.0-p2) (Liao et al. 2014). Differential expression between the L2 females, L2 converted females and L2 hermaphrodites (three comparisons) was assessed using the R package DEseq2 (version 1.18.1, (Love et al. 2014)) following the standard procedure and generating diagnostic plots (as described in the DESeq2 documentation at https:// bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/ DESeq2.html). An adjusted P-value of 0.01 and an absolute log2 fold change (FC) of 2 were used to define differentially expressed (DE) genes. Fold change of DE genes were plotted along the chromosomes using the package karyoploteR after lifting over the gene annotations. Gene ontology (GO) enrichment analyses were conducted on the "down-regulated" and "up-regulated" genes separately for each comparison using the procedure explained above. Homologs of known sex determination genes were identified through reciprocal best hit BLASTp searches (BLAST+, "-evalue 0.01 -max_target_seqs 100, -outfmt 6") using the A. rhodensis and C. elegans proteomes and analyzed manually.
Global protein-coding gene expression of L2 females, L2 hermaphrodites, adult males and mixed stages was examined separately for each chromosome. The transcript abundances obtained using FeatureCounts (from the SubRead Package, version 1.5.0-p2) (Liao et al. 2014) were corrected by library size. The log2 of the global gene expression of each chromosome was plotted using ggplot2 in R. To determine if the genes of the X chromosome were significantly less expressed than those on the autosomes, we randomly sampled the same number of autosomal genes and X genes (600) and compared the sets using a Kruskal-Wallis test followed by Wilcoxon-Mann-Whitney tests between autosome-X pairs. The gene expression across the same chromosome in different replicates of the same condition was confirmed to be similar by performing Kruskal-Wallis tests.

RNAi of the DM domain gene Arh-g5747
The DM domain gene Arh-g5747 was significantly up-regulated in hermaphrodites compared to females in both female and converted female samples at L2. To further investigate this sex-specific expression of Arh-g5747 we targeted the gene for down-regulation using RNA interference (RNAi). Target specific dsRNA was produced using a cDNA template. PCR amplification was performed using the following primers (Forward primer: 59-TAATACGACTCACTATAGGGTCAT-CAACGAGCAGAGCCGAGA-39, reverse primer: 59-TAATACGACT-CACTATAGGGTCCGCCTTCAGTGTTGGAGCT-39) to amplify an 858 bp fragment of the transcript of Arh-g5747. The T7 promoter (shown in italics above) was included at the 59 end of each primer to allow in vitro dsRNA synthesis. RNA extraction and cDNA synthesis were performed on $300 adult hermaphrodite individuals as detailed above, with the exception that samples were subjected to repeated cycles of freeze-thawing instead of bead-beating. RNA was treated with DNase I (Sigma) to remove residual genomic DNA. cDNA synthesis was performed with 0.5 mg of RNA using random primers (Promega) and the MMLV reverse transcriptase enzyme (Promega) following the manufacturer's instructions. The Arh-g5747 cDNA was then PCR amplified using GoTaq Green Master Mix (Promega), using approximately 200 ng of cDNA and 250 nM of each primer. PCR conditions were: an initial denaturation at 94°of 7 min, followed by 30 cycles of 94°for 15 s, 55°for 30 s and 72°for 60 s and a final extension of 10 min at 72°. After verification of the product size by gel electrophoresis, the amplicon was cleaned using the QIAquick PCR Purification kit (Qiagen), according to the manufacturer's protocol and eluted in a final volume of 25 mL. dsRNA was in vitro transcribed by incubating approximately 200 ng of the cDNA template with 2 ml (40 U) of T7 polymerase (Promega), 20 ml of 5x T7 polymerase buffer (Promega), 10 ml of DDT (Promega), 2.5 ml (100 U) of rRNAsin (Promega), 20 ml of 2.5 mM rNTPs (Thermo Fisher Scientific) and RNase free water to a final volume of 50 ml, for 4 h at 37°. The dsRNA product was size verified by gel electrophoresis and cleaned using the RNA clean-up protocol in the RNeasy Mini Kit (Qiagen). A mixture of short interfering RNAs (siRNA) was produced by digesting 5 mg of the Arh-g5747 dsRNA with ShortCut RNase III (NEB) for 20 min at 37°, cleaned by glycogen/ethanol precipitation and eluted in 20 ml of RNase free water, according to the manufacturer's instructions. The RNAi mixture was produced by combining 5 ml of Arh-g5747 siRNA (approximately 100 ng), 4 ul of M9 buffer (Stiernagle 2006) and 1 ml (10% v/v) LipofectamineRNAiMax reagent (Invitrogen) and incubating at 25°for 20 min. Previously, we have shown that inclusion of the transfection reagent Lipofectamine dramatically improves RNAi efficiency in A. rhodensis (Adams et al. 2019). For control injections Arh-g5747 siRNA was omitted from the mixture and replaced with additional M9 buffer.
For RNAi, young hermaphrodites (day 1 of adulthood) were immobilized on dried 2% agarose (w/v) pads in a small drop of undiluted halocarbon oil 700 (Sigma). The required injection mixture was loaded into pre-pulled microcapillary needles (Tritech Research) and microinjected into the gonad arms using an IM-300 Pneumatic Microinjector system with an oil hydraulic Micromanipulator (Narishige) using an injection pressure of 20 psi. Injected worms were rescued by adding a drop of M9 to the slide and moving them separately to a fresh 6 cm NGM plate seeded with E. coli OP50-1. The selfprogeny from the injected hermaphrodites was sexed throughout the life of the mother. The mother was placed on a new plate every 24 h. Sex was determined according to the developmental rate, coloration and morphology of the larvae. Females were larger and whiter than hermaphrodites (dark and thin) due to their faster development. Males displayed a characteristic blunt tail.

Genome characteristics
The scaffolded genome assembly spans 60.6 Mb in 440 scaffolds longer than 1000 bp. This span is smaller than that of C. elegans (100.2 Mb) but similar to that of other rhabditine nematodes (Heterorhabditis bacteriophora, 76.8 Mb; Oscheius tipulae, 59.0 Mb; Caenorhabditis sulstoni, 65.1 Mb). We predicted 11,570 protein coding genes and 833 unique non-coding RNAs. Previously sequenced rhabditine nematodes have been predicted to have more protein coding genes (C. elegans, 20,082; H. bacteriophora, 15,701; O. tipulae, 14,650). We presume that the difference in coding gene content is not due to a large number of missed genes, as the current assembly contains 99.19% of the core eukaryotic genes (predicted by CEGMA). The top BLAST hits of A. rhodensis proteins were more likely to be from parasitic Strongylomorpha species (Ancylostoma ceylanicum, Haemonchus contortus and Necator americanus) than from C. elegans. This is consistent with molecular phylogenies derived from ribosomal RNA and small numbers of protein coding loci (Kanzaki et al. 2017). However, it conflicts with analyses based on larger protein-coding gene datasets, which group A. rhodensis and the free-living O. tipulae closer to Caenorhabditis species than to strongylomorph species. A majority of the proteins (10,449, 90%) was assigned at least one type of annotation (InterProScan signature, GO term, BLAST hit) and 8181 (70.7%) were assigned with at least one GO term. The non-coding RNAs included all the expected major classes (Table S1). We also assembled and annotated the circular, 13,907 base pair mitochondrial genome ( Figure S1). The genome assembly has been submitted to ENA under the accession number ERS3049325 (SAMEA5241922).

Genetic map
We generated 95 F2-derived progeny Lines (F2Ls) from crosses between two polymorphic inbred strains of A. rhodensis (strains APS4 and APS6). We identified 1,052 polymorphic RAD-seq markers that clustered in 7 linkage groups (Table 3), presumably corresponding to the seven chromosomes in A. rhodensis identified by DAPI staining (Tandonnet et al. 2018). We anchored the genomic scaffolds of A. rhodensis to the genetic map to complete the sequence of each linkage group. Of the 1,052 markers and 636 scaffolds (. 200 bp), 1,038 markers ($94%) and 143 scaffolds ($22%) were used to build the chromosomal assembly. The excluded scaffolds were generally short, and either lacked a RAD marker or had a marker that was not able to be placed. The anchored scaffolds represent 95.3% of the span of the scaffold span and contain 93.8% of the predicted proteins (Table 3).
Macrosynteny with C. elegans and identification of the X chromosome A. rhodensis chromosomes are smaller than those of C. elegans (mean 8.3 Mb for A. rhodensis compared to 20.1 Mb for C. elegans), and one chromosome, LG5, is less than half the average size (3.5 Mb). To explore the origins of the changed complement of chromosomes and reduced size, we aligned the A. rhodensis protein-coding gene set to the one of C. elegans (Figure 2). While there was a minor background of between-chromosome translocation, most chromosomes had congruent gene sets. The majority of loci on three A. rhodensis linkage groups (LG1, LG6 and LG7) mapped to single C. elegans chromosomes (V, IV and II, respectively), suggesting one-to-one chromosomal correspondence. For A. rhodensis LG2, LG3 and LG4 we observed mapping to two or more chromosomes. Thus LG2 of A. rhodensis combines segments of C. elegans chromosomes III and X, LG3 combines segments of I and III, and LG4 combines segments of I, III and X. LG5, the smallest A. rhodensis chromosome, mapped almost entirely to the X chromosome of C. elegans, but segments of the C. elegans X chromosome were also found on LG2 and LG4 (Figure 3).
We confirmed that A. rhodensis LG5 was the X chromosome by genotyping F1 hybrid APS4/APS6 males at polymorphic loci spread across the linkage groups (Tandonnet et al. 2018). Hybrid males were heterozygous for all inter-strain polymorphic loci with the exception of those loci located in LG5, which were hemizygous in males (Tandonnet et al. 2018).

Chromosomal rearrangements
To further understand chromosomal rearrangements that took place in the lineage leading to A. rhodensis, we analyzed the synteny relationships of loci conserved between A. rhodensis linkage groups and the chromosomal assemblies of C. elegans, H. contortus, O. tipulae and Pristionchus pacificus (Figures 2 and 3, Figure S2).
A. rhodensis LG1 (Ar_LG1) contained loci that had orthologs on a single chromosome in C. elegans and H. contortus (chromosomes V/5; Ce_V, Hc_5) and on the left arm of P. pacificus chromosome I (Pp_IL) (Figure 2 and Figure S2). Comparing to O. tipulae, the orthologs of loci on Ar_LG1 were on chromosome X (Oti_X) ( Figure S2). A. rhodensis LG6 and LG7 had similar single-chromosome counterparts in the other species analyzed. Orthologs of loci on Ar_LG6 were found on Ce_IV, Ot_IV, Hc_4 and Pp_IV. Orthologs of loci on Ar_LG7 were found on Ce_II, Ot_II, Hc_2 and Pp_II (Figure 2 and Figure S2).
Orthologs of loci on Ar_LG2 were found on two distinct chromosomes in C. elegans, O. tipulae and H. contortus (Ce_III, Hc_3, Ot_III and Ce_X, Hc_X, Ot_V). On Ar_LG2, these loci were partially segregated into blocks with different chromosomal locations in the other species (Figure 4, first column). In P. pacificus, these same blocks of loci had orthologs segregated on the right arm of chromosome 1 (Pp_IR) and on Pp_III. We concluded that Ar_LG2 was the product of fusion and rearrangement of a fragment or fragments of an ancestral chromosome represented by Pp_IR, Ot_V, part of Ce_X and part of Hc_X and an ancestral chromosome now present as parts of Pp_III, Ot_III, Ce_III and Hc_3. The presence of interspersed, extended blocks of loci that appeared to derive from the same ancestral chromosome suggested that the rearrangement was relatively recent, as the processes of intrachromosomal inversion, known to be very rapid in rhabditine nematodes (Stein et al. 2003), had not yet mixed up these blocks of genes.
Analyses of Ar_LG3 and Ar_LG4 identified similar patterns of breakage and fusion. Ar_LG3 contained two sets of distinct blocks of loci, one that had orthologs on Ce_III, Ot_III, Hc_3, and Pp_III, and a second that had orthologs on Ce_I, Ot_I, Hc_1, and Pp_V n ( Figure 4, middle column). Ar_LG4 had three sets of blocks of loci with orthologs on three chromosomes in other species: one set on Ce_III, Hc_3, Ot_III and Pp_III, one on Ce_X, Hc_X, Ot_V and Pp_ IR, and one on Ce_I, Hc_I, Ot_I and Pp_V (Figure 4, last column).
The evolutionary history of the X chromosome The A. rhodensis X chromosome (LG5X) was the smallest chromosome (3.6 Mb) and had the lowest number of protein-coding genes (604, 5.5% of the total). X chromosomes differed markedly between species, but each contained orthologs of loci found on Ar_LG5X. For example, when compared to the most distantly related species, P. pacificus (Figure 5), the majority (82%) of the homologs of genes on Ar_LG5X were found on Pp_X ( Figure 5A, Figure S2C). However, Pp_X was five times larger (16 Mb) and contained 2,998 protein coding genes (11.7% of all predicted P. pacificus genes). C. elegans and H. contortus X chromosomes shared a striking pattern of macrosynteny with A. rhodensis. Mapping homologs from the C. elegans X to A. rhodensis, there were three distinct blocks of synteny on Ar_LG4 ( Figure 5D), and 5 blocks of synteny on Ar_LG2 ( Figure 5D). These same blocks were observed in comparisons with H. contortus ( Figure 5C). Intriguingly, while the C. elegans orthologs   LG4. The circos plots show macrosynteny relationships of (columns) LG2 (yellow), LG3 (pink) and LG4 (blue) of A. rhodensis to (rows) the chromosomal genomes of C. elegans, H. contortus, O. tipulae and P. pacificus, based on the mapping of presumed orthologs between the species. Each line in the circos plots corresponds to a predicted orthologous gene between the two species. In the second row, A. rhodensis is abbreviated 'A. rho.' of Ar_LG5X, Ar_LG4 and Ar_LG2 loci were evenly distributed across Ce_X ( Figure 5E), the mapping of the A. rhodensis orthologs on the H. contortus X was partitioned. Ar_LG4 matches were clustered on the left end of Hc_X and Ar_LG5X matches on the right ( Figure 5E). Ar_LG2 matches were more evenly distributed, although we can note a clustering on both ends of Hc_X. As discussed above, Ar_LG2 and Ar_LG4 may have originated through chromosome breakage and rearrangement. The segregation of Ar_LG5X-like and Ar_LG4-like regions on the H. contortus X may reflect conservation of ancestral synteny that has not been homogenized by within-chromosome rearrangement. The contrast between Ce_X and Hc_X, two chromosomes that otherwise appear highly homologous, suggested that either intrachromosomal rearrangement has been much more active in the lineage leading to C. elegans or that A. rhodensis and H. contortus shared a more recent common ancestor. A. rhodensis orthologs of genes on the O. tipulae X chromosome were found on Ar_LG5X and Ar_LG1 ( Figure 5B). This pattern was different from the one shared by mappings to P. pacificus, H. contortus and C. elegans and may reflect a novel trajectory of X chromosome evolution in the branch leading to O. tipulae.
Contrasting patterns of genome structure between the X chromosome and the autosomes We explored large-scale patterns of genome structure and evolution across the A. rhodensis genome. In C. elegans, conserved genes are more frequently found in the centers of chromosomes and are rarer in autosomal chromosome arms (Wilson 1999; C. elegans Sequencing Consortium 1998). However, in A. rhodensis the gene density and localization of genes with orthologs in Drosophila melanogaster across chromosomes was uniform ( Figure 6A). LG5X had a lower gene density than the autosomes (one protein-coding gene per 5.8 kb on LG5X compared to 5.3 kb 6 0.2 kb on the autosomes), and fewer conserved genes were present on LG5X ( Figure 6A, Table  S2). Most strikingly, none of the nearly 500 tRNA loci were on LG5X ( Figure 6B).
A gene ontology analysis comparing the X vs. the autosomal gene sets found that the GO terms 'translation', 'ribosome', 'nucleic acid binding', 'intracellular membrane-bounded organelle' and 'hydrolase activity' were under-represented on LG5X compared to the autosomes (Table S3). The process governing the 'neuropeptide signaling pathway' was found to be enriched on LG5X (Table S3).
Although the strains used in this study were inbred, we expected to observe a low level of within-strain heterozygosity due to incomplete inbreeding. We found that strain APS6 had higher heterozygosity than APS4, probably because APS6 underwent less inbreeding than APS4 (11 rounds of bottlenecking vs. 50) ( Figure 6C, Table S4). While the overall frequency of variants was different in the two strains, the distribution of these variants across the genome was similar ( Figure 6C), including shared chromosomal regions with higher natural variability. The LG5X chromosome displayed more within-strain variation than the autosomes, probably due to the atypical inheritance of this chromosome.
The genetic map displayed deviation from expected Hardy-Weinberg equilibria in several regions. We found that almost all RAD markers for LG5X were heterozygous across all 95 samples ( Figure 6D). This distorted pattern of X heterozygosity can be explained by the fact that the RAD data were derived from F2 hermaphrodite progenitors left to propagate for 3-10 generations, and there is no X recombination in hermaphrodites (Tandonnet et al. 2018). The few ($10%) markers that were homozygous for the X chromosome are probably the result of recombination that occurred in females during the population expansions originating from the F2s. We also observed a high frequency of homozygous markers for APS6 alleles at the right end of LG3 ( Figure 6D), suggesting that APS6 alleles had been positively selected in the culture conditions tested or possibly that segregation distorters were present. Less extreme deviation from expected equilibrium was also observed at the left end of LG7. These deviations were not explored further.
The different dosage of X chromosomes in females and hermaphrodites compared to males results in a requirement for dosage compensation for X-linked genes. We examined global protein-coding gene expression in L2 females, L2 hermaphrodites, adult males and mixed stages and compared autosomal and X-linked genes (Figure 7). After correction for library size, gene expression from each autosome was found to be similar both between autosomes and across lifecycle stages. However, and unlike C. elegans, genes on A. rhodensis LG5X showed consistently lower expression than those on the autosomes, even in the L2 stage (Wilcoxon Mann-Whitney, p-value ,= 1.0e-11 in all conditions and replicates).
Transcriptomic identification of loci associated with sexual morph development We compared gene expression in developing XX L2 larvae to identify loci that may be associated with the different sexual morphs of A. rhodensis. We generated replicate RNA-seq datasets from L2 fated to become hermaphrodites, L2 fated to become females, and L2 from hermaphrodite-fated nematodes that were converted to females by treatment with DA (converted females). Using standard thresholds (absolute log2(Fold Change) .= 2, FDR ,0.01), we found 2,422 (21%) of the predicted genes were differentially expressed (DE) between L2 hermaphrodites and L2 females. Slightly fewer genes (2,121,18%) were DE between L2 hermaphrodites and L2 converted females. Most of the genes found to be DE between females and hermaphrodites and between converted females and hermaphrodites were the same ( Figure  8A). The genes more expressed in females and converted females compared to hermaphrodites were enriched in GO terms related to translation, protein synthesis, ribosomal function, gonad and embryo development and structural constituents of the cuticle (see File S1 for the list of complete terms). Genes more expressed in hermaphrodites were enriched in few GO terms, with only "structural constituent of cuticle" in common between both comparisons. DE genes were distributed across the A. rhodensis genome, with no enrichment or depletion on LG5X ( Figure 8B). Some A. rhodensis orthologs of C. elegans sex determination genes were DE between female (normal or converted) and hermaphrodite L2s. The known sex determination genes Arh-g5696-gld-1 and Arh-g4999-tra-1 were 200 and 4 times more expressed in females, respectively. The precise roles of tra-1 and gld-1 in the sex determination in A. rhodensis are not yet known. A number of daf (dauer formation) genes (Arh-g6122-daf-11, Arh-g7695-daf-16, Arh-g7696-daf-like) were expressed at higher levels in hermaphrodite L2 compared to female or converted female L2, consistent with the obligate transition through dauer of hermaphrodite-fated nematodes.
The DM (doublesex/mab-3) domain transcription factor Arh-g5747 (dmd-10/11-like) was found to be more than 200 times more expressed in hermaphrodite L2 than in female or converted female L2. To investigate the role of this locus in the decision between female and hermaphrodite sexual fate in A. rhodensis, we downregulated Arh-g5747 by injecting RNAi in young hermaphrodites (first day of adulthood). If Arh-g5747 is required for determining hermaphrodite fate, we would expect to see more female progeny from injected hermaphrodite mothers. Indeed, downregulation of Arh-g5747 in 8 hermaphrodite mothers resulted in more female progeny than control injections performed on 9 hermaphrodites (Wilcoxon Mann-Whitney test, W = 67, p-value = 0.001563, Table 4). Thus, this DM domain locus may drive A. rhodensis hermaphrodite fate, either by inhibiting a female induction signal or through positive upregulation of a hermaphrodite-inducing pathway.
The comparison of the L2 females and the L2 converted females is particularly interesting for identifying genes or mechanisms involved in the hermaphrodite-female decision, upstream of the DA pathway. Gene expression in L2 females and L2 converted females was strikingly similar. Only 55 genes (0.5%) were found to be more highly expressed in females compared to converted females ( Figure 8A and 7B). No genes were found to be significantly less expressed in females compared to converted females. This result is surprising, since the female-inducing treatment (DA) was applied at the L1 stage, when sexual fate has already been decided, and the transcriptome was sampled less than 24 h after DA application. Functional annotation of these DE genes revealed several whose C. elegans homologs are involved in embryogenesis and developmental processes. Three chondroitin proteoglycan genes (Arh-g2548, Arh-g5439, Arh-g2211) were more expressed in female L2, and were also DE between female L2 and hermaphrodite L2. Chondroitin proteoglycans are important for embryonic cell division and vulval morphogenesis. More strikingly, we identified the homologs of the zinc-finger genes mex-1 (required for germ cell formation, and somatic cell differentiation in the early embryo in C. elegans) and pos-1 (essential for proper fate specification of germ cells, intestine, pharynx, and hypodermis in C. elegans). Both these genes had very low expression in converted female L2 and hermaphrodite L2. As these genes are maternally supplied in C. elegans, one possibility is that they are maternal regulators of sexual fate in A. rhodensis, although this hypothesis remains speculative.

DISCUSSION
Auanema rhodensis is a rare example of a three-sexed animal. Here we sequenced the A. rhodensis genome and used a linkage map to construct a chromosomal assembly. At 60 Mb, the A. rhodensis genome is smaller than that of C. elegans (100 Mb), but within the range (55-160 Mb) of other free-living rhabditomorph nematodes. We predict only 11,570 protein coding genes, many fewer than the 23,000 identified in C. elegans, and fewer than would be predicted from the reduction in genome size alone. However, considering that $99% of the core eukaryotic genes (CEGMA prediction) were identified, the reduced gene count in A. rhodensis is unlikely due to a high number of unannotated genes.
It is known that the mating system can influence genome size. In the Caenorhabditis clade, selfing (hermaphrodite) species have smaller genomes than their outcrossing sister species (Yin et al. 2018;Fierst et al. 2015), but these differences are in the order of 10%. However, other free-living and entomopathogenic rhabditomorphs have genomes smaller than C. elegans, and the related animal-parasitic Strongylomorpha have much larger genomes (250 -700 Mb). Detailed understanding of the evolutionary drivers of genome size in this group awaits additional, dense sampling across the Rhabditomorpha.
While there is little shared gene order between nematode species, we identified strong macro-syntenic patterns between A. rhodensis, P. pacificus, H. contortus, O. tipulae and C. elegans. These patterns allow us to propose a preliminary model of the evolution of chromosomes in the Rhabditina, which includes Diplogasteromorpha (P. pacificus), Strongylomorpha (H. contortus), and Rhabditomorpha (A. rhodensis, O. tipulae and C. elegans). It has long been noted that the majority of rhabditine nematodes have a karyotype of n= 6, with an XX:XO sex chromosome system (Walton 1959). While there are deviations from this pattern, including the trioecious A. rhodensis, with n= 7, and the variously parthenogenetic Diploscapter species with n = 1 to n = 9, the phylogenetic perdurance of this karyotype is striking. Using loci identified as orthologs in each species pair, we could identify six putative ancestral macrosynteny groups ( Figure 9) and also map the macrosyntenic changes that may have given rise to present day karyotypes.
We call these ancestral linkage groups Nigon elements in homage to Victor Nigon (Nigon and Félix 2017), a name coined by Matt Rockman (personal communication) in analogy with the Muller units of Drosophila chromosomes. Some of these Nigon elements have been transmitted intact through Rhabditina, while others have undergone rearrangement. Where a rearrangement has resulted in the fusion of (parts of) Nigon elements, in most cases the dynamic processes of intrachromosomal rearrangement, which are very active in rhabditine nematodes (Stein et al. 2003), have acted to mix up the genes originally derived from different units. In other cases, either because the fusions were more recent or because the processes of intra-chromosomal rearrangement are less active, the sets of loci from distinct Nigon elements are found as blocks in the fusion chromosome. We define six Nigon elements, NA, NB, NC, ND, NE and NX, as well as an additional NN unit which we are currently unable to place (Figure 9). It could originally link to NE or NX, or be a unit of its own (which would imply seven Nigon elements in total).
A. rhodensis LG1, LG6 and LG7 represent chromosomal units that have survived largely intact through rhabditine evolution (Nigon elements NE, ND and NB, respectively). However, the A. rhodensis X chromosome (LG5X) was formed from a subset of the loci now found on the C. elegans X chromosome, and A. rhodensis LG2, LG3 and LG4 are the products of major interchromosomal rearrangement events. For example, NA is intact in C. elegans (Ce_I), H. contortus (Hc_1), O. tipulae (Ot_I) and P. pacificus (Pp_V), but underwent fission in the A. rhodensis lineage, with subsequent fusion forming two hybrid chromosomes. A. rhodensis LG3 is formed largely from loci from part of NA (subset a1) and NC (c2), while LG4 includes loci from NA (a2), NC (c3) and NN (n2). Overall, compared to the other four species with chromosomal assemblies (or chromosome-allocated scaffolds), the fission/fusion event(s) may be directly associated with the origin of the novel n = 7 karyotype of this species. It will be informative to explore the origins of other species where n6 ¼ 6, such as species in the genus Diploscapter. In A. rhodensis, these events may have been relatively recent, phylogenetically speaking, as there are still clear blocks of genes of different Nigon element origin within the fusion chromosomes LG2, LG3 and LG4.
Fusions are not unique to A. rhodensis. NE is intact in A. rhodensis (Ar_LG1), C. elegans (Ce_V) and H. contortus (Hc_V) but has fused with NN (n1n2) in P. pacificus to form Pp_I. As noted previously (Rödelsperger et al. 2017), Pp_I is a fusion chromosome incorporating components of Ce_V (NE-derived) and Ce_X (NN-derived), but we note that the continued distinctiveness of the NE-derived and NNderived components within Pp_I suggests that this fusion is phylogenetically recent. The two Nigon element components within Pp_I retain an arms-and-centers long-range structure that is presumably derived from the original separate chromosomes, with high repeat density in the ancestral arms and high gene density in the ancestral centers. It has been proposed that the NE-NN fusion observed in P. pacificus is ancient, based on identification of a NE-NN like junction fragment in the genome sequence of the tylenchine (Clade IV) nematode Bursaphelenchus xylophilus, which is an outgroup to the rhabditine species. This apparent conservation of the junction fragment conflicts with the within-chromosome rearrangements dynamic observed elsewhere in the genome, and may be a chance homoplastic association of NE and NN elements in this species. In O. tipulae, NE has fused with NX to form Ot_X, and the NN unit forms a chromosome on its own (Ot_V).
The X chromosomes of the species analyzed always contain the NX unit either as the sole component of the X (Pp_X and Ar_LG5X), or associated with other Nigon elements: NN in Ce_X and Hc_X, and NE in Ot_X. The complex history of the NN, NX and NE units requires additional analyses, as it is unclear if NN belongs to NE (as found in Pp_I), to NX (as found in Ce_X and Hc_X) or if it is a unit of its own (as found in Ot_V). The number of chromosomally-assembled rhabditine genomes is still too few to fully define the ancestral gene content of Nigon elements.
Meiosis of the A. rhodensis X chromosome seems to follow different patterns, largely depending on the organismal sex and type of gametogenesis. By using five polymorphic markers, we have previously shown that the homologous X chromosomes undergo meiotic recombination in females, but not in hermaphrodites (Tandonnet et al. 2018;Shen and Ellis 2018). The genetic linkage map produced in the present study, using 92 markers for the X chromosome, confirms that the lack of recombination in the X in hermaphrodites is chromosome-wide ( Figure 6D). Lack of recombination of the X is observed during hermaphrodite oogenesis and spermatogenesis, leading to nullo-X oocytes and diplo-X sperm (Tandonnet et al. 2018). Additionally, during Figure 6 Contrasting genomic patterns between the X chromosome and the autosomes. (A) Distribution and conservation of protein-coding genes across A. rhodensis chromosomes. Density of A. rhodensis protein coding genes (upper panel) and conserved genes between A. rhodensis and D. melanogaster (lower panel) along each linkage group using a 200,000 bp window size. Overall, 4,544 conserved genes were identified between A. rhodensis and D. melanogaster. (B) Localization of the annotated non-coding genes (upper panel) and transfer RNAs (tRNAs) (lower panel) along each linkage group using a 300,000 bp window size. No tRNAs were found on the X chromosome (LG5). (C) Patterns of variation across two inbred strains of A. rhodensis. Variant density along each chromosome in 250,000 base windows for the within-strain variants (upper panels) and 100,000 base windows for the between-strain variants (lower panel). (D) Genotype frequencies across A. rhodensis chromosomes. Black and red lines represent the frequencies of RAD sites homozygous for the APS4 or for the APS6 allele, respectively, in the 95 genotyped F2Ls. Blue lines represent heterozygous genotypes. More than 80% of the progeny samples were heterozygous for the X chromosome (LG5X). The black ticks on the x-axis show the positions of the 1052 mapped RAD markers.
outcrossing the X chromosome is always transmitted from father to son (males produce exclusively haplo-X sperm). One of predictions from this atypical inheritance is that the genes on the X will be more exposed to selection. Thus, essential genes will tend to migrate from the X to autosomes, leading to a reduction in size. The X chromosome of A. rhodensis has many distinctive features. It is much smaller than the autosomes, representing only 6% of the genome and containing just over 600 genes. It has no tRNA genes, and fewer conserved genes were found on the X compared to the autosomes. In C. elegans, the X chromosome carries 44% of all tRNA genes and there is no marked exclusion of conserved genes (Wilson 1999). In A. rhodensis, the X is inherited from father to son and is haploid in males. Thus, genes on X chromosomes transmitted between males will be more exposed to natural selection, which will tend to exclude essential genes from the X (Tandonnet et al. 2018). In addition, the lack of recombination of the X chromosome in hermaphrodites will slow down the removal of deleterious mutations, contributing to the exclusion of essential genes on the X. The lower prevalence of essential genes on the X chromosome was reported in C. elegans (Kamath et al. 2003), where a genome-wide RNAi analysis revealed that the X chromosome was depleted of essential genes. The C. elegans X chromosome, which is also haploid in males, would also be more exposed to natural selection than autosomes, although to a lesser degree than A. rhodensis. The heightened exposure to selective forces and lack of recombination of the X would predict a lower diversity on the X due to genetic hitchhiking. Indeed, the presence of a beneficial allele on the X would quickly spread through the population drawing along the rest of the chromosome (selective sweep), and, correspondingly, the negative selection of a deleterious allele on the X would also lead to a decrease of the genetic diversity of the whole X (background selection). However, populations of A. rhodensis are composed of a high proportion of selfing hermaphrodites (estimated around 60% of the adult fraction of the population), in which the X chromosome does not recombine. The XX progeny resulting from a selfer therefore always retain maternal heterozygosity on the X (Tandonnet et al. 2018). Using the inbred strains APS4 and APS6, we found that within-strain and between-strain genetic diversities were higher on the A. rhodensis X chromosome than on the autosomes. In our inbreeding protocol, bottlenecking was performed by isolating a single selfing hermaphrodite every few generations. Thus, the X chromosome will only have recombined in females during the population expansion from each isolated hermaphrodite. As several generations occurred between each hermaphrodite isolation, we expect that the X would become homozygous at a much slower rate than the autosomes. In nature, the genetic diversity of the X chromosome compared to the autosomes (whether higher or lower) could depend on the proportion of the different sexes and on the effect of X-linked mutations.
The X chromosome gene expression was also found to be consistently lower than that of autosomes. In C. elegans hermaphrodites, the X to autosomal gene expression ratio (X:A ratio) varies through development from 0.92 in the L2 to 0.41 fold in adults (Xiong et al. 2010). This change is likely to be associated with the exclusion of genes with germline expression from the C. elegans X chromosome and the growth of the germline in adults (Gama et al. 2002;Strome et al. 2014). Indeed, as the individual (XO or XX) develops, the proportion of germ cells Figure 7 Expression of genes on the A. rhodensis X chromosome is generally lower than those on autosomes. Boxplots of the log2 normalized expression of the genes located on each linkage group of A. rhodensis in different sexes and stages in single replicate libraries. The expression levels were normalized by library size and log2-transformed.
LG5 (in red) is the X chromosome. Boxplots for all libraries are represented in Figure S4. This plot was generated using the R package ggplot2 (Wickham 2016). increases. Because the X chromosome is repressed in germ cells, the X:A ratio steadily declines as the individual develops (Deng et al. 2011). However, in C. elegans mutants lacking germline proliferation, the X:A ratio is close to 1, due to a dosage compensation mechanism equalizing the X chromosome and autosomal gene expression (Deng et al. 2011). In A. rhodensis, the X chromosome expression was found lower than autosomal expression even at the L2 stage. Based on these observations, it is possible that the X chromosome of A. rhodensis lacks a dosage compensation mechanism to equalize the X:A ratio (unlike C. elegans). However, it is to note that the X expression seems similar between XO males and XX animals, indicating that some dosage compensation mechanism is acting to prevent a higher expression in XX animals compared to XO males. Considering the low gene number on the X, one possibility is that a low X:A ratio is viable in A. rhodensis.
Another fundamental question in A. rhodensis biology is the mechanism that controls female vs. hermaphrodite sex determination in XX   Figure 6A) is shown under the lower panel.
animals. Females and hermaphrodites are karyotypically identical, and are thought to be genetically identical. This is because hermaphrodites of a strain inbred for 50 generations (APS4) still produce hermaphrodite and female progeny (Chaudhuri et al. 2015). Transcriptome comparisons between female, hermaphrodite and converted female early larvae (L2) show that the sex differentiation process modulates the expression of many genes, with $20% of all genes found differentially expressed between females (normal and converted females) and hermaphrodites. The considerable difference in transcriptomic profiles at L2 is reflected in their developmental trajectory, with the female-fated larvae undergoing faster development toward adulthood, whereas hermaphrodite-fated larvae arresting at the dauer stage. Additionally, physiological and metabolic differences between females and hermaphrodites are likely to be present. One established example is the production of male-attracting pheromones by females, but not in hermaphrodites (Chaudhuri et al. 2015). Of particular interest, we noted that A. rhodensis orthologs of some genes previously shown to be active in C. elegans sex determination were differentially expressed in females vs. hermaphrodites. A. rhodensis gld-1, for example, was 200-fold overexpressed in females. In C. elegans hermaphrodites, gld-1 is necessary for normal oogenesis and promotes spermatogenesis in hermaphrodites (Jan et al. 1999;Jones and Schedl 1995). However, it has the opposite role in C. briggsae (Beadell et al. 2011;Beadell and Haag 2014), probably due to changes in target transcripts. In C. elegans, GLD-1, and its cofactor FOG-2, regulate the translation of tra-2 and are necessary for hermaphrodite sperm fate (Hu et al. 2019). In C. elegans, wild-type and fog-2 mutants have very similar transcriptomic profiles, which is consistent their role at the translational level (Hu et al. 2019). In A. rhodensis however, no fog-2 ortholog was found and the L2 female and hermaphrodite transcriptomes are strikingly different, indicating that the regulation of sexual fate is different from C. elegans. A. rhodensis tra-1, a master sex determination gene, was four fold overexpressed in females. The zinc finger protein TRA-1 is the terminal regulator of the sex determination pathway in C. elegans, where it promotes female development in somatic tissues (Hodgkin 1987;Zarkower and Hodgkin 1992). Its role in sex determination is n Table 4 Inhibition of the DM transcription factor Arh-g5747 (dmd-10/11) by RNAi in hermaphrodite mothers results in more female progeny conserved in nematodes, making it an interesting target for functional studies (Pires-daSilva and Sommer 2004). We also identified a DM (doublesex/mab-3) domain transcription factor, homologous to C. elegans dmd-10 and dmd-11 that was 200-fold overexpressed in hermaphrodites. DM domain genes regulate sex determination and sexual differentiation processes in a number of organisms (Zarkower 2013), but specific roles for C. elegans dmd-10 and dmd-11 have not yet been elucidated. RNAi knockdown of this locus in hermaphrodites resulted in the production of more female progeny, suggesting a role for this DM domain-coding gene in A. rhodensis sex determination. The near identical expression pattern between converted females and females shows that the conversion of L1 hermaphrodite-fated larvae by exposure to DA is almost complete and that the initial sex decision can be overridden almost fully by hormonal manipulation. The age and sex of an A. rhodensis mother affect the proportion of each sex in its progeny (Chaudhuri et al. 2015) and thus it is likely that maternal effects may directly establish the distinct developmental trajectories of females and hermaphrodites. However, the female vs. hermaphrodite decision could also be modulated by environmental cues acting during embryogenesis and the L1 stage. These maternal and environmental effects could be modulations of what is essentially a random sex determination (RSD) system (Perrin 2016). RSD occurs when fluctuations in the expression of genes at the top of the sex determining cascade or "developmental noise" are enough to canalize sexual fate down contrasting paths. An RSD component in A. rhodensis is plausible as females and hermaphrodites likely share the same genome, all sexual morphs are produced in a single environmental condition and the proportion of each sex produced varies greatly between mothers (although they are from inbred lines).
Auanema nematodes thus offer a fascinating and potentially highly informative model system for depth exploration of the origin of novel traits and their consequences. The genomic and transcriptomic resources we present for A. rhodensis will be critical for future analyses of the origins of new chromosomes in an otherwise stable karyotypic system, the biology of the highly regulated pattern of X chromosome segregation, the dynamics of mating system evolution, and the evolution of sex determination mechanisms. Toward this, we have identified the orthologs of 16 main sex determination genes, including tra-1/2, gld-1, her-1 and fem-1/2 (Table S5), which are clear candidates to investigate the sex determination mechanisms in A. rhodensis. In parallel, we are developing reverse genetic and functional genomic technologies for these species, and these promise routes to rapid validation of hypotheses of gene function (Adams et al. 2019). The A. rhodensis sex determination system may integrate genetic, maternal, environmental and random components, and this nexus of interacting components will also become amenable to manipulation and dissection. Genetic and genomic investigation of additional Auanema and closely related rhabditomorph species will contribute to a complete understanding of the origins and maintenance of this unusual mating system. ACKNOWLEDGMENTS G.D.K. was supported by a BBSRC Ph.D. studentship and S.T. by a Ph.D. training grant from CAPES/CNPq (201116/2014-6). A.P.S. was supported by a grant from National Science Foundation (IOS 1122095), BBSRC (BB/L019884/1) and University of Warwick start-up funds. The authors are very grateful to Stephen Doyle (from the Sanger Institute) for providing early access to the Haemonchus contortus genome and annotation, and to Fabrice Besnard for discussion about O. tipulae linkage groups. We thank staff of Edinburgh Genomics and UT Southwestern facilities for support in sequencing. We are also grateful to two anonymous reviewers for their constructive comments on an earlier version of the manuscript.