The occasional reversal of sex-typical behavior suggests that many of the neural circuits underlying behavior are conserved between males and females and can be activated in response to the appropriate social condition or stimulus. Most poison frog species (Family Dendrobatidae) exhibit male uniparental care, but flexible compensation has been observed in some species, where females will take over parental care duties when males disappear. We investigated hormonal and neural correlates of sex-typical and sex-reversed parental care in a typically male uniparental species, the Dyeing Poison Frog (Dendrobates tinctorius). We first characterized hormone levels and whole brain gene expression across parental care stages during sex-typical care. Surprisingly, hormonal changes and brain gene expression differences associated with active parental behavior in males were mirrored in their non-caregiving female partners. To further explore the disconnect between neuroendocrine patterns and behavior, we characterized hormone levels and neural activity patterns in females performing sex-reversed parental care. In contrast to hormone and gene expression patterns, we found that patterns of neural activity were linked to the active performance of parental behavior, with sex-reversed tadpole transporting females exhibiting neural activity patterns more similar to those of transporting males than non-caregiving females. We suggest that parallels in hormones and brain gene expression in active and observing parents are related to females’ ability to flexibly take over parental care in the absence of their male partners.The occasional reversal of sex-typical behavior suggests that many of the neural circuits underlying behavior are conserved between males and females and can be activated in response to the appropriate social condition or stimulus. Most poison frog species (Family Dendrobatidae) exhibit male uniparental care, but flexible compensation has been observed in some species, where females will take over parental care duties when males disappear. We investigated hormonal and neural correlates of sex-typical and sex-reversed parental care in a typically male uniparental species, the Dyeing Poison Frog (Dendrobates tinctorius). We first characterized hormone levels and whole brain gene expression across parental care stages during sex-typical care. Surprisingly, hormonal changes and brain gene expression differences associated with active parental behavior in males were mirrored in their non-caregiving female partners. To further explore the disconnect between neuroendocrine patterns and behavior, we characterized hormone levels and neural activity patterns in females performing sex-reversed parental care. In contrast to hormone and gene expression patterns, we found that patterns of neural activity were linked to the active performance of parental behavior, with sex-reversed tadpole transporting females exhibiting neural activity patterns more similar to those of transporting males than non-caregiving females. We suggest that parallels in hormones and brain gene expression in active and observing parents are related to females’ ability to flexibly take over parental care in the absence of their male partners.

Tissue collection

D. tinctorius poison frogs were housed in breeding pairs to ensure parental/partner identity. To control for potential effects of parental experience, we allowed all pairs to successfully rear at least one clutch from egg-laying through tadpole transport prior to the study. For the ‘no care’ group, we collected frog pairs between parental bouts when they were not caring for eggs or tadpoles (N=8 pairs for hormones, N=5 pairs for RNAseq). For the ‘egg care’ group, we collected frogs one week after they laid a clutch of eggs (N=7-10 pairs for hormones, N=5 pairs for RNAseq). We collected only frogs that were caring for eggs, as evidenced by the healthy development of embryos. For the ‘tadpole transport’ group, we checked daily for transport behavior and collected tissue when we found frogs actively transporting (N=7-8 pairs for hormones, N=5 pairs for RNAseq). For all behavioral groups we simultaneously collected tissue from parental males and their non-caregiving female partners.

Library construction and sequencing

Flash frozen brains were placed in Trizol (Life Technologies, Grand Island, NY) and RNA was extracted according to manufacturer instructions. RNA quality was confirmed on an Agilent 2100 Bioanalyzer (Santa Clara, CA, USA) where all samples had a RIN (RNA Integrity Number) score >7; one sample (a male transporter) had a RIN score below 5 and was excluded from the further tissue processing and gene expression analysis. Poly-adenylated RNA was isolated from each sample using the NEXTflex PolyA Bead kit (Bioo Scientific, Austin, TX, USA) following manufacturer guidelines. Lack of contaminating ribosomal RNA was confirmed using an Agilent 2100 Bioanalyzer. Strand specific libraries for each sample were prepared using the dUTP NEXTflex RNAseq kit (Bioo Scientific), which includes a magnetic bead-based size selection. Libraries were pooled in equimolar amounts after library quantification using both quantitative PCR with the KAPA Library Quantification Kit (KAPA Biosystems, Wilmington, MA, USA) and the fluorometric Qubit dsDNA high sensitivity assay kit (Life Technologies), both according to manufacturer instructions. Libraries were sequenced on an Illumina HiSeq 2500 at the Bauer Sequencing Core at Harvard University.

Transcriptome construction and annotation

We constructed de novo transcriptomes using a pipeline established in our laboratory (Caty et al., 2019). We amended sequencing errors using Rcorrector (Song and Florea, 2015) and trimmed raw reads using Trim Galore! (Babraham Bioinformatics, Babraham Institute) to remove Illumina adapters and restrict all reads to high-quality sequence. Following developer recommendations, we used a quality score of 33, a stringency of 5, and a minimum read length of 36 bp. Corrected, trimmed reads from all individuals were pooled prior to de novo transcriptome construction with Trinity (Grabherr et al, 2011;  Haas et al, 2013). Our initial assembly contained 2,655,485 isoforms representing 2,272,431 presumptive genes.

We filtered the raw assembly using several approaches. We clustered overlapping contigs using CD-HIT-EST (http://weizhongli-lab.org/cd-hit/) and removed contigs that were smaller than 250 bp following clustering, leaving 1,479,791 total contigs. Preliminary annotations of this and other assemblies from closely related poison frog species revealed a high percentage of contaminant (i.e. non-vertebrate) sequences. To remove non-vertebrate (i.e. bacterial, fungal and other pathogen/parasite) contaminants, we annotated sequences using blastx queries with default parameters and an e-value cutoff of 10^-10 against the SwissProt database (www.uniprot.org) and retained only those contigs with annotations to known vertebrate genes. Our final assembly contained 68,399 contigs representing 35,058 genes. Additional assembly annotation was done using Trinotate (Bryant et al, 2017) and assembly completeness was assessed using BUSCO (Simão et al, 2015). Final assembly completeness was estimated at 87% (37% duplicated; 3.5% fragmented, 9.0% missing). All high-powered computing for transcriptome assembly and filtering was performed on the Odyssey computing cluster supported by the FAS Science Division Research Computing Group at Harvard University.