Introduction

Xenopus laevis is a good model to study molecular mechanisms of gonad development because the structural changes in developing gonads and the master gene determining sex, the W-linked DM-domain gene (dm-w), are well known. The dm-w is located on W chromosome and thus is present only in the genetic females (ZW) (Yoshimoto et al. 2008). At the earliest stage of gonad development, the gonads are undifferentiated and bipotential. The expression of dm-w triggers ovary development, while its absence promotes testis development. It is believed that the DM-W protein blocks the DMRT1 (doublesex and mab-3-related transcription factor 1) involved in male sex determination (Yoshimoto et al. 2010). In addition to the dm-w, many other genes, which act independently or downstream of dm-w, are involved in the development of bipotential gonads into the ovaries or the testes (Piprek et al. 2016). However, the expression and role of many genes involved in gonadal development is still vague. At the initial stage of gonadogenesis (NF50, Nieuwkoop-Faber stage 50), the gonads consist of the gonadal cortex and the medulla. The gonadal cortex contains coelomic epithelium and the germ cells, which adhere to the interior face of the epithelium. The medulla is sterile and contains medullar cells only (Piprek et al. 2016, 2017). At this stage, the sex-determining genes (dm-w and dmrt1) are expressed in the somatic cells of the gonads. In the absence of dm-w, i.e., in the differentiating testis (ZZ), around stage NF53, the cortex and medulla fuse. Subsequently, around stage NF56, the germ cells become enclosed by the somatic cells, which results in the formation of testis cords (Piprek et al. 2017). The typical structure of the testis, i.e., fully differentiated testis cords separated by the interstitium, is established at stage NF62. In contrast, in differentiating ovaries, which express dm-w, the germ cells remain in the cortical position, and at stage NF56, the ovarian cavity forms inside the gonad. Around NF62, the ovaries are fully differentiated, with the oocytes located in the cortex (Piprek et al. 2017; Yoshimoto et al. 2008). This divergent development of the female and male gonads has to be controlled by differential gene expression. A global analysis of Xenopus gonad transcriptome, which we performed in this study, is the step in obtaining a broad database of gene expression pattern in developing male and female Xenopus gonads.

Among vertebrates, the transcriptome of developing gonads has been studied in the mouse (Beverdam and Koopman 2006; Chen et al. 2012; Gong et al. 2013; Jameson et al. 2012; Nef et al. 2005; Small et al. 2005), chicken (Ayers et al. 2015; Scheider et al. 2014), slider Trachemys scripta (Czerwinski et al. 2016), Alligator mississippiensis (Yatsu et al. 2016), and in several teleost fish species (Bar et al. 2016; Lin et al. 2017; Sreenivasan et al. 2008; Sun et al. 2018; Xu et al. 2016). These studies provided valuable insights into the genes involved in gonad development and identified new sex-determining gene candidates.

Among anurans, a transcriptome analysis was performed only in Silurana (Xenopus) tropicalis and only on already sexually differentiated gonads (from stage NF58) (Haselman et al. 2015). Thus, the genes expressed before and during the sexual differentiation of amphibian gonads are still unknown. The aim of our study was to examine the transcriptome of developing Xenopus gonads from the earliest stage of gonad development. We studied the gene expression pattern in four different stages of gonad development: the undifferentiated gonad during the period of sex determination (NF50), gonads at the onset of sexual differentiation (NF53), the differentiating gonads (NF56), and during the developmental progression of differentiated gonads (NF62) (Fig. 1).

Fig. 1
figure 1

Structural changes in developing gonads. a, b At stage NF50, there is no difference in the gonad structure between genetic sexes (ZW and ZZ). Such undifferentiated gonads (arrows) are composed of the somatic cells of coelomic epithelium (ce) covering the gonad, and germ cells (g) located inside; the germ cells are attached to the coelomic epithelium. The somatic cells gather in the gonad center forming gonadal medulla (m). At stage NF53, the first sexual differences appear in the gonad structure; in the differentiating ovaries (c, ZW), the germ cells remain in the peripheral position forming the ovarian cortex, whereas the centrally located medulla remains sterile. In the ZZ (male) gonads at the onset of sexual differentiation (d, the onset of the testis differentiation), the germ cells (g) detach from the coelomic epithelium and move towards the gonad center (medulla, m). At stage NF56, the differentiating ovaries (e) becomes compartmentalized into cortex and medulla; all germ cells (g) are located in the cortex and are attached to the coelomic epithelium; an ovarian cavity forms in the medulla (asterisk). In the differentiating testes (f), the germ cells (g) are dispersed and the cortex and medulla are absent. At stage NF62, the ovaries (g) contain large ovarian cavity (asterisk); the ovarian cortex contains meiotic cells (o). In the testes (h), the germ cells (g) are located within the testis cords (encircled). Scale bar, 25 μm

Full size image

Results and discussion

Sex-specific changes in the level of gene expression

In developing Xenopus laevis gonads (stages NF50, NF53, NF56, and NF62 combined), we detected the expression of 63,084 transcripts in total. We found that while the expression level of the majority of genes was similar between stages and between male and female gonads, a subpopulation of genes showed distinct changes in the expression level between stages and sexes, which suggested that they may play a role in sex determination and/or sexual differentiation (Figs. 2A, B and 3, Tables 1 and 2).

