Rethinking eco-evo studies of gene expression for non-model organisms in the genomic era.

Recent advances in genomic technology, including the rapid development of long-read sequencing technology and single-cell RNA-sequencing methods, are poised to significantly expand the kinds of studies that are feasible in ecological genomics. In this perspective, we review these new technologies and discuss their potential impact on gene expression studies in non-model organisms. Although traditional RNA-sequencing methods have been an extraordinarily powerful tool to apply functional genomics in an ecological context, bulk RNA-seq approaches often rely on de novo transcriptome assembly, and cannot capture expression changes in rare cell populations or distinguish shifts in cell type abundance. Advancements in genome assembly technology, particularly long-read sequencing, and improvements in the scalability of single-cell RNA-sequencing (scRNA-seq) offer unprecedented resolution in understanding cellular heterogeneity and gene regulation. We discuss the potential of these technologies to enable disentangling differential gene regulation from cell type composition differences and uncovering subtle expression patterns masked by bulk RNA-seq. The integration of these approaches provides a more nuanced understanding of the ecological and evolutionary dynamics of gene expression, paving the way for refined models and deeper insights into the generation of biodiversity.

A nuclear genome assembly of an extinct flightless bird, the little bush moa.

We present a draft genome of the little bush moa (Anomalopteryx didiformis)-one of approximately nine species of extinct flightless birds from Aotearoa, New Zealand-using ancient DNA recovered from a fossil bone from the South Island. We recover a complete mitochondrial genome at 249.9× depth of coverage and almost 900 megabases of a male moa nuclear genome at ~4 to 5× coverage, with sequence contiguity sufficient to identify more than 85% of avian universal single-copy orthologs. We describe a diverse landscape of transposable elements and satellite repeats, estimate a long-term effective population size of ~240,000, identify a diverse suite of olfactory receptor genes and an opsin repertoire with sensitivity in the ultraviolet range, show that the wingless moa phenotype is likely not attributable to gene loss or pseudogenization, and identify potential function-altering coding sequence variants in moa that could be synthesized for future functional assays. This genomic resource should support further studies of avian evolution and morphological divergence.

A phylogenetic method linking nucleotide substitution rates to rates of continuous trait evolution.

Genomes contain conserved non-coding sequences that perform important biological functions, such as gene regulation. We present a phylogenetic method, PhyloAcc-C, that associates nucleotide substitution rates with changes in a continuous trait of interest. The method takes as input a multiple sequence alignment of conserved elements, continuous trait data observed in extant species, and a background phylogeny and substitution process. Gibbs sampling is used to assign rate categories (background, conserved, accelerated) to lineages and explore whether the assigned rate categories are associated with increases or decreases in the rate of trait evolution. We test our method using simulations and then illustrate its application using mammalian body size and lifespan data previously analyzed with respect to protein coding genes. Like other studies, we find processes such as tumor suppression, telomere maintenance, and p53 regulation to be related to changes in longevity and body size. In addition, we also find that skeletal genes, and developmental processes, such as sprouting angiogenesis, are relevant.

A Fast, Reproducible, High-throughput Variant Calling Workflow for Population Genomics.

The increasing availability of genomic resequencing data sets and high-quality reference genomes across the tree of life present exciting opportunities for comparative population genomic studies. However, substantial challenges prevent the simple reuse of data across different studies and species, arising from variability in variant calling pipelines, data quality, and the need for computationally intensive reanalysis. Here, we present snpArcher, a flexible and highly efficient workflow designed for the analysis of genomic resequencing data in nonmodel organisms. snpArcher provides a standardized variant calling pipeline and includes modules for variant quality control, data visualization, variant filtering, and other downstream analyses. Implemented in Snakemake, snpArcher is user-friendly, reproducible, and designed to be compatible with high-performance computing clusters and cloud environments. To demonstrate the flexibility of this pipeline, we applied snpArcher to 26 public resequencing data sets from nonmammalian vertebrates. These variant data sets are hosted publicly to enable future comparative population genomic analyses. With its extensibility and the availability of public data sets, snpArcher will contribute to a broader understanding of genetic variation across species by facilitating the rapid use and reuse of large genomic data sets.

Convergence in hearing-related genes between echolocating birds and mammals.

