Identifying transgene insertions in Caenorhabditis elegans genomes with Oxford Nanopore sequencing

Line description

Both BY250 and UA44 strains were created using vectors inserted into the C. elegans Bristol N2 background (Adams, 2022). Specifically, in the BY250 line, neurons express Green Fluorescent Protein (GFP) after the integration of the pPdat-1::GFP construct into the genome (created by Dr. Randy Blakely, Florida Atlantic University;  Nass et al., 2005). BY250 is a well-used model to explore the impact of neurodegeneration following chemical exposure from 6-hydroxydopamine and rotenone (Nass et al., 2005; Ray et al., 2014). The UA44 transgenic line was created at the University of Alabama (Cao et al., 2005) through coinjection of a plasmid with GFP and a second plasmid with human α-synuclein. When the two plasmids are coinjected into C. elegans they form stable extrachromosomal concatamers (Mello et al., 1991). We verified this using phenotype analysis where dopaminergic neurons, highlighted with GFP, displayed α-synuclein dose-dependent neurodegeneration as the animals aged (Cao et al., 2005; Hamamichi et al., 2008). Stably expressing worms were then chromosomally integrated to create UA44 (Qiao et al., 2008; Kritzer et al., 2009).

DNA extraction and sequencing

Nematodes were grown to large population size on NGM plates seeded with E. coli OP50 (Stiernagle, 2006). Plates were washed with M9 media into 15mL conical tubes and rocked on a table rocker for an hour to purge any biological waste. The tubes were then centrifuged to pellet the worms and remove the supernatant, and then the pellet was washed with M9 approximately five times until the supernatant was clear to remove as much bacteria and waste from the worm pellet as possible. The pellet was then moved to a two mL tube and frozen at −20 °C until extraction. Genomic DNA was extracted with a modified phenol-chloroform extraction following Sutton et al. (2021).

We used the ONT SQK-LSK109 ligation sequencing kit for library preparation with a modification, replacing the first AmpureXP bead clean with a treatment of the Short Read Eliminator Kit available from Circulomics Inc. Approximately 700 ng of gDNA were loaded onto a R9.4.1 RevD flow cell on the ONT GridION X5 platform and sequenced for 48 h. We performed standard base calling using Guppy v.4.0.11 and trimmed adaptor sequences from DNA reads using Porechop (− −discard_middle; Wick, 2018). We used Nanoplot to calculate the read statistics after Porechop, and on the subset of reads used in the final Flye assembly (see assembly details; De Coster et al., 2018).

Short-read sequence data

In order to polish the assembled DNA libraries, we downloaded Illumina paired-end DNA libraries from the NCBI Sequence Read Archive (SRA) in April 2020 (BioProject PRJDB2670; NCBI Resource Coordinators, 2017). In order to polish the sequences containing the transgene, we simulated 150 bp, 100x coverage paired-end DNA libraries for the insertion sequences with the software package ART (versionMountRainier; Huang et al., 2012). The simulated data was was added to the respective Illumina data downloaded from SRA and used to polish the completed genome assemblies with the software Pilon (described below; Pilon v1.23; Walker et al., 2014).

Genome assembly, polishing, and scaffolding

We corrected raw ONT DNA libraries using the Canu –correct option, which corrects errors based on best overlaps among reads (Canu v1.9; Koren et al., 2017). The Canu-corrected reads were then assembled using Flye (v2.8.1; Kolmogorov et al., 2019). We used our paired-end short-read dataset containing the respective short-read libraries and simulated transgene insertion libraries to polish our draft assemblies four times with Pilon and eliminate small base pair errors and insertions/deletions (Walker et al., 2014) Caenorhabditis DNA libraries are frequently contaminated with foreign DNA (Fierst & Murdock, 2017; Fierst et al., 2017). To decontaminate our DNA libraries, we created a reference database for taxonomic assignment of contiguous sequences (contigs) with blastn using the NCBI ‘nt’ database (v2.2.31; Altschul et al., 1990). Only contigs aligning to Caenorhabditis were kept. Draft assemblies were then aligned, corrected, and scaffolded according to the genomic position within the C. elegans N2 reference genome (WBcel235) using RagTag (v2.0.1; Alonge et al., 2022). During the RagTag − −correct step we used the −j flag which allowed us to “hide” a list of query sequences with hits to the respective transgenes to prevent the module from removing the insertion due to lack of similarity to the surrounding C. elegans reference sequence. Six small sequence fragments under 1 kb were removed for final statistics and depositing with the NCBI.

Transgene insertion identification

We aligned the vector insertion sequence for each assembled genome strain to the draft genome assembly before and after scaffolding using minimap2 (Li, 2018; Li, 2021). These locations were used to inform the correction step of RagTag to keep the insertion from being removed from the corrected reads. Additionally the insertion vector sequence was aligned to the final genome using blastn, to identify all stop and start locations for the insertion in each genome (Altschul et al., 1990). Annotated figures were made with gggenomes (Hackl, 2022).

