Improving amphibian genomic resources: a multitissue reference transcriptome of an iconic invader

Whole-mitogenome analysis unveils previously undescribed genetic diversity in cane toads across their invasion trajectory

Invasive species offer insights into rapid adaptation to novel environments. The iconic cane toad (Rhinella marina) is an excellent model for studying rapid adaptation during invasion. Previous research using the mitochondrial NADH dehydrogenase 3 (ND3) gene in Hawai'ian and Australian invasive populations found a single haplotype, indicating an extreme genetic bottleneck following introduction. Nuclear genetic diversity also exhibited reductions across the genome in these two populations. Here, we investigated the mitochondrial genomics of cane toads across this invasion trajectory. We created the first reference mitochondrial genome for this species using long-read sequence data. We combined whole-genome resequencing data of 15 toads with published transcriptomic data of 125 individuals to construct nearly complete mitochondrial genomes from the native (French Guiana) and introduced (Hawai'i and Australia) ranges for population genomic analyses. In agreement with previous investigations of these populations, we identified genetic bottlenecks in both Hawai'ian and Australian introduced populations, alongside evidence of population expansion in the invasive ranges. Although mitochondrial genetic diversity in introduced populations was reduced, our results revealed that it had been underestimated: we identified 45 mitochondrial haplotypes in Hawai'ian and Australian samples, none of which were found in the native range. Additionally, we identified two distinct groups of haplotypes from the native range, separated by a minimum of 110 base pairs (0.6%). These findings enhance our understanding of how invasion has shaped the genetic landscape of this species.

A cane toad (Rhinella marina) N-methyltransferase converts primary indolethylamines to tertiary psychedelic amines

Psychedelic indolethylamines have emerged as potential medicines to treat several psychiatric pathologies. Natural sources of these compounds include 'magic mushrooms' (Psilocybe spp.), plants used to prepare ayahuasca, and toads. The skin and parotid glands of certain toads accumulate a variety of specialized metabolites including toxic guanidine alkaloids, lipophilic alkaloids, poisonous steroids, and hallucinogenic indolethylamines such as DMT, 5-methoxy-DMT, and bufotenin. The occurrence of psychedelics has contributed to the ceremonial use of toads, particularly among Mesoamerican peoples. Yet, the biosynthesis of psychedelic alkaloids has not been elucidated. Herein, we report a novel indolethylamine N-methyltransferase (RmNMT) from cane toad (Rhinella marina). The RmNMT sequence was used to identify a related NMT from the common toad, Bufo bufo. Close homologs from various frog species were inactive, suggesting a role for psychedelic indolethylamine biosynthesis in toads. Enzyme kinetic analyses and comparison with functionally similar enzymes showed that recombinant RmNMT was an effective catalyst and not product inhibited. The substrate promiscuity of RmNMT enabled the bioproduction of a variety of substituted indolethylamines at levels sufficient for purification, pharmacological screening, and metabolic stability assays. Since the therapeutic potential of psychedelics has been linked to activity at serotonergic receptors, we evaluated binding of derivatives at 5-HT1A and 5-HT2A receptors. Primary amines exhibited enhanced affinity at the 5-HT1A receptor compared with tertiary amines. With the exception of 6-substituted derivatives, N,N-dimethylation also protected against catabolism by liver microsomes.

Captivity induces large and population-dependent brain transcriptomic changes in wild-caught cane toads (Rhinella marina)

Gene expression levels are key molecular phenotypes at the interplay between genotype and environment. Mounting evidence suggests that short-term changes in environmental conditions, such as those encountered in captivity, can substantially affect gene expression levels. Yet, the exact magnitude of this effect, how general it is, and whether it results in parallel changes across populations are not well understood. Here, we take advantage of the well-studied cane toad, Rhinella marina, to examine the effect of short-term captivity on brain gene expression levels, and determine whether effects of captivity differ between long-colonized and vanguard populations of the cane toad's Australian invasion range. We compared the transcriptomes of wild-caught toads immediately assayed with those from toads captured from the same populations but maintained in captivity for seven months. We found large differences in gene expression levels between captive and wild-caught toads from the same population, with an over-representation of processes related to behaviour and the response to stress. Captivity had a much larger effect on both gene expression levels and gene expression variability in toads from vanguard populations compared to toads from long-colonized areas, potentially indicating an increased plasticity in toads at the leading edge of the invasion. Overall, our findings indicate that short-term captivity can induce large and population-specific transcriptomic changes, which has significant implications for studies comparing phenotypic traits of wild-caught organisms from different populations that have been held in captivity.

