{'fraction_non_alphanumeric': 0.03310900780940728, 'fraction_numerical': 0.029258286248965343, 'mean_word_length': 4.981058975462764, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Less is more: extreme genome complexity reduction with ddRAD using Ion Torrent semiconductor technology

Massively parallel sequencing a small proportion of the whole genome at high coverage enables answering a wide range of questions from molecular evolution and evolutionary biology to animal and plant breeding and forensics. In this study, we describe the development of restriction‐site associated DNA (RAD) sequencing approach for Ion Torrent PGM platform. Our protocol results in extreme genome complexity reduction using two rare‐cutting restriction enzymes and strict size selection of the library allowing sequencing of a relatively small number of genomic fragments with high sequencing depth. We applied this approach to a common freshwater fish species, the Eurasian perch (Perca fluviatilis L.), and generated over 2.2 MB of novel sequence data consisting of ~17 000 contigs, identified 1259 single nucleotide polymorphisms (SNPs). We also estimated genetic differentiation between the DNA pools from freshwater (Lake Peipus) and brackish water (the Baltic Sea) populations and identified SNPs with the strongest signal of differentiation that could be used for robust individual assignment in the future. This work represents an important step towards developing genomic resources and genetic tools for the Eurasian perch. We expect that our ddRAD sequencing protocol for semiconductor sequencing technology will be useful alternative for currently available RAD protocols.


Introduction
As the introduction of next-generation sequencing (NGS) technologies, whole genome sequencing has been carried out on an increasing number of species, triggering a major breakthrough in genetics. For many studies, however, obtaining the complete genome sequence from a large number of individuals still remains prohibitively expensive and is not even necessary. Therefore, several genome complexity reduction methods which allow only a subset of the genome to be sequenced have been recently described . These methods rely on the generation of reduced representation libraries (RRL), generated by restriction enzymes (Altshuler et al. 2000), selective amplification (van Orsouw et al. 2007), targeted amplicon sequencing (Hyten et al. 2010), capture probes (Shen et al. 2013), methylation filtering (Pal-mer et al. 2003), high-C 0 t selection (Yuan et al. 2003) or by sequencing the transcribed proportion of the genome, that is transcriptome (Barbazuk et al. 2007).
One of the most widely used genome complexity reduction methods that have been applied in both model (e.g. humans, laboratory mice and Drosophila) and nonmodel species is restriction-site associated DNA (RAD) sequencing. This method targets DNA sequences flanking specific restriction enzyme cutting sites throughout the genome (Baird et al. 2008;Hohenlohe et al. 2010). To date, RAD sequencing has been applied in various research fields ranging from molecular evolution and evolutionary biology (Emerson et al. 2010) to animal and plant breeding (Yang et al. 2012), conservation (Sharma et al. 2012) and forensics (Ogden et al. 2013). Importantly, RAD sequencing is highly suitable for generating genome-wide genotype data in situations where not much is known about the target genome. As a result, RAD sequencing has been increasingly used in population genetic, biogeographic and phylogenetic studies of nonmodel organisms providing a genome-wide view on the role of different evolutionary forces during population differentiation, adaptation and speciation (Baxter et al. 2011;Bruneaux et al. 2013;Keller et al. 2013).
Earlier RAD library preparation protocols typically consisted of the following steps: restriction enzyme digestion, ligation of the first adapter, physical shearing, end-repair, ligation of the second adapter and size selection. However, more recent studies have simplified, and further developed, the library preparation procedure by eliminating physical shearing and introducing two restriction enzyme digestion with strict size selection (termed double RAD (Bruneaux et al. 2013) or doubledigest RAD (Peterson et al. 2012)). These modifications enable greater flexibility in targeting the optimal number of genomic regions ranging from thousands of consensus sequences required for QTL mapping to hundreds of thousands of contigs necessary for association studies and hitchhiking mapping.
Most of the RAD studies to date have been carried out using Illumina GAII and HiSeq NGS platforms, but to fully exploit the potential of RAD sequencing, other NGS technologies can be used. For example, the new semiconductor sequencing technology (Ion Torrent PGM) represents a simple, fast (hours instead of days) and flexible solution for small laboratories and individual research groups, but currently, only a few studies have reported the utility of Ion Torrent PGM technology for RAD sequencing (Mascher et al. 2013;Kai et al. 2014). This can be partially explained by the relatively small sequencing throughput of PGM technology, which produces tens of millions of short reads compared to, for example, hundreds of millions of sequences generated by the Illumina HiSeq platform. However, for RAD projects which require sequencing of a smaller number of genomic regions, such as QTL mapping and basic population genetic analysis, Ion Torrent PGM technology is more accessible to smaller research groups making it a faster and simpler alternative compared to other NGS platforms. Therefore, it is important to further develop and evaluate RAD sequencing protocols adopted for Ion Torrent PGM that enable high genome complexity reduction by targeting a relatively small number of genomic regions while retaining the high coverage necessary for robust SNP discovery and reliable estimation of various population genetic parameters.
In this study, we describe the development of a twoenzyme RAD sequencing approach for Ion Torrent PGM platform by carrying out extreme genome complexity reduction using rare-cutting restriction enzymes and strict size selection of the library while using DNA pooling as a cost-effective approach to estimate allele frequencies on a genome-wide scale. We applied this approach to the common freshwater fish species, the Eurasian perch (Perca fluviatilis L.). Eurasian perch contributes significantly to both commercial and recreational fisheries in the coastal Estonian waters of the Baltic Sea and in Lake Peipus, the fifth largest lake in Europe (Moran 2003;Pukk et al. 2013). Because of different fishing regulations in these water bodies, the ability to use genetic markers for assigning individuals back to their source population would be of considerable importance for the identification of illegal fishing and fish trade. To develop genomic resources and genetic tools for the Eurasian perch, which could be used in a court of law, we: i) developed a restriction-site associated DNA (RAD) sequencing approach for Ion Torrent PGM platform which results in extreme genome complexity reduction allowing sequencing of a relatively small number of genomic fragments with high sequencing depth; ii) characterized a small proportion of the genome and identified more than a thousand single nucleotide polymorphisms (SNPs) in perch; iii) estimated genetic differentiation between the DNA pools constructed based on individuals collected from freshwater (Lake Peipus) and brackish water (the Baltic Sea) environments; and iv) identified SNPs with the strongest signal of differentiation which could be used for individual assignment in the future.

