
==== Front
Genetics
Genetics
genetics
Genetics
0016-6731
1943-2631
Oxford University Press US

38431281
10.1093/genetics/iyae032
iyae032
Investigation
Molecular Genetics of Development
AcademicSubjects/SCI01180
AcademicSubjects/SCI01140
The contribution of an X chromosome QTL to non-Mendelian inheritance and unequal chromosomal segregation in Auanema freiburgense
Al-Yazeedi Talal School of Life Sciences, University of Warwick, Coventry CV4 7AL, UK

Adams Sally School of Life Sciences, University of Warwick, Coventry CV4 7AL, UK

Tandonnet Sophie School of Life Sciences, University of Warwick, Coventry CV4 7AL, UK

Turner Anisa School of Life Sciences, University of Warwick, Coventry CV4 7AL, UK

Kim Jun Institute of Molecular Biology and Genetics, Seoul National University, Seoul 08826, South Korea

https://orcid.org/0000-0002-6421-1195
Lee Junho Institute of Molecular Biology and Genetics, Seoul National University, Seoul 08826, South Korea

https://orcid.org/0000-0002-0741-7197
Pires-daSilva Andre School of Life Sciences, University of Warwick, Coventry CV4 7AL, UK

Engebrecht J Editor
Corresponding author: School of Life Sciences, University of Warwick, Coventry CV4 7AL, UK. Email: andre.pires@warwick.ac.uk
Conflicts of interest: The author(s) declare no conflicts of interest.

Present address: Center for Applied and Translational Genomics (CATG), Mohammed bin Rashid University of Medicine and Health Sciences, Dubai, United Arab Emirates
Present address: Department of Convergent Bioscience and Informatics, Chungnam National University, Daejeon 34134, South Korea
5 2024
02 3 2024
02 3 2024
227 1 iyae03224 12 2023
15 2 2024
21 3 2024
© The Author(s) 2024. Published by Oxford University Press on behalf of The Genetics Society of America.
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs licence (https://creativecommons.org/licenses/by-nc-nd/4.0/), which permits non-commercial reproduction and distribution of the work, in any medium, provided the original work is not altered or transformed in any way, and that the work is properly cited. For commercial re-use, please contact reprints@oup.com for reprints and translation rights for reprints. All other permissions can be obtained through our RightsLink service via the Permissions link on the article page on our site—for further information please contact journals.permissions@oup.com.

Abstract

Auanema freiburgense is a nematode with males, females, and selfing hermaphrodites. When XO males mate with XX females, they typically produce a low proportion of XO offspring because they eliminate nullo-X spermatids. This process ensures that most sperm carry an X chromosome, increasing the likelihood of X chromosome transmission compared to random segregation. This occurs because of an unequal distribution of essential cellular organelles during sperm formation, likely dependent on the X chromosome. Some sperm components are selectively segregated into the X chromosome's daughter cell, while others are discarded with the nullo-X daughter cell. Intriguingly, the interbreeding of 2 A. freiburgense strains results in hybrid males capable of producing viable nullo-X sperm. Consequently, when these hybrid males mate with females, they yield a high percentage of male offspring. To uncover the genetic basis of nullo-spermatid elimination and X chromosome drive, we generated a genome assembly for A. freiburgense and genotyped the intercrossed lines. This analysis identified a quantitative trait locus spanning several X chromosome genes linked to the non-Mendelian inheritance patterns observed in A. freiburgense. This finding provides valuable clues to the underlying factors involved in asymmetric organelle partitioning during male meiotic division and thus non-Mendelian transmission of the X chromosome and sex ratios.

epistasis
transgression
asymmetric cell division
trioecy
sex determination
sex ratio
BBSRC 10.13039/501100000268 BB/L019884/1 Leverhulme Trust 10.13039/501100000275 RPG-2019-329 Ciência sem Fronteiras 10.13039/501100017564 201116/2014-6 Doctoral Training Program from Natural Environment Research Council Samsung Science and Technology Foundation 10.13039/501100014364 SSTF-BA1501-52
==== Body
pmcIntroduction

Animal reproduction involves the production of gametes with distinct characteristics and patterns of inheritance. Male and female gametes differ not only in size but also in the transmission of cytoplasmic elements. Typically, cytoplasmic components such as mitochondria and bacterial endosymbionts are maternally inherited (but see Hoeh et al. 1991; Zouros et al. 1992; Hurst 1993; Zouros et al. 1994). In contrast, the nuclear genome is usually inherited equally from both parents, ensuring a balanced contribution of genetic material. However, exceptions to this symmetrical nuclear inheritance exist across various species, ranging from maternal-only to paternal-only inheritance of the nuclear genome (for review, see Ross et al. 2022). These variations in inheritance patterns for cytoplasmic elements and the nuclear genome highlight the dynamic nature of reproductive strategies in animals. Understanding these exceptions contributes to our broader understanding of reproductive biology and the diverse mechanisms employed by organisms to propagate their genetic material.

One such variation is the asymmetric transmission of sex chromosomes. In XX:XY sex determination systems, for instance, half of the gametes produced by the male have 1 X chromosome and the other half have 1 Y chromosome. The transmission is asymmetric because only the male can transmit the Y chromosome to the next generation, and females can transmit either of the 2 X chromosomes. Although the expected XX:XY ratio in the offspring for this type of sex determination is 1:1, selfish genetic elements in one of the sex chromosomes may drive a non-Mendelian transmission, resulting in a bias toward female or male offspring (Burt and Trivers 2006).

Similarly, XX:XO sex-determining systems may also result in non-Mendelian transmission of the X chromosome, resulting in sex ratio bias after crosses. In this type of sex-determining system, the male is heterogametic and therefore is expected to produce an equal number of X-bearing and nullo-X gametes. In some organisms, however, the nullo-X sperm are not produced by males, resulting in gametes that only carry an X chromosome (Blackman 1985; Wilson et al. 1997; Tandonnet et al. 2018). Thus, the following generation is biased toward XX animals.

Nematodes are known to have diverse sex determination systems (Pires-Dasilva 2007), which makes them well suited for studying sex ratios due to their short life cycles, prolific reproduction, and ease of husbandry. In some nematode clades, crosses between XX and XO individuals result in highly biased XX offspring (for review, see Van Goor et al. 2021). In only a few cases, the mechanisms underlying this sex ratio bias are known. In nematodes that include Auanema and Strongyloides, for instance, the male-producing nullo-X spermatid is obligately eliminated during spermatogenesis (Shakes et al. 2011; Winter et al. 2017; Tandonnet et al. 2018; Dulovic et al. 2022).

Nematodes of the genus Auanema are trioecious, consisting of XX hermaphrodites, XX females, and XO males. In hermaphrodites and males, an asymmetric division occurs during spermatogenesis, resulting in organelles partitioning to different sides of the dividing cell (Fig. 1; Shakes et al. 2011; Winter et al. 2017; Al-Yazeedi et al. 2022). During spermatogenesis in A. freiburgense, components essential for sperm function, such as mitochondria and organelles containing cytoskeletal proteins necessary for sperm motility, cosegregate with the X chromosome (Winter et al. 2017). Meanwhile, the sister cell lacking an X chromosome and consisting only of autosomes undergoes differentiation into a residual body that is subsequently eliminated together with organelles such as the Golgi complex and endoplasmic reticulum (Shakes et al. 2011; Winter et al. 2017; Al-Yazeedi et al. 2022). Consequently, crosses between males and XX females result in the production of mostly XX offspring (Félix 2004; Shakes et al. 2011; Tandonnet et al. 2018, 2022).

Fig. 1. Models of predominant meiosis in Auanema. a) In wild-type male spermatogenesis, the chromatids of the X chromosome separate in meiosis I. In anaphase II, the cell that receives the X chromosome also inherits mitochondria and the fibroid bodies containing the major sperm protein (MSP). The cell without the X chromosome (nullo-X cell) inherits nonsperm materials and is discarded as a polar body (not to scale). b) Auanema females have canonical meiosis, where the X chromosomes undergo recombination, and there is the production of an X-bearing oocyte and 3 polar bodies. Therefore, most offspring from crosses between females and males are XX. However, if nondisjunction occurs during female oogenesis, it can result in a nullo-X oocyte and male offspring (not pictured in the diagram). c) During hermaphrodite spermatogenesis, diplo-X sperm are produced; the 2 homologous X chromosomes are segregated into 1 daughter cell during anaphase II. d) In hermaphrodite oogenesis, the resulting oocyte does not harbor an X chromosome. We note that these diagrams depict the predominant X segregation patterns. Rare males are produced by female–male crosses and by selfing hermaphrodites. In A. rhodense, males from female–male crosses are derived from the fertilization of rare nullo-X oocytes with X-bearing sperm. Males from selfing hermaphrodites may originate from viable haplo-X sperm fertilizing a nullo-X oocyte. The mechanisms controlling these meiosis patterns are still unknown.

Based on the timing of the X chromosome cosegregation relative to the organelles and the symmetric segregation of organelles in masculinized XX mutants, it has been suggested that the X chromosome may serve as a polarizing signal for this asymmetric cell division (Winter et al. 2017; Al-Yazeedi et al. 2022). In this study, we report the use of recombinant inbred advanced intercross lines (RIAILs), whole-genome sequencing, and quantitative trait locus (QTL) mapping to identify the genetic components controlling asymmetric segregation of organelles in male meiosis in Auanema freiburgense (Kanzaki et al. 2017; Winter et al. 2017). Some of the RIAILs displayed a high rate of male production after outcrossing males. We show that males in these lines produce viable nullo-X sperm, and that this new phenotype maps to a region in the X chromosome. These findings are consistent with the hypothesis that loci in the X chromosome are acting as a polarizing signal during sperm formation, although loci in other chromosomes may also be involved.

Materials and methods

Strains and their maintenance

Auanema freiburgense (previously known as Rhabditis sp. SB372, or Auanema freiburgensis) was first isolated in Freiburg, Germany, by Prof. Walter Sudhaus from a horse dung pile (Kanzaki et al. 2017; Sudhaus 2023). The A. freiburgense SB372 strain underwent 11 generations of bottlenecking (expansion from a single-selfing hermaphrodite) to produce the inbred strain APS7 (Adams et al. 2022). Auanema freiburgense JU1782 strain was isolated from a rotting Petasites stem sampled in Ivry (Val-de-Marne, France; Robles et al. 2021) by Marie-Anne Felix. JU1782 underwent bottlenecking for 10 generations and was renamed inbred strain APS14. Nematodes were maintained on NGM plates seeded with Escherichia coli OP50-1 at 20°C.

Sexual morph identification and female–male cross

In uncrowded conditions, A. freiburgense hermaphrodites produce mostly female and male progeny (Zuco et al. 2018; Robles et al. 2020). Auanema freiburgense dauer larvae invariably develop into self-fertilizing hermaphrodite adults (Kanzaki et al. 2017; Zuco et al. 2018). Therefore, to isolate female and male sexual morphs for crosses, dauer larvae were incubated on NGM OP50-1 plates at a low density (3 dauers per 2-cm-diameter OP50-1 bacterial lawn) until they reached adulthood and began laying eggs (approximately 48 h after collection; Adams et al. 2022). Approximately 36 h after the start of egg laying, L2 stage female larvae were differentiated from males by tail morphology and moved to fresh female-only plates to prevent fertilization. After 24 h, virgin females reached adulthood and were used in crosses. Young adult males were isolated from the original low-density plates.

