The macroevolutionary singularity of snakes

Snakes and lizards (Squamata) represent a third of terrestrial vertebrates and exhibit spectacular innovations in locomotion, feeding, and sensory processing. However, the evolutionary drivers of this radiation remain poorly known. We infer potential causes and ultimate consequences of squamate macroevolution by combining individual-based natural history observations (>60,000 animals) with a comprehensive time-calibrated phylogeny that we anchored with genomic data (5400 loci) from 1018 species. Due to shifts in the dynamics of speciation and phenotypic evolution, snakes have transformed the trophic structure of animal communities through the recurrent origin and diversification of specialized predatory strategies. Squamate biodiversity reflects a legacy of singular events that occurred during the early history of snakes and reveals the impact of historical contingency on vertebrate biodiversity.

Reference genome of an iconic lizard in western North America, Blainville's horned lizard Phrynosoma blainvillii

Genome assemblies are increasingly being used to identify adaptive genetic variation that can help prioritize the population management of protected species. This approach may be particularly relevant to species like Blainville's horned lizard, Phrynosoma blainvillii, due to its specialized diet on noxious harvester ants, numerous adaptative traits for avoiding predation (e.g. cranial horns, dorsoventrally compressed body, cryptic coloration, and blood squirting from the orbital sinuses), and status as Species of Special Concern in California. Rangewide decline since the early 20th century, the basis of its conservation status, has been driven mainly by habitat conversion, over-collecting, and invasion of a non-native ant that displaces its native ant prey base. Here, we report on a scaffold-level genome assembly for P. blainvillii as part of the California Conservation Genomics Project (CCGP), produced using Pacific Biosciences HiFi long reads and Hi-C chromatin-proximity sequencing technology. The de novo assembly has 78 scaffolds, a total length of ~2.21 Gb, a scaffold N50 length of ~352 Mb, and BUSCO score of 97.4%. This is the second species of Phrynosoma for which a reference genome has been assembled and represents a considerable improvement in terms of contiguity and completeness. Combined with the landscape genomics data being compiled by the CCGP, this assembly will help strategize efforts to maintain and/or restore local genetic diversity, where interventions like genetic rescue, translocation, and strategic land preservation may be the only means by which P. blainvillii and other low-vagility species can survive in the fragmented habitats of California.

A chromosome-level genome assembly and annotation of the desert horned lizard, Phrynosoma platyrhinos, provides insight into chromosomal rearrangements among reptiles

BACKGROUND

The increasing number of chromosome-level genome assemblies has advanced our knowledge and understanding of macroevolutionary processes. Here, we introduce the genome of the desert horned lizard, Phrynosoma platyrhinos, an iguanid lizard occupying extreme desert conditions of the American southwest. We conduct analysis of the chromosomal structure and composition of this species and compare these features across genomes of 12 other reptiles (5 species of lizards, 3 snakes, 3 turtles, and 1 bird).

FINDINGS

The desert horned lizard genome was sequenced using Illumina paired-end reads and assembled and scaffolded using Dovetail Genomics Hi-C and Chicago long-range contact data. The resulting genome assembly has a total length of 1,901.85 Mb, scaffold N50 length of 273.213 Mb, and includes 5,294 scaffolds. The chromosome-level assembly is composed of 6 macrochromosomes and 11 microchromosomes. A total of 20,764 genes were annotated in the assembly. GC content and gene density are higher for microchromosomes than macrochromosomes, while repeat element distributions show the opposite trend. Pathway analyses provide preliminary evidence that microchromosome and macrochromosome gene content are functionally distinct. Synteny analysis indicates that large microchromosome blocks are conserved among closely related species, whereas macrochromosomes show evidence of frequent fusion and fission events among reptiles, even between closely related species.

CONCLUSIONS

Our results demonstrate dynamic karyotypic evolution across Reptilia, with frequent inferred splits, fusions, and rearrangements that have resulted in shuffling of chromosomal blocks between macrochromosomes and microchromosomes. Our analyses also provide new evidence for distinct gene content and chromosomal structure between microchromosomes and macrochromosomes within reptiles.

Genome-scale data reveal deep lineage divergence and a complex demographic history in the Texas horned lizard (Phrynosoma cornutum) throughout the southwestern and central US