DNA isolation and ddRAD library preparation
Seventy-six wild Eurasian perch (Perca fluviatilis L.) were used for ddRAD library preparation. Genomic DNA was isolated from fin clips and dried scale samples as described by Pukk et al. (2013). These samples were divided into two pools: the Baltic Sea pool (Turku Bay, n = 33; Matsalu Bay, n = 6 and P€ arnu Bay, n = 11) and Lake Peipus pool (n = 26). DNA quality and concentration were assessed by agarose gel electrophoresis and with a NanoDrop spectrophotometer (Thermo Fisher, Inc.). RAD library preparation protocol broadly followed the methods described by Bruneaux et al. (2013), with some modifications as outlined by Pukk et al. (2014). In short, pooled DNA (800 ng) was digested for 2 h at 37°C in 20 lL reaction volume, simultaneously with two restriction enzymes, 20 U of PstI (restriction site 5 0 CTGCAG 3 0 ) and BamHI (restriction site 5 0 GGATCC 3 0 ) (New England Biolabs) and heat-inactivated for 15 min at 75°C (Fig. 1). The ligation consisted of forward (0.016 pmol) and reverse (0.004 pmol) adapters which were added to 20 lL of restriction reaction together with 5 lL 10 9 T4 DNA Ligase buffer and 1200 U of T4 DNA ligase (M0202S; NEB) (forward adapter, top: 5 0 CCATCT-CATCCCTGCGTGTCTCCGACTCAGXXXXXTGCA 3 0 , forward adapter, bottom: 5 0 XXXXXCTGAGTCGGAGA CACGCAGGGATGAGATGG 3 0 , where XXXXX is a bar-code sequence; reverse adapter, top: 5 0 GATCATCA CCGACTGCCCATAGAGAGG 3 0 , reverse adapter, bottom: 5 0 CCTCTCTATGGGCAGTCGGTGAT 3 0 ). To differentiate between the pools, two following 5 bp barcodes were used (Sea pool: 5 0 AGAAC 3 0 and Lake pool: 5 0 TCGTT 3 0 ) (Fig. 1). The 50 lL ligation reaction was carried out at 22°C for 1 h, heat-inactivated for 30 min at 65°C. Each library was subsequently loaded onto E-Gel â SizeSelect 2% Agarose Gel (Life Technology) to extract DNA fragments of approximately 300 bp length (ranging between 200 and 351 bp). Adapter-ligated products were nick-translated and PCR-amplified in 67.5 lL volume containing 14.1 lL of E-Gel extraction product, 50 lL of Platinum PCR SuperMix High Fidelity (Invitrogen) and 0.25 pmol of Ion Torrent primers A and P1. PCR consisted of 72°C for 20 min, 95°C for 5 min followed by 18 cycles of 95°C for 15 s, 62°C for 15 s with a final extension step at 68°C for 1 min. Libraries were then purified twice using a 1.8-fold volume of Solid-phase reversible immobilization (SPRI) bead solution (Meyer & Kircher 2010) to remove fragments smaller than 100 bp.
In the final RAD libraries, the majority of the fragments (65.4% and 67.9% for Sea and Lake pools, respectively) were between 262 and 299 bp (Agilent 2100 bioanalyser). The concentration of the libraries was measured using a Qubit 2.0 Fluorometer (Invitrogen) and then diluted to 28 pM. Sample emulsion PCR, emulsion breaking and enrichment were performed using an Ion Xpress Template Kit, according to the manufacturer's instructions. The libraries (100 lL) were loaded on Ion 314 (10 lL) and 318 (90 lL) chips and sequenced with an Ion Torrent PGM.