The crosses were conducted with a ratio of 1 female to 1 male. The male was subsequently removed after 24 h, and the mother was transferred to a new plate. The male and nonmale progeny were counted 48–72 h after egg laying, and the counts were combined for all 3 plates for each cross.

Generation of A. freiburgense RIAILs

One hundred A. freiburgense RIAILs were generated by crossing APS7 males with APS14 females. Hybrid F1 progeny resulting from the cross was left to mate (or self-reproduce) in a large mating pool. F1 self-reproducing hermaphrodites and females, each mated with 1 or more males, were isolated to establish the lines. Three to five F2 females were picked from each line and crossed with 2 males from a different line in an inbreeding avoidance scheme. Intercrosses between lines from F2 to F7 were established in such a way that every 2 lines were only crossed once to maximize haplotype breakpoints. After 7 generations of interline crosses, lines were inbred by single-worm descent for 10 generations to bring alleles into homozygosity at most loci. Once the 17th generation was reached, the 100 lines were maintained by sampling animals from a crowded plate to a fresh new plate once a month.

Crossing males from A. freiburgense RIAILs with wild-type APS7 females

In each cross, a male from each RIAIL was crossed with a wild-type APS7 female for ∼24 h. After the cross, males were removed. The fertilized female was moved to a new plate every day until it stopped laying eggs, to synchronize the growth of the F1 progeny. Ratios of male-to-female progeny were scored for each cross. Males from all lines were crossed except line numbers 17, 48, 79, and 107.

DNA extraction from A. freiburgense RIAILs and sequencing

For each line, the nematodes were harvested with water from 5 NGM plates (10 cm in diameter). They were collected in a conical tube and washed 2–3 times with water. After each wash, nematodes were allowed to settle naturally to the bottom of the tube rather than by centrifugation. Five hundred microliters of lysis buffer [100 mM Tris (pH 8.5), 100 mM NaCl, 50 mM EDTA, 1% SDS, and 1% beta-mercaptoethanol] was added, and tubes were frozen at −80°C overnight.

Three cycles of thawing and freezing were performed before adding 2.5 μL of proteinase K (20 mg/mL) to each tube, followed by incubation at 65°C for 3–4 h. DNA extraction was performed with the Gentra Puregene Core Kit (Qiagen) following the manufacturer's instructions.

Sequencing libraries were generated using TruSeq DNA nanogel free at the GenePool facility at the University of Edinburgh. Sequencing libraries were sequenced using the Illumina HiSeq platform to generate 150-bp paired-end reads with an insert size estimation of 350 bp.

DNA extraction from A. freiburgense for PacBio, Illumina mate-pair, and paired-end sequencing

