SNP discovery in European lobster (Homarus gammarus) using RAD sequencing

S3 Supplementary Material – Bioinformatics

Quality filtering

Quality of the raw reads was initially examined using the FastQC software (Babraham Bioinformatics). The sequencing data provided by the Exeter Sequencing Service had already been demultiplexed, which meant fastq files were received for each sample (two files per sample because of paired-end sequencing). The use of phased adapters meant there were one of four different combinations at the start of each forward read: TGCA, ATGCA, CATGCA or GCATGCA (nucleotides marked in red represent the beginning of the SbfI cut site). As a result, the raw data was trimmed at the start of each read, such that each read started with the SbfI cut site (TGCA). Raw reads were then cleaned and truncated to 97 bp using the process_radtagsprogram in the Stacks software v1.45(Catchen et al. 2013). A read was discarded if the score fell below a 90 % probability of being correct (a raw Phred score of 10). Reads were truncated because 97 bp was the shortest read in the dataset and Stacks requires all reads to be the same length. No barcode argument was provided because the data had already been demultiplexed.

Building loci denovo

The formation of RAD loci was carried out by running the wrapper script denovo_map.plin Stacks. The minimum number of reads required to form a stack (m; ustacks), the number of mismatches allowed between stacks to merge them into a locus (M; ustacks) and the number of mismatches allowed between stacks during the construction of the catalog (n; cstacks) parameters were optimised following the methods ofParis et al.(2017). Each parameter was sequentially changed while keeping the other two constant, and the highest number of r80 polymorphic loci was considered the optimum parameter set. The r80 polymorphic loci are identified by adjusting the r parameter in the populations program, whereby a locus must be present in at least 80 % (r = 0.80) of the individuals in a population for it to be processed. All other parameters were set to their default values during the optimisation of m, M and n. Once the optimum values m, M and n were obtained, denovo_map.plwas run again with the maximum number of stacks at a single locus (max_locus_stacks, ustacks) set to 2, meaning only loci with two alleles (biallelic loci) were allowed. This parameter was enforced because, for a diploid species, if more than two alleles are allowed during the denovo assembly of loci, the locus assembly could contain merged paralogs or additional alleles introduced as a result of sequencing error(Catchen et al. 2013; Mastretta-Yanes et al. 2014). Individuals were dropped from the analysis if their coverage was lower than 15x.

SNP discovery

SNPs were discovered by running the populationsprogram with two sets of parameters.The first run of the populations program was executed using all 55 samples with no population map, a minimum allele frequency set to 0.1, a maximum observed heterozygosity set to 0.7 and r was set to 0.8. The rationale for this run was to search for SNPs that can detect genetic differentiation across broad geographical areas. The second run of populations searched for SNPs that could potentially detect genetic differentiation (if any) between geographical locations within the Atlantic Ocean. In this run, the Mediterranean and Skagerraksamples were removed because early genetic analyses indicated differentiation between these samples and all other Atlantic samples (Fig. S3.1); accordingly, a population map for only the Atlantic sampleswas submittedand organised by geography (Table S3.1). Since most populations were only composed of two or three samples, the r parameter was set to 1, meaning a locus had to be present in every individual in that population for it to be processed. Similarly, the p parameter was set to nine (the total number of populations in the population map), which meant that a locus must be present in all populations to be eligible. All other parameters remained the same as the first run.

Figure S3.1: Discriminant analysis of principal components (Jombart et al. 2010). Sample locations were used as priors and the optimal alpha score was used to assess the number of principal components to retain.

Table S3.1: Population map composed of nine putative populations and 40 individuals.

Population / Samples
English Channel / Ald3, Ald4, Jer2, Jer3, Sbs2, Sbs3
West Ireland / Ara1, Ara2, Don2, Don3, Ven1, Ven3, Mul3, Mul4
Southwest England / Bs55, Bs67, Ios5, Ios6, Lo31, Loo18
Southeast Ireland / Cor1, Cor2, Hoo1, Hoo2
North Sea / Cro4, Cro5, Nh8
Irish Sea / Iom1, Iom4, Pem12, Pem14
Orkney / Shetland / Or6, Ork1, She14_3, She14_4
France / Lr4, Lro11, Lro13
Spain / Vg1, Vg3

SNP panel development

To develop the SNP panel, a set of criteria were implemented to filter out uninformative SNPs and SNPs that were inadequate for primer design, whilst retaining those SNPs which were most informative and which were suitable for primer design. Firstly, SNPs derived from populations were ranked and sorted by highest G”st(Meirmans & Hedrick 2011) using the diff_statsfunction in the mmod R package (Winter 2012). The top 300 SNPs were then exported to an Excel file whereby information about each SNP was recorded (e.g.the position of the target SNP in the RAD-tag locus, the RAD-tag locus ID named by Stacks and the ranking according to highest G”st).

Secondly, the locus ID for each SNP was queried using the Stacks SQL user interface to find out which loci had additional non-target SNPs in the RAD-tag. The SNP database was sorted by the fewest number of SNPs per locus, followed by the highest ranking. The ideal RAD-tag locus contained only one SNP as this was advantageous for primer design. However, loci with two or more SNPs were considered, particularly if they were ranked highly, because high-throughput assay designs can often deal with non-target SNPs assuming they are at least 30 bp away from the target SNP.

Thirdly, mini-contigs were assembled by aligning paired-end reads to the original RAD-tag sequence to try and extend the length of each RAD-tag locus. To do this, paired-ends sequences were collated for each locus by executing the sort_read_pairs.pl program in Stacks using a whitelist of locus IDs. Paired-ends reads were then assembled using the exec_velvet.pl program (Velvet wrapper script) with the M parameter set to 100, generating paired-end consensus sequences at least 100 bp long. To build mini-contigs, the paired-end consensus sequences were then aligned to the original RAD-tag sequence (97 bp length) using the alignment tool in Geneious v.10.1.3(Kearse et al. 2012). SNPs with less than 25 bp flanking sequences or non-target SNPs within 30 bp either side of the target SNP were discarded. Finally, the sequences (mini-contigs) for the remaining SNPs were then submitted to Fluidigm’s online portal for in silico assay design. This performs a quality control analysis on each sequence submitted, checking that SNP assays will be compatible with each other (i.e. no primer dimers) and that parameters set by Fluidigm are met (e.g. no regions of large repeats, suitable GC content, adequate flanking sequence, etc.).

References

Catchen J, Hohenlohe PA, Bassham S, Amores A, Cresko WA (2013) Stacks: an analysis tool set for population genomics. Molecular Ecology, 22, 3124–3140.

Jombart T, Devillard S, Balloux F et al. (2010) Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC Genetics, 11, 94.

Kearse M, Moir R, Wilson A et al. (2012) Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics, 28, 1647–1649.

Mastretta-Yanes A, Arrigo N, Alvarez N et al. (2014) Restriction site-associated DNA sequencing, genotyping error estimation and de novo assembly optimization for population genetic inference. Molecular Ecology Resources, 1–14.

Meirmans PG, Hedrick PW (2011) Assessing population structure: FST and related measures. Molecular Ecology Resources, 11, 5–18.

Paris JR, Stevens JR, Catchen JM (2017) Lost in parameter space: a road map for stacks. Methods in Ecology and Evolution, 8, 1360–1373.

Winter DJ (2012) MMOD: An R library for the calculation of population differentiation statistics. Molecular Ecology Resources, 12, 1158–1160.

1