Sex determination in annual fishes: Searching for the master sex-determining gene in Austrolebias charrua (Cyprinodontiformes, Rivulidae).

Evolution of sex determination and differentiation in fishes involves a broad range of sex strategies (hermaphroditism, gonochorism, unisexuality, environmental and genetic sex determination). Annual fishes inhabit temporary ponds that dry out during the dry season when adults die. The embryos exhibit an atypical developmental pattern and remain buried in the bottom mud until the next rainy season. To elucidate genomic factors involved in the sex determination in annual fish, we explored the presence of a candidate sex-specific gene related to the cascade network in Austrolebias charrua. All phylogenetic analyses showed a high posterior probability of occurrence for a clade integrated by nuclear sequences (aprox. 900 bp) from both adults (male and female), with partial cDNA fragments of A. charrua from juveniles (male) and the dsx D. melanogaster gene. The expressed fragment was detected from blastula to adulthood stages showing a sexually dimorphic expression pattern. The isolated cDNA sequence is clearly related to dsx D. melanogaster gene and might be located near the top of the sex determination cascade in this species.


Introduction
In vertebrates different master genes acting on sex determination have been identified. In most mammals, Sry (testis-determining gene) located in the Y chromosome is the transcription factor that triggers the testis determining cascade. The first non-mammalian master gene was discovered in fish (Oryzias latipes dmy). In birds, Dmrt1 (doublesex-mab-3 related transcription factor 1) located in the Z chromosome plays this role while in Xenopus laevis, DM-W a truncated copy of Dmrt1 found in the W chromosome, determines ovary fate in this species (Marín and Baker, 1998;Graves and Peichel, 2010).
Recently, four strong master sex determining candidate genes were identified in fish: amhy in Odontesthes hatchei, gsdf in Oryzias luzonensis, amhr2 in Takifugu rubripes and sdY in Oncorhynchus mykiss. Three of the four sex determining candidate genes (amhy, gsdf and amhr2) code for growth factors and one of their receptors demonstrating that novel actors, other than transcription factors, can be recruited at the top of the sex determination cascade (reviewed by Kikuchi and Hamaguchi, 2013).
Downstream genes involved in sex differentiation regulatory cascades are conserved in vertebrates and inver-tebrates (Smith et al., 1999). Among these genes, a DM gene family, Dmrt (doublesex-mab-3 related transcription factor) is expressed in association with the development of sex-specific organs in all animals studied to date (Kopp, 2012).
The first identified family member: Dmrt1, exhibits an expression pattern mainly involved in postnatal (or post hatching) male gonad development (Hodgkin, 2002;Yamaguchi et al., 2006). In some fish species, its expression pattern is restricted to testes (e.g Oryzias latipes, Kobayashi et al., 2004) while in others it is expressed in both gonads (e.g. Odontesthes bonariensis, Fernandino et al., 2008). In hermaphrodite fish species Dmrt1 expression is related to male differentiation phases (Herpin and Schartl, 2011). The interest in this gene was greater after the discovery that Dmrt1 paralogs have moved up in the regulatory hierarchy from downstream position in gonad differentiation to the top of sex determination cascade in at least three distantly related organisms (O. latipes, Xenopus laevis and chicken Kopp, 2012). In O. latipes the male sex specific Dmrt1 copy, Dmrt1b(Y) (Nanda et al., 2002) or dmy (Matsuda et al., 2002), is located in the Y chromosome with a function equivalent to mammalian Sry (male master sex determining gene).The Dmrt1b(Y) gene is expressed in male embryos before gonadal differentiation (early development: neurula stage, Nanda et al., 2002). It is involved in the specification and maintenance of Sertoli cells fate and inhibits male germ cell division at the beginning of gonadal differentiation . The autosomal paralogue Dmrt1a begins its expression between 20-30 days post hatching during testicular differentiation (Kobayashi et al., 2004). In adult testes both paralogues are expressed but Dmrt1a predominates (Hornung et al., 2007). It has been demonstrated that an insertion of a transposable element is responsible for the regulation of Dmrt1bY expression and also contributed to the establishment of this new regulatory hierarchy (Herpin et al., 2010).
Fishes represent the most basal and diverse vertebrate group enclosing next to 28,000 species (Nelson, 2006). This diversity includes different mechanisms involved in sex determination (genetic and environmental) and sex differentiation (unisexuality, hermaphroditism: synchronic or sequential and gonochorism: undifferentiated or differentiated) (Devlin and Nagahama, 2002). Sex determination systems differ in closely related fish species, even in populations of the same species (Conover and Kynard, 1981;Devlin and Nagahama, 2002). The evolutionary basis of this variability could be explained by ancient genomic duplication events that caused extra gene copies capable to acquire new functions and probably additional plasticity in the sex determination gene networks. The ability to modify sex determination control could be selected in response to environmental disturbances affecting sexual proportions (Mank et al., 2006). This fact could be critical as an adaptation to environmental shifts (e.g. temperature increments) that result in one sex proportion bias affecting population survival. In this context, the possibility to modify sex determination mechanisms could restore sex ratio balance according to the specific habitat condition in threatened species (Volff et al., 2007).
Annual fishes (Cyprinodontiformes; Aplocheiloidei) constitute a freshwater teleosts group with a short lifespan exposed to an extremely variable environment. They inhabit temporary ponds in South America and Africa that dry out during dry season leading to juvenile and adult death. The species survival depends on desiccationresistant embryos that remain hidden in the bottom mud until the next rainy season when they hatch (Myers, 1952;Wourms, 1967, Berois et al., 2012. In particular, Austrolebias charrua (Costa and Cheffe, 2001) is an annual fish species distributed from southern Brazil (Patos-Merín lagoon) to eastern Uruguay (Rocha Department). Chromosomal studies and analyses at DNA level demonstrated high genetic variations (García, 2006) and the existence of ancestral polymorphisms in certain A. charrua populations (García et al., 2009). Aspects related to reproductive diversity are essential to understand the evolution of these mechanisms as well as management of species that are considered potential environmental pollution biomonitors (Devlin and Nagahama, 2002).
Sex differentiation, determination of sexual strategy, and gametogenesis of this species were previously estab-lished (Arezo et al., 2007). There is no evidence about the mechanism of sex determination and no genetic sex markers have been identified so far. Moreover, until now there are no complete genome sequences available. Therefore, the aim of the present work is to focus on the identification of a key gene involved in A. charrua sex determination. As a first step we explored the presence of a candidate master sex-specific gene and its expression during ontogenesis.
Fragments of different length amplified in each individual were eluted with the GFX PCR DNA and Gel Band Purification kit (GE Healthcare Life Sciences) and cloned with the GeneJet 1.2 PCR cloning kit according to manufacturer's instructions (Fermentas). Recombinant plasmids were obtained by DNA minipreparations of individual clones by alkaline lysis (Sambrook et al., 1989). Sequencing reactions were performed on each template using the primers supplied in the cloning kit (pJET1.2 forward and pJET1.2 Reverse) in a Perkin-Elmer ABI Prism 377 Automated Sequencer (MACROGEN, Seoul, Korea). Nucleotide sequences were compared against the National Center for Biotechnology Information (NCBI) protein da-tabase using the BLAST program on the Basic Local Alignment Search Tool (BLAST) network service (Altschul et al., 1990).
Fragments were eluted, cloned, sequenced and compared against the NCBI database as described above. The obtained RNA sequence was analyzed online at the RNAfold WebServer (Gruber et al., 2008;Lorenz et al., 2011). Default parameters at 20°C were selected.

Phylogenetic analyses
In order to resolve the phylogenetic relationships, an alignment with A. charrua isolated genomic fragments as described above and the Dmrt gene family sequences retrieved from GenBank (Table 2) was conducted using ClustalW (Thompson et al., 1997) implemented in MEGA4 (Tamura et al., 2007). Three data set were generated to develop phylogenetic analyses: a) one including four A. charrua isolated genomic fragments (two for male and two for female) vs. Dmrt gene family sequences from different fish groups (Table 2); b) we also carried out an analysis to reveal relationships between the expressed A. charrua partial cDNA fragment and Dmrt1 mRNAs retrieved from GenBank (Table 3); c) a reduced data set including four A. charrua isolated genomic fragments, the expressed A. charrua partial cDNA fragment and the dsx D. melanogaster sequence.