For PacBio long-read sequencing, DNA was extracted from plates containing A. freiburgense APS7 at various stages. To lyse worms, we used 6 mL of 55°C Cell Lysis Solution from The Gentra Puregene® Cell and Tissue Kit (Qiagen) with 0.1 mg/mL proteinase K and 1% β-mercaptoethanol. The solution was directly poured into a 200-µL worm pellet, and worms were lysed in the solution at 55°C for 8 h with occasional inverting. Polynucleotides were purified from the mixture by using phenol–chloroform–isoamyl alcohol (25:24:1 v/v) DNA extraction and ethanol precipitation methods coupled with phase-lock gel to minimize pipetting and DNA shearing. Purified polynucleotides were redissolved in the TE buffer and treated with 10 μg/mL RNase for 2 h. DNA was extracted by the same DNA purification procedure and dissolved in 10 mM Tris-HCl (pH 8.0). Macrogen (South Korea, https://www.macrogen.com/en/main) performed PacBio library preparation and sequencing on the Sequel platform with the continuous long-read sequencing mode.

DNA extraction for Illumina mate-pair and paired-end sequencing of APS7 was performed as described for the RIAIL sequencing, using DNA from dauer larvae collected from 200 plates, as previously described (Pires-Dasilva 2013).

Genome assembly

The assembly of the PacBio data was performed with the Canu software (version 1.6; Koren et al. 2017), using a minimum read length (minReadLength) of 3 kb and setting the corrected error rate (correctedErrorRate) to 0.030, as this resulted in the most contiguous assembly. The estimated genome size (genomeSize) was set to 55 Mb.

Bacterial contamination was identified and removed through BLASTn (version 2.7.1; Camacho et al. 2009) alignments to a contaminant database that included 3,000 bacterial genomes downloaded from the European Nucleotide Archive (ENA) on March 30, 2018. The following BLAST parameters were used: “-task megablast - evalue 1e-06 -outfmt 6 -perc_identity 50.”

We first polished our preliminary draft genome using Quiver with PacBio raw reads aligned to the assembly using pbalign (version 0.3.1). We then further polished the assembly using Pilon with all available Illumina short reads. Illumina short reads were first preprocessed using Skewer (parameters “-n -Q 20 -l 51”). Trimmed reads were then aligned against the assembly using BWA. These alignments were then used to fix base-level inconsistencies between the Illumina read and the genome assembly using Pilon (“--changes --fix bases --chunksize 8000000 --diploid; Walker et al. 2014).

A list of all the programs used in all bioinformatics analyses, their versions, and parameters is included in Supplementary Table 9.

RIAIL genotyping

The assembled genome, polished with Illumina reads, was used as a reference genome to genotype the 100 RIAILs and the APS14 parent samples. The quality of individual samples’ raw reads was checked using FastQC to ensure sequencing adaptors were removed. Each line/strain paired-end DNA sample was aligned against the long-read genome assembly using the BWA aligner (Li and Durbin 2009). Sequence Alignment Map (SAM) files were converted into Binary Alignment Map (BAM) files, then sorted according to the position in the reference genome (Li et al. 2009). Aligned reads in every BAM file were assigned a new read-group tag using the Picard AddOrReplaceReadGroups tool to make it compatible with the Genome Analysis Toolkit (GATK) pipeline (Mckenna et al. 2010). Variants from every alignment bam file with a new read-group tag were identified against the APS7 reference genome using genome analysis toolkit GATK tools, producing a single variant calling file (VCF) containing the variants of each sample against the reference APS7 genome. Low-quality variants were filtered with VariantFiltration provided by GATK tools using default parameters (Mckenna et al. 2010; Depristo et al. 2011; Van Der Auwera et al. 2013; Poplin et al. 2018). Then variants were further filtered using vcftools keeping only variants with a minimum depth of 10 and a minimum genotype quality of 30, removing variants that are missing in more than 25% of RIAILs, removing all the indels, and only keeping biallelic single nucleotide polymorphism (SNP) variants (Danecek et al. 2011). The final 274,394 SNPs were used as markers to construct a genetic linkage map.

Construction of a high-density genetic linkage map

The genetic linkage map was constructed using R/qtl and ASMap R packages (Broman et al. 2003; Taylor and Butler 2017). R/qtl was used to process markers for the preconstruction of the genetic linkage map and ASMap for the genetic linkage map. Markers were imported into R using rqtl as RIL data, crosstype = “riself,” expecting no heterozygous markers in the data set, and heterozygous markers were assigned a missing value. A total of 10,326 markers either missing the genotype of one of the parents or found heterozygous were omitted from the data set, leaving 264,068 markers from 100 samples on 40 scaffolds. Using the rqtl function “drop.markers,” the marker set was further refined by (1) removing markers that were absent in 90% of the data, allowing only 10% of missing values; (2) removing markers that shared similar genotypes to reduce redundancy in the data set; and (3) removing markers with abnormal segregation distortion greater than the genome-wide alpha level of 0.05/number of markers (P < 0.05). After that, the ASMap R package was used to construct a genetic linkage map with the remaining 14,955 markers from 100 samples on 25 scaffolds. The genetic map was constructed using the “mstmap.cross” function provided by the ASMap package using a P-value of 1e−11. The genetic distance between markers was computed, and markers were linked and organized into 7 main linkage groups. An initial genetic linkage map of 7 linkage groups was constructed using 14,884 markers from 100 samples on 22 scaffolds, representing 93.7% of the genome.

To improve X chromosome assembly, another genetic linkage map was constructed specifically for the X chromosome by including all the heterozygous and homozygous markers and removing all autosomal scaffolds. Linkage group 7 from the initial genetic linkage map was identified using synteny mapping as the X chromosome in A. freiburgense. To construct an independent genetic linkage map for the X chromosome, markers were imported into rqtl as “f2 intercross” to include all the heterozygous and homozygous markers. In total, 264,068 markers from 40 scaffolds were imported. Using the rqtl function “drop.markers,” we removed markers that were not present in 75% of the data, allowing only 25% of missing values, and redundancy in the data set was reduced by removing markers that share similar genotypes. Then, markers belonging to autosomal linkage groups (L.1–L.6) were removed, leaving only 21 scaffolds with 3,128 markers from X chromosome scaffolds and unplaced scaffolds. Before linkage map construction, using the ASMap “pullCross” function, markers with more than 10% missing values and those with a high segregation distortion (segregation ratio less than 1:98:1) were removed to be pushed back into the map after construction. The X chromosome genetic map was constructed using the “mstmap.cross” function provided by the ASMap package using a P-value of 1e−12. Using the “pushCross” function, most markers with high segregation distortion were pushed back to the constructed map using a segregation ratio threshold of 0.5:99:0.5. Small linkage groups were subsetted, leaving only the main linkage group representing the X chromosome with 2,632 markers. The new X chromosome genetic linkage map was added to the previously constructed genetic linkage map, replacing linkage group 7 (L.7). The final genetic linkage map of 7 linkage groups contained 16,792 markers from 29 scaffolds and represents 97% of the genome (Datasheet 3).

Synteny analysis and identification of ancestral chromosomal elements (Nigon elements)

Scaffolds were anchored to the genetic map using ALLMAPS software, producing a chromosomal-scale assembly where scaffolds were ordered and oriented into their respective positions within each linkage group (Tang et al. 2015). The chromosomal assembly was aligned to the Caenorhabditis elegans genome using MUMmer with default parameters, and macrosynteny patterns were visualized using Circos (Kurtz et al. 2004; Krzywinski et al. 2009).

To assess the completeness of the genome assembly, a BUSCO analysis was conducted in the genome mode using Augustus gene discovery against the latest nematode database (nematoda_odb10; Supplementary Tables 3 and 8). To examine the patterns of Nigon elements in each chromosomal-scale scaffold, we used the program vis-ALG (https://github.com/pgonzale60/vis_ALG; Gonzalez De La Rosa et al. 2021). Briefly, we used the location of the BUSCO genes identified in A. freiburgense and the association between BUSCO genes and Nigon elements previously determined to paint the chromosomal-scale scaffolds according to the Nigon elements.

X chromosome genotyping and calculation of nullo-X ratio

To follow the X chromosome inheritance, we used an X-linked polymorphic marker(X634) where a HindIII restriction site (AAGCTT) is present in APS7 but not in the APS14 strain (AAACTT). X chromosome genotyping was conducted by PCR amplification of a region spanning the SNP followed by digestion of the product using HindIII. Genomic DNA was extracted from individual males using a modified single-worm PCR method. A single male was frozen in 20 µL of 1× PCR buffer at −80°C for a minimum of 24 h. After thawing, the tissue was lysed and genomic DNA was released by the addition of 0.5 µL of proteinase K (20 mg/mL) and incubation at 65°C for 60 min, followed by 95°C for 15 min, to inactivate the enzyme. Samples were frozen at −20°C for at least 24 h before use in PCR. Each PCR reaction was conducted with 2 µL of DNA, 10 µL of GoTaq Green Master Mix (Promega), 10 µM of the forward primer (UW634_F 5′-AGGGACACGATTGCCTTCTG-3′), and 10 µM of the reverse primer (UW635_R 5′-AATGCCGCGGAGGTCTTTAA-3′) in a final volume of 20 µL. The following cycling conditions were applied: 94°C for 5 min, followed by 30 cycles of 94°C for 15 s, 55°C for 30 s, and 72°C for 1 min. The PCR products were digested by direct addition of 0.5 µL of HindIII (Promega) and incubation for 1 h at 37°C. The genotype of each sample was determined by the agarose gel electrophoresis of the digested sample. The APS7 X allele gave 2 fragments (329 and 234 bp), and the APS14 allele remained undigested (563 bp). Note that the X634 marker region is outside the QTL region (coordinates 438,785–439,347).

To determine the percentage of nullo-X sperm that contributed to the generation of sons, the paternal X chromosome was genotyped using the X-linked marker X634 (Supplementary Table 1). To determine the percentage of viable nullo-X sperm for each strain, we multiplied the percentage of sons resulting from crosses by the proportion of sons that inherited the maternal X chromosomes from the total number of genotyped males (Supplementary Tables 1 and 2). This is based on the understanding that only sons that have a maternal X are derived from the fertilization of a nullo-X sperm.

Structural annotation of the genome

A comprehensive repeat library was produced for the genome assembly using multiple repeat-finding programs. The pipeline was based on the TransposableELMT wrapper script (https://github.com/PlantDr430/TransposableELMT). Repeats were identified using RepeatModeler v2.0.1 (Smit and Hubley 2008-2015), TransposonPSI v08222010 (Haas 2010), LTRfinder v1.0.7 (Xu and Wang 2007), and LTRharvest (Ellinghaus et al. 2008; implemented in GenomeTools v1.6.1; Gremme et al. 2013). To limit the identification of false positives, the LTRharvest output was postprocessed with LTRdigest (Steinbiss et al. 2009; implemented in GenomeTools v1.6.1; Gremme et al. 2013). The resulting libraries were combined, classified using RepeatClassifier v2.0.1 (part of the RepeatModeler v2.0.1 package (Xu and Wang 2007), and redundancy removed using USEARCH (Edgar 2010) based on 90% similarity. The nonredundant custom library was used to soft-mask repeat regions in the assembly with RepeatMasker v4.1.0 (Smit and Hubley 2013-2015).

Gene predictions were made using the ab initio and evidence-driven gene predictors GeneMark-ES (Lomsadze et al. 2005), SNAP (implemented in Maker2; Korf 2004), Maker2 (Holt and Yandell 2011), and Augustus (Stanke and Waack 2003; trained with BUSCO v5.1.3; Seppey et al. 2019).

An A. freiburgense transcriptome assembled using Trinity (Grabherr et al. 2011), the C. elegans protein database (UP000001940_6239), and the UniProt/Swiss-Prot database (uniprot_sprot.fasta) were used as evidence-based inputs for the first round of Maker2. The resulting output was then used to train Augustus and SNAP. The unmasked genome was used as input into GeneMark-ES. Finally, the outputs from the first round of Maker2, Augustus, SNAP, and GeneMark were used as inputs in Maker2. The gene predictions from the second round of Maker2 were used for our analyses. The annotation was later lifted over to the final chromosomal-scale assembly using “flo” using chain files (Pracana et al. 2017).

Functional annotation of the genome

The functional annotation of the genome was conducted using Blast2Go (Götz et al. 2008). A local Blast search of the A. freiburgense predicted proteins was conducted using a searchable database of the Swiss-Prot proteins (swissprot.gz downloaded from NCBI in July 2021) prepared within the Blast2Go software. The same predicted proteins were also analyzed with InterProScan (within the Blast2Go software), and the annotations were merged to give the final functional annotation.

The completeness and quality of the predicted protein-coding genes were assessed using BUSCO on the transcript data set in genome mode with gene discovery via Augustus against the nematode database (nematoda_odb10).

Mitochondrial genome

The mitochondrial scaffold was identified by BLASTn (Version 2.9.0+; Camacho et al. 2009), using the Auanema rhodense mitochondrial genome as a reference database. Only 1 scaffold was identified as the mitochondrion. This scaffold was analyzed and annotated using MITOS (Donath et al. 2019), which revealed a large duplication, probably due to a misassembly due to the mitochondrial genome being circular. We used the pairwise local aligner Water (Smith and Waterman 1981) to identify precisely the junctions of the duplicated region and removed them using the subseq function of seqkit (version 0.16.1; Shen et al. 2016). The curated mitochondrial genome was reintegrated into the assembly to replace the uncurated one. The annotations associated with the uncurated mitochondrial genome were removed and replaced by the MITOS annotations.

Pooling DNA from lines with similar phenotypes into discrete pools and sequencing

Equal amounts of DNA from 10 RIAIL lines with the same phenotype were mixed to create pools of DNA from high-male lines (HM-pools) and low-male lines (LM-pools). Two pools for each category were prepared, each containing 1.5 µg of DNA. The RIAIL lines used for the first and second HM-Pools were 23, 24, 26, 28, 29, 30, 33, 45, 57, and 61 and 35, 42, 46, 55, 56, 59, 63, 65, 72, and 97, respectively. The RIAIL lines used for the first and second LM-Pools were 2, 5, 9, 13, 14, 17, 20, 54, 74, and 77 and 6, 8, 15, 18, 25, 27, 37, 40, 69, and 95, respectively. The quality of DNA was examined by running 1 μL of DNA on 1.8% agarose gel, and concentration was measured using a Qubit fluorometer. Five DNA samples—2 HM-pools, 2 LM-pools, and DNA from the APS14 maternal strain—were sequenced. Sequencing libraries were generated using TruSeq DNA nanogel free at the GenePool facility at the University of Edinburgh on an Illumina HiSeq platform to generate 150-bp paired-end reads with an insert size estimation of 350 bp.

Variant calling for APS14 and next-generation sequencing bulk segregant analysis

The quality of the reads was assessed using FastQC software to get an overview of the reads’ quality (Patel and Jain 2012). Paired-end reads were cleaned using Skewer using the following parameters: -n -Q 20 -l 51 -t 32 -m pe (Jiang et al. 2014). The quality of the reads was reexamined after cleaning using FastQC software. Then, the BWA program was used to align short reads from each pool and APS14 strain separately to the APS7 reference genome (Li and Durbin 2009). Alignment files in SAM format were converted to BAM files using the Samtools “view” command (Li et al. 2009). BAM files from the HM-pools and the LM-pools were merged to produce a single BAM for each phenotype and were subsequently sorted using samtools (Li et al. 2009). Variants from sorted BAM files of both pool samples and the APS14 strain were called individually using “HaplotypeCaller” from GATK (Van Der Auwera et al. 2013). Low-quality variants were filtered out with “VariantFiltration” provided by GATK tools using default parameters (McKenna et al. 2010; DePristo et al. 2011; Van der Auwera et al. 2013; Poplin et al. 2018). Variants were filtered further using bcftools based on genotype quality ≥ 30 and depth ≥ 10 while keeping homozygous sites only (Danecek et al. 2011). SnpEff (v5.0.1; Danecek et al. 2011), SnpEff (v5.0.1; Cingolani, Platts, et al. 2012), and SnpSift (v4.3.1; Cingolani, Patel, et al. 2012) were used to annotate and predict the effects of the genetic variants between the APS14 and APS7 reference genomes. Large structural variants between the APS14 and APS7 genomes were identified using the Parliament2 pipeline (Zarate et al. 2020). Structural variants were called using the following software: Breakdancer (Chen et al. 2009), CNVnator, (Abyzov et al. 2011), Manta (Chen et al. 2016), Lumpy (Layer et al. 2014), and DELLY (Rausch et al. 2012). To avoid false positive calls, we only considered structural variants that were called with more than 1 software.

Pools’ VCFs were merged into 1 file using the GATK “GenotypeGVCFs” command (Van Der Auwera et al. 2013). Joint genotyping using “GenotypeGVCFs” combined all SNP and indel records from both pools to produce the correct genotype likelihood, outputting a single combined VCF (Brouard et al. 2019). Variants were filtered using vcftools (Danecek et al. 2011), removing all sites with missing values and indels and keeping only biallelic sites. In total, 334,230 SNPs were kept. The joint VCF was converted to a table using the “VariantToTable” command provided by the GATK (Van Der Auwera et al. 2013). Regions displaying differences between the HM-pools and LM-pools were identified using the R package “QTLseqr” (Mansfeld and Grumet 2018). QTL-seq analysis for next-generation sequencing bulk segregant analysis (NGS-BSA) was used to calculate the allele frequency difference (SNP-index) from the allele depth at each individual SNP (Takagi et al. 2013). Candidate regions were identified by setting a sliding window size to 1 Mb, and the number of SNPs was counted in that window. Within each sliding window, a tricube-smoothed delta-SNP-index was calculated by constant local regression. A simulation was performed where the delta-SNP-index per bulk was calculated and simulated over 1,000 replications based on the RIAIL F2 population with a bulk size of 20. Confidence intervals at 95 and 99% were estimated using the quantile from the simulation. An alternative approach, G-statistics, was used to identify significant QTLs from BSA (Magwene et al. 2011). G-statistics was calculated genome-wide, and a tricube-smoothed G-statistics (G′) was predicted in a sliding window of 1 Mb. P-values were estimated and adjusted (Benjamini–Hochberg method), and negative log10 was calculated from G′. Candidate regions were identified using a genome-wide false discovery rate (FDR) of 0.1 (Benjamini and Hochberg 1995).

QTL analysis, Fst differentiation analysis, and identification of potential candidate genes

The QTL analysis was conducted with rqtl using all the markers in all chromosomes (Broman et al. 2003; Zuo et al. 2019). SNPs identified in all the RIAILs during the construction of the genetic linkage were lifted to a new coordinate based on the latest genome assembly using LiftOvervcf. In total, 261,860 markers were imported into rqtl to perform a qtl scan, allowing only 10% of missing values where markers have to be present in at least 90% of RIAIL individuals. All the markers with high segregation distortions were retained. Marker segregation distortion is a natural phenomenon, and the inclusion of distorted markers in QTL mapping increases the power of detecting a QTL (Lyttle 1991; Zuo et al. 2019). The interval mapping genotype probability was calculated using the Kosambi map function with a 1 cM step size and an error probability of 0.001. A whole-genome scan on all RIAILs was performed with a single QTL model using Haley–Knott regression, with the percentage of male progeny after an outcross for each RIAIL as a phenotype (Haley and Knott 1992). The result from the first scan was permuted 1,000 times to obtain a genome-wide log of odds (LOD) score significance threshold (Churchill and Doerge 1994). Genome-wide thresholds were obtained at P < 0.05, corresponding to a LOD score of 2.95.

To identify if there is an interaction between the identified QTL and the rest of the genome, we performed a scan using the QTL genotype as a covariate. We initially performed a scan with a single QTL model using Haley–Knott regression with QTL genotype as an additive covariate to detect a putative QTL with the same effect in both QTL genotypes. Another single scan was performed using Haley–Knott regression with QTL genotype as an interactive covariate so that the QTL is allowed to be different in the 2 different genotypes. To test for an interaction between the QTL genotype and the rest of the genome, we obtained the difference between the LOD score with QTL genotype as an interactive covariate and the LOD score with QTL genotypes as an additive covariate. A separate permutation was performed 1,000 times with the QTL genotype as an additive covariate and with the QTL genotype as an interactive covariate using Haley–Knott regression to obtain a significant threshold. Again, the difference in permutation tests concerns the interaction between the QTL genotype and the rest of the genome. We obtained a genome-wide LOD threshold P < 0.05, corresponding to a LOD score of 1.85, where no significant peaks were detected (Supplementary Fig. 9).

Scanning for the likelihood of interacting QTLs using a 2D genome scan for all pairwise combinations of intervals was not possible due to the large number of markers (Dupuis et al. 1995). To obtain the estimated effect of the QTL, a qtl object was created containing the identified QTL from the single scan using the rqtl function (makeqtl). Then, the QTL was added to an additive model using the function (fitqtl) with the formula y ∼ QTL.

The window size Fst differentiation analysis between LM and HM lines was performed using vcftools with a window size of 5,000 bp and a step size of 1,000 pb. SNP-level wcFst and pFst were calculated using vcflib.

Results

Auanema freiburgense RIAILs exhibit a transgressive phenotype

Auanema freiburgense is a species with 3 sexual morphs (XX females, XX hermaphrodites, and XO males; Kanzaki et al. 2017). Previous research has shown that when an A. freiburgense female crosses with a male, the resulting offspring have a skewed sex ratio against males (Kanzaki et al. 2017; Winter et al. 2017). An XX sex ratio bias occurs because the X-bearing sperm of the male is viable, whereas the nullo-X spermatid is not (Winter et al. 2017).

Here, we tested 2 A. freiburgense inbred strains, APS7 and APS14 (Adams et al. 2022), which produced 18% or fewer sons after outcrossing (Datasheet 1). The high ratio of XX offspring confirms previous studies with other strains of the same species (Kanzaki et al. 2017; Winter et al. 2017; Tandonnet et al. 2022). To investigate sex determination in A. freiburgense, we created a genetic linkage map using 100 RIAILs (Fig. 2). Approximately 60% of these recombinant lines exhibited a transgressive phenotype, indicating the presence of traits not observed in either parental strain. When males from these hybrid lines crossed with APS7 females, they produced a higher percentage of sons (high-male RIAILs or HM-RIAILs) than the parental strains (Fig. 3; Datasheet 2).

Fig. 2. Construction of A. freiburgense RIAILs. The breeding scheme involved 7 generations of crosses and 10 generations of bottlenecking by the propagation of individuals derived from the self-fertilization of a single hermaphrodite. The somatic cells are represented at various stages in the construction of RIAILs, showcasing mitochondrial and hypothetical chromosomal backgrounds depicted in different colors for each strain. Mitochondrial representation is color-coded to indicate the maternal origin of the mitochondria. RIAILs were created by crossing an APS14 female with an APS7 male. F1s resulting from the cross were left to mate in a large mating pool. F1 hermaphrodites and females mated with male siblings were isolated on 100 individual plates to produce F2s, each with unique recombination patterns. From the F2 to the F7 generation, lines were intercrossed in an inbreeding avoidance scheme. In this scheme, each line was crossed with another line only once to maximize haplotype breakpoints. Then, each line was selfed by single-worm descent for 10 generations to increase genome homozygosity. By the F17 generation, each line had a unique mixed genome from both parental strains and was homozygous at most loci.

Fig. 3. Percentage of male offspring from RIAIL crosses with wild-type APS7 female. In an outcross with a wild-type female, males from most RIAILs produce more sons than the APS7 and APS14 strains. In total, there are 57 HM lines and 40 LM lines, which include 16 lines where no male offspring from crosses were observed. Additionally, 3 lines were not phenotyped. *Indicates that crossing data is not available.

HM lines produce more viable nullo-X sperm than the parental and LM lines

To monitor the inheritance of the X chromosome after a cross, we used a polymorphic marker on the X chromosome. We found that females and males from APS7 and APS14 produce rare viable nullo-X oocytes and rare viable nullo-X sperm, respectively (Fig. 4; Supplementary Table 2). For instance, in reciprocal crosses between APS7 and APS14, sometimes the sons inherited the X from the mother, indicating fertilization between an X-bearing oocyte and nullo-X sperm. At other times, the sons inherited from the father, indicating fertilization between a nullo-X oocyte and X-bearing sperm (Fig. 4). Overall, the percentage of viable nullo-X sperm produced by APS7 and APS14 males is relatively low (<20%; Supplementary Table 2; Datasheet 1).

Fig. 4. Punnett’s square diagrams depicting the observed progeny ratios between different A. freiburgense strains. The percentages of each progeny type were calculated by (1) determining the ratio of average male offspring (Datasheet 1) and (2) determining the inheritance pattern of the X chromosome in male offspring through genotyping, distinguishing whether the X chromosome originated from the APS7 or APS14 strain (Supplementary Tables 1 and 2). Crosses were performed between APS7 and APS14 males a), between APS7 females and line 19 (HM) males b), and between APS7 females and line 38 (LM) males c). The results of additional crosses are in Datasheet 1 and Supplementary Tables 1 and 2. Oocytes may undergo X chromosome nondisjunction (ND) to generate viable nullo-X female gametes, and spermatocytes in males may undergo symmetric cell division (SCD) to generate viable nullo-X sperm. Note that the OO progeny cannot be counted, as these cases are most likely inviable and not observable.

