A draft genome of a field-collected Steinernema feltiae strain NW

S. feltiae strain NW was collected using trap insect technique with modification (see Kaya and Stock, 1997) from a potato field near Patterson, WA under organic management. Specifically, a trap consisted of 10 wax moth, Galleria mellonella, larvae (Nature’s Way, Ross, OH) in a sealed mesh bag was buried 15 cm deep in the soil, and the trap was retrieved 48 hr later. Nematodes were extracted from dead larvae using White traps (White, 1927). A week later, infective juveniles (IJs) were collected and kept at 12°C. To reduce heterozygosity, 20 IJs were used to establish an inbred colony. These IJs were used to infect a single wax moth, and seven days after infection, IJs were recovered using a White trap (White, 1927). Approximately 200 emerged IJs were used to infect 10 wax moths for two more times. Emerged IJs were washed three times in sterilized tap water and stored in −80°C until DNA purification.

Approximately 10,000 IJs were thawed and frozen twice. High molecular weight genomic DNA was extracted using a phenol chloroform protocol described in the study of Donn et al. (2008). The DNA pellet was re-suspended in 100 µl 10 mM Tris-Cl (pH 8.5) buffer. Washington State University’s Genomic Core Facility (Pullman, WA) performed library preparation and sequencing. Three batches of sequencing with different chemistries were conducted: XL-XL, P4-C2, and P5-C3. In total, 42 SMRT cells (24 cells of XL chemistry, 7 cells of P4-C2 chemistry, and 11 cells of P5-C3 chemistry) were included in the sequencing.

We used bash5tool.py algorithm in the SMRT analysis pipeline (v2.3) to extract subreads from all 42 SMRT cells with minReadScore set to 0.75. Subreads were used for genome assembly. We tested two genome assembly programs initially: Canu (v1.3; Koren et al., 2017) and Celera (v8.3rc2; Myers et al., 2000). For both assemblers, we started with the default setting. Canu outperformed Celera in running time and memory use. Canu also generated an assembly with fewer contigs and greater N50. Thus, we picked Canu for further fine-tuning to improve assembly quality. For the final assembly, we used sprai (v0.9.9.18, http://zombie.cb.k.u-tokyo.ac.jp/sprai) to correct errors for the subreads, and we imported corrected reads into Canu to perform de novo assembly.

Once we obtained an assembly with the best continuity and lowest number of contigs, we removed contigs that had low read support (<20 reads). The read coverage information was indicated in Canu (v1.3; Koren et al., 2017). Often, contigs with low read support were highly fragmented and not informative in building gene models during the annotation process. Next, we polished the assembly by first aligning raw PacBio-h5 files to the assembly using the program pbalign following consensus calling using Quiver (SMRT analysis pipeline (v2.3)). Finally, we queried each contig in the assembly to genome sequences of the endosymbiont X. bovienii (GenBank assembly accession: GCA_000027225.1) using BLASTN (Altschul et al., 1990). Contigs showed homology to any sequence of X. bovienii with bit score >200 were removed in the assembly.

MUMmer (v4.09; Marçais et al., 2018) was used to align our assembly to the genome assembly of SN. Program dnadiff in MUMmer was used to summarize the comparative analysis of the two genomes. Program nucmer provided a 1-to-1 and all-vs-all comparisons of the two genomes, and this analysis was limited to the 100 longest contigs/scaffolds to minimize noise.

Because the genome assembly of strain NW seemed to be divergent from the previously published genome of strain SN, we annotated this genome to present as an alternative resource for the research community. To do this, we modeled and identified potential repeat region of the S. feltiae NW genome based on a repeat library kindly provided by WormBase using RepeatModeler (v1.0.11, http://www.repeatmasker.org/RepeatModeler/). After that, we used RepeatMasker (v4.0.9; Smit et al., 2013-2015) to soft mask the genome. Two RNA-seq data sets (BioProjects accession numbers PRJNA488495 and PRJNA283145, Chang et al., 2019), were aligned to the NW genome assembly using program HISAT2 (v2.1.0; Kim et al., 2015) with the default setting. The resulting bam files from HISAT2 were used as RNA-seq evidence to train gene models in BRAKER2 (Hoff et al., 2015). The generated annotation file in GFF3 format from BRAKER2, along with the protein FASTA file of strain SN (Dillman et al., 2015), were incorporated into the program MAKER (v2.31.10; Cantarel et al., 2008) to synthesize gene models. BUSCO (v3.02; Waterhouse et al., 2017) was used to evaluate the completeness of the genome based on a set of universal single-copy orthologs across Nematoda.

The whole genome assembly was archived in NCBI with WGS accession number MQUG00000000. PacBio subreads were deposited in SRA under BioProject PRJNA353610. Genome in GFF3/GTF format along with predicted mRNA, CDs, and protein sequences were deposited in WormBase ParaSite and will be available from release 15 and forward.

In total, 7,166,944 reads (read length ranging from 500 bp to 39,605 bp, median=1,429 bp, mean=1,681 bp) were generated from 42 SMRT cells across three chemistries. Reads generated in our study were sufficient to cover the genome 145-fold, assuming the genome size of 82 Mb as previously described (Dillman et al., 2015). Such high coverage would adequately overcome the high error rate often associated with long reads generated from PacBio platforms (Rhoads and Au, 2015).

The final assembly used in the downstream analysis was generated in Canu (Koren et al., 2017) with an error rate=0.03, and raw reads covered the genome assembly 99-fold. After removing contigs that were likely contaminated with X. bovienii, and contigs with low coverage support (<20 reads), the final assembly contained 4,876 contigs with N50 (contigs) of 60,433 bp. The assembled genome size of NW is significantly larger (47%) than SN with noticeable improvement regarding the number of contigs and continuity (N50, Table 1). In addition, our assembly of NW contained no ‘gaps’ (i.e., unknown sequences between two contigs of which order and orientation are known; Clum, 2018) in comparison to numerous gaps totaling 2,769,616 bp in the SN genome; this would facilitate annotation of the genome (e.g., Dallery et al., 2017).

We performed a comparative genomic analysis of our NW assembly with the earlier SN assembly using a whole genome alignment approach. Over 90% of the contigs/scaffolds from the two genomes aligned, totaling 94 Mb and 74 Mb from the NW and SN assemblies, respectively. In total, 24 percent (29 Mb) of the NW sequences failed to align to the SN assembly, indicating divergence between the two strains. Genome alignment determined that there were 49,949 1-to-1 best alignments between these two assemblies with an average identity of 97%, totaling 72.5 Mb of the sequences from both genomes. There were 133,667 many-to-many alignments with a similar average identity of 96%. There were 907,706 SNPs detected in common between the two genomes in addition to numerous indels, breakpoints, and allocations (Table 1). The majority of the top 100 longest contigs/scaffolds from the two assemblies aligned with the presence of insertions and rearrangements. We also noted that a few contigs/scaffolds failed to align (absence of dots at the top right of Figure 1). Despite our improvement of the assembly, neither of these two genomes has assembled to the chromosome level.

Figure 1:

Because our genome seemed to be quite divergent from the assembly of strain SN with many structural variations and SNPs, we proceeded to annotate the genome with evidence generated from RNA-Seq alignment. Despite the larger assembled genome size, our annotation predicted rather fewer (12% fewer genes) number of genes/proteins compared to the previous annotation of SN genome (Table 1). However, the total length of amino acids in our annotation is longer than SN, suggesting a longer or greater number of full-length genes. Furthermore, our NW genome had a slightly higher percentage of completeness compared to the SN genome (Table 1).