Echolocation, the detection of objects by means of sound waves, has evolved independently in diverse animals. Echolocators include not only mammals such as toothed whales and yangochiropteran and rhinolophoid bats but also Rousettus fruit bats, as well as two bird lineages, oilbirds and swiftlets. In whales and yangochiropteran and rhinolophoid bats, positive selection and molecular convergence has been documented in key hearing-related genes, such as prestin (SLC26A5), but few studies have examined these loci in other echolocators. Here, we examine patterns of selection and convergence in echolocation-related genes in echolocating birds and Rousettus bats. Fewer of these loci were under selection in Rousettus or birds compared with classically recognized echolocators, and elevated convergence (compared to outgroups) was not evident across this gene set. In certain genes, however, we detected convergent substitutions with potential functional relevance, including convergence between Rousettus and classic echolocators in prestin at a site known to affect hair cell electromotility. We also detected convergence between Yangochiroptera, Rhinolophidea, and oilbirds in TMC1, an important mechanosensory transduction channel in vertebrate hair cells, and observed an amino acid change at the same site within the pore domain. Our results suggest that although most proteins implicated in echolocation in specialized mammals may not have been recruited in birds or Rousettus fruit bats, certain hearing-related loci may have undergone convergent functional changes. Investigating adaptations in diverse echolocators will deepen our understanding of this unusual sensory modality.

PhyloAcc-GT: A Bayesian Method for Inferring Patterns of Substitution Rate Shifts on Targeted Lineages Accounting for Gene Tree Discordance.

An important goal of evolutionary genomics is to identify genomic regions whose substitution rates differ among lineages. For example, genomic regions experiencing accelerated molecular evolution in some lineages may provide insight into links between genotype and phenotype. Several comparative genomics methods have been developed to identify genomic accelerations between species, including a Bayesian method called PhyloAcc, which models shifts in substitution rate in multiple target lineages on a phylogeny. However, few methods consider the possibility of discordance between the trees of individual loci and the species tree due to incomplete lineage sorting, which might cause false positives. Here, we present PhyloAcc-GT, which extends PhyloAcc by modeling gene tree heterogeneity. Given a species tree, we adopt the multispecies coalescent model as the prior distribution of gene trees, use Markov chain Monte Carlo (MCMC) for inference, and design novel MCMC moves to sample gene trees efficiently. Through extensive simulations, we show that PhyloAcc-GT outperforms PhyloAcc and other methods in identifying target lineage-specific accelerations and detecting complex patterns of rate shifts, and is robust to specification of population size parameters. PhyloAcc-GT is usually more conservative than PhyloAcc in calling convergent rate shifts because it identifies more accelerations on ancestral than on terminal branches. We apply PhyloAcc-GT to two examples of convergent evolution: flightlessness in ratites and marine mammal adaptations, and show that PhyloAcc-GT is a robust tool to identify shifts in substitution rate associated with specific target lineages while accounting for incomplete lineage sorting.

Expanding on Our Knowledge of Ecdysozoan Genomes: A Contiguous Assembly of the Meiofaunal Priapulan Tubiluchus corallicola.

Genomic data for priapulans are limited to a single species, restricting broad comparative analyses and thorough interrogation of questions spanning phylogenomics, ecdysozoan physiology, and development. To help fill this void, we present here a high-quality priapulan genome for the meiofaunal species Tubiluchus corallicola. Our assembly combines Nanopore and Illumina sequencing technologies and makes use of a whole-genome amplification, to generate enough DNA to sequence this small meiofaunal species. We generated a moderately contiguous assembly (2,547 scaffolds), with a high level of completeness (metazoan BUSCOs n = 954, single-copy complete = 89.6%, duplicated = 3.9%, fragmented = 3.5%, and missing = 3.0%). We then screened the genome for homologs of the Halloween genes, key genes implicated in the ecdysis (molting) pathway of arthropods, recovering a putative homolog of shadow. The presence of a shadow ortholog in two priapulan genomes suggests that the Halloween genes may not have evolved in a stepwise manner in Panarthropoda, as previously thought, but may have a deeper origin at the base of Ecdysozoa.

Interpreting the pervasive observation of U-shaped Site Frequency Spectra.

The standard neutral model of molecular evolution has traditionally been used as the null model for population genomics. We gathered a collection of 45 genome-wide site frequency spectra from a diverse set of species, most of which display an excess of low and high frequency variants compared to the expectation of the standard neutral model, resulting in U-shaped spectra. We show that multiple merger coalescent models often provide a better fit to these observations than the standard Kingman coalescent. Hence, in many circumstances these under-utilized models may serve as the more appropriate reference for genomic analyses. We further discuss the underlying evolutionary processes that may result in the widespread U-shape of frequency spectra.

Sizing Up the Onychophoran Genome: Repeats, Introns, and Gene Family Expansion Contribute to Genome Gigantism in Epiperipatus broadwayi.

