Skin transcriptional profiles in Oophaga poison frogs

Abstract

Aposematic organisms advertise their defensive toxins to predators using a variety of warning signals, including bright coloration. While most Neotropical poison frogs (Dendrobatidae) rely on crypsis to avoid predators, Oophaga poison frogs from South America advertise their chemical defenses, a complex mix of diet-derived alkaloids, by using conspicuous hues. The present study aimed to characterize the skin transcriptomic profiles of South American Oophaga poison frogs. Our analyses showed very similar transcriptomic profiles for these closely related species in terms of functional annotation and relative abundance of gene ontology terms expressed. Analyses of expression profiles of Oophaga and available skin transcriptomes of cryptic anurans allowed us to propose initial hypotheses for the active sequestration of alkaloid-based chemical defenses and to highlight some genes that may be potentially involved in resistance mechanisms to avoid self-intoxication and skin coloration. In doing so, we provide an important molecular resource for the study of warning signals that will facilitate the assembly and annotation of future poison frog genomes.

Keywords:
Dendrobatids; RNA sequencing; transcriptomes; gene ontology; candidate genes

Introduction

Aposematic organisms advertise their defensive toxins to predators using a variety of warning signals (Harvey et al., 1982; Brown, 2013; Cummings and Crothers, 2013). Among aposematic species, the poison frogs (Dendrobatidae) from the tropical rain forests of Central and South America represent one of the most spectacular examples of warning coloration. While the majority of dendrobatids rely on crypsis to avoid predators, some members of this family are both brightly colored and chemically defended (Grant et al., 2006). From a genetic point of view, these aposematic defenses can be defined as a complex phenotype resulting from the integration (i.e. covariation) of different genetic elements related to conspicuousness, bold behavior, unpalatability, diet specialization, etc. While a long history of research has been devoted to understanding the genetics of warning coloration in arthropods, particularly in Heliconius butterflies (Ramos and Freitas, 1999; Langham and Benkman, 2004), the molecular underpinnings of aposematism in vertebrates, particularly the mechanisms whereby individuals become toxic (or distasteful) remain mostly unknown. In this study, our primary aim was to characterize the skin transcriptomic profiles of Oophaga species as a first attempt to shed light on the molecular genetic bases of aposematic components in poison frogs: the ability to sequester alkaloid-based chemical defenses, warning coloration, and the resistance mechanisms to avoid self-intoxication.

The frogs of this complex inhabit the lowland Pacific rainforests of the Colombian and Ecuadorian Chocó. Previous molecular data from nuclear and mitochondrial markers showed a similar genetic background among lineages in contrast to an extraordinary diversity of morphotypes (Medina et al., 2013). Individuals from different allopatric lineages can be relatively homogeneous, striped, or spotted, and their colors range from bright red, to orange and yellow (Figure 1) (Posso-Terranova and Andrés, 2016). These polymorphic coloration patterns serve as a warning signal of their chemical defenses (Maan and Cummings, 2012), a complex mix of diet-derived alkaloids secreted by the dermal glands (Saporito et al., 2006; Saporito et al., 2012). Although these chemicals may involve metabolism and transport through other tissues, the final accumulation of toxins as well as the production of pigmentary cells is performed in the skin tissue of these organisms even during the adult stage (Bagnara, 1982; Daly, 1995; Wolnicka-Glubisz et al., 2012). Thus, we expected that some of the genes, pathways, and/or gene networks potentially associated with coloration, alkaloid metabolism, transport and storage, could be actively expressed in the skin tissue of these organisms.

Figure 1
Hierarchical clustering (bootstrap=999) and heatmap contrasting the differentially expressed unigenes from skin tissue between Oophaga and cryptic species. Color within the heatmap represent the log2-fold change values. Dark blue indicates upregulation while a lighter coloration indicates downregulation. Rows of the heatmap correspond to the differentially expressed unigenes (n=1,931) detected in this study.

Although the lineages of Oophaga studied here display a wide variety of warning signals (i.e., coloration and patterns), most of them share a black background coloration (Figure 1). In adult frogs, this melanistic coloration is controlled by the lack of iridophores and xanthophores in the dermal chromatophore unit, as well as the size and dispersion pattern of the melanin-containing organelles (melanosomes) within the dermal melanophores (Posso-Terranova and Andrés, 2017). Thus, we also hypothesized that in adult individuals, genes involved in the amount, size and distribution of the melanosomes should be differentially expressed in the skin. Similarly, because the secreted alkaloids of chemically-defended dendrobatids are known to disrupt the normal ion-channel activity and to alter the neurotransmitter-receptor binding capacity of nerve and muscle cells (Daly et al., 1978; Daly et al., 2002; Daly et al., 2003), we hypothesized that the fraction of genes differentially expressed in the skin of poison frogs (when compared to non-aposematic species) may also be enriched in genes related to alkaloid sequestration, transport and autotoxicity avoidance.

To test these hypotheses and to provide an overall transcriptomic characterization, we first assembled the skin transcriptomes of different Harlequin poison frogs lineages: O. histrionica (Berthold, 1846), O. lehmanni (Myers and Daly, 1976), O. sylvatica, (Funkhouser, 1956), O. anchicayensis and O. solanensis (Posso-Terranova and Andrés, 2018). Then, we compared the transcription profiles of these aposematic species with that of a diverse group of cryptic anurans. Functional annotation of unigenes (i.e., assembled transcripts that includes all isoforms from a unique gene) showing a phylogenetic signal (i.e., highly represented and differentially expressed in Oophaga poison frogs) allowed us to propose plausible metabolic pathways (and candidate genes) related with toxicity and color. Overall, this study provides an important molecular resource for the study of aposematism in poison frogs and will facilitate the future assembling and annotation of the complex dendrobatid genomes (Rogers et al., 2018).

Material and Methods

Library preparation and sequencing

Due the conservation status of these frogs (CR, EN, VU; IUCN) and the restrictions imposed by the Ethics Committee, we were allowed to sample and euthanize a very limited number of individuals from natural populations. Seven samples (O. lehmanni n=2; O. anchicayensis, n=2; O. histrionica, O. sylvatica, and O. solanensis, n=1; each) were euthanized in the field with benzocaine gel at 5% by direct application into the mouth cavity. Individual skin samples were deposited in RNAlater® (Life Technologies, Carlbad, CA) for transportation to the laboratory and then, stored at −80 °C until further processing. RNA extractions were performed using Tri-Reagent (Molecular Research Center, Cincinnati, OH). After quality check (Agilent Bioanalyzer Agilent Technologies, Wilmington, DE) samples were used to generate RNA-seq libraries using the Illumina Truseq RNA Sample Prep protocol (Illumina, San Diego, CA). Libraries were cleaned using AMPure XP and sequenced on a single Illumina HiSeq2000 lane (TruSeq SBS v. 3) as follows: four libraries (O. anchicayensis, O. histrionica, O. lehmanni and O. sylvatica) were sequenced using single-end (150-bp) while the remaining three libraries (O. anchicayensis, O. solanensis and O. lehmanni) were run in a single lane of a paired-end module (100-bp, × 2). All animal procedures were approved by the Ethics Committee of Universidad Nacional de Colombia (Acta No.03, July 22^nd, 2015) and were conducted based on the NIH Guide for the Principles of Animal Care. Sampling was conducted according to the research permits granted by Autoridad Nacional de Licencias Ambientales ANLA, Resolución 0255 del 12 de Marzo de 2014 and Resolución 1482 del 20 de Noviembre de 2015. Sequencing was carried out by the Cornell University's BioResource Center.

