Meiotic maps of sockeye salmon derived from massively parallel DNA sequencing

SNP discovery in parents

The parents from a single family, HX13-WL, were selected for SNP discovery. Over 3.5 million RAD tag sequences were generated in both the male and female parents (Table2). These RAD tag sequences as well as the RAD tag sequences from the offspring (see below) were deposited in the NCBI short read archive (Accession: SRA051991.1). These sequences were grouped and counted resulting in more than 250,000 unique sequences in each parent (Table2). A frequency histogram of the number of occurrences per unique sequence (Figure1) shows a sharp peak of sequences that only occur between one and five times in either parent. We filtered these low count sequences and grouped those remaining into loci shared between the two parents (see methods). This resulted in a set of 64,613 shared loci (Figure1, Table2) that contain both monomorphic and polymorphic RAD tags. Of these, 61,183 loci were monomorphic, 120 loci were homozygous in each parent but polymorphic between the parents, 1,596 loci were polymorphic in one parent and homozygous in the other and 1,714 were polymorphic in both parents. Only loci that were polymorphic within one or both parents could be mapped, resulting in 3,430 putative SNPs to be analyzed further (Table2, Figure1). The location of each SNP in the RAD tag was relatively evenly distributed (Figure2), indicating that the SNPs were unlikely to represent sequencing errors due to drop off of base quality near the end of sequence reads[11].

Individual	Filtered reads	Unique sequences	Shared loci	Putative SNPs
HX13 Female	3,995,897	265,492	64,613	3,430
HX13 Male	3,517,798	258,051	64,613	3,430

Filtered reads is the number of sequence reads after each sequence was trimmed to 70 bp; sequences with less than 80% chance of being error free were discarded from further analysis. Unique sequences is the number of unique RAD tag sequences found in the data from each parent using Perl scripts. Shared loci is the number of unique loci, shared between the two parents; this includes both polymorphic and monomorphic loci. Putative SNPs is the SNP containing loci from the shared loci.

Genotyping of offspring

Sequencing of RAD tags was carried out on 96 offspring from family HX13-WL (Table1, Table3). Two rounds of sequencing were necessary to achieve sufficient depth of coverage to accurately call genotypes in the offspring. The first round was carried out on the Illumina GAII, and the second round on the Illumina HiSeq 2000 (Table3). The combined sequencing efforts resulted in between 1,443,900 and 6,672,291 filtered reads per individual with a mean value of 3,754,867 ± 1,007,970. RAD tags from each offspring were aligned to the 3,430 putative SNPs scored in the parents, and coverage counts for each allele were calculated. Coverage values ranged from 0 to 2,246 reads with an average of 49.0 aligned reads. After converting aligned read counts to genotype calls, the majority of offspring (88/96) were successfully genotyped at over 90% of the putative SNPs discovered in the parents (Figure3).

Illumina sequencing technology	Number of individuals per lane	Number of lanes
GAII paired end	6 parents	1
GAII single end	16 progeny	6
HiSeq 2000 single end	24 progeny	3

Parents were sequenced with the GAII (paired end). Offspring were first sequenced with the GAII (single end); the depth of coverage was insufficient for genotyping and a second round of sequencing was done using the HiSeq 2000 (single end). The results of the GAII single end and HiSeq 2000 single end runs were filtered and trimmed to the same lengths and combined for analysis.

In addition to RAD-based SNPs, we genotyped EST-based SNPs using 127 existing 5’-nuclease assays on parents and 45–141 offspring from each of the thirteen families (Table1). After filtering, each family contained between 44 and 57 segregating loci, with a total of 100 loci segregating across all families (Table2). Linkage relationships in these families were merged into the final map.

