Male heterogametic sex determination in Rana dybowskii based on sex-linked molecular markers.

Identifying the mechanism for sex determination in amphibians is challenging. Very little is known about sex determination mechanisms of Rana dybowskii, a species of importance to evolutionary and conservation biology. We screened for sex-linked molecular markers in R. dybowskii in China using target region amplification polymorphism with 2 fixed primers against the sequences of Dmrt1. We found 2 male-linked molecular markers in R. dybowskii, which were 222 bp and 261 bp long. The detection rates of 222 bp marker in males form Xinglong, Huadian, and Dandong were 93.79%, 69.64%, and 13.64%, respectively, while the rate in females from Huadian was 27.50%. Besides, the detection rates of 261 bp marker in the above 3 regions were only observed in males at the rate of 93.79%, 87.50%, and 32.73%, respectively. The inheritance patterns of sex-linked molecular markers showed that the 2 sex-linked molecular markers were heterozygous. Compared to the XY-male parent, progeny from XX-pseudo-male parent possessed lower sex reversal ratio at the same rearing temperature, and the proportion of female froglets from an XX-pseudo-male parent was more than 95% at low rearing temperature (15°C). Our findings suggest that R. dybowskii displays male heterogamety, and the 2 sex-linked molecular markers may have a guiding significance for the protection and utilization of R. dybowskii.

INTRODUCTION

The study of sex‐determination systems in amphibians has been challenging. In sharp contrast with mammals and birds, a few amphibians possess morphologically distinguishable heteromorphic sex chromosomes (Eggert 2004; Chang et al. 2017). Even when these heteromorphic sex chromosomes exist, it is a time‐consuming and resource‐demanding task to identify the heterogametic sex chromosomes using the classic direct method (Miura 1994; Matsuba & Merilä 2006; Pradit et al. 2018). Besides, variations in sex‐determination systems have also been observed in closely related species, populations of the same species, and even in individuals within the same populations due to the complexity and diversity of sex‐determination systems in amphibians (Sumida & Nishioka 2000; Miura 2017; Ito 2018).

Male and female heterogametic systems can generally be distinguished by sex‐linked markers (Berset‐Brändli et al. 2006; Gamble & Zarkower 2014; Gamble et al. 2015; Nemesházi et al. 2020). The amplified fragment length polymorphisms in Palearctic green toads have been used to identify genotypic sex (Stöck et al. 2011) and sex‐linked microsatellite markers could successfully identify an XX/XY sex system in the spiny frog (Yuan et al. 2017). The continuous development of next‐generation sequencing creates new breakthroughs for screening sex‐linked markers. Diversity arrays technology has also been applied to screen sex‐linked markers to identify the heterogamety patterns in the North American green frog (Lambert et al. 2016). A genotyping‐by‐sequencing approach has determined homomorphic sex chromosomes in the moor frog (Brelsford et al. 2017). Jeffries et al. (2018) used restriction‐site associated DNA sequencing (RADseq) to search for sex‐linked markers in 20 frog species, of which 12 species harbored these markers. However, not all attempts using next‐generation sequencing have successfully explored sex‐linked markers to identify sex‐determination systems. Brelsford et al. (2016) analyzed Rana temporaria from a Swiss lowland population with RADseq markers. They mapped a total of 2177 single nucleotide polymorphisms, of which none of them showed significant associations with offspring sex. Similarly, Jeffries et al. (2018) found no sex‐specific markers in 8 species of frog, including Rana chensinensis.

Rana dybowskii is a frog that lives in valley groves of the northeast of China, the Russian Far East, the southeast of East Siberia, the Korean peninsula, and Japan (Amphibia Web 2007). It was earlier thought to be a subspecies of R. chensinensis; however, it has now been classified as an independent species according to the difference in distribution range and morphology (Xie et al. 1999). The heteromorphic sex chromosomes of R. dybowskii from different places in the northeast of China have not been found using cytometry methods such as chromosomal karyotype analysis and Ag‐band (Shao et al. 1999; Zeng et al. 2001). Considering its importance in ecological and economic value, we have previously examined sex‐linked molecular markers of R. dybowskii by random amplification of polymorphic DNA (RAPD), sequence‐related amplified polymorphism (SRAP), and RADseq; however, we did not find evidence of supporting male or female heterogamety (unpublished data).

