The detailed methods of the two experiments have been published in two open-access papers (Bókony et al., 2020; Mikó et al., 2021a). Here we provide a brief description of each experiment, and present detailed methods only for those aspects that are directly relevant for the current study. All experimental procedures were approved by the Ethical Commission of the Plant Protection Institute and carried out according to the permits issued by the Government Agency of Pest County (Department of Environmental Protection and Nature Conservation) and the Budapest Metropolitan Municipality (Department of City Administration, FPH061/2472-4/2017).

The study design is illustrated in Figure 1. For both experiments, we collected freshly spawned agile frog eggs in March 2018 from ponds in woodland habitats in north-central Hungary. We sampled 8 egg masses from each of three ponds for experiment 1, and further 10 egg masses from one of these ponds for experiment 2 (Supplementary Table 1). The eggs were taken into our laboratory, where the two experiments were conducted simultaneously in the same room under artificial photoperiod that mimicked the outdoors dark-light cycles. We raised eggs and tadpoles in reconstituted soft water (RSW) at 19°C water temperature. We started both experiments when the hatchlings reached development stage 25 (Gosner, 1960) by placing haphazardly selected tadpoles into white plastic boxes (14 × 9.5 cm base area) filled with 1 L RSW. Each animal was kept in a separate, individually labeled box throughout the entire experiment. The boxes were arranged in a randomized block design to ensure that all treatments were homogeneously distributed across the shelves in the laboratory.

In both experiments, we exposed the tadpoles to environmentally relevant concentrations of water-polluting chemicals for which sex-related endocrine-disrupting effects had been reported. In experiment 1, we applied two concentrations each of carbamazepine (0.5 and 50 μg/L), a pharmaceutical drug (Galus et al., 2013, 2014), and terbuthylazine (0.003 and 0.3 μg/L), a herbicide frequently used in Europe (Kjeldsen et al., 2013). In experiment 2, we applied two concentrations of chlorpyrifos (0.5 and 5 μg/L), a broad-spectrum insecticide (Bernabò et al., 2011). The applied two concentrations, respectively, for each chemical correspond to the mean (or median) and close-to-maximum values reported from surface waters (Bókony et al., 2020; Mikó et al., 2021a). In both experiments, the control group of tadpoles was kept in clean RSW to which we added ethanol as solvent control (1 μL 96% ethanol to 1 L RSW); all other treatment groups also contained this amount of ethanol as vehicle. In experiment 1, we distributed 480 tadpoles (20 from each family) evenly across five treatment groups with four replicates in each treatment × family combination (4 tadpoles × 24 families × 5 treatments). In experiment 2, we distributed 144 tadpoles evenly across three treatment groups (48 tadpoles per treatment, with 4–6 tadpoles in each treatment × family combination).

We exposed tadpoles to the treatments over the entire duration of their larval development. Twice a week, we renewed each treatment by changing the rearing water and fed the tadpoles ad libitum with chopped, slightly boiled commercial spinach. We collected different kinds of data on tadpole behavior in the two studies (see below), with our focus being on foraging activity in experiment 1 and anti-predatory response in experiment 2. When the tadpoles reached the start of metamorphic climax (appearance of forelimbs at development stage 42), we recorded this date and measured their body mass (±0.1 mg). We raised the post-metamorphic animals in individual rearing boxes, feeding them ad libitum with small crickets until they reached the minimum age that allows reliable identification of phenotypic sex by the gross anatomy of the gonads (Ogielska and Kotusz, 2004; Bernabò et al., 2011). On average 15 weeks (96–136 days) after starting the experiment, we weighed the animals to the nearest 0.01 g and dissected them for phenotypic sexing (see below). During dissection, we removed the entire digestive tract and measured its mass; this value was subtracted from froglet body mass to remove any variance due to food remains in the gut.

In experiment 1, we observed the behavior of each tadpole in their rearing boxes containing ad libitum food, using the “instantaneous sampling” method (Altmann, 1974) four times a week during the larval period, totaling 20 observations per tadpole. Four researchers conducted the observations on 2 days each week, in one morning session and one afternoon session (each lasting ca. one hour) each day. During each session, we scanned all tadpoles once in a fixed order, recorded their instantaneous behavior as inactive, feeding, or swimming, and we also categorized the location of the tadpole within the box as on the bottom, next to the wall (i.e., within one tadpole distance from the wall, but not on the bottom), or in the open (i.e., not near the bottom or wall). In each session, a single observer scanned all tadpoles, and the identity of observers was rotated between sessions.