Both RAD and EST genotypes between parents and offspring were compared for inheritance error (genotypes in the offspring that do not correspond to the parental types), and erroneous genotypes were converted to missing data. Next the genotypes were screened for segregation distortion. The salmonid genome contains an ancient whole genome duplication and consequently contains many duplicate regions[41]. The presence of these duplications has resulted in high false positive rates in previous SNP discovery efforts. In our data set, markers with significant segregation distortion (χ² test, p < 0.05, deviation from 1:1 or 1:2:1 segregation patterns) may represent these paralogous sequence variants (PSVs), true distortion due to other genomic factors, or genotyping errors. Our genotyping method currently does not allow identification of genotyping error versus true distortion, thus we exclude all markers with significant distortion. Out of 3,530 putative SNPs, 1,758 (49%) showed significant segregation distortion and were removed from further analysis. The majority of these loci showed allelic distributions that indicated that they were likely PSVs (majority of all individuals heterozygous). All of the loci removed for segregation distortion were from the RAD tag sequences. Finally, after removing markers with greater than 25% missing data, 1,725 RAD and 5’-nuclease based SNPs were used in linkage analyses.

Meiotic map construction

Linkage mapping was carried out in two steps. First, the genotypes from the 53 EST loci segregating in family HX13-WL were combined with the RAD tag genotypes, and mapping was carried out on the combined RAD-EST dataset from family HX13-WL. Second, mapping was done on the EST loci that were segregating in the remaining families to maximize the number of ESTs that could be included in the final map (Table1). Linkage groups identified among these EST loci were merged with the HX13-WL linkage maps to create male and female consensus maps.

The female RAD-EST map from family HX13-WL contained 1,050 markers distributed among 29 linkage groups, consistent with the expected number of haploid chromosomes for sockeye salmon[42]. The total map length was 4,943.7 cM. Twenty three markers failed to map and were excluded from further analysis. Additionally, this map contained two small fragments, each containing only two markers, which were discarded from further analysis.

The male RAD-EST map initially contained 28 linkage groups and a total length of 4,352.8 cM. Visual inspection of individual linkage groups and comparison to the female map revealed that one male linkage group corresponded to three female linkage groups (female linkage groups LG27, LG28, and LG29). Inspection of the recombination frequency matrix for the corresponding male and female linkage groups revealed that the male linkage group contained two linkage groups joined through an apparently spurious linkage relationship due to a single marker. No corresponding relationship could be observed in the female. This marker was removed, and the two linkage groups were separated. One of the divided male linkage groups matched to female LG28 and the other matched to the remaining two female groups (LG27, 29). Additionally, one female linkage group corresponded to two male linkage groups, LG7a and LG7b. LG7a contained only five markers and may represent a fragment of the larger LG7b. A number of markers within LG7a had LOD >4 with markers in 7b, suggesting that these are part of a single group. After corrections, the male RAD-EST linkage map contained 1,112 markers distributed among 29 linkage groups and eight small fragments containing only two or three markers. These fragments were discarded from further analysis. Sixteen markers failed to map. The male map length was 4,310.3 cM.

There were 29–35 segregating loci per parent in the EST data set from the additional families (Table1). As a result, only small linkage groups were established, consisting of between two and four markers per group with up to nine groups per family. Overlaps of two or more markers between groups among the families allowed some groups to be merged. After merging, 31 linkage relationships in the males and 28 in the females, consisting of between two and five markers each, were established. These fragments were compared to the RAD-based map from family HX13-WL and merged where two or more markers overlapped. In this fashion, thirteen additional EST loci were added to the consensus map, seven additional markers in the female map and six additional markers in the male. The remaining EST groups were excluded from further analysis.

Merging the EST markers did not alter the number of linkage groups. The consensus female map contained 1,057 markers in 29 linkage groups and a total length of 4,895.8 cM and the consensus male map contained 1,118 markers with a total length of 4,220.0 cM across 29 linkage groups (Figure4, Additional files1,2,3,4,5,6 and7).

There were 457 overlapping loci between the male and the female maps (40% of total female markers and 41% of total male markers) (Figure4, Additional files1,2,3,4,5,6 and7). In all cases the overlapping loci were assigned to the same linkage groups. The male map had two linkage groups (LG7a and LG7b) that corresponded to a single linkage group (LG7) in the female map, resulting in 28 unique male groups when compared to the female map. Examination of the recombination matrix for the corresponding female linkage group did not support splitting this group. Additionally, as noted above, male LG27 corresponds to female LG27 and LG29.