The southwestern and central United States serve as an ideal region to test alternative hypotheses regarding biotic diversification. Genomic data can now be combined with sophisticated computational models to quantify the impacts of paleoclimate change, geographic features, and habitat heterogeneity on spatial patterns of genetic diversity. In this study, we combine thousands of genotyping-by-sequencing (GBS) loci with mtDNA sequences (ND1) from the Texas horned lizard (Phrynosoma cornutum) to quantify relative support for different catalysts of diversification. Phylogenetic and clustering analyses of the GBS data indicate support for at least three primary populations. The spatial distribution of populations appears concordant with habitat type, with desert populations in AZ and NM showing the largest genetic divergence from the remaining populations. The mtDNA data also support a divergent desert population, but other relationships differ and suggest mtDNA introgression. Genotype–environment association with bioclimatic variables supports divergence along precipitation gradients more than along temperature gradients. Demographic analyses support a complex history, with introgression and gene flow playing an important role during diversification. Bayesian multispecies coalescent analyses with introgression (MSci) analyses also suggest that gene flow occurred between populations. Paleo-species distribution models support two southern refugia that geographically correspond to contemporary lineages. We find that divergence times are underestimated and population sizes are overestimated when introgression occurred and is ignored in coalescent analyses, and furthermore, inference of ancient introgression events and demographic history is sensitive to inclusion of a single recently admixed sample. Our analyses cannot refute the riverine barrier or glacial refugia hypotheses. Results also suggest that populations are continuing to diverge along habitat gradients. Finally, the strong evidence of admixture, gene flow, and mtDNA introgression among populations suggests that P. cornutum should be considered a single widespread species under the General Lineage Species Concept.

Phylogenetic and morphological influence on habitat choice in moisture-harvesting horned lizards (Phrynosoma spp.)

In previous studies, the superhydrophilic skin of moisture-harvesting lizards has been linked to the morphological traits of the lizards' integument, that is, the occurrence of honeycomb-shaped microstructures. Interestingly, these structures can also cover the skin of lizards inhabiting wet habitats. We therefore tested the influence of the microstructures' main features on the habitat choice and wettability in the genus Phrynosoma. The genus Phrynosoma comprises moisture-harvesting species as well as nonspecialists. Lizards of this genus inhabit large areas of North America with diverse climatic conditions. Remarkably, the differences in the manifestation of microstructures are just as versatile as their surroundings. The phylogeny of the lizards as well as the depth of their ventral microstructures, though independent of each other, correlated with the precipitation in their respective habitat. All other morphological traits, as well as the skin's wettability itself, could not predict the habitat of Phrynosoma species. Hence, it is unlikely that the microstructure influences the wettability, at least directly. Hence, we presume an indirect influence for the following reasons: (a) As the ventral side cannot get wet by rain, but the belly could easily interact with a wet surface, the microstructure might facilitate water absorption from wet soil following precipitation. (b) We found the number of dorsal microstructures to be linked to the occurrence of silt in the habitat. In our study, we observed scales being heavily contaminated, most likely with a mixture of dead skin (after shedding) and silt. As many lizards burrow themselves or even shovel sand onto their backs, deploying the substrate might be a mechanism to increase the skin's wettability.

Phylogenomic Assessment of Biodiversity Using a Reference-Based Taxonomy: An Example With Horned Lizards (Phrynosoma)

Phylogenomic investigations of biodiversity facilitate the detection of fine-scale population genetic structure and the demographic histories of species and populations. However, determining whether or not the genetic divergence measured among populations reflects species-level differentiation remains a central challenge in species delimitation. One potential solution is to compare genetic divergence between putative new species with other closely related species, sometimes referred to as a reference-based taxonomy. To be described as a new species, a population should be at least as divergent as other species. Here, we develop a reference-based taxonomy for Horned Lizards (Phrynosoma; 17 species) using phylogenomic data (ddRADseq data) to provide a framework for delimiting species in the Greater Short-horned Lizard species complex (P. hernandesi). Previous species delimitation studies of this species complex have produced conflicting results, with morphological data suggesting that P. hernandesi consists of five species, whereas mitochondrial DNA support anywhere from 1 to 10 + species. To help address this conflict, we first estimated a time-calibrated species tree for P. hernandesi and close relatives using SNP data. These results support the paraphyly of P. hernandesi; we recommend the recognition of two species to promote a taxonomy that is consistent with species monophyly. There is strong evidence for three populations within P. hernandesi, and demographic modeling and admixture analyses suggest that these populations are not reproductively isolated, which is consistent with previous morphological analyses that suggest hybridization could be common. Finally, we characterize the population-species boundary by quantifying levels of genetic divergence for all 18 Phrynosoma species. Genetic divergence measures for western and southern populations of P. hernandesi failed to exceed those of other Phrynosoma species, but the relatively small population size estimated for the northern population causes it to appear as a relatively divergent species. These comparisons underscore the difficulties associated with putting a reference-based approach to species delimitation into practice. Nevertheless, the reference-based approach offers a promising framework for the consistent assessment of biodiversity within clades of organisms with similar life histories and ecological traits.