Given the observation that the HM lines described above produce a large percentage of sons, we hypothesized that males in those strains were producing more nullo-X sperm than the parental lines. To investigate this, we selected 2 HM lines and 2 LM lines to study the inheritance patterns of the X chromosome in male descendants. These lines were chosen at random, with the requirement that LM lines must yield a minimum number of male offspring to allow for a thorough analysis of X chromosome inheritance. Males from these lines carried the APS14 X genotyping marker and were crossed with an APS7 mother to enable tracking of the X chromosome. Males from the HM lines (lines 12 and 19) produced around 50% of sons, whereas males from the LM lines (lines 2 and 38) generated less than 5% of sons after outcrossing (Fig. 4; Supplementary Table 2). Sons from the LM males inherited the paternal X chromosome in all samples tested. These results indicate that these lines produce viable X-bearing sperm and few or no nullo-X sperm. In contrast, almost all the sons from the HM males inherited the maternal X chromosome, indicating that ∼50% of the sperm in these lines are composed of viable nullo-X sperm (Fig. 4; Supplementary Table 2).

X chromosome independent high-density genetic linkage map improved X chromosome assembly

To map the locus involved in the production of sons after crossing, the genome of A. freiburgense (strain APS7) was sequenced using long reads from Pacific Biosciences (PacBio) and short Illumina paired-end and mate-pair data. The initial genome assembly with these sequencing data resulted in 75 scaffolds (Afr-genome-v1). To improve the genome assembly, the scaffolds were anchored and ordered onto a genetic linkage map. To create this map, we used sequencing data from 100 RIAILs, the parental lines APS7 and APS14, and used Afr-genome-v1 as a reference.

In the related species A. rhodense, the X chromosomes do not recombine in hermaphrodites (Tandonnet et al. 2018). Therefore, we expected that the bottlenecking of the RIAILs by selfing single hermaphrodites would still result in an X chromosome containing significant heterozygosity. To accommodate this unique biology, we employed 3 different approaches to building the genetic linkage map, with parameter settings as follows (see Materials and methods):

Using only homozygous markers for all chromosomes (Supplementary Figs. 1a and 2). The resulting map consisted of 7 linkage groups, with the X linkage group containing fewer markers than the autosomes (Supplementary Table 3).

Using both homozygous and heterozygous markers for all chromosomes. However, the markers clustered in several groups (more than 7; data not shown).

A map was constructed independently for the autosomes and another one for the X chromosome (Table 1; Supplementary Figs. 1b and 2). For this latest approach, we used only homozygous markers to group the autosomal markers and both homozygous and heterozygous markers to group the X chromosome markers. Using this approach, the final length of the X chromosome doubled compared to the first approach (Table 1; Supplementary Fig. 1).

Table 1. Characteristics of the A. freiburgense genome and genetic linkage map.

LG	N markers	AA%	BB%	Length (cM)	Length (Mb)	N genes	% repeat DNA	
L.1	1,079	23.3	76.7	41.0	7.2	1,112	19.7	
L.2	4,404	48.6	51.4	267.4	13.0	1,964	20.4	
L.3	3,337	59.1	40.9	312.7	11.4	1,833	13.6	
L.4	1,638	25.7	74.3	54.9	5.3	896	14.9	
L.5	2,290	27.7	72.3	95.2	8.0	1,442	10.9	
L.6	1,412	28.4	71.6	82.0	5.0	979	8.8	
L.7	2,632	14.5	84.4	186.0	3.3	496	35.5	
Total	16,792	59.0	34.0	10,038.9	53.5	8,759	16.7	
Heterozygous markers are only included for the X chromosome, comprising 1.1% of the X chromosome markers and 0.2% in the overall linkage map.

The final genetic linkage map consists of 7 linkage groups with 16,792 markers from 29 scaffolds, accounting for 97% of the sequenced genome (Fig. 5; Supplementary Fig. 2; Datasheet 3). The linkage map’s overall percentage of missing markers is 6% (Fig. 5; Supplementary Figs. 3 and 4). Around 3% of sequenced DNA (1,689,119 bp in 47 scaffolds) was left unplaced, as it did not cluster with any linkage group (Supplementary Table 3); 15 kb of those unplaced sequences correspond to the mitochondrial genome. When the coding genes from this integrated map-scaffold genome were compared to those of C. elegans, a high degree of synteny was observed between L.7 and the C. elegans X chromosome (Supplementary Fig. 5).

Fig. 5. Parental background of the RIAILs, genotype distribution, and recombination across the genetic linkage map. Genotypes are clustered into 7 linkage groups, with each group appearing in a separate panel on the x-axis. a) Each row in the y-axis represents the genotype of a RIAIL, color coded according to its genetic background. The x-axis shows the physical distance in Mb along each linkage group. b) Genome-wide parental genotype distribution, with the color indicating the proportion of parental genotype occurrence (on the y-axis) per physical distance (on the x-axis). Linkage groups including the X chromosome (L.7) were plotted separately to show the distribution of heterozygous, homozygous, and missing genotypes. Parental genotypes vary consistently per linkage group due to differences in recombination frequencies between linkage groups and within a linkage group. For details about the genotype composition and genotype vs the distance per phenotype of individual RIAILs, see Supplementary Figs. 3 and 4.

Characteristics of the A. freiburgense genome

The resulting chromosomal-scale assembly, which was the result of the integration of the genetic and physical maps, spanned 53.5 Mb (Table 1; Supplementary Table 3). The genome contained 89.2% of complete nematode BUSCO orthogroups, similar to the genome of A. rhodense (Supplementary Table 3). The annotation of repeat sequences (as described in the Materials and methods section) revealed that 10 Mb (18.2%) of the genome was repetitive (Fig. 6; Supplementary Table 4). Using a combination of ab initio predictors and evidence-based methods (including curated proteins and transcriptomic data), we predicted 8,759 protein-coding genes (Fig. 6 and Table 1). To assess the completeness and quality of the predicted protein-coding genes, we used BUSCO (nematoda_odb10 database) on the annotated transcript data set. The annotated transcripts contained 84.8% complete BUSCO orthogroups with 1.5% duplicated BUSCOs. The total percentage of repeat DNA in the table excludes the percentage of repeat DNA in unplaced scaffolds (1.5%).