In this study, we used the target region amplification polymorphism (TRAP) to screen for sex‐linked molecular markers of R. dybowskii and analyzed their geographic variation. We found that 2 male‐linked markers existed in male R. dybowskii. Furthermore, we analyzed the inheritance patterns of the 2 male‐linked markers and sex reversal of the offspring from 3 males and 2 pseudo‐males for the father.

MATERIALS AND METHODS

Marker development

Sampling and DNA extraction

Between 2016 and 2018, R. dybowskii were collected from 3 populations in northeast China (Fig. 1; Table 1) including 198 individuals (145 males and 53 females) from Xinglong, Heilongjiang Province (128°28′E, 46°38′N); 96 individuals (56 males and 40 females) from Huadian, Jilin Province (126°50′E, 43°7.8′N); and 166 individuals (110 males and 56 females) from Dandong, Liaoning (124°58′E, 40°34′N). DNA was extracted from the toes of individuals by phenol‐chloroform method described by Xu et al. (2019) for TRAP‐PCR.

Figure 1

Geographical localization of the 3 populations analyzed in northeast China.

Table 1

Sex‐linked molecular marker frequencies of R. dybowskii from Xinglong (n = 198), Huadian (n = 96) and Dandong (n = 166)

		Proportion of individuals amplified the sex‐linked molecular marker
Population	Phenotypic sex	Dmrt1‐1/EM4	Dmrt1‐2/EM7
Xinglong	Male	93.79%	93.79%
	Female	0.00%	0.00%
Huadian	Male	69.64%	87.50%
	Female	27.50%	0.00%
Dandong	Male	13.64%	32.73%
	Female	0.00%	0.00%

Primer design and polymerase chain reaction

In this study, 20 arbitrary primers (Table 2) for TRAP are taken based on reports by Gao et al. (2008), Hu and Vick (2003), and Guo et al. (2012). Arbitrary primer sequences for TRAP included 3 selective nucleotides at the 3′ end, 4 nucleotides of AT‐ or GC‐rich content in the core region, and 11 or 10 nucleotides as filler sequences at the 5′ end. Two fixed primers (Table 2) for TRAP markers were designed based on 2 R. dybowskii sequences of Dmrt1 (NCBI Accession No.: MW316654 and MW316653). Forty pairs of primers consisting of arbitrary primer and fixed primer were used to amplify sex‐linked molecular markers. PCRs were performed in a total volume of 10 µL including 1 µL 10× Taq Buffer (Tris, KCl, MgCl2, and H2O, Vazyme, China), 0.2 µL of 10 mM dNTPs (dATP, dCTP, dGTP, dTTP), 0.5 U DNA polymerase (Vazyme, China), 0.4 µL forward primer (10 pm/µL), 0.4 µL reverse primer (10 pm/µL), 0.8 µL of extracted DNA (50 ng/µL), and supplemented with ddH2O to 10 µL. PCRs were conducted as follows: 5 min of Taq polymerase activation at 95°C, followed by 5 cycles of denaturation at 94°C for 1 min, annealing at 35°C for 1 min, and elongation at 72°C for 1 min, and 35 cycles of denaturation at 94°C for 1 min, annealing at 51°C for 1 min and elongation at 72°C for 1 min, and ending with a final elongation of 7 min at 72°C. PCR products were analyzed by electrophoresis in 14% polyacrylamide gel and silver stain. Sex‐linked bands were sequenced by Comate Biosciences Co., Ltd. (Comate Biosciences, Changchun, China).

Table 2

Primer sequences and pairs of primers used for TRAP analysis

Primer		NCBI Accession No.	Sequence (5′→3′)
Fixed (forward)	Dmrt1‐1	MW316654	GGCTATTCGTCGCTACTAAAGG
	Dmrt1‐2	MW316653	GGTCATTCCTTGTCCTAATTATCAGT
