Sequence capture using PCR‐generated probes: a cost‐effective method of targeted high‐throughput sequencing for nonmodel organisms

Recent advances in high‐throughput sequencing library preparation and subgenomic enrichment methods have opened new avenues for population genetics and phylogenetics of nonmodel organisms. To multiplex large numbers of indexed samples while sequencing predominantly orthologous, targeted regions of the genome, we propose modifications to an existing, in‐solution capture that utilizes PCR products as target probes to enrich library pools for the genomic subset of interest. The sequence capture using PCR‐generated probes (SCPP) protocol requires no specialized equipment, is highly flexible and significantly reduces experimental costs for projects where a modest scale of genetic data is optimal (25–100 genomic loci). Our alterations enable application of this method across a wider phylogenetic range of taxa and result in higher capture efficiencies and coverage at each locus. Efficient and consistent capture over multiple SCPP experiments and at various phylogenetic distances is demonstrated, extending the utility of this method to both phylogeographic and phylogenomic studies.     

Long‐term sperm storage in the brownbanded bamboo shark Chiloscyllium punctatum

This study investigated the birth of a brownbanded bamboo shark Chiloscyllium punctatum at the Steinhart Aquarium. Genetic analyses suggest this is the longest documented case of sperm storage for any species of shark (45 months).     

Estimating and accounting for genotyping errors in RAD‐seq experiments

In non‐model organisms, evolutionary questions are frequently addressed using reduced representation sequencing techniques due to their low cost, ease of use, and because they do not require genomic resources such as a reference genome. However, evidence is accumulating that such techniques may be affected by specific biases, questioning the accuracy of obtained genotypes, and as a consequence, their usefulness in evolutionary studies. Here, we introduce three strategies to estimate genotyping error rates from such data: through the comparison to high quality genotypes obtained with a different technique, from individual replicates, or from a population sample when assuming Hardy‐Weinberg equilibrium. Applying these strategies to data obtained with Restriction site Associated DNA sequencing (RAD‐seq), arguably the most popular reduced representation sequencing technique, revealed per‐allele genotyping error rates that were much higher than sequencing error rates, particularly at heterozygous sites that were wrongly inferred as homozygous. As we exemplify through the inference of genome‐wide and local ancestry of well characterized hybrids of two Eurasian poplar (Populus) species, such high error rates may lead to wrong biological conclusions. By properly accounting for these error rates in downstream analyses, either by incorporating genotyping errors directly or by recalibrating genotype likelihoods, we were nevertheless able to use the RAD‐seq data to support biologically meaningful and robust inferences of ancestry among Populus hybrids. Based on these findings, we strongly recommend carefully assessing genotyping error rates in reduced representation sequencing experiments, and to properly account for these in downstream analyses, for instance using the tools presented here.     

Applicability of RAD‐tag genotyping for interfamilial comparisons: empirical data from two cetaceans

Restriction‐site‐associated DNA tag (RAD‐tag) sequencing has become a popular approach to generate thousands of SNPs used to address diverse questions in population genomics. Comparatively, the suitability of RAD‐tag genotyping to address evolutionary questions across divergent species has been the subject of only a few recent studies. Here, we evaluate the applicability of this approach to conduct genome‐wide scans for polymorphisms across two cetacean species belonging to distinct families: the short‐beaked common dolphin (Delphinus delphis; n = 5 individuals) and the harbour porpoise (Phocoena phocoena; n = 1 individual). Additionally, we explore the effects of varying two parameters in the Stacks analysis pipeline on the number of loci and level of divergence obtained. We observed a 34% drop in the total number of loci that were present in all individuals when analysing individuals from the distinct families compared with analyses restricted to intraspecific comparisons (i.e. within D. delphis). Despite relatively stringent quality filters, 3595 polymorphic loci were retrieved from our interfamilial comparison. Cetaceans have undergone rapid diversification, and the estimated divergence time between the two families is relatively recent (14–19 Ma). Thus, our results showed that, for this level of divergence, a large number of orthologous loci can still be genotyped using this approach, which is on par with two recent in silico studies. Our findings constitute one of the first empirical investigations using RAD‐tag sequencing at this level of divergence and highlights the great potential of this approach in comparative studies and to address evolutionary questions.     

