Of contigs and quagmires: next‐generation sequencing pitfalls associated with transcriptomic studies

Comparative transcriptomics of cyprinid minnows and carp in a common wild setting: a resource for ecological genomics in freshwater communities

Comparative transcriptomics can now be conducted on organisms in natural settings, which has greatly enhanced understanding of genome–environment interactions. Here, we demonstrate the utility and potential pitfalls of comparative transcriptomics of wild organisms, with an example from three cyprinid fish species (Teleostei:Cypriniformes). We present extensively filtered and annotated transcriptome assemblies that provide a valuable resource for studies of genome evolution (e.g. polyploidy), ecological and morphological diversification, speciation, and shared and unique responses to environmental variation in cyprinid fishes. Our results and analyses address the following points: (i) ‘essential developmental genes’ are shown to be ubiquitously expressed in a diverse suite of tissues across later ontogenetic stages (i.e. juveniles and adults), making these genes are useful for assessing the quality of transcriptome assemblies, (ii) the influence of microbiomes and other exogenous DNA, (iii) potentially novel, species-specific genes, and (iv) genomic rearrangements (e.g. whole genome duplication). The data we present provide a resource for future comparative work in cypriniform fishes and other taxa across a variety of sub-disciplines, including stress response, morphological diversification, community ecology, ecotoxicology, and climate change.

The genome sequence and insights into the immunogenetics of the bananaquit (Passeriformes: Coereba flaveola)

Avian genomics, especially of non-model species, is in its infancy relative to mammalian genomics. Here, we describe the sequencing, assembly, and annotation of a new avian genome, that of the bananaquit Coereba flaveola (Passeriformes: Thraupidae). We produced ∼30-fold coverage of the genome with an assembly size of ca. 1.2 Gb, including approximately 16,500 annotated genes. Passerine birds, such as the bananaquit, are commonly infected by avian malarial parasites (Haemosporida), which presumably drive adaptive evolution of immunogenetic loci within the host genome. In the context of our research on the distribution of avian Haemosporida, we specifically characterized immune loci, including toll-like receptor (TLR) and major histocompatibility complex (MHC) genes. Additionally, we identified novel molecular markers in the form of single nucleotide polymorphisms (SNPs), both genome-wide and within identified immune loci. We discovered nine TLR genes and four MHC genes and identified five other TLR- or MHC- associated genes. Genome-wide, over 6 million high-quality SNPs were annotated, including 568 within TLR genes and 102 in MHC genes. This newly described genome and immune characterization expands the knowledge base for avian genomics and phylogenetics and allows for immune genotyping in the bananaquit, providing tools for the investigation of host-parasite coevolution.

The influence of trematode parasite burden on gene expression in a mammalian host

Background

Parasites can profoundly impact their hosts and are responsible for a plethora of debilitating diseases. To identify global changes in host gene expression related to parasite infection, we sequenced, assembled, and annotated the liver transcriptomes of Balb/cj mice infected with the trematode parasite Schistosoma mansoni and compared the results to uninfected mice. We used two different methodologies (i.e. de novo and reference guided) to evaluate the influence of parasite sequences on host transcriptome assembly.

Results

Our results demonstrate that the choice of assembly methodology significantly impacted the proportion of parasitic reads detected from the host library, yet the presence of non-target (xenobiotic) sequences did not create significant structural errors in the assembly. After removing parasite sequences from the mouse transcriptomes, we analyzed host gene expression under different parasite infection levels and observed significant differences in the associated immunologic and metabolic responses based on infection level. In particular, genes associated with T–helper type 1 (Th–1) and T–helper type 2 (Th–2) were up-regulated in infected mice whereas genes related to amino acid and carbohydrate metabolism were down-regulated in infected mice. These changes in gene expression scale with infection status and likely impact the evolutionary fitness of hosts.

Conclusions

Overall, our data indicate that a) infected mice reduce the expression of key metabolic genes in direct proportion to their infection level; b) infected mice similarly increase the expression of key immune genes in response to infection; c) patterns of gene expression correspond to the pathological symptoms of schistosomiasis; and d) identifying and filtering out non-target sequences (xenobiotics) improves differential expression prediction. Our findings identify parasite targets for RNAi or other therapies and provide a better understanding of the pathology and host immune repertoire involved in response to S. mansoni infections.