Do Alignment and Trimming Methods Matter for Phylogenomic (UCE) Analyses?

Alignment is a crucial issue in molecular phylogenetics because different alignment methods can potentially yield very different topologies for individual genes. But it is unclear if the choice of alignment methods remains important in phylogenomic analyses, which incorporate data from hundreds or thousands of genes. For example, problematic biases in alignment might be multiplied across many loci, whereas alignment errors in individual genes might become irrelevant. The issue of alignment trimming (i.e., removing poorly aligned regions or missing data from individual genes) is also poorly explored. Here, we test the impact of 12 different combinations of alignment and trimming methods on phylogenomic analyses. We compare these methods using published phylogenomic data from ultraconserved elements (UCEs) from squamate reptiles (lizards and snakes), birds, and tetrapods. We compare the properties of alignments generated by different alignment and trimming methods (e.g., length, informative sites, missing data). We also test whether these data sets can recover well-established clades when analyzed with concatenated (RAxML) and species-tree methods (ASTRAL-III), using the full data ($\sim $5000 loci) and subsampled data sets (10% and 1% of loci). We show that different alignment and trimming methods can significantly impact various aspects of phylogenomic data sets (e.g., length, informative sites). However, these different methods generally had little impact on the recovery and support values for well-established clades, even across very different numbers of loci. Nevertheless, our results suggest several "best practices" for alignment and trimming. Intriguingly, the choice of phylogenetic methods impacted the phylogenetic results most strongly, with concatenated analyses recovering significantly more well-established clades (with stronger support) than the species-tree analyses. [Alignment; concatenated analysis; phylogenomics; sequence length heterogeneity; species-tree analysis; trimming].

Phylogenomic data resolve higher-level relationships within South American Liolaemus lizards

Phylogenomic approaches now generate hundreds of loci representative of the whole genome that can be used for phylogenetic analyses. The South American lizard genus Liolaemus is the most species-rich vertebrate radiation from temperate zones (more than 265 described species), yet most higher-level phylogenetic relationships within Liolaemus remain poorly resolved. In this study, we used 584 nuclear loci collected using targeted sequenced capture to estimate the phylogenetic relationships among 26 species representing the two subgenera within Liolaemus (Eulaemus + Liolaemus), and all major groups within Eulaemus. Previous molecular and morphological-based phylogenetic analyses of Eulaemus based on a limited number of characters resolved few higher-level relationships, although one point of agreement is that the early divergence within Eulaemus corresponds to the lineomaculatus section, followed by the diversification of eight main clades that are strongly supported and recognized. Liolaemus probably experienced relatively rapid divergences during parts of its evolutionary history, and a phylogenomic approach was used to resolve the relationships among the major groups. The new analyses presented here support the division of Liolaemus into two subgenera, and resolve relationships among many of the major clades of Eulaemus with strong support. A Bayesian divergence dating analysis using 44 protein-coding genes provides an estimation of the split of the two Liolaemus subgenera of approximately 19,7 ma (95% HPD = 16,94-23,04), while diversification within Eulaemus started at 15,05 ma (95% HPD = 12,94 - 17,59) among the L. lineomaculatus and the L. montanus series by Mid Miocene. A novel phylogenetic network analyses for SNP data identified two hybridizing edges among different groups of Eulaemus at different points in time. Having a solid phylogenetic hypothesis of the main Eulaemus clades opens new opportunities to test a variety of macroevolutionary questions for this unique radiation.

Fully automated sequence alignment methods are comparable to, and much faster than, traditional methods in large data sets: an example with hepatitis B virus