Fig. 2
figure 2

Diagram of changes in the number of genes upregulated and downregulated (≥ 2-fold change) between different stages in ZW gonads (a) and ZZ gonads (b)

Full size image
Fig. 3
figure 3

Diagram of changes in the number of genes with higher expression in ZW or ZZ gonads (≥ 2-fold change)

Full size image
Table 1 Number of genes with up- and downregulated (≥ 2-fold change) expression in ZW and ZZ gonads
Full size table
Table 2 Number of genes with up- and downregulated (≥ 2-fold change) expression in ZW versus ZZ gonads
Full size table

Analysis of gene expression level in the gonads showed that in the genetic females (ZW), the gonads at the onset of sexual differentiation (NF53) had 376 genes with upregulated expression and 1078 genes with downregulated expression in comparison to the undifferentiated gonad during sex determination period (NF50) (Fig. 2, Table 1). In the differentiating ovaries (NF56), only 143 genes were upregulated and 128 genes were downregulated in comparison to NF53 (Table 1). In differentiated ovaries (NF62), there were 918 genes with upregulated expression and 1834 genes with downregulated expression in comparison to NF56 (Table 1).

The genetic male (ZZ) gonads at the onset of sexual differentiation (NF53) had 659 genes with upregulated expression and 436 genes with downregulated expression in comparison to NF50 stage (Fig. 2, Table 1). In differentiating testes (NF56), 340 genes were up-, and 340 downregulated in comparison to NF53 stage. The differentiated testes at stage NF62 had 334 genes with upregulated expression and 831 genes with downregulated expression in comparison to NF56 stage.

Altogether, these data indicate that in both sexes, the transcriptional regulation is more robust during early gonadal development, i.e., at the onset of sexual differentiation of the gonad (NF50-NF53) and in the already differentiated gonads NF56-NF62 than in the differentiating gonads (NF53-NF56).

The comparison of gene expression level in between ZW and ZZ gonads showed significant differences between the sexes and revealed sexually dimorphic pattern of gene expression. At the initial phase of gonad development, i.e., in the undifferentiated gonads during sex determination (NF50), there were 1192 genes (i.e., 3.4%) with sexually dimorphic expression (≥ 2-fold change). Eight hundred twenty genes showed higher expression in ZW (genetic females), and only 372 showed higher expression in ZZ (genetic males) gonads (Fig. 3, Table 2). This indicates that female sex determination in Xenopus involves a robust transcriptional regulation. In contrast, in mice, during the sex determination period (between embryonic day E10.5 and E12.5), a higher number of genes were upregulated in the XY (genetic males) than in the XX (genetic females) gonads (Nef et al. 2005), which suggested that programs of sex determination may be diverse among vertebrates.

Our analysis showed that at NF53, i.e., at the beginning of sexual differentiation of Xenopus gonads, 1083 genes (i.e., 3%) showed sexually dimorphic expression (≥ 2-fold change), which was slightly lower number than at NF50 (during sex determination). One hundred ninety-three genes showed higher expression in ZW gonads, and 890 in ZZ gonads (Fig. 3, Table 2). Thus, at the onset of sexual differentiation, more genes were specifically expressed in ZZ (male) gonads than in ZW (female) gonads in Xenopus, which was opposite to the mouse, where more genes were specifically expressed in XX (female) than XY (male) gonads at the beginning of sexual differentiation (E13.5) (Nef et al. 2005). This again indicates differences in the molecular programs of gonad development among vertebrates.

At NF56, i.e., in the differentiating gonads, only 421 genes (i.e., 1.2%) showed sexually dimorphic expression (≥ 2-fold change). This stage showed the lowest percentage of genes with sexually dimorphic expression among all stages. Seventy-five genes had higher expression in ZW, and 346 in ZZ gonads (Fig. 3, Table 2). Thus, more genes were highly expressed in ZZ gonads (differentiating testes) than in ZW (differentiating ovaries). We previously showed that the testis differentiation in Xenopus is a complex process during which the basement membranes between gonadal cortex and medulla disintegrate, the cortex and medulla fuse, and the germ cells and somatic cells gather to form the testis cords (Piprek et al. 2017). This sequence of profound structural changes certainly requires an involvement of a number of different genes, which is reflected in the high number of genes expressed in ZZ gonads at this stage.

At stage NF62, the sexual dimorphism of gene expression is the most pronounced. At this stage, 3224 genes (i.e., 5%) showed sexually dimorphic expression (≥ 2-fold change). However, only 594 genes showed higher expression in ZW (ovaries), and as many as 2630 in ZZ (testes) gonads. This is the stage when the gonads of both sexes are already differentiated and fully prepared to perform their sex-specific functions, and therefore the sexual dimorphism is evident not only at structural but also at molecular level.

The expression of genes during different stages of ovary development

