Sex Reversal and Performance in Fitness-Related Traits During Early Life in Agile Frogs

Sex reversal is a mismatch between genetic sex (sex chromosomes) and phenotypic sex (reproductive organs and secondary sexual traits). It can be induced in various ectothermic vertebrates by environmental perturbations, such as extreme temperatures or chemical pollution, experienced during embryonic or larval development. Theoretical studies and recent empirical evidence suggest that sex reversal may be widespread in nature and may impact individual fitness and population dynamics. So far, however, little is known about the performance of sex-reversed individuals in fitness-related traits compared to conspecifics whose phenotypic sex is concordant with their genetic sex. Using a novel molecular marker set for diagnosing genetic sex in agile frogs (Rana dalmatina), we investigated fitness-related traits in larvae and juveniles that underwent spontaneous female-to-male sex reversal in the laboratory. We found only a few differences in early life growth, development, and larval behavior between sex-reversed and sex-concordant individuals, and altogether these differences did not clearly support either higher or lower fitness prospects for sex-reversed individuals. Putting these results together with earlier findings suggesting that sex reversal triggered by heat stress may be associated with low fitness in agile frogs, we propose the hypothesis that the fitness consequences of sex reversal may depend on its etiology.


INTRODUCTION
Sex is a fundamental aspect of individual state in all sexually reproducing organisms. Having testes versus ovaries often comes with a diverse set of differences in physiology, morphology, life history, and behavior, including mating and parental strategies often labeled as "sex roles" (Schärer et al., 2012;Immonen et al., 2018). In species with genetic sex determination, where the process of gonad development is triggered by genomic elements, males and females also often differ in their genetic make-up. For example, the chromosome restricted to one sex (e.g., Y in male-heterogametic systems) is inclined to undergo degeneration, which may lead to sex differences in mortality rates and senescence (Marais et al., 2018). In species with environmental sex determination, where the fate of the gonads is decided by external factors such as temperature during early ontogeny, sex ratios and hence population viability may be particularly vulnerable to environmental changes (Mitchell and Janzen, 2010). Thus, the way males and females come to be has crucial implications for population dynamics and thereby biodiversity conservation.
In vertebrate animals, genetic and environmental sex determination have traditionally been thought of as mutually exclusive systems, but research on ectotherms, especially in the past decade, has revealed a growing number of species in which the two systems naturally coexist (Baroiller and D'Cotta, 2016;Holleley et al., 2016). In such species, environmental influences can override the effect of genes during sex determination in early life, resulting in sex reversal whereby individuals develop phenotypic sex discordant with their genetic sex (Alho et al., 2010;Baroiller and D'Cotta, 2016;Holleley et al., 2016;Lambert et al., 2019;Nemesházi et al., 2020). Theoretical studies predict that sex reversal has far-reaching consequences for demography, population persistence, sex chromosome loss and evolutionary transitions between sex-determination systems (Grossen et al., 2011;Quinn et al., 2011;Bókony et al., 2017;Wedekind, 2017;Schwanz et al., 2020;Nemesházi et al., 2021). Studying these consequences empirically is especially important in light of the ongoing rapid human-induced environmental alterations, including climate change and chemical pollution, which may increase sex-reversal frequency in nature .
One of the most urgent questions regarding sex reversal is how sex-reversed individuals compare to concordant males and females in terms of performance in fitness-related traits. Understanding this issue would be highly valuable for at least three reasons. First, it would facilitate forecasting the effects of sex reversal on demography and evolution, because many of these theoretically predicted effects critically depend on the viability and reproductive success of sex-reversed individuals (Grossen et al., 2011;Bókony et al., 2017;Nemesházi et al., 2021). Second, it would provide insight into the ultimate and/or proximate drivers of sex reversal. On the one hand, sex reversal may be an adaptive sex-allocation strategy that allows individuals to develop the sex that is most beneficial under the prevailing environmental conditions (Geffroy and Douhard, 2019), similarly to temperature-dependent sex determination (Schwanz et al., 2016). In this case, sex-reversed individuals may perform at least as well or even better than concordant individuals, as found for fecundity in a reptile (Holleley et al., 2015). On the other hand, sex reversal may be a mechanistic consequence of endocrine disruption due to early life stress, as demonstrated in fish (Senior et al., 2012;Baroiller and D'Cotta, 2016). In this case, sex-reversed individuals may display reduced performance in important life-history traits due to early life stress, as seems to be the case in an anuran amphibian Mikó et al., 2021b). Yet another alternative is that sex reversal may arise by random sex determination (Perrin, 2016): in absence of strong genetic and environmental triggers, sex may be determined by stochastic variability in gene expression levels. This process may coexist with genetic and environmental sex determination, and can explain considerable proportion of phenotypic variance (Perrin, 2016). Thus, random variation (i.e., no systematic difference) between sex-reversed and concordant individuals in fitnessrelated performance might indicate random sex determination. Thirdly, because sex reversal de-couples genetic and phenotypic sex, it allows for evaluating the relative importance of sex-linked genes versus gonadal effects (sex hormones and other sex-specific modifiers that orchestrate sex-biased gene expression) in the development of sex-specific life histories and behaviors. Despite all these reasons for studying the consequences of sex reversal, we have very little empirical information on the fitness of sexreversed individuals, apart from fish in aquaculture (Senior et al., 2012), where sex reversal is artificially induced and thus may not necessarily be ecologically relevant. Researchers have only just begun to investigate the relationship between ecologically relevant sex reversal and individual performance, and most of the little existing knowledge comes from a single reptilian species, where high incubation temperatures produce male-to-female sex-reversed individuals that display a complex combination of male-like, female-like, "supermale" and "superfemale" traits (Holleley et al., 2015;Li et al., 2016;Jones et al., 2020).
Here we address the fitness-related performance of sexreversed individuals by using data from two previous experiments Mikó et al., 2021a) on agile frogs (Rana dalmatina), a European species whose natural populations contain a considerable number of female-to-male sex-reversed adults . Both experiments were designed to test sub-lethal effects of larval exposure to environmentally relevant concentrations of chemical pollutants on early life traits related to fitness. Neither of the chemical treatments affected phenotypic sex ratios significantly Mikó et al., 2021a), but as we report here, sexreversed individuals occurred in both experiments, indicating that some agile frog tadpoles spontaneously undergo sex reversal even in the absence of any sex-reversing chemical or thermal treatment, in accordance with another study (Mikó et al., 2021b). Therefore, these data offered us the opportunity to examine whether genetic and phenotypic sex and the combination thereof (i.e., reversed vs. concordant sex) are associated with differences in life history and behavior during early ontogeny. The larval phase is of critical importance in amphibian life history because mortality during this stage can be extremely high (Riis, 1991), mainly due to predation, pond desiccation, and limited food availability. Thus, larval survival depends to a large extent on the behavioral strategies (predator avoidance, foraging activity) adopted by tadpoles (Skelly, 1994) and the speed of their development (Griffiths, 1997). Also, the rates of development and growth until metamorphosis can have life-long effects on fitness in amphibians (Smith, 1987;Berven, 1990;Altwegg and Reyer, 2003). However, sex differences in larvae are very rarely investigated due to the difficulties of phenotypic sex identification in immature animals  and genetic sexing in amphibians overall . In the agile frog, males typically start reproducing 1 year earlier than females (Riis, 1991;Sarasola-Puente et al., 2011), and larger males are more successful in male-male competition (Vági and Hettyey, 2016). Therefore, we predicted to observe higher growth rate and faster larval development in males. Fast growth and development require high food intake, so we predicted that male tadpoles would spend more time feeding and, in trade-off, take higher predation risk (Urszán et al., 2015). At least some of these sex differences might be sex-chromosome-linked, since the agile frog has an XX/XY sex-chromosome system (Jeffries et al., 2018); this would make the female-to-male sex-reversed individuals resemble concordant females. However, the agile frog sex chromosomes are homomorphic (Jeffries et al., 2018), suggesting limited genetic differentiation between the sexes. By this latter logic, female-to-male sex-reversed individuals may be more likely to resemble concordant males, due to the presence of testes which produce androgen hormones that stimulate the expression of male phenotypic traits (Guarino and Bellini, 1993). We evaluated these predictions by comparing early life development, growth, and behavior among three groups: males and females with concordant genetic and phenotypic sex, and female-to-male sex-reversed individuals.