Arbitrary (reverse)	Me1		TGAGTCCAAACCGGATA
	Me2		TGAGTCCAAACCGGAGC
	Me3		TGAGTCCAAACCGGAAT
	Me4		TGAGTCCAAACCGGACC
	Me5		TGAGTCCAAACCGGAAG
	Me6		TGAGTCCAAACCGGTAA
	Me7		TGAGTCCAAACCGGTCC
	Me8		TGAGTCCAAACCGGTAC
	Me9		TGAGTCCAAACCGGTCA
	Me10		TGAGTCCAAACCGGTCG
	EM1		GACTGCGTACGAATTAAT
	EM2		GACTGCGTACGAATTTGC
	EM3		GACTGCGTACGAATTGAC
	EM4		GACTGCGTACGAATTTGA
	EM5		GACTGCGTACGAATTAAC
	EM6		GACTGCGTACGAATTGCA
	EM7		GACTGCGTACGAATTGAG
	EM8		GACTGCGTACGAATTGCC
	EM9		GACTGCGTACGAATTTCA
	EM10		GACTGCGTACGAATTCAA
Pairs of primers	Combination 1	Dmrt1‐1/EM4
	Combination 2	Dmrt1‐2/EM7

The inheritance patterns of sex‐linked molecular markers and sex reversal

Field sampling and husbandry

Twenty mating pairs were captured in amplexus during the 2019 breeding season from Xinglong, Heilongjiang Province, and all adults were genotyped using these 2 screened sex‐linked molecular markers. Five clutches from 5 fathers (MF1, MF2, MF3, PMF1, and PMF2) were collected, and each clutch was divided into 3 groups, with 240 fertilized eggs in each group. Three parts from each clutch were reared at 15°C, 20°C, and 25°C, respectively, at indoor facilities of the Northeast Agricultural University (Table 3).

Table 3

The genotypic sex ratio of offspring and the ratio sex reversal of offspring raised under different temperature conditions

					Phenotypic sex
Father	Rearing temperature	Code of clutches	Genotypic sex (female/male)		Female/male	Female%	Ratio of sex reversal (%)
Male	15°C	MF1	110/128	912/896 ≈ 1/1 (χ2 = 0.0542, P = 0.8159)	104/134	43.70	5.45
	15°C	MF2	141/80		133/88	60.18	5.67
	15°C	MF3	117/100		116/101	53.46	0.85
	20°C	MF1	109/124		90/144	38.46	18.35
	20°C	MF2	105/115		96/124	43.64	8.57
	20°C	MF3	113/115		76/152	33.33	32.74
	25°C	MF1	72/87		40/119	25.16	44.44
	25°C	MF2	72/85		32/125	20.38	55.56
	25°C	MF3	73/62		23/112	17.04	68.49
Pseudo‐male	15°C	PMF1	168/0	917/0	162/6	96.43	3.57
	15°C	PMF2	172/0		165/7	95.93	4.07
	20°C	PMF1	181/0		164/17	90.61	9.39
	20°C	PMF2	166/0		152/14	91.57	8.43
	25°C	PMF1	101/0		87/14	86.14	13.86
	25°C	PMF2	129/0		115/14	89.15	10.85

Progeny sexing

The offspring were euthanized in 0.2% ethyl3‐aminobenzoate methanesulfonate salt solution (MS222) when tadpoles reached Gosner stages 45 (Gosner 1960). For DNA extraction and amplification of sex‐linked molecular markers as described above, head and hind limbs were soaked in 75% ethanol and preserved at −20°C to identify genotypic sex. Chi‐square test was computed to compare the variation of genotypic sex. Besides, the abdomen was fixed in Bouin's fluid for 24 h, then soaked in 75% ethanol, and preserved at 4°C for gonad scan under a dissecting microscope. Paraffin sections were made from 100 observed gonads as described by Huang and Bai (2003). Then, the testes of the froglets at Gosner stage 45 were classified into 5 and ovaries were classified into 4 categories, based on their morphology (Figs S1 and S2, Supporting Information). Phenotypic sex of the remaining gonads without paraffin section was identified by this classification. Male offspring without the 2 sex‐linked molecular markers were considered to be sex reversals, and their ratios to all individuals without these markers were calculated.

RESULTS

Two male‐linked molecular markers in R. dybowskii

