Novel genetic sex markers reveal unexpected lack of, and similar susceptibility to, sex reversal in free-living common toads in both natural and anthropogenic habitats.

Anthropogenic environmental changes are affecting biodiversity and microevolution worldwide. Ectothermic vertebrates are especially vulnerable because environmental changes can disrupt their sexual development and cause sex reversal, a mismatch between genetic and phenotypic sex. This can potentially lead to sex-ratio distortion and population decline. Despite these implications, there is scarce empirical knowledge on the incidence of sex reversal in nature. Populations in anthropogenic environments may be exposed to sex-reversing stimuli more frequently, which may lead to higher sex-reversal rate or, alternatively, these populations may adapt to resist sex reversal. We developed PCR-based genetic sex markers for the common toad (Bufo bufo) to assess the prevalence of sex reversal in wild populations living in natural, agricultural and urban habitats, and the susceptibility of the same populations to two ubiquitous oestrogenic pollutants in a common garden experiment. We found negligible sex-reversal frequency in free-living adults despite the presence of various endocrine-disrupting pollutants in their breeding ponds. Individuals from different habitat types showed similar susceptibility to sex reversal in the laboratory: all genetic males developed female phenotype when exposed to 1 µg L-1 17α-ethinylestradiol (EE2) during larval development, whereas no sex reversal occurred in response to 1 ng L-1 EE2 and a glyphosate-based herbicide with 3 µg L-1 or 3 mg L-1  glyphosate. The latter results do not support that populations in anthropogenic habitats would have either increased propensity for or higher tolerance to chemically induced sex reversal. Thus, the extremely low sex-reversal frequency in wild toads compared to other ectothermic vertebrates studied before might indicate idiosyncratic, potentially species-specific resistance to sex reversal.

INTRODUCTION

Anthropogenic environmental change confronts wildlife with challenges that call for rapid phenotypic or genetic adaptations. For example, both urban and agricultural land use loads the environment with various chemical pollutants including pesticides, heavy metals, road de‐icers, pharmaceuticals and industrial products, while climate change increases the frequency and intensity of heat waves and other extreme weather events. Understanding how these human‐induced environmental changes affect the ecosystem, and whether and how wildlife can adapt to overcome these challenges, is an important mission of current evolutionary ecology and conservation biology (Tilman et al., 2017).

In ectothermic animals, various environmental stimuli can cause a developmental effect rarely seen in endotherms: sex reversal, whereby individuals exposed to such stimuli during their embryonic or larval life phase develop the sexual phenotype opposite to their genetic sex (Baroiller & D’Cotta, 2016; Flament, 2016; Whiteley et al., 2021). Sex reversal occurs in fish, amphibians and reptiles in nature (Alho et al., 2010; Baroiller & D’Cotta, 2016; Lambert et al., 2019; Nemesházi et al., 2020; Whiteley et al., 2021; Xu et al., 2021), and theoretical studies caution that it may have far‐reaching consequences including skewed sex ratios, sex‐chromosome evolution and even population extinction (Bókony et al., 2017; Grossen et al., 2011; Nemesházi et al., 2021; Perrin, 2009; Schwanz et al., 2020; Wedekind, 2017). Laboratory experiments show that sex reversal can be induced by anthropogenic stressors such as chemical pollution and elevated temperature (Flament, 2016; Lambert et al., 2018; Mikó et al., 2021; Tamschick et al., 2016), and thus we may expect that the contemporary and future increase in the levels of anthropogenic stressors will influence the rates of sex reversal in free‐living populations of ectothermic vertebrates. Whether this influence would be an increased or decreased sex‐reversal frequency in anthropogenic environments is not a trivial question, for the following reasons.

On the one hand, sex reversal may happen more often in areas where the chemical and thermal stimuli triggering it are more pervasive, such as in agricultural areas polluted by pesticides and in urban heat islands. This may be simply a consequence of sex‐reversing stimuli being more frequent in such habitats. Alternatively or additionally, sex reversal might also be an adaptive response to anthropogenic environments, given that adjusting phenotypic sex to environmental conditions can be adaptive (Geffroy & Douhrad, 2019), and sex reversal might be an evolved mechanism for achieving the sexual phenotype that best matches the environment, similarly to environmental sex determination (Schwanz & Georges, 2021). These mechanisms may facilitate the spread of sex reversal in populations persisting in anthropogenic environments, especially because the propensity to develop into one sex or the other may exhibit genetic or epigenetic inheritance (McGaugh & Janzen, 2011; Piferrer & Anastasiadi, 2021).