Multiplex restriction amplicon sequencing: a novel next‐generation sequencing‐based marker platform for high‐throughput genotyping

To enable rapid selection of traits in marker‐assisted breeding, markers must be technically simple, low‐cost, high‐throughput and randomly distributed in a genome. We developed such a technology, designated as Multiplex Restriction Amplicon Sequencing (MRASeq), which reduces genome complexity by polymerase chain reaction (PCR) amplification of amplicons flanked by restriction sites. The first PCR primers contain restriction site sequences at 3’‐ends, preceded by 6‐10 bases of specific or degenerate nucleotide sequences and then by a unique M13‐tail sequence which serves as a binding site for a second PCR that adds sequencing primers and barcodes to allow sample multiplexing for sequencing. The sequences of restriction sites and adjacent nucleotides can be altered to suit different species. Physical mapping of MRASeq SNPs from a biparental population of allohexaploid wheat (Triticum aestivum L.) showed a random distribution of SNPs across the genome. MRASeq generated thousands of SNPs from a wheat biparental population and natural populations of wheat and barley (Hordeum vulgare L.). This novel, next‐generation sequencing‐based genotyping platform can be used for linkage mapping to screen quantitative trait loci (QTL), background selection in breeding and many other genetics and breeding applications of various species.     

Large‐scale connectivity,cryptic population structure,and relatedness in Eastern Pacific Olive ridley sea turtles (Lepidochelys olivacea)

Endangered species are grouped into genetically discrete populations to direct conservation efforts. Mitochondrial control region (mtCR) haplotypes are used to elucidate deep divergences between populations, as compared to nuclear microsatellites that can detect recent structuring. When prior populations are unknown, it is useful to subject microsatellite data to clustering and/or ordination population inference. Olive ridley sea turtles (Lepidochelys olivacea) are the most abundant sea turtle, yet few studies have characterized olive ridley population structure. Recently, clustering results of olive ridleys in the Eastern Tropical Pacific Ocean suggested weak structuring (F_ST = 0.02) between Mexico and Central America. We analyzed mtCR haplotypes, new microsatellite genotypes from Costa Rica, and preexisting microsatellite genotypes from olive ridleys across the Eastern Tropical Pacific, to further explore population structuring in this region. We subjected inferred populations to multiple analyses to explore the mechanisms behind their structuring. We found 10 mtCR haplotypes from 60 turtles nesting at three sites in Costa Rica, but did not detect divergence between Costa Rican sites, or between Central America and Mexico. In Costa Rica, clustering suggested one population with no structuring, but ordination suggested four cryptic clusters with moderate structuring (F_ST = 0.08, p < .001). Across the Eastern Tropical Pacific, ordination suggested nine cryptic clusters with moderate structuring (F_ST = 0.103, p < .001) that largely corresponded to Mexican and Central American populations. All ordination clusters displayed significant internal relatedness relative to global relatedness (p < .001) and contained numerous sibling pairs. This suggests that broadly dispersed family lineages have proliferated in Eastern Tropical Pacific olive ridleys and corroborates previous work showing basin‐wide connectivity and shallow population structure in this region. The existence of broadly dispersed kin in Eastern Tropical Pacific olive ridleys has implications for management of olive ridleys in this region, and adds to our understanding of sea turtle ecology and life history, particularly in light of the natal‐homing paradigm.     

RAD‐sequencing for estimating genomic relatedness matrix‐based heritability in the wild: A case study in roe deer

