Search

Nuclear mitochondrial pseudogenes (numts) may hinder the reconstruction of mtDNA genomes and affect the reliability of mtDNA datasets for phylogenetic and population genetic comparisons. Here, we present the program Numt Parser, which allows for the identification of DNA sequences that likely originate from numt pseudogene DNA. Sequencing reads are classified as originating from either numt or true cytoplasmic mitochondrial (cymt) DNA by direct comparison against cymt and numt reference sequences. Classified reads can then be parsed into cymt or numt datasets. We tested this program using whole genome shotgun-sequenced data from two ancient Cape lions (Panthera leo), because mtDNA is often the marker of choice for ancient DNA studies and the genus Panthera is known to have numt pseudogenes. Numt Parser decreased sequence disagreements that were likely due to numt pseudogene contamination and equalized read coverage across the mitogenome by removing reads that likely originated from numts. We compared the efficacy of Numt Parser to two other bioinformatic approaches that can be used to account for numt contamination. We found that Numt Parser outperformed approaches that rely only on read alignment or Basic Local Alignment Search Tool (BLAST) properties, and was effective at identifying sequences that likely originated from numts while having minimal impacts on the recovery of cymt reads. Numt Parser therefore improves the reconstruction of true mitogenomes, allowing for more accurate and robust biological inferences.

The basal South American notothenioid Eleginops maclovinus (Patagonia blennie or róbalo) occupies a uniquely important phylogenetic position in Notothenioidei as the singular closest sister species to the Antarctic cryonotothenioid fishes. Its genome and the traits encoded therein would be the nearest representatives of the temperate ancestor from which the Antarctic clade arose, providing an ancestral reference for deducing polar derived changes. In this study, we generated a gene- and chromosome-complete assembly of the E. maclovinus genome using long read sequencing and HiC scaffolding. We compared its genome architecture with the more basally divergent Cottoperca gobio and the derived genomes of nine cryonotothenioids representing all five Antarctic families. We also reconstructed a notothenioid phylogeny using 2918 proteins of single-copy orthologous genes from these genomes that reaffirmed E. maclovinus’ phylogenetic position. We additionally curated E. maclovinus’ repertoire of circadian rhythm genes, ascertained their functionality by transcriptome sequencing, and compared its pattern of gene retention with C. gobio and the derived cryonotothenioids. Through reconstructing circadian gene trees, we also assessed the potential role of the retained genes in cryonotothenioids by referencing to the functions of the human orthologs. Our results found E. maclovinus to share greater conservation with the Antarctic clade, solidifying its evolutionary status as the direct sister and best suited ancestral proxy of cryonotothenioids. The high-quality genome of E. maclovinus will facilitate inquiries into cold derived traits in temperate to polar evolution, and conversely on the paths of readaptation to non-freezing habitats in various secondarily temperate cryonotothenioids through comparative genomic analyses.

The basal S. American notothenioidEleginops maclovinus(Patagonia blennie) occupies a uniquely important phylogenetic position in Notothenioidei as the singular closest sister species to the Antarctic cryonotothenioid fishes. Its genome and the traits encoded therein would be nearest representatives of the temperate ancestor from which the Antarctic clade arose, providing an ancestral reference for deducing polar derived changes. In this study, we generated a gene- and chromosome-complete assembly ofE. maclovinusgenome using long read sequencing and HiC scaffolding. We compared its genome architecture with the more basally divergentCottoperca gobioand the derived genomes of nine cryonotothenioids representing all five Antarctic families. We additionally curated its repertoire of circadian rhythm genes, ascertained their functionality by transcriptome sequencing, and compared its pattern of gene retention withC. gobioand the derived cryonotothenioids. Both analyses foundE. maclovinusto share greater conservation with the Antarctic clade, solidifying its evolutionary status as direct sister and best suited ancestral proxy of cryonotothenioids. The high quality genome ofE. maclovinuswill facilitate inquiries into cold derived traits in temperate to polar evolution, as well as inform on the paths of readaptation to non-freezing habitats in various secondarily temperate cryonotothenioids through comparative genomic analyses.