Transcriptomic analysis reveals potential candidate pathways and genes involved in toxin biosynthesis in true toads

Synthesized chemical defenses have broadly evolved across countless taxa and are important in shaping evolutionary and ecological interactions within ecosystems. However, the underlying genomic mechanisms by which these organisms synthesize and utilize their toxins are relatively unknown. Herein, we use comparative transcriptomics to uncover potential toxin synthesizing genes and pathways, as well as interspecific patterns of toxin synthesizing genes across ten species of North American true toads (Bufonidae). Upon assembly and annotation of the ten transcriptomes, we explored patterns of relative gene expression and possible protein-protein interactions across the species to determine what genes and/or pathways may be responsible for toxin synthesis. We also tested our transcriptome dataset for signatures of positive selection to reveal how selection may be acting upon potential toxin producing genes. We assembled high quality transcriptomes of the bufonid parotoid gland, a tissue not often investigated in other bufonid related RNAseq studies. We found several genes involved in metabolic and biosynthetic pathways (e.g. steroid biosynthesis, terpenoid backbone biosynthesis, isoquinoline biosynthesis, glucosinolate biosynthesis) that were functionally enriched and/or relatively expressed across the ten focal species that may be involved in the synthesis of alkaloid and steroid toxins, as well as other small metabolic compounds that cause distastefulness in bufonids. We hope that our study lays a foundation for future studies to explore the genomic underpinnings and specific pathways of toxin synthesis in toads, as well as at the macroevolutionary scale across numerous taxa that produce their own defensive toxins.

Brain transcriptome analysis reveals gene expression differences associated with dispersal behaviour between range-front and range-core populations of invasive cane toads in Australia

Understanding the mechanisms allowing invasive species to adapt to novel environments is a challenge in invasion biology. Many invaders demonstrate rapid evolution of behavioural traits involved in range expansion such as locomotor activity, exploration and risk‐taking. However, the molecular mechanisms that underpin these changes are poorly understood. In 86 years, invasive cane toads (Rhinella marina) in Australia have drastically expanded their geographic range westward from coastal Queensland to Western Australia. During their range expansion, toads have undergone extensive phenotypic changes, particularly in behaviours that enhance the toads’ dispersal ability. Common‐garden experiments have shown that some changes in behavioural traits related to dispersal are heritable. At the molecular level, it is currently unknown whether these changes in dispersal‐related behaviour are underlain by small or large differences in gene expression, nor is known the biological function of genes showing differential expression. Here, we used RNA‐seq to gain a better understanding of the molecular mechanisms underlying dispersal‐related behavioural changes. We compared the brain transcriptomes of toads from the Hawai'ian source population, as well as three distinct populations from across the Australian invasive range. We found markedly different gene expression profiles between the source population and Australian toads. By contrast, toads from across the Australian invasive range had very similar transcriptomic profiles. Yet, key genes with functions putatively related to dispersal behaviour showed differential expression between populations located at each end of the invasive range. These genes could play an important role in the behavioural changes characteristic of range expansion in Australian cane toads.

Identification of Na+/K+-ATPase α/β isoforms in Rhinella marina tissues by RNAseq and a molecular docking approach at the protein level to evaluate α isoform affinities for bufadienolides

