De novo Genome Assembly of Auanema melissensis, a Trioecious Free-Living Nematode

Auanema melissensis (aka, Auanema sp. JU1783) was maintained under the standard culture conditions of Caenorhabditis elegans (Chaudhuri et al., 2011b) at 20°C. Plates were seeded with the streptomycin-resistant Escherichia coli strain OP50-1. Microbial contamination was prevented by adding 50 mg/ml of streptomycin and 10 mg/ml of nystatin to the nematode growth medium (NGM).

To determine if A. melissensis was trioecious, we randomly isolated eggs from the culture plates. The eggs were individually placed into 48-well plates seeded with Escherichia coli OP50-1 and allowed to develop to adulthood in isolation. The sexual morph was determined by morphological characters and the ability or not to reproduce on their own. Hermaphrodites laid eggs in the absence of a mating partner, whereas females could only reproduce if paired with a male. Males were distinguished by their blunt tails.

Dauers of other Auanema species invariably develop into self-reproducing hermaphrodites (Félix, 2004; Chaudhuri et al., 2011a; Zuco et al., 2018; Robles et al., 2021). To determine if this was also the case in A. melissensis, we isolated dauers onto individual plates seeded with E. coli OP50-1 and determined their sexual morph in adulthood, as described above.

To test if the crowding conditions have an effect on the proportion of each sexual morph produced, we counted the number of male, female and hermaphrodite progeny produced by hermaphrodite mothers placed in the presence and absence of crowding cues. The crowding cues were prepared by first washing a Petri dish (ϕ = 60 mm) containing a crowded culture of nematodes with 1 ml of M9 buffer. The resuspended nematode culture was placed in a 50 ml tube and left on a thermomixer for 16–24 hr at 20°C. The liquid culture was then transferred to a 1.5 ml microcentrifuge tube and centrifuged for 45 min, at 15,000 rpm. Once clear from nematodes, the supernatant with the crowding cues was placed in a clean 1.5 ml tube and used immediately.

For the treatment conditions, we used Petri dishes (ϕ = 60 mm) with NGM that were seeded with 50 ml of OP50-1 and supplemented with either 100 ml of supernatant with crowding cues or with M9 buffer (control). The crowding cue or M9 buffer was added directly to the OP50-1 lawn in two installments of 50 ml, letting the liquid dry after each installment.

Dauers (fated to become hermaphrodites) were isolated separately on either the “crowding conditions” or the “control conditions” plates and allowed to develop into adulthood at 20°C. Hermaphrodite mothers were moved to a new plate (under the same conditions) each day and egg collection was carried out >3 d (Fig. S1 in Supplementary Materials). Eggs were placed individually on standard (non-treated) plates and the sexual morph was identified as described above.

Figure S1

We built a generalized linear model (GLM) in R using the “glm” function to test if the proportion of females and hermaphrodites depends on the experimental conditions. The data were fitted to a binomial distribution.

Pictures of ten males, hermaphrodites, and females were taken on day 1 of adulthood with a Zeiss Axio Zoom V16 microscope, using a microscopy camera Axiocam and processed with the Zeiss ZEN2 software. Pictures were stored and edited in TIFF file formats. All pictures were taken at 100× magnification.

Measurements of whole-body adults of each sexual morph were taken using the ImageJ “Measure” function. Nematode length was measured, head to tail, and the scale was defined using image scale bars and the ImageJ “Set Scale” feature.

A time course following ten females and ten hermaphrodites over the first 3 d of adulthood was also conducted using the same microscope settings stated above. Measurements of whole-body adults were taken in two parts: from the tip of the head to the anus and from the anus to the tip of the tail. This avoided measurement errors in pictures where the worm had an extremely curled tail. We analyzed the data by performing a two-way repeated measures ANOVA to evaluate simultaneously the effect of the sexual morph and the age on the total body length variable. A repeated measures ANOVA was necessary as the same individuals were measured for three consecutive days. The normality of the data was assessed by performing Shapiro–Wilk tests. We performed pairwise comparisons (t-tests) to further analyze the differences between age and sex.

