- Research article
- Open Access
Transcriptome display during tilapia sex determination and differentiation as revealed by RNA-Seq analysis
BMC Genomics volume 19, Article number: 363 (2018)
The factors determining sex in teleosts are diverse. Great efforts have been made to characterize the underlying genetic network in various species. However, only seven master sex-determining genes have been identified in teleosts. While the function of a few genes involved in sex determination and differentiation has been studied, we are far from fully understanding how genes interact to coordinate in this process.
To enable systematic insights into fish sexual differentiation, we generated a dynamic co-expression network from tilapia gonadal transcriptomes at 5, 20, 30, 40, 90, and 180 dah (days after hatching), plus 45 and 90 dat (days after treatment) and linked gene expression profiles to both development and sexual differentiation. Transcriptomic profiles of female and male gonads at 5 and 20 dah exhibited high similarities except for a small number of genes that were involved in sex determination, while drastic changes were observed from 90 to 180 dah, with a group of differently expressed genes which were involved in gonadal differentiation and gametogenesis. Weighted gene correlation network analysis identified changes in the expression of Borealin, Gtsf1, tesk1, Zar1, Cdn15, and Rpl that were correlated with the expression of genes previously known to be involved in sex differentiation, such as Foxl2, Cyp19a1a, Gsdf, Dmrt1, and Amh.
Global gonadal gene expression kinetics during sex determination and differentiation have been extensively profiled in tilapia. These findings provide insights into the genetic framework underlying sex determination and sexual differentiation, and expand our current understanding of developmental pathways during teleost sex determination.
Sex-determination mechanisms in mammals and birds are highly conserved. The Sry gene, which arose by gene duplication 170 million years ago (Mya), determines male sex in most mammals . Analysis of human and mouse mutants has identified many additional genes involved in sexual differentiation, including Dmrt1, Foxl2, Rspo1, Sox9, and Wnt4 [2, 3]. Some of these same genes have been implicated as sex determiners in other tetrapods. Most birds share a ZW system of sex determination that arose 130 Mya, which is mediated at least in part by expression of Dmrt1 . A duplicated and truncated Dmrt1 gene on the female-specific W chromosome, Dmw, is required for female sex determination in Xenopus laevis [5, 6]. The fact that many of the same genes are involved in sex differentiation in these distantly related species suggests that the gene regulatory network underlying sex determination and sexual differentiation might be fundamentally similar among all tetrapods.
In contrast, teleost fishes are remarkably diverse in their sex chromosome systems [7,8,9]. Seven master sex-determining genes have been identified in teleosts, namely dmy/dmrt1bY in Oryzias latipes and Oryzias curvinotus [10, 11], gsdfY in Oryzias luzonensis , sox3 in Oryzias dancena , amhy in Odontesthes hatcheri  and Oreochromis niloticus , amhrII in Takifugu rubripes , gdf6Y in Nothobranchius furzeri  and sdY in Oncorhynchus mykiss and several other salmonids [18, 19]. Moreover, sex chromosomes and sex determining genes of the same species were found to be diverse. For example, a tandem duplicate of Amh with a missense SNP on LG23 may contribute to male sex determination in some strains [15, 20], while main sex determining region of other strains has been found on LG1 . The diversity of sex determining genes suggest there are at least some differences between mammals and teleosts in the structure of the gene network controlling sex determination. Indeed, Gsdf is found only in teleosts, and the teleost genomes lack Fgf9, suggesting that the specific function, regulation, and interconnection of the sex network has changed during the evolution of teleosts. Nevertheless, most of the genes Foxl2, Dmrt1, and pathways such as Rspo1/Wnt/TGF-β which are important in mammalian sex determination, still play key roles in teleost sex determination .
Research in Drosophila and Caenorhabditis identified linear pathways of gene interaction leading to sex determination. These pathways are thought to have evolved largely by addition of upstream master controllers, although the pathways may also be vulnerable to takeover at intermediate steps [23, 24]. In contrast, vertebrate pathways have been characterized as balanced networks, in which changes in any of a number of genes can shift the fate of the undifferentiated gonad [25, 26]. The developmental pathways following the initiation of sex determination in fishes may consist of a conserved cascade of genes regulated by sex steroid hormones . Data from natural and induced sex reversal after primary sex determination in various fishes also suggest that sex is not a stepwise, hierarchical trait, but instead a consensus of interconnected gene networks that also incorporate environmental factors [28,29,30].
There are two general approaches to characterize the gene network underlying sex determination and differentiation in a particular evolutionary lineage. The first is to genetically map and identify the top-level sex determiner in each of several species. This is essentially a ‘natural mutant screen’ in which evolution has identified various genes affecting the underlying regulatory network. This approach has been particularly successful in medaka  and tilapia . The second approach is systems biology. Sex-related genes frequently exhibit sexually dimorphic patterns of expression in the developing gonad both before and after overt differentiation of the testis or ovary. If we collect transcriptome data from male and female gonads at various developmental stages, we can use statistical methods to infer the regulatory network from the correlations of gene expression across samples. These two approaches are not mutually exclusive. Gonadal transcriptomic data with different top-level sex determiners might be combined to understand their shared gene regulatory network. Together, these approaches have identified several genes in the pathway of sex determination and sexual differentiation, but our understanding of it is far from complete. Gene co-expression analysis, uses the correlation (or related measures) of gene expression profiles across multiple samples to identify common patterns of regulation . The most widely used method is Weighted Gene Co-expression Network Analysis (WGCNA) . This method discovers groups of genes with similar patterns of gene expression, called modules. Regulatory networks and candidate genes associated with gonad differentiation have been identified using WGCNA in turbot and catfish [33, 34].
Tilapia is the world’s second most farmed fish, and an important source of animal protein around the world. Male tilapias grow faster than females, so the tilapia industry prefers to grow all-male populations. A better understanding of the molecular mechanisms of sex determination and differentiation will facilitate the production of mono-sex progeny for commercial production. The objectives of this study were to explore the molecular mechanism of tilapia sex determination and differentiation and to identify new genes involved in this process. Our approach was to collect a set of gonadal transcriptomes from genetic females and males, and sex-reversed individuals, through the course of gonad development from 5 to 180 days after hatch (dah). WGCNA was used to identify gene co-expression modules that are differentiated between ovary and testis during gonad development. Network neighbors of some canonical genes involved in vertebrate sex differentiation, including Cyp19a1a, Foxl2, Dmrt1, Gsdf, 3beta-Hsd, Cyp11b2, and Amh were identified. Our ultimate goal was to identify putative regulatory links that might be tested by CRISPR modifications.
Fish materials and ethics statement
The founder strain of the Nile tilapia, firstly introduced from Egypt in Africa, was obtained from Prof. Nagahama (Laboratory of Reproductive Biology, National Institute for Basic Biology, Okazaki, Japan). Nile tilapia were maintained and reared in re-circulating aerated freshwater tanks at 26 °C prior to use. The monosex genetic female (XX) and male (XY) tilapia were obtained as described previously . All fish experiments were conducted in accordance with the regulations of the Guide for Care and Use of Laboratory Animals and were approved by the Committee of Laboratory Animal Experimentation at Southwest University.
Gonadal histology, RNA extraction and Illumina library preparation
All fishes were euthanized using an overdose of MS222 and stored in 70% EtOH. Special care was taken to minimize suffering of fish. Gonads were dissected at 5, 20, 30, 40, 90, and 180 dah, fixed in Bouin’s solution for 24 h at room temperature, dehydrated, and embedded in paraffin. All tissue blocks were sectioned at 5 μm and stained with hematoxylin and eosin.
A total of 300, 200, 150, and 120 gonads for each sex were pooled at 5, 20, 30, and 40 dah, respectively. Total RNA was extracted from gonads of XX and XY fishes using the Trizol Reagent (Invitrogen, Carlsbad, CA) according to the manufacturer’s instruction, and purified with DNaseI (RNase-free 5 U/uL) to eliminate genomic DNA contamination. All the tilapia gonadal samples used in sequencing were again genetically sexed using a marker tightly linked to the AMH polymorphism on LG23 . Both agarose gel electrophoresis and a Nanodrop spectrophotometer were used to assess the integrity and concentration of RNA. Poly-T Oligo-attached magnetic beads were used to isolate poly(A) mRNA from total RNA. The first and second strand cDNA were synthesized form the fragmented poly(A) mRNA according to the Illumina sample preparation guide. Sequencing was performed on the HiSeq system. The reads have been deposited in NCBI’s Sequence Read Archive (SRA) database with accession number SUB3254083. Gonadal transcriptomes of normal development at 90 and 180 dah were downloaded from SRA database . Gonadal transcriptomes of 45 dat (days after treatment) and 90 dat were from Fadrozole-treated XX fishes started at 90 dah . These 90 dah-old XX fish were fed with a diet sprayed with 95% ethanol containing Fadrozole (Sigma) at a concentration of 200 μg/g diet for 45 (45 dat) and 90 days (90 dat).
After trimming adapters and low quality sequence with Trimmomatic , the clean reads were aligned to a tilapia reference genome assembly  using Tophat , allowing up to 2 base mismatches per read. NCBI RefSeq annotations were used to guide the Cufflinks assembly. Cuffdiff and Cuffnorm  were used to determine gene expression. A threshold of FPKM > 1 in at least one sample was used to filter the extremely low expression gene. Clustering of samples by global gene expression data (transformed to log2FPKM) was constructed using the heatmap function of Gplot in R. A false discovery rate (FDR) < 0.05 and two-fold difference in expression between XX and XY gonads at the same stage were used as thresholds to identify differentially expressed genes (DEGs).
Gene co-expression network analysis
To further understand the relationships between genes, a network analysis based on gene-to-gene correlations was performed using WGCNA, an R package . The automatic one-step network construction and module detection method was used, which include a signed topological overlap matrix (TOM), a softPower of 18, a minimal module size of 30, and a branch merge cut height of 0.19. The module eigengene (the first principal component of a given module) was calculated and used to test the association of modules with “Genotype”, “Phenotype”, “Developmental stages”, and “key genes (Cyp19a1a, Foxl2, Dmrt1, Gsdf, 3beta-Hsd, Cyp11b2, and Amh)” of 16 samples. A value of 0 was assigned to sex reversal samples, a value of 1 to female samples and a value of − 1 to male samples. Gene significance (GS, the correlation between gene expression and traits), total network connectivity, and module membership (also known as eigengene-based intramodular connectivity), were calculated for all modules. The most relevant connections with key genes (Cyp19a1a, Foxl2, Dmrt1, Gsdf, 3beta-Hsd, Cyp11b2, and Amh), using a high threshold power and differential expression in at least 4 stages, were kept as the candidate genes.
Data validation by in situ hybridization (ISH)
To validate which population of cells in the developing gonads express particular genes, we performed ISH on ovaries and testes from adult tilapia (180 dah) as described previously . Gonads were dissected and fixed in 4% paraformaldehyde in 0.1 M phosphate buffer (pH 7.4, 4% PFA) at 4 °C overnight. After fixation, the tissues were embedded in paraffin and cross sections were cut at 5 μm. Digoxigenin (DIG)-labeled sense and antisense RNA probes were transcribed in vitro from linearized pGEM-T easy plasmids containing-Borealin, Gtsf1, Zar1, C15orf65 or Rbp2 cDNA using a RNA labeling kit (Roche).
Results and discussion
Timeline of gonad development
Histological staining can be used to observe morphological differentiation during the development of tilapia gonads . In this study, we focused on samples beginning at the earliest stage gonads can be dissected (5 dah) through to sexual maturation (180 dah). Histological observations were made of the XX and XY gonads at 5, 20, 30, 40, 90, and 180 dah. Sex determination occurs during the earliest stages of development, when the fate of bi-potential gonadal primordium is directed toward becoming a testis or an ovary. Typically, there is no overt morphological differentiation during the period of sex determination. Between 5 and 20 dah, no morphological differences could be observed in the germ cells of XX and XY gonads (Fig. 1a, b, c, and d). Only a handful of genes show transcriptional differences between male and female gonads at 5 dah. During the period from 5 to 20 dah, the primordial gonad loses its intersexuality and the phenotypic sex is fixed as in other gonochoristic fish . Sexual differentiation begins between 30 and 40 dah, as demonstrated by the initiation of meiosis with appearance of the oocytes. At this stage, the germ cells in the testis have not initiated meiosis (Fig. 1e, f, g, and h). At 90 dah, spermatogenesis began in testis, while an ovarian cavity was observed in the females (Fig. 1i and j). At 180 dah, females develop a mature ovary; in its swollen anterior end, large previtellogenic oocytes are present in a lamellar structure. At the same time, males develop a mature testis characterized by the appearance of germ cell meiosis and numerous spermatozoa in the efferent duct (Fig. 1k and l).
Gonadal transcriptome sequencing of tilapia gonad
To generate a comprehensive profile of the tilapia gonadal transcriptome during sex determination and differentiation, we sampled XX and XY gonads at 5, 20, 30, 40, 90, and 180 dah. Transcriptomes were sequenced to a depth of ~ 90 million reads after trimming, and 78% of these were mapped to the tilapia genome assembly. Of these, 90% reads were mapped to annotated regions and 96.4% of the mapped reads were placed uniquely, indicating the good quality of reads. Overall, 28,586 transcripts were expressed in gonads during at least one of the developmental stages. A threshold FPKM > 1 in at least one sample was used to define genes as robustly expressed [43, 44]. In total, 19,188 genes were retained for subsequent analysis. At later stages of development (40 dah, 90 dah, and 180 dah), thousands of genes were differentially expressed between XX and XY gonads. These genes are enriched for functional annotations related to transcription factor binding and nuclear hormone receptor binding.
Global gene expression diverges at 30 dah
The gene expression patterns of gonads at 5 (with two replicates), 20, 30, 40, 90, and 180 dah, as well as that of 45 and 90 dat, were hierarchically clustered in a heatmap using the Pearson correlation coefficient as distance measure (Fig. 2). Two different groups were clearly identified. One consisted of fishes at earlier stages (Branch I: 5, 20, 30, and 40 dah), and the other consisted fishes at later developmental stages (Branch II: 90, 180 dah). In Branch I, XX and XY gonads are clustered together at both 5 and 20 dah. The first molecular signs of sex differentiation were observed at 5 dah, characterized by increased expression of Dmrt1 in the testis, and increased expression of Cyp19a1a in the ovary, as reported in previous studies [35, 42]. These genes are an indication of some of the earliest steps in the female and male cascade during gonadal sex determination in tilapia. However, most other genes showed no obvious differences between males and females at 20 dah. Global gene expression of male and female gonads started to differentiate at 30 dah, the period when morphological differentiation of the gonads is first observed. Gonadal samples at later stages consistently clustered into two groups corresponding to the XX and XY genotypes. Gene expression of 45 dat sample grouped with that of XX gonads at 90 and 180 dah, while gene expression of 90 dat sample clustered with that of XY gonads at 90 and 180 dah.
Network analysis and association of sex candidate genes with modules
Regulatory networks can be organized into modules of highly interconnected genes, which are implicated in the same biological processes or associated with particular cellular states. Gene co-expression analysis is based on the similarity of gene expression profiles across multiple samples or conditions. In general, adding more datasets obtained from different samples increases the predictive power of this method . In the present study, we performed a network analysis based on 16 samples using WGCNA, which partitions a gene set into modules defined as branches of the co-expressed gene cluster tree. Our analysis resulted in 23 modules (Fig. 3) with sizes ranging from 37 to 7720 genes. For purposes of discussion, each of these modules is assigned a color. Unassigned genes might belong to smaller expression pathways because of the chosen minimal module size of 30 genes. Next, we used ANOVA to test for associations between the WGCNA modules and genotype, phenotype, developmental stages (dah), and the expression of key genes in vertebrate sexual development. Investigation of the relationships between the module eigengenes and traits uncovered correlation coefficients that varied widely from -0.57 to 0.54 with genotype and from − 0.86 to 0.81 with developmental stage. The black and brown modules were significantly associated with the developmental stage, and grouped more than 72% of the genes included in the analysis. This result reflected the dramatic changes in gene expression associated with gonadal morphogenesis. Most of the genes already known to be involved in sex differentiation were assigned to these two large modules, including Cyp19a1a, AmhrII, Sox9b, Fst, Sox3, Gata5, and Lim1. Several genes in the hedgehog (HH) signaling pathway, including Shh, Ihh, and Ptch2, were also placed in brown and black modules.
Genes within a module are frequently known to be functionally related. For instance, Piwil1 and Sox30 from the midnight module are both involved in mRNA binding and regulation of translation [46, 47]. Similarly, another study has provided evidence to suggest that Piwil1 was expressed with a proposed regulatory role in striped bass . The tan module included genes involved in epigenetic regulation throughout gonad development, including histones and methyltransferase. Some homologs, such as Ara and Arb, Sox9a and Sox9b are in the same module. However, other homologous genes, Cyp19a1a and Cyp19a1b, Sox8a and Sox8b were placed in different modules, suggesting a functional difference of these duplicate genes. The modules steelblue and darkturquoise were associated with canonical genes reported to be important in testis development, including Amh, Cyp11b2, Dmrt1, Gsdf, and 3beta-Hsd. Other genes in these modules, such as Tex2, Er2, Fgfr1, Ar, and Cyp2k1, were known to be key regulators of testis differentiation. Genes in the pink module were also preferentially expressed in males, especially at the later stages (90 and 180 dah), suggesting their crucial roles in spermatogenesis. For example, null Klhl10 mice causes haploinsufficiency with meiotic arrest, absence of mature spermatozoa in semen and male infertility . Spata6 was essential for the connecting piece formation and tight head–tail conjunction during spermiogenesis in human . The pink module also includes Tdrd6 and Hira, which are involved in transcription regulation via modification of histones [51, 52]. Most of the genes in the smallest module sienna3, including gene Hsp11 and Svep1, were overexpressed in the ovaries and significantly associated with Foxl2.
Apoptosis associated with the masculinization of the female gonad
The effects of Fadrozole (aromatase inhibitor) on gene expression during sex differentiation have been investigated in several teleost species, including zebrafish , rainbowfish , tilapia , and pejerrey . Regardless of the underlying differences in sex determining mechanism, a shared characteristic of these species is that apoptosis of somatic cells in the gonads plays a key role during testicular differentiation. Degeneration of oocytes, together with the suppression of female pathway genes and activation of male pathway genes, was observed during Fadrozole-induced sex reversal of tilapia . In the present study, module white had the most significant correlation with the trait of phenotype (r = − 0.5). Genes in this module were enriched in functional annotations relevant to PPAR signaling pathwayand oxidoreductase activity. The expression of apoptosis-involved genes precedes the development of testicular structures and degeneration of oocytes in fadrozole-treated fish. These apoptosis-related processes likely play a role in the transdifferentiation of the gonads during sex reversal, as already known from other fishes [56, 57].
Candidate genes involved in tilapia sex differentiation
The top-level genes controlling sex determination vary among vertebrate species, but the downstream genes that control sexual differentiation appear to be relatively conserved . To identify other candidate genes that might be involved in tilapia sex determination and differentiation, we identified sex-biased genes that also showed strong correlations to key genes (Cyp19a1a, Foxl2, Dmrt1, Gsdf, Amh, and 3beta-Hsd) already known to be involved in these processes in other vertebrates (Additional file 1). A hierarchically-clustered heatmap of known sexual development genes and candidate genes showing sexually dimorphic expression in at least 4 developmental stages was constructed additionally to display the expression of these genes (Fig. 4). Interestingly, this subset of genes produced a different clustering result compared with that observed in the overall gene cluster. The gene expression patterns of XX gonads at 20, 30, 40, 90, and 180 dah were hierarchically clustered, while gene expression of XY gonads at these stages grouped together. Namely, these genes displayed obviously different sex-specific gene expression profiles as early as 20 dah. However, similar clustering pattern of Fig. 2 was found while we constructed heatmap using all the genes in Additional file 1, indicating most genes in this list were involved in tilapia sexual differentiation at later stages. Until now, no example of a sex determining gene is still highly expressed in adult fish. For example, the expression of the sex determiner Amhy in Nile tilapia was detected at 5 dah, increased at 20 dah and peaked at approximately 40 dah, and decreased at 90 and 180 dah . Thus, further study of sex determiners expression covering the development from embryo to adults is needed to explore the gene regulatory network during early sex determination.
Knockout of Foxl2 and Cyp19a1a resulted in 100% complete female to male sex reversal , while disruption of AmhRII and Gsdf resulted in male to female sex reversal with normal ovaries in tilapia [15, 59]. Besides these previously reported sex-determining genes, some new candidates genes identified in present study possibly play an important role in tilapia sex determination and sexual differentiation. For example, mutations of Gtsf1 in Drosophila caused de-repression of transposons and loss of ovary follicle layers, resulting in female infertility . Zar1 displayed ovary-specific expression in mice , and it was up-regulated at XX tilapia gonad. Notably, proteomic analysis of striped bass ovary indicated that Zar1 was also predominantly detected in early secondary and vitellogenic growth ovary (equivalent to phase I and phase II oocytes in tilapia) and protein expression dramatically decreased during post-vitellogenesis . A recent study indicated that mutation of Zar1 in zebrafish resulted in early oocyte apoptosis and fully penetrant male development, while the EE2-treated Zar1 homozygous mutants were recovered as females .
Other genes were not sexually dimorphically expressed at 5 and 20 dah, but highly expressed at later stages. These genes, including Pacap, Hsp70, Lgr6, Mphosph6, Svep1, and Tbxas1, might act late in gonad differentiation. Pacap, also known as Adcyap1, is involved in gonadotropin synthesis and release, either alone or in cooperation with GnRH [63, 64]. The coexpression of Pacap, with Cyp11b2, Tbxas1, Amh, Abcb11, and Lgr6 in tilapia testis, indicated its critical role in tilapia endocrine systems. Hsp70 is associated with gonad development and egg quality in different species [64,65,66,67,68]. In the tilapia gonad, Hsp70 was highly expressed in testis during gonad differentiation, and co-expressed with Sox9a, Gata4, and Lhx6. Consistently, a previous study indicated that HSP70 was an interacting partner for SOX9 . Thus, Hsp70 is a good candidate for future studies at investigating sex-related regulatory network in tilapia and possibly other fish. The higher expression of Lgr6 in tilapia testis, together with Lgr4 and Lgr5, might encode orphan 7-transmembrane receptors and mediate Wnt/β-catenin and Wnt/PCP signaling in tilapia gonad differentiation, as indicated in mammals [69,70,71]. Mphosph6 play a critical role in M-phase characteristics of growing oocytes in mouse . In this study, Mphosph6 was highly expressed in tilapia ovaries from 90 to 180 dah, indicating its its role in tilapia oocyte maturation. A previous study showed that both estradiol and TNFα regulated Svep1 expression , but its function in sexual differentiation has not been reported yet. Tbxas1, also called Cyp5, was highly expressed in mouse testis and essential for testis development .
Finally, 17 uncharacterized transcripts were found to be differentially expressed in XX and XY gonads at various time points. Interestingly, 9 of these transcripts were identified as ncRNAs (LOC100696883, LOC102076355, LOC102077364, LOC102077676, LOC102078352, LOC102079282, LOC102080335, LOC102081303, LOC102082652). The presence of these ncRNAs may have significant consequences on tilapia sex determination as ncRNA (MHM) regulation of DMRT1 plays a pivotal role in Gallus gallus sex determination . However, their functions and roles during tilapia sexual development need to be further investigated. This again confirmed that sex determination in vertebrates was far from being a linear pathway built from the bottom up; instead different pathways and multilayered feedback loops worked together in a non-hierarchical network to produce a male or female phenotype . Both in vitro techniques (reporter assays) and in vivo (gene knockout) will be useful to evaluate such interactions and provide greater insight into the functional relationships between genes in this putative network.
Validation of candidate genes by ISH
The gonad is a mixture of different cell types. We used ISH to determine the cell types in which several of these new candidate genes, including Borealin, Gtsf1 and ZarI, were expressed (Fig. 5). By ISH, Borealin was highly expressed in the phase I and II oocytes of the ovary. Borealin was moderately expressed in the phase III oocytes, but was not expressed in the phase IV oocytes. No Borealin expression was found in the testis. Gtsf1 was weakly expressed in the oogonia and phase I oocytes, highly expressed in the phase II oocytes, and moderately expressed in the phase III oocytes of the ovary, while no Gtsf1 expression were detected in the testis. ZarI, located on LG23, was highly expressed in the phase I and II oocytes, and weakly expressed in the phase III oocytes in the ovary, and moderately expressed in the spermatocytes in the testis. Cdn15 and Rpl were expressed in XX and XY gonads at later stages, indicating their possible roles in late gonad development (Additional file 2: Figure S1). Cdn15 was highly expressed in primary spermatocytes and secondary spermatocytes, as well as in the oocytes. Rpl was expressed in spermatogonia and efferent duct, and in the phase II oocytes of the ovary. In contrast, no signal was detected in either testis or ovary using the Borealin, Gtsf1, ZarI, Cdn15 and Rpl sense probes. The fact that these novel genes recapitulate the expression patterns of the key sex genes used to identify relevant gene expression modules provides evidence for the biological relevance of these gene modules.
In this study, we chose gonadal transcriptome data of XX and XY fish at 20, 30, 40, 90 and 180 dah, as well as 45 and 90 dat, to investigate the specific events that trigger sex determination and differentiation. We found that transcriptomic profiles of female and male gonads at 5 and 20 dah exhibited high similarities. Gene co-expression analysis constructed interactive networks of highly correlated genes that modulated sexual differentiation in XX and XY gonads. The global gene expression profiles during gonad development will provide early insights into the potential processes and pathways involved in sexual differentiation. We propose the present candidate genes as a starting point for future mathematical modeling and integration of more genes regarding the regulatory network of sex determination and sexual differentiation. The network can be upgraded later in several aspects, for example, incorporating additional nodes and interactions (i.e. details of TFBS), as well as modeling different cell lineages of the gonad such as the Leydig or theca cells. A functional analysis of the identified candidate genes is required to help elucidate their potential significance in gonadal sex determination and differentiation process.
Marshall Graves JA. The rise and fall of SRY. Trends Genet. 2002;18(5):259–64.
Ottolenghi C, Pelosi E, Tran J, Colombino M, Douglass E, Nedorezov T, Cao A, Forabosco A, Schlessinger D. Loss of Wnt4 and Foxl2 leads to female-to-male sex reversal extending to germ cells. Hum Mol Genet. 2007;16(23):2795–804.
Lavery R, Chassot AA, Pauper E, Gregoire EP, Klopfenstein M, de Rooij DG, Mark M, Schedl A, Ghyselinck NB, Chaboissier MC. Testicular Differentiation Occurs in Absence of R-spondin1 and Sox9 in Mouse Sex Reversals. Plos Genet. 2012;8(12):e1003170.
Smith CA, Roeszler KN, Ohnesorg T, Cummins DM, Farlie PG, Doran TJ, Sinclair AH. The avian Z-linked gene DMRT1 is required for male sex determination in the chicken. Nature. 2009;461(7261):267–71.
Yoshimoto S, Okada E, Umemoto H, Tamura K, Uno Y, Nishida-Umehara C, Matsuda Y, Takamatsu N, Shiba T, Ito M. A W-linked DM-domain gene, DM-W, participates in primary ovary development in Xenopus laevis. P Natl Acad Sci USA. 2008;105(7):2469–74.
Mawaribuchi S, Takahashi S, Wada M, Uno Y, Matsuda Y, Kondo M, Fukui A, Takamatsu N, Taira M, Ito M. Sex chromosome differentiation and the W- and Z-specific loci in Xenopus laevis. Dev Biol. 2016;26(2):393–400.
Devlin RH, Nagahama Y. Sex determination and sex differentiation in fish: an overview of genetic, physiological, and environmental influences. Aquaculture. 2002;208(3–4):191–364.
Barske LA, Capel B. Blurring the edges in vertebrate sex determination. Curr Opin Genet Dev. 2008;18(6):499–505.
Pan Q, Anderson J, Bertho S, Herpin A, Wilson C, Postlethwait JH, Schartl M, Guiguen Y. Vertebrate sex-determining genes play musical chairs. C R Biol. 2016;339(7-8):258–62.
Matsuda M, Nagahama Y, Shinomiya A, Sato T, Matsuda C, Kobayashi T, Morrey CE, Shibata N, Asakawa S, Shimizu N, et al. DMY is a Y-specific DM-domain gene required for male development in the medaka fish. Nature. 2002;417(6888):559–63.
Nanda I, Kondo M, Hornung U, Asakawa S, Winkler C, Shimizu A, Shan ZH, Haaf T, Shimizu N, Shima A, et al. A duplicated copy of DMRT1 in the sex-determining region of the Y chromosome of the medaka, Oryzias latipes. P Natl Acad Sci USA. 2002;99(18):11778–83.
Myosho T, Otake H, Masuyama H, Matsuda M, Kuroki Y, Fujiyama A, Naruse K, Hamaguchi S, Sakaizumi M. Tracing the emergence of a novel sex-determining gene in Medaka, Oryzias luzonensis. Genetics. 2012;191(1):163.
Takehana Y, Matsuda M, Myosho T, Suster ML, Kawakami K, Shin-I T, Kohara Y, Kuroki Y, Toyoda A, Fujiyama A et al. Co-option of Sox3 as the male-determining factor on the Y chromosome in the fish Oryzias dancena. Nat Commun. 2014;5:4157.
Hattori RS, Murai Y, Oura M, Masuda S, Majhi SK, Sakamoto T, Fernandino JI, Somoza GM, Yokota M, Strussmann CA. A Y-linked anti-Mullerian hormone duplication takes over a critical role in sex determination. P Natl Acad Sci USA. 2012;109(8):2955–9.
Li MH, Sun YL, Zhao JE, Shi HJ, Zeng S, Ye K, Jiang DN, Zhou LY, Sun LN, Tao WJ, et al. A Tandem Duplicate of Anti-Mullerian Hormone with a Missense SNP on the Y Chromosome Is Essential for Male Sex Determination in Nile Tilapia, Oreochromis niloticus. Plos Genet. 2015;11(11):e1005678.
Kamiya T, Kai W, Tasumi S, Oka A, Matsunaga T, Mizuno N, Fujita M, Suetake H, Suzuki S, Hosoya S, et al. A Trans-Species Missense SNP in Amhr2 Is Associated with Sex Determination in the Tiger Pufferfish, Takifugu rubripes (Fugu). Plos Genet. 2012;8(7):e1002798.
Reichwald K, Petzold A, Koch P, Downie BR, Hartmann N, Pietsch S, Baumgart M, Chalopin D, Felder M, Bens M, et al. Insights into sex chromosome evolution and aging from the genome of a short-lived fish. Cell. 2015;163(6):1527–38.
Yano A, Guyomard R, Nicol B, Jouanno E, Quillet E, Klopp C, Cabau C, Bouchez O, Fostier A, Guiguen Y. An immune-related gene evolved into the master sex-determining gene in rainbow trout, Oncorhynchus mykiss. Curr Biol. 2012;22(15):1423–8.
Yano A, Nicol B, Jouanno E, Quillet E, Fostier A, Guyomard R, Guiguen Y. The sexually dimorphic on the Y-chromosome gene (sdY) is a conserved male-specific Y-chromosome sequence in many salmonids. Evol Appl. 2013;6(3):486–96.
Sun YL, Jiang DN, Zeng S, Hu CJ, Ye K, Yang C, Yang SJ, Li MH, Wang DS. Screening and characterization of sex-linked DNA markers and marker-assisted selection in the Nile tilapia (Oreochromis niloticus). Aquaculture. 2014;433:19–27.
Gammerdinger WJ, Conte MA, Acquah EA, Roberts RB, Kocher TD. Structure and decay of a proto-Y region in Tilapia, Oreochromis niloticus. Bmc Genomics. 2014;15:975.
Graham P, Penn JKM, Schedl P. Masters change, slaves remain. BioEssays. 2003;25(1):1–4.
Hodgkin J. Sex determination compared in Drosophila and Caenorhabditis. Nature. 1990;344(6268):721–8.
Marin I, Baker BS. The evolutionary dynamics of sex determination. Science. 1998;281(5385):1990–4.
Capel B, Tanaka M. Forward to the special issue on sex determination. Dev Dynam. 2013;242(4):303–6.
Ross AJ, Capel B. Signaling at the crossroads of gonad development. Trends Endocrin Met. 2005;16(1):19–25.
Nakamura M. The mechanism of sex determination in vertebrates-are sex steroids the key-factor? J Exp Zool Part A. 2010;313A(7):381–98.
Paul-Prasanth B, Bhandari RK, Kobayashi T, Horiguchi R, Kobayashi Y, Nakamoto M, Shibata Y, Sakai F, Nakamura M, Nagahama Y. Estrogen oversees the maintenance of the female genetic program in terminally differentiated gonochorists. Sci Rep. 2013;3:2862.
Sun LN, Jiang XL, Xie QP, Yuan J, Huang BF, Tao WJ, Zhou LY, Nagahama Y, Wang DS. Transdifferentiation of differentiated ovary into functional testis by long-term treatment of aromatase inhibitor in Nile tilapia. Endocrinology. 2014;155(4):1476–88.
Herpin A, Schartl M. Plasticity of gene-regulatory networks controlling sex determination: of masters, slaves, usual suspects, newcomers, and usurpators. EMBO Rep. 2015;16(10):1260–74.
Eisen MB, Spellman PT, Brown PO, Botstein D. Cluster analysis and display of genome-wide expression patterns. P Natl Acad Sci USA. 1998;95(25):14863–8.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. Bmc Bioinformatics. 2008;9:559.
Robledo D, Ribas L, Cal R, Sánchez L, Piferrer F, Martínez P, Viñas A. Gene expression analysis at the onset of sex differentiation in turbot (Scophthalmus maximus). BMC Genomics. 2015;16:973.
Zeng Q, Liu S, Yao J, Zhang Y, Yuan Z, Jiang C, Chen A, Fu Q, Su B, Dunham R, Liu Z. Transcriptome display during testicular differentiation of channel catfish (Ictalurus punctatus) as revealed by RNA-Seq analysis. Biol Reprod. 2016;95:19.
Tao WJ, Yuan J, Zhou LY, Sun LN, Sun YL, Yang SJ, Li MH, Zeng S, Huang BF, Wang DS: Characterization of Gonadal Transcriptomes from Nile Tilapia (Oreochromis niloticus) Reveals Differentially Expressed Genes. Plos One 2013, 8(5).
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20.
Conte MA, Gammerdinger WJ, Bartie KL, Penman DJ, Kocher TD. A high quality assembly of the Nile Tilapia (Oreochromis niloticus) genome reveals the structure of two sex determination regions. BMC Genomics. 2017;18(1):341.
Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14(4):R36.
Trapnell C, Hendrickson DG, Sauvageau M, Goff L, Rinn JL, Pachter L. Differential analysis of gene regulation at transcript resolution with RNA-seq. Nat Biotechnol. 2013;31(1):46.
Kobayashi T, Kajiura-Kobayashi H, Nagahama Y. Differential expression of vasa homologue gene in the germ cells during oogenesis and spermatogenesis in a teleost fish, tilapia, Oreochromis niloticus. Mech Develop. 2000;99(1–2):139–42.
Nakamura M, Kobayashi T, Chang XT, Nagahama Y. Gonadal sex differentiation in teleost fish. J Exp Zool. 1998;281(5):362–72.
Ijiri S, Kaneko H, Kobayashi T, Wang DS, Sakai F, Paul-Prasanth B, Nakamura M, Nagahama Y. Sexual dimorphic expression of genes in gonads during early differentiation of a teleost fish, the Nile tilapia Oreochromis niloticus. Biol Reprod. 2008;78(2):333–41.
Ramskold D, Wang ET, Burge CB, Sandberg R. An abundance of ubiquitously expressed genes revealed by tissue transcriptome sequence data. PLoS Comput Biol. 2009;5(12):e1000598.
Shin H, Shannon CP, Fishbane N, Ruan J, Zhou M, Balshaw R, Wilson-McManus JE, Ng RT, McManus BM, Tebbutt SJ. Variation in RNA-Seq transcriptome profiles of peripheral whole blood from healthy individuals with and without globin depletion. PLoS One. 2014;9(3):e91041.
D'haeseleer P, Liang SD, Somogyi R. Genetic network inference: from co-expression clustering to reverse engineering. Bioinformatics. 2000;16(8):707–26.
Pevny LH, LovellBadge R. Sox genes find their feet. Curr Opin Genet Dev. 1997;7(3):338–44.
Girard A, Sachidanandam R, Hannon GJ, Carmell MA. A germline-specific class of small RNAs binds mammalian Piwi proteins. Nature. 2006;442(7099):199–202.
Reading BJ, Williams VN, Chapman RW, Williams TI, Sullivan CV. Dynamics of the striped bass (Morone saxatilis) ovary proteome reveal a complex network of the translasome. J Proteome Res. 2013;12(4):1691–9.
Yan W, Ma L, Burns KH, Matzuk MM. Haploinsufficiency of kelch-like protein homolog 10 causes infertility in male mice. P Natl Acad Sci USA. 2004;101(20):7793–8.
Yuan SQ, Stratton CJ, Bao JQ, Zheng HL, Bhetwal BP, Yanagimachi R, Yan W. Spata6 is required for normal assembly of the sperm connecting piece and tight head-tail conjunction. P Natl Acad Sci USA. 2015;112(5):E430–9.
Iwasaki YW, Siomi MC, Siomi H. PIWI-interacting RNA: its biogenesis and functions. Annu Rev Biochem. 2015;84:405–33.
Tagami H, Ray-Gallet D, Almouzni G, Nakatani Y. Histone H3.1 and H3.3 complexes mediate nucleosome assembly pathways dependent or independent of DNA synthesis. Cell. 2004;116(1):51–61.
Villeneuve L, Wang RL, Bencic DC, Biales AD, Martinovic D, Lazorchak JM, Toth G, Ankley GT. Altered gene expression in the brain and ovaries of zebrafish (Danio rerio) exposed to the aromatase inhibitor fadrozole: microarray analysis and hypothesis generation. Environ Toxicol Chem. 2009;28(8):1767–82.
Shanthanagouda AH, Nugegoda D, Patil JG. Effects of bisphenol a and fadrozole exposures on cyp19a1 expression in the Murray rainbowfish, Melanotaenia fluviatilis. Arch Environ Contam Toxicol. 2014;67(2):270–80.
Yamamoto Y, Hattori RS, Kitahara A, Kimura H, Yamashita M, Strussmann CA. Thermal and endocrine regulation of gonadal apoptosis during sex differentiation in pejerrey Odontesthes bonariensis. Sex Dev. 2013;7(6):316–24.
Rodriguez-Mari A, Canestro C, BreMiller RA, Nguyen-Johnson A, Asakawa K, Kawakami K, Postlethwait JH. Sex Reversal in Zebrafish fancl Mutants Is Caused by Tp53-Mediated Germ Cell Apoptosis. Plos Genet. 2010;6(7):e1001034.
He Y, Shang X, Sun JH, Zhang L, Zhao W, Tian YH, Cheng HH, Zhou RJ. Gonadal apoptosis during sex reversal of the Rice field eel: implications for an evolutionarily conserved role of the molecular chaperone heat shock protein 10. J Exp Zool Part B. 2010;314b(4):257–66.
Zhang X, Li M, Ma H, Liu X, Shi H, Li M, Wang D. Mutation of foxl2 or cyp19a1a results in female to male sex reversal in XX Nile Tilapia. Endocrinology. 2017;158(8):2634–47.
Jiang DN, Yang HH, Li MH, Shi HJ, Zhang XB. Wang DS: gsdf is a downstream gene of dmrt1 that functions in the male sex determination pathway of the Nile tilapia. Mol Reprod Dev. 2016;
Ohtani H, Iwasaki YW, Shibuya A, Siomi H, Siomi MC, Saito K. DmGTSF1 is necessary for Piwi-piRISC-mediated transcriptional transposon silencing in the Drosophila ovary. Genes Dev. 2013;27(15):1656–61.
Wu XM, Viveiros MM, Eppig JJ, Bai YC, Fitzpatrick SL, Matzuk MM. Zygote arrest 1 (Zar1) is a novel maternal-effect gene critical for the oocyte-to-embryo transition. Nat Genet. 2003;33(2):187–91.
Miao LY, Yuan Y, Cheng F, Fang JS, Zhou F, Ma WR, Jiang Y, Huang XH, Wang YC, Shan LJ, et al. Translation repression by maternal RNA binding protein Zar1 is essential for early oogenesis in zebrafish. Development. 2017;144(1):128–38.
Radleff-Schlimme A, Leonhardt S, Wuttke W, Jarry K. Evidence for PACAP to be an autocrine factor on gonadotrope cells. Ann N Y Acad Sci. 1998;865:486–91.
Winters SJ, Moore JP. PACAP, an autocrine/paracrine regulator of Gonadotrophs. Biol Reprod. 2011;84(5):844–50.
Lin LL, Janz DM. Effects of binary mixtures of xenoestrogens on gonadal development and reproduction in zebrafish. Aquat Toxicol. 2006;80(4):382–95.
He Y, Luo MJ, Yi MH, Sheng Y, Cheng YB, Zhou RJ, Cheng HH: Identification of a Testis-Enriched Heat Shock Protein and Fourteen Members of Hsp70 Family in the Swamp Eel. Plos One 2013, 8(6).
Mu W, Wen H, Li J, He F. Cloning and expression analysis of a HSP70 gene from Korean rockfish (Sebastes schlegeli). Fish Shellfish Immunol. 2013;35(4):1111–21.
Sullivan CV, Chapman RW, Reading BJ, Anderson PE. Transcriptomics of mRNA and egg quality in farmed fish: some recent developments and future directions. Gen Comp Endocrinol. 2015;221:23–30.
Harley VR. The molecular action of testis-determining factors SRY and SOX9. Novart Fdn Symp. 2002;244:57–67.
Glinka A, Dolde C, Kirsch N, Huang YL, Kazanskaya O, Ingelfinger D, Boutros M, Cruciat CM, Niehrs C. LGR4 and LGR5 are R-spondin receptors mediating Wnt/beta-catenin and Wnt/PCP signalling. EMBO Rep. 2011;12(10):1055–61.
Carmon KS, Gong X, Lin QS, Thomas A, Liu QY. R-spondins function as ligands of the orphan receptors LGR4 and LGR5 to regulate Wnt/beta-catenin signaling. P Natl Acad Sci USA. 2011;108(28):11452–7.
Wickramasinghe D, Ebert KM, Albertini DF. Meiotic competence acquisition is associated with the appearance of M-phase characteristics in growing mouse oocytes. Dev Biol. 1991;143(1):162–72.
Glait-Santar C, Benayahu D. Regulation of SVEP1 gene expression by 17beta-estradiol and TNFalpha in pre-osteoblastic and mammary adenocarcinoma cells. J Steroid Biochem Mol Biol. 2012;130(1–2):36–44.
Mould AW, Duncan R, Serewko-Auret M, Loffler KA, Biondi C, Gartside M, Kay GF, Hayward NK. Global expression profiling of sex cord stromal tumors from Men1 heterozygous mice identifies altered TGF-beta signaling, decreased Gata6 and increased Csf1r expression. Int J Cancer. 2009;124(5):1122–32.
Roeszler KN, Itman C, Sinclair AH, Smith CA. The long non-coding RNA, MHM, plays a role in chicken embryonic development, including gonadogenesis. Dev Biol. 2012;366(2):317–26.
Capel B. Vertebrate sex determination: evolutionary plasticity of a fundamental switch. Nat Rev Genet. 2017;
We are grateful to Dr. Beide Fu (Institute of Hydrobiology, Chinese Academy of Sciences, Wuhan) and Dr. Yifei Shi (University of Maryland, USA) for data analysis.
This work was supported by grants 31502147, 31630082, and 31572609 from the National Natural Science Foundation of China; grant 20130182130003 from the Specialized Research Fund for the Doctoral Program of Higher Education of China; grant cstc2013kjrc-tdjs80003 and cstc2014jcyjA80001 from the Natural Science Foundation Project of Chongqing, Chongqing Science and Technology Commission, grant XDJK2018B025, XDJK2017B007, and XDJK2016A003 from the Fundamental Research Funds for the Central Universities (Ministry of Education of China).
Availability of data and materials
The datasets generated during the current study are publicly available at NCBI – SRA (http://www.ncbi.nlm.nih.gov/sra) accession: SUB3254083.
Ethics approval and consent to participate
All fish experiments were conducted in accordance with the regulations of the Guide for Care and Use of Laboratory Animals and were approved by the Committee of Laboratory Animal Experimentation at Southwest University.
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Candidate genes involved in tilapia sex determination and differentiation. (XLSX 65 kb)
Figure S1. Cellular localization of C15orf65 and Rbp2 in tilapia testis and ovary by ISH. C15orf65 was highly expressed in primary spermatocytes and secondary spermatocytes (A), as well as in the oocytes (B). Rbp2 was expressed in spermatogonia and efferent duct (C), and in the phase II oocytes of the ovary (D). PSC, primary spermatocytes; SSC, secondary spermatocytes; OC, oocytes; I-IV, phase I to phase IV oocytes; SG, spermatogonia; ED, efferent duct; Arrowheads indicate the positive signal. (TIF 40771 kb)