Fig. 6. The density of genes, repeats, and GC content across the 7 chromosomal-scale scaffolds of A. freiburgense. Gene density, repeats, and GC content were plotted along the A. freiburgense chromosomal-scale assembly using a window size of 200 kb.

Next, we examined the chromosomal-scale genome of A. freiburgense to examine the evolution of chromosomal gene content in nematodes. Previous research has identified 7 possible ancestral chromosomal units of Rhabditida nematodes, called Nigon elements (Tandonnet et al. 2019; Gonzalez De La Rosa et al. 2021). They represent the reconstructed genic content of the hypothesized 7 chromosomes of the ancestor of all Rhabditida nematodes. The Nigon elements in A. freiburgense are mostly fragmented into different chromosomes (Fig. 7; Supplementary Table 5). Notably, Nigon B is fragmented into 3 different chromosomes in A. freiburgense (although it is intact in A. rhodense; Fig. 7a). In contrast, Nigon N, which is present in 2 chromosomes in A. rhodense, is mainly in 1 chromosome in A. freiburgense (Fig. 7b). When analyzed in context with other nematodes, A. freiburgense chromosomes underwent recent fission and fusion events of Nigon elements (Fig. 7b; Supplementary Fig. 5 and Table 5).

Fig. 7. Gene painting according to their Nigon element association. a) Nigon elements in A. freiburgense. b) Model of chromosome evolution in Rhabditida.

Identifying a candidate region on the X chromosome underlying the rates of production of males

NGS-BSA is a high-throughput strategy to identify potential QTL regions associated with traits of interest (Li and Xu 2022). We used NGS-BSA as a first approach using bulk DNA from a subset of the RIAILs (Supplementary Fig. 6 and Table 6), performed delta-SNP-index and G′ analysis, and calculated P-values and Q-values based on the tricube-smoothed G′ values (see Materials and methods; Magwene et al. 2011; Takagi et al. 2013). Two peaks of high delta-SNP-index and G′ values were associated with HM rates based on a FDR of 0.1 (Supplementary Figs. 7 and 8). The first region is on L.5, about 6.3 Mb long, and the second is on L.7 (chrX), approximately 1.0 Mb long (Supplementary Fig. 8 and Tables 7 and 8). It is important to acknowledge the small size of bulks from the total RIAILs in the BSA analysis, together with the increased risk of false positives associated with the FDR. This underscores the need for a cautious interpretation of the observed associations. Therefore, we performed a QTL analysis using a total of 261,860 genome-wide SNPs from all RIAILs.

A whole-genome QTL scan employing a single QTL model, using Haley–Knott regression and the rate of production of males from RIAIL crosses as a phenotype, revealed a prominent QTL peak on chromosome 7 (Haley and Knott 1992). We identified a candidate region spanning ∼420 kb (positions 1,949,767 and 2,375,391 bp) on the X chromosome (chr7) based on genome-wide P-value thresholds < 0.05 corresponding to a LOD score of 2.95 obtained through 1,000 times permutation of the first scan (Fig. 8a). The QTL peak is located at 2,135,500 bp with a LOD score of 3.17, and 83% of lines contain the homozygous APS14 genotype at the QTL peak and 17% contain the AP7 genotype (Supplementary Figs. 9a and 12). The mean percentage of males in lines with the homozygous APS14 genotype at the QTL peak is 33.1%, compared to just 9.3% in lines with the homozygous APS7 genotype (Supplementary Fig. 9b). Only 2 lines out of the total 57 HM lines (3.5%) have the APS7 genotype at the QTL peak, while the majority contain the APS14 genotype (Supplementary Fig. 12a). On the other hand, LM lines contain 13/40 (32.5%) APS7 genotypes at the QTL peak (Datasheet 2). To determine if there is an interaction between the identified QTL and the rest of the genome, we performed a scan using the QTL genotype as a covariate. However, no significant QTL peaks were detected (Supplementary Fig. 9c–e).

Fig. 8. Genome-wide QTL scan and genome-wide population differentiation between HM and LM RIAILs. a) A whole-genome QTL scan with a single QTL model, using Haley–Knott regression and the rate of male production as a phenotype, identified a significant QTL region spanning ∼420 kb between 1,949,767 and 2,375,391 bp on the X chromosome (chr7). LOD scores resulting from the QTL scan were plotted on the y-axis for each segregating marker across the genome on the x-axis. a) The genome-wide P-value threshold <0.05 corresponding to a LOD score of 2.95 was obtained through 1,000 times permutation of the first scan (horizontal dotted line). Population differentiation using an Fst analysis between HM lines and LM lines revealed a region of high differentiation aligning with the QTL region. b) Genetic differentiation between HM and LM lines was calculated using 5 kb weighted-windowed Fst with a step size of 1 kb and plotted on the y-axis against the genomic position in Mb (x-axis). c) Weir and Cockerham wcFst for each SNP was calculated and plotted on the y-axis against the genomic position in Mb (x-axis). d) P-value based on the LRT to quantify allele frequency differences between populations (pFst) is plotted on the y-axis against the genomic position in Mb (x-axis). d) Significant threshold at P-value 0.05 is plotted in dotted horizontal line. The smoothed line using a GAM on all population differentiation analyses between the HM line and LM line peaks exactly at the identified QTL region b–d).

Population differentiation using an Fst analysis between HM lines and LM lines revealed a region of high differentiation aligning with the QTL region. Genetic differentiation between HM and LM lines was calculated using 5 kb-windowed Fst with a step size of 1 kb, using Weir and Cockerham wcFst, and a P-value based on the likelihood ratio test (LRT) to quantify allele frequency differences between populations pFst (Fig. 8b–d). The smoothed line using a generalized additive model (GAM) on all population differentiation analyses between the HM line and LM line peaks exactly at the identified QTL region (Fig. 8b–d).

Identification of genes for polarity establishment in QTL region

The variant calling analysis (Van Der Auwera et al. 2013) detected 354,705 variants between APS7 and APS14, resulting in an average density of 1 variant per 155 bp. Most of these variants (84.1%) were SNPs (Supplementary Fig. 10). Of the identified SNPs, 0.1% were nonsense, and 5% were missense SNPs. In addition to SNPs, 8.6% of all variants are insertions, and 7.2% are deletions in the APS14 genome compared to the APS7 reference genome (Supplementary Fig. 11 and Table 7). The QTL region illustrates a mixing of haplotype blocks from both parental backgrounds (Supplementary Fig. 12). Within the QTL region, the APS14 genome, compared to the APS7 reference, contains 290 insertions, 269 small deletions, and 6 large structural deletions (Supplementary Table 7; Datasheets 4 and 5).

The QTL region contains 76 predicted protein-coding genes and 205 repeat motifs (Datasheets 6 and 7). SnpEff is a variant annotation and effect prediction tool to annotate and predict the effects of genetic variants, such as SNPs, insertions, deletions, and structural variants, on genes and proteins (Cingolani, Platts, et al. 2012). Using this tool, 181 missense variants with high LOD scores were detected within the QTL region on the X chromosome (Fig. 9; Datasheets 4–8). Forty-seven of the 76 (70.1%) predicted protein-coding genes have at least 1 missense variation, causing an amino acid change (Datasheet 6).

Fig. 9. QTL region significantly differentiated nonsynonymous SNPs. The LOD score of SNPs in the QTL region (y-axis) was plotted against the QTL bp (x-axis) position a). The −log(P-value) of QTL SNP differentiation was plotted (y-axis) against the QTL bp position in the x-axis. The horizontal dotted line in a) and b) indicates the significance threshold. Vertical lines in a) and b) indicate the genomic location of candidate genes involved in maintaining polarity and controlling cell division. c) The top half of the plot shows the APS14 variants (relative to the reference genome of APS7) with large structural variant deletions and small indels. The bottom half of the plot shows the APS7 annotation of repeats and genes.

To identify candidate genes within the QTL region involved in establishing the asymmetry in A. freiburgense spermatogenesis, we used Gene Ontology (GO) enrichment terms. In this region, we identified a few genes that are known to be involved in cell polarity in other systems: 2 copies of genes similar to serine/threonine-protein kinase sax-1 (Afr_13764 and Afr_13767; Zallen et al. 2000), a C2 domain-containing protein 2 (Afr_13773; Johnson et al. 2023), and a gene similar to the transcription factor pax-6 (Afr_13785; Asami et al. 2011). A full description of GO terms and GO IDs for genes located in the QTL region is provided in Datasheet 6. Nonsynonymous APS14 SNPs with a significant LOD score were identified in the C2 domain-containing protein 2, and only synonymous SNPs were identified in the pax-6 homolog (Datasheets 6 and 8). However, we did not detect any SNPs with significant LOD scores or differentiation in the gene regions of the duplicated sax-1-like genes.

Discussion

Auanema populations consist of selfing hermaphrodites, females, and males (Félix 2004; Kanzaki et al. 2017). Female vs hermaphrodite development is determined intergenerationally according to maternal age or exposure to social cues (Chaudhuri et al. 2011, 2015; Zuco et al. 2018; Tandonnet et al. 2019; Robles et al. 2021; Adams et al. 2022; Tandonnet et al. 2022), while males are determined chromosomally (Shakes et al. 2011).

In Auanema, the process of spermatogenesis in XO males and XX hermaphrodites involves the elimination of nullo-X spermatids (Tandonnet et al. 2018). Males produce predominantly X-bearing sperm (Shakes et al. 2011; Winter et al. 2017), resulting in non-Mendelian transmission of the X chromosome and thus giving rise to mostly XX offspring when mating with females (Félix 2004). When mating with hermaphrodites, males sire only sons (Tandonnet et al. 2018). This occurs due to meiotic nondisjunction of the X chromosome in hermaphrodite oocytes (Tandonnet et al. 2018), which leads to the production of nullo-X oocytes. As a result, the resulting embryos lack an X chromosome, leading to the development of male individuals when fertilized by the X-bearing sperm of the male (Tandonnet et al. 2018).

Unlike in C. elegans spermatogenesis, where nonsperm components are deposited into a residual body without DNA (Ward et al. 1981), Auanema eliminates those cytoplasmic components together with autosomes. Similar phenomena are observed in other organisms, such as sciarid flies, which undergo asymmetric segregation of chromosomes during meiosis and eliminate an entire set of chromosomes into a residual body (for review, see Gerbi 1986; Goday and Esteban 2001). Scale insects also remove sets of heterochromatic chromosomes during spermatogenesis (Bongiorni et al. 2004).

The removal of DNA observed in Auanema and other organisms represents a type of programmed DNA elimination (PDE), which has independently evolved in different taxa (Wang and Davis 2014; Dedukh and Krasikova 2022; Drotos et al. 2022; Kloc et al. 2022). PDE processes may influence sex ratios in various species. In Nasonia wasps, for instance, PDE removes sets of chromosomes to convert female embryos into males (Nur et al. 1988; Werren and Stouthamer 2003). The nullo-X sperm of the nematode Strongyloides spp., responsible for male development, are eliminated during spermatogenesis (Dulovic et al. 2022). To generate males, specific portions of the X chromosome are targeted for degradation in XX embryos (Nemetschke et al. 2010). Hermaphrodites of the nematode Rhabdias also expel one of their X chromosomes into a residual body during spermatogenesis (Runey et al. 1978). The mechanisms by which specific chromosomes are targeted for PDE are still not well understood.