On the other hand, sex reversal may be costly in terms of fitness. For example, reproductive performance of fish is reduced by sex reversal (Pandian & Sheela, 1995; Senior et al., 2012), as well as by intersex, a form of imperfect sex reversal when the gonads contain both male and female tissues (Fuzzen et al., 2015; Harris et al., 2011). Due to their sex‐chromosome genotype, sex‐reversed individuals may be unable to produce daughters or sons (Wedekind, 2017), and therefore may be selected against by sex‐ratio selection (Schwanz & Georges, 2021; see also figure S11a,d in Nemesházi, Kövér, et al., 2021). Also, sex‐reversed individuals may perform poorly in traits that influence survival (Mikó et al., 2021; Nemesházi et al., 2020) or sexually selected traits (Nemesházi, Kövér, et al., 2021). In such situations, we can expect resistance to sex reversal to be adaptive in environments where sex‐reversing stressors are pervasive. As a result of such adaptation, populations exposed to sex‐reversing environments might maintain the same or similar frequency of sex reversal as unexposed populations.

Assessing sex‐reversal frequencies in wild populations has been hindered by the difficulty of diagnosing sex reversal in nonmodel organisms. Due to the high evolutionary lability and homomorphy of sex chromosomes in ectothermic vertebrates, genetic sexing methods are available only for a small fraction of species (e.g., Alho et al., 2010; Baroiller & D’Cotta, 2016; Lambert et al., 2019; Nemesházi et al., 2020; Tamschick et al., 2016; Whiteley et al., 2021; Xu et al., 2021). In two such species, recently developed genetic sex markers have been used to investigate whether sex reversal is more prevalent in anthropogenic habitats, and they reported contradictory answers: yes in one frog species (Nemesházi et al., 2020) but no in another (Lambert et al., 2019). Furthermore, no study, to our knowledge, has yet tested whether animal populations living in anthropogenic habitats have increased or reduced inherent propensity for sex reversal.

In this study, we first aimed to produce a reliable molecular marker set for diagnosing genetic sex in the common toad (Bufo bufo), an anuran amphibian widespread in Eurasia that occupies a wide range of habitats from pristine woodlands to anthropogenic areas (Agasyan et al., 2009). This species has a female‐heterogametic (ZZ/ZW) sex‐chromosome system (Dufresnes et al., 2020), and is liable to chemically induced sex reversal (Hayes, 1998). Then, using our novel marker set, we investigated whether the frequency of sex reversal in toads differed between natural, agricultural and urban habitats. Finally, we performed a common garden experiment to test whether toads originating from these three types of habitat differ in their susceptibility to sex reversal induced by chemical pollutants.

We focused on the sex‐reversing effects of two endocrine‐disrupting chemical (EDC) compounds with high prevalence in surface water in agricultural and urban areas, respectively: glyphosate, the most used herbicide worldwide (Brovini et al., 2021), and 17α‐ethinylestradiol (EE2), a common ingredient of contraceptives that pollutes natural water bodies via wastewater (Bhandari et al., 2015). Both these EDCs may cause male‐to‐female sex reversal based on their effects on oestrogenic enzymatic activities, female‐skewed sex ratios and intersex gonads (Bhandari et al., 2015; Howe et al., 2004; Lanctôt et al., 2014; Tamschick et al., 2016). As both chemicals have been in use for about half a century, we can expect resistance to have potentially evolved in populations chronically exposed to these pollutants. Similar, rapid evolutionary changes due to anthropogenic habitat alterations have been documented in various taxa, including evolved tolerance to lethal effects of pollutants (Brans et al., 2021; Cothran et al., 2013; Johnson & Munshi‐South, 2017; Marques da Cunha et al., 2019; Reid et al., 2016). Here we test for altered susceptibility to sex reversal (a sublethal EDC effect) in common toad populations living in anthropogenic habitats.

METHODS

Data collection

We captured 352 adult toads during the spawning seasons of 2016 and 2017 at 14 breeding sites in north–central Hungary, which represented three habitat types: natural, agricultural and urban areas, with four or five sites per habitat type (Table S1, Figures S1 and S2). The habitats were categorized based on land cover within a 500‐m‐wide belt zone around each breeding pond (see Supporting Information: section 1). We identified the phenotypic sex of adults by sexual characteristics: nuptial pads in males (N = 216) and presence of eggs in females (N = 136). We took a DNA sample from each individual (buccal swab or tissue sample) and stored it in 96% ethanol.

