Genome-Wide SNP Analysis of Male and Female Rice Field Frogs, Hoplobatrachus rugulosus , Supports a Non-Genetic Sex Determination System

: Sex determination systems (SDSs) in anurans are diverse and have undergone independent evolutionary transitions among species. The mode of sexual reproduction of the rice field frog ( Hoplobatrachus rugulosus )—an economically viable, edible amphibian species—is not well known. Previous studies have proposed that threshold temperature conditions may determine sex in these frogs. To elucidate the SDS in H. rugulosus , we karyotyped 10 male and 12 female frogs, and performed fluorescence in situ hybridization combined with sequencing analyses using DArTseq™. Our results revealed a highly conserved karyotype with no sex chromosome heteromorphism, and the sequencing analyses did not identify any consistent sex-linked loci, supporting the hypothesis of temperature-dependent sex determination. The results of this study, and others, on SDSs in the rice field frog and related species also provide support for the theory that heteromorphic sex chromosomes may lead to an evolutionary trap that prevents variable SDSs. These findings add im-portant information to the body of knowledge on H. rugulosus and are likely to have a significant impact on the productivity and economic success of rice field frog farming.


Introduction
Sex determination systems (SDSs) in anurans involve both genotypic sex determination (GSD) and environmental sex determination (ESD) systems [1][2][3][4][5][6]. The investigation of GSD is difficult because many anurans have homomorphic sex chromosomes that cannot be easily identified using cytogenetic approaches [7,8]. Advanced high-throughput molecular methods have, therefore, been applied to determine genotype by sequencing, These frogs were originally obtained from a Northern Thai population that represents a single clade [27,47]. Adults were sampled with a standard weight of 200-250 g and length of 12.7-15.2 cm. The sex of each individual was determined based on morphological characteristics of which male and female are sexually distinct around 1 year and mature for mating [48][49][50] (Figure S1). The animals were euthanized in 0.25% (w/v) MS-222 (tricaine methanesulfonate). Blood samples were collected for DNA extraction, bone marrow for mitotic chromosome preparation. The sex of each individual was also confirmed by internal examination of gonadal morphology. All animal care and experimental procedures were approved by the Animal Experiment Committees of Chulalongkorn University (Protocol Review No. 1623002) and Kasetsart University (approval no. ACKU63-SCI-011), Thailand, and were conducted in accordance with the Regulations on Animal Experiments at Chulalongkorn University and Kasetsart University. Whole genomic DNA was extracted following the standard salting-out protocol as described previously, with slight modifications for extraction from different tissues [51]. The high-molecular-weight DNA samples were stored at −20 °C until required for DArTseq library construction, as described previously [10][11][12][13].

DArT Sequencing, Genotyping and Analysis
The DArTseq methodology for sequencing and genotyping by SNP loci was applied according to the protocol described by Jaccoud et al. (2001) [52]. Variability among SNP loci generates presence/absence polymorphisms in restriction sites, which are called presence-absence (PA) markers. Genotyping of multiple loci was performed using DArTseq™ (Diversity Arrays Technology Pty Ltd., Canberra, ACT, Australia) for SNP loci and in silico DArT to determine candidate sex-specific loci between male and female individuals. Approximately 100 ng of DNA from each sample was used for the development of DArTseq arrays. The DNA samples were subjected to digestion and ligation reactions, as described previously [10][11][12][13]53]. Sequences were processed using proprietary DArTseq analytical pipelines [54]. The outputs generated by DArTsoft14 were filtered based on reproducibility values (>3.5), average count for each sequence (sequencing depth > 5), balance of average counts for each SNP allele (>0.9), and call rate (>0.8) (proportion of samples for which the marker was scored as described previously [10][11][12][13]. Sex-linked or specific loci were obtained from the analysis of SNPs and PA markers. For an XX/XY sex determination system, the SNP and PA loci sequenced for 60%, 70%, 80%, 90%, and 100% of males were included in a separate data set. Loci that passed the 100% filtering criterion were designated as perfectly sex-linked, whereas those passing the 60-90% thresholds were considered moderately sex-linked loci, as described previously [10][11][12][13]. An opposite but similar approach was used to target loci based on a ZZ/ZW system. The Hamming distance was calculated to determine the number of combined loci between male and female individuals for pairwise differences in SNP and PA loci using the "rdist" function in R version 3.5.1. The Hamming distance represents the number of pairwise differences between all individuals across all loci. The Cochran-Armitage trend test (CATT) was performed to examine the genetic association between each locus and phenotypic sex from SNP and PA loci using the "catt" function of R version 3.5.1 with the HapEstXXR package. The CATT results were similar to those of a chi-square test that assessed whether the proportion of different genotypes followed the null expectation. Polymorphic information content (PIC), which is an index for evaluation of the informativeness of SNP and PA loci, was calculated for each locus and ranged from 0 (fixation of one allele) to 0.5 (frequencies of both alleles are equal) [10][11][12][13]55,56]. The probability of candidate sex-linked loci showing random associations with sex in a small sample size was estimated using the formula Pi = 0.5 n , where P is the probability for a given locus, i is sex-linked, 0.5 is the probability that either a female is homozygous or a male is heterozygous at a given locus, and n is the number of individuals sequenced at the locus, as described previously [10][11][12][13]