Experimental Procedures
The detailed methods of the two experiments have been published in two open-access papers 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(Galus et al., , 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 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.

Behavioral Observations and Video Analysis
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 antipredator 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 predatornaïve .
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.

Phenotypic and Genetic Sexing
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 The table excludes 15 individuals for which we had no data on sex (7 escaped before dissection and DNA sampling; 8 died before dissection and could not be sexed genetically). *Phenotypic sex could not be diagnosed for the individuals that died before dissection. † Three phenotypic males had unknown genetic sex due to marker disagreement (Rds3 genotype was female whereas Rds1 genotype was male). ‡ These two odds ratios were taken from Fisher's exact tests; the rest from a binomial mixed model (overall effect of treatment: Wald test, χ 2 = 1.24, df = 4, P = 0.871).
phenotypically sexed individuals, using a recently developed molecular marker set  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 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).

Statistical Analyses
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 threecategory factor (concordant male, concordant female, or sexreversed), 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 prestartle 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.

Sex Ratios and Sex Reversal
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 sexreversed ( 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 (χ 2 = 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).

Development and Growth
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 Significant differences (i.e., CI excluding 0 or, in case of odds ratios and hazard ratios, 1) are marked with bold text. The comparisons marked with an asterisk are calculated from the effect sizes given in Table 3. § Concordant individuals include XX females and XY males. † Genetic females include sex-reversed individuals and concordant females. ‡ Phenotypic males include sex-reversed individuals and concordant males.
males and concordant females (Table 3), whereas in sexreversed 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 sexreversed individuals metamorphosed relatively late and with small mass ( Figure 3A); most of them had testicular oogonia (Supplementary Figure 2).

Tadpole Behavior
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).

