Dataset DOI: 10.5061/dryad.ttdz08m9p

Description of the data and file structure

Partner cues and individual variation underlie sex-reversed parental care in poison frogs

This dataset contains all raw video/audio scoring tabulations, gene count and annotation data, and relevant trial and sampling information used in the characterization of parental behaviors and associated gene expression profiles of Ranitomeya imitator, contextualized in experimental procedures and results. Nine files (described in detail in Files/Variables) collectively serve as inputs for three scripts, Behavioral_Analysis_Moss_et_al_PROC_B.R, Hormone_Analysis_Moss_et_al_PROC_B.R, and Gene_Expression_Analysis_Moss_et_al_PROC_B.R, which include all data visualizations and statistical analyses presented in the manuscript.

Files and variables

File: Male_Call_Rates.csv

Description: contains all blind scoring data and relevant trial information pertaining to the quantification of male call rates over three time periods for a subset of behavioral trials recorded in the presence (N=10 trials) and absence (N=10 trials) of females.

Variables

• Trial Number: A numeric designation given to the trial for the purposes of blind scoring, with trial types and time periods numbered in a random order.
• Focal Tank: An alphanumeric designation given to the tank in which the trial was conducted.
• Trial Type: One of two possible trial types described in the Methods of the manuscript, wherein males were present in the tank: Control = behavioral baseline, wherein the female partner was present. Female Removal = the female was removed to observe the behavior of single males.
• Stage: One of three possible time periods within a trial, described in the Methods of the manuscript: Pre-Hatching = the last full (12-hr) day before hatch day. Pre-Transport = the day of hatching up to the moment of tadpole pickup as identified from the video analysis. Transport = the time that the tadpole was in transit (i.e., pick-up to drop-off) as identified from the video analysis.
• Call Rate (Calls/Hour): Male call rate during the time period, quantified by dividing the total number of calls tallied by the duration of the trimmed and filtered track (in hours).    

File: Scoring_Info.csv

Description: contains all video scoring data and trial information pertaining to the 71 behavioral trials conducted between October 2021 and June 2023 using established pairs of R. imitator.

Variables