Electronic supplementary material

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

In silico hybridization enables transcriptomic illumination of the nature and evolution of Myxozoa

Background

The Myxozoa, a group of oligocellular, obligate endoparasites, has long been poorly understood in an evolutionary context. Recent genome-level sequencing techniques such as RNA-seq have generated large amounts of myxozoan sequence data, providing valuable insight into their evolutionary history. However, sequences from host tissue contamination are present in next-generation sequencing reactions of myxozoan tissue, and differentiating between the two has been inadequately addressed. In order to shed light on the genetic underpinnings of myxozoan biology, assembled contigs generated from these studies that derived from the myxozoan must be decoupled from transcripts derived from host tissue and other contamination. This study describes a pipeline for categorization of transcripts asmyxozoan based on similarity searching with known host and parasite sequences, explores the extent to which host contamination is present in previously existing myxozoan datasets, and implements this pipeline on a newly sequenced transcriptome of Myxobolus pendula, a parasite of the common creek chub gill arch.

Methods

The insilico hybridization pipeline uses iterative BLAST searching and database-driven e-value comparison to categorize transcripts as deriving from host, parasite, or other contamination. Functional genetic analysis of M. pendula was conducted using further BLAST searching, Hidden Markov Modeling, and sequence alignment and phylogenetic reconstruction.

Results

Three RNA libraries of encysted M. pendula plasmodia were sequenced and subjected to the method. Nearly half of the final set of contiguous assembly sequences (47.3 %) was identified as putative myxozoan transcripts. Putative contamination was also identified in at least 1/3^rd of previously published myxozoan transcripts. The set of M. pendula transcripts was mined for a range of biologically insightful genes, including taxonomically restricted nematocyst structural proteins and nematocyst proteins identified through mass tandem spectrometry of other cnidarians. Several novel findings emerged, including a fourth myxozoan minicollagen gene, putative myxozoan toxin proteins,and extracellular matrix glycoproteins.

Conclusions

This study serves as a model for the handling of next-generation myxozoan sequence. The need for careful categorization was demonstrated in both previous and new sets of myxozoan sequences. The final set of confidently assigned myxozoan transcripts can be mined for any biologically relevant gene or gene family without spurious misidentification of host contamination as a myxozoan homolog. As exemplified by M. pendula, the repertoire of myxozoan polar capsules may be more complex than previously thought, with an additional minicollagen homolog and putative expression of toxin proteins.

Electronic supplementary material

The online version of this article (doi:10.1186/s12864-015-2039-6) contains supplementary material, which is available to authorized users.

Antarctic krill population genomics: apparent panmixia, but genome complexity and large population size muddy the water

Antarctic krill (Euphausia superba; hereafter krill) are an incredibly abundant pelagic crustacean which has a wide, but patchy, distribution in the Southern Ocean. Several studies have examined the potential for population genetic structuring in krill, but DNA-based analyses have focused on a limited number of markers and have covered only part of their circum-Antarctic range. We used mitochondrial DNA and restriction site-associated DNA sequencing (RAD-seq) to investigate genetic differences between krill from five sites, including two from East Antarctica. Our mtDNA results show no discernible genetic structuring between sites separated by thousands of kilometres, which is consistent with previous studies. Using standard RAD-seq methodology, we obtained over a billion sequences from >140 krill, and thousands of variable nucleotides were identified at hundreds of loci. However, downstream analysis found that markers with sufficient coverage were primarily from multicopy genomic regions. Careful examination of these data highlights the complexity of the RAD-seq approach in organisms with very large genomes. To characterize the multicopy markers, we recorded sequence counts from variable nucleotide sites rather than the derived genotypes; we also examined a small number of manually curated genotypes. Although these analyses effectively fingerprinted individuals, and uncovered a minor laboratory batch effect, no population structuring was observed. Overall, our results are consistent with panmixia of krill throughout their distribution. This result may indicate ongoing gene flow. However, krill's enormous population size creates substantial panmictic inertia, so genetic differentiation may not occur on an ecologically relevant timescale even if demographically separate populations exist.