Transcriptome assemblies and functional annotation

Initial read quality trimming, filtering and removal of adapters was performed using FLEXBAR v.2.5 (Dodt et al., 2012). All retained reads were ≥ 50bp with an average quality of ≥ 30 and less than two uncalled (ambiguous) bases (Table S1). For each library, a de novo transcriptome was assembled using TRINITY v.2.1.0 (Haas et al., 2013) with default settings (kmer=25, minimum contig length=48) and keeping only the longest transcript per cluster for subsequent analyses. Varying these parameters (i.e. default settings) did not result in assemblies with longer N50 values. To account for differences in the quality and sequencing strategy (single vs. paired-end) among individual transcriptomes, we constructed a composite de novo reference transcriptome. Briefly, we followed (Gould et al., 2015) and combined 20% randomly selected reads (>33 million in total) from each of the four single-end libraries (O. anchicayensis, O. histrionica, O. lehmanni and O. sylvatica). Then, a TRINITY assembly was performed using a minimum contig length of 350 bp. To filter out highly similar contigs that may potentially represent alternatively spliced transcripts, we implemented the error correction module of iAssembler v1.3.2 (Zheng et al., 2011) with default parameters (maximum length of end clips=30 bp, minimum overlap length=40 bp, minimum percent sequence identity=95%).

Gene annotation was conducted using a sequential BLASTX search to both, the available Xenopus transcriptomes and the NCBI non-redundant (nr) database. The composite transcripts were first compared with that of the Xenopus databases retaining annotations with E-values ≤ 10^−5. Unannotated contigs were then submitted for BLASTX to the nr protein database for possible identification. Gene ontology (GO) annotation and term mapping was performed using Blast2Go v 5.2.0 with default significance cutoffs (Conesa et al., 2005).

To estimate the completeness of each transcriptome, we implemented the expected gene content of Benchmarking Universal Single-Copy Orthologs (BUSCO v3) (Waterhouse et al., 2018) as implemented in gVolante v1.1.0 (Nishimura et al., 2017). This widely used metrics for transcriptome assembly (Holding et al., 2018) assesses the quality of any given transcriptome by estimating the completeness of core vertebrate genes predicted to be ubiquitous in eukaryotes. The pre-processing of the reads and all BLASTX analyses were run in the Bugaboo Dell Xeon cluster of the western Canada's WestGrid computer facilities (www.westgrid.ca).

Transcriptome profiles

Analyses on the distribution of reads and the relative abundance of different transcripts within each individual sample (i.e., lineage) were carried out by mapping RNA-seq reads back to the composite de novo assembly. In order to identify highly represented unigenes, we first implemented eXpress (probabilistic assignment of ambiguously mapping sequenced fragments) (Roberts and Pachter, 2013) to estimate the effective number of reads (ER) that mapped to the contigs in the reference transcriptome after adjusting for read number and length biases. Then, to compare the proportion of reads that mapped to a transcript in each Oophaga RNA-library (n=7), we estimated the number of transcripts per million (TPM), a measure of RNA abundance that allows the comparison between samples (sum of all TPMs in each sample are the same) (Wagner et al., 2012). To visually explore for highly represented transcripts, we performed a Principal Component Analysis (PCA) of the reference contigs dataset and using each RNA-library TPMs values as independent variables (n=7). Then, we constructed dispersion plots of the ER values and percentile plots of the TPM distribution to select the overall top expressed transcripts (i.e., 2%) for their characterization.

To further characterize the skin transcriptome profile of Oophaga poison frogs, we compared the transcription levels of annotated unigenes in this clade with those observed in several cryptic (i.e. non-aposematic) anurans, including a toad (Bufo gargarizans), and five frog species (Pelophylax nigromaculatus, Polypedates megacephalus, Rana catesbeiana, Rana sylvatica, and Xenophrys sangzhiensis, Table S2) (Huang et al., 2016; Eskew et al., 2018). RNA-seq data from these species was obtained from the Sequence Read Archive platform (SRA) and represented the only skin-based RNA-seq available to us for comparison purposes. Although transcriptomic datasets are available for phylogenetically closer and cryptic species (i.e Colostethus) (Santos et al., 2018), these datasets were not included in our analysis because they are not derived from skin tissue and may bias our results. Read quality, contamination screening, duplicate removal, and quality trimming steps were performed in FLEXBAR v.2.5 as previously described. The filtered raw reads from individual libraries (Oophaga, n=7; cryptic species n=6) were mapped to our composite de novo reference transcriptome using Bowtie2 v2.3.4.1 (Langmead and Salzberg, 2012). To accommodate for sequence divergence among the taxa included in our study (75–98 similarity), we followed (Nicholls et al., 2017) and ran the alignment algorithm allowing for soft clipping (—local) and with a maximum penalty value of three (—mp 3). Then, we implemented eXpress (probabilistic assignment of ambiguously mapping sequenced fragments) (Roberts and Pachter, 2013) to estimate the effective number of reads (ER) that mapped to the contigs in the reference transcriptome after adjusting for read number and length biases. Only unigenes with average of >10 ER/library were included in differential expression analyses (Nanoth Vellichirammal et al., 2014).

To identify genes showing a phylogenetic signal in their expression levels, we estimated fold-changes in expression levels between aposematic (i.e., Oophaga) and cryptic species. To do so, we followed (Schurch et al., 2016) and implemented the R package DEseq (Anders and Huber, 2010) to select unigenes with an adjusted p-value of <0.05. To visualize patterns of expression, we constructed a multidimensional scaling plot (Euclidean distances) of the ER data using Past v. 3.18 (Hammer et al., 2001). A heatmap of the expression differences between Oophaga and cryptic anuran lineages (n_genes=1,931) was generated using the ‘heatmap.2’ function of the ‘gplots’ package, with Euclidean distances and complete linkage for clustering (bootstrap=999) (Warnes et al., 2016). BLAST analysis and GO term mapping for 1,931 genes included in this analysis was performed as described above.

Data availability

The datasets generated and/or analyzed during the current study will be available at Dryad and/or GenBank before peer reviewed publication.

Results

Transcriptome assemblies and functional annotation