In experiment 2, we video-recorded tadpole behavior 1 month after starting the experiment and 3 days thereafter. As we had a limited number of cameras, we video-recorded only a subsample (83%) of individuals (40 from each chemical treatment group). On each occasion, we transferred the tadpoles into new containers identical to their rearing boxes, but with no food (to facilitate automatic tracking). On the first occasion, the box contained RSW with the same chemical treatment the individual was reared in, whereas on the second occasion all tadpoles were moved into clean RSW for video recording. Each time, we recorded the tadpoles’ behavior for 20 min, and then exposed them to a startling stimulus by abruptly pouring 40 ml RSW into their water. Typically, the tadpoles reacted by a short burst of swimming (“escape”), followed by a period of motionlessness (“freezing”). To allow for measuring these responses we continued recording for another 20 min after the stimulus. As a main objective of experiment 2 was to assess anti-predator responses, half of the tadpoles within each chemical treatment group received clean RSW as startling stimulus on both occasions, whereas the other half received chemical cues indicating predation risk. The chemical cues were prepared as described in Hettyey et al. (2016), using water from the tank of European perch (Perca fluviatilis) that had been feeding on agile frog tadpoles. As the European perch is a native predator in our region, agile frog tadpoles respond to chemical cues indicating its presence by decreasing their activity, even if they are predator-naïve (Hettyey et al., 2016).

From the video recordings, we collected data on tadpole activity using the automatic tracking software ToxTrac (Rodriguez et al., 2018). All tracking results were manually checked for quality; all tracking data used in the current paper were error-free. We calculated the following variables from the first 20 min of each recording (i.e., before the addition of predator cues): total distance moved as a measure of locomotor activity, proportion of area used as a measure of exploration rate, and proportion of time spent near the wall (within a 50-pixel wide stretch from the side of the box; the tadpoles had an average snout-to-vent length of ca. 40 pixels in the videos) as a measure of risk aversion. For quantifying the startle response, the second 20 min of each video recording (i.e., after the addition of predator cues) were analyzed manually by a single observer, who recorded the following variables: the intensity of immediate reaction to the startle (“startle response”), categorized on a 0–3 scale (0: no movement, 1: slight movement, 2: swimming away apparently calmly, 3: swimming around fervently); the duration of escape, measured as the time spent moving continuously from the startle stimulus until the first stop; and the duration of freezing, measured as the time spent motionless after the escape until the first movement thereafter. Note that the duration of escape and freezing were quantifiable only in those individuals whose startle response had a non-zero intensity score.

We euthanized the animals by immersion into 6.66 g/L MS-222 solution buffered to neutral pH. After dissection, we examined the gonads under an Olympus SZX12 stereomicroscope at 16 × magnification, and categorized phenotypic sex as male (testes), female (ovaries), or uncertain. There were only two individuals in the latter category, and one of them also had uncertain genetic sex (see below); we excluded this animal from all statistics presented in this paper. Due to a low level of mortality in both studies, we had data on phenotypic sex for 439 individuals in experiment 1 and 135 individuals in experiment 2 (Table 1).

Throughout both experiments we monitored survival daily, and stored a tissue sample from each animal that died before dissection, as well as from all dissected animals, in 96% ethanol. We extracted DNA using Geneaid Genomic DNA Extraction Kit for animal tissue (Thermo Fisher Scientific), following the manufacturer’s protocol, except that digestion time was 2 h. We identified genetic sex of all animals, and sex reversal in phenotypically sexed individuals, using a recently developed molecular marker set (Nemesházi et al., 2020) which has been validated for agile frog populations in our study area, including the three populations sampled for the present study. In short, first we tested all froglets for marker Rds3 (≥95% sex linkage; primers: Rds3-HRM-F and Rds3-HRM-R) using high-resolution melting (HRM), and we accepted an individual to be concordant male or concordant female if its Rds3 genotype was in accordance with its phenotypic sex. Those individuals that appeared sex-reversed by Rds3 were also tested for marker Rds1 (≥89% sex linkage; primers: Rds1-F, Rds1-R and Rds1-Y-R) using PCR, and were accepted to be sex-reversed only if both markers confirmed sex reversal. Those individuals that were not phenotypically sexed due to early mortality were screened for both Rds1 and Rds3. Individuals with discrepant genotyping results (i.e., contradiction between Rds1 and Rds3) were considered to be of unknown genetic sex. This approach yielded data on genetic sex for 465 individuals in experiment 1 and 140 individuals in experiment 2 (Table 1).

For those individuals that we identified as sex-reversed, we also examined the gonads histologically to make sure that the mismatch was not due to incorrect categorization of phenotypic sex. We also investigated gonad histology in the individual for which phenotypic sex was uncertain based on gross anatomy but the genetic sex was unambiguous. We fixed the dissected gonads (not separated from the kidneys, because of small gonad size) in neutral-buffered 10% formalin, and prepared histological sections as described in our earlier papers (Nemesházi et al., 2020; Mikó et al., 2021b). For one individual, the sections failed to include gonadal tissue. For all other individuals examined, the gonads were clearly identifiable by histology either as testes (n = 14) or ovotestes (testes containing a few oogonia, see Supplementary Figure 1; n = 6, including the individual whose phenotypic sex had been uncertain based on gross anatomy). Therefore, in the analyses we treated these individuals as phenotypic males (female-to-male sex-reversed individuals).

