A toolkit for rapid gene mapping in the nematode Caenorhabditis briggsae

Background The nematode C. briggsae serves as a useful model organism for comparative analysis of developmental and behavioral processes. The amenability of C. briggsae to genetic manipulations and the availability of its genome sequence have prompted researchers to study evolutionary changes in gene function and signaling pathways. These studies rely on the availability of forward genetic tools such as mutants and mapping markers. Results We have computationally identified more than 30,000 polymorphisms (SNPs and indels) in C. briggsae strains AF16 and HK104. These include 1,363 SNPs that change restriction enzyme recognition sites (snip-SNPs) and 638 indels that range between 7 bp and 2 kb. We established bulk segregant and single animal-based PCR assay conditions and used these to test 107 polymorphisms. A total of 75 polymorphisms, consisting of 14 snip-SNPs and 61 indels, were experimentally confirmed with an overall success rate of 83%. The utility of polymorphisms in genetic studies was demonstrated by successful mapping of 12 mutations, including 5 that were localized to sub-chromosomal regions. Our mapping experiments have also revealed one case of a misassembled contig on chromosome 3. Conclusions We report a comprehensive set of polymorphisms in C. briggsae wild-type strains and demonstrate their use in mapping mutations. We also show that molecular markers can be useful tools to improve the C. briggsae genome sequence assembly. Our polymorphism resource promises to accelerate genetic and functional studies of C. briggsae genes.


Background
Comparative analysis of developmental and behavioral processes in closely related species is a powerful approach to understand the mechanisms of evolution. It facilitates identification of molecular components that are conserved over millions of years due to their role in specifying common features as well as those that are variable because they confer species-specific features. The model organism Caenorhabditis elegans (a nematode) and its congener, C. briggsae, are particularly suitable for such investigations. Their many experimental advantages include rapid growth, small size, transparency, ease of culture and genetic manipulations, and the availability of fully sequenced genomes [1][2][3].
C. briggsae is phenotypically almost indistinguishable from C. elegans and has a similar (hermaphroditic) repro-ductive mode. The last common ancestor of these two species lived about 30 million years ago [4], and despite the rapid molecular evolution typical of the family Rhabditidae, more than half (~52%) of the C. elegans genome aligns with the C. briggsae genome assembly [2]. This includes two-thirds of all C. briggsae genes (13,107 or 67.8%) with reciprocal orthologs in C. elegans [5]. Thus C. elegans-C. briggsae comparative genomic and genetic studies promise powerful new tools for the identification of genes and pathways and the study of both conservation and divergence.
Like C. elegans, C. briggsae has six chromosomes that display extensive conservation of synteny, but not exact colinearity relative to C. elegans [6]. While C. briggsae shares many of the experimental advantages of C. elegans, it has the further advantage of increased natural variability for single nucleotide polymorphisms (SNPs) and insertion-deletions (indels) [7,8]. This elevated natural variation potentially enhances its use for genotype-phe-notype association studies, and is also very useful for the mapping aspects of forward genetics projects.
Among the tools needed to facilitate forward genetics in C. briggsae, a set of easily scored DNA polymorphisms is especially important. Experimentally validated polymorphisms can serve as useful markers for mapping mutations that cause visible phenotypes. Additionally, these markers can be integrated with the phenotypebased genetic linkage map (e.g., dpy and unc mutants [21]) to further enhance their utility. Integration of polymorphisms and phenotype-based maps increases map density and anchors the relative locations of molecular and phenotypic markers. With this goal in mind we have discovered a large set of genome-wide polymorphisms (SNPs and indels) in wild-type strains, using AF16 as a reference strain and four other natural isolate strains: HK104, HK105, VT847, and PB800.
The indels were placed into three classes: small (7-49 bp), medium (50-2,000 bp), and large (>2 kb). We have focused on medium and small indels (212 and 7,530, respectively), which offer the greatest utility as genetic markers. In the case of SNPs (23,829), we found that 4,700 modify restriction enzyme sites (termed snip-SNPs) and therefore can be easily detected as restriction fragment length polymorphisms (RFLPs). We established assay conditions for bulk segregant analysis (BSA) and used these to experimentally validate 14 snip-SNPs, 28 medium and 32 small indels. The validated polymorphisms were used to genetically map known mutations causing visible phenotypes thus demonstrating the effectiveness of the polymorphisms in linkage mapping studies. We also developed single animal-based PCR assay to determine map distance. Five mutations were successfully localized to sub-chromosomal regions by 3 or more indels, greatly facilitating the search for each candidate gene. These results demonstrate the utility of our mapping toolkit in genetic linkage and gene identification studies.