Illumina sequencing produced an average number of reads per sample of 118.3 million for paired-end and 37.7 million for single-end libraries. After a stringent read trimming involving the removal of low-quality sequences, duplicated reads and reads containing adapter/primer sequences, we retained over 81% of the initial sequencing data. Paired-end libraries produced a fairly consistent number of reads (Table S1).

After the removal/merging of highly similar contigs (paired-end:7,5% – 9,6%; single-end: 6.0% – 13.4%), the 150-bp based transcriptomes recovered a large number of contigs ranging from 35,287 (O. anchicayensis) to 107,381 (O. sylvatica). Paired-reads libraries (100-bp; 2X) generated transcriptomes with a relatively similar number of contigs (40,398 for O. solanensis and 60,494 for O. lehmanni). N50 values and average transcript length (AL) were lower for paired-end libraries (N50= 498-561bp; AL= 441.78 ± 379 – 475.93 ± 436bp) than for assemblies produced with single-end libraries (N50=667-1579bp; AL=538.9 ± 1001bp – 668.6 ± 819 bp) (Table 1).

Table 1
Summary statistics and quality assessment estimators of the eight de novo assembled transcriptomes generated in this study (Oophaga species n=7 transcriptomes; composite reference n=1).

	Paired-end libraries			Single-end libraries				Reference transcriptome
	O. solanensis	O. anchicayensis	O. lehmanni	O. sylvatica	O. lehmanni	O. anchicayensis	O. histrionica
File size (Mbytes)	22.23	32.06	35.7	93.9	77.3	23.16	85.97	31.6
# contigs	43,690	63,520	66,937	123,916	117,082	37,552	111,823	31,498
# contigs after filtration (removal/merging of highly similar contigs)	40,398	57,516	60,494	107,381	104,523	35,287	96,892	31,498
# bases	19,853,219	28,061,593	31,856,993	81,266,530	66,367,430	20,238,091	74,768,159	31,968,628
# bases after filtration	17,164,733	23,756,585	26,287,312	58,224,462	50,845,151	15,278,290	53,582,406	31,968,628
Merged contigs	7.5%	9.5%	9.6%	13.3%	10.7%	6.0%	13.4%	0.0%
Average transcript length (AL)	454.41 ± 405	441.78 ± 379	475.93 ± 436	655.8 ± 799	566.8 ± 648	538.9 ± 1579	668.6 ± 819	1,014.91 ± 948.43
Maximum length	6,115	7,915	8,400	15,102	12,906	54,265	15,136	16,041
N50	528	498	561	1079	809	667	1,101	1,316
Contigs with BLAST hits	46.5%	47.7%	50.8%	42.1%	41.6%	53.4%	44.5%	62.6%

Comparative BLAST analysis indicated that our assembled transcriptomes recovered a significant proportion of the X. laevis and X. tropicalis reference transcripts (50,592 and 41,042 unigenes respectively). Paired-end transcriptomes showed a lower number of significant BLAST hits (18,263 – 29,670) than single-end libraries (18,526 – 40,651) (Figure S1). In general, individual de-novo transcriptomes were of good quality and recovered 51-79% of the core eukaryotic universal single-copy orthologs (Waterhouse et al., 2018) (Table S3). The final de novo composite skin transcriptome reconstructed from all the RNA reads of the Oophaga species yielded 31,498 contigs greater than 350 bp. The overall assembly incorporated 94% of all initial reads and the level of fragmentation was low with half of all base pairs clustered into contigs of 1,316 bp in length or greater. The maximum contig length was 16,041 bp and the AL was 1014.91 ± 948.43 (Figure S2A). Nucleotide-based BLAST analyses (BLASTX) revealed that ~63% of the contigs (n=19,732) show significant similarities with either annotated gene products and/or known protein domains (E-value ≤ 10-5) (Figure S2B) and only a small fraction of unigenes (3.4%) showed significant homology to the same annotated transcript. The percentage of transcripts in our reference assembly with the highest BLAST scores were found with X. tropicalis (74%), X. laevis (2%) and the green sea turtle Chelonia mydas (0.5%). The composite skin transcriptome was of higher quality when compared with individual assemblies, having 88.54% of the ortholog count of BUSCO (Waterhouse et al., 2018) (Table S3).

Gene Ontology (GO) assignments were used to classify the functions of the predicted unigenes based on contigs with significant BLASTX (E-value ≤ 10-5). Based on GO level II, unigenes were assigned to 27 biological processes (BP), 22 cell components (CC) and 15 molecular functions (MF) (Figure S2C). Some unigenes were associated with multiple GO annotations because a single sequence may be annotated in any or all categories, giving more GO annotations than sequences annotated (Xie et al., 2002). Within BP, ~47% of the annotations were assigned to basic cellular, metabolic processes and biological regulation. The remaining unigenes were involved in a broad range of BP such as response to stimulus (7%), response to stress (7%), localization (6%), developmental process (5%), signal transduction (4%), biogenesis (4%), immune response (2%), reproductive process (1.7%) and cellular adhesion (1%). Within the CC category, other than the essential cell constituents (44%), the membrane components were highly represented in the transcriptome (19%). Within MF, most of the unigenes were assigned to binding and catalytic activities (72%) followed by transporter activity (9%), molecular transducers (7%) and molecular function regulators (4%) (Figure S2C).

Transcriptome profiles

Dispersion plots of ER values in Oophaga RNA libraries showed a common pattern of over-represented unigenes across RNA-libraries. That is, the same group of unigenes showed similar patterns of over-representation no matter RNA library's origin (i.e. different Oophaga species) or composition (paired vs. single-read libraries) (Figure S3). The TPM represents a measure of RNA abundance and hence, it provides a general overview of gene expression levels in a particular sample. Interestingly, TPM distribution plots (percentiles) showed a distinct spike in expression levels towards percentile 98% for all our RNA-libraries (Figure S4A). After the conservative filtration of the 2% top expressed unigenes, we selected a subset of 1,437 unigenes with particular higher levels of expression (higher TPM values and PCA outliers, Figure S4B) in Oophaga skin tissue. Homology blast analysis of these unigenes revealed that 76% (n=1,092) show significant similarities with annotated gene products and/or known protein domains distributed mainly across amphibian (frogs) species (X. tropicalis=30%, X. laevis=14% and Rana catesbiana=4%).

After adjusting for library size and length bias, the maximum effective number of reads (ER) in Oophaga species ranged from 2,461 (O. anchicayensis single-end) to 5,813 (O. sylvatica – single-end) while in for the cryptic species varied from 82,509 (Rana catesbeiana) to 452,018 (Bufo gargarizans) (higher coverage in cryptic species, Table S2). This indicated that after the optimization of our alignment parameters, a considerable amount of reads from each RNA-seq library actually mapped to the composite de novo transcriptome. After the removal of contigs with low expression profiles (<10 read counts), our RNA-seq dataset was composed of reads mapping to 20,771 genes. Multidimensional scaling plots (stress <0.05) revealed a clear separation by groups (Oophaga vs. cryptic species, p <0.001 Bonferroni corrected) with higher variation among the cryptic species, as expected given the broad taxonomic range of the members of this group (Figure 2).