Quality assessment

The quality of the assemblies was assessed and compared to the C. elegans N2 reference genome (WBcel235) using QUAST (v5.1.0rc1; Gurevich et al., 2013). BUSCO (v5.3.2; database: Nematoda_odb10) was used to assess the completeness of our genomes using a unique set of 3,131 genes expected to be conserved in a single copy within Nematoda (Simão et al., 2015; Manni et al., 2021). Both quality assessments were measured before and after scaffolding with RagTag (v2.0.1; Alonge et al., 2022).

Annotation

Gene annotation was completed using the protein-coding sequences from C. elegans N2 (WBcel235) as a reference with the software Liftoff which is able to map annotated protein-coding sequences from one assembled genome to another when the two strains are closely related (Shumate & Salzberg, 2021). After annotation, the coding sequences (CDS) were extracted from the assembled genome sequence using the ‘extract_sequences’ function from the AGAT package (Dainat, 2020), and a BED coordinates file was created from the GFF using the ‘gff2bed’ function from BEDOPS (Neph et al., 2012). We used the WormBase Parasite BioMart (Howe et al., 2017) to extract identifying features associated with genes that were present in the C. elegans N2 reference protein-coding annotations but missing in our genome annotations. These features included orthologous and paralogous genes, chromosomal locations, gene ontology (GO) terms (The Gene Ontology Consortium, 2000) and Interproscan protein domain annotations (Zdobnov & Apweiler, 2001; Jones et al., 2014; Finn et al., 2017).

Synteny analysis

The mRNA sequences and BED coordinates files generated during the annotation step were used to assess synteny with MCscan (python version) from the JCVI utility library (v1.1.17; Tang et al., 2008; Tang et al., 2017). We compared the UA44 and BY250 genome sequences to the standard C. elegans N2 (WBcel235) genome limiting our comparison to the six largest linkage groups representing the six Caenorhabditis chromosomes.

Sequence data summary

After sequencing, the Oxford Nanopore DNA library for the UA44 strain contained 7.44 Gb of sequence reads with an N50 (the size median) of 4,732 bp with an average (mean) sequencing depth of 74.2x (Table 1). The BY250 library contained 12.45 Gb of sequence reads with an N50 of 6,116 bp with an average depth of 124.12x. The reads were then corrected with Canu –correct (Koren et al., 2017). The Canu-corrected UA44 DNA library was 3.88 Gb with an N50 of 8,082 bp with an average depth of 38.69x (Table 1). The Canu-corrected BY250 library was 4.10 Gb with an N50 of 15,109 bp and an average depth of 40.84x. Additional ONT library statistics can be found in Table S1.

Assembly

The Canu-corrected reads were used to assemble both genomes with Flye (v2.8.1; Kolmogorov et al., 2019). The UA44 draft assembly contained 89 contigs with a length of 102.8 Mb (Table 2), and the BY250 draft assembly contained 118 contigs with a length of 102.5 Mb. After scaffolding with RagTag (v2.0.1; Alonge et al., 2022) and removing contigs under 1 kb, the UA44 assembly contained 55 contigs, and the N50, an assessment to measure the length of the scaffold that represents 50% of the genome, improved from 3.6 Mb to 18 Mb. The BY250 assembly contained 46 contigs after scaffolding, and the N50 improved from 1.8 Mb to 17.9 Mb. The scaffolded N50s of 18 Mb and 17.9 Mb for UA44 and BY250 respectively are comparable to the N50 of 18 Mb found in the reference assembly WBcel235. Additionally, the L90 for UA44 improved from 28 contigs to six, and from 57 contigs to six for BY250. With final L90s of six contigs and L50s of three contigs, our assemblies match those of the WBcel235 genome, indicating that correction and scaffolding improved the UA44 genome to near chromosome level. L90 represents the number of contigs that make up the first 90% of the length of the genome, and L50 represents the number of contigs that make up 50% of the length. Both genomes also show a consistent GC% of 35.47% matching the 35.44% of the reference. The final UA44 genome assembly represented 99.65% of the standard reference genome, and the final percentage match for BY250 was 99.44% (Table 2). These metrics indicate that our genome assemblies with reference-based scaffolding are very complete and near chromosome-level.

To assess the completeness of our genomes, we used BUSCO (v5.3.2) to search our genomes for a set of 3,131 orthologous genes present in nematodes (Nematode odb v10; Simão et al., 2015; Manni et al., 2021). The final assembly for UA44 showed a 98.9% BUSCO score matching the score of the reference (Table 3). UA44 had 21 missing genes, which is also the same as the C. elegans reference. BY250 had a final BUSCO score of 98.6% with 24 missing genes, only three more than the reference. Over 97.9% of the genes were present in single copy in both genomes, very similar to the 98.3% of single copy genes in the reference. Our gene completeness results show that our genomes compare very well to the standard published C. elegans reference genome.

Annotations