Chromosome Preparation
Mitotic chromosomes were obtained from bone marrow cells using an air-drying method. Briefly, after an abdominal cavity injection of 0.01% colchicine (Sigma, St. Louis, MI, USA) (0.7 mL per 100 g of frog weight) for 2 h, frogs were anesthetized by surgical anesthesia, and femurs and tibias were collected from each individual. The bone heads were cut just enough to insert a 23-gauge needle into the marrow cavity, and cells were flushed out into a conical centrifuge tube using a 1 cc syringe filled with 0.075 M KCl. After hypotonic treatment of bone marrow in 0.075 M KCl for 30 min at room temperature, the cells were collected by filtration using gauze and then fixed with 3:1 methanol/acetic acid. A drop of cell suspension was placed onto a clean glass slide and air-dried. The slides were stored at −80 °C until use. For karyotyping with conventional Giemsa staining, the chromosome slides were stained with 4% Giemsa solution (pH 7.2) for 10 min.

Microsatellite Repeat Motifs, Telomeric (TTAGGG)n FISH Mapping
An amplification of microsatellite repeat motifs is often observed with Y or W sex chromosomes in vertebrates [57][58][59][60][61][62]. To extensively identify sex chromosomes in jade perch, we performed telomeric repeats and 19 microsatellite repeat motifs using FISH mapping. Chromosomal locations of telomeric (TTAGGG)n sequences and 19 microsatellite repeat motifs were determined using FISH as previously described [59,63,64].  6, and (AAAAT)6. Fluorescence hybridization signals were captured using a cooled charge-coupled device camera mounted on a Nikon Eclipse 80 microscope and processed using NIS-Elements BR 3.2, software (Nikon Corporation, Tokyo, Japan).

Determination of Sex System and Identification of Sex-Linked Loci in Rice Field Frog
We sequenced 70,269 SNP loci and an additional 96,789 PA loci. PIC values ranged from 0.00 to 0.50 with an average of 0.33 for SNPs, and 0.01 to 0.50 with an average of 0.34 for PA markers. To determine whether GSD (XX/XY or ZZ/ZW) or ESD was the SDS in the rice field frog, we compared a number of SNP and PA loci by filtering using a gradually changing set of criteria. After filtering using male: female ratios of 40:60, 30:70, 20:80, 10:90, and 0:100, no sex-linked or sex-specific significant loci were associated with phenotypic sex in either the ZZ/ZW or the XX/XY GSD system ( Figure 1 and Tables 1 and 2).