Figure 2
Multidimensional scaling (MDS) plot of 13 RNA-sequencing samples from skin tissue of Oophaga (n=7, red dots) and cryptic species (n=6, blue dots). RNA-seq data shown here represent the effective read counts (ER) for 31,498 skin-specific contigs. Color ellipses represent a 95% confidence interval.

About 9% of the highly expressed genes (n= 1,931) also showed a phylogenetic signal (i.e. significant fold changes in expression levels between Oophaga and all other anurans). Most of these genes (n=1,656; 85.8%) show higher levels of expression in Oophaga. The clustering analysis (Figure 1) revealed that the toad (B. gargarizans) has the closest transcription profile to Oophaga species. One possible explanation is the presence of alkaloid-based toxins in its skin glands (Qi et al., 2011).

Homology BLAST analyses revealed that 65% (n=1,252) show significant similarities with annotated gene products and/or known protein domains (Table S4; Figure S5A) distributed mainly across amphibian (frogs) species (X. tropicalis=32%, X. laevis=19% and Rana catesbeiana=4%) (Figure S5B). Based on GO mapping level II, these unigenes were assigned to 24 BP, 20 CC and 13 MF (Figure S5C). Multilevel GO term classification assigned the highly expressed unigenes to 16 BP, 5 CC and 9 MF (Figure S5C). Within BP, 13% of unigenes were associated with oxidation-reduction processes, followed by response to stimulus (10%) and regulation of macromolecule metabolic processes (7.3%). Within CC, the integral components of the cell membrane were highly represented (39%), followed by protein complex (27%) and nuclear components (20%). Within MF, most of the differentially expressed unigenes were assigned to metal ion binding activity (21%), followed by protein binding activity and transmembrane transporters (17%), oxidoreductase activity (12%) and structural molecule activity (10%) (Figure S5C). The full functional annotation of these genes allowed us to propose plausible mechanisms that may be involved in the ability to sequester alkaloid-based chemical defenses, warning coloration and the auto-resistance mechanisms to avoid self-intoxication (Table S5).

Discussion

Studies focusing on highly represented and differential gene expression in related lineages have provided valuable insights into the molecular underpinnings of phenotypic variation (Stranger et al., 2007; Gamazon et al., 2015). These analyses rely on the premises that lineages show contrasting phenotypes of interest, and that most differentially expressed genes relate to the trait of interest rather than the overall of genomic differentiation between lineages. Here, we examined shared patterns of gene expression in a small clade of closely related aposematic frogs (Oophaga lineages). Because all members of this monophyletic group show contrasting hues and accumulate alkaloid-derived toxins in their dermal glands (Silverstone, 1973, 1975), we assumed that the genes related with these traits are likely to 1) be highly expressed in the skin of Oophaga poison frogs and 2) showed contrasting expression patterns between aposematic and cryptic and less (or not) poisonous species. Using expression levels in a multivariate framework, we first identified ortholog genes that are highly represented in Oophaga skin tissue and then, we compared with those differentially expressed in the skin of contrasting groups and that show similar patterns of expression among Oophaga species. We finally used a functional annotation analyses to highlight a series of genes that are possibly associated with a series of aposematic traits.

The characterization of our skin transcriptomic profiles allowed us to identify genes that are potentially associated with alkaloid transportation and resistance to autotoxicity in poison frogs (Dendrobatidae). Because the alkaloids secreted by these frogs are structurally very similar to those present in plants (Daly et al., 1994, 2002, 2005, 2009), it is reasonable to assume that the transportation and accumulation of alkaloids in these frogs may be carried out by similar systems to those described in plants in which alkaloids are transferred from source to sink organs (Hashimoto and Yamada, 1994; Hartmann, 1999; Yazaki, 2005). If so, at least, two alternative membrane mechanisms may help to explain the transport and storage of unmodified diet-derived alkaloid by the specialized cells of skin secretory glands. First, alkaloids in a lipophilic state may freely pass through the cell membranes by simple diffusion and accumulate in acidic secretory lysosomes if they become protonated to form hydrophilic cations. This ion-trap mechanism is not energy-dependent and does not necessarily require the expression of any transporters. Alternatively, the transportation of alkaloids may be managed by proton-antiport carrier systems in an energy-requiring manner (Otani, et al., 2005; Shitan and Yazaki, 2007; Carqueijeiro et al., 2013). Under this alternative model, diet derived alkaloids may be taken up by ABC transporters (Sakai et al., 2002; Shitan et al., 2003) and may accumulate in the secretory lysosomes of the skin gland cells by a cation exchanger antiporter (CAX) system, dependent on the pH gradient generated by “vacuolar” type ATPases (V-H⁺-ATPases) and/or pyrophosphatases (V-H⁺-PPases). To find support for these two potential models that might explain the accumulation and posterior secretion of dietary alkaloids, we investigated the common transcriptional profile across Oophaga species. For each of these species, our differential expression and transcript abundance analysis revealed that at least 15 homologs of different type II CAX are highly expressed in skin tissue (Table S5). A result suggesting the existence of an active membrane transportation-accumulation mechanism of alkaloids in harlequin poison frogs. Consistent with this hypothesis, we also found putative ABC transporters and at least one V-H⁺-ATPase on the skin transcriptome. All together, these results suggest that active lysosomal exocytosis may play a key role in the secretion of alkaloids in these frogs.

The accumulation of toxic compounds implies that organisms must avoid self-intoxication (auto-resistance). While the membrane transportation mechanisms described above reduce auto-toxicity by compartmentalizing the sequestered alkaloids, other non-alternative excluding mechanisms are likely to contribute to auto-toxicity resistance. In harlequin poison frogs, the major toxic components present in skin tissue are histrionicoxins (HTX) (Daly et al., 2005). These alkaloids are known to cause temporary paralysis or even death by inactivating or blocking voltage-gate ion channels (Daly et al., 1985). Resistance to this type of cytotoxic compounds usually arises through either, an increased expression level of P450 enzymes (CYPs) that metabolize the toxin, or through target insensitivity via mutations that reduce the toxin's ability to bind to the ion channel itself (Chaudhary and Willett, 2006). Our analyses revealed that five CYP homologs are amongst the most highly expressed genes in the skin of harlequin poison frogs (Table S5). This result suggests that the oxidative biotransformation of lipophilic alkaloids to hydrophilic compounds (Sigel et al., 2007) is an important auto resistance mechanism in these frogs. Although the overproduction of the CYP-enzymes has been interpreted as evidence for metabolic detoxification (Santos et al., 2016), there is an alternative explanation for the high expression levels of CYP proteins in skin of harlequin poison frogs. While many of the alkaloids sequestered in the wild are accumulated unchanged in the dermal glands, at least two species of dendrobatids, stereo-selectively hydroxylate the toxic pumiliotoxin PTX-(+)-251D to convert it into a much more potent toxin, the allopumiliotoxin aPTX-(+)-267A (Daly et al., 2003). Both PTX-(+)-251D and aPTX-(+)-267A have been isolated from the skin of O. histrionica, and CYPs are known to be involved in the stereo-selective hydroxylation of alkaloids (Giddings et al., 2011). Thus, it is also possible that the CYPs expressed in the skin of harlequin poison frogs may enhance the antipredator potency of ingested PTXs.