Estimating the evolutionary potential of quantitative traits and reliably predicting responses to selection in wild populations are important challenges in evolutionary biology. The genomic revolution has opened up opportunities for measuring relatedness among individuals with precision, enabling pedigree‐free estimation of trait heritabilities in wild populations. However, until now, most quantitative genetic studies based on a genomic relatedness matrix (GRM) have focused on long‐term monitored populations for which traditional pedigrees were also available, and have often had access to knowledge of genome sequence and variability. Here, we investigated the potential of RAD‐sequencing for estimating heritability in a free‐ranging roe deer (Capreolous capreolus) population for which no prior genomic resources were available. We propose a step‐by‐step analytical framework to optimize the quality and quantity of the genomic data and explore the impact of the single nucleotide polymorphism (SNP) calling and filtering processes on the GRM structure and GRM‐based heritability estimates. As expected, our results show that sequence coverage strongly affects the number of recovered loci, the genotyping error rate and the amount of missing data. Ultimately, this had little effect on heritability estimates and their standard errors, provided that the GRM was built from a minimum number of loci (above 7,000). Genomic relatedness matrix‐based heritability estimates thus appear robust to a moderate level of genotyping errors in the SNP data set. We also showed that quality filters, such as the removal of low‐frequency variants, affect the relatedness structure of the GRM, generating lower h² estimates. Our work illustrates the huge potential of RAD‐sequencing for estimating GRM‐based heritability in virtually any natural population.     

Mapping migration in a songbird using high‐resolution genetic markers

Neotropic migratory birds are declining across the Western Hemisphere, but conservation efforts have been hampered by the inability to assess where migrants are most limited—the breeding grounds, migratory stopover sites or wintering areas. A major challenge has been the lack of an efficient, reliable and broadly applicable method for measuring the strength of migratory connections between populations across the annual cycle. Here, we show how high‐resolution genetic markers can be used to identify genetically distinct groups of a migratory bird, the Wilson's warbler (Cardellina pusilla), at fine enough spatial scales to facilitate assessing regional drivers of demographic trends. By screening 1626 samples using 96 highly divergent single nucleotide polymorphisms selected from a large pool of candidates (~450 000), we identify novel region‐specific migratory routes and timetables of migration along the Pacific Flyway. Our results illustrate that high‐resolution genetic markers are more reliable, precise and amenable to high throughput screening than previously described intrinsic marking techniques, making them broadly applicable to large‐scale monitoring and conservation of migratory organisms.     

Targeted multiplex next‐generation sequencing: advances in techniques of mitochondrial and nuclear DNA sequencing for population genomics

Next‐generation sequencing (NGS) is emerging as an efficient and cost‐effective tool in population genomic analyses of nonmodel organisms, allowing simultaneous resequencing of many regions of multi‐genomic DNA from multiplexed samples. Here, we detail our synthesis of protocols for targeted resequencing of mitochondrial and nuclear loci by generating indexed genomic libraries for multiplexing up to 100 individuals in a single sequencing pool, and then enriching the pooled library using custom DNA capture arrays. Our use of DNA sequence from one species to capture and enrich the sequencing libraries of another species (i.e. cross‐species DNA capture) indicates that efficient enrichment occurs when sequences are up to about 12% divergent, allowing us to take advantage of genomic information in one species to sequence orthologous regions in related species. In addition to a complete mitochondrial genome on each array, we have included between 43 and 118 nuclear loci for low‐coverage sequencing of between 18 kb and 87 kb of DNA sequence per individual for single nucleotide polymorphisms discovery from 50 to 100 individuals in a single sequencing lane. Using this method, we have generated a total of over 500 whole mitochondrial genomes from seven cetacean species and green sea turtles. The greater variation detected in mitogenomes relative to short mtDNA sequences is helping to resolve genetic structure ranging from geographic to species‐level differences. These NGS and analysis techniques have allowed for simultaneous population genomic studies of mtDNA and nDNA with greater genomic coverage and phylogeographic resolution than has previously been possible in marine mammals and turtles.     

Population genomic analyses from low‐coverage RAD‐Seq data: a case study on the non‐model cucurbit bottle gourd