We found that in ZW gonads at stage NF53, in comparison to stage NF50, 376 genes had upregulated expression. The list of genes is presented in Suppl. Table 1, and chosen genes are presented in Table 3. Functional analysis grouped these genes in several distinct categories shown in Table 4. Among the upregulated genes, monoacylglycerol O-acyltransferase 2 gene 1 (mogat2.1) is involved in synthesis of diacylglycerol (DAG) that acts as a messenger lipid in cell signaling (Toker 2005); retinol-binding protein 2 (rbp2) is involved in retinoic acid regulation; extracellular proteins: collagen 2 and collagen 9, cysteine protease cathepsin K, epithelium-specific intermediate filaments: keratin 14 and keratin 19, estrogen receptor 1 (esr1), and synuclein gamma. At this early stage, the germ and somatic cells proliferate, and somatic cells start gathering in the gonad center forming medulla (Fig. 1A, C). Collagens accumulate between the gonad cortex and medulla (Piprek et al. 2017). Importantly, around stage NF50, a sex determination period takes place and gene expression analysis suggest that DAG, retinol, and estradiol may be involved in Xenopus sex determination.

Table 3 Chosen genes up- and downregulated in ZW (genetic females) gonads at NF53 in relation to NF50 stage
Full size table
Table 4 Number of genes assigned to functional groups up- and downregulated in ZW (genetic female) gonads
Full size table

We also found that in ZW gonads at stage NF53, there were 1078 genes with a downregulated expression in comparison to stage NF50. All these genes are listed in Suppl. Table 2, and chosen genes are presented in Table 3. Functional analysis grouped these genes in four categories shown in Table 4. Among these downregulated genes, there were signaling protein chordin (chrd), retinol-binding protein (rbp4), several protease inhibitors serpins, signaling proteins wnt10b and igf3 (insulin-like growth factor 3), transcription factors dmrt2, and mafb (Table 3).

In developing ZW gonad at stage NF56, in comparison to stage NF53, there were 143 genes with upregulated expression (Suppl. Table 3, and chosen genes are presented in Table 5). Functional analysis grouped these genes in three categories shown in Table 4. One of important genes upregulated in this period is a neurotrophin receptor a-1 (p75NTRa) (Table 2); its role in gonad development has never been studied; however, its upregulation suggests that neurotrophins (ligands of this receptor) can play a role in ovarian differentiation. We also found that in ZW gonad at stage NF56, in comparison to stage NF53, there were 128 genes with downregulated expression (Suppl. Table 4, and chosen genes are presented in Table 5). Functional analysis grouped these genes in several categories shown in Table 4. At NF56 stage, more genes responsible for reorganization of extracellular matrix and epithelial differentiation in ZW gonads are expressed than at stage NF53. Between stages NF53 and NF56, the medulla cells disperse, which results in the formation of the cavity in the ovary center (Fig. 1E). The mechanism of this event is not known and would be interesting to study how the neurotrophins, extracellular matrix, and epithelial differentiation are involved in this process.

Table 5 Chosen genes up- and downregulated in ZW (genetic females) gonads at NF56 in relation to NF53 stage
Full size table

In developing ZW gonad at stage NF62, in comparison to stage NF56, there were 918 genes with upregulated expression (Suppl. Table 5, and chosen genes are presented in Table 6). Functional analysis grouped these genes in the many categories (Table 4). Among known genes upregulated in the ovaries at stage NF62 are genes involved in meiosis and oocyte development, such as poly(A)-binding protein, oocyte-specific pou5f3.3, zygote arrest 1, zona pellucid proteins (zp2, zpd, zpy1), sycp3 (synaptonemal complex protein 3),and lhx8 (LIM homeobox 8). This reflects the onset of meiosis at stage NF62 and appearance of first oocytes (Fig. 1G). Also, more genes involved in the regulation of development, such as genes encoding the following: vegt protein, growth differentiation factor (gdf1), foxh1, foxr1, wnt11b, ddx25, and the survivin which prevents apoptosis, were upregulated at stage NF62 than at stage NF56.

Table 6 Chosen genes up- and downregulated in ZW (genetic females) gonads at NF62 in relation to NF56 stage
Full size table

In developing ZW gonad at stage NF62, in comparison to stage NF56, there were 1834 genes with downregulated expression (Suppl. Table 6, and chosen genes are presented in Table 6). Functional analysis grouped these genes into several categories (Table 4). Also, many (24) pathways were downregulated, including metabolic pathways, steroid hormone biosynthesis, retinol metabolism, PPAR signaling pathway, and adipocytokine signaling pathway (Table 4). Among known genes downregulated in the ovaries at stage NF62 are the following genes: retinol-binding protein 4 (rbp4), rdh16 (retinol dehydrogenase 16), several serpins, emx1.2 (empty spiracles homeobox 1), igf3 (insulin-like growth factor 3), igfbp1-a (insulin-like growth factor–binding protein 1), gata2 (GATA binding protein 2), and chordin. This indicates that retinol pathway and insulin-like growth factor pathway are downregulated at a later stage of ovarian development (NF62), and that these two pathways may be important for earlier stages of ovarian development. The PPAR signaling pathway and adipocytokine signaling pathway are involved in fat tissue differentiation (Ogunyemi et al. 2013) and are probably important for the development of corpora adiposa (fat tissue) at the anterior edges of the developing gonads at stages before NF62. Thus, after the fat tissue had been formed, these pathways are downregulated at stage NF62.