A number of linkage groups contain minor order discrepancies between the male and the female maps (Additional files1,2,3,4,5,6 and7). The markers contained these groups consistently in even when the LOD score is varied. Removal of each marker and repositioning using try.seq in both the male and the female was unable to resolve the order variation in these groups. In no cases did shared markers map to different linkage groups between the male and the female. Linkage groups containing marker order discrepancies include LG14, LG20, LG22, and LG27.

Recombination in the male was reduced over what was observed in the female. The male map contained 61 more markers than the female but was 675.8 cM shorter overall. As not all markers are shared between the male and female, it is inappropriate to make a direct comparison between corresponding linkage groups. Rather, we compared the distance between the most proximal and distal shared pair within each group. The ratio for individual linkage groups ranged from 0.71:1 to 1.71:1. Additionally, the male had more extensive clusters of non-recombinant markers across linkage groups (Table4, Figure4, Figure5 Additional files1,2,3,4,5,6 and7).

Linkage Group	Female cM	Male cM	Female: Male
1	42.7	25	1.708
2	119.5	123	0.97154
3	79.3	96.9	0.81837
4	153.5	160.7	0.9552
5	72	76.6	0.93995
6	91.3	72.8	1.25412
7a	6.9	6.3	1.09524
7b	57.1	57.6	0.99132
8	71.2	81.3	0.87577
9	197.4	224.1	0.88086
10	155.4	126.2	1.23138
11	98.8	99.5	0.99296
12	113.2	123	0.92033
13	144.4	131.4	1.09893
14	126.9	138.9	0.91361
15	151.3	120.2	1.25874
16	80.9	84.1	0.96195
17	39.8	44.3	0.89842
18	120.9	115	1.0513
19	51.5	50.5	1.0198
20	156.8	212.8	0.73684
21	72.1	77.8	0.92674
22	189	176.2	1.07264
23	113.1	159.6	0.70865
24	104.2	115.9	0.89905
25	76.9	74.5	1.03221
26	53.3	46.8	1.13889
27	90.7	82.4	1.10073
28	71.3	91.3	0.78094
29	22.7	22.2	1.02252
Total	2924.1	3016.9	0.96924

The distance between the most proximal and distal of the shared markers was calculated for the male and female on each linkage group and a ratio was compared from these distances. Male LG27 corresponds to female LG27 and LG29; the distances between the shared markers from LG27-F and LG29-F are both included. LG7 is split in the male and 7a and 7b are the shared markers from each of the male linkage groups.

Comparison to rainbow trout

Paired-end sequencing

Paired-end sequencing of genomic DNA from six parents, to provide longer contigs for annotation, produced between 3,517,798 and 3,995,897 filtered reads per individual for assembly. Out of 6,860 potential contigs, 5,722 were successfully assembled with a minimum length of 150 bp. Average contig length was 211 ± 30 bp, with an average coverage value of 38.25 ± 33.76 reads.