Despite the major role that detoxification enzymes may have in the resistance to diet acquired toxins, mutations conferring constitutive resistance to alkaloids are also likely to be involved in auto-resistance in chemically defended dendrobatids. In insects, single amino-acid substitution in the voltage-gated sodium channels (Na⁺K⁺-ATPases) are responsible for resistance to host-plant phytochemicals (Labeyrie and Dobler, 2004; Aardema et al., 2012; Dobler et al., 2012). In amphibians, the only documented example of this type of target insensitivity is from distantly-related lineages of alkaloid-defended frogs, 24 species of Neotropical dendrobatids and the aposematic Madagascar frog Mantella aurantiaca (Clark et al., 2005). In all of these cases, they carried different amino acid replacements in the inner pore of the voltage-gated sodium channel Nav1.4 that are not found in other frog (Tarvin et al., 2016). Our analyses revealed eight genes encoding voltage-gated ion channel proteins (VGIC), one of which (8966), encodes the gamma-1 subunit of a sodium channel.

Aposematism in harlequin frogs is a complex phenotype that results from the integration of different elements including unpalatability and conspicuousness. Thus, another of our goals was to identify genes potentially related with warning coloration. Many aposematic organisms, including the Oophaga poison frogs studied here exhibit color patterns that show strong color and/or luminance contrast such between black and red, orange or yellow (Posso-Terranova and Andrés, 2017). It has been proposed that such patterning increases aversion learning as well as conspicuousness and distinctiveness from palatable prey (Bowdish and Bultman, 1993; Gamberale-Stille and Guilford, 2003; Sherratt and Beatty, 2003; Prudic et al., 2007; Aronsson and Gamberale-Stille, 2009, 2012). In anurans, coloration relates to both the structure and pigment composition of the dermal chromatophore units in which three types of pigment cells (xanthophores, iridophores and melanophores) are laid one on another (Bagnara et al., 1968). In dark/black skin areas of harlequin poison frogs, no xanthophores or iridophores are present and the distal fingers of the melanophores are filled with melanin granules (melanosomes) obscuring the dermis (Posso-Terranova and Andrés, 2017). Thus, genes involved in the amount, size and distribution of the melanosomes are likely to play a significant role in the coloration pattern of these frogs.

In one of the few detailed studies of coloration in anurans (Bagnara and Taylor, 1970), the authors suggest a common origin for of all the pigment granules found in the cells of the chromatophore: a primordial organelle (vesicle) derived from the rough endoplasmic reticulum (RER). According to this model, in the formation of melanosomes, the pre-melanosomes are derived from cisternae of the RER, which then fuse with vesicles derived from the Golgi complex containing tyrosinase enzymes (Bagnara and Taylor, 1970; Bagnara et al., 1973, 1978; Bagnara, 1982; Palumbo et al., 1997). A highly expressed gene identified here (tyrosinase regulator, 25062) may contribute to regulate this mechanism, which in turn might translate into hue differences. Interestingly, among Oophaga species, we have found populations that are characterized by light-brown background coloration as opposed to black (Posso-Terranova and Andrés 2017, 2018). One tantalizing possibility is that this phenotypic difference is indeed associated with the expression differences of these tyrosinase enzymes.

Obviously, there are other molecular and cellular mechanisms that might be associated to the difference between light and dark background coloration. In fishes and amphibians, dark hues are known to be produced by the interaction between high levels of melanocyte-stimulating hormone (α-MSH) and several variants of its transmembrane receptor (MC1R) through the dispersion of melanosomes within the melanophore (by increasing cAMP intracellular levels) (Sugimoto, 2002; Logan et al., 2006). Thus, it is possible that structural or expression differences in the MC1R might contribute to dark phenotypes in Oophaga frogs. While there is a strong evidence that different mutations at MC1R cause either light or dark phenotypes in many mammals, birds and reptiles (Everts et al., 2000; Gross et al., 2009; Gangoso et al., 2011; Baião and Parker, 2012; Corso et al., 2012), the only two studies conducted in frogs are inconclusive (Herczeg et al., 2010; Chikako, 2012). A detailed inspection of the coding sequences recovered for this gene revealed the presence of different length isoforms, making MC1R a promising candidate gene candidate to explain the differences in background coloration in poison frogs (Posso-Terranova and Andrés, 2017). Finally, our study revealed another two highly and differentially expressed genes (14447-Melanoregulin-; 2038-Melanocortin phosphoprotein), which may also contribute to dark hues. In this case, the predicted products of these genes are key proteins that mediate the melanosome transport and distribution in epidermal cells through the formation of a tripartite protein complex (Ohbayashi et al., 2012). The disruption of the transport protein complex results in pigmentary dilution and lighter phenotypes by the clustering of melanosomes around the nucleus (Van Gele et al., 2009).

Overall, this study demonstrates the utility of using RNA-sequencing with non-model organisms to identify loci that may be of adaptive importance. Altogether, these data enabled us to provide a first global study of Oophaga poison frogs transcriptomes and propose potential mechanisms for alkaloid sequestration, auto-resistance to toxic compounds, and variation in coloration and patterns. It is important to recognized that due to sampling restrictions (see Methods), we are most likely missing some of the genes responsible for these traits. It is also probable that the genes associated with coloration, alkaloid metabolism, transport and storage are differentially expressed not only in skin, but also in other tissues in the organism. Despite the potential caveats, two recent molecular studies have shown that some of the genes identified in this study are indeed potentially related with aposematic traits in dendrobatids (Tarvin et al., 2016; Posso-Terranova and Andrés, 2017).

Despite the critical impact of the genetic basis of coloration in the evolution and diversification of aposematic phenotypes, the genetic architecture of coloration in poison frogs (Dendrobatidae) remains virtually unexplored. Here, we report probably the first RNA-seq study in Oophaga poison frogs, a model system for understanding the relationship between toxicity, diet, and coloration. The skin-expressed genes that we have identified here provide an initial working hypothesis to further unravel the molecular genetics mechanisms to sequester alkaloid-based chemical defenses, warning coloration, and the auto-resistance mechanisms to avoid self-intoxication. Hence, further analysis aiming to compare the amino acid substitutions in the expressed cDNAs will provide insights about the structure and function of each of these genes. Finally, our comparative transcriptome data provide an important new resource to better understand the evolution of warning signals in nature.