Karyotype
More than 20 Giemsa-stained metaphase spreads were examined for each rice field frog. Metaphase analysis revealed chromosome numbers to be 2n = 26, comprising one pair of large metacentrics (first), one pair of large submetacentrics (second), two pairs of medium-sized submetacentrics (third and fourth), one pair of medium-sized metacentrics (fifth), four pairs of small-sized metacentrics (sixth, seventh, eighth, and ninth), and four pairs of small-sized submetacentrics (tenth, eleventh, twelfth, and thirteenth) (Figure 2).

Chromosomal Locations of the Telomeric (TTAGGG)n Sequences and Microsatellite Repeat Motifs
The results from FISH analysis revealed hybridization signals indicating the presence of TTAGGG repeats at the telomeric ends of all chromosomes, and interstitial signals were observed in seven chromosome pairs ( Figure S2). Hybridization signals for microsatellite repeat motifs of (AGAT)8 were detected in the subterminal region of the short arm of chromosome 1 in all males and females ( Figure 2). However, no signals were observed for the

Discussion
Many changes in the sex determination mode of anurans have been reported, making them a compelling focus for studies of dynamic SDSs [65,66]. Most Dicroglossidae, Ranidae, Mantellidae, and Rhacophoridae family members have highly conserved karyotypes with diploid chromosome numbers (2n) ranging from 22 to 26, with GSD. Our chromosome analyses of the rice field frog revealed that it has a highly conserved karyotype with 2n = 26, in line with previous reports [67][68][69][70][71][72], although slight karyotypic variation with different numbers of metacentric and submetacentric chromosomes was observed among different origins of specimens. This type of variation was also reported in the bullfrog (H. tigerinus, 2n = 26), which is a closely related species [73], possibly resulting from the presence of small inversions, translocations, or heterochromatin propagation in the lineages of the rice field frog and the bullfrog. Comparative chromosome mapping with FISH and functional cDNA or bacterial artificial chromosome probes is required to deduce the process of chromosomal rearrangements in these lineages [60,61,74]. No heteromorphism between male and female karyotypes was found in the rice field frog. These results suggested that rice field frogs may contain cryptic sex chromosomes. We, therefore, applied DArTseq™ technology to a large number of SNP and PA loci to enable the prediction of SDSs and sex-linked loci for the rice field frog using a sample size of 22 individuals (10 males and 12 females). No SNP or PA markers were identified for male-or female-linked loci across all specimens examined, indicating that the sex determination of the rice field frog is likely to be non-GSD. However, it should be noted that false-positive signals might be expected in such specimens because of their diverse genetic backgrounds [75]. The DArTseq™ technique has a few inherent deficiencies, such as low genome coverage and lack of prior information regarding the association of markers with targeted genes [76]. The phenotypic sex of anurans is easily affected by environmental factors such as temperature, endocrine-disrupting chemical pollutants, and genetic disturbances, such as hybridization through breeding between genetically different populations [77]. Previous studies have reported that the hatching success rate and survival rate of rice field frogs are significantly affected by temperature [30]. The proportion of males in a group rose to over 80% at 29 °C to 34 °C, suggesting that the gonads of rice field frog tadpoles are biased toward males at high temperatures [46]. The results from our genome-wide SNP analyses support the hypothesis of a TSD in rice field frogs. A similar pattern was also observed in Rana chensinensis, David, 1875 [78], Hong Kong rice-paddy frogs (Fejervarya multistriata, Hallowell, 1861) [79], and giant spiny frogs (Quasipaa spinose, David, 1875) [80], suggesting that high temperature promotes a bias toward males in most ESD anurans. Studies have shown that gonadal differentiation, locomotion, and growth in the rice field frog are not completely controlled by genes, while environmental factors such as temperature and hormones also affect gonadal differentiation to determine phenotypic sex [81]. However, the molecular mechanisms underlying TSD are difficult to identify. Genetic variants that alter the function of TSD genes are associated with gonadal phenotypes and allelic polymorphisms can affect biochemical pathways in response to temperature. Regulatory variants may also change the level of gene expression influenced by temperature [82]. Alternatively, if a single gene on a recombining chromosome determines sex, this gene is likely to be missed in a genome SNP analysis using reduced-representation approaches such as DArTseq™. Moreover, sex hormone dependence such as steroid hormones or androgens can also induce sex reversal in many frog species. In the rice field frog, the sex ratios treated with letrozole at 29 °C and 34 °C were significantly biased toward males, and male ratio increased as letrozole concentration increased [46]. It will be interesting to discover which of the three alternatives is most applicable to determine SDS in the rice field frog by conducting actual investigations.
Phylogenetic comparison of published SDSs in Dicroglossidae, Ranidae, Mantellidae, and Rhacophoridae, as well as our current data, raise questions about how and why temperature (a stochastic environmental factor), rather than GSD, influences the fate decision toward female or male differentiation in the rice field frog. Together, these findings support the hypothesis that sex chromosomes may form an evolutionary trap with respect to SDSs. Sex chromosomes undergo cycles of turnover by default unless the tipping point of differentiation is crossed. This establishes a heteromorphic sex chromosome trap, whereas homomorphic sex chromosomes retain the ability to turn over [83]. Anurans possess homomorphic sex chromosomes that appear to be evolutionarily young owing to their frequent turnover [84]. The transition from GSD to TSD or turnover to different GSD systems requires traversing a group of fitness-related genes, where individuals are produced carrying suboptimal or lethal WW or YY genotypes. However, several species appear to have escaped the sex chromosome evolutionary trap and evolved independently, such as pleurodonts and their sister group, corytophanids, which harbor different, partial-linkage sex chromosomal groups within their internal lineage [85,86]. The tendency for recurrence among sex chromosomal groups often results in homoplasy, and nearly homomorphic sex chromosomes are difficult to identify as X/Z and Y/W counterparts when using the Cbanding approach during chromosome analysis. Here, no C-positive heterochromatin was found in rice field frog chromosomes (data not shown). By contrast, large C-positive heterochromatins of the entire W sex chromosome in the bullfrog (H. tigerinus) were observed, while no heteromorphic sex chromosomes were found by conventional Giemsa staining [73]. Applications of advanced omics technology and molecular cytogenetics are necessary to further elucidate the status of sex chromosomes in many anurans. Highthroughput transcriptome sequencing of male and female gonad tissues may be another powerful approach to determine any key genes/differential expression patterns of genes that can be targeted by either sex in the rice field frog to elucidate the sex determination system in this frog species. Studying doses compensation events in rice field frogs (H. rugulosus) and other anurans may shed light on the pattern of sexual differentiation in this frog species. While there is a solid understanding of the evolutionary significance of TSD, the mechanistic basis of this SDS is still an enigma.