Molecular tools and bumble bees: revealing hidden details of ecology and evolution in a model system

Bumble bees are a longstanding model system for studies on behaviour, ecology and evolution, due to their well-studied social lifestyle, invaluable role as wild and managed pollinators, and ubiquity and diversity across temperate ecosystems. Yet despite their importance, many aspects of bumble bee biology have remained enigmatic until the rise of the genetic and, more recently, genomic eras. Here, we review and synthesize new insights into the ecology, evolution and behaviour of bumble bees that have been gained using modern genetic and genomic techniques. Special emphasis is placed on four areas of bumble bee biology: the evolution of eusociality in this group, population-level processes, large-scale evolutionary relationships and patterns, and immunity and resistance to pesticides. We close with a prospective on the future of bumble bee genomics research, as this rapidly advancing field has the potential to further revolutionize our understanding of bumble bees, particularly in regard to adaptation and resilience. Worldwide, many bumble bee populations are in decline. As such, throughout the review, connections are drawn between new molecular insights into bumble bees and our understanding of the causal factors involved in their decline. Ongoing and potential applications to bumble bee management and conservation are also included to demonstrate how genetics- and genomics-enabled research aids in the preservation of this threatened group.

Immunomics of the koala (Phascolarctos cinereus)

The study of the koala transcriptome has the potential to advance our understanding of its immunome--immunological reaction of a given host to foreign antigens--and to help combat infectious diseases (e.g., chlamydiosis) that impede ongoing conservation efforts. We used Illumina sequencing of cDNA to characterize genes expressed in two different koala tissues of immunological importance, blood and spleen. We generated nearly 600 million raw sequence reads, and about 285 million of these were subsequently assembled and condensed into ~70,000 subcomponents that represent putative transcripts. We annotated ~16% of these subcomponents and identified those related to infection and the immune response, including Toll-like receptors (TLRs), RIG-I-like receptors (RLRs), major histocompatibility complex (MHC) genes, and koala retrovirus (KoRV). Using phylogenetic analyses, we identified 29 koala genes in these target categories and report their concordance with currently accepted gene groups. By mapping multiple sequencing reads to transcripts, we identified 56 putative SNPs in genes of interest. The distribution of these SNPs indicates that MHC genes (34 SNPs) are more diverse than KoRV (12 SNPs), TLRs (8 SNPs), or RLRs (2 SNPs). Our sequence data also indicate that KoRV sequences are highly expressed in the transcriptome. Our efforts have produced full-length sequences for potentially important immune genes in koala, which should serve as targets for future investigations that aim to conserve koala populations.

Transcriptomic response to heat stress among ecologically divergent populations of redband trout

Background

As ectothermic organisms have evolved to differing aquatic climates, the molecular basis of thermal adaptation is a key area of research. In this study, we tested for differential transcriptional response of ecologically divergent populations of redband trout (Oncorhynchus mykiss gairdneri) that have evolved in desert and montane climates. Each pure strain and their F1 cross were reared in a common garden environment and exposed over four weeks to diel water temperatures that were similar to those experienced in desert climates within the species’ range. Gill tissues were collected from the three strains of fish (desert, montane, F1 crosses) at the peak of heat stress and tested for mRNA expression differences across the transcriptome with RNA-seq.

Results

Strong differences in transcriptomic response to heat stress were observed across strains confirming that fish from desert environments have evolved diverse mechanisms to cope with stressful environments. As expected, a large number of total transcripts (12,814) were differentially expressed in the study (FDR ≤ 0.05) with 2310 transcripts in common for all three strains, but the desert strain had a larger number of unique differentially expressed transcripts (2875) than the montane (1982) or the F1 (2355) strain. Strongly differentiated genes (>4 fold change and FDR ≤ 0.05) were particularly abundant in the desert strain (824 unique contigs) relative to the other two strains (montane = 58; F1 = 192).

Conclusions