Acknowledgments

We would like to thank the former members of Rick Harrison's lab (requiescat in pace) at Cornell University for their valuable discussions and comments on early versions of the manuscript, which was greatly improved by the thoughtful comments of Dr. Andrew Crawford at Universidad de los Andes (Colombia). This study was funded by COLCIENCIAS (Departamento Administrativo de Ciencia, Tecnología e Innovación, Colombia, 529-2011) grant to A.P.-T. and an NSERC-Discovery grant to J.A.A.

References

• Aardema ML, Zhen Y and Andolfatto P (2012) The evolution of cardenolide-resistant forms of Na+,K+-ATPase in Danainae butterflies. Mol Ecol 21:340-349.
• Anders S and Huber W (2010) Differential expression analysis for sequence count data. Genome Biol 11:R106.
• Aronsson M and Gamberale-Stille G (2009) Importance of internal pattern contrast and contrast against the background in aposematic signals. Behav Ecol 20:1356-1362.
• Aronsson M and Gamberale-Stille G (2012) Evidence of signaling benefits to contrasting internal color boundaries in warning coloration. Behav Ecol 24:349-354.
• Bagnara JT (1982) Development of the spot pattern in the leopard frog. J Exp Zool 224:283-287.
• Bagnara JT and Taylor JD (1970) Differences in pigment-containing organelles between color forms of the red-backed salamander, Plethodon cinereus. Z Zellforsch Mikrosk Anat 106:412-417.
• Bagnara JT, Taylor JD and Hadley ME (1968) The dermal chromatophore unit. J Cell Biol 38:67-79.
• Bagnara JT, Taylor JD and Prota G (1973) Color changes, unusual melanosomes, and a new pigment from leaf frogs. Science 182:1034-1035.
• Bagnara JT, Frost SK and Matsumoto J (1978) On the development of pigment patterns in mmphibians. Am Zool 18:301-312.
• Baião PC and Parker PG (2012) Evolution of the melanocortin-1 receptor (MC1R) in boobies and gannets (Aves, Suliformes). J Hered 103:322-329.
• Berthold AA (1846) Üeber verschiedene neue oder seltene Reptilien aus Neu-Granada und Crustaceen aus China. Nachr Ges Wiss Göttingen 3:43.
• Bowdish TI and Bultman TL (1993) Visual cues used by mantids in learning aversion to aposematically colored prey. Am Midl Nat 129:215-222.
• Brown J (2013) The evolution of parental care, aposematism and color diversity in Neotropical poison frogs. Evol Ecol 27:825-829.
• Carqueijeiro I, Noronha H, Duarte P, Gerós H and Sottomayor M (2013) Vacuolar transport of the medicinal alkaloids from Catharanthus roseus is mediated by a proton-driven antiport. Plant Physiol 162:1486-1496.
• Chaudhary A and Willett KL (2006) Inhibition of human cytochrome CYP1 enzymes by flavonoids of St. John's wort. Toxicology 217:194-205.
• Chikako M (2012) Geographic variations of melanocortine 1 receptor gene (MC1R) in the common frog (Rana temporaria) in Northern Europe. Amphibia-Reptilia 33:105-111.
• Clark VC, Raxworthy CJ, Rakotomalala V, Sierwald P and Fisher BL (2005) Convergent evolution of chemical defense in poison frogs and arthropod prey between Madagascar and the Neotropics. Proc Natl Acad Sci U S A 102:11617-11622.
• Conesa A, Götz S, García-Gómez JM, Terol J, Talón M and Robles M (2005) Blast2GO: A universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 21:3674-3676.
• Corso J, Gonçalves G and Thales dF (2012) Sequence variation in the melanocortin-1 receptor (MC1R) pigmentation gene and its role in the cryptic coloration of two South American sand lizards. Genet Mol Biol 35:81-87.
• Cummings M and Crothers L (2013) Interacting selection diversifies warning signals in a polytypic frog: an examination with the strawberry poison frog. Evol Ecol 27:693-710.
• Daly JW (1995) The chemistry of poisons in amphibian skin. Proc Natl Acad Sci U S A 92:9-13.
• Daly JW, Brown GB, Mensah-Dwumah M and Myers CW (1978) Classification of skin alkaloids from neotropical poison-dart frogs (Dendrobatidae). Toxicon 16:163-188.
• Daly JW, McNeal ET, Overman LE and Ellison DH (1985) A new class of cardiotonic agents: structure-activity correlations for natural and synthetic analogs of the alkaloid pumiliotoxin B (8-hydroxy-8-methyl-6-alkylidene-1-azabicyclo[4.3.0]nonanes). J Med Chem 28:482-486.
• Daly JW, Martin Garraffo H, Spande TF, Jaramillo C and Stanley Rand A (1994) Dietary source for skin alkaloids of poison frogs (Dendrobatidae). J Chem Ecol 20:943-955.
• Daly JW, Kaneko T, Wilham J, Garraffo HM, Spande TF, Espinosa A and Donnelly MA (2002) Bioactive alkaloids of frog skin: Combinatorial bioprospecting reveals that pumiliotoxins have an arthropod source. Proc Natl Acad Sci U S A 99:13996-14001.
• Daly JW, Garraffo HM, Spande TF, Clark VC, Ma J, Ziffer H and Cover JF (2003) Evidence for an enantioselective pumiliotoxin 7-hydroxylase in dendrobatid poison frogs of the genus Dendrobates Proc Natl Acad Sci U S A 100:11092-11097.
• Daly JW, Spande TF and Garraffo HM (2005) Alkaloids from amphibian skin:? A tabulation of over eight-hundred compounds. J Nat Prod 68:1556-1575.
• Daly JW, Ware N, Saporito RA, Spande TF and Garraffo HM (2009) N-Methyldecahydroquinolines: An unexpected class of alkaloids from Amazonian poison frogs (Dendrobatidae). J Nat Prod 72:1110-1114.
• Dobler S, Dalla S, Wagschal V and Agrawal AA (2012) Community-wide convergent evolution in insect adaptation to toxic cardenolides by substitutions in the Na,K-ATPase. Proc Natl Acad Sci U S A 109:13040-13045.
• Dodt M, Roehr J, Ahmed R and Dieterich C (2012) FLEXBAR—Flexible barcode and adapter processing for next-generation sequencing platforms. Biology 1:895-905.
• Eskew EA, Shock BC, LaDouceur EEB, Keel K, Miller MR, Foley JE and Todd BD (2018) Gene expression differs in susceptible and resistant amphibians exposed to Batrachochytrium dendrobatidis. R Soc Open Sci 5:170910.
• Everts RE, Rothuizen J and van Oost BA (2000) Identification of a premature stop codon in the melanocyte-stimulating hormone receptor gene (MC1R) in Labrador and Golden retrievers with yellow coat colour. Anim Genet 31:194-199.
• Funkhouser J (1956) New frogs from Ecuador and southwestern Colombia. Zoologica 41:73-80.
• Gamazon ER, Wheeler HE, Shah KP, Mozaffari SV, Aquino-Michaels K, Carroll RJ, Eyler AE, Denny JC, Consortium GT, Nicolae DL et al. (2015) A gene-based association method for mapping traits using reference transcriptome data. Nat Genet 47:1091-1098.
• Gamberale-Stille G and Guilford T (2003) Contrast versus colour in aposematic signals. Anim Behav 65:1021-1026.
• Gangoso L, Grande JM, Ducrest AL, Figuerola J, Bortolotti GR, Andres JA and Roulin A (2011) MC1R-dependent, melanin-based colour polymorphism is associated with cell-mediated response in the Eleonora's falcon. J Evol Biol 24:2055-2063.
• Giddings LA, Liscombe DK, Hamilton JP, Childs KL, DellaPenna D, Buell CR and O’Connor SE (2011) A stereoselective hydroxylation step of alkaloid biosynthesis by a unique cytochrome P450 in Catharanthus roseus J Biol Chem 286:16751-16757.
• Gould B, McCouch S and Geber M (2015) De novo transcriptome assembly and identification of gene candidates for rapid evolution of soil Al tolerance in Anthoxanthum odoratum at the long-term park grass experiment. PLoS One 10:e0124424.
• Grant T, Frost DR, Caldwell JP, Gagliardo R, Haddad CFB, Kok PJR, Means DB, Noonan BP, Schargel WE and Wheeler WC (2006) Phylogenetic systematics of dart-poison frogs and their relatives (Amphibia: Athesphatanura: Dendrobatidae). Bull Amer Mus Nat Hist 299:1-262.
• Gross JB, Borowsky R and Tabin CJ (2009) A novel role for MC1R in the parallel evolution of depigmentation in independent populations of the Cavefish Astyanax mexicanus PLoS Genet 5:e1000326.
• Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, Couger MB, Eccles D, Li B, Lieber M et al. (2013) De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat Protocols 8:1494-1512.
• Hammer Ø, Harper D and Ryan P (2001) PAST: Paleontological statistics software package for education and data analysis, https://palaeo-electronica.org/2001_1/past/issue1_01.htm
  » https://palaeo-electronica.org/2001_1/past/issue1_01.htm