Quality control, de novo assembly, mapping and SNP calling
The sequence data acquired using Ion Torrent semiconductor technology were analysed as follows. All the raw sequence reads from the Ion 314 and 318 chips were merged into a single FASTQ file and subsequently split into two (Sea and Lake), based on barcode sequences using fastx_barcode_splitter.pl implemented in FASTX -Toolkit version 0.0.13 (http://hannonlab.cshl.edu/fas-tx_toolkit/index.html). All three files were then trimmed using a FASTQ/A Trimmer with a FASTX -Toolkit. First, one bp was removed from the 3 0 end of all merged raw reads to avoid integration of an extra nucleotide at the end of the restriction enzyme (BamHI) cutting site. Second, after dividing all the reads into Sea and Lake data sets, 10 bp were trimmed from the 5 0 end to remove the barcodes. Subsequently, three filtering criteria were used: (i) reads with lengths of 30-300 bp were retained; (ii) reads were removed if 80% or more of a sequence had a quality score (QS) < 13 and (iii) using a sliding window method with the quality threshold set to 20 and minimum read length to 20 bp. Custom Python scripts (available on Dryad) and computer cluster at the Computing Centre of Finland (CSC; http://www.csc.fi/) were used for all analyses. Following filtering, data quality was checked using PRIN-SEQ-LITE v.0.17.3 (Schmieder & Edwards 2011), after which de novo assembly was performed using MIRA version 3.9.15 (Chevreux et al. 1999) with job specified as genome, de novo and accurate in manifest file using Ion Torrent specific settings. After de novo assembly, the following filtering was carried out to exclude contigs that most likely contain multiple loci. First, using TABLET version 1.13.07.31 (Milne et al. 2013), all the contigs with more than a 2.5% mismatch were discarded. This threshold was set based on the overall distribution of mismatch percentage from the whole data set. Similarly to mismatch percentage, we evaluated the distribution of sequencing depth of all contigs and excluded contigs with more than 5009 coverage from the subsequent analysis to avoid repetitive sequences. To avoid mis-assembled contigs due to chimeric reads, the remaining contigs were split into two or three shorter contigs, depending upon the number of uncut restriction sites found. The contigs consisting of 40 or more nucleotides were selected as reference for mapping in MIRA. For mapping, two files containing population-specific reads were used. The longer reads containing uncut motifs were split from the point of restriction site into shorter (> 20 bp) reads. The job mentioned in the manifest file for MIRA was genome, mapping and accurate, and the chimeric read clipping option was also used. The BAM file produced by MIRA was then used as an input for Fig. 1 The adapters and primers for Ion Torrent sequencing technology. Genomic DNA (gDNA) is simultaneously digested with two restriction enzymes (PstI and BamHI) producing fragments with 5 0 TGCA '3 and 5 0 GATC '3 overhangs, respectively. The 5-bp barcode in forward adapter is marked XXXXX.
SAMTOOLS v.0.1.18 (Li et al. 2009) for SNP discovery. Potential SNPs were accepted only if they met the following stringent selection criteria: read alignment depth of ≤ 5009; ≥ 5 reads per SNP allele; ≤ 2% of polymorphic sites per contig; quality threshold ≥ 20; and ≤ 2% of mismatch. To avoid incorrect SNP calling caused by frequent homopolymer sequencing errors of the Ion Torrent PGM platform, SNPs associated with three or more consecutive homopolymers were discarded.

Genetic differentiation
G ST values for SNPs were calculated using formulae described by Bruneaux et al. (2013). To visualize the effect of coverage on genetic differentiation, a LOESS smoothing model in PAST (Hammer et al. 2001) was used (smoothing value was set at 0.15). As G ST values were strongly affected by coverage, we also used a chisquare test to identify loci showing the strongest genetic differentiation signal between Sea and Lake pools (G€ unther & Coop 2013). For the 25 top SNP loci that showed the strongest genetic differentiation between Sea and Lake pools based on the chi-square test, we also calculated the 95% confidence limits of G ST estimates by resampling with replacement (1000 times) the observed allele counts and recalculating G ST values using POPTOOLS 3.2.5. (Hood 2011). This allowed quantification of the effect of coverage (number of reads) on genetic differentiation estimates; small and/or an uneven number of reads per contig and per population were expected to provide relatively inaccurate G ST estimates, while genetic differentiation estimates are more precise for SNPs, showing high coverage in both population pools.

SNP validation
For SNP validation, primers for PCR amplification were developed for 24 putative SNP markers, which were selected from the pre-reanalysed data set according to their chi-square value (12 markers showing high differentiation between Lake and Sea pools) or the presence of a PstI restriction site to test the potential presence of chimeras (n = 12). For primer development, an online PRIMER3 v.0.4.0 (Koressaar & Remm 2007;Untergasser et al. 2012) program was used and primers were selected only if they met the following selection criteria: any' and 3 0 close to 0; melting temperature 57-60°C and SNP position at least 40 bp from left or right primer to ensure successful sequencing results (primer sequences available on Dryad). Sanger sequencing of seven individuals (four from Lake Peipus and three from P€ arnu Bay) was carried out at Estonian Biocentre (http://vvv.ebc.ee) using 31309l Genetic Analyzer (Applied Biosystems). For sequence alignment and analysis, CLUSTALW multiple alignment option in BIOED-IT v7.2.5 (Hall 1999) was used when possible, otherwise sequences were aligned manually.

Contig annotation and categorization
RAD contigs were annotated using a BLAST v2.2.28 (Altschul et al. 1990) (BLASTXagainst a Uniprot database and BLASTNagainst a nucleotide database) with the E-value threshold set at < 19 e À04 . Gene names and Gene Ontology (GO) terms were retrieved from the Uniprot database according to the BLASTX top hit results. Repetitive genomic elements were detected using REPEATMASKER v4.0.3 software (Smit et al. 1996(Smit et al. -2010.

RAD tag sequencing, de novo assembly and mapping
The sequencing of the RAD libraries with 314 and 318 chips generated 126 327 and 4 045 509 raw reads, respectively. After filtering, 16 886 contigs covering 2 134 779 bp were retained (N50 contig size 82.56 bp; mean coverage 45.199) and used as a reference for mapping. From the combined reads of both pools (1 377 041 for Sea and 1 158 807 for Lake pool), MIRA used 1 259 177 sequences (2 208 684 bp; 711 194 and 547 983 reads from the Sea and Lake pool, respectively) for mapping, with an N50 contig size of 148 bp and a mean coverage of 74.69 (Table S1, Supporting information).

SNP identification and classification
When all available sequence reads (1 259 177) were used against the assembled reference sequence for SNP discovery, 7325 potential SNPs were identified. After stringent filtering, a total of 1259 putative SNPs remained (Table S2, Supporting information), corresponding to a SNP density of one per 1754 bp. The most common SNP variant consisted of G/A or C/T. The distribution of the SNPs along the contigs showed that the frequency of identified polymorphisms decreased after the first 100-150 bp (Fig. 2), indicating that an elevated sequencing error rate towards the end of the reads most likely does not affect the identified 1259 SNPs. The mean overall coverage so-called SNPs was 81.79 (mean coverage of Sea and Lake pools: 46.19 and 35.59, respectively), and 842 of SNPs (67%) had ten or more reads per DNA pool. In total, 32.9% of the filtered SNP-containing contigs had significant BLASTN (23.3%; 293) and/or BLASTX (18.8%;237) hits. Based on the BLASTX results, 168 and 69 SNPs were located in noncoding and coding regions, respectively, the latter corresponding to 40 synonymous and 29 nonsynonymous substitutions.

Genetic differentiation between Sea and Lake pools
Genetic differentiation, measured as G ST , varied considerably (0-0.8933) among 1 228 SNP loci, which contained alleles counted for both pooled samples (Fig. 3). G ST values were strongly affected by coveragemost of the SNPs with very high level of differentiation showed rather low coverage (< 509), while loci with higher coverage showed similar or slightly higher levels of differentiation than individually genotyped microsatellite loci (for 17 microsatellites average G ST is 0.016; 95% CI: 0.007-0.025; 24 individuals analysed per population; (Pukk et al. 2014)). Chi-square test identified loci with high levels of differentiation and low sampling error (e.g. high coverage in both pools) between Sea and Lake pools.

SNP validation
From 24 putative SNP markers selected for validation from pre-reanalysed data set (12 highly differentiated loci and 12 loci containing seemingly uncut PstI restriction site), twelve yielded positive PCR amplification products. Eleven of them represented highly differentiated markers while only a single fragment containing seemingly undigested PstI restriction site was successfully amplified. After cleaning, trimming and aligning the sequences, the presence of both homozygote and heterozygote genotypes could be confirmed in 6 of 12 markers while only homozygotes were observed at four SNP loci and the sequence quality was too low for reliable SNP calling at two loci (Table S2).

Discussion
In this study, we described an efficient genome complexity reduction protocol that uses two rare-cutting restriction enzymes and stringent size selection to make a RAD library compatible for sequencing with Ion Torrent semiconductor technology. We obtained a total of 4 M reads, which were de novo assembled into~17 000 contigs predicted to cover 0.18% of the P. fluviatilis genome (assuming a genome size of 1.19 GB, based on Animal Genome Size Database; http://www.genomesize.com).  Compared to earlier RAD sequencing studies carried out on a wide range of species, our RAD protocol enables very efficient genome complexity reduction ( Fig. 4; Table  S3, Supporting infomation). Therefore, our RAD protocol for the Ion Torrent is expected to be useful for studies which require sequencing of a smaller number of genomic regions, such as mapping QTL loci influencing variation in traits of interest, phylogenetic and population genomic analysis. Using two 6-cutter restriction enzymes (PstI and BamHI) in combination with stringent size selection, we were able to target a small proportion of the genome while retaining high sequencing depth allowing reliable SNP calling despite the higher sequencing errors of the Ion Torrent platform compared to currently available light-based sequencing approaches (Loman et al. 2012;Merriman et al. 2012;Bragg et al. 2013). Unlike earlier RAD sequencing studies, that typically suffer the risk of high DNA loss due to physical shearing and purification steps during agarose gel size selection, our protocol included a fast and simple, yet reproducible, narrow size selection strategy using an E-Gel â SizeSelect TM Agarose Gel. This enabled quick separation and collection of size fractions suitable for a particular sequencing platform, significantly reducing the overall library preparation time. Using our RAD protocol, the mean distance between two contigs in the Eurasian perch is more than 75 Kb, and only a small number of published RAD studies have resulted in a similar level of genome complexity reduction ( Fig. 4; Gagnaire et al. 2013;Hale et al. 2013;Ogden et al. 2013;Richards et al. 2013). All these studies have used a single restriction enzyme (SbfI) with an 8-bp recognition site (5 0 CCTGCAGG 3 0 ) and size selection of DNA fragments from 300 bp to 600-700 bp. Importantly, 8-cutter SbfI creates the identical four base pair overhangs (5 0 TGCA 3 0 ) compared to 6-cutter PstI used in this study and therefore, the same forward adapters can be used in relation to both restriction enzymes to target the optimal number of RAD loci to address a particular research question. As such, our protocol enables fast and flexible RAD sequencing and efficient genome complexity reduction with high coverage using semiconductor sequencing technology.
Eurasian perch, as the most common and widely distributed member of the Percidae family, contributes significantly to commercial and recreational fisheries in Estonia. Because of different fishing regulations, the ability to use genetic markers for assigning individuals back to their source population would be of considerable importance. Also, perch has a great potential for freshwater aquaculture (Rossi et al. 2007), where the use of molecular markers could be beneficial. However, not much is known about the P. fluviatilis genome as the genomic resources of the perch currently include a small number of mitochondrial and nuclear sequences. In this study, we obtained over 2.2 MB of novel sequence data for perch and fifteen per cent of all the consensus sequences gave significant BLASTX hits. As expected, the largest proportion of those hits came from the Nile tilapia O. niloticus (Perciformes). In total, over 1200 putative SNPs were identified and loci with large allele frequency differences between freshwater (Lake Peipus) and brackish water DNA pools (the Baltic Sea) were found, which could be used for the development of genotyping assays to correctly assign fish to populations of origin in the future.
Compared to the microsatellite markers (Pukk et al. 2014), DNA pooling resulted in similar or higher levels of genetic differentiation estimates for the high coverage SNPs (Fig. 3). On the other hand, the SNP markers with relatively low coverage (< 20-509) showed inflated G ST values, demonstrating the importance of high coverage for accurate estimation of genetic differentiation from DNA pools. However, SNP validation of a small subset of loci using Sanger sequencing indicated that not all putative SNP markers represent genuine polymorphisms and therefore, great care should be taken during development of SNP panels for identifying illegal fishing and fish trade. The other critical issue related to allelotyping is the importance of adding equal amount of DNA from each individual to the DNA pool. However, when using large number of individuals for generation of DNA Log10 number of contigs Fig. 4 Varying levels of genome complexity reduction in published RAD studiesgenome size (MB) versus number of contigs (log10 scale). Empty circles correspond to earlier studies from 2010 to 2013 that refer to (Baird et al. (2008) and contained all necessary information (genome size, number of contigs, platform used and fragment size selection). Filled circle corresponds to this study. pools, accurate allele frequency estimation can be achieved even with relatively unequal contributions of each individual forming the pool (Gautier et al. 2013a,b). The results of this study also revealed that majority of the loci with apparently uncut PstI and/or BamHI sites represent most likely chimeric sequences. Therefore, we recommend that future studies either exclude putative chimeric sequences before assembly or alternatively carry out in silico cutting of the fragments before subsequent analysis steps. Further optimization of the library protocols and developing alternative strategies to reduce the occurrence of chimeras would be useful. In summary, this study represents an important step towards developing genomic resources and genetic tools for the Eurasian perch, which could be used in a court of law. Our RAD protocol enables very high genome complexity reduction by targeting a small number of genomic regions while retaining the high coverage necessary for robust SNP discovery and reliable estimation of various population genetic parameters. We expect that the described RAD sequencing protocol for semiconductor sequencing technology will be particularly useful for quick SNP discovery, testing different size selection procedures, restriction enzymes and other critical steps during library preparation. However, Ion Torrent PGM is currently not the cheapest available NGS technology in the market, which combined with rather low throughput and relatively high sequencing error rate makes it less suitable for larger genotyping-by-sequencing projects. Therefore, it is important that both the strengths and limitations of the approach are taken into account when choosing appropriate genome complexity reduction method and sequencing platform (Arnold et al. 2013;Andrews & Luikart 2014;Puritz et al. 2014).

Data Accessibility
Raw sequence data, Sanger sequence alignments, Python scripts and primer sequences are available on Dryad, doi:10.5061/dryad.s2405. Table 1 SuppInfo.xlsb contains the annotated assembled contigs. Table 2 SuppInfo.xlsb contains the filtered SNP data set.

Supporting Information
Additional Supporting Information may be found in the online version of this article:    Baird et al. (2008).