Genome assemblies are growing at an exponential rate and have proved indispensable for studying evolution but the effort has been biased toward vertebrates and arthropods with a particular focus on insects. Onychophora or velvet worms are an ancient group of cryptic, soil dwelling worms noted for their unique mode of prey capture, biogeographic patterns, and diversity of reproductive strategies. They constitute a poorly understood phylum of exclusively terrestrial animals that is sister group to arthropods. Due to this phylogenetic position, they are crucial in understanding the origin of the largest phylum of animals. Despite their significance, there is a paucity of genomic resources for the phylum with only one highly fragmented and incomplete genome publicly available. Initial attempts at sequencing an onychophoran genome proved difficult due to its large genome size and high repeat content. However, leveraging recent advances in long-read sequencing technology, we present here the first annotated draft genome for the phylum. With a total size of 5.6Gb, the gigantism of the Epiperipatus broadwayi genome arises from having high repeat content, intron size inflation, and extensive gene family expansion. Additionally, we report a previously unknown diversity of onychophoran hemocyanins that suggests the diversification of copper-mediated oxygen carriers occurred independently in Onychophora after its split from Arthropoda, parallel to the independent diversification of hemocyanins in each of the main arthropod lineages.

Genomic signature of sexual reproduction in the bdelloid rotifer Macrotrachella quadricornifera.

Bdelloid rotifers, common freshwater invertebrates of ancient origin and worldwide distribution have long been thought to be entirely asexual, being the principal exception to the view that in eukaryotes the loss of sex leads to early extinction. That bdelloids are facultatively sexual is shown by a study of allele sharing within a group of closely related bdelloids of the species Macrotrachella quadricornifera, supporting the view that sexual reproduction is essential for long-term success in all eukaryotes.

HieRFIT: a hierarchical cell type classification tool for projections from complex single-cell atlas datasets.