Total 40 TRAP primer pairs were tested in the pre‐screening experiment. Only 2 of them could produce sex‐linked bands in R. dybowskii from Xinglong and were selected in the following study (Table 2). TRAP primer combination 1 (Dmrt1‐1/EM4) amplified a consistent band of 222 bp in 14 out of 16 males, completely absent in all 16 females (Fig. 2a, named MSM‐222). Similarly, the other sex‐linked marker of 261 bp was amplified with TRAP primer combination 2 (Dmrt1‐2/EM7) and found in 14 out of 16 males, completely absent in all 16 females too (Fig. 2b, named MSM‐261). The remaining 2 males lacking the 2 sex‐linked markers were considered as pseudo‐male, which is the induction of genotypic female‐to‐phenotypic male sex reversal. The complete sequences of male sex‐linked markers produced by TRAP primer combination 1 (Dmrt1‐1/EM4) and combination 2 (Dmrt1‐2/EM7) were 222 bp and 261 bp in length, respectively, as shown by DNA sequencing (Supplement 1). BLAST searches of GenBank showed that the sequence amplified by Dmrt1‐1/EM4 primers was similar to satellite DNA from Rana arvalis (NCBI Accession No.: HF564834), Rana pyrenaica (NCBI Accession No.: AM900346), and Rana iberica (NCBI Accession No.: AJ555239), and the corresponding hit was 49%, 46%, and 43%, respectively. Besides, no hit was observed for the sequence amplified by Dmrt1‐2/EM7 (261 bp) primers.

Figure 2

Two male‐linked molecular TRAP markers isolated from 14 female individuals of R. dybowskii in Xinglong population. (a) Male‐linked molecular marker (222 bp) obtained with the TRAP primer combination 1 (Dmrt1‐1/EM4). (b) Male‐linked molecular marker (261 bp) obtained with the TRAP primer combination 2 (Dmrt1‐2/EM7). The individual with ‘*’ were pseudo‐male.

Geographic variation of sex‐linked molecular TRAP marker in R. dybowskii

R. dybowskii from 3 populations were collected for validation of the male‐linked molecular markers described above and analyzed for geographic variations (Table 1). The detection rates of 222 bp sex‐linked molecular markers from TRAP primer combination 1 in males from Xinglong, Huadian, and Dandong were 93.79% (136/145), 69.64% (39/56), and 13.64% (15/110), respectively, while the rate in females from Huadian was 27.50% (11/40). Unlike the markers from TRAP primer combination 1, we could not detect the sex‐linked molecular marker from TRAP primer combination 2 in any females, regardless of their region of origin. In contrast, the frequency of the sex‐linked molecular marker from TRAP primer combination 2 in males from Xinglong, Huadian, and Dandong was 93.79% (136/145), 87.50% (49/56), and 32.73% (36/110), respectively.

The inheritance patterns of sex‐linked molecular markers

Clutches of 3 fathers with the 2 sex‐linked molecular markers from Xinglong were used to analyze the inheritance patterns of the 2 sex‐linked molecular markers. The DNA from 1808 progeny of 3 fathers was analyzed for the 2 sex‐linked molecular markers (Table 3); among these, 912 harbored the markers while the remaining 896 did not. A chi‐square test revealed no significant sex variation (χ2 = 0.0542, df = 1, P = 0.8159). Like fathers, 2 sex‐linked molecular markers were simultaneously amplified by partial sons (individuals with phenotypic male). The remaining sons and all daughters (individuals with phenotypic female) did not amplify these sex‐linked molecular markers. This indicated that the male parent was heterozygous.

Sex reversal of offspring from male and pseudo‐male parent

The genotypic female‐to‐phenotypic male sex reversal was observed in R. dybwoskii offspring reared at different temperatures from both males (with the 2 male‐linked molecular markers) and pseudo‐males (missing male‐linked molecular markers). The sex reversals ratio of the offspring from 3 male parents, shown in Table 3, are 0.85–5.67% at 15°C, 8.57–32.74% at 20°C, and 44.44–68.49% at 25°C. The result shows that there is a positive trend with the rearing temperature. Similar patterns occurred in the 2 families from pseudo‐male parents. Compared to male parents, progeny from pseudo‐male parents possessed a lower ratio of sex reversal at the same rearing temperature, which was more prominent in group of 20°C and 25°C (8.43–9.39% vs 8.57–32.74% at 20°C and 10.85–13.86% vs 44.44–68.49% at 25°C). Interestingly, the proportion of female froglets from a pseudo‐male parent was more than 95% at 15°C.