Na⁺/K⁺-ATPase (NKA) function is inhibited by Bufadienolides (BD), a group of cardiotonic steroids (CTS) primarily produced by anurans of the Bufonidae family, such as Rhinella marina. This study characterized the presence of α and β NKA subunit isoforms in R. marina via RNAseq in four tissues: oocytes, skin, heart, and skeletal muscle. Transcripts encoding three α-like isoforms (α₁, α₂, α₃) and three β-like isoforms (β₁, β₂, β₄) were identified. The amino acid sequence of α₁-like isoform shared 99.4% identity with the α₁ isoform previously published for R. marina. Sequences for α₂, α₃, and β₄ from R. marina were previously unavailable. The first extracellular loop in the α₂-like isoform in R. marina showed similar substitutions to those found in their susceptible homologues in other taxa (L/Q111T and S119T); in contrast, this same loop in α₃-like isoform showed similar substitutions (Q111L and G120R) to those reported for toad-eating animals such as snakes, which suggests relatively lower affinity for CTS. Docking results showed that all three α-like isoforms identified in R. marina transcriptomes have low affinity to CTS compared to the susceptible α₁ isoform of Sus scrofa (pig), with α₁-like isoform being the most resistant. The tissue-specific RNAseq results showed the following expression of NKA α-like and β-like subunit isoforms: Oocytes expressed α₁ and β₁; skin α₁, β₁, and low levels of β₂; heart α₁, α₃, and β₁; skeletal muscle α₁, β₄, with low levels of α₂, α₃, and β₁. R. marina could be used as an important model for future structural, functional and pharmacological studies of NKA and its isoforms.

The rise and fall of globins in the amphibia

The globin gene repertoire of gnathostome vertebrates is dictated by differential retention and loss of nine paralogous genes: androglobin, neuroglobin, globin X, cytoglobin, globin Y, myoglobin, globin E, and the α- and β-globins. We report the globin gene repertoire of three orders of modern amphibians: Anura, Caudata, and Gymnophiona. Combining phylogenetic and conserved synteny analysis, we show that myoglobin and globin E were lost only in the Batrachia clade, but retained in Gymnophiona. The major amphibian groups also retained different paralogous copies of globin X. None of the amphibian presented α^D-globin gene. Nevertheless, two clades of β-globins are present in all amphibians, indicating that the amphibian ancestor possessed two paralogous proto β-globins. We also show that orthologs of the gene coding for the monomeric hemoglobin found in the heart of Rana catesbeiana are present in Neobatrachia and Pelobatoidea species we analyzed. We suggest that these genes might perform myoglobin- and globin E-related functions. We conclude that the repertoire of globin genes in amphibians is dictated by both retention and loss of the paralogous genes cited above and the rise of a new globin gene through co-option of an α-globin, possibly facilitated by a prior event of transposition.

The physiology of movement

Movement, from foraging to migration, is known to be under the influence of the environment. The translation of environmental cues to individual movement decision making is determined by an individual's internal state and anticipated to balance costs and benefits. General body condition, metabolic and hormonal physiology mechanistically underpin this internal state. These physiological determinants are tightly, and often genetically linked with each other and hence central to a mechanistic understanding of movement. We here synthesise the available evidence of the physiological drivers and signatures of movement and review (1) how physiological state as measured in its most coarse way by body condition correlates with movement decisions during foraging, migration and dispersal, (2) how hormonal changes underlie changes in these movement strategies and (3) how these can be linked to molecular pathways. We reveale that a high body condition facilitates the efficiency of routine foraging, dispersal and migration. Dispersal decision making is, however, in some cases stimulated by a decreased individual condition. Many of the biotic and abiotic stressors that induce movement initiate a physiological cascade in vertebrates through the production of stress hormones. Movement is therefore associated with hormone levels in vertebrates but also insects, often in interaction with factors related to body or social condition. The underlying molecular and physiological mechanisms are currently studied in few model species, and show -in congruence with our insights on the role of body condition- a central role of energy metabolism during glycolysis, and the coupling with timing processes during migration. Molecular insights into the physiological basis of movement remain, however, highly refractory. We finalise this review with a critical reflection on the importance of these physiological feedbacks for a better mechanistic understanding of movement and its effects on ecological dynamics at all levels of biological organization.