DISCUSSION
In this study, we confirmed that female-to-male sex reversal occurs in agile frogs at a relatively low frequency (6.4%) in the absence of thermal stress, and demonstrated that it was independent of chemical treatments representing ecologically relevant concentrations of carbamazepine, terbuthylazine, and chlorpyrifos. We are confident that these mismatches between genetic and phenotypic sex were indeed sex reversals, because our phenotypic sexing was backed up by histological analysis, and our genetic sexing method is based on two sex-linked markers that are located relatively far from each other on the sex chromosomes. The latter makes it highly unlikely that we would misdiagnose rare mutation or recombination events as sex reversal, which can FIGURE 2 | Development, growth, and behavior before the addition of predator cues in sex-reversed individuals (XX males), concordant (XX) females and concordant (XY) males. In each box plot, the thick middle line and the box represent the median and interquartile range, respectively; whiskers extend to the most extreme data points within 1.5 × interquartile range from the box. happen when only one marker is used (Toli et al., 2016). Because the animals were raised in the laboratory at benign temperatures with ad libitum food and no predators, and their sex development was not altered by the chemical treatments, it seems likely that these instances of sex reversal occurred independently of any obvious environmental stressor. Furthermore, the sex-reversed individuals in this study were similar to their sex-concordant siblings in almost all morphological, life-history and behavioral traits that we examined. Taken together, these results may be explained by two, not mutually exclusive ideas. First, the theory of random sex determination (Perrin, 2016) postulates that, in the lack of strong genetic and environmental effects on sex, developmental noise (i.e., random fluctuations in the expression of sex-determining genes) decides the sexual fate of individuals. Second, the threshold model of sex determination (Quinn et al., 2011;Nemesházi et al., 2021) assumes that phenotypic sex depends on whether the amount of "male signal" (i.e., expression of male-producing developmental signals, which can be influenced by both genotype and environment) exceeds the individual's threshold for male development, a trait encoded by FIGURE 3 | Relationships between mass at metamorphosis and duration of larval development (A) and between feeding rate and tadpole age (B) in sex-reversed individuals (XX males; purple, filled circles), concordant females (XX females; red, open circles), and concordant males (XY males; blue, triangles). The lines and surrounding polygons, respectively, represent the slopes and their confidence intervals given in Table 3. To facilitate graph readability, panel A was cropped to exclude an outlier point (a concordant male with 373.4 mg mass at metamorphosis and 78 days of larval development; note that this individual was included in all analyses and also in the estimation of the curve fitted here). Significant effects (i.e., CI excluding 0 or, in case of odds ratios and hazard ratios, 1) are marked with bold text.
genetic elements. Thus, for individuals who happen to have genes encoding high "male signal" levels and/or low threshold levels, even a small elevation of environmentally induced "male signal" expression may result in female-to-male sex reversal. This theory is supported in the present study by the non-random distribution of sex reversals among agile frog families, and by the high rate of sex reversal in those families where we detected only genetic females and no genetic males at all. The latter fits the threshold theory because the agile frog has an XX/XY sex determination system, so families containing 100% female offspring suggest that the sire in those families had been a female-to-male sexreversed individual (i.e., an XX male, mating with a concordant XX female) who may have passed on his alleles encoding high propensity for sex reversal to his offspring. Such a combination of genetic variation and random environmental noise might explain at least some occurrences of sex reversal in natural populations, especially where sex-reversal rate is not correlated with environmental factors such as the level of urbanization (Lambert et al., 2019), climate (Castelli et al., 2021), or elevation (Phillips et al., 2020).
The fact that the sex-reversed individuals in this study did not differ from concordant individuals in growth and development stands in stark contrast with our findings from another experiment on agile frogs (Mikó et al., 2021b), where sex reversal was induced by a six-days heat treatment during larval development. Heat treatment resulted in high rates of femaleto-male sex reversal, but also reduced survival, development, growth, and fat reserves (Mikó et al., 2021b). Thus, in that experiment, sex reversal was strongly associated with signs of developmental stress and poor fitness prospects, similarly to what has been reported about various fishes (Senior et al., 2012;Baroiller and D'Cotta, 2016). Combining those findings with our current results, we speculate that the fitness of sexreversed individuals may depend on the etiology of sex reversal. When it arises by stochastic variation in the biochemical processes of sex determination or in individual sensitivity to environmental effects on sex, it might not be systematically accompanied by changes in fitness-related traits. In contrast, when sex reversal is triggered by strong environmental effects and/or high physiological stress, it might be associated with FIGURE 4 | Intensity of response to the startling stimulus with and without predator cues in sex-reversed individuals (XX males), concordant (XX) females, and concordant (XY) males in experiment 2 (0: no movement, 1: slight movement, 2: swimming away apparently calmly, 3: swimming around fervently).
poor health or reduced performance in life-history traits. This association may arise by the same stressor affecting both sex and fitness-related traits, perhaps mediated by stress-induced glucocorticoid hormone effects (Geffroy and Douhard, 2019) or cellular calcium-redox regulation (Castelli et al., 2020). For example, a meta-analysis concluded that the poor fitness of fish that underwent chemically induced sex reversal was not due to sex reversal per se, but was the result of the chemical treatments themselves (Senior et al., 2012). Additionally, the association between sex reversal and fitness might be exacerbated by sex reversal itself directly affecting some fitness-related traits (Mikó et al., 2021b) or by making the offspring of sex-reversed individuals more sensitive to environmentally induced sex reversal (Piferrer and Anastasiadi, 2021). If one sex can do better than the other under stressful conditions, environment-induced sex reversal may serve as an adaptive sex-allocation strategy (Geffroy and Douhard, 2019). However, in agile frogs that spawn in early spring and develop in cool waters, high temperatures during larval development might not have been frequent enough in their evolutionary past for such an adaptive strategy to evolve. Nevertheless, recent findings indicate that sex-reversed agile frogs occur more frequently in anthropogenic habitats , and phenotypic sex ratios have become more malebiased in some amphibian species since the start of contemporary climate change (Bókony et al., 2017), suggesting that sex reversals might be shifting from mostly spontaneous or stochastic to increasingly stress-induced incidences. These speculations would deserve further empirical testing.
In the few instances where we found significant differences between sex-reversed and sex-concordant individuals in the present study, the former stood out by having lower locomotor activity and responding to disturbance with shorter escape duration when predator cues were present. Furthermore, the natural trade-off between larval development speed and growth rate seemed to be lacking in sex-reversed individuals. These findings support neither higher nor lower performance in terms of overall fitness for sex-reversed animals. First, while low activity may constrain foraging success, sex-reversed tadpoles were feeding at least as often as concordant individuals. Second, while short escape duration may lower the probability of being noticed by predators, it may be disadvantageous for escaping predators if they are already in pursuit. Third, although both early metamorphosis and large metamorphic mass are considered beneficial for amphibians in general (Smith, 1987;Berven, 1990;Altwegg and Reyer, 2003), sex-reversed individuals in our study tended to perform either well or poorly in both traits. Those sex-reversed individuals that did well in these traits tend to conform to the theory of heritable, random variation in the propensity for sex reversal, because most of them originated from two families with only XX genotypes (suggesting a sexreversed sire; see above). Almost all sex-reversed individuals that did poorly in both development and growth appeared unsuccessful also in executing sex reversal completely, as they had oogonia in their testes. This supports the above idea that sex-reversed individuals may represent a heterogeneous group whose life history and health might depend on the etiology of sex reversal.
When comparing males and females (either genetically or phenotypically), we found no difference in development and growth, and only a few differences in behavior. Concordant females were the only group that reacted to predator cues by less intense startle response and longer freezing. This may indicate lower risk taking by females, which may agree with the behavior of adult agile frogs observed in nature, where females were reported to forage less in open areas than males (Cicort-Lucaciu et al., 2011). Phenotypic males, including sexreversed individuals and concordant males, did not show the same responses to predator cues as females did, suggesting that FIGURE 5 | Duration of escape and freezing after the startle stimulus with and without predator cues in sex-reversed individuals (XX males), concordant (XX) females, and concordant (XY) males in experiment 2, excluding those individuals that did not respond to the stimulus by moving. In each box plot, the thick middle line and the box represent the median and interquartile range, respectively; whiskers extend to the most extreme data points within 1.5× interquartile range from the box. sex differences in these aspects of risk-taking behavior may not be genetically determined, but rather may develop after sex determination, e.g., by sex hormones. However, because none of the male-female differences in our study were statistically significant despite the relatively large sample for both sexes, we conclude that most of the divergent life histories and behaviors making up sex roles in agile frogs do not seem to arise in their larval life. In this species, males search and compete actively for females at high densities (Lodé et al., 2004), whereas at lower densities males maintain territories and females appear to choose males by their call characteristics (Lesbarrères et al., 2008), but both parents abandon the eggs after spawning. It would be very interesting to perform similar studies with species where either the male or the female parent takes the risky job of providing care to the offspring, as the developmental determinants of sex roles and therefore the effects of sex reversal may vary greatly between traditional and sex-role reversed systems. Amphibians and fish offer excellent model systems for such studies given their great diversity in mating and parental-care systems (Mank et al., 2005;Vági et al., 2019), but notably, sex reversal can also be experimentally induced in birds and mammals to study the development of sex roles (Adkins-Regan and Wade, 2001; Renfree et al., 2014).
Since sex reversal occurs relatively rarely under natural circumstances, most of our existing knowledge about ecologically relevant sex reversal comes from studies that include relatively small numbers of sex-reversed individuals in each population, year or treatment group (e.g., Li et al., 2016;Lambert et al., 2019;Jones et al., 2020;Nemesházi et al., 2020). The present study is no exception to this constraint. However, the fact that sexreversed individuals do not make up a large proportion of current populations does not mean that they are merely a curiosity: they may be powerful catalyzers of evolutionary change Nemesházi et al., 2021). Therefore, we call out for many more empirical studies to solidify our understanding of the evolutionary-ecological significance of sex reversal, and to extend it from a few species to a broad spectrum of ectothermic vertebrates faced with the challenges of rapid human-induced environmental change.