Hybrid zones are dynamic systems where natural selection, sexual selection and other evolutionary forces can act on reshuffled combinations of two distinct genomes. The movement of hybrid zones, individual traits, or both are of particular interest for understanding diversification processes. In a hybrid zone involving two lek-breeding birds (genusManacus), secondary sexual plumage traits ofM. vitellinus,especially bright yellow collars and olive bellies, have introgressed asymmetrically ∼50 km across the genomic center of the zone into populations ofM. candei. Males with yellow collars are preferred by females and are more aggressive than parentalM. candei, suggesting that sexual selection was responsible for the introgression of male traits. To assess the spatial and temporal dynamics of this hybrid zone, we carried out a longitudinal reassessment using historical (1989 - 1994) and contemporary (2017 - 2020) transect samples to survey both morphological and genetic variation. Genome-scale SNP data and several male phenotypic traits show that the genomic center of the zone has remained spatially stable, whereas the olive belly color of malevitellinushas continued to introgress over time, suggesting that sexual selection may be an ongoing evolutionary force in this system.

Library preparation protocols for most sequencing technologies involve PCR amplification of the template DNA, which open the possibility that a given template DNA molecule is sequenced multiple times. Reads arising from this phenomenon, known as PCR duplicates, inflate the cost of sequencing and can jeopardize the reliability of affected experiments. Despite the pervasiveness of this artifact, our understanding of its causes and of its impact on downstream statistical analyses remains essentially empirical. Here, we develop a general quantitative model of amplification distortions in sequencing datasets, which we leverage to investigate the factors controlling the occurrence of PCR duplicates. We show that the PCR duplicate rate is determined by the ratio between library complexity and sequencing depth, and that amplification noise (including in its dependence on the number of PCR cycles) only plays a minor role for this artifact. We confirm our predictions using new and published RAD-seq libraries and provide a method to estimate library complexity and amplification noise in any dataset containing PCR duplicates. We discuss how amplification-related artifacts impact downstream analyses, and in particular genotyping accuracy. The proposed framework unites the numerous observations made on PCR duplicates and will be useful to experimenters of all sequencing technologies where DNA availability is a concern.

Restriction enzymes have been one of the primary tools in the population genetics toolkit for 50 years, being coupled with each new generation of technology to provide a more detailed view into the genetics of natural populations. Restriction site-Associated DNA protocols, which joined enzymes with short-read sequencing technology, have democratized the field of population genomics, providing a means to assay the underlying alleles in scores of populations. More than 10 years on, the technique has been widely applied across the tree of life and served as the basis for many different analysis techniques. Here, we provide a detailed protocol to conduct a RAD analysis from experimental design to de novo analysis-including parameter optimization-as well as reference-based analysis, all in Stacks version 2, which is designed to work with paired-end reads to assemble RAD loci up to 1000 nucleotides in length. The protocol focuses on major points of friction in the molecular approaches and downstream analysis, with special attention given to validating experimental analyses. Finally, the protocol provides several points of departure for further analysis.

Body plans often evolve through changes in the number of repeated parts or serial homologs. Using the butterfly Bicyclus anynana, RiveraColón et al. studied the genetics underlying heritability... The underlying genetic changes that regulate the appearance and disappearance of repeated traits, or serial homologs, remain poorly understood. One hypothesis is that variation in genomic regions flanking master regulatory genes, also known as input–output genes, controls variation in trait number, making the locus of evolution almost predictable. Another hypothesis implicates genetic variation in up- or downstream loci of master control genes. Here, we use the butterfly Bicyclus anynana, a species that exhibits natural variation in eyespot number on the dorsal hindwing, to test these two hypotheses. We first estimated the heritability of dorsal hindwing eyespot number by breeding multiple butterfly families differing in eyespot number and regressing eyespot numbers of offspring on midparent values. We then estimated the number and identity of independent genetic loci contributing to eyespot number variation by performing a genome-wide association study with restriction site-associated DNA sequencing from multiple individuals varying in number of eyespots sampled across a freely breeding laboratory population. We found that dorsal hindwing eyespot number has a moderately high heritability of ∼0.50 and is characterized by a polygenic architecture. Previously identified genomic regions involved in eyespot development, and novel ones, display high association with dorsal hindwing eyespot number, suggesting that homolog number variation is likely determined by regulatory changes at multiple loci that build the trait, and not by variation at single master regulators or input–output genes.

