Skip to main content

Analyzing the gonadal transcriptome of the frog Hoplobatrachus rugulosus to identify genes involved in sex development

Abstract

Background

The tiger frog (Hoplobatrachus rugulosus) is listed as a national Class II protected species in China. In the context of global warming, the sex ratio of amphibians will be affected, and the development of the population will be limited. Therefore, considering the potential for a decrease in the number of amphibians, studying sex evolution and molecular regulation of gonadal development in H. rugulosus, phenomenon that are currently unclear, is of great significance.

Results

Here, H. rugulosus was used to explore the mechanisms regulating gonadal development in amphibians. Illumina HiSeq 3000 was used to sequence the gonadal transcriptome of male and female H. rugulosus at two growth stages to identify genes related to gonadal development and analyze expression differences in the gonads. This analysis indicated that cyp17α, hsd3β, hsd11β1, cyp19α, and hsd17β12 perform vital functions in sex development in amphibians. Specifically, the expression of cyp3α, cyp17α, hsd3β, hsd11β1, sox2, sox9, sox30, soat, cyp19α, hsd17β12, and hspα1s was correlated with gonadal development and differentiation in H. rugulosus, as determined using the quantitative reverse transcriptase-polymerase chain reaction.

Conclusion

Significant differences were found in the gonadal gene expression levels in H. rugulosus of both sexes, and we identified a steroid hormone synthesis pathway in this species and analyzed related gene expression, but the changes during sex differentiation were still unclear. To our knowledge, this report presents the first analysis of the H. rugulosus gonadal transcriptome and lays the foundation for future research.

Peer Review reports

Introduction

Sex development in amphibians is a complicated process. As species that undergo genetic sex determination, the genotypic sex of amphibians is determined at the fertilization stage, but the phenotypic sex is influenced by the environment [1,2,3] and ultimately depends on the results of gonadal development. Therefore, the gonads represent the most appropriate organ for studying the sex of amphibians. Gonadal development in amphibians begins with the migration of primordial germ cells (PGCs) to the genital ridges, and then the PGCs and genital ridges co-develop into the primitive gonads, which further differentiate into testes or ovaries [4,5,6]. In most studies, the period from tadpole growth to the end of metamorphosis describes gonadal development [7,8,9], although a few studies have explored the complete progression from gonadogenesis to differentiation of the testes or ovaries and then to sexual maturation [10,11,12]. For example, Mali et al. [12] completely described the events that occurred at different stages of gonadal development in Microhyla ornata and Hylarana malabarica frogs and compared differences in their gonadal development during growth. Thus, previous research has established a relatively deep understanding of the processes involved in gonadal development in amphibians.

Modern biotechnology has made it possible to explore the molecular mechanisms underlying sex determination in amphibians. Many researchers have focused on genotypic sex determination systems and sex chromosomes of amphibians [13,14,15,16]. Some of the genes involved in sex development in amphibians, such as dmrt1 [17,18,19], sox3 [20], sox9 [21, 22], dax1 [23], sf1 [24], mis [25], and amh [19], were discovered over the course of a few years. Researchers have identified the relevant genes through molecular-biology–based approaches and gained a preliminary understanding of their functions and mechanisms of action [19, 20]. Further, the only sex-determining gene dm-w has been found in Xenopus laevis [26]. Some reports have also shown the differential expression of some sex-related genes during gonadal development in amphibians [17, 27, 28]. However, the overall differences in the related gene expression levels during gonadal development in amphibians remain unknown.

During primary gonadal differentiation in amphibians, steroid hormones secreted by the gonads can control the development of accessory structures, and by extension, secondary sexual characteristics, which ultimately affects the sex phenotype [29]. Because of the special sex development in amphibians, the effects of steroid hormones on sex differentiation have also received extensive attention. Early researchers treated amphibians with exogenous steroid hormones and found that the phenotypic sex was reversed [30,31,32], which demonstrated the importance of steroid hormones in phenotypic sex differentiation in amphibians. In addition, previous findings have also indicated that the effects of the same steroid hormone are not uniform in different species; even in the same species, the effects can vary with different doses [14, 33, 34]. Several studies have also revealed that steroid hormones play important roles in early gonadal development [14, 15, 33, 35]. However, because the molecular mechanism of action is unclear, the functions of steroid hormones during gonadal development are controversial [14, 35]. Previous studies showed that genes associated with steroid hormone synthesis, such as cyp11a1, star, hsd3b, cyp17, hsd17b, and cyp19 are expressed before gonadal differentiation [36]. Differential expression of cyp19 has been reported in Pleurodeles waltl tadpoles of different genotypic sexes when the original gonads would soon develop into testis or ovaries [37]; this demonstrated the regulatory role of steroid hormones in gonadal differentiation. However, relevant research is still limited, and the differential expression of these steroid-hormone-synthesis–related genes during the entire gonadal development needs to be further studied.

Substantial research has been conducted on the molecular mechanisms underlying sex development in amphibians, although sex-related gene expression changes occurring during gonadal development are unclear, as are the regulation of sexual steroid hormones in amphibians and associated molecular pathways during gonadal development. Therefore, many problems remain to be solved, and further studies are needed. The main obstacle to resolving these questions is the lack of genetic and genomic information for amphibian species, although this lack of information can now be overcome using next-generation high-throughput sequencing technologies [38, 39].

Hoplobatrachus rugulosus (Wiegmann, 1834) is a member of Dicroglossidae within the Anura order in the Amphibia class. This species is listed in Appendix II of CITES as a national Class II protected species in China and represents a common frog in farmlands in southern China [40]. To strengthen the protection of this species, many scholars have begun studying the food habits, acoustic characteristics, growth, development, and physiology of this species [41,42,43,44,45,46]. Currently, the increase in global temperature owing to the greenhouse effect is modifying the sex of amphibians [14]. Thus, in the context of climatic warming, the sex ratio of amphibians will be affected and the development of this population will be limited. Therefore, considering the potential for a decreasing number of amphibians, based on the particularities of amphibian sex determination, studying sex evolution and the molecular regulation of gonadal development in H. rugulosus is of utmost significance. In this study, we generated millions of sequence reads from the gonadal transcriptome of male and female H. rugulosus at two growth stages, for the first time, using the Illumina HiSeq 3000 platform. A non-redundant set of transcripts was generated, various analyses were performed, and the expression of sex-related– and steroid-hormone-synthesis–related genes was further examined to evaluate the expression differences in gonadal glands. The results of this study provide a molecular foundation for future research.

Results

Anatomical observations and histological analysis of testes and ovaries

Male and female frogs that were reared for 3 months after complete metamorphosis exhibited differences in body structure and size; the vocal sac structures were visible in the males, but the gonads were immature (Fig. 1a, b). After 15 months, the males and females exhibited evident sexual dimorphism, the frogs had grown to become breeding individuals, and the gonads were fully mature (Fig. 1c, d); therefore, male and female gonad samples were obtained from 3- and 15-month-old H. rugulosus by dissection, and they were then observed and weighed. The ovaries exhibited apparent changes in size and morphology as the development progressed, age increased, volumes of the ovaries increased, and the color changed from pinkish white to black. In contrast, only the size of the testes changed, with little differences in the morphological characteristics. The mass of the gonads in 3-month–old females (hereafter referred to as 3F) was larger than that in 3-month–old males (hereafter 3 M; one-way ANOVA; F1, 4 = 37.0, P < 0.01), and the gonad mass increased with the age of H. rugulosus (Fig. 1e). The mass of the gonads in 15-month–old females (hereafter 15F) was greater than that in 15-month–old males (hereafter 15 M; one-way ANOVA; F1, 4 = 34.4, P < 0.01). We also prepared sections for histological analysis, which revealed the presence of growing follicles in the 3F gonads under a microscope. At 15 months of age, the follicles became larger and tended to mature. However, primary spermatocytes appeared in the 3 M gonads, and sperm cells were already present. Numerous sperm cells were already present in the 15 M gonads, and sexual maturation was achieved (Fig. 2).