An overview of all our analyses is given in Supplementary Table 2. First, we tested whether sex-reversal rate was independent of chemical treatment by analyzing the phenotypic sex of genetic females in a generalized linear mixed-effects model with binomial error and logit link, including treatment type as fixed factor, and family nested in experiment as random factors. Two treatment groups lacking sex-reversed individuals had to be excluded from this analysis because such separation in logistic models results in unreliable estimates; thus, we used Fisher’s exact tests to compare sex-reversal rate in these two groups with the respective control group.

We tested whether survival rate differed between genetic males and genetic females using a Cox’s proportional hazards model, including family nested in experiment as random effects, and treating the dissected individuals as censored observations (i.e., these animals survived until the end of the study). In all remaining models (see below), we included sex as a three-category factor (concordant male, concordant female, or sex-reversed), thus individuals with missing data on either genetic or phenotypic sex were excluded.

We analyzed the duration of larval development (number of days from starting the experiment at stage 25 until the start of metamorphosis at stage 42), body mass at metamorphosis, and juvenile body mass (measured at dissection ca. 2 months after metamorphosis; excluding gut mass) by pooling the data of the two experiments, taking into account chemical treatment as a fixed factor and family nested in experiment as random factors. Duration of larval development was analyzed with a Cox’s proportional hazards model, whereas both mass variables were analyzed with linear mixed-effects models, allowing for heteroscedasticity among the three sex categories. In the model of juvenile body mass, we also included age (number of days from starting the experiment until dissection) as a covariate. To investigate the trade-off between development and growth in tadpoles in the three sex categories, we added the interaction between the duration of larval development and sex into the model of mass at metamorphosis; in this model we allowed for heteroscedasticity also among families because residual diagnostics indicated that a few families exhibited outliers in the relationship between development and growth.

We analyzed tadpole behavior separately for the two experiments. For experiment 1, we analyzed two variables. We compared the proportion of observations in the open among the three sex categories with a Fisher’s exact test, because separation in the data precluded the use of logistic models (thus, in this statistical test we could not control for potential effects of other predictors such as chemical treatment). To analyze the proportion of individuals feeding, we used a generalized linear mixed-effects model with binomial error and logit link. For experiment 2, we analyzed the pre-startle distance moved, exploration rate, and time spent near the wall with linear-mixed effects models, allowing for heteroscedasticity among the three sex categories. We analyzed the intensity of startle response using a cumulative-link mixed model with logit link function, and the durations of escape and of freezing with Cox’s proportional hazards models, excluding those tadpoles that did not react with movement. All models of behavioral variables included individual nested in family as random factors, chemical treatment as a fixed factor, and the fixed effects of date (a covariate in experiment 1, expressing the number of days from starting the experiment; and a two-category factor in experiment 2) and time of day (a two-category factor in experiment 1, and a covariate in experiment 2, expressing the order of video recordings which were done in 7 consecutive bouts). Additionally, the model of experiment 1 included shelf height and observer identity as fixed factors and water temperature as a covariate (for detailed explanation on these covariates, see Bókony et al., 2020). To investigate the change of feeding rate over time in the three sex categories, we added the interaction between date and sex into the model of feeding rate (we expected that young tadpoles with undifferentiated gonads would behave similarly, whereas feeding rate may diverge between sexes in later stages as the gonads become differentiated). For experiment 2, the model of exploration rate included the total distance moved as a covariate, because we aimed to investigate the percentage of area used independently of locomotor activity. The three models of post-startle variables also included the total distance moved as a covariate, because post-startle activity may depend on pre-startle activity. Furthermore, these latter three models included stimulus type (i.e., whether the stimulus water contained predator cues or not) as a fixed factor, and its interaction with sex to test whether the effect of predator cues on behavior differed between sex categories.