Restriction enzymes have been one of the primary tools in the population genetics toolkit for fifty years, being coupled with each new generation of technology to provide a more detailed view into the genetics of natural populations. Restriction site-Associated DNA protocols, which joined enzymes with short-read sequencing technology have democratized the field of population genomics, providing a means to assay the underlying alleles in scores of populations. More than ten years on, the technique has been widely applied across the tree of life and served as the basis for many different analysis techniques. Here, we provide a detailed protocol to conduct a RAD analysis from experimental design to de novo analysis – including parameter optimization – as well as reference-based analysis, all in Stacks version 2, which is designed to work with paired-end reads to assemble RAD loci up to 1000 nucleotides in length. The protocol focuses on major points of friction in the molecular approaches and downstream analysis, with special attention given to validating experimental analyses. Finally, the protocol provides several points of departure for further analysis.

Morphologically similar species, that is cryptic species, may be similar or quasi-similar owing to the deceleration of morphological evolution and stasis. While the factors underlying the deceleration of morphological evolution or stasis in cryptic species remain unknown, decades of research in the field of paleontology on punctuated equilibrium have originated clear hypotheses. Species are expected to remain morphologically identical in scenarios of shared genetic variation, such as hybridization and incomplete lineage sorting, or in scenarios where bottlenecks reduce genetic variation and constrain the evolution of morphology. Here, focusing on three morphologically similar Stygocapitella species, we employ a whole-genome amplification method (WGA) coupled with double-digestion restriction-site associated DNA sequencing (ddRAD) to reconstruct the evolutionary history of the species complex. We explore population structure, use population-level statistics to determine the degree of connectivity between populations and species, and determine the most likely demographic scenarios which generally reject for recent hybridization. We find that the combination of WGA and ddRAD allowed us to obtain genomic-level data from microscopic eukaryotes (∼1 millimetre) opening up opportunities for those working with population genomics and phylogenomics in such taxa. The three species share genetic variance, likely from incomplete lineage sorting and ancient admixture. We speculate that the degree of shared variation might underlie morphological similarity in the Atlantic species complex.

1. The restriction site-associated DNA (RADseq) family of protocols involves digesting DNA and sequencing the region flanking the cut site, thus providing a cost and time-efficient way for obtaining thousands of genomic markers. However, when working with non-model taxa with few genomic resources, optimization of RADseq wet-lab and bioinformatic tools may be challenging, often resulting in allele dropout—that is when a given RADseq locus is not sequenced in one or more individuals resulting in missing data. Additionally, as datasets include divergent taxa, rates of dropout will increase since restriction sites may be lost due to mutation. Mitigating the impacts of allele dropout is crucial, as missing data may lead to incorrect inferences in population genetics and phylogenetics. 2. Here, we demonstrate a simple pipeline for the optimization of RADseq datasets which involves partitioning datasets into subgroups, namely by reducing and analysing the dataset at a population or species level. By running the software Stacks at a subgroup level, we were able to reliably identify and remove individuals with high levels of missing data (bad apples) likely stemming from artefacts in library preparation, DNA quality or sequencing artefacts. 3. Removal of the bad apples generally led to an increase in loci and decrease in missing data in the final datasets. 4. The biological interpretability of the data, as measured by the number of retrieved loci and missing data, was considerably increased.

Restriction-site Associated DNA sequencing (RADseq) has become a powerful and versatile tool in modern population genomics, enabling large-scale genomic analyses in otherwise inaccessible biological systems. With its widespread use, different variants on the protocol have been developed to suit specific experimental needs. Researchers face the challenge of choosing the optimal molecular and sequencing protocols for their experimental design, an often-complicated process. Strategic errors can lead to improper data generation that has reduced power to answer biological questions. Here we present RADinitio, simulation software for the selection and optimization of RADseq experiments via the generation of sequencing data that behaves similarly to empirical sources. RADinitio provides an evolutionary simulation of populations, implementation of various RADseq protocols with customizable parameters, and thorough assessment of missing data. Using the software, we test its efficacy using different RAD protocols across several organisms, highlighting the importance of protocol selection on the magnitude and quality of data acquired. Additionally, we test the effects of RAD library preparation and sequencing on allelic dropout, observing that library preparation and sequencing often contributes more to missing alleles than population-level variation.