Aligning sequences for phylogenetic analysis (multiple sequence alignment; MSA) is an important, but increasingly computationally expensive step with the recent surge in DNA sequence data. Much of this sequence data is publicly available, but can be extremely fragmentary (i.e., a combination of full genomes and genomic fragments), which can compound the computational issues related to MSA. Traditionally, alignments are produced with automated algorithms and then checked and/or corrected "by eye" prior to phylogenetic inference. However, this manual curation is inefficient at the data scales required of modern phylogenetics and results in alignments that are not reproducible. Recently, methods have been developed for fully automating alignments of large data sets, but it is unclear if these methods produce alignments that result in compatible phylogenies when compared to more traditional alignment approaches that combined automated and manual methods. Here we use approximately 33,000 publicly available sequences from the hepatitis B virus (HBV), a globally distributed and rapidly evolving virus, to compare different alignment approaches. Using one data set comprised exclusively of whole genomes and a second that also included sequence fragments, we compared three MSA methods: (1) a purely automated approach using traditional software, (2) an automated approach including by eye manual editing, and (3) more recent fully automated approaches. To understand how these methods affect phylogenetic results, we compared resulting tree topologies based on these different alignment methods using multiple metrics. We further determined if the monophyly of existing HBV genotypes was supported in phylogenies estimated from each alignment type and under different statistical support thresholds. Traditional and fully automated alignments produced similar HBV phylogenies. Although there was variability between branch support thresholds, allowing lower support thresholds tended to result in more differences among trees. Therefore, differences between the trees could be best explained by phylogenetic uncertainty unrelated to the MSA method used. Nevertheless, automated alignment approaches did not require human intervention and were therefore considerably less time-intensive than traditional approaches. Because of this, we conclude that fully automated algorithms for MSA are fully compatible with older methods even in extremely difficult to align data sets. Additionally, we found that most HBV diagnostic genotypes did not correspond to evolutionarily-sound groups, regardless of alignment type and support threshold. This suggests there may be errors in genotype classification in the database or that HBV genotypes may need a revision.

Cryptic diversity in the Mexican highlands: Thousands of UCE loci help illuminate phylogenetic relationships, species limits and divergence times of montane rattlesnakes (Viperidae: Crotalus)

With the continued adoption of genome-scale data in evolutionary biology comes the challenge of adequately harnessing the information to make accurate phylogenetic inferences. Coalescent-based methods of species tree inference have become common, and concatenation has been shown in simulation to perform well, particularly when levels of incomplete lineage sorting are low. However, simulation conditions are often overly simplistic, leaving empiricists with uncertainty regarding analytical tools. We use a large ultraconserved element data set (>3,000 loci) from rattlesnakes of the Crotalus triseriatus group to delimit lineages and estimate species trees using concatenation and several coalescent-based methods. Unpartitioned and partitioned maximum likelihood and Bayesian analysis of the concatenated matrix yield a topology identical to coalescent analysis of a subset of the data in bpp. ASTRAL analysis on a subset of the more variable loci also results in a tree consistent with concatenation and bpp, whereas the SVDquartets phylogeny differs at additional nodes. The size of the concatenated matrix has a strong effect on species tree inference using SVDquartets, warranting additional investigation on optimal data characteristics for this method. Species delimitation analyses suggest up to 16 unique lineages may be present within the C. triseriatus group, with divergences occurring during the Neogene and Quaternary. Network analyses suggest hybridization within the group is relatively rare. Altogether, our results reaffirm the Mexican highlands as a biodiversity hotspot and suggest that coalescent-based species tree inference on data subsets can provide a strongly supported species tree consistent with concatenation of all loci with a large amount of missing data.

A genomic evaluation of taxonomic trends through time in coast horned lizards (genus Phrynosoma)