366
A. charrua sex determination   We also used two model-based approaches, i.e., maximum-likelihood (ML) and Bayesian inference (BI) implemented in BEAST v.1.5.4 (Drummond and Rambaut, 2007). For the ML and Bayesian analyses, the best-fitted nucleotide substitution models for the data set were determined in Modeltest v.3.7 (Posada and Crandall, 1998) based on the Akaike Information Criterion (Akaike, 1974), which simultaneously compares multiple nested or non-nested models.
For data set a) the best fit out of the 56 models was obtained with the general time-reversible model (GTR + G, Rodriguez et al., 1990), which includes six different nucleotide substitution types and variable substitution rates among sites drawn from a gamma distribution (G). The gamma distribution shape parameter in the data set was 0.7311. For data set b) the best fitted model selected was K81uf+G model (K81+G: Kimura 3-parameter plus Gamma, Kimura 1981) with a gamma distribution shape parameter of 0.6519. In dataset c), HKY85 (Hasegawa et al. 1985) resulted in the best fitted model of molecular evolution.
For data set a) the likelihood scores estimated for the selected model GTR +G was used as the prior settings for the ML analysis (-ln L = 32940.4721). A heuristic ML search (again with 100 replicates of stepwise addition and TBR branch swapping) implemented in PAUP* 4.0b10 (Swofford, 2002) was used. The robustness of the nodes was determined after 500 bootstrapping replicates as implemented in PhyML 3.0, according to the algorithm developed by Guindon et al. (2010). In this case, the NNI (a fast nearest neighbour edge interchange search) swapping algorithm option was implemented. For data set b) the likelihood scores estimated for the selected model K81+G were used as the prior settings for the ML analysis (-ln L = -15985.54790). Finally, for the c) data set the likelihood scores estimated for the selected model HKY85 were used as the prior settings for the ML analysis (-ln L = -723.286). All trees were rooted by means of an outgroup criterion using more distantly related sequences in each data set.
Bayesian posterior probabilities of the trees were calculated using the BEAST v.1.5.4 program (Drummond and Rambaut, 2007). BEAST performs Bayesian statistical inferences of parameters, such as divergence times, using MCMC as a framework. Input files were generated using Beauti v.1.5.4 (Drummond and Rambaut, 2007), assuming uncorrelated log-normal trees and a constant population size as prior information. This prior tree is the most suitable for trees describing the relationships between individuals in the same population/species (Drummond and Rambaut, 2007). The nucleotide substitution model and its parameter values were selected based on the aforementioned Modeltest v.3.7 (Posada and Crandall, 1998) results. We carried out two independent runs of 10 million generations. Trees and parameters were sampled every 1,000 iterations, with a burn-in of 10%. Results for each run were visualized using the Tracer v.1.5 program (Rambaut and Drummond, 2009) to ensure stationarity and convergence. Each analysis was repeated many times to optimize the operators of the parameters until no suggestion message appeared in the log file. Posterior probabilities and the maximum credibility tree were calculated using the TreeAnnotator v.1.5.4 program (Drummond and Rambaut, 2007).