Typological characters of the male tail, which are considered to represent species-specific characters, were observed and micrographed using a light microscope Eclipse Ni (Nikon) facilitated with DIC optics and a digital camera MC170 HD (Leica) attached to the microscope. Dauer larvae were individually cultured on small (ϕ = 40 mm) Petri dishes and allowed to develop into self-reproducing hermaphrodites. Newly emerged F1 males were picked up using a stainless steel insect pin (Insect pin #00, Shiga Kontyu), and observed and micrographed using the silicon grease method (Kanzaki, 2013). Micrographs were edited with PhotoShop Elements 2021 (Adobe) to construct the figure plate.

The type material of each sexual morph was prepared according to the Natural History Museum of London specifications. Ten individuals of each sexual morph were placed in a microcentrifuge tube with 200 ml of 80% ethanol.

An Illumina whole genome sequencing (WGS) pair-end library of insert-size 200 bp and three mate-pair Illumina WGS libraries of insert-sizes 3 kb, 5 kb, and 8 kb were sequenced using a HiSeq2000 sequencer (Table S1 in Supplementary Materials). Raw reads were preprocessed using Skewer, version 0.2.2 (Jiang et al., 2014) to remove poor-quality sequence regions (<20) and reads smaller than 51 bases (Table S2 in Supplementary Materials). A preliminary genome was assembled to estimate the libraries’ insert sizes and assess contamination. We used Blobtools, version 1.1.1 (Laetsch and Blaxter, 2017) to visualize and assess contamination levels (Fig. S2 in Supplementary Materials). Our libraries contained little contaminants (<2% of the trimmed reads), which predominantly corresponded to E. coli. We assumed it was the strain OP50-1 as this strain was used as the food source of the nematode culture. The reads mapping to the genome of OP50-1 (GCA_000176815.1, ASM17681v1) were removed and not used for the final genome assembly (Table S2 in Supplementary Materials). Reads were error-corrected using Fiona, version 0.2.10 (Schulz et al., 2014) using 60 Mb as an estimated genome size.

Figure S2

An Illumina paired-end RNA-seq library from a mixed population was sequenced using an Illumina NovaSeq 6000 sequencer (Table S1 in Supplementary Materials). Mature RNA was selected using an oligo-dT method. Raw reads were preprocessed using Skewer, version 0.2.2 (Jiang et al., 2014) to remove poor-quality sequences.

Repeats were identified and soft-masked using a combination of programs: TransposonPSI (http:// transposonpsi.sourceforge.net/), LTR_harvest (Ellinghaus et al., 2008) and LTR_Digest (Steinbiss et al., 2009), LTR_finder (version 1.0.7) (Xu and Wang, 2007), RepeatModeler (version 2.0.1) (Flynn et al., 2020) (https://www.repeatmasker.org/RepeatModeler/) and RepeatMasker (version 4.0.9-p2) (https://www.repeatmasker.org/).

The masked genome was then annotated using BRAKER2 version 2.1.6 (Brůna et al., 2021), which uses evidence-based data in conjunction with the ab initio predictors GeneMark (Brůna et al., 2020) and Augustus version 3.4.0 (Stanke et al., 2006). The first run of BRAKER2 was performed using transcriptomic data (RNA-seq reads and de-novo assembled transcriptome). The RNA-seq data generated in this study were pre-processed using Skewer (version 0.2.2) (Jiang et al., 2014) and the trimmed reads were used directly as gene hints by BRAKER2. We also assembled a de novo transcriptome from the trimmed RNA-seq reads using Trinity (version 2.10.0) (Grabherr et al., 2011). The second run of BRAKER2 used the curated protein database “metazoa_odb10” database from orthoDB (Grabherr et al., 2011; Kriventseva et al., 2019) as the evidence-based data. The results of both runs were merged using TSEBRA (Gabriel et al., 2021).

The phylogenetic position of A. melissensis was determined in relation to other Clade V nematode species: C. elegans (WBcel235, GCF_000002985.6), Oscheius tipulae (GCA_013425905.1), Pristionchus pacificus (el paco v4, GCA_000180635.4), A. rhodensis (GCA_947366455), and A. freiburgensis (Talal Al Yazeedi, pers. comm.), using BUSCO single-copy orthologs determined through the Busco Phylogenomics script of Jamie McGowan (https://github.com/jamiemcg/BUSCO_phylogenomics). Briefly, after running BUSCO version 5.2.2 (Simão et al., 2015) on all genomes (Table S3 in Supplementary Materials), single-copy BUSCOs found in at least four species were used to construct alignments using MUSCLE (Edgar, 2004). The alignments were trimmed with TrimAl (Capella-Gutiérrez et al., 2009) and used to construct a phylogenetic tree with either a supertree or a supermatrix approach. For the supertree approach, maximum likelihood trees were constructed for each BUSCO family using the IQ-Tree (Minh et al., 2020). The resulting BUSCO phylogenies were concatenated with ASTRAL, version 5.7.8 (Zhang et al., 2018) to generate the final supertree species’ phylogenetic relationships. For the supermatrix approach, we constructed the phylogenetic tree using Fasttree2 (option “-pseudo,” version 2.1.11) on the concatenated trimmed alignments.

Auanema melissensis n. sp. (Figs. 1, 2A, and 3) (=Auanema sp. JU1783) typologically fits the generic character of Auanema (Kanzaki et al., 2017). Therefore, only body size and some diagnostic characteristics are described here. Whole body length differed between the sexual morphs at day 1 of adulthood: males were small (437.2 ± 18.4 mm, n = 10) and females (854.0 ± 77.0 mm, n = 10) were found to be smaller than hermaphrodites (963.1 ± 64.0 mm, n = 10) (Wilcoxon Rank-Sum exact test, W = 11; P-value = 0.002089; Table S4 and Fig. S3 in Supplementary Materials). A time course following 10 females and 10 hermaphrodites >3 d revealed that the size difference between the non-male sexual morphs was maintained later in adulthood (repeated measures ANOVA F(1, 8) = 26.580, P < 0.05, Fig. 2).

Figure 1

Figure 2

Figure 3

Figure S3

We also observed a difference in pigmentation between females and hermaphrodites, which accentuated during the time course: hermaphrodites seemed to accumulate pigments in the gut whereas females may effectively lose them (Fig. 2B). The bursa is anteriorly open, supported by eight pairs of papillae (rays) arranged as follows <GP1, GP2, (GP3, CO), GP4, GP5d, GP6, (GP7d, GP8) phasmid>, where the distances between GP1-GP2, GP2-GP3 and GP6-GP7 are similar to each other, the distance between GP5d-GP6 is clearly larger and that between GP3-GP4 and GP4-GP5d is shorter; GP7d and GP8 are close to each other (Fig. 3).

In addition to its generic characters, A. melissensis n. sp. is characterized by its male tail characters described above. The new species described here is typologically identical to A. rhodensis, i.e., the male tail characters are shared by these two species (Kanzaki et al., 2017). These two species can be distinguished only by molecular phylogenetic status and mating experiments.

The species was originally isolated from a rotting starfruit (Averrhoa carambola) on September 21, 2009 on the Indian Ocean island of La Réunion in Saint Benoit (Melissa domain). The species epithet is after its type locality, Melissa domain.

Crosses between A. melissensis females and A. rhodensis or A. freiburgensis males (as well as the reciprocal crosses) did not result in progeny (Table 1), indicating that A. melissensis can be considered a distinct species from the two other known trioecious Auanema nematodes.

Similar to the other Auanema nematodes, all dauer larvae of A. melissensis developed into hermaphrodites (n = 48). To determine if the female versus hermaphrodite decision in A. melissensis is mediated by overcrowding cues as in A. freiburgensis (Zuco et al., 2018; Robles et al., 2021), we exposed hermaphrodite mothers to crowding cues. These cues were extracted from plates with high population densities of A. melissensis. Mothers exposed to the crowding cues produced a higher proportion of hermaphrodite progeny than under the control conditions (GLM P-value <0.01; Fig. 4, Table S5 in Supplementary Materials).

Figure 4

The genome of A. melissensis was sequenced using a combination of pair-end and mate-pair Illumina short reads (see Section “Material and Methods”). The best assembly obtained was 59.7 Mb long and estimated to be 88.9% complete using the program BUSCO (version 5.2.2) with the Nematoda_odb10 database (Table 2). This BUSCO score was very close to the score obtained for A. freiburgensis (89.1%) and A. rhodensis (89.8%) (Table S3 in Supplementary Materials), which indicates that A. melissensis’ genome contains most genes. Statistics on the A. melissensis genome are summarized in Table 2. Repeat annotation estimated that 4.8 Mb (8%) of the genome was repetitive (Table S6 in Supplementary Materials).

We used the BUSCO output from A. melissensis to look at the profile of Nigon elements across the scaffolds. Most scaffolds with BUSCO genes corresponded to only one Nigon element (190/224, or 84.8%). Of the scaffolds containing BUSCO genes of several Nigons, 26 (11.6%), 6 (2.6%), and 2 (0.8%) had a mix between 2, 3, and 4 Nigon elements, respectively. The data was provided as an Excel data file in the attachment.

Using the Nigon element concept, it was possible to identify putative scaffolds of the X chromosome. We considered putative X scaffolds containing at least three orthologs pertaining to the Nigon X element (Table 3).

The phylogenetic position of A. melissensis was determined relative to two other Auanema species, A. rhodensis (Tandonnet et al., 2019) and A. freiburgensis (Talal Al Yazeedi, pers. comm.), as well as O. tipulae (Gonzalez de la Rosa et al., 2021), C. elegans (WBcel235) and P. pacificus (el paco v4, GCA_000180635.4).

Auanema melissensis is placed as a sister taxon to A. freiburgensis, with A. rhodensis as an outgroup (Fig. 5). This phylogenetic position is the same as that previously reported, which was constructed based on rDNA and RNA polymerase II sequences (Kiontke et al., 2007; Kanzaki et al., 2017; Winter et al., 2017).

Figure 5