MOTIVATION: The emergence of single-cell RNA sequencing (scRNA-seq) has led to an explosion in novel methods to study biological variation among individual cells, and to classify cells into functional and biologically meaningful categories. RESULTS: Here, we present a new cell type projection tool, Hierarchical Random Forest for Information Transfer (HieRFIT), based on hierarchical random forests. HieRFIT uses a priori information about cell type relationships to improve classification accuracy, taking as input a hierarchical tree structure representing the class relationships, along with the reference data. We use an ensemble approach combining multiple random forest models, organized in a hierarchical decision tree structure. We show that our hierarchical classification approach improves accuracy and reduces incorrect predictions especially for inter-dataset tasks which reflect real-life applications. We use a scoring scheme that adjusts probability distributions for candidate class labels and resolves uncertainties while avoiding the assignment of cells to incorrect types by labeling cells at internal nodes of the hierarchy when necessary. AVAILABILITY AND IMPLEMENTATION: HieRFIT is implemented as an R package, and it is available at (https://github.com/yasinkaymaz/HieRFIT/releases/tag/v1.0.0). SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

The genome of the stable fly, Stomoxys calcitrans, reveals potential mechanisms underlying reproduction, host interactions, and novel targets for pest control.

BACKGROUND: The stable fly, Stomoxys calcitrans, is a major blood-feeding pest of livestock that has near worldwide distribution, causing an annual cost of over $2 billion for control and product loss in the USA alone. Control of these flies has been limited to increased sanitary management practices and insecticide application for suppressing larval stages. Few genetic and molecular resources are available to help in developing novel methods for controlling stable flies. RESULTS: This study examines stable fly biology by utilizing a combination of high-quality genome sequencing and RNA-Seq analyses targeting multiple developmental stages and tissues. In conjunction, 1600 genes were manually curated to characterize genetic features related to stable fly reproduction, vector host interactions, host-microbe dynamics, and putative targets for control. Most notable was characterization of genes associated with reproduction and identification of expanded gene families with functional associations to vision, chemosensation, immunity, and metabolic detoxification pathways. CONCLUSIONS: The combined sequencing, assembly, and curation of the male stable fly genome followed by RNA-Seq and downstream analyses provide insights necessary to understand the biology of this important pest. These resources and new data will provide the groundwork for expanding the tools available to control stable fly infestations. The close relationship of Stomoxys to other blood-feeding (horn flies and Glossina) and non-blood-feeding flies (house flies, medflies, Drosophila) will facilitate understanding of the evolutionary processes associated with development of blood feeding among the Cyclorrhapha.

Deeply Altered Genome Architecture in the Endoparasitic Flowering Plant Sapria himalayana Griff. (Rafflesiaceae).

Despite more than 2,000-fold variation in genome size, key features of genome architecture are largely conserved across angiosperms. Parasitic plants have elucidated the many ways in which genomes can be modified, yet we still lack comprehensive genome data for species that represent the most extreme form of parasitism. Here, we present the highly modified genome of the iconic endophytic parasite Sapria himalayana Griff. (Rafflesiaceae), which lacks a typical plant body. First, 44% of the genes conserved in eurosids are lost in Sapria, dwarfing previously reported levels of gene loss in vascular plants. These losses demonstrate remarkable functional convergence with other parasitic plants, suggesting a common genetic roadmap underlying the evolution of plant parasitism. Second, we identified extreme disparity in intron size among retained genes. This includes a category of genes with introns longer than any so far observed in angiosperms, nearing 100 kb in some cases, and a second category of genes with exceptionally short or absent introns. Finally, at least 1.2% of the Sapria genome, including both genic and intergenic content, is inferred to be derived from host-to-parasite horizontal gene transfers (HGTs) and includes genes potentially adaptive for parasitism. Focused phylogenomic reconstruction of HGTs reveals a hidden history of former host-parasite associations involving close relatives of Sapria's modern hosts in the grapevine family. Our findings offer a unique perspective into how deeply angiosperm genomes can be altered to fit an extreme form of plant parasitism and demonstrate the value of HGTs as DNA fossils to investigate extinct symbioses.

Error, noise and bias in de novo transcriptome assemblies.

De novo transcriptome assembly is a powerful tool, and has been widely used over the last decade for making evolutionary inferences. However, it relies on two implicit assumptions: that the assembled transcriptome is an unbiased representation of the underlying expressed transcriptome, and that expression estimates from the assembly are good, if noisy approximations of the relative abundance of expressed transcripts. Using publicly available data for model organisms, we demonstrate that, across assembly algorithms and data sets, these assumptions are consistently violated. Bias exists at the nucleotide level, with genotyping error rates ranging from 30% to 83%. As a result, diversity is underestimated in transcriptome assemblies, with consistent underestimation of heterozygosity in all but the most inbred samples. Even at the gene level, expression estimates show wide deviations from map-to-reference estimates, and positive bias at lower expression levels. Standard filtering of transcriptome assemblies improves the robustness of gene expression estimates but leads to the loss of a meaningful number of protein-coding genes, including many that are highly expressed. We demonstrate a computational method, length-rescaled CPM, to partly alleviate noise and bias in expression estimates. Researchers should consider ways to minimize the impact of bias in transcriptome assemblies.

Evaluating single-cell cluster stability using the Jaccard similarity index.

MOTIVATION: One major goal of single-cell RNA sequencing (scRNAseq) experiments is to identify novel cell types. With increasingly large scRNAseq datasets, unsupervised clustering methods can now produce detailed catalogues of transcriptionally distinct groups of cells in a sample. However, the interpretation of these clusters is challenging for both technical and biological reasons. Popular clustering algorithms are sensitive to parameter choices, and can produce different clustering solutions with even small changes in the number of principal components used, the k nearest neighbor and the resolution parameters, among others. RESULTS: Here, we present a set of tools to evaluate cluster stability by subsampling, which can guide parameter choice and aid in biological interpretation. The R package scclusteval and the accompanying Snakemake workflow implement all steps of the pipeline: subsampling the cells, repeating the clustering with Seurat and estimation of cluster stability using the Jaccard similarity index and providing rich visualizations. AVAILABILITYAND IMPLEMENTATION: R package scclusteval: https://github.com/crazyhottommy/scclusteval Snakemake workflow: https://github.com/crazyhottommy/pyflow_seuratv3_parameter Tutorial: https://crazyhottommy.github.io/EvaluateSingleCellClustering/.

Studying Natural Selection in the Era of Ubiquitous Genomes.

A major goal of comparative genomics research is modeling changes in DNA sequences between species to understand the evolutionary forces acting on species differences. Application of these models to a number of species over the past decade has revealed some commonalities across organisms, most notably a consistent role of positive selection in shaping the molecular evolution of the immune system. However, models of DNA sequence evolution also have important limitations that are increasingly being recognized, including issues with data quality and biases caused by simplifying assumptions. While new approaches have begun to address these challenges, ultimately, additional data, such as resequencing data from populations, will provide more power to fully understand the unique evolutionary forces acting on different species. In this review, I summarize the conclusions of recent genome-wide studies of selection, highlight some important challenges to applying these methods to large data sets, and discuss ways forward for the field.

Short paired-end reads trump long single-end reads for expression analysis.

BACKGROUND: Typical experimental design advice for expression analyses using RNA-seq generally assumes that single-end reads provide robust gene-level expression estimates in a cost-effective manner, and that the additional benefits obtained from paired-end sequencing are not worth the additional cost. However, in many cases (e.g., with Illumina NextSeq and NovaSeq instruments), shorter paired-end reads and longer single-end reads can be generated for the same cost, and it is not obvious which strategy should be preferred. Using publicly available data, we test whether short-paired end reads can achieve more robust expression estimates and differential expression results than single-end reads of approximately the same total number of sequenced bases. RESULTS: At both the transcript and gene levels, 2 × 40 paired-end reads unequivocally provide expression estimates that are more highly correlated with 2 × 125 than 1 × 75 reads; in nearly all cases, those correlations are also greater than for 1 × 125, despite the greater total number of sequenced bases for the latter. Across an array of metrics, differential expression tests based upon 2 × 40 consistently outperform those using 1 × 75. CONCLUSION: Researchers seeking a cost-effective approach for gene-level expression analysis should prefer short paired-end reads over a longer single-end strategy. Short paired-end reads will also give reasonably robust expression estimates and differential expression results at the isoform level.

Genomics of host-pathogen interactions: challenges and opportunities across ecological and spatiotemporal scales.

Evolutionary genomics has recently entered a new era in the study of host-pathogen interactions. A variety of novel genomic techniques has transformed the identification, detection and classification of both hosts and pathogens, allowing a greater resolution that helps decipher their underlying dynamics and provides novel insights into their environmental context. Nevertheless, many challenges to a general understanding of host-pathogen interactions remain, in particular in the synthesis and integration of concepts and findings across a variety of systems and different spatiotemporal and ecological scales. In this perspective we aim to highlight some of the commonalities and complexities across diverse studies of host-pathogen interactions, with a focus on ecological, spatiotemporal variation, and the choice of genomic methods used. We performed a quantitative review of recent literature to investigate links, patterns and potential tradeoffs between the complexity of genomic, ecological and spatiotemporal scales undertaken in individual host-pathogen studies. We found that the majority of studies used whole genome resolution to address their research objectives across a broad range of ecological scales, especially when focusing on the pathogen side of the interaction. Nevertheless, genomic studies conducted in a complex spatiotemporal context are currently rare in the literature. Because processes of host-pathogen interactions can be understood at multiple scales, from molecular-, cellular-, and physiological-scales to the levels of populations and ecosystems, we conclude that a major obstacle for synthesis across diverse host-pathogen systems is that data are collected on widely diverging scales with different degrees of resolution. This disparity not only hampers effective infrastructural organization of the data but also data granularity and accessibility. Comprehensive metadata deposited in association with genomic data in easily accessible databases will allow greater inference across systems in the future, especially when combined with open data standards and practices. The standardization and comparability of such data will facilitate early detection of emerging infectious diseases as well as studies of the impact of anthropogenic stressors, such as climate change, on disease dynamics in humans and wildlife.

Evolutionary Dynamics of Sex Chromosomes of Paleognathous Birds.

Standard models of sex chromosome evolution propose that recombination suppression leads to the degeneration of the heterogametic chromosome, as is seen for the Y chromosome in mammals and the W chromosome in most birds. Unlike other birds, paleognaths (ratites and tinamous) possess large nondegenerate regions on their sex chromosomes (PARs or pseudoautosomal regions). It remains unclear why these large PARs are retained over >100 Myr, and how this retention impacts the evolution of sex chromosomes within this system. To address this puzzle, we analyzed Z chromosome evolution and gene expression across 12 paleognaths, several of whose genomes have recently been sequenced. We confirm at the genomic level that most paleognaths retain large PARs. As in other birds, we find that all paleognaths have incomplete dosage compensation on the regions of the Z chromosome homologous to degenerated portions of the W (differentiated regions), but we find no evidence for enrichments of male-biased genes in PARs. We find limited evidence for increased evolutionary rates (faster-Z) either across the chromosome or in differentiated regions for most paleognaths with large PARs, but do recover signals of faster-Z evolution in tinamou species with mostly degenerated W chromosomes, similar to the pattern seen in neognaths. Unexpectedly, in some species, PAR-linked genes evolve faster on average than genes on autosomes, suggested by diverse genomic features to be due to reduced efficacy of selection in paleognath PARs. Our analysis shows that paleognath Z chromosomes are atypical at the genomic level, but the evolutionary forces maintaining largely homomorphic sex chromosomes in these species remain elusive.