DISCUSSION

Polymerase chain reaction (PCR)‐based TRAP was used to detect genetic variation at the DNA level (Hu & Vick 2003). Here, we used TRAP to screen for a sex‐linked marker of R. dybowskii. Dmrt1 is either sex‐linked or possibly important in the sexual differentiation of several Hylidae, Rana temporaria, Rana rugosa, and Bufo viridis (Uno et al. 2008; Brelsford et al. 2013, 2016; Lambert et al. 2016). Considering that Dmrt1 could be a candidate of the sex‐determining gene in frogs (Nakamura 2012; Brelsford et al. 2017; Miura 2017), we employed 2 fixed primers from the Dmrt1 sequences of R. dybowskii and 2 selected arbitrary primers, which consist of GC‐rich “core” regions for annealing with an intron corresponding to each fixed primer. This study was one of the few researches to use TRAP‐PCR to develop sex‐linked molecular markers of animals without a genetic linkage map or complete genome sequencing. In this study, the sequence of the sex‐linked molecular markers was not hit to the exons of Dmrt1 except the region of fixed primers; however, it was possible that the 2 sex‐linked molecular markers were across intron and exon of Dmrt1, or in the intron of Dmrt1, considering that arbitrary primers had regions for annealing with an intron. Due to the relatively short molecular marker length of the 2 markers, up‐and downstream sequences of them should be amplified to analysis their relationship with Dmrt1 in R. dybowskii for further research.

Most amphibians have homomorphic sex chromosomes, which make it challenging to identify sex‐determination systems, despite all amphibians having a genetic component to sex determination. Sex‐linked markers had been widely used in inferring the sex‐determination systems of frogs (Lambert et al. 2016; Brelsford et al. 2017; Jeffries et al. 2018). Two male‐linked molecular markers of R. dybowskii were screened in this research. In this context, inheritance patterns of male‐linked molecular markers between parents and offspring suggested that R. dybowskii probably possessed the XX/XY sex‐determination system. These discoveries provided new clues for further studies of a sex‐determining gene in the genus Rana.

Geographic variations in sex determination are common in the same amphibian (Miura et al. 1998; Rodrigues et al. 2014) and reptile (Ewert et al. 2005), which appears in 2 aspects based on sex‐linked molecular markers. The first one is that sex is related to different sex‐linked molecular markers because of population diversity. In this context, Cano et al. (2011) found 7 markers of R. temporaria linked to phenotypic sex. Among them, 6 were mapped to a linkage group in wild individuals in Northern Finland and Sweden, the seventh marker mapped to another linkage group in wild individuals in Southern Sweden. The present results showed that the male‐linked molecular markers of 222 bp appeared in male individuals from 3 populations and a small portion of females from Huadian. The sex‐linked molecular markers of the 261 bp band appeared to be stricter, and it occurred only in male individuals from 3 populations. In all these populations, males in the Xinglong population shared the same proportion of the 2 male‐linked molecular markers. The reasons for these variations mainly include sex chromosome recombination during evolution (Cano et al. 2011; Lambert et al. 2016). Also, the associations of a sex‐linked marker with phenotypic sex were variant among different populations. Although the association between phenotypic sex and sex‐linked microsatellite markers of linkage group 2 has been found to be perfect in a northern Swedish population of R. temporaria, sex‐linked microsatellite markers show decreasing levels of sex linkage with decreasing latitude in Sweden (Rodrigues et al. 2013, 2014, 2016). Similarly, different proportions of males from 3 populations possessed the 2 male‐linked markers in our study, and males from higher latitude population (Xinglong) owned the highest proportion. The possible main reason for this was due to the variation of the environmental factors such as temperature (Wallace & Wallace 2000; Alho et al. 2010), anthropogenic land cover (Nemesházi et al. 2020), or chemical exposure (Pettersson & Berg 2007; Hayes et al. 2010), besides sex chromosome recombination. In this study, interpretating the nature of geographical variation of the 2 male‐linked markers in R. dybowskii among populations along an 800‐km transect across northeast China was advantageous for revealing clues into the origin, evolution, and migration of R. dybowskii.

