ddRAD-Seq and RAD-Seq for Population Genetics and Phylogenomics: From Library Preparation to Data Analysis

Reduced-representation sequencing solves a fundamental economic problem: whole-genome sequencing of hundreds of individuals is cost-prohibitive for most population genetics studies, yet genotyping only a handful of markers misses the genome-wide signal needed for robust demographic inference and phylogenetic reconstruction. Restriction-site-associated DNA sequencing (RAD-seq) and its derivatives occupy the sweet spot — they sequence a consistent, reproducible subset of the genome adjacent to restriction enzyme cut sites, yielding thousands to tens of thousands of SNP markers distributed across the genome at a fraction of the cost of whole-genome sequencing.

This article provides a practical guide to the ddRAD-seq workflow — from enzyme selection and library preparation through bioinformatic analysis — with emphasis on study design decisions that affect data quality, analytical power, and budget.

Figure 1: RAD-Seq Technology Comparison — RAD, ddRAD, 2b-RAD, and GBS

Reduced-Representation Sequencing Explained — RAD, ddRAD, 2b-RAD, and GBS

The RAD family has diverged into several distinct methods, each with different trade-offs in library preparation complexity, locus repeatability, and marker density. Understanding these differences is essential for selecting the right method — the choice affects not only cost and throughput but also which genomic regions are sampled, how much missing data to expect, and which bioinformatic pipelines are appropriate.

Original RAD-seq, as described by Baird et al. (2008), uses a single restriction enzyme followed by random mechanical shearing to generate fragments for sequencing. While flexible, the shearing step introduces variability in which loci are recovered across libraries — a problem when samples are prepared in different batches. The method also requires multiple purification steps, increasing hands-on time and DNA loss from low-input samples. These limitations motivated the development of double-digest RAD-seq (ddRAD-seq), which replaces mechanical shearing with a second restriction enzyme plus a precise size-selection window. By using one rare-cutting enzyme (e.g., PstI, EcoRI, or SbfI) paired with one common-cutting enzyme (e.g., MspI or MseI), ddRAD-seq defines fragment boundaries enzymatically, substantially improving the consistency of locus recovery across samples. The addition of a narrow size-selection window — typically 300 to 500 base pairs — ensures that only fragments within a predictable size range enter the sequencing library, reducing stochastic variation in locus representation. A 2025 empirical comparison across three enzyme combinations in safflower confirmed that ddRAD-seq consistently outperformed single-digest RAD-seq on raw read count, alignment rate, coverage depth, and SNP yield (Pathania et al., Scientific Reports, 2025). For population genetic studies where F-statistics, ADMIXTURE, and phylogenetic inference depend on genotypes called at shared loci, this improved repeatability makes ddRAD-seq the standard choice for non-model organisms. ddRAD-seq services support the full workflow from enzyme selection through library construction, sequencing, and bioinformatic analysis.

2b-RAD represents a more recent innovation that uses type IIB restriction enzymes (e.g., BsaXI, AlfI) to produce uniform 33 to 36 base-pair fragments from every restriction site in the genome. Because all fragments are the same length, 2b-RAD eliminates the size-selection step entirely, achieving the highest locus repeatability of any RAD method with tag recurrence rates above 95 percent between replicate libraries (Chambers et al., Ecology and Evolution, 2023). The trade-off is read length: at 33 to 36 base pairs, 2b-RAD tags are too short for robust de novo assembly, making the method most suitable when a reference genome is available. For high-precision applications with a reference genome — QTL mapping, GWAS, genomic selection — 2b-RAD sequencing offers a streamlined, highly reproducible alternative.

GBS (genotyping-by-sequencing), while conceptually related, uses a single restriction enzyme and relies on PCR for complexity reduction. It is the most economical option for very large cohorts — breeding panels of thousands of individuals — but produces higher missing-data rates and lower per-locus depth than ddRAD. For ultra-high-throughput projects where per-sample cost is the dominant constraint, genotyping by sequencing services support production-scale genotyping with imputation pipelines to mitigate missing data. For a broader overview of how these methods fit into the larger genotyping landscape — including microsatellites, SNP arrays, and GWAS pipelines — see our Genotyping and Genetic Diversity Services overview.