This study demonstrated patterns of acclimation (i.e., phenotypic plasticity) within strains and evolutionary adaptation among strains in numerous genes throughout the transcriptome. Key stress response genes such as molecular chaperones (i.e., heat shock proteins) had adaptive patterns of gene expression among strains, but also a much higher number of metabolic and cellular process genes were differentially expressed in the desert strain demonstrating these biological pathways are critical for thermal adaptation to warm aquatic climates. The results of this study further elucidate the molecular basis for thermal adaptation in aquatic ecosystems and extend the potential for identifying genes that may be critical for adaptation to changing climates.

Electronic supplementary material

The online version of this article (doi:10.1186/s12864-015-1246-5) contains supplementary material, which is available to authorized users.

Recent insertion/deletion (reINDEL) mutations: increasing awareness to boost molecular-based research in ecology and evolution

Today, the comparative analysis of DNA molecules mainly uses information inferred from nucleotide substitutions. Insertion/deletion (INDEL) mutations, in contrast, are largely considered uninformative and discarded, due to our lacking knowledge on their evolution. However, including rather than discarding INDELs would be relevant to any research area in ecology and evolution that uses molecular data. As a practical approach to better understanding INDEL evolution in general, we propose the study of recent INDEL (reINDEL) mutations – mutations where both ancestral and derived state are seen in the sample. The precondition for reINDEL identification is knowledge about the pedigree of the individuals sampled. Sound reINDEL knowledge will allow the improved modeling needed for including INDELs in the downstream analysis of molecular data. Both microsatellites, currently still the predominant marker system in the analysis of populations, and sequences generated by next-generation sequencing, a promising and rapidly developing range of technologies, offer the opportunity for reINDEL identification. However, a 2013 sample of animal microsatellite studies contained unexpectedly few reINDELs identified. As most likely explanation, we hypothesize that reINDELs are underreported rather than absent and that this underreporting stems from common reINDEL unawareness. If our hypothesis applies, increased reINDEL awareness should allow gathering data rapidly. We recommend the routine reporting of either the absence or presence of reINDELs together with standardized key information on the nature of mutations when they are detected and the use of the keyword “reINDEL” to increase visibility in both instances of successful and unsuccessful search.

Transcriptomic characterization of the immunogenetic repertoires of heteromyid rodents

Background

When populations evolve under disparate environmental conditions, they experience different selective pressures that shape patterns of sequence evolution and gene expression. These may be manifested in genetic and phenotypic differences such as a diverse immunogenetic repertoire in species from tropical latitudes that have greater and/or different parasite burdens than more temperate species. To test this idea, we compared the transcriptomes of one tropical species (Heteromys desmarestianus) and two species from temperate latitudes (Dipodomys spectabilis and Chaetodipus baileyi) from the Heteromyidae. We did so in a search for positive selection on sequences and/or differential expression, while controlling for phylogenetic history in our choice of species.

Results

We identified 127,812 contigs and annotated 34,878 of these, identifying immune genes associated with interleukins, cytokines, and the production of mast cells. We identified 632 genes that were upregulated in H. desmarestianus (8.7% of genes tested) and 492 (6.7%) that were downregulated. Gene ontology terms including “immune response” were associated with 31 (4.9%) of the 632 upregulated genes. We found preliminary evidence for positive selection on three genes (Palmitoyltransferase ZDHHC5 Ubiquitin-conjugating enzyme E2 N, Krueppel-like factor 10, and Spindle and kinetochore-associated protein 1) along the H. desmarestianus lineage.

Conclusions

Overall our findings pinpoint genes in species from disparate environments that are on different evolutionary trajectories in terms of expression levels and/or nucleotide sequence. Our data indicate there are significant differences in the expression of genes among the spleen transcriptomes of these species and that a number of these differentially expressed genes do not show the same pattern of differential expression in another tissue type. This points to the possibility of expression differences between these species specific to the spleen transcriptome.

Electronic supplementary material

The online version of this article (doi:10.1186/1471-2164-15-929) contains supplementary material, which is available to authorized users.

Models of selection, isolation, and gene flow in speciation