Increased Adaptive Variation Despite Reduced Overall Genetic Diversity in a Rapidly Adapting Invader

Invasive species often evolve rapidly following introduction despite genetic bottlenecks that may result from small numbers of founders; however, some invasions may not fit this “genetic paradox”. The invasive cane toad (Rhinella marina) displays high phenotypic variation across its introduced Australian range. Here, we used three genome-wide datasets to characterize their population structure and genetic diversity. We found that toads form three genetic clusters: 1) native range toads, 2) toads from the source population in Hawaii and long-established areas near introduction sites in Australia, and 3) toads from more recently established northern Australian sites. Although we find an overall reduction in genetic diversity following introduction, we do not see this reduction in loci putatively under selection, suggesting that genetic diversity may have been maintained at ecologically relevant traits, or that mutation rates were high enough to maintain adaptive potential. Nonetheless, toads encounter novel environmental challenges in Australia, and the transition between genetic clusters occurs at a point along the invasion transect where temperature rises and rainfall decreases. We identify environmentally associated loci known to be involved in resistance to heat and dehydration. This study highlights that natural selection occurs rapidly and plays a vital role in shaping the structure of invasive populations.

Ecotoxicity of nanomaterials in amphibians: A critical review

Nanomaterials (NMs) have been used in a growing number of commercial products, and their rapid expansion could lead to their release into the aquatic environments. However, there is limited knowledge about the impact of NMs in the biota, especially the amphibians. The present study revised the historical use of amphibian species as a model system for nanoecotoxicological studies and summarized the data available in the scientific literature about the genotoxic, mutagenic, histopathological, embryotoxic and reproductive effects of NMs in different groups of amphibians. The interaction, bioaccumulation, mode of action (MoA) and ecotoxicity of NMs on amphibians were also revised. The nanoecotoxicological studies were conducted with 11 amphibian species, being eight species of the order Anura and three species of the order Caudata. Xenopus laevis was the most studied species. The studies were conducted mainly with inorganic NMs (72%) compared to organic ones. The nanoecotoxicity depends on NM behavior and transformation in the environment, as well as the developmental stages of amphibians. The known effects of NMs in amphibians were mainly reported with reactive oxygen species (ROS) production, oxidative stress, and genotoxic effects. Results emphasize the need for further studies testing the ecotoxicity of different NMs, concentrations and exposure periods at environmentally relevant approaches. Furthermore, standard protocols for nanoecotoxicological tests using amphibians are required. Revised data showed that amphibians are suitable organisms to assess the environmental impact of NMs and indicated significant research gaps concerning the ecotoxicity of NMs on freshwater ecosystems and recommendations for future researches.

De Novo Assembly and Annotation of the Larval Transcriptome of Two Spadefoot Toads Widely Divergent in Developmental Rate

Amphibians are highly vulnerable and diverse vertebrates for which we still have modest genomic resources. Amphibian larvae are key components of continental wetlands, where they have strong influences on energy fluxes, nutrient cycling, and community structure. Amphibian larvae are highly sensitive to environmental conditions and can often alter their physiology, behavior and even morphology in response to the local conditions experienced, although we still know relatively little about the transcriptomic changes that enable such plasticity. Here we contribute the larval transcriptomes of two spadefoot toad species with divergent developmental rates and degree of developmental plasticity in response to pond drying.

Genomic Takeover by Transposable Elements in the Strawberry Poison Frog