• Male ID: ID of the male subject (“Ri.####”)    
• Female ID: ID of the female subject (“Ri.####”)
• Pair: Compound variable representing the unique combination of Male and Female ID.
• Trial Type: One of six possible trial types described in the Methods of the manuscript: Control = behavioral baseline, wherein male and female behavior was observed in the absence of any manipulation. Male Removal = the male was removed to observe the behavior of single females. Female Removal = the female was removed to observe the behavior of single males. Visual = males were removed, and females were presented with a 3D-printed male dummy painted to match their mate’s coloration and pattern. Playback = males were removed, and females were presented with playback of their mate’s advertisement calls. Playback + Visual = males were removed, and females were presented with both a 3D dummy and playback.
• Oviposition Date: Date on which focal clutch for the trial was discovered (Day 0), recorded during daily checks.
• Egg Position: Compound variable representing the number of eggs in the clutch and its location in the tank (L = left egg canister, R = right egg canister).
• Hatch Date: First date on which an egg(s) in the focal clutch was discovered hatched (Day 8), recorded during daily checks.
• Days Before Hatching: Hatch Date minus Oviposition Date, in days, if both are known. If one or the other is not known, this field is reported as a range.
• Video to Score (last full day before hatch): The date used to score pre-hatching behavior, if video was available. If the day before Hatch Date fell outside the recording window (<day 9), this field is NA.
• Tad Transported: A binary variable reporting the trial outcome: N = the tadpole was not transported within 6 days of hatching; Y = the tadpole was transported to a pool within 6 days of hatching.
• Transport Date: Date on which a tadpole from the focal clutch was first discovered in a pool or video revealed a tadpole deposition event, whichever came first. If the tadpole was not transported (Tad Transported = N), this field is NA.
• Transport Day: Transport Date minus Oviposition Date, in days, if both are known. Otherwise, NA.
• Transport on Video: A binary variable reporting whether the event of tadpole transport was captured on video: N = transport was not captured, possibly due to behaviors occurring outside the light cycle or to other technical problems; Y = transport was captured on video. For trials where the tadpole was not transported (Tad Transported = N), this field is NA.
• Latency to Pickup (Days): Transport Date minus Hatch Date, in days, if both are known and if transport occurred. Otherwise, NA.
• Pickup Time (from first entry of egg canister): The last time the focal frog was seen entering the egg canister before it emerged with a tadpole on its back, scored using video. If the tadpole was not transported (Tad Transported = N), this field is NA.
• Dropoff Time (from entry of deposition pool): The last time the focal frog entered a deposition pool before emerging without a tadpole, scored using video. If the tadpole was not transported (Tad Transported = N), this field is NA.
• Transit Time (from pickup to dropoff): Pickup Time – Dropoff Time (in hh:mm:ss). If the tadpole was not transported (Tad Transported = N), this field is NA.
• No. of failed pool visits before deposition: Number of times a transporting frog was observed visiting a pool without depositing a tadpole, prior to deposition, scored using video. If the tadpole was not transported (Tad Transported = N), this field is NA.

File: Behavior_on_last_full_day.csv

Description: contains all video scoring data and relevant trial information pertaining to the quantification of pre-hatching parental behaviors, recorded on the last full 12-hr light cycle before hatch day.

Variables

• Pair: Compound variable representing the unique combination of Male and Female ID.
• Trial Type: One of six possible trial types described in the Methods of the manuscript: Control = behavioral baseline, wherein male and female behavior was observed in the absence of any manipulation. Male Removal = the male was removed to observe the behavior of single females. Female Removal = the female was removed to observe the behavior of single males. Visual = males were removed, and females were presented with a 3D-printed male dummy painted to match their mate’s coloration and pattern. Playback = males were removed, and females were presented with playback of their mate’s advertisement calls. Playback + Visual = males were removed, and females were presented with both a 3D dummy and playback.
• Date: Date from which video was scored to quantify pre-hatching parental behaviors.
• ID: Identity of the focal individual performing the behavior, either the male or female.
• Behavior: Coded variable corresponding to one of four possible behaviors: Egg care (EC; visiting the focal clutch, defined as full body inside of the canister), and Left Pool/Right Pool (LP/RP); visiting deposition pools, defined as full body inside of the pool).
• Start: Timestamp in video at which a single instance of a behavior was first observed (in hh:mm:ss).
• End: Corresponding timestamp in video at which the present instance of a behavior concluded (in hh:mm:ss).

File: Rimi_plasma_hormones.csv

Description: contains all MyAssay data and relevant sampling information related to assays of cortisol concentrations in 24 plasma samples collected from frogs at sacrifice.

Variables

• sample id: An alphanumeric designation given to the tube containing the plasma sample to be assayed.
• tank id: An alphanumeric designation given to the tank from which the frog was sampled.
• frog id: Identity of the sampled frog.
• sex: Sex of the sampled frog.
• group: One of four sampling groups to which the frog was assigned: Control, or egg care (i.e., males (N=3) and females (N=5) performing egg care only, sampled 9 days post-oviposition), Sex-Typical (i.e., male transporting (N=4) and female observing (N=4), sampled 0–1 days post-hatching), Sex-Reversed-Flexible (i.e., female transporting in absence of male (N=5), sampled 0–1 days post-hatching), and Sex-Reversed-Inflexible (i.e., non-transporting female in absence of male (N=3), sampled 3 days post-hatching).
• approx. plasma vol (uL): The volume, in uL, of plasma recovered from the whole blood sample following refrigerated centrifugation in a heparinized tube.
• plate date: Date on which samples were assayed using competitive ELISA immunoassay kits.
• plasma vol used cortisol: The volume of plasma resuspended in assay buffer for quantification, determined by the starting volume of plasma.
• diss reagent vol: The volume of dissociation reagent combined with the plasma prior to resuspension, which was of equal volume to the plasma volume used.
• assay buffer vol: The volume of assay buffer added to the plasma:dissociation reagent mixture to bring all samples to a near-equal concentration.
• total vol: The final volume of the plasma:dissociation reagent:assay buffer mixture, either 200 or 400 μL, determined by the starting volume of plasma.
• dilution factor: The dilution factor of the plasma:dissociation reagent:assay buffer mixture, either 80 or 100, determined by the starting volume of plasma.
• EIA notes: miscellaneous notes recorded for samples while performing the assay.
• EIA cv: Coefficient of variability calculated among sample replicates and demonstrated as variation among Optical Density (OD) readouts post-assay.
• EIA avg conc (pgmL): Average concentration of cortisol in pg/mL, calculating by averaging replicates.
• plasma (pgmL): Plasma concentration of cortisol in pg/mL, calculated by multiplying the EIA avg conc by the dilution factor.
• plasma (ngmL): Plasma concentrations of cortisol in ng/mL, calculated by multiplying the plasma (pg/mL) by 0.0001.

File: StudyDesign.csv

Description: contains all relevant information about the samples selected for RNASeq analysis.

Variables

• Sample ID: Name given to the sample at time of collection.
• Tank ID: An alphanumeric designation given to the tank from which the frog was sampled.
• Frog ID: Identity of individual from which sample was collected.
• Sex: Sex of the sampled individual.
• Group: Behavioral group to which individual was assigned, either control (Egg Care (i.e., males (N=5) and females (N=5) performing egg care only, sampled 9 days post-oviposition), Sex-Typical (i.e., male transporting (N=5) and female observing (N=5), sampled 0–1 days post-hatching), Sex-Reversed-Flexible (i.e., female transporting in absence of male (N=5), sampled 0–1 days post-hatching), or Sex-Reversed-Inflexible (i.e., non-transporting female in absence of male (N=5), sampled 3 days post-hatching).
• Behavior: Behavior of focal individual at the time of sample collection.
• State: Social state of the pair, either Sex-Typical (both male and female present) or Sex-Reversed (male removed).
• Set: Sampling set to which sample contributed (1 to 6) to generate approximately time-matched behavioral groups.
• Lay Date: Date of oviposition for the focal clutch.
• Focal Clutch: Number of eggs and location (left (L) or right (R) egg canister) of the focal clutch
• Sampling Date: Date sample was collected.
• Sampling Day: Day sample was collected relative to the oviposition date (Sampling Date – Lay Date).
• Time Dissection Started: Time the dissection started.
• Time Dissection Finished: Time the dissection finished.
• Plasma Volume: Volume, in microliters, of total plasma obtained from trunk blood after spinning down red blood cells.
• Notes: Any notes recorded at the time of sampling, typically only recorded for one individual in the pair if two frogs were sampled simultaneously.
• RNA Extraction ID: Temporary identifier used to label tubes during RNA Extraction.
• Final Volume (uL): Total volume of eluent + RNA obtained following RNA Extraction.
• RNA Concentration (ng/uL): Concentration of RNA in extracted sample, as determined on a Qubit.

File: Trinity.gene.counts.matrix.txt

Description: is a matrix with 30 columns corresponding to the number of samples and 859,656 rows corresponding to the number of genes comprising the de novo transcriptome. Each is labeled according to the unique individual or gene identifier. Cells are populated by raw gene counts generated by Trinity.

File: trinotate_report.tsv

Description: is a file containing annotations for 1,229,709 transcripts produced via the Trinity-Trinotate pipeline. Trinotate makes use of a number of different well-referenced methods for functional annotation, each represented by the column headers in the output. For each field, ‘.’ denotes no hit for the transcript from the respective database.

Variables

• X.gene_id: The identifier for the gene to which the transcript belongs (if gene-level information is available).
• transcript_id: The unique identifier for the transcript, often derived from the input transcriptome assembly (e.g., Trinity IDs).
• sprot_Top_BLASTX_hit: Details of the top hit from a BLASTX search against a protein database (e.g., SwissProt), including the accession, description, E-value, and percent identity.
• infernal: noncoding RNA Annotations
• prot_id: The identifier for the protein to which the transcript belongs (if protein-level information is available).
• prot_coords: The coordinates, in base pairs, of the protein-coding region within the parent transcript.
• sprot_Top_BLASTP_hit: Details of the top hit from a BLASTP search against a protein database, typically for predicted protein sequences from the ORFs.
• Pfam: Information about identified protein domains from the Pfam database, including domain ID, description, and E-value.
• SignalP: Indication of whether a signal peptide was predicted for the corresponding protein according to the SignalP database.
• TmHMM: Indication of whether transmembrane helices were predicted for the corresponding protein according to the TmHMM database.
• Eggnog: Functional categories and descriptions from the eggNOG database.
• Kegg: KEGG pathway annotations, linking transcripts to metabolic and signaling pathways.
• gene_ontology_BLASTX: Gene Ontology terms corresponding to BLASTX nucleotide hits, providing information on biological process, molecular function, and cellular component.
• gene_ontology_BLASTP: Gene Ontology terms corresponding to BLASTP preptide hits, providing information on biological process, molecular function, and cellular component.
• transcript: optional identifier
• peptide: optional identifier

File: Missing_Annotations_Blast2Go_Table.csv

Description: contains additional gene annotations obtained by searching unannotated sequences of interest (i.e., 923 transcripts differentially expressed between behavioral groups) against NCBI’s non-redundant vertebrate nucleotide database with an e-value cutoff of 1.0E-3 using BLAST2GO software. Column headers were modified from the raw Blast2Go output to match Trinotate column headers. See above for definitions of specific fields.

Variables

• gene_id: The identifier for the gene to which the transcript belongs (if gene-level information is available).
• transcript_id: The unique identifier for the transcript, often derived from the input transcriptome assembly (e.g., Trinity IDs).
• sprot_Top_BLASTX_hit: Details of the top hit from a BLASTX search against a protein database (e.g., SwissProt), including the accession, description, E-value, and percent identity.
• infernal: noncoding RNA Annotations
• prot_id: The identifier for the protein to which the transcript belongs (if protein-level information is available).
• prot_coords: The coordinates, in base pairs, of the protein-coding region within the parent transcript.
• sprot_Top_BLASTP_hit: Details of the top hit from a BLASTP search against a protein database, typically for predicted protein sequences from the ORFs.
• Pfam: Information about identified protein domains from the Pfam database, including domain ID, description, and E-value.
• SignalP: Indication of whether a signal peptide was predicted for the corresponding protein according to the SignalP database.
• TmHMM: Indication of whether transmembrane helices were predicted for the corresponding protein according to the TmHMM database.
• Eggnog: Functional categories and descriptions from the eggNOG database.
• Kegg: KEGG pathway annotations, linking transcripts to metabolic and signaling pathways.
• gene_ontology_BLASTX: Gene Ontology terms corresponding to BLASTX nucleotide hits, providing information on biological process, molecular function, and cellular component.
• gene_ontology_BLASTP: Gene Ontology terms corresponding to BLASTP preptide hits, providing information on biological process, molecular function, and cellular component.
• transcript: optional identifier
• peptide: optional identifier

Code/software

All statistical analyses were performed in R Studio v. 2024.04.2+764 (R Core Development Team, 2020). Analyses are divided into three sections according to the type of data: Behavioral Analysis, Hormone Analysis, and Gene Expression Analysis. These are described generally below, with step-by-step annotations included within the scripts.

Behavioral Analysis: This script is divided into three subsections with header lines demarcated using ### HEADER ####

Male Call Rates Across Parenting Cycle and Between Social Conditions (beginning on L15): We first read in the male call rate data and summarize the mean and standard error for call rates by trial type and stage for the purpose of visualization. Raw data (call rates per hour) are then converted to integer form and differences between trial types and stages are modeled using a generalized linear mixed effects model with a poisson link to evaluate significance. Post-hoc comparisons are conducted using the emmeans package.

Sex and Treatment Differences in Transport Success (beginning on L46): We first read in the scoring info and summarize the proportion of successes and failures for each trial type for the purpose of visualization. Raw data (Tad Transported, Y/N) are then converted to binary format (1/0) and differences between trial categories (sex typical, or male transporting, and sex reversed, or female transporting) and trial types are modeled using bayeseian mixed effects models with a bernoulli link to evaluate specific hypotheses. Order effects are incorporated into the models by creating an order variable for each individual frog based on trial type. For specific pairwise contrasts, we concluded if the 95% credible interval (CI) of the difference contained zero, there was no evidence for a difference in transport success. Individual repeatability of transport success is assessed with the function rptBinary in the package rptR for all subjects and for females separately.

Sex Differences in Transport Performance or Prehatching Behavior (beginning L 211): Working from the previous dataset, we first create a subsetted dataset including only trials where successful transport was observed and evaluate sex differences in latency to transport and time in transit by fitting survival curves and gamma-distributed generalized linear mixed models, respectively. We next read in the behavioral scoring data from the last full day before hatching and calculate the duration of each instance of a prehatching behavior. These data are then summarized to give the total time spent on each behavior and the number of times the behavior was observed over the approximately 12-hour period, and null data (behaviors that were not observed for a focal individual) are imputed to account for zeros. Individual repeatability of prehatching behaviors is assessed with either the rptPoisson or the rptGaussian function for all subjects. Next, variation in prehatching behaviors across sexes and trial types is visualized as violin plots and modeled as either a linear mixed effects model or a generalized linear mixed effects model with poisson links. Finally, to explore the relationship between prehatching and posthatching parental behaviors, we use boxplots and binomial generalized linear mixed models to assess the predictive power of prehatching behaviors on female transport success and use linear mixed effects model or a generalized linear mixed effects model with poisson links to assess the predictive power of prehatching behaviors on latency to transport and tadpole transport time in trials where transport occurred.

This script requires the following packages: plyr, dplyr, stringr, survival, ggplot2, survminer, brms, tidybayes, lme4, lmerTest, emmeans, cowplot, lubridate, rptR, and data.table

Hormone Analysis: This script reads in the hormone data, plots it as a function of behavioral group, and statistically compares cortisol concentrations between behavioral groups using an ANOVA followed by Tukey post hoc contrasts.

This script requires the following packages: tidyverse, multcomp

Gene Expression Analysis: This script is divided into three subsections with header lines demarcated using ### HEADER ####
Prepare Datasets (beginning L23): We first read in the Trinity gene count matrix, associated gene annotations, and study design, and clean up headers, collapse duplicate annotations, and create sub-datasets in preparation for differential expression analysis.
Differential Expression Analysis (beginning L84): We next carry out the differential expression (DE) analysis on unnormalized, raw counts using DeSeq2. The first model considers all adults and compares gene expression between sexes and behavioral groups. This DeSeq output is used for PCA analysis and visualization of the separation of sexes and behavioral groups (Figure 5). Subsequently, separate DE analyses are performed for males and females to compare gene expression of parenting individuals to sex-matched controls (egg care). Each DE analysis includes a filtering step for low count genes, excluding genes that reported fewer than 10 transcripts in at least half of samples. Genes DE at alpha<0.1 in any comparison are saved to tables for functional enrichment analysis.

Functional Enrichment Analysis (beginning L301): Finally, we perform functional enrichment analysis of differentially expressed genes using the TopGO package. Briefly, GO terms annotated for each gene in the gene annotations table are consolidated across various database hits and extracted for genes of interest (i.e., genes differentially expressed at alpha<0.1 in tested contrasts, above). We then follow the TopGO workflow to perform functional enrichment analysis of various gene sets for GO terms belonging to each of the three gene ontology categories: biological process, molecular function, and cellular component. We retain GO terms that are significantly enriched at p<0.01. The final step is to combine enriched GO terms from each of the three categories in a single table for each gene set and sort them by significance.

This script requires the following packages: plyr, dplyr, data.table, magrittr, qdapRegex, DESeq2, pcaExplorer, factoextra, VennDiagram, RRHO, RRHO2, topGO, GOplot, genefilter, devtools, smart, and ggalt