Annotations were mapped to the genome using Liftoff (Shumate & Salzberg, 2021). Results are shown in Table 4. Of the 20,191 genes present in the reference, 20,100 and 20,099 were lifted over to UA44 and BY250 respectively representing over 99.5% of the reference genes. Only 91 and 92 genes were not able to be mapped during liftover to UA44 and BY50. Of those 91 and 92 missing genes, 77 were shared between the UA44 and BY250 annotation. 46% of the genes missing in BY250 were located on Chromosome V in the C. elegans N2 reference annotations, 20% on Chromosome X, 10% on Chromosome II, 9% on Chromosome IV, 8% on Chromosome III and 8% on Chromosome I. UA44 missing genes had a similar chromosomal distribution (within 1–2%) to BY250.

Inspection of missing genes on Chromosome V of BY250 revealed two common patterns. First, 13 missing genes had no orthologous or paralogous genes, Gene Ontology (GO) (The Gene Ontology Consortium, 2000) or Interproscan protein domain annotations (Zdobnov & Apweiler, 2001; Jones et al., 2014; Finn et al., 2017). These genes were also very small for protein-coding genes, ranging in size from 20–36 nucleotides. Computational protein-coding genome annotation has a high false positive rate (Bruna et al., 2021) and it is possible these gene sequences were not accurate.

The second pattern we identified was missing genes that were verified in C. elegans and the result of complex patterns of molecular evolution. For example, the Seven TransMembrane Receptor gene str-120 is part of a gene family that contains 21 orthologues across Caenorhabditis and 256 paralogues within the C. elegans genome. Similarly, the EB1 C-terminal domain-containing protein ebp-3 is part of a gene family that contains 158 orthologues but just two paralogous genes in C. elegans. Caenorhabditis genomes show rapid rates of gene family expansion and shrinkage through tandem duplications and deletions (Adams et al., 2023), leaving closely related complete and partial gene sequences that can be challenging for alignment and annotation (Figs. S1B, S2B). Overall, our coding sequence annotations (>99.7%) and mRNA annotations (>99.6%) were highly consistent with the reference genome. Additionall y > 99.8% of the annotated genes were contained on the six main linkage groups (I, II, III, IV, IV, X) or the mitochondria.

Transgene insertion

Prior to scaffolding, the vector insertion sequence for UA44 aligned to Contig59 which then scaffolded into Chromosome IV. For BY250 the vector insertion sequence aligned to Contig86 which then scaffolded into Chromosome I. The BY250 and UA44 transgenes mapped to expected genomic regions based on outcomes from genetic crosses with these transgenes (Gaeta et al., 2023). These locations were used to exclude Contig59 and Contig86 from the RagTag (Alonge et al., 2022) correction step of scaffolding for their respective genomes.

We also used blastn to align the vector sequence to the genomes (Altschul et al., 1990). For UA44, the vector insertion was found in four complete copies on Chromosome IV. The locations of the 4,924 base pair vector were: IV:8767911-8763011 (−4,900 bp), IV:8772965-8768045 (−4,920 bp), IV:8777841-8772966 (−4,875 bp), and IV:8784656-8789078 (4,422 bp) (Table 5). Three of the locations were negatively oriented, while one was positively oriented (Fig. 1, created with gggenomes; Hackl, 2022). Additional partial hits of the insertion vector sequence to the genome were found, however, none of these partial hits contained the α-synuclein gene (Table S2). The four-copy insertion of α-synuclein into Chromosome IV of UA44 may increase the overall expression of the genes present in the insertion. Parkinson’s severity in humans shows copy-number dependence including within-family genomic triplication of SNCA (the α-synuclein locus) (Singleton et al., 2003; reviewed in Cognata et al., 2017). Importantly, annotations of the surrounding areas indicate that the four-copy insertion did not affect genes in that region (Fig. 1, created with gggenomes Hackl, 2022). Some of the partial vector insertions either overlapped or were adjacent to genic sequences (Table S2). These intersections may be potentially capable of altering the regulation of proteins or gene expression.

For the BY250 transgene, only two BLAST (Altschul et al., 1990) hits were found. The primary insertion location was on chromosome I:11058905–11058020 and was 885 bp of the 1,113 bp insertion vector sequence (Table 5). A second partial insertion hit was also found on chromosome I:11059035–11058951, but was only 84 base pairs long. The BY250 insertion found in single copy on chromosome I did not intersect any surrounding genes (Fig. 2, created with gggenomes Hackl, 2022).

Synteny

Coding sequences and gene BED coordinates were used to align each genome to the C. elegans N2 (WBcel235) reference genome using MCscan (Python version) to assess synteny (Tang et al., 2008; Tang et al., 2017). Macrosynteny visualization shows our genomes covering the seven chromosomes of the C. elegans N2 reference (Fig. 3). Macrosynteny plots for each assembly aligned to the reference along with syntenic maps and alignment depth comparisons are available in the Figs. S1 and S2.