Figure 2: ddRAD-Seq Wet-Lab Workflow — From DNA Extraction to Sequencing-Ready Library

ddRAD-Seq Wet-Lab Workflow — From DNA to Sequencing-Ready Library

The ddRAD-seq wet-lab workflow consists of five core steps, each with decision points that affect the quality, cost, and informativeness of the final dataset. Getting the wet lab right is more important than optimizing the bioinformatics downstream — no pipeline can rescue libraries built from degraded DNA or digested with poorly chosen enzymes.

DNA quality and quantity

High-molecular-weight DNA is non-negotiable. The restriction digestion step requires intact recognition sites; degraded DNA with nicks and breaks produces fewer fragments that include both enzyme cut sites, reducing the number of loci recovered. A minimum of 500 nanograms of DNA at a concentration of at least 20 nanograms per microliter is recommended, measured by fluorometry (Qubit) rather than spectrophotometry (NanoDrop), since the latter overestimates concentration in the presence of RNA or degraded nucleic acids. For challenging samples — museum specimens, herbarium vouchers, non-invasive hair or fecal samples — DNA integrity should be assessed by gel electrophoresis or a TapeStation before committing to library preparation. Samples with a DNA integrity number (DIN) below 6.0 typically yield unacceptably low library complexity and are better directed toward targeted amplicon or whole-genome approaches that tolerate fragmentation.

Enzyme selection

The choice of restriction enzymes determines which genomic compartments are sampled and how many loci are recovered. GC-rich enzyme recognition sites (e.g., PstI with its CTGCAG motif) tend to enrich for hypomethylated, gene-rich regions in many plant and vertebrate genomes, while AT-rich sites skew toward intergenic regions — a pattern confirmed by Galla-Camps et al. (2024) who analyzed 80 genome assemblies across plants, protostomes, and deuterostomes and found that enzyme recognition site GC content significantly biases locus distribution (BMC Genomics, 2024). For a typical vertebrate genome of 1 to 2 gigabases, EcoRI (rare cutter) plus MspI (common cutter) produces 30,000 to 80,000 loci within a 300 to 500 base-pair window. In silico digestion using tools such as SimRAD or ddgRADer against a reference genome — or a close relative's assembly — should precede wet-lab work to predict fragment counts and size distributions for candidate enzyme pairs. Spending an afternoon on in silico testing routinely prevents weeks of troubleshooting failed libraries.

Size selection precision

This is the step where ddRAD-seq libraries most commonly fail. A loose size-selection window admits fragments that are too short (dominated by adapter dimers at approximately 120 to 130 base pairs) or too long (inefficiently clustered on Illumina flow cells), both of which reduce usable reads per sample. The BluePippin system, which uses pulsed-field electrophoresis for automated, programmable size selection, provides the most reproducible results with a coefficient of variation below 5 percent for the selected size range. Manual gel excision is less expensive but introduces operator-dependent variability that can produce batch effects visible in downstream PCA plots. Regardless of the method, the size-selection window should be aligned with the sequencing read length: for paired-end 150 base-pair sequencing, a 300 to 450 base-pair insert window ensures that paired reads do not overlap excessively and that R2 quality is maintained.

Library QC checkpoints

Three QC measurements should be recorded for every library before pooling: (1) concentration by fluorometry, (2) fragment size distribution by TapeStation or Bioanalyzer — expect a single sharp peak at the selected size window plus adapter length, and (3) the absence of a secondary peak at approximately 120 to 130 base pairs indicating adapter dimers. Libraries with adapter dimer content above 5 percent of total molarity should be re-purified by an additional round of bead cleanup. For studies spanning multiple 96-well plates, include at least one replicate sample across plates to quantify technical variation attributable to library batch — a practice that costs one extra library but can distinguish biological from technical signal when unexpected structure appears in downstream analyses.

Pooling and sequencing

Equimolar pooling across samples is critical for even read distribution. A commonly cited target is 1 to 5 million reads per sample, which provides sufficient depth for genotype calling when 30,000 to 80,000 loci are targeted. Paired-end 150 base-pair sequencing on Illumina platforms is standard; single-end reads can suffice for reference-guided analysis but limit de novo assembly options. Recent protocol improvements include the use of quick-acting ligases compatible with restriction enzyme buffers — eliminating buffer-exchange steps — combined with single-step amplification-and-barcoding PCR, which reduce hands-on time while improving library yield across diverse sample types.