The study of meiotic drive and PDE mechanisms has often been hindered by the lack of genetic tools available for many organisms that exhibit these processes. However, Auanema has emerged as a promising model system for investigating these mechanisms due to its advantageous features, including easy husbandry, short life cycle, large brood sizes, and the recent development of genetic tools and resources (Adams et al. 2019; Tandonnet et al. 2019; Kranse et al. 2021; Tandonnet et al. 2022).

In this study, we aimed to improve our understanding of the mechanisms of X chromosome drive in A. freiburgense, generating genetic resources in the form of a chromosomal-scale assembly of its genome. To enhance the contiguity of the assembly, we constructed independent genetic linkage maps for the autosomes and the X chromosome. By accounting for differences in recombination events between these chromosomes during the inbreeding phase of RIAILs, we were able to improve the contiguity of the X chromosome assembly. Our efforts resulted in a 47.3% increase in the size of the X chromosome compared to the initial genetic linkage map. This improved assembly will provide valuable insights into the unique biology and mechanisms underlying PDE in A. freiburgense. We observed an inverse relationship between the genetic distance in cM and the chromosome size in bp between the independently constructed genetic linkage maps for the X chromosome (LG7) and the autosomes especially (LG4, LG5, and LG7). This disparity potentially arises from the difference in segregation distortion handling methodologies in the construction of the 2 independent genetic linkage maps and the inclusion of heterozygous markers exclusively for the X chromosome genetic linkage map. The level of heterozygosity in the X chromosome was less than previously observed in RILs of A. rhodense (Tandonnet et al. 2018). This could indicate that A. freiburgense hermaphrodites do not undergo the noncanonical meiosis typical of A. rhodense (Fig. 1) or may reflect differences in population dynamics between the 2 species during the inbreeding process. Auanema rhodense hermaphrodites are produced concurrently with females and male siblings (Chaudhuri et al. 2011). In contrast, A. freiburgense populations must become crowded before hermaphrodites are produced (Félix 2004; Robles et al. 2021). Thus, inbreeding populations of A. freiburgensis will have likely already gone through female sister/male brother mating (with canonical X meiosis) before hermaphrodites are produced in each generation, potentially reducing heterozygosity in the RIAILs.

We serendipitously discovered that by crossing 2 inbred strains of A. freiburgense, males of hybrid lines can produce more sons than the parental strains. By using genetic markers specific to the X chromosome, we found that males in these transgressive lines produce viable nullo-X sperm (Fig. 4; Supplementary Table 1), which is absent in the parental lines. We hypothesize that HM lines display alterations of the asymmetric spermatogenesis to either a C. elegans-like system with the formation of a central residual body or a completely symmetric spermatogenesis, without a central residual body formation (Fig. 10). As a result, when these males mate with females, there is an increased proportion of sons.

Fig. 10. Proposed models of male spermatogenesis in HM- and LM-RIAILs. a) Male cross-progeny resulting from a cross between APS7 females and LM lines inherited the paternal X chromosome, indicating that nullo-X sperm are discarded during the spermatogenesis of males from those lines. b) Nullo-X sperm could be produced by the symmetrical distribution of sperm components between X-bearing and nullo-X sperm during anaphase II. The symmetric segregation could be total, where subcellular compartments, including sperm components, segregate equally, or a C. elegans-like anaphase II, where a central residual body is formed for discarded materials.

Since the same allele combination was present in the parental line APS14, we can discard a model based on dominance or exposure to recessive alleles (Rieseberg et al. 1999). Instead, the “HM” phenotype is likely to be the result of epistatic interactions between a combination of homozygous alleles of more than 1 locus. Although the BSA indicated that a second locus associated with high production of males is located in chromosome 5, we could not confirm this by QTL mapping using whole-genome sequencing of all RIAILs. The success and precision of NGS-BSA rely on the high number of individuals used per bulk and the coverage of sequencing. Even though we sequenced bulks at high coverage, there were only 10 lines per bulk. We reasoned that the genetic variation in the small number of lines selected per bulk did not represent the genetic variation in all of the RIAILs sharing a similar phenotype. As a result, another candidate region that is probably not associated with the phenotype was detected. To improve confidence in NGS-BSA, future experiments should increase the number of lines per bulk by pooling all the HM-RIAILs in 1 bulk and the LM-RIAILs in another bulk to capture all the genetic variations within the RIAILs.

In this study, we established that DNA elimination (i.e. elimination of autosomes from nullo-X spermatids) during sex determination is driven by epistatic interactions between specific allele combinations. This process leads to the polarization of the spermatocyte cytoplasm, resulting in the production of viable sperm and nonviable spermatids. The presence of the APS14 genome within the QTL region is a key factor in the transgressive phenotype observed, with 83% of lines exhibiting a homozygous APS14 genotype at the QTL peak marker, compared to 17% with the APS7 genotype at the QTL peak marker, as illustrated in Supplementary Fig. 9a. Among the 57 HM lines studied, only 2 lines (representing 3.5%) possess the APS7 genotype at the peak marker, with the majority displaying the APS14 genotype. On the other hand, LM lines contain 13/40 (32.5%) APS7 genotypes at the QTL peak marker (Datasheet 2). The mechanism of nullo-X elimination by a polarizing cue from the X chromosome is still unclear. However, in the ∼420 kb mapped QTL region, we identified genes that are potentially interesting for further investigation, including genes coding for serine/threonine-protein kinase sax-1, a C2 domain-containing protein 2 type, and pax-6. Mutants of sax-1 in C. elegans exhibit defects in neuronal cell shape and polarity (Zallen et al. 2000). At this stage, we cannot confirm with confidence if the sax-1 gene is also duplicated in APS14 or is present as a single copy. C2 domain-containing protein 2 (Afr_13773; Johnson et al. 2023) and the transcription factor pax-6 (Afr_13785; Asami et al. 2011) are involved in the polarization of neuronal cells in mammalian systems.

It is possible that additional structural variants, not covered in this study, mediate the polarization of cytoplasmic components during meiosis. Therefore, an investigation of the QTL region by deep sequencing with long reads of the APS14 strain and RIAILs could shed more light on the differences in coding sequences, regulating elements, and structural differences between the strains.

Supplementary Material

iyae032_Supplementary_Data

iyae032_Peer_Review_History

Data availability