In each model of development, growth, and behavior, we tested the effect of sex by pre-planned comparisons (Ruxton and Beauchamp, 2008). Specifically, we extracted three linear contrasts from each model: genetic males versus genetic females (the latter including sex-reversed individuals), phenotypic females versus phenotypic males (the latter including sex-reversed individuals), and sex-reversed versus concordant individuals (the latter including males and females). We provide these contrast estimates with 95% confidence intervals (CI) as non-standardized measures of effect size (Nakagawa, 2004; Nakagawa and Cuthill, 2007); we interpret CIs excluding zero (or one, in case of odds ratios and hazard ratios) as statistically significant. Further, to test whether the sexes differed in the trade-off between development and growth, we estimated the slope of relationship between the duration of larval development and body mass at metamorphosis for each sex category, and then applied the above three linear contrasts to the slopes. We used the same approach to compare the slope of change over time in feeding rate among the sexes. To test whether the sexes differed in the effect of predator cues on post-startle behaviors, for each sex category we estimated the difference in each behavioral variable between animals startled with versus without predator cues, and again we applied the above three linear contrasts to these predator-effect estimates.

All analyses were performed in the R computing environment v4.0.3 (R Core Team., 2020). Although our sample sizes are unbalanced due to the small number of sex-reversed individuals, we used mixed models throughout, which provide a flexible and powerful tool for appropriately analyzing unbalanced and heteroscedastic data (Pinheiro and Bates, 2000). For each analysis, we checked that the statistical requirements of the model were met by visually inspecting relevant graphs of residuals. Our data with a detailed, annotated R script are available as given in the Data Availability Statement.

In the total sample, there was significant female bias in genetic sex ratio (344/605 = 56.9% females, 95% CI = 0.53–0.61) but not in phenotypic sex ratio (305/574 = 53.1% females, 95% CI = 0.49–0.57). Out of the 571 individuals for which both genetic and phenotypic sex was identifiable, 21 were female-to-male sex-reversed (Table 1), 16 of which originated from the pond sampled for both experiments (9.8% sex-reversal rate), and 5 from the two ponds sampled only for experiment 1 (2.4 and 3.1% sex-reversal rate, respectively). The sex-reversed individuals came from 8 (out of 34) different families (5–75% of individuals sex-reversed per family; Supplementary Table 1), which does not conform to a homogeneous distribution of sex reversal among families (χ² = 233, df = 33, p < 0.001). In 3 out of the 34 families we found no genetic males at all; two of these families exhibited 30 and 75% sex-reversal rate, respectively (Supplementary Table 1). The frequency of sex reversal among genetic females was independent of chemical treatments (6.4% overall; Table 1). Survival rate until dissection did not depend on genetic sex (18 females and 16 males died; hazard ratio: 0.84, 95% CI: 0.60–2.38).

We found no significant difference between any combination of genetic and phenotypic sex in the length of larval development and body mass at metamorphosis or at dissection (Table 2 and Figure 2). However, the relationship between the length of larval development and mass at metamorphosis varied significantly with sex (Table 2 and Figure 3A): animals that metamorphosed later had higher body mass in both concordant males and concordant females (Table 3), whereas in sex-reversed individuals the slope of this relationship did not differ significantly from zero (Table 3). Some sex-reversed individuals metamorphosed relatively early and with large mass (Figure 3A); most of them came from those families where we detected no genetic males at all (Supplementary Figure 2). Other sex-reversed individuals metamorphosed relatively late and with small mass (Figure 3A); most of them had testicular oogonia (Supplementary Figure 2).

In experiment 1, none of the 5 sex-reversed individuals were ever observed in the open; 0.37 and 0.39% of observations of concordant females and concordant males, respectively, were in the open. This difference was not significant (Fisher’s exact test: p = 0.900). Also, feeding frequency did not vary significantly with sex (Table 2). The proportion of tadpoles feeding increased slightly as the tadpoles aged and this increase appeared greatest in sex-reversed individuals and smallest in concordant males (Figure 3B); however, the slopes had wide confidence intervals (Table 3) and did not differ significantly between any combination of genetic and phenotypic sex (Table 2).

In experiment 2, the total distance moved was significantly shorter in sex-reversed tadpoles than in concordant individuals (Table 2 and Figure 2), but exploration rate and time spent near the wall did not vary significantly with sex (Table 2 and Figure 2). Similarly, we found no significant differences in the intensity of startle response and the durations of escape and freezing between any combination of genetic and phenotypic sex (Table 2). However, the presence of predator cues modified these behaviors in sex-dependent ways (Table 3). Among concordant females, those that received predator cues responded less intensely to the disturbance than those that received clean water (Table 3 and Figure 4); there was no such difference in concordant males or in sex-reversed individuals (Table 3 and Figure 4). The duration of escape was shorter in sex-reversed individuals if they received predator cues than when they did not (Table 3 and Figure 5); there was no such difference in concordant individuals (Table 3 and Figure 5). The duration of freezing after the escape reaction was longer in concordant females in the presence of predator cues than in clean water (Table 3 and Figure 5); the similar trends in concordant males and sex-reversed individuals were not statistically significant (Table 3 and Figure 5). Out of all these sex differences in the effects of predator cues, only the one for escape duration was statistically significant (Table 2).