Fig. 1
figure1

Anatomic observations of gonads from Hoplobatrachus rugulosus and comparison of their masses at two ages. a Left and right gonads from a 3-month-old female. b Left and right gonads from a 3-month-old male. c Left and right ovaries from a 15-month-old female. d Left and right gonads from a 15-month-old male. e Comparisons of the masses of gonads from males and females at two developmental stages. The scale bar represents 1 cm in the figure

Fig. 2
figure2

Histological observations of gonads from both sexes of Hoplobatrachus rugulosus at two ages. a Ovary from a 3-month-old female. b Ovary from a 15-month-old female. c Testis from a 3-month-old male. (d) Testis from a 15-month-old male

Sequence analysis and assembly

cDNA libraries were constructed using testes and ovaries to study the transcriptome of H. rugulosus at two ages. The analysis began with preprocessing and ended with comparisons of differentially expressed genes (DEGs) and functional annotations between different samples. After strict quality assessments and data filtering, 91,416,254 (96.53%) and 69,155,150 (95.07%) high-quality reads were generated from the 3 M and 3F groups, and 80,734,092 (91.31%) and 81,576,146 (91.23%) high-quality reads were generated from the 15 M and 15F groups, respectively, and these reads were selected for further analysis. The GC percentage of the reads in the 3 M group was 45.92%, that in the 3F group was 47.92%, that in the 15 M group was 45.92%, and that in the 15F group was 46.19% (Table S1). In addition, the respective Q20 and Q30 values of the DNA bases in the 3 M group after filtration were 94.36 and 86.51%, respectively, whereas those in the 3F group were 93.18 and 84.64%, those in the 15 M group were 91.89 and 82.48%, and those in the 15F group were 91.84 and 82.31% (Table S2). Therefore, the overall sequencing quality was good, and the data met the requirements for subsequent analysis. The Trinity program [47] was used to splice short RNA-seq sequences to obtain the reference genes, and 49,149 non-redundant unigenes were obtained. The average length obtained was 1593 base pairs (bp) and the N50 length was 2455 bp, which was indicative of a good-quality assembly (Table S3). The number of genes expressed in each sample was counted. Through this analysis, 48,910 unigenes were obtained from the 3 M samples, accounting for 99.51% of the total, 42,502 (86.48%) unigenes were obtained from the 3F samples, 48,588 (98.86%) unigenes were obtained from the 15 M samples, and 45,289 (92.15%) unigenes were obtained from the 15F samples (Table S4).

Functional annotation and expression analysis

The obtained unigenes were analyzed using Nr, Swiss-Prot, Kyoto Encyclopedia of Genes and Genomes (KEGG), Clusters of Orthologous Groups (COG), and EuKaryotic Orthologous Groups (KOG) databases for functional annotation, and 20,942 unigenes were matched in at least one database. Among them, the most unigenes (20,925) matched with the Nr database, and the fewest unigenes (11,982) matched with the KOG database (Table S5). The annotation information is represented in a Venn diagram in Fig. S1. The E values of the best comparisons for all unigenes in the four major databases were counted and divided into five ranges (10–20 < E value ≤10–5; 10–50 < E value ≤10–20; 10–100 < E value ≤10–50; 10–150 < E value ≤10–100; 0 ≤ E value ≤10–150), and the number of genes in each range was counted. Nr is a well-known protein database, which enables the use of BLASTx to compare assembled unigene sequences with the Nr database to obtain the sequence with the best alignment result (lowest E value) of each unigene in the Nr database as the corresponding homologous sequence. In cases of equal E values, the first value is used. The number of homologous sequences was statistically compared to each species to identify the species to which the homologous sequences belong, and the top three species were all amphibians (Fig. S2). The homologous sequences with the highest similarity were from Nanorana parkeri, belonging to the same family as H. rugulosus, which validated our comparisons. In the KOG database (Fig. S3), 11,982 unigenes were grouped into 25 categories, and “signal transduction mechanisms” accounted for the highest proportion (5209, 43.47%), followed closely by “general function prediction only” (3814, 31.83%). ln contrast, “nuclear structure” (80, 0.67%) and “cell motility” (52, 0.43%) accounted for the lowest proportions, which suggested that most of the gonadal functional genes of H. rugulosus were used for signal transduction and basic cellular activities.

Next, we calculated and analyzed the gene expression levels in each sample, and generated a bar chart (Fig. S4A) showing the statistical correlations of DEGs among the samples. Pairwise comparisons of DEGs were performed using the four samples, revealing the following results. Comparison of the 3F and 3 M groups revealed the presence of 29,166 upregulated genes and 8579 downregulated genes in the 3 M group (Table S6). Comparison of the 15F and 15 M groups revealed the presence of 24,660 upregulated genes and 10,879 downregulated genes in the 15 M group (Table S7). Comparison of the 3F and 15F groups revealed the presence of 12,786 upregulated genes and 3074 downregulated genes in the 15F group (Table S8). Finally, comparison of the 3 M and 15 M groups revealed the presence of 3357 upregulated genes and 5847 downregulated genes in the 15 M group (Table S9). These differences might reflect the fact that the 15 M and 3 M samples were derived from testes, whereas the 15F and 3F samples were from ovaries, which represent different tissue types, thus resulting in increased numbers of DEGs. In addition, genes expressed in the testes were more upregulated and less downregulated than those expressed in the ovaries at both 3 and 15 months of age. Comparison of the ovary samples at 3 and 15 months of age indicated that more genes were upregulated and fewer genes were downregulated at 15 months, whereas the opposite results were obtained when analyzing the testes samples. The scatter plot in Fig. S4B reveals greater correlations between samples obtained from females (3F and 15F) and males (3 M and 15 M), compared with samples obtained from different sexes. The volcanic plots of the four samples more intuitively show the expression levels of DEGs in gonads from different sexes at two developmental stages (Fig. S4C).

Enrichment analysis of DEGs

We carried out Gene Ontology (GO)/KEGG pathway enrichment analysis of the DEGs between the different groups, and GO functions and KEGG pathways associated with the DEGs were determined. GO analysis is mainly divided into three categories, molecular function, cellular components, and biological processes. Through GO functional-enrichment analysis, we classified DEGs obtained by pairwise comparisons, and selected GO terms (and the corresponding GO numbers) that might be related to gonadal development. Comparison of the 3F group with the 3 M group revealed 23 GO terms related to gonadal development or function, among which the DEGs were significantly enriched in biological processes such as gonadal development (GO:0008406), development of primary sexual characteristics (GO:0045137), sex differentiation (GO:0007548), and developmental process involved in reproduction (GO:0003006), among others (Fig. 3). Comparison of the 15F group with the 15 M group revealed 25 GO terms related to gonadal development or function, although no DEGs were significantly enriched. Comparison of the 3F group with the 15F group revealed 20 GO terms related to gonadal development or function, with significant enrichments in biological processes like the development of primary sexual characteristics (GO:0045137), sex differentiation (GO:0007548), gonad development (GO:0008406), development of primary female sexual characteristics (GO:0046545), female gonad development (GO:0008585), female sex differentiation (GO:0046660), developmental process involved in reproduction (GO:0003006), reproductive structure development (GO:0048608), and reproductive system development (GO:0061458). Finally, upon comparing the 3 M group with the 15 M group, we found that 17 GO terms were related to gonadal development or function, but no significant enrichment was found in terms of the DEGs. Through significant pathway enrichment, the major biochemical metabolic pathways and signal transduction pathways in which DEGs are involved can be determined.