Conclusions
The results from our study support a TSD SDS in rice field frogs, as this species is considered economically viable and edible in several countries [5,7]. The many frog farms could consider improving production efficiency by manipulating frog sexes via temperature. A thorough examination of SDSs across Dicroglossidae, Ranidae, Mantellidae, and Rhacophoridae family members using the same approach is required to elucidate the tempo and mechanism of evolutionary transitions between modes of sex determination in anurans.

Supplementary Materials:
The following are available online at www.mdpi.com/article/10.3390/d13100501/s1, Figure S1: Morphological characteristics between (a) male and (b) female. Arrows indicate vocal sacs in male, Figure S2: FISH patterns of the telomeric (TTAGGG)n sequence on DAPI-stained metaphase chromosome spreads of male and female rice field frog (Hoplobatrachus rugulosus, Wiegmann, 1834).
Funding: This research was financially supported in part by the Graduate Scholarship Program of the Graduate School, Kasetsart University, Thailand awarded to T.P. and K.S. and by a Postdoctoral Researcher Grant award at Kasetsart University awarded to S.F.A. and K.S., the Office of the Ministry of Higher Education, Science, Research and Innovation; and the Thailand Science Research and Innovation through The Kasetsart University Reinventing University Program 2021 awarded to T.P. and K.S., and the e-ASIA Joint Research Program (e-ASIA JRP) (no. P1851311) awarded to K.S. and W.S. No funding source was involved in the study design, collection, analysis and interpretation of the data, writing the report and the decision to submit the article for publication.