We sequenced the genome of the strawberry poison frog, Oophaga pumilio, at a depth of 127.5× using variable insert size libraries. The total genome size is estimated to be 6.76 Gb, of which 4.76 Gb are from high copy number repetitive elements with low differentiation across copies. These repeats encompass DNA transposons, RNA transposons, and LTR retrotransposons, including at least 0.4 and 1.0 Gb of Mariner/Tc1 and Gypsy elements, respectively. Expression data indicate high levels of gypsy and Mariner/Tc1 expression in ova of O. pumilio compared with Xenopus laevis. We further observe phylogenetic evidence for horizontal transfer (HT) of Mariner elements, possibly between fish and frogs. The elements affected by HT are present in high copy number and are highly expressed, suggesting ongoing proliferation after HT. Our results suggest that the large amphibian genome sizes, at least partially, can be explained by a process of repeated invasion of new transposable elements that are not yet suppressed in the germline. We also find changes in the spliceosome that we hypothesize are related to permissiveness of O. pumilio to increases in intron length due to transposon proliferation. Finally, we identify the complement of ion channels in the first genomic sequenced poison frog and discuss its relation to the evolution of autoresistance to toxins sequestered in the skin.

Immune and environment-driven gene expression during invasion: An eco-immunological application of RNA-Seq

Host-pathogen associations change rapidly during a biological invasion and are predicted to impose strong selection on immune function. It has been proposed that the invader may experience an abrupt reduction in pathogen-mediated selection ("enemy release"), thereby favoring decreased investment into "costly" immune responses. Across plants and animals, there is mixed support for this prediction. Pathogens are not the only form of selection imposed on invaders; differences in abiotic environmental conditions between native and introduced ranges are also expected to drive rapid evolution. Here, we use RNA-Seq to assess the expression patterns of immune and environmentally associated genes in the cane toad (Rhinella marina) across its invasive Australian range. Transcripts encoding mediators of costly immune responses (inflammation, cytotoxicity) showed a curvilinear relationship with invasion history, with highest expression in toads from oldest and newest colonized areas. This pattern is surprising given theoretical expectations of density dynamics in invasive species and may be because density influences both intraspecific competition and parasite transmission, generating conflicting effects on the strength of immune responses. Alternatively, this expression pattern may be the result of other evolutionary forces, such as spatial sorting and genetic drift, working simultaneously with natural selection. Our findings do not support predictions about immune function based on the enemy release hypothesis and suggest instead that the effects of enemy release are difficult to isolate in wild populations, especially in the absence of information regarding parasite and pathogen infection. Additionally, expression patterns of genes underlying putatively environmentally associated traits are consistent with previous genetic studies, providing further support that Australian cane toads have adapted to novel abiotic challenges.

RNA Sequencing of Peripheral Blood Revealed that the Neurotropic TRK Receptor Signaling Pathway Shows Apparent Correlation in Recovery Following Spinal Cord Injury at Small Cohort

Spinal cord injury (SCI) can be lethal; however, the precise mechanisms underlying healing are unclear, limiting the development of effective therapies. In this study, the molecular mechanisms involved in SCI were investigated. Clinical peripheral blood samples from normal individuals and patients with incomplete SCI (ISCI) and complete SCI (CSCI) were analyzed by RNA-Seq. The expression levels of EPHA4, CDK16, BAD, MAP2 Normal 2, EGR, and RHOB differed significantly between the SCI group and normal individuals, and these results were verified by q-PCR. A gene ontology (GO) enrichment analysis showed that differentially expressed genes were mostly enriched for the neurotrophin TRK receptor signaling pathway. We verified the expression of neurotrophic factors and found that they were all expressed most highly in the SCI group. The results of this study demonstrate that neurotrophic factors are highly expressed after SCI and the neurotrophin TRK receptor signaling pathway may be involved in the initiation of nerve system regeneration.

Improving amphibian genomic resources: a multitissue reference transcriptome of an iconic invader

Background

Cane toads (Rhinella marina) are an iconic invasive species introduced to 4 continents and well utilized for studies of rapid evolution in introduced environments. Despite the long introduction history of this species, its profound ecological impacts, and its utility for demonstrating evolutionary principles, genetic information is sparse. Here we produce a de novo transcriptome spanning multiple tissues and life stages to enable investigation of the genetic basis of previously identified rapid phenotypic change over the introduced range.

Findings

