A common goal in evolutionary biology is to discern the mechanisms that produce the astounding diversity of morphologies seen across the tree of life. Aposematic species, those with a conspicuous phenotype coupled with some form of defense, are excellent models to understand the link between vivid color pattern variations, the natural selection shaping it, and the underlying genetic mechanisms underpinning this variation. Mimicry systems in which species share a conspicuous phenotype can provide an even better model for understanding the mechanisms of color production in aposematic species, especially if comimics have divergent evolutionary histories. Here we investigate the genetic mechanisms by which mimicry is produced in poison frogs. We assembled a 6.02 Gbp genome with a contig N50 of 310 Kbp, a scaffold N50 of 390 Kbp, and 85% of expected tetrapod genes. We leveraged this genome to conduct gene expression analyses throughout development of four color morphs of R. imitator and two color morphs from both R. fantastica and R. variabilis which R. imitator mimics. We identified a large number of pigmentation and patterning genes differentially expressed throughout development, many of them related to melanophores/melanin, iridophore development, and guanine synthesis. We also  identify the pteridine synthesis pathway (including genes such as qdpr and xdh) as a key driver of the variation in color between morphs of these species, and identify several plausible candidates for coloration in vertebrates (e.g., cd36, ep-cadherin, perlwapin). Finally, we hypothesize that keratin genes (e.g., krt8) are important for producing different structural colors within these frogs.

Gene expression sample preparation:

                Samples were prepared differently for the mimic (R. imitator) and the model species (R. fantastica and R. variabilis). During the course of our work, we discovered that there were multiple groups approaching the same questions using collected samples from different species, but at slightly different timepoints. In light of this, we chose to combine our efforts into a single manuscript in an attempt at making broader inferences. We acknowledge this, and as a result, the data in this manuscript are analyzed in a manner concordant with these differences.

Ranitomeya imitator:

The initial breeding stock of Ranitomeya imitator was purchased from Understory Enterprises, LLC (Chatham, Canada). Frogs used in this project represent captive-bred individuals sourced from the following wild populations going clockwise from top left in Figure 1: Tarapoto (green-spotted), Sauce (orange-banded), Varadero (redheaded), and Baja Huallaga (yellow-striped). Tadpoles were sacrificed for analyses at 2, 4, 7, and 8 weeks of age. We sequenced RNA from a minimum of three individuals at each time point from the Sauce, Tarapoto, and Varadero populations (except for Tarapoto at 8 weeks), and two individuals per time point from the Huallaga population. Individuals within the same time points were sampled from different family groups (Appendix Table 1).

Tadpoles were anesthetized with 20% benzocaine (Orajel), then sacrificed via pithing. Whole skin was removed and stored in RNA later (Ambion) at -20°C until RNA extraction. Whole skin was lysed using a BeadBug (Benchmark Scientific, Sayreville, NJ, USA), and RNA was then extracted using a standardized Trizol protocol. RNA was extracted from the whole skin using a standardized Trizol protocol, cleaned with DNAse and RNAsin, and purified using a Qiagen RNEasy mini kit. RNA Libraries were prepared using standard poly-A tail purification with Illumina primers, and individually barcoded using a New England Biolabs Ultra Directional kit as per the manufacturer’s protocol. Individually barcoded samples were pooled and sequenced using 50 bp paired end reads on three lanes of the Illumina HiSeq 2500 at the New York Genome Center.

Ranitomeya fantastica and Ranitomeya variabilis:

We set up a captive colony in Peru (see Appendix) consisting of between 6 and 10 wild collected individuals per locality. We raised the tadpoles on a diet consisting of a 50/50 mix of powdered spirulina and nettle, which they received 5 times a week. Tadpoles were raised individually in 21oz plastic containers, within outside insectaries covered with 50% shading cloth, and water change was performed with rainwater.  Three tadpoles per stage (1, 2, 5, 7, and 8 weeks after hatching; see Appendix Table 1) were fixed in an RNAlater (Ambion) solution. To do so, tadpoles were first euthanized in a 250 mg/L benzocaine hydrochloride bath, then rinsed with distilled water before the whole tadpole was placed in RNAlater and stored at 4°C for 6h before being frozen at -20°C for long-term storage. Before RNA extraction, tadpoles were removed from RNA later and the skin was dissected off. Whole skin was lysed using a Bead Bug, and RNA was then extracted using a standardized Trizol protocol. RNA libraries were prepared using standard poly-A tail purification, prepared using Illumina primers, and individually dual-barcoded using a New England Biolabs Ultra Directional kit. Individually barcoded samples were pooled and sequenced on four lanes of an Illumina HiSeq X at NovoGene (California, USA). Reads were paired end and 150 base pairs in length.

Differential gene expression:

                We indexed our new Ranitomeya imitator genome using STAR version2.7.10a_alpha_220601 (Dobin et al. 2013). We aligned our reads to our genome using STAR version2.7.10a_alpha_220601 (Dobin et al. 2013), allowing 10 base mismatches (--outFilterMismatchNmax 10), a maximum of 5 multiple alignments per read (--outFilterMultimapNmax 5), and discarding reads that mapped at less than 50% of the read length (--outFilterScoreMinOverLread 0.5). We then counted aligned reads using htseq-count version 2.0.2 (Anders, Pyl, and Huber 2015).

References:

Anders, Simon, Paul Theodor Pyl, and Wolfgang Huber. 2015. “HTSeq—a Python Framework to Work with High-Throughput Sequencing Data.” Bioinformatics  31 (2): 166–69.

Dobin, Alexander, Carrie A. Davis, Felix Schlesinger, Jorg Drenkow, Chris Zaleski, Sonali Jha, Philippe Batut, Mark Chaisson, and Thomas R. Gingeras. 2013. “STAR: Ultrafast Universal RNA-Seq Aligner.” Bioinformatics  29 (1): 15–21.