Results
The PCR amplifications employing the Dmrt1bY specific oligonucleotides (Nanda et al., 2002) detected two fragments of 1000 and 900 bp in the A. charrua genome. They were cloned, sequenced and named M N1A and M N1B in male, and H N2A and H N2B in female respectively (Figure 1). The alignment of these sequences and highly related ones retrieved from GenBank (Table 2) revealed the Arezo et al. 367  Phylogenetic analysis based on ML and BI phylogenetic tree reconstruction ( Figure 2) showed a topology with three main supported clades integrated by: 1) autosomic Dmrt1 sequences from O. latipes, O. hatcheri and Monopterus albus, Dmrt1bY from O. latipes, dmy from O. curvinotus and the isolated genomic sequences from A. charrua (70% bootstrap support); 2) Dmrt2 sequences from Danio rerio, Xiphophorus maculatus and Takifugu rubripes (100% boostrap support) and 3) Dmrt3 sequences from D. rerio and T. rubripes; Dmrt4 sequences from X. maculatus and T. rubripes genes and Dmrt5 sequences from O. latipes, X. maculatus, D. rerio and T. rubripes (63% boostrap support). Tetraodon nigroviridis Dmrt3 and Dmrt4 resulted in most divergent sequences of this data set which collapsed in a basal polytomy with the aforementioned clades.
The RT-PCR amplifications with the same oligonucleotides (Nanda et al., 2002) revealed an expressed 205 bp cDNA fragment in embryo stages from dispersed phase (late blastula) to pre-hatching embryos. Later on (juveniles and adults stages), the fragment was only detected in male individuals (Figure 3). This sequence showed 8% of identity (e-value = 0.67) with Dmrt1 bY from O. latipes strain HNI (accession number AY129241.1) and dmy (accession number NM001104680.1). The BLAST search was restricted to Oryzias database. The corresponding translated fragment sequence included one stop codon. Remarkably there is an absence of a conserved DM domain in the partial cDNA fragment (205 bp) isolated in A. charrua male.
All phylogenetic analyses (MP, ML and BI) of the isolated partial cDNA sequence with Dmrt1 mRNA sequences from different fish species Mus musculus Dmrt1 as well as D. melanogaster dsx gene were carried out (Figure 4). The tree topology shows two main clades (64% bootstrap support). The major one is subdivided into four subclades embracing: 1) Dmrt1 mRNAs from some Oryzias species ( A. charrua sex determination   ing versions of Dmrt1 from Clarias gariepinus (bootstrap support of 94%). The mRNA Dmrt1 sequence from Mus musculus is the most divergent and not well supported basal taxon which collapsed in a basal polytomy with the remaining aforementioned clades. All tree topology reconstruction methods revealed that the isolated partial cDNA sequence is more closely related to the small fragments amplified from male (M N1B) and female (H N2B) genomes (HPP > 96) and straightforward support the relationships of these sequences with the dsx D. melanogaster gene (HPP > 97) ( Figure 5).

Discussion
Two fragments (M N1A and M N1B) were isolated from the A. charrua male genome. This pattern was similar to O. latipes male where a fragment of 1289 nucleotides corresponds to the Dmrt1 gene and other of 965 nucleotides to the male specific duplicated gene version, Dmrt1bY (Nanda et al., 2002). It is worthy of notice that the same pattern was also found in the A. charrua female genome. Nucleotide sequence comparison of the larger fragments, M N1A and H N2A (1000 bp), showed 42% similarity while the smaller fragments (900 bp), M N1B and H N2B are 99% similar. Thus, the smaller fragment is present in both sexes of the A. charrua genomes. In fish, five classes of Dmrt genes have been described: Dmrt 1, 2, 3, 4 and 5 (Huang et al., 2002). A moderate similarity was detected among the A. charrua genomic fragments and the Dmrt1 sequences. Therefore these results suggest a non conclusive relationship between the isolated fragments and the Dmrt1 from different fish groups.
The absence of a conserved DM domain in the partial cDNA fragment (205 bp) isolated in A. charrua males could be explained by restricted sequence information. The fragment length covers about the 11% of complete Dmrt1bY and Dmrt1 mRNAs reported in the GenBank. An- Arezo et al. 369  (Table 3). Numbers above the branches (left to right) correspond to the Bayesian posterior probability for clades obtained using BEAST 1.5.4v , followed by ML and MP boostrap support (> 60%) respectively, recovered in 1000 replicates. other possibility could be that the DM domain is not present in the A. charrua fragment, as was shown for the Dmrt1c isoform from C. gariepinus and C. batrachus (Raghuveer and Senthilkumaran, 2009) and for the Dmrt1 gonadal isoform from mouse (Lu et al., 2007). Phylogenetic relationships among DM domain containing genes are not obvious in vertebrates, as different taxa show little sequence conservation outside this domain (Volff et al., 2003). Phylogenetic analyses grouped the A. charrua expressed cDNA fragment with Dmrt1 alternative splicing versions c and d from O. latipes (Beloniformes) instead of Dmrt1 sequences from more closely related species of Cyprinodontiformes (K. marmoratus, X. maculatus and P. reticulata). This discordance may be explained by ancestral polymorphism maintenance or homoplasy. Evolution of homoplastic characters results from adaptations of different lineages in response to selective pressures in similar environments (Futuyma, 2005;Wake et al., 2011). Atheriniformes, Beloniformes and Cyprinodontiformes belong to the monophyletic Atherinomorpha Series. Beloniformes and Cyprinodontiformes are considered sister groups, while Atheriniformes occupies a basal position (Parenti, 2005;Nelson, 2006;Setiamarga et al., 2008). Thus, this discordance may be explained by ancestral polymorphism maintenance, as in the topology the A. charrua fragment grouped with O. latipes Dmrt1 c and d alternative splicing versions and Drosophila dsx. Furthermore, this clade is basally related with Odontesthes (Atheriniformes) Dmrt1 sequences. Discordance between Dmrt1 sequence relationships was also evidenced in phylogenetic analyses performed with Dmrt1 from O. bonariensis (Atheriniformes) and other fish species (Fernandino et al., 2006).
Although homology relationships among members of the Dmrt family have not been clearly established in metazoans, it is interesting to consider that the Dmrt1 gene is the closest related member in terms of structure and function to Drosophila dsx and Caenorhadbitis elegans mab-3 (Herpin and Schartl, 2011). In this sense, all phylogenetic analyses in the present work straightforward grouped Drosophila dsx with A. charrua expressed cDNA partial fragment and amplified sequences from male (M N1B) and female (H N2B) genomes showing high posterior probability of occurrence for such clade.
The expression pattern of cDNA partial sequence during the ontogeny of the studied species showed that this sequence is detected from dispersed phase (late blastula). This suggests that its expression depends on the zygotic genome activation (mid-blastula transition). Since no sex specific molecular markers are available to genotype the sex of the embryos, another possibility is that all the embryos analyzed at early stages (1 to 128 cells and early blastulas) were 100% females. In 1 to 128 cell embryos this probability is, however, only 0.004% and in early blastula stage it is 0.008%. The probability of simultaneous occurrence of two or more independent events is calculated as the product of probabilities of each independent event (Canavos, 1988). We assumed the probability for a given embryo to be male or female to be 0.5. For this reason, we consider that the absence of fragment amplification in these developmental stages indicates that it is not maternally supplied and that its expression depends on zygotic genome activation.
mykiss Marchand et al., 2000;K. marmoratus, Kanamori et al., 2006;Halichoeres tenuispinis, Jeong et al., 2009). In D. rerio and Gadus morhua Dmrt1 expression was detected in pre-hatching stages. However, these embryos are in somitogenesis (Schulz et al., 2007;Johnsen and Andersen, 2012). Therefore, the finding that this cDNA partial sequence is expressed very early during A. charrua development and its closely phylogenetic relationships with the dsx_D. melanogaster gene, could be indicative that this sequence could play an essential role during the sex determination process and might be located near the top of the sex determination cascade in this species.
Taken together, these results suggest that both isolated A. charrua genomic sequences and the expressed cDNA partial sequence are more closely related to dsx from D. melanogaster than the two alternative splicing Dmrt1sequences from O.latipes. The observed male-specific expression pattern supports the hypothesis that this sequence probably belongs to an alternative splicing Dmrt1 gene version in A. charrua, as this gene is present in both sexes and represents the most conserved downstream gene member implicated in male sex development during vertebrate evolution. Nevertheless, it is also possible that the isolated cDNA sequence belongs to a class of abundant non-coding RNAs with potential regulatory function. Regulatory RNAs from 100 to 100,000 nucleotides are referred to as long non-coding RNAs (lncRNAs) because they have typical structures. It has furthermore been documented that splicing isoforms encode a protein, whereas other specific splicing isoforms encode regulatory lncRNAs that show tissue specificity and dynamic expression during development (Amaral et al., 2008(Amaral et al., , 2011. The results obtained with RNAfold analysis (Gruber et al., 2008;Lorenz et al., 2011) evidenced that the 205 bp cDNA sequence is capable of generating a stable secondary structure (-84.39 Kcal/mol at 20°C) ( Figure 6). This predicted secondary structure could be formed by lncRNAs, as mRNAs are usually not structured. Nevertheless, it is essential to obtain more sequence information (e.g. stagespecific transcriptomes) to elucidate A. charrua expressed cDNA sequence identity. To further evaluate if this RNA is effectively translated (characterization of expression pattern in temporal and spatial scales) the generation of specific antibodies would be necessary.
For sex determination studies in a given species it is important to consider the criterion proposed by Valenzuela et al. (2003) which emphasizes the presence of sex chromosomes as strong evidence for genetic sex determination (GSD). Sex chromosomes were, however, not identified in A. charrua (García, 2006), in agreement with data obtained from most studied teleosts. Cytogenetically differentiated sex chromosomes are sporadic in different fish taxa, suggesting a recent and polyphyletic origin of sex chromosomes in this vertebrate group (Mank et al., 2006). Fish sex chromosomes are difficult to evidence by classic cytogenetic techniques due to their putative recent origin (ho- Arezo et al. 371 Figure 6 -RNA structure analysis obtained with RNAfold. A) Minimum free energy predicted structure (-84.39 Kcal/mol). B) Centroid secondary structure. In both images color intensity denotes base-pairing probabilities. For unpaired regions color is related to the probability of being unpaired. C) Mountain plot of minimum free energy structure, the thermodynamic ensemble of RNA, centroid structures and entropy for each position. momorphic chromosomes). Ultrastructural analyses of synaptonemic complexes could reveal tiny unpaired regions allowing sex chromosome pair recognition in the heterogametic sex (Devlin and Nagahama, 2002). This approach made it possible to identify sex chromosomes in O. niloticus (Carrasco et al., 1999). For this reason we cannot discard the presence of homomorphic sex chromosomes in A. charrua. Even in the absence of sex chromosomes, GSD may be demonstrated by the construction of a microsatellite-based linkage map to define genome regions involved in sex determination. This approach allowed the identification of the male sex determining non-recombinant region in the African annual fish Nothobranchius furzeri (Valenzano et al., 2009), similar to those identified in other phylogenetically related fish species as O. latipes (Matsuda et al., 2002;Nanda et al., 2002), Gasterosteus aculeatus (Peichel et al., 2004), X. maculatus (Schultheis et al., 2009) and P. reticulata (Tripathi et al., 2009).