Restriction site‐associated DNA sequencing (RAD‐Seq), a next‐generation sequencing‐based genome ‘complexity reduction’ protocol, has been useful in population genomics in species with a reference genome. However, the application of this protocol to natural populations of genomically underinvestigated species, particularly under low‐to‐medium sequencing depth, has not been well justified. In this study, a Bayesian method was developed for calling genotypes from an F₂ population of bottle gourd [Lagenaria siceraria (Mol.) Standl.] to construct a high‐density genetic map. Low‐depth genome shotgun sequencing allowed the assembly of scaffolds/contigs comprising approximately 50% of the estimated genome, of which 922 were anchored for identifying syntenic regions between species. RAD‐Seq genotyping of a natural population comprising 80 accessions identified 3226 single nuclear polymorphisms (SNPs), based on which two sub‐gene pools were suggested for association with fruit shape. The two sub‐gene pools were moderately differentiated, as reflected by the Hudson's F_ST value of 0.14, and they represent regions on LG7 with strikingly elevated F_ST values. Seven‐fold reduction in heterozygosity and two times increase in LD (r²) were observed in the same region for the round‐fruited sub‐gene pool. Outlier test suggested the locus LX3405 on LG7 to be a candidate site under selection. Comparative genomic analysis revealed that the cucumber genome region syntenic to the high F_ST island on LG7 harbors an ortholog of the tomato fruit shape gene OVATE. Our results point to a bright future of applying RAD‐Seq to population genomic studies for non‐model species even under low‐to‐medium sequencing efforts. The genomic resources provide valuable information for cucurbit genome research.     

Helping decision making for reliable and cost‐effective 2b‐RAD sequencing and genotyping analyses in non‐model species

High‐throughput sequencing has revolutionized population and conservation genetics. RAD sequencing methods, such as 2b‐RAD, can be used on species lacking a reference genome. However, transferring protocols across taxa can potentially lead to poor results. We tested two different IIB enzymes (AlfI and CspCI) on two species with different genome sizes (the loggerhead turtle Caretta caretta and the sharpsnout seabream Diplodus puntazzo) to build a set of guidelines to improve 2b‐RAD protocols on non‐model organisms while optimising costs. Good results were obtained even with degraded samples, showing the value of 2b‐RAD in studies with poor DNA quality. However, library quality was found to be a critical parameter on the number of reads and loci obtained for genotyping. Resampling analyses with different number of reads per individual showed a trade‐off between number of loci and number of reads per sample. The resulting accumulation curves can be used as a tool to calculate the number of sequences per individual needed to reach a mean depth ≥20 reads to acquire good genotyping results. Finally, we demonstrated that selective‐base ligation does not affect genomic differentiation between individuals, indicating that this technique can be used in species with large genome sizes to adjust the number of loci to the study scope, to reduce sequencing costs and to maintain suitable sequencing depth for a reliable genotyping without compromising the results. Here, we provide a set of guidelines to improve 2b‐RAD protocols on non‐model organisms with different genome sizes, helping decision‐making for a reliable and cost‐effective genotyping.     

Development of a universal double‐digest RAD sequencing approach for a group of nonmodel,ecologically and economically important insect and fish taxa

The generation of genome‐scale data is critical for a wide range of questions in basic biology using model organisms, but also in questions of applied biology in nonmodel organisms (agriculture, natural resources, conservation and public health biology). Using a genome‐scale approach on a diverse group of nonmodel organisms and with the goal of lowering costs of the method, we modified a multiplexed, high‐throughput genomic scan technique utilizing two restriction enzymes. We analysed several pairs of restriction enzymes and completed double‐digestion RAD sequencing libraries for nine different species and five genera of insects and fish. We found one particular enzyme pair produced consistently higher number of sequence‐able fragments across all nine species. Building libraries off this enzyme pair, we found a range of usable SNPs between 4000 and 37 000 SNPS per species and we found a greater number of usable SNPs using reference genomes than de novo pipelines in STACKS. We also found fewer reads in the Read 2 fragments from the paired‐end Illumina Hiseq run. Overall, the results of this study provide empirical evidence of the utility of this method for producing consistent data for diverse nonmodel species and suggest specific considerations for sequencing analysis strategies.     

SNP discovery using Paired‐End RAD‐tag sequencing on pooled genomic DNA of Sisymbrium austriacum (Brassicaceae)