SNP Discovery
We performed SNP discovery in four C. briggsae strains by aligning paired shotgun sequence reads to the AF16based reference sequence (cb25 assembly, [22]). These sequences were obtained by capillary gel electrophoresis at Washington University Genome Center (see Methods). To build on previous SNP discovery efforts [6], we applied the ssahaSNP algorithm, which detects SNPs and small indels based on SSAHA alignments to a reference sequence (see Methods). Compared to AF16, ssahaSNP detected 23,829 unique SNP loci in HK104 DNA, or one substitution per 163 bp on average (Table 1 and additional file 1). Consistent with C. briggsae clade structure [8], SNP density was slightly lower in strains HK105 (1/ 168 bp) and PB800 (1/197 bp) and much lower in strain VT847 (1/475 bp). In HK104, the most common substitu- Total aligned base pairs include redundant matches due to sequence overlaps (between 15% and 25%) in sequence data. The SNP density is based on the number of uniquely aligned base pairs. tion by far was A(T) to G(C), which accounted for 57.1% of all substitutions ( Figure 1A).

RFLP genotyping assays and validation of snip-SNPs
We screened the SNPs predicted for HK104 for variants that altered the recognition site of a restriction enzyme, and thus might be amenable to restriction fragment length polymorphism (RFLP) genotyping. To make this a practical resource, we limited the analysis to 30 restriction enzymes from REBASE [23] that are reliable and inexpensive. Of 23,829 HK104 SNPs, some 4,700 (19.72%) were predicted to alter the recognition site of at least one of the 30 restriction enzymes. To develop restriction fragment length polymorphism (RFLP) assays from these snip-SNPs, we designed PCR primers with a standard protocol and performed in silico digests of the resulting amplicons to infer the banding patterns for each strain. RFLPs not easily distinguishable on a gel, or SNPs on ultracontigs not yet included in the genetic map, were removed. Finally, we used assembly AGP information and BLAST alignment to obtain coordinates for each snip-SNP on the cb3 sequence assembly. Our set contained 1,987 predicted RFLP assays from 1,362 snip-SNPs (some SNPs alter multiple RE sites) positioned on both the genetic and physical maps (see additional file 2). Another snip-SNP, bdP3, was identified in a separate study (see Methods and additional files 1 and 2).
We selected a total of 20 RFLP assays (between 3 and 4 for each chromosome) based on HindIII, DraI and SalI snip-SNPs for validation in AF16 and HK104 parental DNA ( Table 2). Roughly a third of the assays (6) failed PCR in one or both strains in repeated attempts. Although we did not investigate the issue of PCR failure, it is possible that redesigning primers (by moving them out or in) and testing different PCR conditions may produce desired products in some cases. All of the 14 assays successfully gave rise to strain-specific RFLP banding patterns, validating the predicted snip-SNP ( . Interestingly, two of these assays (cb55670 and cb20723) exhibited HK104 fragments that varied from in silico predictions, another possible consequence of unknown variants in this highly divergent strain. Consistent with C. briggsae clade structure [8], VT847 was not polymorphic (from AF16) for the snip-SNPs we examined.

Indel Discovery
The ssahaSNP program [24] was also able to detect insertion/deletion (indel) polymorphisms of 1-50 bp ( Figure  1B). We used the parse_indel utility to extract 7,530 candidate indels (4,686 deletions and 2,844 insertions compared to AF16) for the HK104 strain. Most insertions and deletions detected by ssahaSNP were single base pair events; the largest was 49 bp. To identify larger insertion/ deletion events, we developed a customized algorithm called BreakPointRead that detects indels based on BLAST alignments between read sequences and the reference genome. First, the algorithm identifies "breakpoint reads" with alignment gaps of 10 bp or larger compared to the reference sequence. Next, it analyzes the gap size and alignment orientations to infer the nature (insertion, deletion, inversion, etc.) and approximate size of sequence variation. When we applied BreakPointRead to the HK104 sequence traces, it identified 689 breakpoint reads suggestive of 635 underlying variants. We filtered the results to obtain insertion/deletion events between 50 and 2,000 bp. The resulting set contained 212 predicted indels (144 deletions and 68 insertions), the largest of which was a 1,707 bp deletion on chromosome IV ( Figure  1C and additional file 3).

Medium indels
We used the BreakPointRead algorithm to develop PCR fragment length polymorphism (PLP) assays for medium indels to facilitate high-throughput gene mapping. Of the 212 putative indels, we selected 40 for validation in AF16 and HK104 parental DNA ( Figure 2B, C). Two other indels (bdP1 and bdP4) were also chosen for a similar analysis (see Methods and additional file 3). Ten of the 42 PLP assays failed PCR in one or both strains ( Table 3). As mentioned earlier, some of these errors may be eliminated by redoing PCR using new primers. One assay (cb-m16) resulted in PCR products that were quite large (>1,500 bp) in both strains and therefore could not be accurately resolved on the gel. Of the remaining 31 assays that were successfully amplified, 29 confirmed the pres- The snip-SNPs are arranged by chromosome and location. The corresponding ultracontigs are also listed. The expected and observed DNA fragments refer to products based on in silico analysis and actual experiments, respectively. The bdP3 location is unknown because the corresponding contig is unassembled. The validation type column marks snip-SNPs that were consistent with prediction (*) and those that differed significantly (**). RE: Restriction enzyme used to digest PCR amplified products.
ence of polymorphism, 22 of which were similar in type and size as predicted (Table 3).

Small indels
We developed PLP assays for small (<50 bp) indels in AF16 and HK104 that were identified by ssahaSNP. To allow for gel resolution, we excluded indels smaller than 7 bp. This resulted in 436 assays that had 7-49 bp band size differences between AF16 and HK104 (see additional file 4). We tested 45 indels (between 4 and 9 for each chromosome) by PCR and found that except one (bhP44), for which HK104 amplification failed, all others could be successfully amplified (Table 4). A total of 32 indels showed bands of predicted sizes (Table 4, Figure 2D). Of the remaining 12, 1 showed no indel (i.e., identical PCR products in AF16 and HK104), 4 had multiple products (either due to PCR error, incorrect in silico predictions, or misassembly), and 7 showed PCR products that were inconsistent and unreliable (Table 4). Altogether we experimentally confirmed 75 polymorphisms (14 snip-SNPs, 29 medium indels, and 32 small indels) on all six chromosomes ( Figure 3, Table 5). The utility of these 'working' markers in genetic studies is demonstrated by successful mapping of several mutations that cause visible phenotypes. In two cases these mapping experiments also helped improve the genome sequence assembly. Specifically, the bhP42 contig fpc4184 was placed near the center of chromosome 1 and bhP18 contig fpc4010 was reassigned to the left arm of chromosome 3 (see below).  The indels are arranged by chromosome and location. The corresponding ultracontigs are also listed. The "unknown" locations refer to unassembled contigs. The predicted indel sizes are based on in silico analysis of AF16 and HK104 amplicons. The validation type column marks indels that were consistent with prediction (*) and those that differed significantly (**).

Linkage mapping studies using snip-SNPs and indels
The snip-SNPs and indels were used to map a set of 12 mutations with easily recognizable phenotypes (e.g., Uncoordinated or Unc and Dumpy or Dpy) that were previously isolated in different laboratories (Table 6 and Methods). Except lin(bh25) and unc(sy5415), all other loci were uniquely assigned to linkage groups by standard 2factor mapping using known mutations that serve as reference (Table 6, also see [21]). The dpy(s1272), unc(s1270), and unc(sa997) are reference markers for LGIII, LGIV, and LGV, respectively. The remaining auto-somal loci are linked to C. elegans orthologs Cbr-lin-11 (LGI) and Cbr-unc-4 (LGII). Not only did the polymorphism-based mapping agree with phenotypic markerbased mapping (see dpy(sy5001) and dpy(sy5148) in Figure 4A, B), it also helped to identify linkage groups of lin(bh25) (LGI, see Figure 4D and Table 6) and unc(sy5415) (LGV, see Table 6). In each of these cases a single cross with HK104 provided enough genomic DNA and usually one PCR per chromosome was sufficient to establish the linkage (using bulk-segregant approach, BSA). We also quantified DNA band intensities to deter-   mine linkages as unitless linkage values (ULVs) (see Methods). As expected, for unlinked loci the ULV was one. In the case of dpy(sy5148) the ULV for chromosome 2 indel (bhP21) was 2.7 suggesting a strong linkage (Figure 4C). Overall, these results demonstrate that polymorphism-based mapping can be used to quickly map new mutations in C. briggsae.
In addition to linking mutations to chromosomes we also investigated whether polymorphisms could be used in more precise mapping i.e., placing mutations in spe-cific chromosomal regions (left, right arms, or middle). We reasoned that by narrowing down genetic intervals of mutations it should be possible to identify potential candidates, including C. elegans orthologs, thereby facilitating gene cloning by RNAi and transgene rescue approaches. To this end we used three medium indels on chromosome X to map dpy(sy5001). The ULVs for dpy(sy5001) suggest weak linkages to indels cb-m204 (left arm) and cb-m136 (right arm) and tight linkage to the middle indel cb-m197 ( Figure 5A). Similar result was also obtained with the small indel bhP26 that is located close to cb-m197 and is strongly linked to dpy(sy5001) (data not shown).
Besides dpy(sy5001), we also mapped 4 autosomal mutations to sub-chromosomal regions using small indels. The results showed that dpy(s1272) is most strongly linked to bhP18, unc(sa972) to bhP14 and bhP18, and lin(bh20) to bhP40 (all on chromosome 3) ( Figure  5B). The unc(sy5422) appears to be located on the right arm of chromosome 4 (closer to the middle) since it shows strongest linkage to bhP9 ( Figure 5C).
The C. briggsae sequence assembly 'cb3' had placed the bhP18 contig fpc4010 on the right arm of chromosome 3 (~36.5 mu). We found that this location was inconsistent with ULVs for dpy(s1272) and unc(sa972) since both mutations are also linked to the left arm indel bhP14 (Figure 5B). This suggested that there could be a possible error in the sequence assembly. In a separate study Zhao et al. [25] used the SNP-based oligonucleotide array to map dpy(s1272) and found that the region corresponding to bhP18 is misassembled. Therefore, we have assigned bhP18, as well as dpy(s1272) and unc(sa972), to the very end of the left arm of chromosome 3. Additional mapping using polymorphisms and visible markers will resolve the extent of misassembly.
In addition to the above BSA approach we also analyzed single F2 mutants segregating from a cross (termed single recombinant analysis or SRA) to determine recombination distances between polymorphisms and mutations. For this we used a multivulva mutation lin(sy5353) and an Unc mutation unc(sy5506). The lin(sy5353) mutation is located on chromosome 1 since it is strongly linked to three small indels bhP1, bhP7 and bhP42 (B.P. Gupta, unpublished results). From a cross between lin(sy5353) and HK104 we picked 23 F2 mutant animals (46 chromosomes) and analyzed their DNA for the presence of bhP1. A single recombinant chromosome was recovered suggesting a recombination frequency of ~2% between lin(sy5353) and bhP1. In the case of unc(sy5506) mutation, located on chromosome X, we analyzed 20 F2 mutant animals for the presence of indel bhP26. A total of 4 recombinant chromosomes were recovered (see additional file 5) suggesting that the two loci are 10% apart. These results demonstrate that SRA mapping protocol can be used in C. briggsae to localize mutations to subchromosomal regions and narrow down their genetic interval.

Discussion and Conclusions
We took a bioinformatics approach to identify polymorphisms in the C. briggsae genome and experimentally val-  idated a small set of these to facilitate mapping of mutations. Comparison of AF16 (reference strain) to four other isolates (HK104, HK105, VT847 and PB800) revealed that HK104 is most polymorphic since it has the highest density of SNPs among all the strains. Altogether we identified ~31,300 polymorphisms (23,800 SNPs and 7,500 medium and small indels) between AF16 and HK104 that promise to be a valuable resource for mutation mapping and genome evolution studies. Roughly 20% of the SNPs are predicted to alter restriction enzyme sites (snip-SNPs) that could be detected by PCR followed by restriction digestion and agarose gel electrophoresis.
A total of 107 polymorphisms (20 snip-SNPs, 42 medium indels, and 45 small indels) that were experimentally tested, 66 (12 snip-SNPs, 22 medium indels, and 32 small indels) showed DNA fragments identical (or close) to in silico predictions (Table 5). Another 9 cases (2 snip-SNPs and 7 medium indels) were significantly different but nonetheless showed the presence of underlying variants. In 15 cases no polymorphism could be detected. Thus, excluding PCR failures (total 17), the success rate of correctly predicted polymorphisms was 73% (69-86% range) ( Table 5). This suggests that both ssahaSNP and BreakPointRead algorithms work equally efficiently regardless of the type of polymorphism in question. A similar study in C. elegans [26] showed that greater than 95% of the polymorphisms predicted by the Polybayes program [27] are true. It remains to be seen whether the lower success rate in C. briggsae is due to intrinsic differences between the programs alone or if the quality of sequence data and assembly are additional contributing factors.
We used snip-SNPs and indels to map 12 mutations with visible phenotypes, and found that polymorphismbased mapping agreed with phenotypic marker-based results. Furthermore, it helped map two mutations, lin(bh25) and unc(sy5415), for which no prior genetic linkage data was available. Five mutations were also localized to sub-chromosomal regions. Thus our mapping resource can be used to rapidly map new mutations in C. briggsae. It is also relatively easy to validate additional polymorphisms if one needs a greater resolution. It should be pointed out that Hillier et al. [6] have validated another set of 9 snip-SNPs by sequencing during the process of C. briggsae genome sequence assembly (see additional file 6). Given the high density of such markers (>2,000), it should be possible to map a mutation within a small genetic interval to facilitate molecular cloning (e.g., see [28]).
In addition to mapping mutations, SNPs and indels could also be used to improve the genetic linkage map of C. briggsae. The current C. briggsae sequence assembly, cb3, incorporates 90.2% (91.2 Mb) of the genome united into six chromosomes [6]. The remaining 9.8% of sequences are tentatively associated with chromosomes. These unmapped regions could be integrated into chromosomes by polymorphism-based recombination mapping. We have successfully used this approach to place the contig fpc4184 in the vicinity of fpc3441 (chromosome 1) based on the recombination distance of 5% between bhP42 (fpc4184) and bhP1 (fpc3441) (Figure 3) (A. Seetharaman, P. Cumbo, B. Nagagireesh and B. P. Gupta, manuscript submitted). In the other case, we have reassigned the bhP18 contig fpc4010 to the left arm of chromosome 3 based on its tight linkage to dpy(s1272)  and unc(sa972). Additional snip-SNPs and indels could further refine the locations of these contigs.

In silico predictions of polymorphisms SNPs
SNP discovery was performed on 13,632 shotgun sequence traces from strains HK104, VT847, HK105, and PB800. The ssahaSNP program (version SSAHA2) [24] was used to call SNPs due to its robust and efficient performance; only polymorphisms with quality scores above the minimum threshold were accepted. We also tested Polyphred (v5.04) [30] and PolyBayes (v3.0) [27] programs but found that only ssahaSNP could efficiently handle the entire read set and reference genome sequences as input. For a reference sequence the cb25 genome sequence assembly, which is based on strain AF16 and organized into ultra (fingerprint) contigs, was obtained from Wormbase. Flanking sequences for predicted SNPs were repeat-masked to lower case using the RepeatMasker program (v3. 1.5) [31] with a customized C. briggsae repeat library.
The HK104 SNPs were positioned on the cb25 sequence assembly during SNP discovery. To position them on the newer cb3 sequence assembly, which is by chromosome, we obtained the assembly AGP files from Wormbase. SNP positions were inferred based on the coordinates and orientation of their cb25 ultracontig. SNPs on cb3-unmapped ultracontigs were mapped by WU-BLAST v2.0 (Gish, W., personal communication) alignment of their flanking sequences. Some 699 SNPs could not be positioned on the cb3 assembly by either method.

Indels
Candidate AF16-HK104 indels were extracted from HK104 sequence traces using the parse_indel utility of ssahaSNP. In the HK104 set, the largest indel event identified by ssahaSNP was 49 bp. To identify larger insertion/ deletion variants we implemented BreakPointRead, a custom algorithm that detects structural variations (insertions, deletions, inversions, and copy number variants) spanned by individual sequence traces. Traces were aligned to the reference genome (cb25) using WU-BLAST v2.0 (Gish, W., personal communication), and screened for alignments with "gaps" of > = 10 bp. The alignment patterns of such "breakpoint reads" were used to infer the type and size of polymorphism. Predicted insertions and deletions were set aside for assay development.

bdP polymorphisms
The bdP polymorphisms described in this study (snip-SNP bdP3 and medium indels bdP1 and bdP4) were identified in the laboratory of SEB. The snip-SNP bdP3 was earlier used in a study involving ray pattern variation in C. briggsae [32].

Development of RFLP and PCR Length Polymorphism (PLP) assays
SNPs were screened for substitutions that altered the recognition sequence of restriction enzymes using the Bio::Restriction::Analysis library of BioPerl [33]. The analysis was limited to 30 restriction enzymes from REBASE [23] known to be reliable and inexpensive. PCR assays were designed (amplicon sizes of 500 to 1000 bp, primer Tm's of 54-56 °C) using a local installation of the primer3 program [34]. In silico fragment analysis of the PCR products was performed to predict band sizes for AF16 and HK104; assays with more than 4 bands in either strain were removed.
In the case of small indels (7-49 bp), primers were selected to generate AF16 amplicon sizes within the range of 200 and 400 bp. For medium indels (50-2,000 bp), primers flanking each indel and specifying an AF16 amplicon size of 300-800 bp were selected.

PCR
In all experiments the genomic DNA from F2 worms (derived from a cross between AF16 and HK104 animals) was used as a PCR template. In some control experiments genomic DNA from F1 heterozygous animals was also used. PCR results that gave rise to unexpected or no products were repeated at least twice. In some cases we also tested different annealing temperatures. Those that consistently failed were termed as "PCR failure".

Mutation mapping
We picked 12 strains for linkage mapping studies ( Table  6). The mutations were obtained from EMS (ethyl methane sulfonate) mutagenesis screens in an AF16 genetic background in various laboratories. The strains were outcrossed several times (3 or more). For mapping, mutant hermaphrodites were crossed with HK104 males and the genomic DNA from 20 F2 animals (wild type and mutant separately) was prepared as described in the previous section. The linkage was determined by PCR using protocols established for control experiments.

Linkage and ULV analysis
To determine the linkage of a mutation to a chromosome, we initially relied on the visual inspection of DNA band intensities on Ethidium bromide-stained agarose gels. Subsequently, in indel-based mapping experiments, we calculated linkages as unitless values (ULVs) for an unbiased analysis. The mean intensities of DNA bands were measured by NIH ImageJ software (version 1.41o; [35]) using Measure tool under Analyze menu. For each muta-tion a ratio of band intensities in the "mutant" lane was calculated by dividing the mean intensities of the AF16 bands by the mean intensities of the HK104 bands. This ratio was termed as the ULV. As expected, ULVs were one for unlinked mutations and higher for linked mutations.

Genetic positions of polymorphisms
The genetic positions of snip-SNPs and indels in this study correspond to nearest SNPs that were experimentally validated (D.C.K. and R.D.M., unpublished). These 'verified' SNPs (400 in total) were genotyped in RILs derived from two independent crosses between AF16, HK104 and VT847 (AF16 × HK104 and AF16 × VT847). The details are available on the C. briggsae SNP Research Facility website [36].