• Hartmann T (1999) Chemical ecology of pyrrolizidine alkaloids. Planta 207:483-495.
• Harvey PH, Bull JJ, Pemberton M and Paxton RJ (1982) The evolution of aposematic coloration in distasteful prey: A family model. Am Nat 119:710-719.
• Hashimoto T and Yamada Y (1994) Alkaloid biogenesis: Molecular aspects. Annu Rev Plant Biol 45:257-285.
• Herczeg G, Matsuba C and Merilä J (2010) Sequence variation in the melanocortin-1 receptor gene (Mc1r) does not explain variation in the degree of melanism in a widespread amphibian. Ann Zool Fennici 47:37-45.
• Holding ML, Margres MJ, Mason AJ, Parkinson CL and Rokyta DR (2018) Evaluating the performance of De novo assembly methods for venom-gland transcriptomics. Toxins 10:249.
• Huang L, Li J, Anboukaria H, Luo Z, Zhao M and Wu H (2016) Comparative transcriptome analyses of seven anurans reveal functions and adaptations of amphibian skin. Sci Rep 6:24069.
• Labeyrie E and Dobler S (2004) Molecular adaptation of Chrysochus leaf beetles to toxic compounds in their food plants. Mol Biol Evol 21:218-221.
• Langham GM and Benkman C (2004) Specialized avian predators repeatedly attack novel color morphs of Heliconius butterflies. Evolution 58:2783-2787.
• Langmead B and Salzberg SL (2012) Fast gapped-read alignment with Bowtie 2. Nat Methods 9:357.
• Logan DW, Burn SF and Jackson IJ (2006) Regulation of pigmentation in zebrafish melanophores. Pigment Cell Res 19:206-213.
• Maan ME and Cummings ME (2012) Poison frog colors are honest signals of toxicity, particularly for bird predators. Am Nat 179:E1-E14.
• Medina I, Wang IJ, Salazar C and Amézquita A (2013) Hybridization promotes color polymorphism in the aposematic harlequin poison frog, Oophaga histrionica Ecol Evol 3:4388-4400.
• Myers CW and Daly JW (1976) Preliminary evaluation of skin toxins and vocalizations in taxonomic and evolutionary studies of poison-dart frogs (Dendrobatidae). Bull Amer Mus Nat Hist 157:175-262.
• Nanoth Vellichirammal N, Zera AJ, Schilder RJ, Wehrkamp C, Riethoven JJM and Brisson JA (2014) De novo transcriptome assembly from fat body and flight muscles transcripts to identify morph-specific gene expression profiles in Gryllus firmus PLoS One 9:e82129.
• Nicholls SM, Aubrey W, de Grave K, Schietgat L, Creevey CJ and Clare A (2017) Probabilistic recovery of cryptic haplotypes from metagenomic data. bioRxiv:117838
• Nishimura O, Hara Y and Kuraku S (2017) gVolante for standardizing completeness assessment of genome and transcriptome assemblies. Bioinformatics 33:3635-3637.
• Ohbayashi N, Maruta Y, Ishida M and Fukuda M (2012) Melanoregulin regulates retrograde melanosome transport through interaction with the RILP–p150 glued complex in melanocytes. J Cell Sci 125:1508-1518.
• Otani M, Shitan N, Sakai K, Martinoia E, Sato F and Yazaki K (2005) Characterization of vacuolar transport of the endogenous alkaloid berberine in Coptis japonica Plant Physiol 138:1939-1946.
• Palumbo A, di Cosmo A, Gesualdo I and Vincent J (1997) Subcellular localization and function of melanogenic enzymes in the ink gland of Sepia officinalis Biochem J 323:749-756.
• Posso-Terranova A and Andrés JA (2016) Complex niche divergence underlies lineage diversification in Oophaga poison frogs. J Biogeogr 43:2002-2015.
• Posso-Terranova A and Andrés JA (2017) Diversification and convergence of aposematic phenotypes: Truncated receptors and cellular arrangements mediate rapid evolution of coloration in harlequin poison frogs. Evolution 71:2677-2692.
• Posso-Terranova A and Andrés J (2018) Multivariate species boundaries and conservation of harlequin poison frogs. Mol Ecol 27:3432-3451.
• Prudic KL, Skemp AK and Papaj DR (2007) Aposematic coloration, luminance contrast, and the benefits of conspicuousness. Behav Ecol 18:41-46.
• Qi F, Li A, Inagaki Y, Kokudo N, Tamura S, Nakata M and Tang W (2011) Antitumor activity of extracts and compounds from the skin of the toad bufo Bufo gargarizans Cantor. Int Immunopharmacol 11:342-349.
• Ramos RR and Freitas A (1999) Population biology and wing color variation in Heliconius erato phyllis (Nymphalidae). J Lepid Soc 53:11-21.
• Roberts A and Pachter L (2013) Streaming fragment assignment for real-time analysis of sequencing experiments. Nat Meth 10:71-73.
• Rogers RL, Summers K, Wu Y, Guo C, Zheng J, Xun X, Xiong Z, Yang H, Zhou L, Zhang G et al. (2018) Genomic takeover by transposable elements in the strawberry poison frog. Mol Biol Evol 35:2913-2927.
• Sakai K, Shitan N, Sato F, Ueda K and Yazaki K (2002) Characterization of berberine transport into Coptis japonica cells and the involvement of ABC protein. J Exp Bot 53:1879-1886.
• Santos JC, Tarvin RD and O’Connell LA (2016). A review of chemical defense in poison frogs (Dendrobatidae): Ecology, pharmacokinetics, and autoresistance. In: Schulte AB, Goodwin ET and Ferkin HM (eds) Chemical Signals in Vertebrates 13. Springer International Publishing, Cham, pp 305-337.
• Santos JC, Tarvin RD, O’Connell LA, Blackburn DC and Coloma LA (2018) Diversity within diversity: Parasite species richness in poison frogs assessed by transcriptomics. Mol Phylogenet Evol 125:40-50.
• Saporito R, Donnelly M, Garraffo H, Spande T and Daly J (2006) Geographic and seasonal variation in alkaloid-based chemical defenses of Dendrobates pumilio from Bocas del Toro, Panama. J Chem Ecol 32:795-814.
• Saporito R, Donnelly M, Spande T and Garraffo HM (2012) A review of chemical ecology in poison frogs. Chemoecology 22:159-168.
• Schurch NJ, Schofield P, Gierlinski M, Cole C, Sherstnev A, Singh V, Wrobel N, Gharbi K, Simpson GG, Owen-Hughes T et al. (2016) How many biological replicates are needed in an RNA-seq experiment and which differential expression tool should you use? RNA 22:839-851.
• Sherratt T and Beatty C (2003) The evolution of warning signals as reliable indicators of prey defense. Am Nat 162:377-389.
• Shitan N and Yazaki K (2007) Accumulation and membrane transport of plant alkaloids. Curr Pharm Biotechnol 8:244-252.
• Shitan N, Bazin I, Dan K, Obata K, Kigawa K, Ueda K, Sato F, Forestier C and Yazaki K (2003) Involvement of CjMDR1, a plant multidrug-resistance-type ATP-binding cassette protein, in alkaloid transport in Coptis japonica Proc Natl Acad Sci U S A 100:751-756.
• Sigel A, Sigel H and Sigel RK (2007). The ubiquitous roles of cytochrome P450 proteins: metal ions in life sciences. John Wiley & Sons, New York.
• Silverstone PA (1973) Observations on the behavior and ecology of a Colombian poison-arrow frog, the Kõkoé-Pá (Dendrobates histrionicus Berthold). Herpetologica 29:295-301.
• Silverstone PA (1975) A revision of the poison-arrow frogs of the genus Dendrobates Wagler. Am Nat 21:1-55.
• Stranger BE, Forrest MS, Dunning M, Ingle CE, Beazley C, Thorne N, Redon R, Bird CP, de Grassi A, Lee C et al. (2007) Relative impact of nucleotide and copy number variation on gene expression phenotypes. Science 315:848-853.
• Sugimoto M (2002) Morphological color changes in fish: regulation of pigment cell density and morphology. Microsc Res Tech 58:496-503.
• Tarvin RD, Santos JC, O’Connell LA, Zakon HH and Cannatella DC (2016) Convergent substitutions in a sodium channel suggest multiple origins of toxin resistance in poison frogs. Mol Biol Evol 33:1068-1081.
• Van Gele M, Dynoodt P and Lambert J (2009) Griscelli syndrome: A model system to study vesicular trafficking. Pigment Cell Melanoma Res 22:268-282.
• Wagner GP, Kin K and Lynch VJ (2012) Measurement of mRNA abundance using RNA-seq data: RPKM measure is inconsistent among samples. Theory Biosci 131:281-285.
• Warnes GR, Bolker B, Bonebakker L, Gentleman R, Liaw WHA, Lumley T, Maechler M, Magnusson A, Moeller S and Venables MS (2016) gplots: Various R programming tools for plotting data. R Package Version 3.0.1.
• Waterhouse RM, Seppey M, Simão FA, Manni M, Ioannidis P, Klioutchnikov G, Kriventseva EV and Zdobnov EM (2018) BUSCO Applications from quality assessments to gene prediction and phylogenomics. Mol Biol Evol 35:543-548.
• Wolnicka-Glubisz A, Pecio A, Podkowa D, Kolodziejczyk LM and Plonka PM (2012) Pheomelanin in the skin of Hymenochirus boettgeri (Amphibia: Anura: Pipidae). Exp Dermatol 21:537-540.
• Xie H, Wasserman A, Levine Z, Novik A, Grebinskiy V, Shoshan A and Mintz L (2002) Large-scale protein annotation through Gene Ontology. Genome Res 12:785-794.
• Yazaki K (2005) Transporters of secondary metabolites. Curr Opin Plant Biol 8:301-307.
• Zheng Y, Zhao L, Gao J and Fei Z (2011) iAssembler: A package for de novo assembly of Roche-454/Sanger transcriptome sequences. BMC Bioinformatics 12:1-8.