Fig. 3
figure3

Gene Ontology (GO) term enrichment analysis. a, b GO terms (related to biological processes) significantly enriched (P < 0.05) for the differentially expressed genes (DEGs) identified when comparing the 3F group with the 3 M group (a) or the 3F group with the 15F group (b). 3F: 3-month-old female, 3 M: 3-month-old male, 15F: 15-month-old female, 15 M: 15-month-old male

To determine active biological pathways during the sexual development of H. rugulosus, we carried out pathway enrichment analysis of the DEGs. Comparison of the 3F group with the 3 M group revealed that the DEGs were enriched in 340 KEGG pathways. Comparison of the 15F group with the 15 M group revealed DEGs that were enriched in 339 KEGG pathways. Comparison of the 3F with the 15F group revealed DEGs that were enriched in 334 KEGG pathways. Finally, comparison of the 3 M group with the 15 M group revealed DEGs that were enriched in 321 KEGG pathways. Some of the signaling pathways identified are essential for development, including steroid hormone biosynthesis (ko00140), ovarian steroidogenesis (ko04913), Gnrh signaling pathway (ko04912), cortisol synthesis and secretion (ko04927), estrogen signaling pathway (ko04915), steroid biosynthesis (ko00100), oocyte meiosis (ko04114), and progesterone-mediated oocyte maturation (ko04914). Further information regarding these predictive pathways might be useful for investigating the roles played by these DEGs during the gonadal development of H. rugulosus.

Quantitative reverse transcriptase PCR (qRT-PCR)

After analyzing the overall transcriptome data, we filtered out the important DEGs associated with sex development and steroid hormone synthesis, and these results were verified by qRT-PCR. We compared the qRT-PCR results for 11 genes related to sex development with the transcriptomics results (Fig. 4). Cyp3α, cyp17α, hsd3β, hsd11β1, sox2, sox9, and sox30 showed higher expression levels in males than in females, whereas soat, cyp19α, hsd17β12, and hspa1s showed the opposite results (one-way ANOVA, all P < 0.05). soat, cyp3α, cyp17α, and sox2 showed higher expression levels in females at 15 months of age than at 3 months of age, whereas the sox30 gene showed the opposite result, and hsd3β gene expression was lower in the 15 M group than that in the 3 M group (one-way ANOVA, all P < 0.05; Fig. 4). In addition, most KEGG pathway-related genes were related to steroid hormone synthesis. After sorting and simplifying these pathways, we obtained the pathway map shown in Fig. 5. We used coloring to mark the sex-related genes in the pathway map for different groups to present the reads per kilobase per million mapped reads (RPKM) data more intuitively.

Fig. 4
figure4

qRT-PCR analysis of 17 differentially expressed genes (DEGs) in male/female gonads from Hoplobatrachus rugulosus at two ages. Gapdh expression was detected as a reference gene to normalize the qRT-PCR data. Data are expressed as the mean + standard error. *P < 0.05; **P < 0.01; ***P < 0.001

Fig. 5
figure5