For half a century population genetics studies have put type II restriction endonucleases to work. Now, coupled with massively-parallel, short-read sequencing, the family of RAD protocols that wields these enzymes has generated vast genetic knowledge from the natural world. Here, we describe the first software natively capable of using paired-end sequencing to derive short contigs from de novo RAD data. Stacks version 2 employs a de Bruijn graph assembler to build and connect contigs from forward and reverse reads for each de novo RAD locus, which it then uses as a reference for read alignments. The new architecture allows all the individuals in a metapopulation to be considered at the same time as each RAD locus is processed. This enables a Bayesian genotype caller to provide precise SNPs, and a robust algorithm to phase those SNPs into long haplotypes, generating RAD loci that are 400-800 bp in length. To prove its recall and precision, we tested the software with simulated data and compared reference-aligned and de novo analyses of three empirical data sets. Our study shows that the latest version of Stacks is highly accurate and outperforms other software in assembling and genotyping paired-end de novo data sets.

The underlying genetic changes that regulate the appearance and disappearance of repeated traits, or serial homologs, remain poorly understood. One hypothesis is that variation in genomic regions flanking master regulatory genes, also known as input-output genes, controls variation in trait number, making the locus of evolution almost predictable. Other hypotheses implicate genetic variation in up-stream or downstream loci of master control genes. Here, we use the butterflyBicyclus anynana, a species which exhibits natural variation in eyespot number on the dorsal hindwing, to test these two hypotheses. We first estimated the heritability of dorsal hindwing eyespot number by breeding multiple butterfly families differing in eyespot number, and regressing eyespot number of offspring on mid-parent values. We then estimated the number and identity of independent genetic loci contributing to eyespot number variation by performing a genome-wide association study with restriction site-associated DNA Sequencing (RAD-seq) from multiple individuals varying in number of eyespots sampled across a freely breeding lab population. We found that dorsal hindwing eyespot number has a moderately high heritability of approximately 0.50. In addition, multiple loci near previously identified genes involved in eyespot development display high association with dorsal hindwing eyespot number, suggesting that homolog number variation is likely determined by regulatory changes at multiple loci that build the trait and not by variation at single master regulators or input-output genes.Data accessibilityTheBicyclus anynanaPstI RAD-tag sequencing data is available via the Genbank BioProject PRJNA509697. Genotype VCF files will be made available through Figshare upon acceptance.

The underlying genetic changes that regulate the appearance and disappearance of repeated traits, or serial homologs, remain poorly understood. One hypothesis is that variation in genomic regions flanking master regulatory genes, also known as input–output genes, controls variation in trait number, making the locus of evolution almost predictable. Another hypothesis implicates genetic variation in up- or downstream loci of master control genes. Here, we use the butterfly Bicyclus anynana, a species that exhibits natural variation in eyespot number on the dorsal hindwing, to test these two hypotheses. We first estimated the heritability of dorsal hindwing eyespot number by breeding multiple butterfly families differing in eyespot number and regressing eyespot numbers of offspring on midparent values. We then estimated the number and identity of independent genetic loci contributing to eyespot number variation by performing a genome-wide association study with restriction site-associated DNA sequencing from multiple individuals varying in number of eyespots sampled across a freely breeding laboratory population. We found that dorsal hindwing eyespot number has a moderately high heritability of ∼0.50 and is characterized by a polygenic architecture. Previously identified genomic regions involved in eyespot development, and novel ones, display high association with dorsal hindwing eyespot number, suggesting that homolog number variation is likely determined by regulatory changes at multiple loci that build the trait, and not by variation at single master regulators or input–output genes.

The underlying genetic changes that regulate the appearance and disappearance of repeated traits, or serial homologs, remain poorly understood. One hypothesis is that variation in genomic regions flanking master regulatory genes, also known as input-output genes, controls variation in trait number, making the locus of evolution almost predictable. Other hypotheses implicate genetic variation in up-stream or downstream loci of master control genes. Here, we use the butterfly Bicyclus anynana, a species which exhibits natural variation in eyespot number on the dorsal hindwing, to test these two hypotheses. We first estimated the heritability of dorsal hindwing eyespot number by breeding multiple butterfly families differing in eyespot number, and regressing eyespot number of offspring on mid-parent values. We then estimated the number and identity of independent genetic loci contributing to eyespot number variation by performing a genome-wide association study with restriction site-associated DNA Sequencing (RAD-seq) from multiple individuals varying in number of eyespots sampled across a freely breeding lab population. We found that dorsal hindwing eyespot number has a moderately high heritability of approximately 0.50. In addition, multiple loci near previously identified genes involved in eyespot development display high association with dorsal hindwing eyespot number, suggesting that homolog number variation is likely determined by regulatory changes at multiple loci that build the trait and not by variation at single master regulators or input-output genes.