Supplementary material

The following online material is available for this article:

Table S1 - Number of reads and library type for each individual RNA-seq experiment.

Table S2 - RNA-sequence datasets from cryptic species used for the differential expression analysis.

Table S3 - Quality assessment report for each individual and the composite reference transcriptome.

Table S4 - Summary statistics of the differentially expressed unigenes.

Table S5 - Transcripts with potential important functions in the alkaloid sequestration system, auto-resistance to toxic compounds, and variation in coloration in harlequin poison frogs.

Figure S1 - Pie charts representing the number and proportion of contigs with significant BLAST hits (E<1.0E-5) for each individual transcriptome from Oophaga species.

Figure S2 - Contig length distribution of the de-novo composite reference transcriptome of Oophaga species, pie chart representing the number and proportion of contigs with significant BLAST hits, and Gene ontology categories.

Figure S3 - Distribution plots of the effective number of reads from Oophaga RNA-seq experiments.

Figure S4 - Percentile plot of the estimated transcripts per million in the composite reference transcriptome and principal component analysis plot of the reference contigs dataset.

Figure S5 - Pie chart representing the total number of differentially expressed unigenes, BLAST hit species distribution of highly expressed unigenes, and GO categories distribution.

Publication Dates

• Publication in this collection
  16 Nov 2020
• Date of issue
  2020

History

• Received
  17 Dec 2019
• Accepted
  06 June 2020