In 2017, we transferred 89 pairs of the above‐mentioned captured adults to captivity and allowed them to spawn there, as described in Bókony et al. (2018). Depending on the availability of females and their willingness to spawn in captivity, we had 1–15 egg strings (families) from each of 11 sites out of the 14 sites sampled for adult DNA (36, 16 and 37 families from natural, agricultural and urban sites, respectively; Table S1). When the tadpoles hatching from the captive‐laid eggs reached the free‐swimming stage (developmental stage 25; Gosner, 1960), we haphazardly selected six individuals from each family, distributed them among six treatments (control, solvent control and four EDC treatments; N = 534: Table S2; see below), and raised them for ~5 months after metamorphosis as described in Ujhegyi and Bókony (2020). Because the methodological details of these procedures have already been published, we repeat them in the Supporting Information (section 2) of the present paper, and only the most important aspects of the experiment are described below.

The control group was kept in clean, filtered, reconstituted soft water, and served as control for the glyphosate treatments, in which a glyphosate‐based herbicide formulation (Glyphogan Classic; Monsanto Europe S.A.; containing 41.5% [w/w] glyphosate and 15.5% [w/w] polyethoxylated tallow amines) was added to the rearing water to maintain a nominal concentration of either 3 µg L−1 or 3 mg L−1 glyphosate. The solvent‐control group, in which the rearing water contained 1 µl L−1 ethanol, served as control for the EE2 treatments, in which the nominal concentration was either 1 ng L−1 or 1 µg L−1 EE2, obtained by dissolving EE2 powder (Sigma E4876) in 96% ethanol and adding 1 µl of this solution to each litre of rearing water. Actual EDC concentrations were close to the nominal concentrations (Ujhegyi & Bókony, 2020). Both EDCs are documented to occur in our actual study ponds (Bókony et al., 2018). The lower and higher concentrations we used for each EDC represent the typical and maximum concentrations, respectively, detected in surface waters (Avar et al., 2016; Bhandari et al., 2015; Bókony et al., 2018; Brovini et al., 2021). The treatments lasted throughout the entire larval period for each individual, and were renewed twice a week at each water change.

When the toadlets (N = 417) reached the age by which their gonads are completely differentiated (Ogielska & Kotusz, 2004), they were killed by using MS‐222, and we identified whether each individual had testes or ovaries by dissection. We stored the gonads in 10% buffered formalin and later examined them histologically (for detailed methods, see Nemesházi et al., 2020). In a few cases where we could not unambiguously categorize the gonads as testes or ovaries based on gross anatomy and histology, we treated the phenotypic sex as uncertain. We stored the body of dissected toadlets in 96% ethanol until extracting DNA from a foot sample (see section 3 in Supporting Information). Metadata for both adults and juveniles are publicly available on FigShare (Nemesházi, Sramkó, et al., 2021).

All captures and experimental procedures were 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 experiments were further approved by the Ethical Commission of the Plant Protection Institute (ATK NÖVI).

Marker development and validation

For developing and validating genetic sex markers, we used a reduced‐representation genomic library approach on toadlets with known phenotypic sex from the control group of the common garden experiment. Since these animals had not been exposed to any stimuli that are expected to cause sex reversal, they are likely to have phenotypic sex concordant with their genetic sex. First, we selected 24 nonsibling individuals sexed by gonad morphology (11 males, 13 females), representing 10 out of the 11 capture sites each by 1–4 toadlets. After DNA extraction (see section 3 in Supporting Information), we applied restriction‐site associated DNA sequencing (RADseq) to identify sex‐specific markers using the approach of Feron et al. (2021), which statistically examines RAD‐tags as being significantly associated with a priori sex.

To generate RADseq data, we adopted the original RADseq protocol of Baird et al. (2008) and used dual‐barcoded modified Illumina adapters. Next generation sequencing libraries were prepared using 210 ng genomic DNA (gDNA) as starting material from each isolate that was digested with the rare cutter restriction enzyme SbfI‐HF (New England Biolabs). Custom P1 adapters were ligated at sticky ends, and then isolates with a different adapter were pooled and sheared using a Bioruptor Pico machine (five cycles of 30 s “on” 30 s “off”). Libraries were size selected using the SPRI Select Kit (Beckman Coulter) to contain 300–600‐bp fragments only. Custom P2 adapters were ligated to the fragments and sublibraries were pooled equimolarly. PCR (polymerase chain reaction) enrichment used the Phusion High‐Fidelity PCR Master Mix (New England Biolabs) and a decreased number of PCR cycles (14) compared to the original protocol of Baird et al. (2008). After a final size selection the quality and quantity of the library were checked on a Bioanalyzer (Agilent Technologies) device, then the library was sequenced on an Illumina HiSeq platform with 150‐bp paired‐end sequencing option at a commercially available service provider (Novogene).