Single nucleotide polymorphisms SNPs are rapidly replacing anonymous markers in population genomic studies, but their use in non model organisms is hampered by the scarcity of cost‐effective approaches to uncover genome‐wide variation in a comprehensive subset of individuals. The screening of one or only a few individuals induces ascertainment bias. To discover SNPs for a population genomic study of the Pyrenean rocket (Sisymbrium austriacum subsp. chrysanthum), we undertook a pooled RAD‐PE (Restriction site Associated DNA Paired‐End sequencing) approach. RAD tags were generated from the PstI‐digested pooled genomic DNA of 12 individuals sampled across the species distribution range and paired‐end sequenced using Illumina technology to produce ~24.5 Mb of sequences, covering ~7% of the specie's genome. Sequences were assembled into ~76 000 contigs with a mean length of 323 bp (N₅₀ = 357 bp, sequencing depth = 24x). In all, >15 000 SNPs were called, of which 47% were annotated in putative genic regions based on homology with the Arabidopsis thaliana genome. Gene ontology (GO) slim categorization demonstrated that the identified SNPs covered extant genic variation well. The validation of 300 SNPs on a larger set of individuals using a KASPar assay underpinned the utility of pooled RAD‐PE as an inexpensive genome‐wide SNP discovery technique (success rate: 87%). In addition to SNPs, we discovered >600 putative SSR markers.     

Combining genome‐wide association study and FST‐based approaches to identify targets of Borrelia‐mediated selection in natural rodent hosts

Recent advances in high‐throughput sequencing technologies provide opportunities to gain novel insights into the genetic basis of phenotypic trait variation. Yet to date, progress in our understanding of genotype–phenotype associations in nonmodel organisms in general and natural vertebrate populations in particular has been hampered by small sample sizes typically available for wildlife populations and a resulting lack of statistical power, as well as a limited ability to control for false‐positive signals. Here we propose to combine a genome‐wide association study (GWAS) and F_ST‐based approach with population‐level replication to partly overcome these limitations. We present a case study in which we used this approach in combination with genotyping‐by‐sequencing (GBS) single nucleotide polymorphism (SNP) data to identify genomic regions associated with Borrelia afzelii resistance or susceptibility in the natural rodent host of this Lyme disease‐causing spirochete, the bank vole (Myodes glareolus). Using this combined approach we identified four consensus SNPs located in exonic regions of the genes Slc26a4, Tns3, Wscd1 and Espnl, which were significantly associated with the voles’ Borrelia infectious status within and across populations. Functional links between host responses to bacterial infections and most of these genes have previously been demonstrated in other rodent systems, making them promising new candidates for the study of evolutionary host responses to Borrelia emergence. Our approach is applicable to other systems and may facilitate the identification of genetic variants underlying disease resistance or susceptibility, as well as other ecologically relevant traits, in wildlife populations.     

Double‐digest RAD sequencing using Ion Proton semiconductor platform (ddRADseq‐ion) with nonmodel organisms

Research in evolutionary biology involving nonmodel organisms is rapidly shifting from using traditional molecular markers such as mtDNA and microsatellites to higher throughput SNP genotyping methodologies to address questions in population genetics, phylogenetics and genetic mapping. Restriction site associated DNA sequencing (RAD sequencing or RADseq) has become an established method for SNP genotyping on Illumina sequencing platforms. Here, we developed a protocol and adapters for double‐digest RAD sequencing for Ion Torrent (Life Technologies; Ion Proton, Ion PGM) semiconductor sequencing. We sequenced thirteen genomic libraries of three different nonmodel vertebrate species on Ion Proton with PI chips: Arctic charr Salvelinus alpinus, European whitefish Coregonus lavaretus and common lizard Zootoca vivipara. This resulted in ~962 million single‐end reads overall and a mean of ~74 million reads per library. We filtered the genomic data using Stacks, a bioinformatic tool to process RAD sequencing data. On average, we obtained ~11 000 polymorphic loci per library of 6–30 individuals. We validate our new method by technical and biological replication, by reconstructing phylogenetic relationships, and using a hybrid genetic cross to track genomic variants. Finally, we discuss the differences between using the different sequencing platforms in the context of RAD sequencing, assessing possible advantages and disadvantages. We show that our protocol can be used for Ion semiconductor sequencing platforms for the rapid and cost‐effective generation of variable and reproducible genetic markers.     