Using approximately 1.9 billion reads from developing tadpoles and 6 adult tissue-specific cDNA libraries, as well as a transcriptome assembly pipeline encompassing 100 separate de novo assemblies, we constructed 62 202 transcripts, of which we functionally annotated ∼50%. Our transcriptome assembly exhibits 90% full-length completeness of the Benchmarking Universal Single-Copy Orthologs data set. Robust assembly metrics and comparisons with several available anuran transcriptomes and genomes indicate that our cane toad assembly is one of the most complete anuran genomic resources available.

Conclusions

This comprehensive anuran transcriptome will provide a valuable resource for investigation of genes under selection during invasion in cane toads, but will also greatly expand our general knowledge of anuran genomes, which are underrepresented in the literature. The data set is publically available in NCBI and GigaDB to serve as a resource for other researchers.

Draft genome assembly of the invasive cane toad, Rhinella marina

Background

The cane toad (Rhinella marina formerly Bufo marinus) is a species native to Central and South America that has spread across many regions of the globe. Cane toads are known for their rapid adaptation and deleterious impacts on native fauna in invaded regions. However, despite an iconic status, there are major gaps in our understanding of cane toad genetics. The availability of a genome would help to close these gaps and accelerate cane toad research.

Findings

We report a draft genome assembly for R. marina, the first of its kind for the Bufonidae family. We used a combination of long-read Pacific Biosciences RS II and short-read Illumina HiSeq X sequencing to generate 359.5 Gb of raw sequence data. The final hybrid assembly of 31,392 scaffolds was 2.55 Gb in length with a scaffold N50 of 168 kb. BUSCO analysis revealed that the assembly included full length or partial fragments of 90.6% of tetrapod universal single-copy orthologs (n = 3950), illustrating that the gene-containing regions have been well assembled. Annotation predicted 25,846 protein coding genes with similarity to known proteins in Swiss-Prot. Repeat sequences were estimated to account for 63.9% of the assembly.

Conclusions

The R. marina draft genome assembly will be an invaluable resource that can be used to further probe the biology of this invasive species. Future analysis of the genome will provide insights into cane toad evolution and enrich our understanding of their interplay with the ecosystem at large.

Benchmarking differential expression analysis tools for RNA-Seq: normalization-based vs. log-ratio transformation-based methods

BACKGROUND

Count data generated by next-generation sequencing assays do not measure absolute transcript abundances. Instead, the data are constrained to an arbitrary "library size" by the sequencing depth of the assay, and typically must be normalized prior to statistical analysis. The constrained nature of these data means one could alternatively use a log-ratio transformation in lieu of normalization, as often done when testing for differential abundance (DA) of operational taxonomic units (OTUs) in 16S rRNA data. Therefore, we benchmark how well the ALDEx2 package, a transformation-based DA tool, detects differential expression in high-throughput RNA-sequencing data (RNA-Seq), compared to conventional RNA-Seq methods such as edgeR and DESeq2.

RESULTS

To evaluate the performance of log-ratio transformation-based tools, we apply the ALDEx2 package to two simulated, and two real, RNA-Seq data sets. One of the latter was previously used to benchmark dozens of conventional RNA-Seq differential expression methods, enabling us to directly compare transformation-based approaches. We show that ALDEx2, widely used in meta-genomics research, identifies differentially expressed genes (and transcripts) from RNA-Seq data with high precision and, given sufficient sample sizes, high recall too (regardless of the alignment and quantification procedure used). Although we show that the choice in log-ratio transformation can affect performance, ALDEx2 has high precision (i.e., few false positives) across all transformations. Finally, we present a novel, iterative log-ratio transformation (now implemented in ALDEx2) that further improves performance in simulations.

CONCLUSIONS

Our results suggest that log-ratio transformation-based methods can work to measure differential expression from RNA-Seq data, provided that certain assumptions are met. Moreover, these methods have very high precision (i.e., few false positives) in simulations and perform well on real data too. With previously demonstrated applicability to 16S rRNA data, ALDEx2 can thus serve as a single tool for data from multiple sequencing modalities.