Figure 3: Bioinformatics Pipeline — Stacks vs. ipyrad, De Novo vs. Reference-Guided Analysis

Bioinformatics Pipeline — Stacks, ipyrad, and the De Novo vs. Reference-Guided Decision

The raw output of a ddRAD-seq run is a set of demultiplexed FASTQ files. Transforming these into a filtered, analysis-ready SNP matrix requires a bioinformatic pipeline that handles quality filtering, locus assembly or alignment, variant calling, and population-aware filtering. Two pipelines dominate the ddRAD literature: Stacks (version 2) and ipyrad (version 0.9.x).

Stacks 2 uses a modular architecture: process_radtags demultiplexes and filters reads; ustacks assembles loci within each sample at a user-defined mismatch threshold (M, default 2–4); cstacks builds a catalog of consensus loci; sstacks matches samples to the catalog; gstacks calls variants; and populations exports filtered genotype matrices. Stacks supports both de novo (denovo_map.pl) and reference-guided (ref_map.pl) pipelines, and its paired-end-aware assembly — which concatenates forward and reverse reads from the same fragment into longer contigs — is a distinguishing advantage for ddRAD data. The parameters that most affect results are M (within-sample mismatches), n (between-sample mismatches during catalog building), and m (minimum read depth to form a stack, default 3). Raising M and n to 5–8 can increase SNP yield without inflating missing data, but optimal values depend on the study species' diversity and should be tested empirically on a subset of samples.

ipyrad takes a seven-step approach: demultiplexing, read filtering, within-sample clustering at a tunable similarity threshold (clust_threshold, default 0.85), joint estimation of heterozygosity and error, consensus calling, cross-sample clustering, and output formatting. ipyrad's branching architecture allows parallel exploration of multiple parameter sets without recomputing shared steps — useful for sensitivity testing. The critical parameters are clust_threshold (higher values increase locus recovery but potentially inflate paralog inclusion), min_samples_locus, and mindepth_majrule.

The de novo vs. reference-guided decision. When a reference genome is available, reference-guided pipelines offer higher sensitivity and produce coordinates comparable across studies. However, a comparative study running the same ddRAD data through both Stacks reference-based and ipyrad reference assembly illustrates the practical reality: when the same ddRAD data are run through both pipelines, the overlap of called SNPs can be surprisingly small — studies across multiple taxa have reported that only a fraction of variants are shared between pipelines, with Stacks and ipyrad each recovering thousands of pipeline-specific SNPs in addition to a core set of cross-validated markers. This finding, replicated across study systems, argues strongly for cross-validating key biological conclusions against pipeline choice. For researchers establishing ddRAD analysis workflows, bioinformatics services can provide both Stacks and ipyrad pipelines with parameter sensitivity testing and cross-validation reporting.

Filtering and missing data management. Standard post-calling filters include: minor allele frequency at 0.01 to 0.05, per-locus call rate at 70 to 80 percent, and Hardy-Weinberg equilibrium deviation (p > 0.001) to flag genotyping errors. Paralog filtering deserves particular attention: because ddRAD-seq sequences short fragments around restriction sites, paralogous regions with conserved sites can be mistakenly assembled as a single locus, producing inflated heterozygosity and spurious admixture signals. The HDplot method, which combines per-SNP heterozygosity with allele read-ratio deviation using the AD (allele depth) field from Stacks VCF output, can identify and remove paralogous loci. LD pruning before analyses that assume marker independence — ADMIXTURE, PCA — should use sliding windows with an r-squared threshold calibrated by minor allele frequency bin to avoid preferentially removing rare variants.

Figure 4: Population Genetics Applications — Diversity, Structure, and Demography

Population Genetics Applications — Diversity, Structure, and Demography from RAD Data

A filtered ddRAD-seq SNP matrix supports a comprehensive suite of population genetic analyses. The analytical toolkit is mature; the principal challenge is interpreting results in light of reduced-representation data limitations — particularly missing data patterns, ascertainment bias from enzyme choice, and the sensitivity of biological conclusions to parameter choices.