After removing duplicate sequences, including the second allele in heterozygous loci, 3,124 contigs remained. A BLASTX search of these contigs against the Swiss-Prot database identified 669 sequences that match a database sequence with an e-value of 10^-4 or less. After assigning GOslim terms, the sequences span a range of biological processes and cellular functions (Figure6A, Additional files1,2,3,4 and5 ). Of the 669 annotations, 115 could not be assigned GOslim terms. While there were no specific patterns in the distribution of GOslim terms, the two most common processes are transport and transcription. Among the sequences with assigned GOslim terms the most common biological processes are cellular transcription and transport, followed by multicellular organismal development. In cellular function, DNA binding and protein binding were the most common terms with 18%, and 17% respectively (Figure6B). Close inspection of the DNA binding category revealed numerous transposable elements. Additional transposable elements were found in the annotated sequences lacking GOslim terms. In total, out of 669 annotated sequences, 83 were transposable elements or related to transposition (12%). These were the most common annotations found in our data (the next most frequent is zinc finger proteins (17/669).

Animals

Thirteen pairs of sockeye salmon from Lake Aleknagik, Alaska, were collected during August, 2009, and used to create 13 full-sib families. Axillary fin clips were taken from all adults and preserved in ethanol for DNA extraction. All embryos were incubated, first in Alaska, and later at the hatchery facility at the University of Washington. After hatch, juveniles were reared in recirculating aquaria for one month at which time 100 individuals from each family were sampled as whole fry and stored in RNAlater (Ambion, Inc., Foster City, CA). Remaining individuals from each family were stored in ethanol. All animal care and use was carried out using methods approved by the International Animal Care and Use Committee (IACUC), under approved protocol 4229–01.

DNA preparation

DNA extractions were carried out on all parents and offspring using DNeasy-96 kits (Qiagen, Valencia, CA) according to the manufacturers’ protocol. DNA was extracted from tail clips from 45–96 offspring from each family and fin-clips from all parents (Table1). Concentration of the extracted DNA was assessed using the Quant-iT PicoGreen dsDNA Assay (Life Technologies, Carlsbad, California) on a Victor D plate reader (PerkinElmer, Waltham, Massachusetts).

5’-nuclease SNP Genotyping and Parentage test

All individuals were genotyped for 127 previously determined EST-based SNPs[64, 69–71] following procedures for 5’ nuclease assays as described in Seeb et al.[72]. Genotype calling was carried out using the SDS 2.3 software (Life Technologies, Carlsbad California) or the BioMark 3.0.2 software (Fluidigm, South San Francisco, California). Monomorphic assays and assays which failed to amplify in the families were excluded from further analysis. The remaining 5’-nuclease genotypes were successfully filtered for Mendelian inheritance error and segregation distortion as described below, and markers uninformative for mapping were removed. Occasional family-to-family contamination had been previously observed in our fish incubation and rearing system.

RAD library construction and sequencing

RAD library preparation was carried out on 96 offspring from one family and six parents (Table2) using the methods similar to those previously described[39, 40]. The parents were prepared for paired end sequencing to provide longer contigs for annotation and 5’nuclease assay design (see below); the offspring were prepared for single end sequencing. The methods are identical except for the size of the sheared library extracted from the gel, 400–800 bp for single end sequencing and 150–400 bp for paired end sequencing. The restriction enzyme SbfI was used to digest the DNA, and SbfI specific Illumina linkers each containing a unique barcode were ligated to each digested DNA sample as described in Miller et al.[11]. Individual samples were pooled into libraries containing up to 16 individuals. Each library assessed for quality and concentration using a Bioanalyzer DNA 1000 kit (Agilent Technologies, Santa Clara, California), and final concentration was determined by the Bioanalyzer software (Agilent Technologies, Santa Clara, California).

Ten nanomoles of each library in EB (Qiagen, Valencia, California) and 0.1% Tween 20 were sequenced at the University of Oregon High Throughput Sequencing Facility. The parents were sequenced using paired-end technology at 2 × 80 bp on the Illumina Genome Analyzer II (GAII). All offspring were first single end sequenced at 80 bp on the GAII. Another round of sequencing was required to achieve satisfactory depth of coverage; the second round was conducted on the Illumina HiSeq2000 (Table3), which produces 100 bp reads.

SNP discovery

All sequence analysis for parents and offspring was carried out using Perl scripts (Perl scripts used are available from Miller et al.[11] ) and the publically available software programs Bowtie version 0.12.7[73] and Novoalign version 2.07 (http://www.novocraft.com). First the sequences were filtered for quality. The last ten base pairs from each sequence were removed from GAII sequences (paired end and single end); subsequently any reads with less than 80% chance of being error free were discarded from further analysis. Sequences obtained from the HiSeq2000 were trimmed to 70 bp by removing the last 30 bp (to match the GAII sequences), and any reads with less than 80% chance of being error free were discarded. The 80% threshold was determined using the Phred quality score of each read. The quality value at each nucleotide was converted to a probability score and summed across the read length and if the sum fell below 80% the read was discarded. All quality filtered data was de-multiplexed into separate sequence files for each individual using the barcodes contained in the Illumina linker ligated during library preparation (see above)[11]. These procedures were applied to all GAII single end, GAII paired end, and HiSeq2000 data. The data for each offspring from the GAII and the HiSeq2000 was combined for all further analysis. SNP discovery was then carried out using only parental sequence data.

Using Perl scripts, the set of unique RAD tag sequences for each parent was identified, and the count of occurrences for each unique sequence within each parent was obtained. After labeling each sequence with the parent of origin, the set of sequences from both parents were combined for alignment in Novoalign. The combined set of RAD tag sequences were self-aligned in Novoalign. The software was set to carry out an exhaustive search, returning up to 40 alignments per unique sequence. We set a maximum alignment score of 125. The results of this alignment were filtered to maximize true (non-paralogous) loci shared between parents. These loci included monomorphic and SNP containing loci, as well as some paralogs. The loci were filtered using methods adapted from Miller et al.[11]. In order to filter the sequences, the script compares each alignment score to a threshold value, compares the number of reads at each alignment, and finally compares the number of alignments within and between each parent. The Perl script to identify loci takes into account alignment score, coverage, and the number of identified alignments both between and within parents. For this study a minimum alignment score of 30, with no more than one alternative allele in either parent was used. The results of this search were filtered to include only the sequences containing putative SNPs, and these were output in a FASTA file for use as a reference sequence for Bowtie alignment.

SNP genotyping in offspring

The sequences containing putative SNPs in the parents were used as a reference set in a Bowtie alignment, formatted as “allele 1” and “allele 2” for each putative SNP. Quality filtered sequences from all offspring and the two parents were aligned to this reference, with Bowtie set to only align reads which matched the reference perfectly. The output of this alignment was filtered using custom Perl scripts, counting the number of matches to each allele for each individual. These allele counts were exported for conversion into genotype calls in Excel.

SNP genotype calls from the allele counts were made using methods similar to those in Nielsen et al.[59]. A threshold of 10 reads from both alleles (sum of both alleles) was required for genotyping; any locus with less than ten reads was called as missing. This threshold was determined experimentally by testing various threshold values on the data set and examining the resulting calls for inheritance error, segregation distortion, and the number of missing markers, and subsequently comparing these values between thresholds. During SNP discovery each allele was designated as allele 1 or allele 2. For all loci, allele 1 was designated as the reference allele, thus in some cases the alternative allele may have higher read counts than the reference allele. For each locus, a heterozygous genotype was called if the read count for the non-reference allele was between 28%-80% of the total count for both alleles. Otherwise a homozygous genotype was called. Certain combinations of allele counts may still result in a miscalled allele, however, such a miscall in the parents would result in a mismatch in alleles between parents and offspring. All data was screened for such inheritance error (see below) and markers with an excess of these errors were removed from the analysis.

Next, parent and offspring genotype calls were compared to check for inheritance error, where the genotypes of the progeny do not match the expected genotypes based on the parental genotypes. If found, these genotype calls were converted to missing data. The informative genotype calls from existing sockeye 5’-nuclease assays were added to the data set. The data was checked for segregation distortion from 1:1 or 1:2:1 ratio using a χ²test. Markers exhibiting significant levels of segregation distortion were removed from the analysis. Finally, markers for which one parent or more than 25% of offspring genotype calls were missing were removed, and the genotypes were converted into OneMap format[56, 74].

Linkage mapping

Linkage analysis was carried out using the R package OneMap version2.0-1[74]. OneMap uses a maximum likelihood analysis to assign phase to identical heterozygous progeny genotypes, allowing the inclusion of these markers in the linkage map[56]. Male and female linkage maps were generated separately. All linkage map ordering was carried out using the Kosambi map function[8, 75, 76]. Data was imported in outcross format. Next, the recombination fraction between all pairs of markers was calculated using the rf.2pts algorithm. Linkage groups were formed using the group command in two stages. First, linkage groups were formed with LOD = 12 and a maximum recombination fraction of 0.01. These groups were ordered using either the compare function or the order.seq function depending on the number of markers. Groups with up to seven markers were ordered with the compare function, groups containing more than seven were ordered with order.seq. This identified markers which overlap in position. One marker from each position with the fewest missing genotypes was selected to use in the next step of map construction. Next, linkage groups were formed with LOD = 6 and the default recombination fraction. Order.seq was used to place markers in their initial order. After initial ordering, ripple.seq was used to test the marker order, with a sliding window of four markers, and the order was changed when necessary. The marker orders were also inspected visually using rf.graph.table, which plots a heat map of the LOD score and recombination frequency, for ordering errors too far apart to detect with ripple.seq and made manual changes when necessary. Finally, the markers removed in the LOD = 12 analysis were individually added back to each linkage group using the try.seq and ripple.seq functions to finalize marker order. The male and female linkage maps were compared for marker order among the shared markers. Where discrepancies were found the individual marker was removed using drop.seq, and then added back to the map using try.seq to examine optimal marker position. In cases where no consensus order could be established this is noted in the text.

Linkage relationships were obtained among the segregating 5’-nuclease loci in each of the 12 families and in both the RAD tag data and the 5’-nuclease assays in the sequences family. The individual male and female maps from each family were combined where possible using MERGE, part of the LINKMFeX package[77].

Comparison to rainbow trout RAD tags

All putative loci discovered in the HX13-WL parents were compared to the rainbow trout RAD tags in Miller et al.[11] using BLASTN (BLAST version 2.2.25). For unmapped loci, only full length (59 bp) hits with no more than one mismatch were considered. For loci used for either map, the stringency was stepped down and the map locations were compared. If hits occurred in one to two linkage groups, consistent with syntenic matches, they were kept. In this way, hits longer than 56 bp and containing no gaps and no more than four mismatches were considered matches. Microsoft Excel was used to identify shared tags that were present in both rainbow trout and sockeye salmon linkage maps. Tag positions were compared between the rainbow trout and sockeye linkage maps and overlapping map hits were counted.

Paired End assembly

Paired-end sequence data from the parents was assembled to generate longer sequences for use in annotation of RAD tags and as a resource for 5’ nuclease assay design. Paired-end assembly was carried out using custom Perl scripts and the publically available program Velvet version 1.1.06[78]. The Perl script matches each sequence in the FASTA file generated during SNP discovery, containing the 6,860 SNP alleles, against the forward end (paired-end 1) paired-end file and pulls both matching sequences and their pairs and writes them into a temporary data file. Both alleles of each SNP discovered were used as the template for the paired-end sequence assembly. The complete set of paired, filtered reads from all six individuals was compared to the SNP-containing FASTA file. The temporary file is then run through Velvet, and the resulting output is added to the results file. This procedure is repeated iteratively until all the sequences have been compared. Velvet is set with a k-mer value of 21 and a minimum contig length of 150 bp. All other Velvet settings are set to default values. For any contig not meeting the 150 bp minimum, a placeholder value is inserted to indicate the missing sequence. After assembly, duplicate contigs were removed.

Contigs from the paired-end assembly were subjected to a BLAST search against the Swiss-Prot database using BLASTX. An e-value of less than 10^-4 was required to accept the annotation. Microsoft Access and the complete set of GO terms for the Swiss-Prot database (NCBI version 2011_11) were used to associate sequence annotation with their GOSlim term. Terms were assigned for both biological process and cellular function.

The 5’-nuclease assays were designed as described previously ([67, 69, 70]). Five assays were designed on assembled paired end sequence, and five assays were designed from RAD tags alone.