Many marine ecologists aspire to use genetic data to understand how selection and demographic history shape the evolution of diverging populations as they become reproductively isolated species. I propose combining two types of genetic analysis focused on this key early stage of the speciation process to identify the selective agents directly responsible for population divergence. Isolation-with-migration (IM) models can be used to characterize reproductive isolation between populations (low gene flow), while codon models can be used to characterize selection for population differences at the molecular level (especially positive selection for high rates of amino acid substitution). Accessible transcriptome sequencing methods can generate the large quantities of data needed for both types of analysis. I highlight recent examples (including our work on fertilization genes in sea stars) in which this confluence of interest, models, and data has led to taxonomically broad advances in understanding marine speciation at the molecular level. I also highlight new models that incorporate both demography and selection: simulations based on these theoretical advances suggest that polymorphisms shared among individuals (a key source of information in IM models) may lead to false-positive evidence of selection (in codon models), especially during the early stages of population divergence and speciation that are most in need of study. The false-positive problem may be resolved through a combination of model improvements plus experiments that document the phenotypic and fitness effects of specific polymorphisms for which codon models and IM models indicate selection and reproductive isolation (such as genes that mediate sperm-egg compatibility at fertilization).

The genome sequence of a widespread apex predator, the golden eagle (Aquila chrysaetos)

Biologists routinely use molecular markers to identify conservation units, to quantify genetic connectivity, to estimate population sizes, and to identify targets of selection. Many imperiled eagle populations require such efforts and would benefit from enhanced genomic resources. We sequenced, assembled, and annotated the first eagle genome using DNA from a male golden eagle (Aquila chrysaetos) captured in western North America. We constructed genomic libraries that were sequenced using Illumina technology and assembled the high-quality data to a depth of ∼40x coverage. The genome assembly includes 2,552 scaffolds >10 Kb and 415 scaffolds >1.2 Mb. We annotated 16,571 genes that are involved in myriad biological processes, including such disparate traits as beak formation and color vision. We also identified repetitive regions spanning 92 Mb (∼6% of the assembly), including LINES, SINES, LTR-RTs and DNA transposons. The mitochondrial genome encompasses 17,332 bp and is ∼91% identical to the Mountain Hawk-Eagle (Nisaetus nipalensis). Finally, the data reveal that several anonymous microsatellites commonly used for population studies are embedded within protein-coding genes and thus may not have evolved in a neutral fashion. Because the genome sequence includes ∼800,000 novel polymorphisms, markers can now be chosen based on their proximity to functional genes involved in migration, carnivory, and other biological processes.

A pipeline for the de novo assembly of the Themira biloba (Sepsidae: Diptera) transcriptome using a multiple k-mer length approach

BACKGROUND

The Sepsidae family of flies is a model for investigating how sexual selection shapes courtship and sexual dimorphism in a comparative framework. However, like many non-model systems, there are few molecular resources available. Large-scale sequencing and assembly have not been performed in any sepsid, and the lack of a closely related genome makes investigation of gene expression challenging. Our goal was to develop an automated pipeline for de novo transcriptome assembly, and to use that pipeline to assemble and analyze the transcriptome of the sepsid Themira biloba.

RESULTS

Our bioinformatics pipeline uses cloud computing services to assemble and analyze the transcriptome with off-site data management, processing, and backup. It uses a multiple k-mer length approach combined with a second meta-assembly to extend transcripts and recover more bases of transcript sequences than standard single k-mer assembly. We used 454 sequencing to generate 1.48 million reads from cDNA generated from embryo, larva, and pupae of T. biloba and assembled a transcriptome consisting of 24,495 contigs. Annotation identified 16,705 transcripts, including those involved in embryogenesis and limb patterning. We assembled transcriptomes from an additional three non-model organisms to demonstrate that our pipeline assembled a higher-quality transcriptome than single k-mer approaches across multiple species.

CONCLUSIONS

The pipeline we have developed for assembly and analysis increases contig length, recovers unique transcripts, and assembles more base pairs than other methods through the use of a meta-assembly. The T. biloba transcriptome is a critical resource for performing large-scale RNA-Seq investigations of gene expression patterns, and is the first transcriptome sequenced in this Dipteran family.