Amphibians are known to have genotypic sex determination (GSD) (Wallace et al. 1999); however, gender differentiation of frogs is determined not only by genetic factors but also affected by environmental factors such as rearing temperature, exterior hormones, and radiation (Ito 2018), leading to individuals for whom phenotypic sex does not match genotypic sex (such as XX‐males, XY‐females, ZZ‐females, and ZW‐males). This mismatch between phenotypic and genotypic sex is caused by sex reversal (Miura 1994; Wallace et al. 1999; Alho et al. 2010; Ito 2018). Often, sex reversal is inferred in laboratory experiments by biasing the sex ratio (Li et al. 2001; Phuge 2017; Tang et al. 2020). However, due to a lack of heterotypic sex chromosomes for lower vertebrates, it is difficult to identify their genotypic sex, resulting in the barrier to obtaining information on sex reversals in natural populations. With the development of molecular biology, sex reversal can be diagnosed by the discordance between the phenotypic sex and the majority of sex‐linked markers (Wallace et al. 1999; Alho et al. 2010; Lambert et al. 2016; Flament 2016). Three Chinese giant salamanders reversed from genetic female to physiological male were found in a group exposed to elevated temperature, and 13 individuals reversed from genetic male to physiological female were obtained in a 17β‐estradiol exposed group by sex‐linked markers (Hu et al. 2019). Nemesházi et al. (2020) developed a genetic sexing method based on sex‐linked SNP to study the sex reversal in Rana dalmatina, showing only female‐to‐male sex‐reversed adults. Among the wild‐caught R. dalmatina, XX/male ratio increased significantly with total anthropogenic land cover and agricultural activities (Nemesházi et al. 2020). According to the sex‐linked markers we screened, there were only XX‐males and no XY‐females in R. dybowskii were found in our research, as similar to the previous studies on R. temporaria (Alho et al. 2010) and R. dalmatina (Nemesházi et al. 2020). The sex reversal rate of R. dybowskii increased with the rearing temperature, to some degree, which may explain 2 sex‐linked markers showing decreasing levels of sex linkage with decreasing latitude over 800 km.

The effect of pseudo‐males on the phenotypic sex of offspring is different in studies conducted on fish. All‐female Oreochromis niloticus were obtained by mating pseudo‐males (XX) with genetic females (XX) (Pandit et al. 2015). In contrast, Chen et al. (2014) found that the offspring from pseudo‐males (ZW) of Cynoglossus semilaevis can spontaneously develop into functional pseudo‐males without environmental stimuli, which may be caused by the stable epigenetic inheritance pattern of the Z chromosome in the pseudo‐male gene. In this study, the proportion of female froglets from pseudo‐male parent was more than 95% at 15°C, which was similar to other studies that found offspring of R. temporaria lacking a Y haplotype showed an extreme bias towards females (Alho et al. 2010). The desiccated oviduct of the female R. dybowskii, Oviductus Ranae, is a valuable Chinese crude drug and is recorded in the Pharmacopoeia of the People's Republic of China, and we mainly obtain Oviductus Ranae by catching wild and semi‐wild R. dybowskii (Zhang et al. 2019). In confronting the current 2019‐nCoV disease epidemic, the Ministry of Agriculture and Rural Affairs (MARA) of China informed R. dybowskii shall be managed by the Department in charge of Fisheries in accordance with the requirements of aquatic animals in an attempt to curb the risk from diseases that jump from wildlife to humans. Consequently, our result could be of interest for all‐female R. dybowskii cultivation to avoid catching wild R. dybowskii and the protection of wild resources.

AUTHOR CONTRIBUTIONS

Yuan Xu, Jiayu Liu, Xiujuan Bai, and Zhiheng Du designed the study; Yuan Xu, Jiayu Liu, and Hang Su carried out the sampling and analyses; Shengwei Di, Fangyong Ning, Shiquan Cui, Lijuan Wang, Jianming Liu, and Chuanshuai Ren were involved in revising it critically for important intellectual content; and Yuan Xu and Zhiheng Du contributed to the writing of the paper which was checked and approved by all the authors.

Figure S1 The shapes of testis under anatomical microscope. (▲: testis, *: kidney)

Click here for additional data file.

Figure S2 The shapes of ovary under anatomical microscope. (◆: ovary, *: kidney)

Click here for additional data file.

Supplement 1 The sequence of sex‐linked molecular marker

Click here for additional data file.