Raw Illumina reads were demultiplexed and filtered using process_radtags from the stacks version 2.2 pipeline (Rochette et al., 2019). Adapter content was additionally checked and removed by using fastp version 0.20.1 (Chen et al., 2018). Using the “forward” (R1) reads only, we screened the data set for sex‐specific reads of the 24 samples with known sex by radsex version 1.1.2 (Feron et al., 2021). We identified significantly sex‐linked markers by setting the significance threshold to the False Discovery Rate adjustment (Benjamini et al., 2001) of the p = .05 threshold adjusted for the number of tests (N = 168 combinations of number of males and number of females in which the given marker is present; sex linkage tested by Pearson's χ2 test of independence with Yates’ correction for continuity). We set the minimum read depth to one, thus allowing discovery of the maximum number of potentially sex‐linked markers. Read depth and distribution of sex‐specific markers were checked with the sgtr package (Feron et al., 2021) in r (R Core Team, 2014). Since radsex can only process reads that totally overlap (practically restricting this step to the “forward” reads), we retained sequence information of the “reverse” (R2) reads by picking the reads identified by radsex as significant markers and assembling “contigs” using the corresponding paired‐end reads by using the GNU/Linux utility “grep.” First, we extracted all exact sequence matches of the forward reads to the output sequences of radsex signif, then searched for the read pairs by their unique read identifiers. These short reads were clustered into contigs by cd‐hit version 4.8.1 (Fu et al., 2012) with a sequence similarity of 1.0 and pear version 0.9.6 (Zhang et al., 2014), as implemented in ddocent version 2.7.8 (Puritz et al., 2014), which is designed for the de novo processing of RADseq data sets.

Using the ncbi primer designing tool (https://www.ncbi.nlm.nih.gov/tools/primer‐blast), we designed sequencing primers for PCR amplification of potentially sex‐linked loci from the assembled RAD loci in order to obtain Sanger sequences of them. After PCR‐optimization (conditions for each sequencing primer pair are available in Table S3), loci that gave a bright PCR‐product band in the expected size ranges on a 2% agarose gel were further processed. Such PCR products of three female and one male laboratory‐raised juvenile toads from the control group were cut and purified from the gel using the NucleoSpin Gel and PCR Clean‐up Kit (Macherey‐Nagel) and were sequenced on a 3130xl Genetic Analyzer (Thermo Fisher Scientific) by a commercially available service provider (BIOMI Kft.). We initially sequenced more females than males to obtain multiple copies from both sex chromosomes (i.e., a total of five Z and three W copies from the three ZW females and one ZZ male; note: female‐associated alleles may not exclusively be found on the W chromosome, but for simplicity hereafter we will refer to these as W alleles). Sequences from each locus were manually checked using the staden software package (Bonfield et al., 1995) and subsequently aligned in mega (version 7.0.26). We focused on those loci where we both obtained unambiguous sequences and found sex‐linked sequence differences, which were either single nucleotide polymorphisms (SNPs) or insertions and deletions (InDels). These loci were sequenced in five males and five females in total for the purpose of designing diagnostic sexing primers. For four markers, we developed optimized sexing primers as described in Supporting Information section 4 (see Figure S4 for an illustration of W‐allele sequencing). When a Z/W InDel difference was large enough for detection on agarose gel, we designed sexing primers that would bind to both sex chromosomes (Z/W primers) and yield fragments of different length. For other loci, we added a digestion step where the W product was cut into two fragments by a restriction enzyme to be sex‐specific. When the above methods failed, we included a third, W‐specific, primer in the sexing PCRs to obtain a clearly distinguishable W product along with the products of the Z/W sexing primers.

We subsequently tested the developed sexing method for each sex marker (c2, c5, c12 and c16; see sexing primers in Table 1) in 46 males and 36 females with unambiguous sexual phenotype from the control group, all being nonsiblings and representing all 11 study sites that were used for the common garden experiment (including those used for sequencing and primer design). Additionally, we searched the common toad genome (NCBI GenBank identifier: aBufBuf1.1) in order to identify the sex chromosomes based on our sex markers (see Supporting Information section 5).

TABLE 1

Sexing primers

Locus	Primer ID	Primer sequence a	Restriction enzyme	Product size	WW detection b
c2	BbS2r‐F	TCCCCACACAAAGGAGAATGG	—	Z&W: 170 bp; W: 316 bp	No
	Bb_c2‐W‐F	TGTTCTATGCACTATGTGG
	Bb_c2‐R	ATCATCGAAGGGAAGAGCCG
c5 c	BbS5‐F2	CACAGCCCCTTCTTTGCTAAC	—	Z: 406 bp; W: 376 & 238 bp	Yes
	Bb_c5‐R	CTGGACGT A TGTTCTCCACG
	Bb_c5‐W‐R	GGGCCAATTTTTTGGAGAAG
c12 d	Bb_c12‐F	G T TCGGTCCCTCCTGAACG	TaiI	Z: 402 bp; W: 287 & 112 bp	Yes, if W is digested completely
	BbS12r‐R	GCCTAACCCGATGAAGCCG
c16	BbS16‐F3	GAACAGGGCGCCCACAC	—	Z: 218 bp;	Yes
	BbS16‐R	ACTCCAATCTCCAGAACGGC		W: 258 & ca. 290 bp

Nucleotides highlighted are SNP positions and are either W (underlined italic) or Z versions (underlined bold). These primers nevertheless bind to both sex chromosomes.

Although we did not find WW individuals in our data set, theoretically this genotype could be detected with c5, c12 and c16.

Among the two W‐specific products, the 376‐bp fragment is amplified by BbS5‐F2 and Bb_c5‐R, while the 238‐bp fragment is amplified by BbS5‐F2 and Bb_c5‐W‐R. WW identification is possible after ~180 min of electrophoresis at 100 V on a 2% agarose gel.

W‐products represent the post‐restriction fragment sizes.

Identification of sex reversals

We used our novel sexing primers to distinguish sex‐reversed and sex‐concordant individuals among both free‐living adults and their juvenile offspring raised in our common garden experiment. Preparation of DNA samples for subsequent analyses is described in section 3 in the Supporting Information. Sexing PCRs were performed on a Life ECO TC‐96/G/H(b)C (Bioer) or a Biometra Tone 96G (Analytic Jena) instrument using one of the following two touch‐down protocols. PCRs of c5, c12 and c16 were performed as follows: 2 min denaturation at 94°C followed by 13 cycles of 30 s denaturation at 94°C, 30 s annealing gradually decreasing from 68 to 64°C (−0.3°C per cycle) and 30 s elongation at 72°C, followed by 22 more cycles with the same settings but constant 64°C annealing temperature, and a final 10 min extension step at 72°C. Amplification of c2 differed from the above PCR profile: the annealing temperature decreased from 66 to 60°C in three cycles (touch‐down phase), followed by 32 cycles with an annealing temperature of 60°C. PCRs of all markers were performed in a total volume of 10 µl containing 2 µl FIREPol Master Mix (5×, ready to load; Solis BioDyne), 5 µl unquantified DNA for swab samples or 1 µl unquantified DNA for tissue samples, and varying amount of PCR primers and nuclease‐free water. In the two‐primer PCRs (c12 and c16), 0.4 pmol of each PCR primer (i.e., one forward and one reverse) was used. In the three‐primer PCRs (c2 and c5), a mixture of 0.05 pmol BbS2r‐F, 0.45 pmol Bb_c2‐W‐F and 0.4 pmol Bb_c2‐R performed best for c2, while 0.8 pmol BbS5‐F2, 0.175 pmol Bb_c5‐W‐R and 0.15 pmol Bb_c5‐R gave the best result for c5. To perform W‐specific digestion for c12, we subsequently added 1.43 µl Tango Buffer (10×; Thermo Scientific), 0.72 µl TaiI restriction enzyme (10 U µl−1; Thermo Scientific) and 2.85 µl nuclease‐free water to the PCR product, resulting in a 15 µl final volume. Digestion was performed at 65°C for 2 hr. With each marker, genetic sex was identified by electrophoresis on 2% agarose (peqGOLD Electran; VWR Peqlab) gel (0.5% TBE buffer; Thermo Scientific) stained with ECO Safe (Pacific Image Electronics).

Because both the Z‐linked and the W‐linked PCR products of c16 were similarly bright on the agarose gel, and genotyping required only a simple PCR with one sexing‐primer pair (i.e., no third primer or enzymatic restriction was necessary), we decided to use this marker for screening all those individuals from the common garden experiment that had not been sexed during the marker testing phase as well as the wild‐caught adults. If the c16 genotype of an individual did not match its phenotype, we genotyped the individual for c12 as well, to ensure correct assignment of the genetic sex. Individuals with uncertain sexual phenotype were sexed for at least two additional markers besides c16.

RESULTS

Based on 24 individuals, radsex identified 17 significantly sex‐linked markers (Figure 1), out of which 13 showed a female‐biased pattern (as expected under ZW/ZZ sex determination), but two of the latter were suspected to be paralogue sequences. During marker development, we concentrated on the remaining 11 RAD loci and designed sequencing primer pairs for each. Of these, nine primer pairs produced bright PCR products of the expected fragment size, and we obtained unambiguous sequences from seven loci (Table S3). We found sex‐linked InDel or SNP differences in the sequences of four loci hereafter referred to as c2, c5, c12 and c16 sex markers (sequences were registered at NCBI GenBank under the following accession numbers: OK507208–OK507215). NCBI genome blast search (see also Supporting Information section 5) indicated that c12 was located on chromosome 5, whereas c2 and c16 were localized on two different, unplaced scaffolds: the former is currently suggested to belong to chromosome 6, but no such information is available on the latter. For c5, genome blast showed highly similar sequences on several different chromosomes, including multiple locations on chromosome 5.

FIGURE 1

Tile plot showing the number of radsex markers found in different sexes. The shade of each tile refers to the number of markers that were found in a given number of males and females. Thick black frames around tiles show significantly sex‐linked markers, and the numbers within these tiles indicate the exact number of markers

Final sexing PCR primers for our four sex‐linked markers, product sizes and further details are shown in Table 1. Sexual genotypes based on each marker matched the sexual phenotype in all 82 nonsibling individuals chosen for marker validation, yielding 100% reliability for sexing with each of our four, newly devised markers.

We successfully genotyped 349 wild‐caught adults from 14 breeding ponds, while PCRs failed in three individuals (Table S1). We found 135 concordant ZW females and 213 concordant ZZ males, and a single sex‐reversed individual, a phenotypic male from an agricultural site (“Határrét”; see Figure S1) which was diagnosed as genetic female (ZW) by three out of four markers (c16, c2, and c12), while repeated PCRs with marker c5 gave ambiguous results.

We successfully genotyped all 417 toadlets that survived until phenotypic sexing in the common garden experiment. We detected no sex reversal in any of the treatment groups, except for the higher concentration of EE2 (Figure 2). However, detection of sex reversal might have been hindered by the high mortality rate in the treatment group of high concentration of the glyphosate‐based herbicide (see more details in Ujhegyi & Bókony, 2020). In the high‐concentration EE2 treatment, all genetic males developed into phenotypic females, regardless of their original habitat type (Figure 2). The ovaries of the male‐to‐female sex‐reversed individuals (ZZ females) were anatomically and histologically indistinguishable from the ovaries of concordant females (Figure S3). Additionally, six individuals with uncertain sexual phenotype were found to be genetically males (ZZ). Four of these toadlets, all originating from urban ponds, had intersex gonads (Figure 3); three of them had been treated with the lower concentration of the glyphosate‐based herbicide and one with the lower concentration of EE2 (Figure 2). In the remaining two toadlets we could not unambiguously ascertain if the gonads were intersex or normal testes (Figure 3).

FIGURE 2

Proportion of concordant and sex‐reversed individuals in each treatment group by habitat type of the parents’ capture site. Bar widths are proportional to sample size, which varied between two and 35 due to differences in survival (see Table S2). GLY, glyphosate; EE2, 17α‐ethinylestradiol

FIGURE 3

Ambiguous gonads in juvenile common toads. (a) The gonad on the left is a normal testis and the gonad on the right is an ovary. This individual originated from an urban pond (Pesthidegkút) and was treated with 1 ng L−1 EE2. (b) The testes are abnormally shaped and Bidder's organ on the right has an ovary‐like structure. This individual originated from an urban pond (Pilisvörösvár) and was treated with 3 µg L−1 glyphosate. (c) A small ovary‐like structure (white arrow) between the testis and Bidder's organ. Two individuals treated with 3 µg L−1 glyphosate had this morphology; both originated from urban ponds (Göd, Pesthidegkút). (d) Histological image of the gonad from one of the individuals showing the gross anatomy in (c). A single oogonium (black arrow) is found in the testicular tissue. (e) Abnormal gross anatomy of testes and Bidder's organs in an individual from the solvent control group, originating from an urban pond (Pilisszentiván). The histological section of this individual was lost, so its phenotypic sex was categorized as uncertain. (f) Histological image of the gonad from the other individual whose phenotypic sex was categorized as uncertain. The cells shown by arrows may be testicular oogonia, or may belong to Bidder's organ. Gross anatomy showed normal testes. This individual originated from a natural pond (János‐tó) and was raised in the control group

DISCUSSION

The novel sex markers developed in our study confirmed that Hungarian populations of the common toad are female heterogametic (Figures S5 and S6), echoing recent findings from Switzerland (Dufresnes et al., 2020). Identification of the common toad's sex chromosomes has remained unresolved so far (Dufresnes et al., 2020). Genome blast showed that at least three of the four new sex markers are located on different scaffolds, but the specific position of these scaffolds is unknown, except for one which is located on chromosome 5. One of the unplaced scaffolds, however, was suggested to belong to chromosome 6. Some anuran species feature complex sex determination that can include multiple sex chromosomes (Gazoni et al., 2018; Roco et al., 2015). In such cases, sex chromosomes can form a meiotic chain, and therefore multiple loci of different chromosomes may be inherited in linkage. Because each of our sex markers showed 100% sex linkage, they are probably all located on the sex chromosome(s). Based on the available data, we propose chromosome pair 5 (most similar to chromosome 6 in Xenopus tropicalis) as a candidate sex chromosome pair in the common toad, but the exact number of sex chromosomes operating in this species has yet to be determined. Irrespective of the exact location on the chromosomes, our new marker system enables accurate identification of the genetic sex in the common toad. While the majority of amphibian species to which genetic sexing methods have been established feature XX/XY sex determination (Alho et al., 2010; Lambert et al., 2019; Nemesházi et al., 2020; Xu et al., 2021), our new marker set provides a cheap and easy‐to‐use method for future studies aiming to understand sex‐reversal mechanisms in an anuran species with ZW/ZZ sex determination.

Our field study on several hundreds of adult toads found only a single case of sex reversal: three out of our four sex markers confirmed that the W chromosome was present in the DNA sample of one adult male captured at an agricultural site. The fourth marker (c5) gave inconclusive result, cautioning against use of this marker for genetic sexing in future studies. Given the lack of further swab samples from this animal, we cannot completely exclude the possibility of contamination. Nevertheless, we had the highest number of phenotypic male samples from the pond where this individual was captured (Figure S1), and thus finding a single sex reversal at this site is compatible with the idea of an existing but very low frequency of sex reversal in the studied common toad populations. This almost complete lack of sex reversal is surprising, because we found many EDCs in the studied ponds, with higher concentrations in anthropogenic areas, as published in earlier papers (Bókony et al., 2018, 2021). Furthermore, we found a considerable number of female‐to‐male sex‐reversed agile frogs (Rana dalmatina) in some of these ponds in the same years (Figure S1; Nemesházi et al., 2020). Thus, the lack of sex reversal in toads cannot be explained by the general lack of sex‐reversing effects in the studied sites. Instead, this result may suggest that toad populations living in more polluted areas might have evolved resistance to sex reversal, thereby showing the same undisrupted sex development as their conspecifics in natural habitats. However, this interpretation is not supported by the results of our common garden experiment, because the effects of sex‐reversing EDC treatments on the offspring of the studied toads did not depend on their original habitat type. Instead, they either all showed no sex reversal at low concentration of both EDCs, or showed a 100% male‐to‐female sex reversal in the presence of high EE2 concentration (Figure 2). We found only a slight indication of habitat dependence of EDC susceptibility, suggesting that toads originating from anthropogenic habitats may be more, not less, susceptible to disrupted sex development: only urban toadlets displayed intersex gonads in a few cases when treated with ecologically realistic, low EDC concentrations. It remains to be tested if other EDC compounds or other concentrations within the range of the realistically low and close‐to‐maximum values that we applied here would reveal habitat‐dependent sex‐reversal probabilities in toads or any other species liable to sex reversal. Nevertheless, because most EDCs found in amphibian breeding habitats have oestrogenic potential (Bókony et al., 2018), our treatments provide a good overall representation of the oestrogenic EDC effects probably present in the field.

As a possible explanation for our results, the survival rate of sex‐reversed juveniles might be low in the wild, resulting in a low prevalence of sex‐reversed individuals among adults. The environmental stimuli that cause sex reversal may have other developmental effects that might reduce survival; for example, heat stress in agile frogs increases both sex‐reversal rate and mortality (Mikó et al., 2021). However, a meta‐analysis found no significant relationship between chemically induced sex reversal and mortality in aquaculture fish (Senior et al., 2012). Similarly, in our present study, toadlet survival was not reduced in the treatment group that showed 100% male‐to‐female sex reversal (Table S2), although this might not be representative of their survival chances in the wild. Mortality can be especially high during the first winter hibernation (Üveges et al., 2016), which was not assessed in the present study. Moreover, sex‐reversed individuals may show different behaviour (Li et al., 2016; Senior et al., 2015). If their altered behaviour also affects their microhabitat use or results in changed activity during the breading season, these individuals might be harder to find by conventional capturing methods. There is currently very little information on the survival and behaviour of sex‐reversed individuals in nature (Wild et al., 2022), so testing the above ideas will require further research.

As an alternative explanation that is not mutually exclusive with the above hypotheses, we speculate that the common toad may be relatively resistant to sex reversal, regardless of habitat type. In all other anuran species studied so far for sex reversal in free‐living populations, female‐to‐male sex reversal was found in noticeable numbers (Alho et al., 2010; Lambert et al., 2019; Nemesházi et al., 2020; Xu et al., 2021), and a few cases of male‐to‐female sex reversal were also indicated (Lambert et al., 2019). There are several differences between the common toad and the previously studied anuran species, which might contribute to the apparent difference in sex‐reversal frequencies found in their wild populations. First, all species studied so far belong to the family Ranidae, whereas the common toad is a member of Bufonidae; and different phylogenetic lineages may show different sensitivity for certain sex‐reversing conditions (Chardard et al., 2004; Hayes, 1998; Orton & Tyler, 2015; Tamschick et al., 2016). Second, toads produce defensive toxins from cholesterol, the precursor of steroid hormones (Daly, 1995), and they have been selected for resistance to autotoxicity (Moore et al., 2009). This might have conferred them tolerance to other chemical perturbations which mimic the effects of steroid hormones (including oestradiols and “stress hormones”), similarly to the cross‐resistance provided by tolerance to certain pesticides in other anurans (Hua et al., 2014). Third, sex reversal may be triggered by not only chemical but also thermal stimuli, and different species may have adapted to different temperatures. If sex reversal in free‐living amphibians occurs mostly due to extreme temperatures (Lambert et al., 2018; Mikó et al., 2021), the lack of sex reversal in common toads might be explained by their higher tolerance to heat. In line with this idea, the breeding season starts ~1 month later in spring for common toads than for agile frogs in our study region; accordingly, we found female‐to‐male sex reversal in agile frogs (Nemesházi et al., 2020) but not in common toads among free‐living adults, and we found the same difference between the two species in an experimental study of heat‐induced sex reversal (Ujszegi et al., 2021). Similarly, evolution of different temperature thresholds for sex reversal was suggested to explain the finding that in a reptile, Pogona vitticeps, sex reversal is absent in the hottest part of the species’ range (Castelli et al., 2021).

What makes different populations and species more or less susceptible to sex reversal is an important question for evolutionary ecology as well as for conservation biology. Possible reasons include constraints such as the degree of sex‐chromosome heteromorphy (Miura et al., 2016) and direct or indirect selection pressures. For example, artificial selection for increased fecundity in females can indirectly affect male sensitivity to oestrogenic disruption of testis development and spermatogenesis (Spearow et al., 1999). Such selection pressures may force populations to evolve or plastically modulate any element involved in the biochemical pathway that translates environmental stimuli into sex (Castelli et al., 2020), such as by mutation of genes encoding hormone receptors (Castañeda Cortés et al., 2019; Hamilton et al., 2020). Thus, the vulnerability of phenotypic sex development may be shaped by multiple forces, which might explain why researchers have had mixed success in finding clear‐cut relationships of sex‐reversal rate with environmental factors such as climate (e.g., Castelli et al., 2021 vs. Dissanayake et al., 2021) and urbanization (Lambert et al., 2019 vs. Nemesházi et al., 2020) or with taxonomy (Senior & Nakagawa, 2013). Even when a clear correlation is present, the underlying mechanisms are difficult to ascertain: for example, oestrogenic pollution in river stretches is associated with a high frequency of intersex in fish but not with polymorphisms in genes involved in responses to EDCs (Hamilton et al., 2020). Our present results with common toads add to this complex picture, emphasizing the need for further research on sex reversal in a wide diversity of species. Building on our accumulated understanding from laboratory experiments on how environmental perturbations affect sex and from theoretical models on how sex reversal may impact population dynamics and evolution, the time is ripe for empirical studies on the causes and consequences of sex reversal in wild populations in the Anthropocene.

AUTHOR CONTRIBUTIONS

V.B. designed and supervised the study. Fieldwork was done by B.Ü. and V.B. Animals in the laboratory were raised by E.N., N.U. and V.B.; dissection and phenotypic sexing was performed by N.U. RADseq was performed by G.S., L.L. and L.S. All further DNA work was conducted by E.N., E.B. and N.V., supervised by E.N. E.N. and V.B. wrote the manuscript. All authors proofread the manuscript and gave final approval for publication.

BENEFIT‐SHARING STATEMENT

Benefits Generated: Benefits from this research accrue from the sharing of our data and results on public databases as described above.

CONFLICT OF INTEREST STATEMENT

The authors declare no conflict of interest.

Supplementary Material

Click here for additional data file.