Genetic diversity and differentiation. Standard diversity indices — observed heterozygosity (Ho), expected heterozygosity (He), nucleotide diversity (pi), and the inbreeding coefficient (Fis) — are computed per population from VCF output. Because ddRAD targets a non-random subset of the genome, absolute pi values should be interpreted as relative comparisons across populations within a study rather than as unbiased genome-wide estimates. Population differentiation is measured by pairwise Fst (Weir and Cockerham's estimator), which for ddRAD data is robust to moderate missing data when per-population sample sizes exceed 6 to 8 individuals. Analysis of molecular variance (AMOVA) partitions genetic variance among hierarchically defined population groups. Mantel tests or distance-based redundancy analysis (dbRDA) test for isolation-by-distance patterns against geographic or environmental distance matrices. Population evolution and genomics services support the full pipeline from variant calling through diversity estimation, structure analysis, and selection scan detection.

Population structure. ADMIXTURE and sNMF estimate individual ancestry proportions for a specified number of ancestral populations (K), with cross-validation error guiding the choice of K. PCA and DAPC provide complementary, model-free visualizations. With 10,000 to 30,000 ddRAD SNPs, both methods reliably detect structure at Fst values as low as 0.01 to 0.02 when sample sizes are adequate and missing data are below 20 percent. For fine-scale structure, fineRADstructure uses co-ancestry matrices derived from haplotype-level information in RAD loci to resolve relationships that allele-frequency-based methods miss.

Demographic history. Stairway Plot 2 infers historical effective population size trajectories from the site frequency spectrum. daDi and moments fit explicit demographic models — isolation-with-migration, secondary contact, expansion — to the joint site frequency spectrum, testing alternative hypotheses via likelihood ratio tests. A caution specific to RAD data: demographic methods assume unlinked, neutrally evolving SNPs, but ddRAD loci are physically clustered around restriction sites. Thinning to one SNP per locus and filtering high-LD regions is essential before site-frequency-spectrum-based inference.

Real inquiry patterns. Researchers contacting CD Genomics with ddRAD-seq projects typically fall into three categories: (1) population structure and gene flow — for example, 100 beetle individuals from 10 localities requiring 10,000 to 30,000 SNPs to resolve fine-scale differentiation; (2) phylogeographic reconstruction — 60 Calochortus individuals across 6 taxa, using coalescent species tree methods to resolve relationships that organellar markers could not; and (3) conservation and invasion genetics — 40 invasive species samples spanning the introduced range, identifying introduction sources and dispersal corridors. Each scenario presents different demands on DNA quality, marker density, and analytical method, underscoring the value of method-to-question matching discussed in our genotyping services overview.

Figure 5: Phylogenomics with RAD Data — Concatenation vs. Coalescent Methods

Phylogenomics with RAD Data — Concatenation, Coalescence, and Species Tree Inference

ddRAD-seq has become a primary data source for phylogenomic studies in non-model clades, displacing Sanger-sequenced markers that lacked sufficient information to resolve recently diverged lineages or nodes affected by incomplete lineage sorting (ILS). However, RAD data present distinct challenges for phylogenetics that do not arise with whole-genome or targeted-enrichment datasets.

Concatenation vs. coalescent methods. Concatenation — assembling all RAD loci into a supermatrix and analyzing it with maximum likelihood (IQ-TREE with ModelFinder and ultrafast bootstrap) — is computationally straightforward and well powered when gene tree discordance is low. However, concatenation ignores the reality that individual RAD loci have independent genealogical histories. When ILS is high — common in rapid radiations — concatenation can produce highly supported but incorrect trees. Coalescent-based methods such as ASTRAL and SVDquartets estimate the species tree while accounting for gene tree heterogeneity. Studies in amphibian systems using ddRAD data have found that concatenated and coalescent species trees are typically concordant at deeper nodes but may diverge at shallow nodes where incomplete lineage sorting is prevalent, illustrating why both methods should be reported.

Locus filtering for phylogenetics. Standard population genetic filters — MAF, HWE — are often inappropriate for phylogenetic datasets spanning species boundaries. Phylogenetic filtering instead focuses on: locus completeness (data for at least 50 to 70 percent of taxa), removal of loci with anomalously high heterozygosity (potential paralogs), and identification of loci with outlier phylogenetic signal using tools such as PhyParts and Quartet Sampling that map gene tree conflict across the species tree.

Hybridization detection. RAD data's high marker density enables detection of hybridization through D-statistics (ABBA-BABA tests), TreeMix (population splits with migration edges), and HyDe (phylogenetic invariants for hybrid speciation). Because ILS and genuine gene flow can produce similar allele-sharing patterns, running multiple detection methods and reporting concordant results strengthens the case for genuine hybridization. A practical design recommendation: for studies where hybridization is a central question, increase per-population sampling to 10 to 15 individuals — the additional individuals improve allele frequency estimates that power D-statistics and TreeMix far more than adding loci to smaller population samples.

Moving from discovery to validation. When ddRAD data identify SNPs associated with species boundaries or adaptive divergence, extending findings to additional individuals using targeted genotyping is often cost-effective. Microsatellite genotyping and targeted SNP approaches for validation are discussed in our guide on marker-based population analysis.

Cost, Sample Requirements, and Experimental Design

Budget, not technology, is the binding constraint in most ddRAD-seq studies. Understanding the cost structure — and where corners can be cut versus where they cannot — enables informed decisions before funds are committed.

DNA quality — the irreducible minimum

As discussed above, degraded DNA is the single most common cause of ddRAD-seq project failure. Beyond the minimum input requirements covered in the wet-lab section, one additional consideration applies at the budgeting stage: for field-collected samples stored in ethanol, extraction methods that include an RNase treatment and a final ethanol precipitation step typically outperform column-based kits for DNA integrity and yield. For museum or herbarium specimens with severely degraded DNA, ddRAD-seq is unlikely to succeed without substantial optimization; whole-genome sequencing with gentler library preparation or targeted amplicon approaches may be more appropriate.

SNP yield expectations

The number of SNPs recovered depends on nucleotide diversity, enzyme pair, size-selection window, and sequencing depth. For most diploid eukaryotes with moderate diversity, a well-designed ddRAD library yields 10,000 to 50,000 SNPs after filtering. Highly inbred species or those that experienced recent bottlenecks may yield fewer than 5,000; species with exceptionally high diversity can yield over 100,000. A pilot study of 8 to 16 individuals is strongly recommended before scaling — it provides empirical SNP yield data that informs sequencing depth decisions and validates that the chosen enzyme pair produces adequate loci. Spending 10 to 15 percent of the total budget on a pilot is the most cost-effective investment a ddRAD-seq project can make: it answers whether DNA quality supports library preparation, whether the enzyme pair delivers the expected locus count, what the realized SNP yield and missing-data rate are, and whether there are unexpected sources of technical variation requiring protocol adjustment before production.

Cost structure at different scales

ddRAD-seq costs divide into library preparation and sequencing. At small scale (fewer than 50 samples), per-sample library costs dominate; invest in the best data for the question. At medium scale (50 to 200 samples), batch processing brings per-sample costs into the $15 to $25 range for library preparation plus $10 to $20 for sequencing at 2 million reads per sample. At large scale (200 to 500 samples), dual-indexing schemes and bulk reagent purchasing achieve total per-sample costs of $25 to $40. For studies exceeding 500 samples, genotyping by sequencing (GBS) with its simpler single-enzyme protocol offers lower per-sample costs at the expense of higher missing-data rates, though imputation can partially compensate. For researchers developing ddRAD-seq projects, ddRAD-seq services include pilot-scale library preparation and sequencing with full QC reporting, allowing methods to be validated before committing to the full cohort.

Connecting discovery to phenotype. ddRAD-seq excels as a discovery tool for population structure, diversity, and phylogenomic inference. When studies progress from discovery to trait mapping, the same SNP resources can be leveraged for association analysis. Our guide on GWAS experimental design with GBS discusses how reduced-representation SNP data connect genotype to phenotype, completing the arc from marker discovery through functional analysis.

FAQ

What is the difference between RAD-seq, ddRAD-seq, and 2b-RAD?

Original RAD-seq uses one restriction enzyme plus random shearing, introducing variability in locus recovery. ddRAD-seq uses two enzymes plus precise size selection, improving cross-sample repeatability. 2b-RAD uses type IIB enzymes to produce uniform 33–36 bp fragments from all restriction sites, achieving the highest repeatability but requiring a reference genome for short-read mapping.

How much DNA do I need for ddRAD-seq?

A minimum of 500 nanograms at 20 nanograms per microliter, measured by fluorometry. More is better — 1 to 2 micrograms provides margin for repeated steps. Degraded DNA with a DIN below 6.0 typically yields unacceptably low library complexity.

Do I need a reference genome for ddRAD-seq?

No. ddRAD-seq works well de novo using clustering-based pipelines such as Stacks denovo_map.pl or ipyrad. The 400 to 500 base-pair contigs assembled from paired-end reads are long enough for de novo locus identification and flanking primer design. Having a reference genome improves variant calling sensitivity and annotation but is not required, which is why ddRAD dominates non-model organism research.

How many SNPs can I expect from ddRAD-seq?

For most diploid eukaryotes with moderate diversity, expect 10,000 to 50,000 SNPs after standard filtering. Highly inbred species may yield fewer than 5,000; species with exceptionally high diversity can yield over 100,000. A pilot study with 8 to 16 individuals is the most reliable way to obtain empirical estimates for your species and enzyme pair.

Which pipeline should I use — Stacks or ipyrad?

There is no universally superior pipeline. Stacks excels for paired-end ddRAD with its native contig assembly and provides the AD field needed for paralog filtering via HDplot. ipyrad's branching architecture enables efficient parameter sensitivity testing on HPC clusters. The most robust approach is to run both on a subset of samples and cross-validate — SNPs detected by both pipelines are high-confidence.

What is the advantage of ddRAD-seq over whole-genome resequencing?

Cost and reference independence. ddRAD-seq delivers 10,000 to 50,000 genome-wide SNPs at $25 to $40 per sample including library and sequencing, making it feasible for the 100 to 500 individual sample sizes needed for robust population genetic inference. WGS provides more markers and detects structural variants but at several times the per-sample cost, and de novo assembly of WGS data without a reference genome remains computationally intensive.

References:

1. Pathania A, Sharma R, et al. Comparative analysis of RAD-seq methods for SNP discovery and genetic diversity assessment in oil seed crop safflower. Scientific Reports. 2025;15:22600. https://doi.org/10.1038/s41598-025-06706-2
2. Chambers EA, Tarvin RD, Santos JC, Ron SR, Betancourth-Cundar M, Hillis DM, Matz MV, Cannatella DC. 2b or not 2b? 2bRAD is an effective alternative to ddRAD for phylogenomics. Ecology and Evolution. 2023;13(3):e9842. https://doi.org/10.1002/ece3.9842
3. Galla-Camps M, Carreras C, Pascual M, Pegueroles C. Genome composition and GC content influence loci distribution in reduced representation genomic studies. BMC Genomics. 2024;25:410. https://doi.org/10.1186/s12864-024-10312-3
4. Scariolo F, Draga S, et al. A pioneering genotypic and phylogenetic characterisation of Cichorium crops through a genome-scale sequencing for future breeding innovations. BMC Plant Biology. 2025;25:860. https://doi.org/10.1186/s12870-025-06876-1
5. Perez-Bello P, Panero I, De Paoli E, Casolo V, Attorre F, Cambria VE, Strumia S, Santangelo A, Bonomi C, Fabrini G, Marroni F. ddRAD sequencing of the endangered species Primula palinuri Petagna reveals high levels of inter-population diversity. Scientific Reports. 2025;15:15245. https://doi.org/10.1038/s41598-025-98334-z
6. Toker TP, Ulusoy D, Dogan B, Kasapoglu S, Hakan F, Reddy UK, Kordrostami M, Yol E. Genomic insights into Mediterranean pepper diversity using ddRADSeq. PLOS ONE. 2025;20(3):e0318105. https://doi.org/10.1371/journal.pone.0318105
7. Baird NA, Etter PD, Atwood TS, Currey MC, Shiver AL, Lewis ZA, Selker EU, Cresko WA, Johnson EA. Rapid SNP discovery and genetic mapping using sequenced RAD markers. PLOS ONE. 2008;3(10):e3376. https://doi.org/10.1371/journal.pone.0003376

For research use only, not intended for clinical diagnosis, treatment, or individual health assessments.