Determining the boundaries between species and deciding when to describe new species are challenging practices that are particularly difficult in groups with high levels of geographic variation. The coast horned lizards (Phrynosoma blainvillii, Phrynosoma cerroense and P. coronatum) have an extensive geographic distribution spanning many distinctive ecological regions ranging from northern California to the Cape Region of Baja California, Mexico, and populations differ substantially with respect to external morphology across much of this range. The number of taxa recognized in the group has been reevaluated by herpetologists over 20 times during the last 180 years, and typically without the aid of explicit species delimitation methods, resulting in a turbulent taxonomy containing anywhere from one to seven taxa. In this study, we evaluate taxonomic trends through time by ranking 15 of these species delimitation models (SDMs) using coalescent analyses of nuclear loci and SNPs in a Bayesian model comparison framework. Species delimitation models containing more species were generally favoured by Bayesian model selection; however, several three-species models outperformed some four- and five-species SDMs, and the top-ranked model, which contained five species, outperformed all SDMs containing six species. Model performance peaked in the 1950s based on marginal likelihoods estimated from nuclear loci and SNPs. Not surprisingly, SDMs based on genetic data outperformed morphological taxonomies when using genetic data alone to evaluate models. The de novo estimation of population structure favours a three-population model that matches the currently recognized integrative taxonomy containing three species. We discuss why Bayesian model selection might favour models containing more species, and why recognizing more than three species might be warranted.

Phylogenomics of a rapid radiation: is chromosomal evolution linked to increased diversification in north american spiny lizards (Genus Sceloporus)?

Background

Resolving the short phylogenetic branches that result from rapid evolutionary diversification often requires large numbers of loci. We collected targeted sequence capture data from 585 nuclear loci (541 ultraconserved elements and 44 protein-coding genes) to estimate the phylogenetic relationships among iguanian lizards in the North American genus Sceloporus. We tested for diversification rate shifts to determine if rapid radiation in the genus is correlated with chromosomal evolution.

Results

The phylogenomic trees that we obtained for Sceloporus using concatenation and coalescent-based species tree inference provide strong support for the monophyly and interrelationships among nearly all major groups. The diversification analysis supported one rate shift on the Sceloporus phylogeny approximately 20–25 million years ago that is associated with the doubling of the speciation rate from 0.06 species/million years (Ma) to 0.15 species/Ma. The posterior probability for this rate shift occurring on the branch leading to the Sceloporus species groups exhibiting increased chromosomal diversity is high (posterior probability = 0.997).

Conclusions

Despite high levels of gene tree discordance, we were able to estimate a phylogenomic tree for Sceloporus that solves some of the taxonomic problems caused by previous analyses of fewer loci. The taxonomic changes that we propose using this new phylogenomic tree help clarify the number and composition of the major species groups in the genus. Our study provides new evidence for a putative link between chromosomal evolution and the rapid divergence and radiation of Sceloporus across North America.

Electronic supplementary material

The online version of this article (doi:10.1186/s12862-016-0628-x) contains supplementary material, which is available to authorized users.

Short Tree, Long Tree, Right Tree, Wrong Tree: New Acquisition Bias Corrections for Inferring SNP Phylogenies

Single nucleotide polymorphisms (SNPs) are useful markers for phylogenetic studies owing in part to their ubiquity throughout the genome and ease of collection. Restriction site associated DNA sequencing (RADseq) methods are becoming increasingly popular for SNP data collection, but an assessment of the best practises for using these data in phylogenetics is lacking. We use computer simulations, and new double digest RADseq (ddRADseq) data for the lizard family Phrynosomatidae, to investigate the accuracy of RAD loci for phylogenetic inference. We compare the two primary ways RAD loci are used during phylogenetic analysis, including the analysis of full sequences (i.e., SNPs together with invariant sites), or the analysis of SNPs on their own after excluding invariant sites. We find that using full sequences rather than just SNPs is preferable from the perspectives of branch length and topological accuracy, but not of computational time. We introduce two new acquisition bias corrections for dealing with alignments composed exclusively of SNPs, a conditional likelihood method and a reconstituted DNA approach. The conditional likelihood method conditions on the presence of variable characters only (the number of invariant sites that are unsampled but known to exist is not considered), while the reconstituted DNA approach requires the user to specify the exact number of unsampled invariant sites prior to the analysis. Under simulation, branch length biases increase with the amount of missing data for both acquisition bias correction methods, but branch length accuracy is much improved in the reconstituted DNA approach compared to the conditional likelihood approach. Phylogenetic analyses of the empirical data using concatenation or a coalescent-based species tree approach provide strong support for many of the accepted relationships among phrynosomatid lizards, suggesting that RAD loci contain useful phylogenetic signal across a range of divergence times despite the presence of missing data. Phylogenetic analysis of RAD loci requires careful attention to model assumptions, especially if downstream analyses depend on branch lengths.