Another interesting gene expressed at the onset of gonadogenesis (NF50), showing upregulation at NF53 and downregulated at NF62 is chordin (chrd). Several studies showed that this gene is crucial for early organogenesis (dorsalization, gastrulation, and head development (Pappano et al. 1998; Bachiller et al. 2000), but its role in gonad development is unknown. Overall, our gene expression analysis showed that the later development of the ovary (NF62) is a very transcriptionally active period (many genes become upregulated and downregulated between NF56 and NF62), which may be related to the initialization of meiosis and oocyte formation during this developmental period.

The expression of genes during different stages of testis development

Our analysis showed that in the genetic male (ZZ) gonads at stage NF53, i.e., at the beginning of sexual differentiation, there were 659 genes with upregulated expression in comparison to the stage NF50 gonad (Suppl. Table 7, and chosen genes are presented in Table 7). Functional analysis grouped these genes into several categories (Table 8). There were the following genes with known function: igf3 (insulin-like growth factor 3), rbp4 (retinol-binding protein 4), vtn (vitronectin), several serpins, esr2 (estrogen receptor 2), several components of extracellular matrix (collagen 9, matrilin 2), and extracellular matrix (timp3) enzymes. A role of these genes in the early phase of ZZ gonad development is not known, and it would be interesting to study if retinol and/or igf3 are involved in male sex determination in Xenopus. Upregulation of PPAR and adipocytokine signaling pathways, characteristic for fat tissue, possibly reflects the onset of the development of the fat bodies at the anterior edge of the gonad.

Table 7 Chosen genes up- and downregulated in ZZ (genetic males) gonads at NF53 in relation to NF50 stage
Full size table
Table 8 Number of genes assigned to functional groups up- and downregulated in ZZ (genetic male) gonads
Full size table

Our analysis also showed that in the genetic male (ZZ) gonads at stage NF53, there were 436 genes with downregulated expression (Suppl. Table 8, and chosen genes are presented in Table 7). Functional analysis grouped these genes in the categories shown in Table 8.

Comparison of gene expression level in the ZZ gonads between stage NF56 and NF53 showed that at stage NF56, there were 340 genes with upregulated expression (Suppl. Table 9, and chosen genes are presented in Table 9). Functional analysis grouped these genes in categories shown in Table 8. Some of these upregulated genes are rbp2 (retinol-binding protein 2), receptor of prostaglandin E (ptger3), stromal cell-derived factor 2-like 1 (sdf2l1), and neurotrophin receptor (p75NTRa). Further, studies are necessary to establish what is the exact role of the prostaglandin E, retinol, and neurotrophins in testis differentiation. Importantly, around NF53-NF56, the cortex and medulla fuse in differentiating testes, and the germ cells lose their connection with the superficial coelomic epithelium and disperse in the whole testis (Fig. 1F). There were also 340 genes downregulated at stage NF56 ZZ gonad in comparison to stage NF53 gonad (Suppl. Table 10, and chosen genes are presented in Table 9). Functional analysis grouped these genes into several categories (Table 8).

Table 9 Chosen genes up- and downregulated in ZZ (genetic males) gonads at NF56 in relation to NF53 stage
Full size table

Comparison of gene expression level in the ZZ gonads between stages NF62 and NF56 showed that at stage NF62 gonad, there were 334 genes with the upregulated expression (Suppl. Table 11, and chosen genes are presented in Table 10). Functional analysis grouped these genes into several categories (Table 8). Around stage NF56-NF62, cells group into the testis cords (Fig. 1H). Genes involved in this process are not known, and presumably, the genes upregulated at this stage may be responsible for the formation of testis cords. There were also 831 genes downregulated in ZZ gonad at stage NF62 in comparison to stage NF56 (Suppl. Table 12, and chosen genes are presented in Table 10, and the gene categories, which were analyzed functionally are shown in Table 8).

Table 10 Chosen genes up- and downregulated in ZZ (genetic males) gonads at NF62 in relation to NF56 stage
Full size table

Genes with sexual dimorphism of expression in ZW and ZZ gonads in different developmental stages

The master sex-determining gene in Xenopus the dm-w was discovered in 2008 (Yoshimoto et al. 2008), but the molecular machinery of sex determination is certainly very complex and contains many other genes. We previously published the expression profile of known genes involved in sex determination and sexual differentiation in the Xenopus gonads (Piprek et al. 2018). We showed that the gata4, sox9, dmrt1, amh, fgf9, ptgds, pdgf, fshr, and cyp17a1 had upregulated expression in testes, while dm-w, fst, foxl2, and cyp19a1 had upregulated expression in the ovary (Piprek et al. 2018).

Here, we compared gene expression level between ZW and ZZ gonads at different stages of gonad development. These analyses showed that at stage NF50 (undifferentiated gonads during sex determination period), there were 820 genes with upregulated expression in ZW gonad (Suppl. Table 13, and chosen genes are shown in Table 12). Functional analysis grouped these genes into several categories (Table 11). Many genes upregulated in this period are uncharacterized. Among known genes upregulated in ZW gonad at stage NF50 is chordin (chrd). Chordin is a secreted protein responsible for several developmental processes such as dorsalization, head development, and gastrulation (Sasai et al. 1994; Pappano et al. 1998; Bachiller et al. 2000); our study indicates that it may play a crucial role in female sex determination (Table 12, Suppl. Table 13). Other genes upregulated in ZW gonad at NF50 are two protease inhibitors, serpin A3 and serpin I2, extracellular glycoprotein vitronectin, metalloproteinases mmp7 and adam27, retinol-binding protein rbp4, signaling molecules wnt10b, wnt11b, and igf3, helicase ddx25, and transcription factors foxa2 and lhx8. A role of these factors in sex determination in Xenopus is unknown and requires further study.

Table 11 Number of genes assigned to functional groups expressed at higher level in ZW and ZZ gonads
Full size table
Table 12 Chosen genes upregulated in ZW (genetic females) in relation to ZZ (genetic males) gonads at NF50 stage [higher gene expression level in ZW than in ZZ gonads]
Full size table

There were 372 genes with higher expression level in the ZZ (genetic males) gonads at stage NF50 (Suppl. Table 14, and chosen genes are shown in Table 13, and the functional groups are shown in Table 11). Among these upregulated genes are known genes such as epithelium markers keratin 5, 12, and 14, coiled-coil domain containing 50 (ccdc50) that acts as an effector in EGF signaling and negative regulator of NF-kB factor (Tsukiyama et al. 2012), signaling molecules: wnt3a, wnt7b, growth differentiation factor 3 (gdf3), fibroblast growth factor–binding protein 1 (fgfbp1), proteases cathepsin K and H, extracellular matrix molecules lumican, collagen IX and I, and decorin. A role of these genes in male sex determination and early testis development remains unknown.

Table 13 Chosen genes downregulated in ZW in relation to ZZ gonads at NF50 stage [higher gene expression level in ZZ than in ZW gonads]
Full size table

There are 193 genes with a higher expression in ZW (genetic females) gonad at stage NF53 (the onset of sexual differentiation of gonads) (Suppl. Table 15, and chosen genes are shown in Table 14). Functional analysis did not link these genes to any specific pathway. Among these upregulated genes, there are the following known genes: retinol-binding protein 2 (rbp2), protease calpain 8, synuclein gamma with unknown function, cell adhesion gene claudin 6, metalloproteinases mmp1 and adam21, and galectin-la involved in cell adhesion and signaling.

Table 14 Chosen genes upregulated in ZW in relation to ZZ gonads at NF53 stage [higher gene expression level in ZW than in ZZ gonads]
Full size table

There were 890 genes with higher expression in ZZ (genetic males) gonad at stage NF53 (the onset of sexual differentiation of gonads) (Suppl. Table 16, and chosen genes are shown in Table 15). Functional analysis grouped these genes into several categories (Table 11). The upregulated known genes are coiled-coil domain containing 50 (ccdc50), retinol-binding protein 4 (rbp4), signaling molecules igf1 and igf3, estrogen receptor 2 (esr2), transcription factors, Kruppel-like factor 9 (klf9), Kruppel-like factor 15 (klf15), and foxo1 (forkhead box O1), enzyme glycerophosphodiester phosphodiesterase 1 (gde1) responsible for synthesis of signaling molecule lysophosphatidic acid (LPA), cell adhesion proteins gap junction protein alpha 3 (gja3), occluding (ocln), and extracellular matrix component vitronectin (vtn).

Table 15 Chosen genes downregulated in ZW in relation to ZZ gonads at NF53 stage [higher gene expression level in ZZ than in ZW gonads]
Full size table

There were 75 genes with higher expression in ZW (genetic females) gonad at stage NF56 (Suppl. Table 17, and chosen genes are shown in Table 16, and the functional groups are shown in Table 11). Among known genes are retinoic binding protein 4 and vitronectin.

Table 16 Chosen genes upregulated in ZW versus ZZ gonads at NF56 stage [higher gene expression level in ZW than in ZZ gonads]
Full size table

There were 346 genes with higher expression in ZZ (genetic males) gonad at stage NF56 (Suppl. Table 18, and chosen genes are shown in Table 17, and the functional groups are shown in Table 11). Among known genes are keratin 14 and 15, cell molecule gap junction protein, alpha (gja3), endophilin B2 (sh3glb2) and coiled-coil domain containing 50 (ccdc50).

Table 17 Chosen genes downregulated in ZW in relation to ZZ gonads at NF56 stage [higher gene expression level in ZZ than in ZW gonads]
Full size table

There were 594 genes with higher expression in ZW (genetic females) gonad at stage NF62 (Suppl. Table 19, and chosen genes are shown in Table 18, and the functional groups are shown in Table 11). Many genes expressed at this stage such as zona pellucida glycoprotein 4 (zp4) and zona pellucida C glycoprotein (xlzpc) are involved in ovarian follicles and oocytes formation and development. Other genes with upregulated expression at this stage were enzyme arachidonate 12-lipoxygenase 12R type (alox12b) responsible for metabolism of a signal compound—arachidonic acid (ARA), signaling factors such as growth differentiation factor 1 (gdf1), Wnt11b, cell adhesion molecules claudin 6 and connexin 38, transcription factors foxr1 and foxh1, and survivin—an inhibitor of apoptosis.

Table 18 Chosen genes upregulated in ZW versus ZZ gonads at NF62 stage [higher gene expression level in ZW than in ZZ gonads]
Full size table

There were 2630 genes with upregulated expression in ZZ (genetic males) gonad at stage NF62 (Suppl. Table 20, and chosen genes are shown in Table 19). Functional analysis grouped these into many categories (Table 11). Among known genes with upregulated expression were factors involved in signaling and signaling pathways: igf1, desert hedgehog (dhh), sonic hedgehog (shh), indian hedgehog (ihh), wnt3a, wnt8b, wnt7b, Janus kinase 2 (jak2), frizzled receptor 4 and 10 (fzd4, fzd10), cellular retinoic acid–binding protein 2 (crabp2), SMAD family member 4 (smad4); proteases: serine protease 3 (prss3), cathepsin H (ctsh), peptidase inhibitor—serpini2; transcription factors: LIM homeobox 1 (lhx1), homeobox a9, d10, and d13 (hoxa9, hoxd10, hoxd13), foxf1, foxa2, gata2; extracellular matrix components: collagen III (col3a1), collagen I (col1a1), fibrillin 3 (fbn3); extracellular matrix enzymes: mmp2, mmp16, cell adhesion molecule 3 (cadm3); and intermediate filaments: keratin 15 and nestin (nst).

Table 19 Chosen genes downregulated in ZW in relation to ZZ gonads at NF62 stage [higher gene expression level in ZZ than in ZW gonads]
Full size table

Genes identified here that showed sexual dimorphism of expression can be categorized into several functional groups: (1) signaling molecules: chordin (upregulated in ♀), wnt3a (upregulated in ♂), wnt7b (♂), wnt8b (♂), wnt10b (♀), wnt11b (♀), igf1 (♂), igf3 (♀ and ♂), gdf1 (♀), gdf3 (♂), ccdc50 (effector in EGF pathway) (♂), including hedgehog factors (♂): dhh, shh, ihh; (2) retinoic binding proteins: rbp2 (♀), rbp4 (♀ and ♂); (3) enzymes involved in signaling: enzyme glycerophosphodiester phosphodiesterase 1 (gde1) responsible for synthesis of signaling molecule lysophosphatidic acid (LPA) (♂), enzyme arachidonate 12-lipoxygenase 12R type (alox12b) responsible for metabolism of a signal compound—arachidonic acid (♀); (4) receptors of wnt signaling: fzd4 (♂), fzd10 (♂); (5) proteases: cathepsin H (♂), cathepsin K (♂), calpain 8 (♀); (6) protease inhibitors: serpin A3 (♀), serpin I2 (♀ and ♂); (7) transcription factors: foxa2 (♀), foxf1 (♂), foxh1 (♀), foxo1 (♂), foxr1 (♀), lhx1 (♂), lhx8 (♀), gata2 (♂), Kruppel-like factor 9 (klf9) (♂), Kruppel-like factor 15 (klf15) (♂); (8) helicase: ddx25 (♀); (9) cell adhesion molecules: occludin (♂), claudin 6 (♀),galectin-a (♀); (10) extracellular matrix components (mainly in ♂): collagens 1,3,9 (♂), vitronectin (♂), decorin (♂), lumican (♂), fibrillin 3 (♂); (11) extracellular matrix enzymes: mmp1 (♀), mmp2 (♂), mmp7 (♀), mmp16 (♂), adam21 (♀), adam27 (♀); (12) oocyte-specific proteins (♀): zp4, xlzpc; (13) epithelium-specific intermediate filaments (♂): keratins 5, 12, 14, 15.

The changes in the level of the expression of several genes listed above indicate that EGF signaling and lysophosphatidic acid (LPA) signaling may be involved in testis differentiation, arachidonic acid signaling may be involved in ovarian differentiation, while the wnt signaling, insulin-like growth factor signaling, and retinol signaling may be involved in gonad development in both sexes.

Interestingly, from the moment of sexual differentiation (after stage NF53), the genes encoding cytoplasmic and nuclear proteins are upregulated in ZW gonads (developing ovaries), while the genes encoding cell membrane proteins are upregulated in ZZ gonads (developing testes) (Fig. 4). The same trend was noted during gonad development in Silurana tropicalis (Haselman et al. 2015). This indicates that there are important molecular differences between developing ovaries and testes.

Fig. 4
figure 4

Subcellular distribution of gene products (obtained from the Ingenuity Pathway Analysis)

Full size image

Comparison of sex-specifically expressed genes in developing gonads of Xenopus and other vertebrates

We compared Xenopus microarray data to the published microarray data of developing gonads in other vertebrates: mouse (Jameson et al. 2012), chicken (Ayers et al. 2015), a red-eared slider Trachemys scripta (Czerwinski et al. 2016), American alligator (Yatsu et al. 2016)—both species with temperature-dependent sex determination, and zebrafish (Sreenivasan et al. 2008). The comparison is shown in Tables 20, 21, and 22.

Table 20 Comparison of sex-specifically expressed genes in developing gonads of Xenopus and mouse
Full size table
Table 21 Comparison of sex-specifically expressed genes in developing gonads of Xenopus and chicken
Full size table
Table 22 Comparison of sex-specifically expressed genes in developing gonads of Xenopus and red-eared slider (Trachemys scripta), American alligator, and zebrafish
Full size table

The transcriptome of developing mouse gonad did not show the expression of Wnt3, Wnt7, Wnt8, Wnt10, Wnt11, and chordin (Jameson et al. 2012), which were expressed in Xenopus developing gonads. The Igf1 was expressed in XX (genetic females) mouse gonads at a higher level than in XY gonads (Jameson et al. 2012); however, in Xenopus, this gene was expressed in ZZ developing gonads (genetic males). In mouse, in contrast to Xenopus (data presented in this study), the developing gonads did not express the Igf3, Gdf1, and Gdf3 (Jameson et al. 2012). The Ccdc50 was expressed in the developing mouse gonads but did not show sexual dimorphism of expression (Jameson et al. 2012). In Xenopus, this gene had an upregulated expression in ZZ gonads. Among hedgehog growth factors, in developing mouse gonads, only the dhh was expressed (Jameson et al. 2012). In Xenopus, gonads dhh and also shh and ihh were expressed. In mice, the Rbp1 (in XX) and Rbp4 (in XY gonads) were expressed (Jameson et al. 2012). In Xenopus, the rbp2 was expressed in ZW and rbp4 in ZZ and ZW gonads. Gde1 gene was expressed in developing mouse gonads; however, it did not show sexual dimorphism of expression (Jameson et al. 2012) In Xenopus, this gene had an upregulated expression in ZZ gonads. Alox12b gene was not expressed in the developing mouse gonads (Jameson et al. 2012) but was upregulated in Xenopus ZW gonads. A subpopulation of fzd receptors was expressed in the developing mouse gonads. In Xenopus, fzd4 and fzd10 had an upregulated expression in developing ZZ gonads. The calpain 8 (Capn8) was not expressed in developing mouse gonads (Jameson et al. 2012) but was upregulated in Xenopus ZW gonads. The serpins were not expressed in developing mouse gonad (Jameson et al. 2012), but they were expressed in Xenopus developing gonads. In developing mouse gonads, several cathepsins (Cts) were expressed; however, only cathepsin H (ctsh) was upregulated in XY gonads (Jameson et al. 2012), and this gene was also upregulated in ZZ Xenopus gonads. Among forkhead box factors, only Foxo1 was expressed in XY developing mouse gonads (Jameson et al. 2012) and in ZZ Xenopus gonads. Similarly, Lhx1 was expressed in XY developing mouse gonads (Jameson et al. 2012) and ZZ Xenopus gonads. Considering proteins of extracellular matrix, only collagen 9 and metalloproteinase Mmp2 were expressed in a similar manner in XY developing mouse gonads (Jameson et al. 2012) and ZZ Xenopus gonads.

Analysis of transcriptome of developing chicken gonads showed that calpain 5 (Capn5), Gpr56, and Fgfr3 were upregulated in ZW (female) gonads, which suggested that they may be involved in sexual differentiation (Ayers et al. 2015). Calpain 5 was expressed in developing Xenopus gonads, but not in a sex dimorphic manner. We showed the upregulation of calpain 8 in ZW (females) Xenopus gonads, which suggests a role of this group of proteases in sexual differentiation of vertebrate gonads. However, calpain 5 or 8 was not expressed in developing mouse gonads (Jameson et al. 2012). Gpr56 was upregulated in XY mouse and ZW chicken gonads (Ayers et al. 2015; Jameson et al. 2012), but it was not expressed in Xenopus developing gonads. Fgfr3 showed sexual dimorphism of expression in developing chicken gonads (upregulated in ZW) (Ayers et al. 2015) and was also expressed, equally in both sexes, in mouse (Jameson et al. 2012) and Xenopus gonads.

Analysis of transcriptome of a red-eared slider (T. scripta) developing gonads showed that Vwa2, Fdxr, Nov, Kdm6b, Rbm20, and Pcsk6 were upregulated in the male-producing temperature, while Fank1, Avil, Twist1, and Hspb6 were upregulated in the female-producing temperature (Czerwinski et al. 2016). Fdxr2 and Hspb6 were also upregulated in ZW (male) developing gonads of Xenopus, but the sexual dimorphism in the level of expression was not statistically significant. Twist1 gene was slightly upregulated in ZZ gonads of Xenopus, but the sexual dimorphism in the level of expression was also not significant. We detected the expression of Nov and Pcsk6 in Xenopus gonads but these genes did not show a sexual dimorphism of expression. Among Kdms genes, we detected only the expression of kdm6a but it did not show sexual dimorphism. We did not detect the expression of Vwa2, Rbm20, Frank1, or Avil in developing Xenopus gonads.

In American alligator, the expression of Wnt11 was shown at male-producing temperature, which induces the development of the testes (Yatsu et al. 2016). We detected the expression of this gene in ZW (female) developing gonads in Xenopus. Analysis of transcriptome of zebrafish developing gonads showed that the estrogen receptor 2 (esr2) was upregulated in developing testes (Sreenivasan et al. 2008). The ZZ developing Xenopus gonads also upregulated the expression of this gene.

This comparison indicates that there is a profound difference in the pattern of gene expression and sexual dimorphism of gene expression between Xenopus and other vertebrates. Only few genes indicated above show a similar pattern of expression between Xenopus and other vertebrates. This shows how complex and fast-evolving is a molecular regulation of gonad development.

Conclusion

In this study, we revealed genes representing many functional groups, which showed sexual dimorphism of expression in developing Xenopus gonads. Some of these genes are probably involved in sex determination and sexual differentiation of the gonads. We also detected a sexual dimorphism of expression of many uncharacterized and unnamed genes. These genes should be characterized and studied further to discover if they are involved in sex determination and sexual differentiation. Comparative analysis of genes expressed in developing gonads of different classes of vertebrates showed striking inter-specific differences. Only few genes showed similarities of expression pattern between the species. This indicates how little we know and how complex, diversified, and evolutionary malleable are molecular mechanisms driving gonad development in vertebrates.

Material and methods

Animals

Tadpoles of the African clawed frog (Xenopus laevis Daudin, 1802) were raised in 10-L aquaria (30 tadpoles per 10 L) at 22 °C, fed daily with powder food Sera Micron (Sera), and staged according to Nieuwkoop and Faber (1956). The tadpoles at four stages (NF50, NF53, NF56, and NF62) were anesthetized with 0.1% MS222 solution, and the gonads were manually dissected under the dissecting microscope. All individuals used in the experiments were handled according to Polish legal regulations concerning the scientific procedures on animals (Dz. U. nr 33, poz. 289, 2005) and with the permission from the First Local Commission for Ethics in Experiments on Animals.

Sex determination by PCR

The genetic sex of each tadpole was determined using PCR detection of female-specific dm-w gene. DNA was isolated from tadpole tails using NucleoSpin Tissue Kit (Macherey-Nagel, 740952.240C). The dm-w gene (W-linked female-specific marker) and dmrt1 gene (positive control) were used to determine ZZ or ZW status of tested animals. PCR was performed as previously described (Yoshimoto et al. 2008). Following pairs of primers were used: for dm-w, 5′-CCACACCCAGCTCATGTAAAG-3′ and 5′-GGGCAGAGTCACATATACTG-3′, and for dmrt1, 5′-AACAGGAGCCCAATTCTGAG-3′ and 5′-AACTGCTTGACCTCTAATGC-3′.

Histological analysis

Bouin’s solution-fixed and paraffin-embedded samples were sectioned at 4 μm. Sections were deparafinated, rehydrated, and stained with hematoxylin and picroaniline according to Debreuill’s procedure (Piprek et al. 2012). Sections were viewed under the Nikon Eclipse E600 microscope.

RNA isolation

Total RNA was isolated using Trizol and purified with Direct-zol RNA kit according to the manufacturer’s protocol (Zymo Research, R2061). The total RNA was quantified using NanoDrop 2000, and RIN (RNA Integrity Number) was assessed with Bioanalyzer 2100. All samples used in the study had RIN above 8. In order to obtain a sufficient amount of RNA, the samples from 10 individuals were pooled in each experiment as previously described (Piprek et al. 2018). Total RNA in RNase-free water was frozen at − 80 °C until further use.

Microarray analysis

Microarray analysis was performed as previously described (Piprek et al. 2018). Total RNA was labeled with fluorescent dyes using Agilent One-Color Quick Amp Labeling Protocol. RNA isolated from ZW gonads were labeled with Cy3, and RNA from ZZ gonads with Cy5. Fluorescently labeled RNA samples were mixed with Agilent Hi-RPM Hybridization Buffer, and hybridized at 65 °C for 17 h in HybArray12 hybridization station (Perkin Elmer). RNA from ZW and ZZ were mixed together and hybridized to the same chip. The RNA isolated from the gonads in different stages of development was labeled with the same fluorochrome (either Cy3 or Cy5) and hybridized individually to the separate chips. Samples were washed in Gene Expression Wash Buffer 1 (6X SSPE, 0.005% N-lauroylsarcosine; at RT) and Gene Expression Wash Buffer (0.06X SSPE, 0.005% N-lauroylsarcosine; at RT) for 1 min each and immersed in a solution of acetonitrile. Air-dried slides (custom-commercial Agilent-070330 X. laevis Microarray slides) were scanned in the Agilent Technologies G2505C Microarray Scanner at a 5-μm resolution. The microarray experiment was repeated three times.

Data processing

Data processing was performed as previously described (Piprek et al. 2018). TIF files obtained in microarray scanner were processed using Agilent Feature Extraction software version 10.5.1.1. Control and non-uniform features were removed; remaining values for each unique probe sequence were averaged. Log base 2 intensities were median centered between arrays. Differential gene expression was filtered using a statistical significance threshold (FDR < 0.05) and a fold change threshold (2-fold). The data were published in Gene Expression Omnibus (accession number GSE105103). Functional analysis and gene ontology were carried out using DAVID 6.8 (https://david.ncifcrf.gov/tools.jsp) and IPA (Ingenuity Pathway Analysis, Qiagen). First, we compared the level of gene expression between gonads in different stages of development within each sex. The gene expression level at each stage of gonad development was compared to the gene expression level at the previous developmental stage, i.e., the stage NF53 was compared to the stage NF50, the stage NF56 was compared to the stage NF53, and the stage NF62 was compared to the stage NF56. In each comparison, the level of gene expression in the younger stage of gonad development was arbitrarily designated as the reference level of expression. The results of these analyses gave us an overview of the pattern of gene expression in consecutive stages of gonad development. Subsequently, we compared the level of gene expression between genetic female (ZW) versus male (ZZ) gonads at each studied developmental stage.