Differential expression of sex-related genes in male/female gonads from Hoplobatrachus rugulosus at two ages. The original pathways (ko00140, ko04913, ko04927, ko04915, ko00100, and ko04024) from KEGG database (https://www.kegg.jp/kegg/; Permission: 210365)

Discussion

High-throughput transcriptome sequencing is widely used in various molecular studies [48]. Compared to traditional sequencing technology, high-throughput sequencing technology using the Illumina platform offers unparalleled advantages such as high throughput, high accuracy, and low cost; thus, the platform has become a powerful tool for full transcriptome research. Our study represents the first attempt to study the gonadal transcriptome of H. rugulosus. In total, 49,149 unigenes with an average length of 1593 bp were obtained, and the transcriptomes of the gonadal glands at different periods in different sexes were compared. The results are of great significance for the further study of gonadal development and sex differentiation in H. rugulosus.

In this study, the four groups were compared in pairs. We identified far more DEGs when comparing both sexes than when comparing both ages. These results might have been obtained because although both the testis and the ovary develop from the same original gonads, the testes and ovaries (two different sex organs) form owing to the differential expression of sex-related genes during gonadal differentiation [19]. Previously, several sex-related genes were found in amphibians [19,20,21,22, 24], and the expression levels of these were different in the testes and ovaries. For example, the amh gene of the TGF family [49, 50] is expressed in differentiated testicular Sertoli cells, but it is also expressed in the ovaries. Therefore, it was proposed that its function is not limited to testicular development [19]. The dmrt1 gene is involved in male sex development [51, 52] and is highly expressed in male germ cells as well as somatic cells during amphibian development, but it is also expressed at low levels in some amphibian ovaries [19]. The sox9 gene, a member of the SRY-related HMG-box (sox) family, is expressed in both the testes and ovaries, but the expression level in the testes is higher than that in the ovaries [22]. The results of these studies identified significant differences in gonadal gene expression among amphibians of different sexes, consistent with our current results with H. rugulosus.

In addition, more DEGs were upregulated and fewer were downregulated in the testes than those in the ovaries, at both stages of development, in previous reports for other species. Similar results were found in studies of gene differences in the testes and ovaries of Portunus trituberculatus [53] and Oratosquilla oratoria [54], suggesting that more genes were highly expressed in the testes than in the ovaries and that they were involved in testicular development or the regulation of biological processes. However, few studies on this topic have been conducted in amphibians to date, and the gene expression levels in gonads from different sexes remain to be further explored. When comparing pairs of two ages, more genes were upregulated and fewer were downregulated in the 15F ovaries than in the 3F ovaries, whereas the opposite was true in the testes at two developmental stages. Therefore, we suspected that during development, as the ovaries of females progressed from an immature to mature stage, more genes were involved in ovarian development. However, in the testes, more genes were actively expressed during the early stage of development, and these genes were downregulated during the progression from immature to mature stages. There is a lack of data regarding gonadal gene expression during development in amphibians; thus, our research provides a more in-depth understanding of this aspect of development and provides a reference for further study of the dynamic gene expression changes that occur during gonadal development in amphibians.

GO enrichment analysis of DEGs between the different groups showed that terms such as gonadal development, development of primary sexual characteristics, sex differentiation, and developmental process involved in reproduction were enriched for different genes between the 3F and 3 M groups and between the 3F and 15F groups. These results suggest that these DEGs played important roles as potential regulatory factors during the early gonadal development of females. In addition, by performing pathway enrichment analysis, we generated pathway maps for gonadal development-related pathways, such as steroid hormone biosynthesis, ovarian steroidogenesis, the Gnrh signaling pathway, cortisol synthesis and secretion, the estrogen-signaling pathway, steroid biosynthesis, oocyte meiosis, and progesterone-mediated oocyte maturation. Most of these pathways were related to the synthesis of gonadal steroid hormones, which confirmed the important role of steroid hormones in the development and differentiation of gonads in H. rugulosus, whereas the formation of gonads is largely dependent on steroid hormones [55]. It was also previously found that testosterone (T) and estradiol (E2) in amphibian tadpoles are naturally synthesized before gonad differentiation, and the gonads directionally differentiate after binding to the corresponding sex hormone receptors to form complexes [15, 16], which provides a foundation for investigating the molecular mechanisms regulating steroid hormone expression in amphibians. To more clearly analyze the molecular mechanisms underlying sex differentiation in H. rugulosus, we simplified the pathway, screened sex-related GO terms enriched in some important DEGs in the pathway and verified the results by qRT-PCR.

After qRT-PCR verification of the genes screened for sex development and steroid synthesis, the expression trends of 17 genes were consistent with the results obtained upon transcriptome analysis, among which 11 genes were selected for further analyses. The remaining six genes were expressed differently in other pathways instead of the steroid synthesis signaling pathway, and there was no significant difference in the expression of upstream and downstream genes; their specific role and function are unknown, and therefore, we did not analyze them. Among them, cyp3α, cyp17α, hsd3β, hsd11β1, sox2, sox9, and sox30 showed higher expression levels in males than in females, whereas the opposite results were found for soat, cyp19α, hsd17β12, and hspα1s. Figure 5 shows that cyp17α, hsd3β, hsd11β1, cyp19α, and hsd17β12 play important roles in the steroid hormone synthesis pathway, and the functions of cyp17α and cyp19α have been confirmed in numerous studies [6, 16, 55, 56]. It is known that cytochrome P450 17-hydroxylase/17, 20 lyase (cyp17α) can promote the conversion of progesterone into dehydroisoandrosterone and that cytochrome P450 aromatase (cyp19α) can transform T into E2 (Fig. 5). Indeed, cyp17α gene expression is upregulated in male tadpoles before sex determination and maintained a high expression level, and cyp19α is highly expressed in the undifferentiated gonads of female tadpoles [36]. Data from several studies have shown that increased expression of aromatase is a major indicator of ovarian differentiation in non-mammalian vertebrates, such as birds, reptiles, and various amphibians [22]. Therefore, the importance of the cyp19α gene in amphibian sex differentiation is evident.

Based on these data, in combination with our current results, we propose that during amphibian sex differentiation, transcription factors promote cyp17α expression in H. rugulosus at the early stage, some of which develop into males owing to low cyp19α expression and high cyp17α expression, whereas some have high gonadal cyp19α gene expression, which promotes the differentiation of ovaries and the development of females. Figure 5 shows that follicle-stimulating hormone (fsh) and luteinizing hormone (lh) regulate the production of corresponding hormones by binding to specific receptors (fshr and lhcgr, respectively), and resulting in the activation of protein kinase A and promotion of sox9 expression. Previous findings have shown that sox9 is associated with the development of male gonads [21, 57]. Our results showed that sox9 expression was higher in male gonads than that in female gonads, which also demonstrated a correlation between sox9 and male sex development. Other genes with significant differences have not been reported in the literature for amphibians. However, sox2, sox30, and sox9 all belong to the sox family, and sox30 gene expression has been related to gonadal development in other species [58]. Therefore, its function might be similar to that of sox9, and further experimental studies are needed on this topic. In reptiles, the sex of most species is related to temperature [59,60,61], and heat shock protein (hsp) genes have also been considered relevant candidate genes, owing to their temperature-sensitive expression [62]. Hsp70 family genes are necessary for translocation and protein folding [63]. Previous results have shown that the hsp70α mRNA expression is significantly higher in ovaries than that in testes at the same age [62]. In this study, hspα1s (a member of the hsp70 family) was highly expressed in female gonads, and previous data indicated that the sex differentiation process in amphibians is also affected by temperature [14]. These findings suggest the importance of hsp70 in gonadal development in species for which sex is affected by temperature.

In addition, soat, cyp3α, cyp17α, sox2, and other genes showed higher expression in females at 15 months of age than at 3 months of age, whereas the opposite results were found for sox30. These findings indicate that soat, cyp3α, cyp17α, sox2, and other genes might play an important role in mature female gonads. As our study was conducted only at 3 and 15 months of age, we were unable to show a continuous change in development, and whether these genes play increasingly important roles during the continuous sex differentiation of females still needs to be tested. The sox30 gene is thought to be related to gonadal development, spermatogonial differentiation, and spermatogenesis [58], supporting our observations that sox30 did not play an obvious role in female gonads, and we hypothesized that its expression decreased as the female gonads matured. However, as our sample was only analyzed at two time points, and as we cannot infer the entire developmental process, and this speculation needs to be further tested. In addition, hsd3b expression was lower in males at 15 months of age than at 3 months of age. Some scholars have reported that the sex-related gene sf1 is dimorphologically expressed in amphibians of different sexes but that its expression is upregulated in females and downregulated in males, with increasing age [24]. These data suggest that the genes studied here might have also underwent dimorphism changes between both sexes and that their expression levels changed according to different needs during ontogenesis. However, few reports have described changes in the expression levels of sex-related genes during development, which highlights a need for further study and verification in the future.

Conclusions

Our experimental results supported the following conclusions. (1) Significant differences were found in the gonadal gene expression levels in H. rugulosus of both sexes, and compared with those in ovaries, more genes were highly expressed in the testes and were involved in testicular development or in regulating testes-related biological process. (2) The important roles played by steroid hormones in H. rugulosus gonad differentiation were confirmed, and the major regulatory pathways of steroids during sex development were integrated through enrichment with GO and KEGG analyses. Specifically, the crucial functions of cyp17α, hsd3β, hsd11β1, cyp19α, and hsd17β12 were identified, although our results do not yet indicate which transcription factors promoted the expression of these genes. (3) qRT-PCR revealed that the expression of cyp3α, cyp17α, hsd3β, hsd11β1, sox2, sox9, sox30, soat, cyp19α, hsd17β12, and hspα1s differed between both sexes and ages, and therefore, we speculate that these genes are closely related to gonadal development and differentiation in H. rugulosus. However, our paper is limited by certain shortcomings. As we only focused on two time points in the developmental process, we could not observe the continuous changes in gene expression during gonadal development, and thus, during development, it is unknown if more genes are involved in ovarian development as the females mature gradually. Moreover, in the testes, more genes were upregulated during the early stage of development, but it is unclear if the expression of these genes decreased gradually with testicular development. Further verification is needed.

Materials and methods

Ethics statement

The experimental procedures used in this study complied with the current laws related to animal welfare and research in China and were specifically approved by the Animal Research Ethics Committee of Lishui University (permit number ARECLSU201606001).

Sample collection and preparation

The female and male frogs used in this study were collected from our froggery at Lishui University, China. Three frogs of each sex were collected at 3 and 15 months of age, which covered two gonad development stages. In each group, the coefficient of variation of snout-vent length and the body mass of selected frogs ranged from 2.9 to 8.9% and from 3.5 to 12.7%, respectively, and the details are shown in the Table S10. After transfer to our laboratory, each frog was individually reared in a 500 × 360 × 200 mm plastic tank placed in a room for 1 week, where the temperature varied naturally, often within the range of 26 °C to 30 °C. After anaesthetization with 0.05% MS-222 (Sigma), the ovary and testis tissues were collected. Next, samples from frogs of the same sex and age were mixed together and transferred to 2 ml RNase-free plastic tubes. These frogs were derived from the same clutch of tadpoles, and they grew up in the same rearing environment, which could eliminate the differences within the same group caused by genetic and environmental backgrounds. Samples were immediately frozen in liquid nitrogen and stored at − 80 °C.

RNA extraction, library construction, and sequencing

Ovary (or testis) tissues from frogs at each age were subjected to RNA extraction. Total RNA was extracted using the Trizol Reagent Kit (Sangon Biotech, Shanghai, China) following the manufacturer’s protocol. RNA quality (RIN values ranged from 7.5 to 8.7) was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and confirmed using RNase-free agarose gel electrophoresis. Equal amounts of high-quality RNA from each ovary or testis sample (three samples were mixed for per age group) were then pooled for RNA-seq analysis. Oligo (dT) beads were used to purify mRNA from extracted total RNA samples, and the mRNAs were fragmented into short fragments using fragmentation buffer (NEB #E7490, New England BioLabs). The fragmented mRNA was used as the template, and the first strand cDNA was synthesized by random hexamers; the second strand cDNA was then synthesized by adding the buffer, dNTPs, RNase H, and DNA polymerase I. After purification with a QIAquick PCR Purification Kit and elution in EB buffer, the obtained products were suspended into End Repair Mix for end reparation and adenylation of the 3′ ends. Following addition of the sequencing adaptors, fragments with the targeted size (200 bp) were recovered by agarose gel electrophoresis (A620014, Agarose, Sangon Biotech, Shanghai), after which PCR amplification was performed to create the final cDNA library. The PCR products were sequenced using the Illumina HiSeq 3000 sequencing platform. Preparation and sequencing of the cDNA library were implemented by Genedenovo Biotechnology Co., Ltd. (Guangzhou, China).

De novo assembly of sequencing reads

The original imaging data produced were converted to sequence data using Base Calling software (bcl2fastq v2.20.0.422, Illumina), which we used to call raw data or raw reads, and cleanup was performed using Fastp software v0.18.0 [64]. In addition, the sequencing depth was 8G, paired-end sequencing was used, and the sequencing strategy was PE150. Not all clean reads obtained were valid, and reads containing adaptor, repetitive, or low-quality sequences could affect the assembly and subsequent analysis. Thus, we used Fastp to filter the clean reads by deleting adaptor sequences, reads where the proportion of Ns was > 10%, and low-quality reads (Q20 accounting for > 40% of the total reads). After recovering the filtered high-quality clean reads, the sequencing qualities were evaluated using Fastp. Subsequently, Trinity software v2.8.4 [47] was used for de novo transcriptome assembly. Reads with a 30 bp length of overlap were connected to form longer fragments with Trinity software (k-mer = 31), and these N-free assembled reads were assembled to generate unigenes.

Functional annotation and DEGs

The basic functional annotations of unigenes included protein functional annotations, pathway annotations, and functional annotation with the COG and KOG databases (http://www.ncbi.nlm.nih.gov/COG). First, the unigene sequences were aligned to protein databases such as the Nr (date of access: 2020-05-14), SwissProt (date of access: 2020-05-14), KEGG (version number: Release 93.0), and COG/KOG (date of access: 2015-07-24) databases (e value < 0.00001) using BLASTx with default parameters (June, 2019), the protein with the highest sequence similarity to a given unigene was selected, and the protein functional annotation information of the unigene was obtained. Then, the RPKM method [65] was used to calculate the expression levels of unigenes for four samples, using the following formula: RPKM = (1,000,000 × C)/([N × L]/1000). To calculate the expression level of unigene A, C represents the number of reads for unigene A, N is the total number of reads that uniquely aligned to all genes, and L is the number of bases in unigene A. The RPKM method can be used to eliminate the influence of the gene length and sequence-quantity differences on the calculated gene expression level. The calculated gene expression level can be directly used to compare the gene expression differences in different samples, using the edgeR program (http://www.bioconductor.org/packages/release/bioc/html/edgeR.html) for statistical analysis. The false-discovery rate (FDR) and log2 fold-change (FC) were obtained for each unigene, and DEGs were identified using screening criteria of FDR < 0.05 and |log2 FC| > 1.

GO enrichment analysis

We performed GO functional analysis of the DEGs identified among the samples [66]. The GO database employs three ontologies to describe the molecular functions, cellular components, and biological processes of genes, and GO analysis was performed using the DEGs from the edgeR analysis. The calculation was performed using Eq. 1, as follows:

$$ P=1-\sum \limits_{i=0}^{m-1}\frac{\left(\begin{array}{l}M\\ {}i\end{array}\right)\left(\begin{array}{l}N-M\\ {}n-i\end{array}\right)}{\left(\begin{array}{l}N\\ {}n\end{array}\right)} $$

where N is the number of genes with the GO annotation, n is the number of DEGs in N, M is the number of all genes that are annotated to certain GO terms, and m is the number of DEGs in M. Using an FDR of 0.05 as the threshold, the P value was calculated after the FDR was corrected, and then, a hypergeometric test was applied using the R v3.2.1 to identify the GO entries that were significantly enriched in DEGs, when compared with the whole genome background.

KEGG pathway enrichment analysis

The Pathway database is the main component of the public KEGG database (https://www.kegg.jp/kegg/) [67]. In organisms, different genes can exert biological effects in a coordinated manner, and pathway-based analysis is helpful for further understanding the biological functions of genes. KEGG pathway analysis was used to analyze datasets to identify significantly enriched genes, and a hypergeometric test was used to identify pathways that were significantly enriched for DEGs compared with the whole genomic background. Eq. 1 was also used for KEGG enrichment analysis. Specifically, N is the number of genes with a KEGG annotation, n is the number of DEGs in N, M is the number of specific pathways associated with all gene annotations, and m is the number of DEGs in M. Using FDR = 0.05 as the threshold, the P value was calculated after the FDR was corrected. Pathways that met this condition were defined as those with significant enrichment for DEGs.

qRT-PCR analysis

In this study, we also screened for DEGs important for sex development and steroid hormone synthesis and verified the expression levels of these genes by qRT-PCR [68]. For the qRT-PCR experiments, Trizol reagent (Sangon Biotech, China) was used to extract the total RNA, and the PrimeScript™ RT Reagent Kit with gDNA Eraser (RR047A, Takara, Japan) were used to reverse transcribe the total RNA into cDNA. qRT-PCR was performed using TB Green™ Premix Ex Taq™ (RR420A, Takara, Japan) and carried out in a CFX96 Touch Real-Time PCR Detection System (Bio-Rad, USA). The following conditions were used for amplification: 95 °C for 30 s, followed by 40 cycles of 95 °C for 5 s and 60 °C for 30 s. The sequences of the primers used for qRT-PCR are shown in Table S11. The amplicons were subjected to melt curve analysis to check the specificity of the amplified products. The relative expression level of each gene was calculated using the 2−(ΔΔCt) method [69], using the gene encoding glyceraldehyde-3-phosphate dehydrogenase (gapdh) as a reference gene, and the gene repeated on each board.

Availability of data and materials

The data was presented in the manuscript and the supporting materials. The raw reads data was submitted to the Short Read Archive (SRA) under the accession number SRR12516455, SRR12516456, SRR12516457 and SRR12516458 and BioProject accession number PRJNA659277.

Abbreviations

3 M:

3-Month–old males

15 M:

15-Month–old males

3F:

3-Month-old females

15F:

15-Month–old females

COG:

Clusters of orthologous groups

DEGs:

Differentially expressed genes

cyp17α :

Cytochrome P450 17-hydroxylase/17, 20 lyase

cyp19α :

Cytochrome P450 aromatase

E2:

Estradiol

FC:

Fold-change

FDR:

False-discovery rate

fsh :

Follicle-stimulating hormone

fshr :

Follicle-stimulating hormone receptor

gapdh :

Glyceraldehyde-3-phosphate dehydrogenase

GO:

Gene ontology

hsd :

Hydroxysteroid dehydrogenase

hsp :

Heat shock protein

KEGG:

Kyoto encyclopedia of genes and genomes

KOG:

EuKaryotic orthologous groups

lh :

Luteinizing hormone

lhcgr :

Luteinizing hormone/choriogonadotropin receptor

PGCs:

Primordial germ cells

qRT-PCR:

Quantitative reverse transcriptase PCR

sox :

SRY-related HMG-box

T:

Testosterone

References

  1. 1.

    Chardard DC, Penrad-Mobayed M, Chesnel A, et al. Thermal sex reversals in amphibians. In: Valenzuela N, Lance VA, editors. Temperature-Dependent Sex Determination in Vertebrates. Washington, DC: Smithsonian Books; 2004. p. 59–67.

    Google Scholar 

  2. 2.

    Eggert C. Sex determination: the amphibian models. Reprod Nutr Dev. 2004;44(6):539–49. https://doi.org/10.1051/rnd:2004062.

    Article  PubMed  Google Scholar 

  3. 3.

    Tompsett AR, Wiseman S, Higley E, Giesy JP, Hecker M. Effects of exposure to 17aethynylestradiol during larval development on growth, sexual differentiation, and abundances of transcripts in the liver of the wood frog (Lithobates sylvaticus). Aquat Toxicol. 2013;126:42–51. https://doi.org/10.1016/j.aquatox.2012.10.003.

    CAS  Article  PubMed  Google Scholar 

  4. 4.

    Flament S, Chardard D, Chesnel A, Dumond H. Sex Determination and Sexual Differentiation in Amphibians. In: Norris DO, Lopez KH, editors. Hormones and Reproduction of Vertebrates. UK: Academic Press; 2011. p. 1–19.

    Google Scholar 

  5. 5.

    Caneguim BH, Beltrame FL, da Luz JS, Valentini SR, Cerri PS, Sasso-Cerri E. Primordial germ cells (spermatogonial stem cells) of bullfrogs express sex hormone-binding globulin and steroid receptors during seasonal spermatogenesis. Cells Tissues Organs. 2013;197(2):136–44. https://doi.org/10.1159/000341517.

    CAS  Article  PubMed  Google Scholar 

  6. 6.

    Flament S. Sex reversal in amphibians. Sex Dev. 2016;10(5–6):267–78. https://doi.org/10.1159/000448797.

    CAS  Article  PubMed  Google Scholar 

  7. 7.

    Falconi R, Dalpiaz D, Zaccanti F. Ultrastructural aspects of gonadal morphogenesis in Bufo bufo (Amphibia Anura) 1. Sex differentiation. J Exp Zool. 2004;301A(5):378–88. https://doi.org/10.1002/jez.a.20069.

    Article  Google Scholar 

  8. 8.

    Mendoza-Cruz E, Moreno-Mendoza N, Zambrano L, Villagrán-SantaCruz M. Development and gonadal sex differentiation in the neotenic urodele: Ambystoma mexicanum. Zoomorphology. 2017;136(4):497–509. https://doi.org/10.1007/s00435-017-0361-z.

    Article  Google Scholar 

  9. 9.

    Phuge SK. Effect of fromestane on gonadal sex differentiation and sex ratio in the frog, Euphlyctis cyanophlyctis, with undifferentiated type of gonadal differentiation. J Herperol. 2018;52(2):171–5. https://doi.org/10.1670/17-019.

    Article  Google Scholar 

  10. 10.

    Haczkiewicz K, Ogielska M. Gonadal sex differentiation in frogs: how testes become shorter than ovaries. Zool Sci. 2013;30(2):125–34. https://doi.org/10.2108/zsj.30.125.

    Article  Google Scholar 

  11. 11.

    Phuge SK, Gramapurohit NP. Gonadal sex differentiation, development up to sexual maturity and steroidogenesis in the skipper frog, Euphlyctis cyanophlyctis. Gen Comp Endocrinol. 2013;181:65–71. https://doi.org/10.1016/j.ygcen.2012.08.017.

    CAS  Article  PubMed  Google Scholar 

  12. 12.

    Mali PV, Gramapurohit NP. Pattern of gonadal differentiation and development up to sexual maturity in the frogs, Microhyla ornate and Hylarana malabarica: a comparative study. J Exp Zool Part A. 2015;323(9):666–78. https://doi.org/10.1002/jez.1958.

    Article  Google Scholar 

  13. 13.

    Uno Y, Nishida C, Oshima Y, Yokoyama S, Miura I, Matsuda Y, et al. Comparative chromosome mapping of sex-linked genes and identification of sex chromosomal rearrangements in the Japanese wrinkled frog (Rana rugosa Ranidae) with ZW and XY sex chromosome systems. Chromosom Res. 2008;16(4):637–47. https://doi.org/10.1007/s10577-008-1217-7.

    CAS  Article  Google Scholar 

  14. 14.

    Nakamura M. Sex determination in amphibians. Semin Cell Dev Biol. 2009;20(3):271–82. https://doi.org/10.1016/j.semcdb.2008.10.003.

    Article  PubMed  Google Scholar 

  15. 15.

    Nakamura M. The mechanism of sex determination in vertebrates-are sex steroids the key-factor? J Exp Zool Part A. 2010;313(7):381–98.

    Article  Google Scholar 

  16. 16.

    Nakamura M. Is a sex-determining gene(s) necessary for sex-determination in amphibians? Steroid hormones may be the key factor. Sex Dev. 2013;7(1–3):104–14. https://doi.org/10.1159/000339661.

    CAS  Article  PubMed  Google Scholar 

  17. 17.

    Shibata K, Takase M, Nakamura M. The Dmrt1 expression in sex-reversed gonads of amphibians. Gen Comp Endocr. 2002;127(3):232–41. https://doi.org/10.1016/S0016-6480(02)00039-4.

    CAS  Article  PubMed  Google Scholar 

  18. 18.

    Osawa N, Oshima Y, Nakamura M. Molecular cloning of Dmrt1and its expression in the gonad of Xenopus. Zool Sci. 2005;22(6):681–7. https://doi.org/10.2108/zsj.22.681.

    CAS  Article  Google Scholar 

  19. 19.

    Piprek RP, Pecio A, Laskowska-Kaszub K, Kubiak JZ, Szymura JM. Sexual dimorphism of AMH, DMRT1 and RSPO1 localization in the developing gonads of six anuran species. Int J Dev Biol. 2013;57(11–12):891–5. https://doi.org/10.1387/ijdb.130192rp.

    CAS  Article  PubMed  Google Scholar 

  20. 20.

    Oshima Y, Naruse K, Nakamura Y, Nakamura M. Sox3: a transcription factor for Cyp19 expression in the frog Rana rugosa. Gene. 2009;445(1–2):38–48. https://doi.org/10.1016/j.gene.2009.05.011.

    CAS  Article  PubMed  Google Scholar 

  21. 21.

    Takase M, Noguchi S, Nakamura M. Two Sox9 messenger RNA isoforms: isolation of cDNAs and their expression during gonadal development in the frog Rana rugosa. FEBS Lett. 2000;466(2–3):249–54. https://doi.org/10.1016/S0014-5793(00)01078-4.

    CAS  Article  PubMed  Google Scholar 

  22. 22.

    Abramyan J, Feng CW, Koopman P. Cloning and expression of candidate sexual development genes in the cane toad (Bufo marinus). Dev Dynam. 2009;238(9):2430–41. https://doi.org/10.1002/dvdy.22055.

    Article  Google Scholar 

  23. 23.

    Sugita J, Takaseb M, Nakamura M. Expression of DAX1 during gonadal development of frog. Gene. 2001;208:67–74.

    Article  Google Scholar 

  24. 24.

    Mayer LP, Overstreet SL, Dyer CA, Propper CR. Sexually dimorphic expression of steroidogenic factor 1 (SF-1) in developing gonads of the American bullfrog, Rana catesbeiana. Gen Comp Endocr. 2002;127(1):40–7. https://doi.org/10.1016/S0016-6480(02)00019-9.

    CAS  Article  PubMed  Google Scholar 

  25. 25.

    Al-Asaad I, Chardard D, di Clemente N, Picard JY, Dumond H, Chesnel A, et al. Müllerian inhibiting substance in the caudate amphibian Pleurodeles waltl. Endocrinology. 2013;154(10):3931–6. https://doi.org/10.1210/en.2013-1229.

    CAS  Article  PubMed  Google Scholar 

  26. 26.

    Yoshimoto S, Okada E, Umemoto H, Tamura K, Uno Y, Nishida-Umehara C, et al. A W-linked DM-domain gene, DM-W, participates in primary ovary development in Xenopus laevis. Proc Natl Acad Sci U S A. 2008;105(7):2469–74. https://doi.org/10.1073/pnas.0712244105.

    Article  PubMed  PubMed Central  Google Scholar 

  27. 27.

    Suda M, Kodama M, Oshima Y, Yamamoto K, Nakamura Y, Tanaka S, et al. Up-regulation of FSHR expression during gonadal sex determination in the frog Rana rugosa. Gen Comp Endocr. 2011a;172(3):475–86. https://doi.org/10.1016/j.ygcen.2011.04.011.

    CAS  Article  PubMed  Google Scholar 

  28. 28.

    Haselman JT, Olmstead AW, Degitz SJ. Global gene expression during early differentiation of Xenopus (Silurana) tropicalis gonad tissues. Gen Comp Endocr. 2015;214:103–13. https://doi.org/10.1016/j.ygcen.2014.06.009.

    CAS  Article  PubMed  Google Scholar 

  29. 29.

    Mali PV, Gramapurohit NP. Are sex steroids essential for gonadal differentiation of the ornate frog, Microhyla ornata? Gen Comp Endocr. 2016;233:63–72. https://doi.org/10.1016/j.ygcen.2016.05.019.

    CAS  Article  PubMed  Google Scholar 

  30. 30.

    Zhang CY, Witschi E. Genic control and hormonal reversal of sex differentiation in Xenopus. P Soc Exp Biol Med. 1956;93(1):140–1.

    Article  Google Scholar 

  31. 31.

    Nishioka M, Miura I, Saitoh K. Sex chromosomes of Rana rugosa with special reference to local differences in sex-determining mechanism. Sci Rep Lab Amphibian Biol. 1993;12:55–81.

    Google Scholar 

  32. 32.

    Martyniuk CJ, Bissegger S, Langlois VS. Current perspectives on the androgen 5 alpha-dihydrotestosterone (DHT) and 5 alpha-reductases in teleost fishes and amphibians. Gen Comp Endocr. 2013;194:264–74. https://doi.org/10.1016/j.ygcen.2013.09.019.

    CAS  Article  PubMed  Google Scholar 

  33. 33.

    Piprek RP, Pecio A, Kubiak JZ, Szymura JM. Differential effects of testosterone and 17β-estradiol on gonadal development in five anuran species. Reproduction. 2012;144(2):257–67. https://doi.org/10.1530/REP-12-0048.

    CAS  Article  PubMed  Google Scholar 

  34. 34.

    Stephanie T, Beata RK, Maria O, Andreas L, Petros L, Frauke H, et al. Sex reversal assessments reveal different vulnerability to endocrine disruption between deeply diverged anuran lineages. Sci Rep. 2016;6(1):23825.

    Article  Google Scholar 

  35. 35.

    Hayes TB. Steroids as potential modulators of thyroid hormone activity in anuran metamorphosis. Amer Zool. 1997;37(2):185–94. https://doi.org/10.1093/icb/37.2.185.

    CAS  Article  Google Scholar 

  36. 36.

    Maruo K, Suda M, Yokoyama S, Oshima Y, Nakamura M. Steroidogenic gene expression during sex determination in the frog Rana rugosa. Gen Comp Endocrinol. 2008;158(1):87–94. https://doi.org/10.1016/j.ygcen.2008.04.019.

    CAS  Article  PubMed  Google Scholar 

  37. 37.

    Chardard D, Desvages G, Pican C, Dournon C. Aromatase activity in larval gonads of Pleurodeles watl (Urodele Amphibia) during normal sex differentiation and during sex reversal by thermal treatment effect. Gen Comp Endocrinol. 1995;99(1):100–7. https://doi.org/10.1006/gcen.1995.1089.

    CAS  Article  PubMed  Google Scholar 

  38. 38.

    Ansorge WJ. Next-generation DNA sequencing techniques. New Biotechnol. 2009;25(4):195–203. https://doi.org/10.1016/j.nbt.2008.12.009.

    CAS  Article  Google Scholar 

  39. 39.

    Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, et al. De novo transcript sequence reconstruction from RNA-seq using the trinity platform for reference generation and analysis. Nat Protoc. 2013;8(8):1494–512. https://doi.org/10.1038/nprot.2013.084.

    CAS  Article  Google Scholar 

  40. 40.

    Fei L, Ye CY, Jiang JP. Colored atlas of Chinese amphibians and their distributions. Chengdu: Sichuan Science and Technology Publishing House; 2012. p. 572–3.

    Google Scholar 

  41. 41.

    Wei L, Lin ZH, Ma XM, Zhao LH, Ma XH. Acoustic characteristics of the tiger frog, Hoplobatrachus rugulosus, during the breeding season. Zool Res. 2011;32(4):456–60. https://doi.org/10.3724/SP.J.1141.2011.04456.

    Article  Google Scholar 

  42. 42.

    Fan XL, Lei HZ, Lin ZH. Ontogenetic shifts in selected body temperature and thermal tolerance of the tiger frog, Hoplobatrachus chinensis. Acta Ecolpgica Sinica. 2012;32(17):5574–80.

    Article  Google Scholar 

  43. 43.

    Hu ZY, Lai YP, Chen WJ. Comparison of blood cells of Paa spinosa, Rana rugulosa and Rana nigromaculata. Sichuan J Zool. 2015;24(1):5–8. https://doi.org/10.3969/j.issn.1000-7083.2005.01.002.

  44. 44.

    Wang Y, Shao C, Wang XT. A study on relative fatness of tiger frog (Hoplobatrachus rugulosus) in Zhejiang province. Ecol Sci. 2008;27(6):490–3.

    Google Scholar 

  45. 45.

    Lin ZH, Ji X. Sexual dimorphism in morphological traits and food habits in tiger frogs, Hoplobatrachus rugulosus in Lishui, Zhejiang. Zool Res. 2005;26(3):255–62.

    Google Scholar 

  46. 46.

    Liu L, Liu CW, Lin D, Liang WF. Comparison of genetic diversity between frog Rana tigrina rugulosa and frog R .tigrina cantor. Fisheries Sci. 2005;24(12):31–3.

    Google Scholar 

  47. 47.

    Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-seq data without a reference genome. Nat Biotechnol. 2011;29(7):644–52. https://doi.org/10.1038/nbt.1883.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  48. 48.

    Morozova O, Hirst M, Marra MA. Applications of new sequencing technologies for transcriptome analysis: annual review of genomics and human genetics. Annu Rev Genomics Hum Genet. 2009;10(1):135–51. https://doi.org/10.1146/annurev-genom-082908-145957.

    CAS  Article  PubMed  Google Scholar 

  49. 49.

    Cate RL, Mattaliano RJ, Hession C, Tizard R, Farber NM, Cheung A, et al. Isolation of the bovine and human genes for müllerian inhibiting substance and expression of the human gene in animal cells. Cell. 1986;45(5):685–98. https://doi.org/10.1016/0092-8674(86)90783-X.

    CAS  Article  PubMed  Google Scholar 

  50. 50.

    Behringer RR, Finegold MJ, Cate RL. Müllerian-inhibiting substance function during mammalian sexual development. Cell. 1994;79(3):415–25. https://doi.org/10.1016/0092-8674(94)90251-8.

    CAS  Article  PubMed  Google Scholar 

  51. 51.

    Raymond CS, Shamu CE, Shen MM, Seifert KJ, Hirsch B, Hodgkin J, et al. Evidence for evolutionary conservation of sex-determining genes. Nature. 1998;391(6668):691–5. https://doi.org/10.1038/35618.

    CAS  Article  PubMed  Google Scholar 

  52. 52.

    Smith CA, Roeszler KN, Ohnesorg T, Cummins DM, Farlie PG, Doran TJ, et al. The avian Z-linked gene DMRT1 is required for male sex determination in the chicken. Nature. 2009;461(7261):267–71. https://doi.org/10.1038/nature08298.

    CAS  Article  PubMed  Google Scholar 

  53. 53.

    Meng XL, Liu P, Jia FL, Li J, Gao BQ. De novo transcriptome analysis of Portunus trituberculatus ovary and testis by RNA-Seq: identification of genes involved in gonadal development. PLoS One. 2015;10(6):e0128659. https://doi.org/10.1371/journal.pone.0128659.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  54. 54.

    Yan HW, Cui X, Shen XF, Wang LS, Jiang LN, Liu HY, et al. De novo transcriptome analysis and differentially expressed genes in the ovary and testis of the Japanese mantis shrimp Oratosquilla oratoria by RNASeq. Comp Biochem Phys D. 2018;26:69–78.

    CAS  Google Scholar 

  55. 55.

    Suda M, Uno Y, Fujii J, Matsuda Y, Nakamura M. Isolation and characterization of the CYP17A1 gene and its processed pseudogene in Rana rugosa. Comp Biochem Phys B. 2011b;160(1):54–61. https://doi.org/10.1016/j.cbpb.2011.05.008.

    CAS  Article  Google Scholar 

  56. 56.

    Oshima Y, Kato T, Wang D, Murakami T, Matsuda Y, Nagahama Y, et al. Promoter activity and chromosomal location of the Rana rugosa P450 aromatase (CYP19) gene. Zool Sci. 2006;23(1):79–85. https://doi.org/10.2108/zsj.23.79.

    CAS  Article  Google Scholar 

  57. 57.

    Trukhina AV, Lukina NA, Nekrasova AA, Smirnov AF. Sex inversion and epigenetic regulation in vertebrates. Russ J Genet. 2015;51(3):290–7.

    CAS  Article  Google Scholar 

  58. 58.

    Han F, Dong Y, Liu WB, Ma XX, Shi RS, Chen HQ, et al. Epigenetic regulation of Sox30 is associated with testis development in mice. PLoS One. 2014;9(5):e97203. https://doi.org/10.1371/journal.pone.0097203.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  59. 59.

    Pieau C, Dorizzi M, Richard-Mercier N. Temperature-dependent sex determination and gonadal differentiation in reptiles. Cell Mol Life Sci. 1999;55(7):887–900. https://doi.org/10.1007/s000180050342.

    CAS  Article  PubMed  Google Scholar 

  60. 60.

    Sarre SD, Georges A, Quinn A. The ends of a continuum: genetic and temperature-dependent sex determination in reptiles. BioEssays. 2004;26(6):639–45. https://doi.org/10.1002/bies.20050.

    Article  PubMed  Google Scholar 

  61. 61.

    Valenzuela N. Sexual development and the evolution of sex determination. Sex Dev. 2008;2(2):64–72. https://doi.org/10.1159/000129691.

    CAS  Article  PubMed  Google Scholar 

  62. 62.

    Kohno S, Katsu Y, Urushitani H, Ohta Y, Iguchi T, Guillette LJ. Potential contributions of heat shock proteins to temperature-dependent sex determination in the american alligator. Sex Dev. 2010;4(1-2):73–87. https://doi.org/10.1159/000260374.

    CAS  Article  PubMed  Google Scholar 

  63. 63.

    Gething MJ, Sambrook J. Protein folding in the cell. Nature. 1992;355(6355):33–45. https://doi.org/10.1038/355033a0.

    CAS  Article  PubMed  Google Scholar 

  64. 64.

    Chen SF, Zhou YQ, Chen YR, Gu J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34(17):i884–90. https://doi.org/10.1093/bioinformatics/bty560.

  65. 65.

    Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008;5(7):621–8. https://doi.org/10.1038/nmeth.1226.

    CAS  Article  PubMed  Google Scholar 

  66. 66.

    Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. Nat Genet. 2000;25(1):25–9. https://doi.org/10.1038/75556.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  67. 67.

    Kanehisa M, Furumichi M, Sato Y, Ishiguro-Watanabe M, Tanabe M. KEGG: integrating viruses and cellular organisms. Nucleic Acids Res. 2021;49(D1):D545–51. https://doi.org/10.1093/nar/gkaa970.

    CAS  Article  PubMed  Google Scholar 

  68. 68.

    Holzapfel B, Wickert L. Die quantitative Real-Time-PCR (qRT-PCR). Biol Unserer Zeit. 2007;37(2):120–6. https://doi.org/10.1002/biuz.200610332.

    CAS  Article  Google Scholar 

  69. 69.

    Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCT method. Methods. 2001;25(4):402–8. https://doi.org/10.1006/meth.2001.1262.

    CAS  Article  PubMed  Google Scholar 

Download references

Funding

This work was supported by grants from National Natural Science Foundation of China (31500308) and Zhejiang Provincial Natural Science Foundation of China (LQ16C040001), and Public Welfare Technology Application Research Project of Lishui City (2020GYX06).

Author information

Affiliations

Authors

Contributions

Y.T., J.Y.C., and G.H.D. did the statistical analysis and wrote the paper. G.H.D. and Z.H.L. designed and supervised the study, Y.T. and J.Y.C. captured and maintained animals. Y.T., J.Y.C. and G.H.D. collected data for the study. G.H.D. prepared figures and/or tables. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Guo-Hua Ding.

Ethics declarations

Ethics approval and consent to participate

The experimental procedures used in this study complied with the current laws related to animal welfare and research in China and were specifically approved by the Animal Research Ethics Committee of Lishui University (permit number ARECLSU201606001).

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Tang, Y., Chen, JY., Ding, GH. et al. Analyzing the gonadal transcriptome of the frog Hoplobatrachus rugulosus to identify genes involved in sex development. BMC Genomics 22, 552 (2021). https://doi.org/10.1186/s12864-021-07879-6

Download citation

Keywords

  • Assembly
  • Gene
  • Gonadal development
  • Hoplobatrachus rugulosus
  • Molecular pathway
  • Transcriptome