We deposited all sequencing data at the European Nucleotide Archive (ENA) under BioProject numbers PRJEB55706 (Illumina genomic data), PRJEB60474 and PRJEB50372 (Illumina transcriptomic data), and PRJNA640723 (PacBio data; Supplementary Table 10). Furthermore, we uploaded the genome assembly, annotation, and predicted protein data set to NCBI with accession PRJNA947217 and to figshare (https://doi.org/10.6084/m9.figshare.24885039.v1). Linkage groups L.1–L.7 are named chr1 to chr7 in the submitted genome. The RIAILs are available upon request from the authors. Datasheets 1–8 are published at https://doi.org/10.6084/m9.figshare.24901743.

Supplemental material available at GENETICS online.

Funding

AP-dS and SA were supported by grants from BBSRC (BB/L019884/1) and Leverhulme Trust (RPG-2019-329). ST was funded by a full PhD scholarship from the program Ciência sem Fronteiras (Conselho Nacional de Desenvolvimento Científico e Tecnológico agency, process number 201116/2014-6). AT was funded by the Doctoral Training Program from Natural Environment Research Council (NERC CENTA). JL was supported by Samsung Science and Technology Foundation SSTF-BA1501-52.
==== Refs
Literature cited

Abyzov  A, Urban  AE, Snyder  M, Gerstein  M. 2011. CNVnator: an approach to discover, genotype, and characterize typical and atypical CNVs from family and population genome sequencing. Genome Res. 21 (6 ):974–984. doi:10.1101/gr.114876.110.21324876
Adams  S, Pathak  P, Kittelmann  M, Jones  ARC, Mallon  EB, Pires-daSilva  A. 2022. Sexual morph specialisation in a trioecious nematode balances opposing selective forces. Sci Rep. 12 (1 ):6402. doi:10.1038/s41598-022-09900-8.35431314
Adams  S, Pathak  P, Shao  H, Lok  JB, Pires-daSilva  A. 2019. Liposome-based transfection enhances RNAi and CRISPR-mediated mutagenesis in non-model nematode systems. Sci Rep. 9 (1 ):483. doi:10.1038/s41598-018-37036-1.30679624
Al-Yazeedi  T, Xu  EL, Kaur  J, Shakes  DC, Pires-daSilva  A. 2022. Lagging X chromatids specify the orientation of asymmetric organelle partitioning in XX spermatocytes of Auanema rhodensis. Genetics  222 (4 ):iyac159. doi:10.1093/genetics/iyac159.36255260
Asami  M, Pilz  GA, Ninkovic  J, Godinho  L, Schroeder  T, Huttner  WB, Götz  M. 2011. The role of Pax6 in regulating the orientation and mode of cell division of progenitors in the mouse cerebral cortex. Development  138 (23 ):5067–5078. doi:10.1242/dev.074591.22031545
Benjamini  Y, Hochberg  Y. 1995. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Series B Stat Methodol. 57 :289–300. doi:10.1111/j.2517-6161.1995.tb02031.x.
Blackman  RL . 1985. Spermatogenesis in the aphid Amphorophora tuberculata (Homoptera, Aphididae). Chromosoma  92 (5 ):357–362. doi:10.1007/BF00327467.
Bongiorni  S, Fiorenzo  P, Pippoletti  D, Prantera  G. 2004. Inverted meiosis and meiotic drive in mealybugs. Chromosoma  112 (7 ):331–341. doi:10.1007/s00412-004-0278-4.15095094
Broman  KW, Wu  H, Sen  S, Churchill  GA. 2003. R/qtl: QTL mapping in experimental crosses. Bioinformatics  19 (7 ):889–890. doi:10.1093/bioinformatics/btg112.12724300
Brouard  JS, Schenkel  F, Marete  A, Bissonnette  N. 2019. The GATK joint genotyping workflow is appropriate for calling variants in RNA-seq experiments. J Anim Sci Biotechnol. 10 (1 ):44. doi:10.1186/s40104-019-0359-0.31249686
Burt  A, Trivers  R. 2006. Genes in Conflict: the Biology of Selfish Genetic Elements. Cambridge (MA): Belknap Press of Harvard University Press.
Camacho  C, Coulouris  G, Avagyan  V, Ma  N, Papadopoulos  J, Bealer  K, Madden  TL. 2009. BLAST+: architecture and applications. BMC Bioinformatics  10 (1 ):421. doi:10.1186/1471-2105-10-421.20003500
Chaudhuri  J, Bose  N, Tandonnet  S, Adams  S, Zuco  G, Kache  V, Parihar  M, von Reuss  SH, Schroeder  FC, Pires-daSilva  A. 2015. Mating dynamics in a nematode with three sexes and its evolutionary implications. Sci Rep. 5 (1 ):17676. doi:10.1038/srep17676.26631423
Chaudhuri  J, Kache  V, Pires-daSilva  A. 2011. Regulation of sexual plasticity in a nematode that produces males, females, and hermaphrodites. Curr Biol. 21 (18 ):1548–1551. doi:10.1016/j.cub.2011.08.009.21906947
Chen  K, Wallis  JW, McLellan  MD, Larson  DE, Kalicki  JM, Pohl  CS, McGrath  SD, Wendl  MC, Zhang  Q, Locke  DP, et al  2009. BreakDancer: an algorithm for high-resolution mapping of genomic structural variation. Nat Methods. 6 (9 ):677–681. doi:10.1038/nmeth.1363.19668202
Chen  X, Schulz-Trieglaff  O, Shaw  R, Barnes  B, Schlesinger  F, Källberg  M, Cox  AJ, Kruglyak  S, Saunders  CT. 2016. Manta: rapid detection of structural variants and indels for germline and cancer sequencing applications. Bioinformatics  32 (8 ):1220–1222. doi:10.1093/bioinformatics/btv710.26647377
Churchill  GA, Doerge  RW. 1994. Empirical threshold values for quantitative trait mapping. Genetics  138 (3 ):963–971. doi:10.1093/genetics/138.3.963.7851788
Cingolani  P, Patel  VM, Coon  M, Nguyen  T, Land  SJ, Ruden  DM, Lu  X. 2012. Using Drosophila melanogaster as a model for genotoxic chemical mutational studies with a new program, SnpSift. Front Genet. 3 :35. doi:10.3389/fgene.2012.00035.22435069
Cingolani  P, Platts  A, Wang  LL, Coon  M, Nguyen  T, Wang  L, Land  SJ, Lu  X, Ruden  DM. 2012. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly (Austin). 6 (2 ):80–92. doi:10.4161/fly.19695.22728672
Danecek  P, Auton  A, Abecasis  G, Albers  CA, Banks  E, DePristo  MA, Handsaker  RE, Lunter  G, Marth  GT, Sherry  ST, et al  2011. The variant call format and VCFtools. Bioinformatics  27 (15 ):2156–2158. doi:10.1093/bioinformatics/btr330.21653522
Dedukh  D, Krasikova  A. 2022. Delete and survive: strategies of programmed genetic material elimination in eukaryotes. Biol Rev Camb Philos Soc. 97 (1 ):195–216. doi:10.1111/brv.12796.34542224
DePristo  MA, Banks  E, Poplin  R, Garimella  KV, Maguire  JR, Hartl  C, Philippakis  AA, del Angel  G, Rivas  MA, Hanna  M, et al  2011. A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nat Genet. 43 (5 ):491–498. doi:10.1038/ng.806.21478889
Donath  A, Juhling  F, Al-Arab  M, Bernhart  SH, Reinhardt  F, Stadler  PF, Middendorf  M, Bernt  M. 2019. Improved annotation of protein-coding genes boundaries in metazoan mitochondrial genomes. Nucleic Acids Res. 47 (20 ):10543–10552. doi:10.1093/nar/gkz833.31584075
Drotos  KHI, Zagoskin  MV, Kess  T, Gregory  TR, Wyngaard  GA. 2022. Throwing away DNA: programmed downsizing in somatic nuclei. Trends Genet. 38 (5 ):483–500. doi:10.1016/j.tig.2022.02.003.35227512
Dulovic  A, Koch  I, Hipp  K, Streit  A. 2022. Strongyloides spp. eliminate male-determining sperm post-meiotically. Mol Biochem Parasitol. 251 :111509. doi:10.1016/j.molbiopara.2022.111509.35985494
Dupuis  J, Brown  PO, Siegmund  D. 1995. Statistical methods for linkage analysis of complex traits from high-resolution maps of identity by descent. Genetics  140 (2 ):843–856. doi:10.1093/genetics/140.2.843.7498758
Edgar  RC . 2010. Search and clustering orders of magnitude faster than BLAST. Bioinformatics  26 (19 ):2460–2461. doi:10.1093/bioinformatics/btq461.20709691
Ellinghaus  D, Kurtz  S, Willhoeft  U. 2008. LTRharvest, an efficient and flexible software for de novo detection of LTR retrotransposons. BMC Bioinformatics  9 (1 ):18. doi:10.1186/1471-2105-9-18.18194517
Félix  MA . 2004. Alternative morphs and plasticity of vulval development in a rhabditid nematode species. Dev Genes Evol. 214 (2 ):55–63. doi:10.1007/s00427-003-0376-y.14730447
Gerbi  SA . 1986. Unusual chromosome movements in sciarid flies. In: Hennig  W, editor. Germ Line—Soma Differentiation. Heidelberg (Berlin): Springer. p. 71–104.
Goday  C, Esteban  MR. 2001. Chromosome elimination in sciarid flies. Bioessays  23 (3 ):242–250. doi:10.1002/1521-1878(200103)23:3<242::AID-BIES1034>3.0.CO;2-P.11223881
Gonzalez de la Rosa  PM, Thomson  M, Trivedi  U, Tracey  A, Tandonnet  S, Blaxter  M. 2021. A telomere-to-telomere assembly of Oscheius tipulae and the evolution of rhabditid nematode chromosomes. G3 (Bethesda). 11 (1 ):jkaa20. doi:10.1093/g3journal/jkaa020.
Götz  S, Garcia-Gomez  JM, Terol  J, Williams  TD, Nagaraj  SH, Nueda  MJ, Robles  M, Talon  M, Dopazo  J, Conesa  A. 2008. High-throughput functional annotation and data mining with the Blast2GO suite. Nucleic Acids Res. 36 (10 ):3420–3435. doi:10.1093/nar/gkn176.18445632
Grabherr  MG, Haas  BJ, Yassour  M, Levin  JZ, Thompson  DA, Amit  I, Adiconis  X, Fan  L, Raychowdhury  R, Zeng  Q, et al  2011. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 29 (7 ):644–652. doi:10.1038/nbt.1883.21572440
Gremme  G, Steinbiss  S, Kurtz  S. 2013. GenomeTools: a comprehensive software library for efficient processing of structured genome annotations. IEEE/ACM Trans Comput Biol Bioinform. 10 (3 ):645–656. doi:10.1109/TCBB.2013.68.24091398
Haas  B . 2010. TransposonPSI: An Application of PSI-Blast to Mine. (Retro-)Transposon ORF Homologies.  https://transposonpsi.sourceforge.net/.
Haley  CS, Knott  SA. 1992. A simple regression method for mapping quantitative trait loci in line crosses using flanking markers. Heredity (Edinb). 69 (4 ):315–324. doi:10.1038/hdy.1992.131.16718932
Hoeh  WR, Blakley  KH, Brown  WM. 1991. Heteroplasmy suggests limited biparental inheritance of Mytilus mitochondrial DNA. Science  251 (5000 ):1488–1490. doi:10.1126/science.1672472.1672472
Holt  C, Yandell  M. 2011. MAKER2: an annotation pipeline and genome-database management tool for second-generation genome projects. BMC Bioinformatics  12 (1 ):491. doi:10.1186/1471-2105-12-491.22192575
Hurst  LD . 1993. The incidences. Mechanisms and evolution of cytoplasmic sex ratio distorters in animals. Biol Rev. 68 (1 ):121–194. doi:10.1111/j.1469-185X.1993.tb00733.x.
Jiang  H, Lei  R, Ding  SW, Zhu  S. 2014. Skewer: a fast and accurate adapter trimmer for next-generation sequencing paired-end reads. BMC Bioinformatics  15 (1 ):182. doi:10.1186/1471-2105-15-182.24925680
Johnson  B, Iuliano  M, Lam  T, Biederer  T, De Camilli  P. 2023. A complex of the lipid transport ER proteins TMEM24 and C2CD2 with band 4.1 at cell-cell contacts. bioRxiv. doi:10.1101/2023.12.06.570396.
Kanzaki  N, Kiontke  K, Tanaka  R, Hirooka  Y, Schwarz  A, Müller-Reichert  T, Chaudhuri  J, Pires-daSilva  A. 2017. Description of two three-gendered nematode species in the new genus Auanema (Rhabditina) that are models for reproductive mode evolution. Sci Rep. 7 (1 ):11135. doi:10.1038/s41598-017-09871-1.28894108
Kloc  M, Kubiak  JZ, Ghobrial  RM. 2022. Natural genetic engineering: a programmed chromosome/DNA elimination. Dev Biol. 486 :15–25. doi:10.1016/j.ydbio.2022.03.008.35321809
Koren  S, Walenz  BP, Berlin  K, Miller  JR, Bergman  NH, Phillippy  AM. 2017. Canu:scalable and accurate long-read assembly via adaptive k-mer weighting and repeat separation. Genome Res. 27 (5 ):722–736. doi:10.1101/gr.215087.116.28298431
Korf  I . 2004. Gene finding in novel genomes. BMC Bioinformatics  5 (1 ):59. doi:10.1186/1471-2105-5-59.15144565
Kranse  O, Beasley  H, Adams  S, Pires-daSilva  A, Bell  C, Lilley  CJ, Urwin  PE, Bird  D, Miska  E, Smant  G, et al  2021. Toward genetic modification of plant-parasitic nematodes: delivery of macromolecules to adults and expression of exogenous mRNA in second stage juveniles. G3 (Bethesda). 11 :jkaa058. doi:10.1093/g3journal/jkaa058.33585878
Krzywinski  M, Schein  J, Birol  I, Connors  J, Gascoyne  R, Horsman  D, Jones  SJ, Marra  MA. 2009. Circos: an information aesthetic for comparative genomics. Genome Res. 19 :1639–1645. doi:10.1101/gr.092759.109.19541911
Kurtz  S, Phillippy  A, Delcher  AL, Smoot  M, Shumway  M, Antonescu  C, Salzberg  SL. 2004. Versatile and open software for comparing large genomes. Genome Biol. 5 :R12. doi:10.1186/gb-2004-5-2-r12.14759262
Layer  RM, Chiang  C, Quinlan  AR, Hall  IM. 2014. LUMPY: a probabilistic framework for structural variant discovery. Genome Biol. 15 :R84. doi:10.1186/gb-2014-15-6-r84.24970577
Li  H, Durbin  R. 2009. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics  25 :1754–1760. doi:10.1093/bioinformatics/btp324.19451168
Li  H, Handsaker  B, Wysoker  A, Fennell  T, Ruan  J, Homer  N, Marth  G, Abecasis  G, Durbin  R. 2009. The Sequence Alignment/Map format and SAMtools. Bioinformatics  25 :2078–2079. doi:10.1093/bioinformatics/btp352.19505943
Li  Z, Xu  Y. 2022. Bulk segregation analysis in the NGS era: a review of its teenage years. Plant J. 109 :1355–1374. doi:10.1111/tpj.15646.34931728
Lomsadze  A, Ter-Hovhannisyan  V, Chernoff  YO, Borodovsky  M. 2005. Gene identification in novel eukaryotic genomes by self-training algorithm. Nucleic Acids Res. 33 :6494–6506. doi:10.1093/nar/gki937.16314312
Lyttle  TW . 1991. Segregation distorters. Ann Rev Ecol Evol Syst. 25 :511–557. doi:10.1146/annurev.ge.25.120191.002455.
Magwene  PM, Willis  JH, Kelly  JK. 2011. The statistics of bulk segregant analysis using next generation sequencing. PLoS Comput Biol. 7 :e1002255. doi:10.1371/journal.pcbi.1002255.22072954
Mansfeld  BN, Grumet  R. 2018. QTLseqr: an R package for bulk segregant analysis with next-generation sequencing. Plant Genome. 11 . doi:10.3835/plantgenome2018.01.0006.
McKenna  A, Hanna  M, Banks  E, Sivachenko  A, Cibulskis  K, Kernytsky  A, Garimella  K, Altshuler  D, Gabriel  S, Daly  M, et al  2010. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 20 :1297–1303. doi:10.1101/gr.107524.110.20644199
Nemetschke  L, Eberhardt  AG, Hertzberg  H, Streit  A. 2010. Genetics, chromatin diminution, and sex chromosome evolution in the parasitic nematode genus Strongyloides. Curr Biol. 20 :1687–1696. doi:10.1016/j.cub.2010.08.014.20832309
Nur  U, Werren  JH, Eickbush  DG, Burke  WD, Eickbush  TH. 1988. A “selfish” B chromosome that enhances its transmission by eliminating the paternal genome. Science  240 :512–514. doi:10.1126/science.3358129.3358129
Patel  RK, Jain  M. 2012. NGS QC Toolkit: a toolkit for quality control of next generation sequencing data. PLoS One  7 :e30619. doi:10.1371/journal.pone.0030619.22312429
Pires-daSilva  A . 2007. Evolution of the control of sexual identity in nematodes. Semin Cell Dev Biol. 18 :362–370. doi:10.1016/j.semcdb.2006.11.014.17306573
Pires-daSilva  A . 2013. Pristionchus pacificus Protocols. Bethesda, USA: WormBook. p. 1–20.
Poplin  R, Ruano-Rubio  V, DePristo  MA, Fennell  TJ, Carneiro  MO, Van der Auwera  GA, Kling  DE, Gauthier  LD, Levy-Moonshine  A, Roazen  D, et al  2018. Scaling accurate genetic variant discovery to tens of thousands of samples. bioRxiv. 201178. doi:10.1101/201178.
Pracana  R, Priyam  A, Levantis  I, Nichols  RA, Wurm  Y. 2017. The fire ant social chromosome supergene variant Sb shows low diversity but high divergence from SB. Mol Ecol. 26 :2864–2879. doi:10.1111/mec.14054.28220980
Rausch  T, Zichner  T, Schlattl  A, Stutz  AM, Benes  V, Korbel  JO. 2012. DELLY: structural variant discovery by integrated paired-end and split-read analysis. Bioinformatics  28 :i333–i339. doi:10.1093/bioinformatics/bts378.22962449
Rieseberg  LH, Archer  MA, Wayne  RK. 1999. Transgressive segregation, adaptation and speciation. Heredity (Edinb). 83 (Pt 4 ):363–372. doi:10.1038/sj.hdy.6886170.10583537
Robles  P, Turner  A, Zuco  G, Adams  S, Paganopolou  P, Winton  M, Hill  B, Kache  V, Bateson  C, Pires-daSilva  A. 2021. Parental energy-sensing pathways control intergenerational offspring sex determination in the nematode Auanema freiburgensis. BMC Biol. 19 :102. doi:10.1186/s12915-021-01032-1.34001117
Robles  P, Turner  A, Zuco  G, Paganopolou  P, Hill  B, Kache  V, Bateson  C, Pires-daSilva  A. 2020. Chromatin-associated effectors of energy-sensing pathways mediate intergenerational effects. bioRxiv. 2020.2008.2031.275727. 10.1101/2020.08.31.275727.
Ross  L, Mongue  AJ, Hodson  CN, Schwander  T. 2022. Asymmetric inheritance: the diversity and evolution of non-Mendelian reproductive strategies. Annu Rev Ecol Evol Syst. 53 :1–23. doi:10.1146/annurev-ecolsys-021822-010659.
Runey  WM, Runey  GL, Lauter  FH. 1978. Gametogenesis and fertilization in Rhabdias ranae Walton 1929: I. The parasitic hermaphrodite. J Parasitol. 64 :1008–1014. doi:10.2307/3279712.570219
Seppey  M, Manni  M, Zdobnov  EM. 2019. BUSCO: assessing genome assembly and annotation completeness. Methods Mol Biol. 1962 :227–245. doi:10.1007/978-1-4939-9173-0_14.31020564
Shakes  DC, Neva  BJ, Huynh  H, Chaudhuri  J, Pires-daSilva  A. 2011. Asymmetric spermatocyte division as a mechanism for controlling sex ratios. Nat Commun. 2 :157. doi:10.1038/ncomms1160.21245838
Shen  W, Le  S, Li  Y, Hu  F. 2016. SeqKit: a cross-platform and ultrafast toolkitfor FASTA/Q file manipulation. PLoS One  11 :e0163962. doi:10.1371/journal.pone.0163962.27706213
Smit  AFA, Hubley  R. 2008-2015. RepeatModeler Open-1.0.  http://www.repeatmasker.org.
Smit  AFA, Hubley  R. 2013-2015. RepeatMasker-4.0.  www.repeatmasker.org.
Smith  TF, Waterman  MS. 1981. Identification of common molecular subsequences. J Mol Biol. 147 :195–197. doi:10.1016/0022-2836(81)90087-5.7265238
Stanke  M, Waack  S. 2003. Gene prediction with a hidden Markov model and a new intron submodel. Bioinformatics  19 (Suppl 2) ):ii215–ii225. doi:10.1093/bioinformatics/btg1080.14534192
Steinbiss  S, Willhoeft  U, Gremme  G, Kurtz  S. 2009. Fine-grained annotation and classification of de novo predicted LTR retrotransposons. Nucleic Acids Res. 37 :7002–7013. doi:10.1093/nar/gkp759.19786494
Sudhaus  W . 2023. An update of the catalogue of paraphyletic ‘Rhabditidae’ (Nematoda) after eleven years. Soil Org. 95 :95–116. doi:10.25674/so95iss1id312.
Takagi  H, Abe  A, Yoshida  K, Kosugi  S, Natsume  S, Mitsuoka  C, Uemura  A, Utsushi  H, Tamiru  M, Takuno  S, et al  2013. QTL-seq: rapid mapping of quantitative trait loci in rice by whole genome resequencing of DNA from two bulked populations. Plant J. 74 :174–183. doi:10.1111/tpj.12105.23289725
Tandonnet  S, Farrell  MC, Koutsovoulos  GD, Blaxter  ML, Parihar  M, Sadler  PL, Shakes  DC, Pires-daSilva  A. 2018. Sex- and gamete-specific patterns of X chromosome segregation in a trioecious nematode. Curr Biol. 28 :93–99.e93. doi:10.1016/j.cub.2017.11.037.29276124
Tandonnet  S, Haq  M, Turner  A, Grana  T, Paganopoulou  P, Adams  S, Dhawan  S, Kanzaki  N, Nuez  I, Félix  M-A, et al  2022. De novo genome assembly of Auanema melissensis, a trioecious free-living nematode. J Nematol. 54 :20220059. doi:10.2478/jofnem-2022-0059.36879950
Tandonnet  S, Koutsovoulos  GD, Adams  S, Cloarec  D, Parihar  M, Blaxter  ML, Pires-daSilva  A. 2019. Chromosome-wide evolution and sex determination in the three-sexed nematode Auanema rhodensis. G3 (Bethesda). 9 :1211–1230. doi:10.1534/g3.119.0011.30770412
Tang  H, Zhang  X, Miao  C, Zhang  J, Ming  R, Schnable  JC, Schnable  PS, Lyons  E, Lu  J. 2015. ALLMAPS: robust scaffold ordering based on multiple maps. Genome Biol. 16 :3. doi:10.1186/s13059-014-0573-1.25583564
Taylor  J, Butler  D. 2017. R package ASMap: efficient genetic linkage map construction and diagnosis. J Stat Softw.  79 :29. doi:10.18637/jss.v079.i06.
Van der Auwera  GA, Carneiro  MO, Hartl  C, Poplin  R, del Angel  G, Levy-Moonshine  A, Jordan  T, Shakir  K, Roazen  D, Thibault  J, et al  2013. From FastQ data to high confidence variant calls: the Genome Analysis Toolkit best practices pipeline. Curr Protoc Bioinformatics. 43 :11 10 11–11 10.33. doi:10.1002/0471250953.bi1110s43.
Van Goor  J, Shakes  DC, Haag  ES. 2021. Fisher vs. the worms: extraordinary sex ratios in nematodes and the mechanisms that produce them. Cells  10 :1793. doi:10.3390/cells10071793.34359962
Walker  BJ, Abeel  T, Shea  T, Priest  M, Abouelliel  A, Sakthikumar  S, Cuomo  CA, Zeng  Q, Wortman  J, Young  SK, et al  2014. Pilon: an integrated tool for comprehensive microbial variant detection and genome assembly improvement. PLoS One  9 :e112963. doi:10.1371/journal.pone.0112963.25409509
Wang  J, Davis  RE. 2014. Programmed DNA elimination in multicellular organisms. Curr Opin Genet Dev. 27 :26–34. doi:10.1016/j.gde.2014.03.012.24886889
Ward  S, Argon  Y, Nelson  GA. 1981. Sperm morphogenesis in wild-type and fertilization-defective mutants of Caenorhabditis elegans. J Cell Biol. 91 :26–44. doi:10.1083/jcb.91.1.26.7298721
Werren  JH, Stouthamer  R. 2003. PSR (paternal sex ratio) chromosomes: the ultimate selfish genetic elements. Genetica  117 :85–101. doi:10.1023/A:1022368700752.12656576
Wilson  ACC, Sunnucks  P, Hales  DF. 1997. Random loss of X chromosome at male determination in an aphid, Sitobion near fragariae, detected using an X-linked polymorphic microsatellite marker. Genet Res. 69 :233–236. doi:10.1017/S0016672397002747.
Winter  ES, Schwarz  A, Fabig  G, Feldman  JL, Pires-daSilva  A, Müller-Reichert  T, Sadler  PL, Shakes  DC. 2017. Cytoskeletal variations in an asymmetric cell division support diversity in nematode sperm size and sex ratios. Development  144 :3253–3263. doi:10.1242/dev.153841.28827395
Xu  Z, Wang  H. 2007. LTR_FINDER: an efficient tool for the prediction of full-length LTR retrotransposons. Nucleic Acids Res. 35 :W265–W268. doi:10.1093/nar/gkm286.17485477
Zallen  JA, Peckol  EL, Tobin  DM, Bargmann  CI. 2000. Neuronal cell shape and neurite initiation are regulated by the Ndr kinase SAX-1, a member of the Orb6/COT-1/warts serine/threonine kinase family. Mol Biol Cell. 11 :3177–3190. doi:10.1091/mbc.11.9.3177.10982409
Zarate  S, Carroll  A, Mahmoud  M, Krasheninina  O, Jun  G, Salerno  WJ, Schatz  MC, Boerwinkle  E, Gibbs  RA, Sedlazeck  FJ. 2020. Parliament2: accurate structural variant calling at scale. Gigascience  9 :giaa145. doi:10.1093/gigascience/giaa145.33347570
Zouros  E, Freeman  KR, Ball  AO, Pogson  GH. 1992. Direct evidence for extensive paternal mitochondrial DNA inheritance in the marine mussel Mytilus. Nature  359 :412–414. doi:10.1038/359412a0.1357555
Zouros  E, Oberhauser Ball  A, Saavedra  C, Freeman  KR. 1994. An unusual type of mitochondrial DNA inheritance in the blue mussel Mytilus. Proc Natl Acad Sci U S A. 91 :7463–7467. doi:10.1073/pnas.91.16.7463.8052604
Zuco  G, Kache  V, Robles  P, Chaudhuri  J, Hill  B, Bateson  C, Pires-daSilva  A. 2018. Sensory neurons control heritable adaptation to stress through germline reprogramming. bioRxiv. 406033. doi:10.1101/406033.
Zuo  JF, Niu  Y, Cheng  P, Feng  JY, Han  SF, Zhang  Y-H, Shu  G, Wang  Y, Zhang  Y-M. 2019. Effect of marker segregation distortion on high density linkage map construction and QTL mapping in soybean (Glycine max L.). Heredity (Edinb). 123 :579–592. doi:10.1038/s41437-019-0238-7.31152165