Microhaplotypes provide increased power from short‐read DNA sequences for relationship inference

The accelerating rate at which DNA sequence data are now generated by high‐throughput sequencing instruments provides both opportunities and challenges for population genetic and ecological investigations of animals and plants. We show here how the common practice of calling genotypes from a single SNP per sequenced region ignores substantial additional information in the phased short‐read sequences that are provided by these sequencing instruments. We target sequenced regions with multiple SNPs in kelp rockfish (Sebastes atrovirens) to determine “microhaplotypes” and then call these microhaplotypes as alleles at each locus. We then demonstrate how these multi‐allelic marker data from such loci dramatically increase power for relationship inference. The microhaplotype approach decreases false‐positive rates by several orders of magnitude, relative to calling bi‐allelic SNPs, for two challenging analytical procedures, full‐sibling and single parent–offspring pair identification. We also show how the identification of half‐sibling pairs requires so much data that physical linkage becomes a consideration, and that most published studies that attempt to do so are dramatically underpowered. The advent of phased short‐read DNA sequence data, in conjunction with emerging analytical tools for their analysis, promises to improve efficiency by reducing the number of loci necessary for a particular level of statistical confidence, thereby lowering the cost of data collection and reducing the degree of physical linkage amongst markers used for relationship estimation. Such advances will facilitate collaborative research and management for migratory and other widespread species.     

Strengths and potential pitfalls of hay transfer for ecological restoration revealed by RAD‐seq analysis in floodplain Arabis species

Achieving high intraspecific genetic diversity is a critical goal in ecological restoration as it increases the adaptive potential and long‐term resilience of populations. Thus, we investigated genetic diversity within and between pristine sites in a fossil floodplain and compared it to sites restored by hay transfer between 1997 and 2014. RAD‐seq genotyping revealed that the stenoecious floodplain species Arabis nemorensis is co‐occurring with individuals that, based on ploidy, ITS‐sequencing and morphology, probably belong to the close relative Arabis sagittata, which has a documented preference for dry calcareous grasslands but has not been reported in floodplain meadows. We show that hay transfer maintains genetic diversity for both species. Additionally, in A. sagittata, transfer from multiple genetically isolated pristine sites resulted in restored sites with increased diversity and admixed local genotypes. In A. nemorensis, transfer did not create novel admixture dynamics because genetic diversity between pristine sites was less differentiated. Thus, the effects of hay transfer on genetic diversity also depend on the genetic make‐up of the donor communities of each species, especially when local material is mixed. Our results demonstrate the efficiency of hay transfer for habitat restoration and emphasize the importance of prerestoration characterization of microgeographic patterns of intraspecific diversity of the community to guarantee that restoration practices reach their goal, that is maximize the adaptive potential of the entire restored plant community. Overlooking these patterns may alter the balance between species in the community. Additionally, our comparison of summary statistics obtained from de novo‐ and reference‐based RAD‐seq pipelines shows that the genomic impact of restoration can be reliably monitored in species lacking prior genomic knowledge.     

How quantitative is metabarcoding: A meta‐analytical approach

Metabarcoding has been used in a range of ecological applications such as taxonomic assignment, dietary analysis and the analysis of environmental DNA. However, after a decade of use in these applications there is little consensus on the extent to which proportions of reads generated corresponds to the original proportions of species in a community. To quantify our current understanding, we conducted a structured review and meta‐analysis. The analysis suggests that a weak quantitative relationship may exist between the biomass and sequences produced (slope = 0.52 ± 0.34, p < 0.01), albeit with a large degree of uncertainty. None of the tested moderators, sequencing platform type, the number of species used in a trial or the source of DNA, were able to explain the variance. Our current understanding of the factors affecting the quantitative performance of metabarcoding is still limited: additional research is required before metabarcoding can be confidently utilized for quantitative applications. Until then, we advocate the inclusion of mock communities when metabarcoding as this facilitates direct assessment of the quantitative ability of any given study.