DATA AVAILABILITY STATEMENT
The datasets analyzed for this study, along with an annotated R script, can be found in the Figshare repository (https://doi.org/ 10.6084/m9.figshare.14974050).

ETHICS STATEMENT
The animal study was reviewed and approved by Ethical Commission of the Plant Protection Institute, Centre for Agricultural Research, Hungary.

AUTHOR CONTRIBUTIONS
VB conceived the idea, conducted the statistical analyses, wrote the first draft, and supervised the project. VB, EN, NU, ZM, and AH collected the data. RE analyzed the video recordings. EN and NV isolated DNA and performed PCRs. ZG and OH performed HRM analyses, which were then evaluated by EN. VB, OH, and AH acquired funding. EN supervised all molecular work. All authors contributed to writing the article and approved the submitted version.

FUNDING
This study was funded by the National Research, Development and Innovation Office of Hungary (NKFIH, grants 115402 and 135016 to VB, 124375 to AH, and 124708 to OH). The authors were supported by the János Bolyai Research Scholarship of the Hungarian Academy of Sciences (to VB, AH, and OH), the ÚNKP-20-5 New National Excellence Program of the Ministry for Innovation and Technology from the source of the National Research, Development and Innovation Fund ("Bolyai + Scholarship" to VB and AH), the Ministry of Human Capacities (National Program for Talent of Hungary, NTP-NFTÖ-17-B-0317 to EN), and by the Young Researcher program of the Hungarian Academy of Sciences (to NU). None of the funding sources had any influence on the study design, collection, analysis, and interpretation of data, writing of the manuscript, or decision to submit it for publication.