BACKGROUND
sex was thought to have arisen in a single evolutive event in the last common ancestor of all eukaryotes, since sexual reproduction is almost universal and exclusive of this group  <cit> .
considering its consequences over the lifespan of an organism and its influence on population demography, it is thought that the sex determination  mechanism should be under strong selection forces  <cit> .
however, sex can be established by many different and fast-evolving mechanisms  <cit> , indicating that sd triggers have emerged several times throughout evolution  <cit> .
within vertebrates, different sex determining systems have been described.
in therian mammals, with a xx/xy chromosome system, sex depends on the presence of the sry gene, a paralogue of sox <dig>  on the y chromosome  <cit> , while in birds  with a zz/zw chromosome system, the dmrt <dig> gene with a double dosage is required for testis development  <cit> .
also, in xenopus laevis the dm-w gene, a paralogue of dmrt <dig>  is responsible for sd  <cit> .
when dmy, belonging to the dm family of transcription factors like dmrt <dig> and dm-w, was found to be the sex determining gene  of the fish oryzias latipes, a biased and recurrent recruitment of specific sdgs or families throughout evolution was suggested  <cit> .
however, later findings in fish do not seem to support this hypothesis.
fish, with around  <dig>  species  <cit> , is the most diverse group of vertebrates and its study has broadened our knowledge on sd.
fish diversity is also reflected by the variety of reproductive strategies: unisexuality, different types of hermaphroditism and gonochorism; and also by the diversity of sd mechanisms  <cit> .
in the last years, an important effort has been made in order to identify the sdg in several model and aquaculture fish species.
different productive traits are sex-associated in farm fish such as growth rate, color, taste and flesh quality; hence, the interest of the industry in producing monosex stocks  <cit> .
detailed information at gene level is available for only a limited number of fish species.
five different master sdgs have been identified so far: dmy /dmrt1by in oryzias latipes and in o. curvinotus  <cit> , gsdf in o. luzonensis  <cit> , amhy in odontesthes hatchery  <cit> , amhr <dig> in takifugu rubripes, t. pardalis and t. poecilonotus  <cit> , and sdy in salmonid family  <cit> .
recently, a distant cis-regulatory element of sox <dig> necessary for male determination in o. dancena, a species with a xx/xy sd system, has also been identified  <cit> , and dmrt <dig> has been suggested as the sdg in cynoglossus semilaevis  <cit> .
however, little information is available, not only on the sdg, but also on the initial molecular pathways related to sexual differentiation.
traditionally, sd has been related to the switching mechanism of a hierarchical genetic network that causes the activation of downstream genes involved in gonad differentiation  leading to the differentiation of testes or ovaries  <cit> .
thus, concerning whether the first difference between sexes is the expression of a gene/s or the strength of an environmental factor, sd can be genetic or environmental, although both ways can coexist  <cit> .
in the classical view of sd and gd, the downstream genes of the cascade were assumed to be highly conserved, and only the genes at the top of the cascade would change by gene duplication  or by allelic diversification, establishing a new sd mechanism  <cit> .
nowadays, the conservation of the downstream cascade has been questioned  <cit>  and a new view, which considers sex as a threshold phenotype in which both genetic and environmental factors can act alone or in combination and, importantly, at different times during the period of gd is gaining support  <cit> .
in this new view, sd encompasses not only the initial trigger, be it genetic, environmental or both, but also the whole gd process, and different factors such as cell proliferation and hormone levels would be involved in determining a threshold which would give rise to a testis or an ovary, thus fitting to a threshold quantitative trait  <cit> .
turbot is one of the most important species cultured in europe, being galicia  the main production region since the eighties.
production and quality of farmed fish rely on a deep knowledge of biological functions, especially those related to reproduction, growth and disease resistance.
in this context, the production of monosex stocks to exploit sex-associated dimorphisms related to productive traits, especially growth and sexual maturation, has been a sought in finfish aquaculture by different approaches  <cit> .
turbot shows one of the strongest sexual growth dimorphisms amongst marine species and females can reach up to 50 % bigger size than males  <cit> , thus industry is interested in the production of all-female populations.
in the last years, an important effort has been devoted to understanding sd and gd in this species.
analysis on mitotic and meiotic chromosomes revealed the absence of an heteromorphic sex chromosome pair related to sex  <cit> .
the major sd region was located on the linkage group   <dig> at  <dig>  cm of sma-usc <dig> marker   <cit> , but other minor sex-related quantitative trait locus  were detected at lg <dig>  lg <dig> and lg <dig>  <cit> .
in that study, a zz/zw chromosome system was established in accordance with the sex ratios of progenies obtained from hormonally sex-reversed parents  <cit> .
temperature also showed a minor influence on sex ratios in this species  <cit> .
close to the sex-associated marker sma-usc <dig> several candidate genes were identified , but eventually discarded because no association to sex was detected at the species level, so the sd gene remains unidentified in turbot  <cit> .
considering the lack of information on the sd mechanism of turbot, we determined the expression patterns of a suite of  <dig> genes shown to be involved in gd in other species at the initial critical stages of sex differentiation using a large amount of fish and sampling times in turbot.
in a previous study, we analyzed reproduction in this species through the use of oligo-microarrays, spanning a larger age period but in a lower number of samples and finding several genes involved in ovary or testis development .
in this work, our aim was to study gene expression at a very specific time point coinciding with the onset of sex differentiation.
we also evaluated the effect of temperature to clarify whether it has a major role in turbot sd and its possible interaction with genetic factors.
our results enabled us to establish the correspondence between gene expression profiles and histological sex and to identify the first sex-associated genetic cues at the initial stages of gonad development in turbot.
RESULTS
sampling and sexing
a timeline of turbot gonad differentiation is shown in fig.
<dig>  all of the  <dig> turbot samples used in this study were genetically sexed using the sma-usc <dig> marker and, additionally, the  <dig>   <dig> and  <dig> dpf samples were histologically sexed.
eighty-nine females and eighty-five males could be genetically sexed because the sma-usc <dig> marker was informative, the remaining six samples being removed from this analysis since they could not be sexed.
a 7 % sexing discrepancy was observed between the genetic and histological information in the samples obtained at  <dig>   <dig> and  <dig> dpf.
sma-usc <dig> is expected to present a 2 % sex-genotyping error  <cit> , so in our data the 5 % sex allocation error may represent sex-reversal events due the action of a secondary qtl or environmental factors.
however, given the reasonable accuracy of genetic sexing, the sex of samples below  <dig> dpf obtained through smausc-e <dig> genotyping was considered for further analyses.fig.
<dig> turbot gonad differentiation.
main histological and physiological events along sexual differentiation in turbot.
turbot images were obtained from http://larvalbase.org/and http://fishbase.se/
males and females did not show length differences both in the whole dataset and at each age  .
however, significant differences were found at every age between the three temperatures, except at  <dig> and  <dig> dpf, where 18 °c-reared animals did not differ in length from those reared at 23 °c .table  <dig> length comparison between temperatures for each age group
co-localization of targeted genes with sex-related qtls
a main sd qtl in linkage group  <dig>  and three minor ones in lg <dig>  lg <dig> and lg <dig> were previously reported in turbot  <cit> .
after establishing the relationship between the turbot map  and the turbot genome ,  <dig> genes in the subset analyzed here  could be located in lgs harboring a sd qtl .
five genes were found in lg <dig> and two of them, sox <dig> and fxr <dig>  co-localized with the main sd qtl; ar <dig> co-localized with the sex qtl in lg8; four genes were placed in lg <dig> and two of them, sox9a and sox <dig>  within the confidence interval of the sd qtl.fig.
<dig> turbot sex qtls and target genes.
four turbot linkage groups are shown.
estimated location of the target genes is shown in red.
grey shaded lg areas represent the position of the sd qtls  <cit>
global expression patterns
samples and genes were hierarchically clustered in a heatmap using the pearson correlation coefficient as distance measure .
for each sample, rearing temperature, age and sex are shown in the heatmap.
some samples are grouped according to sex or age and can be associated with particular groups of coexpressed genes.
the samples of  <dig> and  <dig> dpf are clustered in two groups  and they are characterized by the high expression of sox <dig>  fxr <dig>  wnt <dig>  hsp <dig>  ptges <dig>  t4_ <dig> and dmrt <dig> , but also by the nearly null expression of tdrd <dig>  vasa, cyp19a1a, foxl <dig> and gsdf, involved in gonad maturation and female differentiation .
these samples are grouped by age independently of temperature or sex, which do not seem to represent relevant factors on the diagnostic genes expressed at these ages.
two different groups of older fish, one essentially made of females and another essentially made of males, could clearly be identified.
the female group  is mainly associated with the up-regulation of two different clusters of genes: one cluster containing cyp19a1a, foxl <dig>  vasa, tdrd <dig> and gsdf, genes not expressed in undifferentiated individuals, as outlined before; and another cluster containing dnmt <dig>  dact <dig>  sox <dig>  rdh <dig> and ctnnb <dig> .
the male group  is associated with the expression of sox <dig>  amh, ar <dig>  fshb, cyp11a .
these “male” genes are also highly expressed in a mix of males and females of around  <dig> dpf and mostly reared at low temperatures .
some “female” genes  are also expressed to a lower extent in the male samples, suggesting a role in gonad development irrespective of sex .
as previously mentioned, some blocks connected to rearing temperature can also be seen, but, in general, it does not seem to be a determining factor for sample clustering.fig.
<dig> global expression patterns.
a heatmap of target genes and all gonad samples.
gene names are shown in the bottom of the figure while gene hierarchical cluster is shown in the top.
log fold change expression values representation ranges from red  to light green .
sample names are not shown, instead each sample is represented by the three colors at the left of the figure which indicate sex , age  and temperature .
yellow, black or blue circles highlight expression patterns characteristic of undifferentiated, female or male individuals respectively.
b samples were grouped according to the fold change expression values of the target genes by a principal component analysis.
samples labeled as “female” and “male” and colored in red and purple, respectively, represent gonad samples which were both genetically and histologically sexed because developing testis and ovaries could be distinguished from each other.
samples labeled as “genfemale” and “genmale” and colored in olive green and light blue, respectively, are gonad samples which were only genetically sexed because histologically the gonads were still undifferentiated.
a 66 % of the samples of each group are expected to be placed in their respective circles.
the arrows with the name of the genes at the end represent how each gene contributes to the two principal analysis components represented in the figure
pca analysis on the  <dig> samples  revealed that histologically sexed males and females , could be in most cases clearly discriminated by their differential expression.
however, younger individuals , appeared fully overlapped in the pca, likely because they are still undifferentiated.
a 66 % prediction ellipse for each group is shown in fig.
3b indicating that if new individuals were added to our analysis from a certain group, 66 % of them would expect to be placed inside the corresponding ellipse.
some of the genetically sexed individuals  are also found in the ellipses drawn for phenotypic males or females which, also considering the previous heatmap results, indicates that turbot gd might start before  <dig> dpf.
interestingly, about half a dozen genetic males were included in the female circle.
the arrows indicate the weight of each gene on the two first principal components.
among the analyzed genes there is a large group seemingly contributing to female differentiation , while the presumed male-related genes are fewer and not so markedly pointing towards male differentiation .
also, as suggested in the heatmap, some genes are clearly related to undifferentiated individuals .
sex differences
fold change  expression values for those differentially expressed genes between males and females  were analyzed in relation to age  .fig.
<dig> gene fold change values along gonad development.
cyp19a1a, vasa, amh, foxl <dig>  tdrd <dig> and gsdf fold change values for each sample plotted according to both its length, in cm, and its age, in days post fertilization.
female samples are shown in magenta and male samples in blue.
in the fc/length figure for each gene non-linear trend lines were calculated by loess regression and genes with significative differences between sexes at any age point present a pink background if the gene is up-regulated in females or a blue one if it is up-regulated in males.
genes without sex differences have a white background.
in the fc/age figure, error bars represent the standard error of the mean, also an asterisk marks those age points were the differences in expression between males and females are significant
fig.
<dig> gene fold change values along gonad development.
sox <dig>  dnmt <dig>  dact <dig>  rdh <dig>  ctnnb <dig>  sf <dig> and piwi <dig> fold change values for each sample plotted according to both its length, in cm, and its age, in days post fertilization.
female samples are shown in magenta and male samples in blue.
in the fc/length figure for each gene non-linear trend lines were calculated by loess regression and genes with significative differences between sexes at any age point present a pink background if the gene is up-regulated in females or a blue one if it is up-regulated in males.
genes without sex differences have a white background.
in the fc/age figure, error bars represent the standard error of the mean, also an asterisk marks those age points were the differences in expression between males and females are significant
fig.
<dig> gene fold change values along gonad development.
fxr <dig>  dmrt <dig>  zar <dig>  wnt <dig>  sox9a, sox <dig> fold change values for each sample plotted according to both its length, in cm, and its age, in days post fertilization.
female samples are shown in magenta and male samples in blue.
in the fc/length figure for each gene non-linear trend lines were calculated by loess regression and genes with significative differences between sexes at any age point present a pink background if the gene is up-regulated in females or a blue one if it is up-regulated in males.
genes without sex differences have a white background.
in the fc/age figure, error bars represent the standard error of the mean, also an asterisk marks those age points were the differences in expression between males and females are significant
sex differences were first observed at  <dig> dpf, when cyp19a1a  and vasa   which were more expressed in females, and amh  in males .
these three genes also increased their expression from  <dig> to  <dig> dpf, in both sexes , which is also observed for foxl <dig>  tdrd <dig>  gsdf  and sox <dig> and rdh <dig> .
gsdf expression increased dramatically even earlier, from  <dig> to  <dig> dpf .
foxl <dig> did not show significant differences in expression in females until  <dig> dpf, mainly because its expression decreased slowly in males .
tdrd <dig> pattern of expression resembled that of vasa, being more highly expressed in females at  <dig> dpf and onwards.
other genes increased their expression at  <dig> dpf  in females, while their expression remained low or even decreased in males.
fold change values of females vs. males at  <dig> dpf were above  <dig> for sox <dig> , dnmt <dig> , dact <dig> , rdh <dig> , ctnnb <dig> , sf <dig>  and piwi <dig>   and these differences increased from  <dig> to  <dig> dpf, when the fcs were between  <dig> and  <dig> for the seven genes.
fxr <dig> was also up-regulated in females at  <dig> dpf  and showed a 47 % increase at  <dig> dpf , but presented higher expression levels in undifferentiated individuals, irrespective of sex.
a similar pattern was also observed for wnt <dig>  dmrt <dig> and zar <dig> , genes which at some point during sex differentiation, 105– <dig> dpf , showed a higher expression in females, but its expression decreased from  <dig> to  <dig> dpf .
in contrast to that observed in females, in our study there were very few genes in the subset assayed whose expression were higher in males and, even in these cases, the differences between males and females were low.
one of them was sox9a for which sex differences increased from  <dig> dpf  up to  <dig> dpf  .
another sox family gene, sox <dig>  highly expressed in undifferentiated individuals in both sexes, was up-regulated in males at  <dig> dpf , but the difference remained constant at  <dig> dpf .
two additional genes, fshb and cyp11a, showed mean expression values slightly higher in males, but not significant .
some genes like ptges <dig>  hh <dig>  hsp <dig> or t4_ <dig> did not show sex differences .
among these, a gene of the sox family, sox <dig>  showed some groups of outliers whose expression was not explained either by sex or by length/age.
other two genes of this family, sox <dig> and sox <dig>  did not present any clear expression pattern along development or by sex, and androgen receptor  <dig>  ar <dig> did not show sex dimorphic expression either.
discriminant analysis
a discriminant analysis considering the earliest dimorphic expressed genes -cyp19a1a, amh and vasa– correctly classified 100 % of the genetic males and 82 % of the genetic females at  <dig> dpf, representing as a whole 91 % of individuals correctly classified .
the remaining 18 % of the genetic females were grouped as males.
these three genes can constitute a fairly efficient genetic tool for early sexing of turbot.
furthermore, from  <dig> dpf onwards , the expression of cyp19a1a alone is capable of perfectly discriminating males and females without error.fig.
<dig> discriminant analysis 3d plot.
ninety days post fertilization samples were plotted in a three dimensions graph according to their fold change values for cyp19a1a, amh and vasa.
female samples are colored in magenta and males in blue
network analysis
to further understand the relationships between genes we performed a network analysis based on gene-to-gene correlations .
a tight cluster with several female up-regulated genes  was found with all their genes inter-connected.
also, sox <dig> and amh constituted a small male-like cluster together with fshb and cyp11a.
the two clusters were connected through two genes: fxr <dig>  gene located at the main sex determining region of turbot, and gsdf.
the absence of some genes in the network  suggests that they did not show significant relationships with any other gene in our analysis and for the chosen correlation threshold.
this does not mean that they do not have any role in sex differentiation, since our study analyzed the expression of a limited number of genes .
if more genes were added, it is possible that these genes could have shown connection to the network through them.fig.
<dig> network representation.
weighed correlation network performed with the fold change expression values of the genes is shown.
genes are represented as blue circles if they are up-regulated in males at any age, magenta if they are up-regulated in females, or dark green if no differences were found.
lines connecting genes indicate significant correlations, red lines are positive correlations and blue lines are negative correlations
temperature effects on gene expression
we found a higher proportion of females at both 15 °c and 18 °c than at 23 °c, where the male:female proportion was close to 1: <dig> .
however, these differences were not significant .
genotyping of the sd marker in the main sd qtl in lg <dig> strongly suggested that some genetic males developed as females: a total of  <dig> genetic males were classified as females by histology.
this discrepancy can either be caused by 1) temperature effects on sd or 2) interaction of a secondary qtl  with the main sd qtl in lg <dig>  further, we analyzed the effects of temperature on gene expression in males and females separately and, since turbot length was different between temperatures in almost every development stage, we checked if the detected temperature differences in gene expression were independent of length  or not .
among those genes with length-independent temperature effects on expression , only sox <dig> showed temperature effects which are not sex dependent.
this gene showed higher expression at 15 and 23 °c, although the difference between  <dig> and 23 °c was not significant in females.
among those genes showing sex-specific temperature effects amh, sox9a and cyp11a were more highly expressed at low temperatures in females, while sox <dig> and dmrt <dig> showed the opposite pattern with higher expression at 23 °c in females.
on the contrary, ctnnb <dig>  piwi <dig>  sf <dig> and sox <dig> were up-regulated at low temperatures in males, and the four genes showed a very similar pattern.fig.
<dig> sex proportions and temperature.
percentage of male and female turbot, histologically sexed, at  <dig>   <dig> and 23 °c.
also, the percentage of phenotypic females which are genetic males is shown.
no genetic females developed as males.
sample size =  <dig> fish per temperature
fig.
<dig> temperature effects on gene expression.
mean fold change gene expression values at  <dig>   <dig> and 23 °c in the whole dataset are shown for males  and females .
error bars represent standard deviation.
significate differences between temperatures are indicated by *, **, *** or ****.
black asterisks indicate that the detected difference is independent of fish length, while red asterisks indicate that fish length has an influence on the temperature differences
finally, since  <dig> dpf is the first age where we find sex dimorphic expression and can be considered the onset of sex differentiation, we decided to test for temperature differences at this age which were independent of both length and sex, since most genes did not show dimorphic gene expression at this age.
a single gene showed significant differences between temperatures at  <dig> dpf: ctnnb <dig> .
ctnnb <dig> highest expression was observed at 18 °c, while its lowest expression was at 23 °c.
discussion
recently, sex determination and differentiation have begun to be perceived as a single modular process rather than a cascade, where sex behaves in many instances as a threshold-like character  <cit> .
under this view, gonad fate depends on several factors acting coordinately, including many genes and possibly also environmental variables.
in this sense, and given the interest in obtaining monosex stocks in many finfish aquaculture species, understanding the different gene patterns during early sex differentiation, supposedly the moment where gonad fate can be more easily controlled or altered, is of great importance to manipulate sex determination in order to control sex ratios.
in this study we have analyzed the expression of  <dig> genes previously connected to sex differentiation in other species, some of them studied in fish regarding sex differentiation for the first time.
a total of  <dig> genes were found to show dimorphic expression at some point during early sex differentiation in turbot.
the influence of temperature was also assessed, finding differences between temperatures for  <dig> genes.
this study has broadened our knowledge of gene expression patterns during early sex determination in turbot in particular and in fish in general.
early sex differentiation and primodial germ cells
although morphological gonad differences between sexes were not detected before  <dig> dpf, the first molecular signs of sex differentiation were observed between  <dig> and  <dig> dpf , characterized by an expression increase of gsdf and tdrd <dig>  but also by the expression increase and differential expression of cyp19a1a, amh and vasa, which allowed discriminating males  and females  at  <dig> dpf.
vasa, tdrd <dig> and gsdf are genes related to primordial germ cell development.
tdrd <dig> proteins were detected in the primordial germ cells of zebrafish  at  <dig> dpf and were involved in both oocyte and sperm development  <cit> .
tdrd <dig> was found to interact with vasa, which is also a highly specific marker of germ cells  required for their development  <cit>  and conserved along several invertebrate and vertebrate taxa  <cit> .
while vasa and tdrd <dig> are specific germ cell markers, gsdf is a teleost-specific member of the tgf-β superfamily which has been reported to be expressed in the somatic cells surrounding the primordial germ cells in rainbow trout and promoting their proliferation  <cit> .
gsdf has shown higher levels of expression in testis in zebrafish, three-spot wrasse  and coelacanth  .
furthermore, gsdf has been proposed as the male sex determining gene in anoploma fimbria  <cit> , and a copy named gsdfy has been found to be the sex determinant gene in oryzias luzonensis, with a higher expression in males 10 days after hatching  <cit> .
gsdf has also been found to be directly up-regulated by dmy, male sdg of oryzias latipes  <cit> , and sox <dig>  male sdg of oryzias dancena  <cit>  during the first stages of gd.
although gsdf does not seem to have such a male-like function in turbot, since it did not present sex dimorphic expression, it seems to be important for gonad development in both sexes since its expression greatly increased between  <dig> and  <dig> dpf.
even more, this is the only gene which showed a significant expression increase from  <dig> to  <dig> dpf in our study, consistent with a function as germ cell inductor since its expression increase precedes that of vasa or tdrd <dig>
the expression patterns of vasa, tdrd <dig> and gsdf in turbot suggested that primordial germ cells start proliferating between  <dig> and  <dig> dpf in both sexes, either faster in females or suffering a certain delay in males, as suggested by vasa and tdrd <dig> expression levels.
another germ-line specific gene, piwi <dig>  was investigated in our study, but it showed a different expression pattern and its induction was limited to females and delayed until  <dig> dpf.
in zebrafish, piwi <dig> has been found to play a crucial role in meiosis  <cit>  and perhaps its different pattern of expression in turbot may be related to the onset of meiosis in female germ cells.
the amount of primordial cells is recognized as one of the initial differences between male and female gonads in some fish species like zebrafish  <cit> , medaka  <cit>  and stickleback   <cit> , although not in others such as loach   <cit>  and goldfish   <cit> .
further, germ cell proliferation has been found to be associated with sd in several fish species.
in medaka, germ cell proliferation is inhibited in males when the sex determining gene, dmy, is expressed at the sd stage before testis differentiation  <cit> .
when dmy is not active in xy embryos, germ cells proliferate and enter meiosis as in xx embryos.
surprisingly, amh does not present a sex dimorphic expression in medaka during gd  <cit> .
although the müllerian ducts are not present in modern teleosts  <cit> , amh orthologs have been described in fish and characterized as gonad specific key factors for male sex differentiation  <cit> .
yet, even if amh does not present sex dimorphic expression in medaka, it has been found to control germ cell proliferation in this species, and mutations on its receptor  lead to excessive proliferation of germ cells which caused male-to-female sex reversal  <cit> , although phenotypic female xy gonads still expressed dmy.
furthermore, if these amhrii mutants are depleted of germ cells, testis development takes place  <cit> .
thus, amh seems to be a repressor of germ cell proliferation in medaka and to play a major role in sex determination.
this is also the case in fugu , where a single snp in the coding region of amhrii is likely responsible for sd  <cit> .
this snp encodes a protein with a reduced function and is fixed in females, which are not sensitive to amh  <cit> .
also, in the patagonian pejerrey, a copy of amh, amhy, has been found to be the sex determining gene and its action has been suggested to regulate germ cell proliferation, being upstream to the autosomal amh and relegating the former to a function in testicular maturation and/or spermatogenesis  <cit> .
interestingly, in the female gonad of zebrafish the maintenance of cyp19a1a expression, but not its activation, has been related to the presence of the primordial germ cells  <cit> , pointing towards a model where amh is responsible for the control of germ cell proliferation while germ cells aid to maintain cyp19a1a expression.
amh and cyp19a1a are genes involved in sex differentiation across all vertebrate taxa  <cit> , and have been reported as male and female like genes, respectively, in several fish species  <cit> .
in our experiment, cyp19a1a showed higher expression in females and amh in males at  <dig> dpf.
a threshold expression of amh could be controlling sexual fate.
if amh does not reach the required expression level, primordial germ cells will proliferate and maintain cyp19a1a levels while amh levels decrease.
on the contrary, if amh expression reaches a certain threshold, germ cells stop proliferating and cyp19a1a expression decreases.
although further data is required, this hypothesis seems to be consistent with the findings in turbot and other fish species, and could be a common mechanism controlling the balance between male and female gonad differentiation in fish.
female sex differentiation
turbot gd is in progress at 105– <dig> dpf, either towards males or females, and sex was easily identified by cyp19a1a expression alone at these developmental stages.
several female up-regulated genes were detected during this period of gd , while male gene expression pattern was very similar to undifferentiated fish for the assayed genes, excluding a few classical male-like genes .
foxl <dig> is a transcription factor that activates cyp19a1a transcription  <cit>  and both genes are strictly co-expressed in mammals  <cit> .
the expression pattern of foxl <dig> was similar to that previously described for cyp19a1a, vasa, tdrd <dig> and gsdf.
its expression in turbot is consistent with an activation of cyp19a1a; however, the later decrease of cyp19a1a while foxl <dig> expression is still high in males suggests other roles for foxl <dig> not related to cyp19a1a activation at early stages of development.
foxl <dig> expression has also been described in the male gonad of tilapia , southern catfish  and goldfish  <cit> .
also, cells with cyp19a1a expression without foxl <dig> expression have been reported in medaka  <cit> .
although the authors did not exclude that those cells had earlier foxl <dig> expression, at least foxl <dig> did not seem to be essential for cyp19a1a expression maintenance  <cit> .
foxl <dig> seems to be involved in the initiation of the female gonad development cascade through the activation of cyp19a1a expression in fish, as previously suggested  <cit> , however foxl <dig> up-regulation alone is not sufficient for the activation of cyp19a1a and it might have additional roles important for gonad development both in males and females.
ctnnb <dig> is another female-like gene with a critical function in female differentiation in mammals, antagonizing sox <dig> and blocking testis development, thus promoting ovarian development  <cit> .
ctnnb <dig> is the key downstream component of the canonical wnt signaling pathway and our results are in agreement with a conserved female function of this pathway, which has also been shown in zebrafish  <cit>  and rainbow trout  <cit> .
however, in turbot, the wnt pathway seems to be at least partially independent of wnt <dig>  the highest expression of this gene was found at undifferentiated stages, consistent with a role in early gonad development in both sexes also observed in mammals, although with a different function since in mammals it is responsible for the development of müllerian ducts  <cit> .
wnt <dig> expression dropped at  <dig> dpf irrespective of sex, and later was found scarcely up-regulated in females.
wnt <dig> is a key female gene in mammals which antagonizes fgf <dig> and down-regulates sox <dig> expression  <cit> .
wnt <dig> not so clear female-pattern during gd in turbot is likely related to the absence of fgf <dig> in teleosts  <cit> .
no dimorphic wnt <dig> expression has been observed in rainbow trout   <cit> , zebrafish  <cit>  or the more distant rana rugosa  <cit> .
wnt <dig> does not show a conserved function in female sd throughout evolution, and the results in our study suggest that it is not involved in the expression of ctnnb <dig> in the female gonad development and therefore other wnt proteins could be responsible for activating the wnt pathway, which seems to have a conserved female prominent function.
six genes, dnmt <dig>  rdh <dig>  sox <dig>  dact <dig>  ctnnb <dig> and sf <dig> showed a similar increasing expression pattern in females from  <dig> to  <dig> dpf.
these genes showed high pair-wise correlation values and also the highest negative correlations with amh amongst all the assayed genes.
hence, amh down-regulation may be required for the activation of genes involved in female gonad development.
the coactivation of these genes in female differentiation does not seem to be related to a specific pathway considering their functional diversity, but with the activation of several concomitant pathways at the beginning of ovarian development, and their up-regulation is possibly connected to the down-regulation of amh.
male-like genes
besides amh, we found two up-regulated genes in males from  <dig> dpf onwards: sox9a and sox <dig>  furthermore, the heatmap and network analyses suggested that cyp11a and fshb are related to male development as well, but they did not display a dimorphic expression pattern at any stage.
sox9a is an essential player in sex differentiation and its male-like nature seems to be rather conserved throughout evolution, but its relevance seems to be variable  <cit> .
in mammals, this gene is directly activated by sry and is responsible and sufficient for fating the male gonad  <cit> , also activating amh transcription  <cit> .
in birds, sox <dig> is co-expressed with amh and induced by the sd gene dmrt <dig>  <cit> .
however, in medaka, sox <dig> is not required for testis development  <cit> .
in turbot, sox9a dimorphic expression is found later than that of amh, and its expression is more stable along the assayed stages, suggesting a less important role in gd in this species.
genes in the main sd region
fxr <dig> is located in the turbot genome very close to sma-usce <dig>  <cit> , the highest sex-associated marker in this species, within the main sd region at lg <dig>  <cit> , thus representing a potential sd candidate gene.
in this experiment, fxr <dig> is highly expressed before the first gd signs and precedes the expression of cyp19a1a and amh, although at this time its expression is not sexually dimorphic.
fxr <dig> is a rna-binding protein and is an autosomal paralogue of fmrp , important for normal female reproductive function and cognition development in humans.
further, fxr <dig> it has been related to female gametogenesis in pigs  <cit>  and xenopus laevis  <cit> .
in turbot, fxr <dig> position in the functional network, connecting male and female clusters, and its higher expression at low temperatures, suggest that further analysis on this gene in younger individuals looking for a putative dimorphic pattern between sexes or other type of dimorphic cue are desirable.
on the contrary sox <dig>  also located in the main sd region of turbot  <cit> , did not present a dimorphic expression pattern and also showed a relatively steady expression along all the assayed stages, suggesting no role of this gene in turbot sd.
temperature effects
a higher proportion of phenotypic females were obtained at  <dig> and 18 °c, however sex ratio differences between temperatures were not significant.
this effect was not due to genetic segregation distortion, since genotyping of the sd marker strongly suggested that some genetic males developed as females.
although proper male-to-female sex reversal at low temperatures has not been described in fish  <cit> , it has been described in reptiles, with genetic males developing as females  <cit> .
the influence of temperature on sex ratio in turbot was previously reported  <cit> , but it was family-dependent and not always in the same way: in this study, two families presented a higher proportion of females than expected at 23 °c while another family presented more females at 15 °c.
in any case, if temperature effects exist in turbot, they seem to be limited and family dependent, however further work should be done to evaluate this issue.
another possible explanation for genetic males developing as phenotypic females is the existence of a family-specific genetic interaction between a secondary sd qtl and the main sd qtl in lg <dig>  which could be responsible for genetic males developing as phenotypic females.
independently of the mechanism, families with higher proportions of females could be very interesting for aquaculture, since it would lead to higher growth rates given the important sex size dimorphism of turbot.
nevertheless, temperature effects on gene expression are very important to explain possible sex ratio shifts and other effects on sex differentiation.
we found temperature effects on the expression of some of the genes assayed: amh, ctnnb <dig>  cyp11a, dmrt <dig>  piwi <dig>  sf <dig>  sox <dig>  sox <dig>  sox9a and sox <dig>  among these, only sox <dig> temperature effects were sex independent.
this gene showed higher expression both at  <dig> and 23 °c, which may indicate some kind of stress response.
sox <dig> is a transcription factor regulating several genes and it is also involved in the maintenance of stem-cell identity  <cit> .
among the genes which presented sex dependent temperature effects, ctnnb <dig>  piwi <dig>  sf <dig> and sox <dig> showed higher expression at low temperatures in males.
among these, ctnnb <dig> is remarkable because it not only occupies a downstream position in the wnt pathway  <cit> , which needs to be up-regulated for the development of an ovary in zebrafish  <cit> , but also because it showed temperature-dependent expression differences in our study at  <dig> dpf, the onset of sd.
genes found upstream in the wnt signalling pathway are likely regulated by temperature and responsible for this increase in ctnnb <dig> expression in turbot.
consistently with our results, elevated ctnnb <dig> expression has been reported connected to low temperatures in rats  <cit>  and tilapia  <cit> , albeit in other tissues.
a recent study in oyster found a biased sex ratio towards females related to higher ctnnb <dig> expression at lower rearing temperatures  <cit> .
thus, ctnnb <dig> and the wnt pathway are good candidates for future studies aimed at investigating temperature effects on sex ratios in turbot and possibly other fish.
other two genes showed the opposite pattern in females, with higher expression at higher temperatures: dmrt <dig> and sox <dig>  dmrt <dig> was upregulated during gonad development and also expressed in germ cells in the swamp eel   <cit> , while in rana rugosa it was expressed in the developing gonad during sd without any dimorphic pattern, suggesting a function both in testicular and ovarian differentiation  <cit> .
sox <dig> has been associated with ovarian development in dicentrarchus labrax  <cit> , although it did not present dimorphic expression in turbot and so, apparently, it is not related to female differentiation in this species.
yet, this gene seemed to have a peak of expression at  <dig> dpf in both sexes, so it could have some function in early gonad development in turbot for both males and females.
in the swamp eel, sox <dig> was also expressed both in testis and ovary  <cit> .
temperature effects are gene and sex specific.
given the labile nature of sd in fish related to specific morphogenetic thresholds, several genes could be responsible for sex ratio shifts.
as seen in this study, several genes involved in sex show expression differences due to temperature, and so these genes are potential candidates for sex ratio alterations.
CONCLUSIONS
turbot sex differentiation is ongoing at  <dig> dpf and sex can be distinguished by the expression levels of three genes when fish are 5–6 cm length: cyp19a1a, amh and vasa; while later females are easily discriminated by the expression of cyp19a1a.
the first molecular signs of sex differentiation are the dimorphic expression of these three genes and an increase in the expression of vasa, gsdf and tdrd <dig>  connected with primordial germ cells, suggesting their proliferation from  <dig> to  <dig> dpf and an important role in sex differentiation.
the primary sex determining gene of turbot remains unknown, since neither sox <dig> nor fxr <dig>  genes located in the main sd region of turbot, showed expression patterns clearly consistent with this role.
our data suggest that female development has more complex machinery and is strongly regulated, suggesting the involvement of both methylation and splicing mechanisms.
furthermore, we have observed that temperature affects the expression of several genes and suggest that the wnt/β-catenin pathway could be a likely candidate to explain possible temperature-induced sex ratio shifts.
the present study, analyzing the expression pattern of many genes related to sex differentiation, has revealed that turbot sex differentiation is a complex process with many factors involved.
these results are more compatible with a view of sex determination as a network where the activation or repression of several genes can affect gonad fate.
this view of sex determination as a threshold character could help us to understand temperature effects during sex differentiation.

 sex was thought to have arisen in a single evolutive event in the last common ancestor of all eukaryotes, since sexual reproduction is almost universal and exclusive of this group. considering its consequences over the lifespan of an organism and its influence on population demography, it is thought that the sex determination mechanism should be under strong selection forces. however, sex can be established by many different and fast-evolving mechanisms, indicating that sd triggers have emerged several times throughout evolution. within vertebrates, different sex determining systems have been described. in therian mammals, with a xx/xy chromosome system, sex depends on the presence of the sry gene, a paralogue of sox on the y chromosome, while in birds with a zz/zw chromosome system, the dmrt gene with a double dosage is required for testis development. also, in xenopus laevis the dm-w gene, a paralogue of dmrt is responsible for sd. when dmy, belonging to the dm family of transcription factors like dmrt and dm-w, was found to be the sex determining gene of the fish oryzias latipes, a biased and recurrent recruitment of specific sdgs or families throughout evolution was suggested. however, later findings in fish do not seem to support this hypothesis. fish, with around species, is the most diverse group of vertebrates and its study has broadened our knowledge on sd. fish diversity is also reflected by the variety of reproductive strategies: unisexuality, different types of hermaphroditism and gonochorism; and also by the diversity of sd mechanisms. in the last years, an important effort has been made in order to identify the sdg in several model and aquaculture fish species. different productive traits are sex-associated in farm fish such as growth rate, color, taste and flesh quality; hence, the interest of the industry in producing monosex stocks. detailed information at gene level is available for only a limited number of fish species. five different master sdgs have been identified so far: dmy /dmrt1by in oryzias latipes and in o. curvinotus, gsdf in o. luzonensis, amhy in odontesthes hatchery, amhr in takifugu rubripes, t. pardalis and t. poecilonotus, and sdy in salmonid family. recently, a distant cis-regulatory element of sox necessary for male determination in o. dancena, a species with a xx/xy sd system, has also been identified, and dmrt has been suggested as the sdg in cynoglossus semilaevis. however, little information is available, not only on the sdg, but also on the initial molecular pathways related to sexual differentiation. traditionally, sd has been related to the switching mechanism of a hierarchical genetic network that causes the activation of downstream genes involved in gonad differentiation leading to the differentiation of testes or ovaries. thus, concerning whether the first difference between sexes is the expression of a gene/s or the strength of an environmental factor, sd can be genetic or environmental, although both ways can coexist. in the classical view of sd and gd, the downstream genes of the cascade were assumed to be highly conserved, and only the genes at the top of the cascade would change by gene duplication or by allelic diversification, establishing a new sd mechanism. nowadays, the conservation of the downstream cascade has been questioned and a new view, which considers sex as a threshold phenotype in which both genetic and environmental factors can act alone or in combination and, importantly, at different times during the period of gd is gaining support. in this new view, sd encompasses not only the initial trigger, be it genetic, environmental or both, but also the whole gd process, and different factors such as cell proliferation and hormone levels would be involved in determining a threshold which would give rise to a testis or an ovary, thus fitting to a threshold quantitative trait. turbot is one of the most important species cultured in europe, being galicia the main production region since the eighties. production and quality of farmed fish rely on a deep knowledge of biological functions, especially those related to reproduction, growth and disease resistance. in this context, the production of monosex stocks to exploit sex-associated dimorphisms related to productive traits, especially growth and sexual maturation, has been a sought in finfish aquaculture by different approaches. turbot shows one of the strongest sexual growth dimorphisms amongst marine species and females can reach up to 50 % bigger size than males, thus industry is interested in the production of all-female populations. in the last years, an important effort has been devoted to understanding sd and gd in this species. analysis on mitotic and meiotic chromosomes revealed the absence of an heteromorphic sex chromosome pair related to sex. the major sd region was located on the linkage group at  cm of sma-usc marker, but other minor sex-related quantitative trait locus were detected at lg lg and lg. in that study, a zz/zw chromosome system was established in accordance with the sex ratios of progenies obtained from hormonally sex-reversed parents. temperature also showed a minor influence on sex ratios in this species. close to the sex-associated marker sma-usc several candidate genes were identified, but eventually discarded because no association to sex was detected at the species level, so the sd gene remains unidentified in turbot. considering the lack of information on the sd mechanism of turbot, we determined the expression patterns of a suite of genes shown to be involved in gd in other species at the initial critical stages of sex differentiation using a large amount of fish and sampling times in turbot. in a previous study, we analyzed reproduction in this species through the use of oligo-microarrays, spanning a larger age period but in a lower number of samples and finding several genes involved in ovary or testis development. in this work, our aim was to study gene expression at a very specific time point coinciding with the onset of sex differentiation. we also evaluated the effect of temperature to clarify whether it has a major role in turbot sd and its possible interaction with genetic factors. our results enabled us to establish the correspondence between gene expression profiles and histological sex and to identify the first sex-associated genetic cues at the initial stages of gonad development in turbot.  sampling and sexing a timeline of turbot gonad differentiation is shown in fig.  all of the turbot samples used in this study were genetically sexed using the sma-usc marker and, additionally, the and dpf samples were histologically sexed. eighty-nine females and eighty-five males could be genetically sexed because the sma-usc marker was informative, the remaining six samples being removed from this analysis since they could not be sexed. a 7 % sexing discrepancy was observed between the genetic and histological information in the samples obtained at and dpf. sma-usc is expected to present a 2 % sex-genotyping error, so in our data the 5 % sex allocation error may represent sex-reversal events due the action of a secondary qtl or environmental factors. however, given the reasonable accuracy of genetic sexing, the sex of samples below dpf obtained through smausc-e genotyping was considered for further analyses.fig. turbot gonad differentiation. main histological and physiological events along sexual differentiation in turbot. turbot images were obtained from http://larvalbase.org/and http://fishbase.se/ males and females did not show length differences both in the whole dataset and at each age. however, significant differences were found at every age between the three temperatures, except at and dpf, where 18 °c-reared animals did not differ in length from those reared at 23 °c.table length comparison between temperatures for each age group co-localization of targeted genes with sex-related qtls a main sd qtl in linkage group and three minor ones in lg lg and lg were previously reported in turbot. after establishing the relationship between the turbot map and the turbot genome, genes in the subset analyzed here could be located in lgs harboring a sd qtl. five genes were found in lg and two of them, sox and fxr co-localized with the main sd qtl; ar co-localized with the sex qtl in lg8; four genes were placed in lg and two of them, sox9a and sox within the confidence interval of the sd qtl.fig. turbot sex qtls and target genes. four turbot linkage groups are shown. estimated location of the target genes is shown in red. grey shaded lg areas represent the position of the sd qtls global expression patterns samples and genes were hierarchically clustered in a heatmap using the pearson correlation coefficient as distance measure. for each sample, rearing temperature, age and sex are shown in the heatmap. some samples are grouped according to sex or age and can be associated with particular groups of coexpressed genes. the samples of and dpf are clustered in two groups and they are characterized by the high expression of sox fxr wnt hsp ptges t4_ and dmrt, but also by the nearly null expression of tdrd vasa, cyp19a1a, foxl and gsdf, involved in gonad maturation and female differentiation. these samples are grouped by age independently of temperature or sex, which do not seem to represent relevant factors on the diagnostic genes expressed at these ages. two different groups of older fish, one essentially made of females and another essentially made of males, could clearly be identified. the female group is mainly associated with the up-regulation of two different clusters of genes: one cluster containing cyp19a1a, foxl vasa, tdrd and gsdf, genes not expressed in undifferentiated individuals, as outlined before; and another cluster containing dnmt dact sox rdh and ctnnb. the male group is associated with the expression of sox amh, ar fshb, cyp11a. these “male” genes are also highly expressed in a mix of males and females of around dpf and mostly reared at low temperatures. some “female” genes are also expressed to a lower extent in the male samples, suggesting a role in gonad development irrespective of sex. as previously mentioned, some blocks connected to rearing temperature can also be seen, but, in general, it does not seem to be a determining factor for sample clustering.fig. global expression patterns. a heatmap of target genes and all gonad samples. gene names are shown in the bottom of the figure while gene hierarchical cluster is shown in the top. log fold change expression values representation ranges from red to light green. sample names are not shown, instead each sample is represented by the three colors at the left of the figure which indicate sex, age and temperature. yellow, black or blue circles highlight expression patterns characteristic of undifferentiated, female or male individuals respectively. b samples were grouped according to the fold change expression values of the target genes by a principal component analysis. samples labeled as “female” and “male” and colored in red and purple, respectively, represent gonad samples which were both genetically and histologically sexed because developing testis and ovaries could be distinguished from each other. samples labeled as “genfemale” and “genmale” and colored in olive green and light blue, respectively, are gonad samples which were only genetically sexed because histologically the gonads were still undifferentiated. a 66 % of the samples of each group are expected to be placed in their respective circles. the arrows with the name of the genes at the end represent how each gene contributes to the two principal analysis components represented in the figure pca analysis on the samples revealed that histologically sexed males and females, could be in most cases clearly discriminated by their differential expression. however, younger individuals, appeared fully overlapped in the pca, likely because they are still undifferentiated. a 66 % prediction ellipse for each group is shown in fig. 3b indicating that if new individuals were added to our analysis from a certain group, 66 % of them would expect to be placed inside the corresponding ellipse. some of the genetically sexed individuals are also found in the ellipses drawn for phenotypic males or females which, also considering the previous heatmap results, indicates that turbot gd might start before dpf. interestingly, about half a dozen genetic males were included in the female circle. the arrows indicate the weight of each gene on the two first principal components. among the analyzed genes there is a large group seemingly contributing to female differentiation, while the presumed male-related genes are fewer and not so markedly pointing towards male differentiation. also, as suggested in the heatmap, some genes are clearly related to undifferentiated individuals. sex differences fold change expression values for those differentially expressed genes between males and females were analyzed in relation to age.fig. gene fold change values along gonad development. cyp19a1a, vasa, amh, foxl tdrd and gsdf fold change values for each sample plotted according to both its length, in cm, and its age, in days post fertilization. female samples are shown in magenta and male samples in blue. in the fc/length figure for each gene non-linear trend lines were calculated by loess regression and genes with significative differences between sexes at any age point present a pink background if the gene is up-regulated in females or a blue one if it is up-regulated in males. genes without sex differences have a white background. in the fc/age figure, error bars represent the standard error of the mean, also an asterisk marks those age points were the differences in expression between males and females are significant fig. gene fold change values along gonad development. sox dnmt dact rdh ctnnb sf and piwi fold change values for each sample plotted according to both its length, in cm, and its age, in days post fertilization. female samples are shown in magenta and male samples in blue. in the fc/length figure for each gene non-linear trend lines were calculated by loess regression and genes with significative differences between sexes at any age point present a pink background if the gene is up-regulated in females or a blue one if it is up-regulated in males. genes without sex differences have a white background. in the fc/age figure, error bars represent the standard error of the mean, also an asterisk marks those age points were the differences in expression between males and females are significant fig. gene fold change values along gonad development. fxr dmrt zar wnt sox9a, sox fold change values for each sample plotted according to both its length, in cm, and its age, in days post fertilization. female samples are shown in magenta and male samples in blue. in the fc/length figure for each gene non-linear trend lines were calculated by loess regression and genes with significative differences between sexes at any age point present a pink background if the gene is up-regulated in females or a blue one if it is up-regulated in males. genes without sex differences have a white background. in the fc/age figure, error bars represent the standard error of the mean, also an asterisk marks those age points were the differences in expression between males and females are significant sex differences were first observed at dpf, when cyp19a1a and vasa which were more expressed in females, and amh in males. these three genes also increased their expression from to dpf, in both sexes, which is also observed for foxl tdrd gsdf and sox and rdh. gsdf expression increased dramatically even earlier, from to dpf. foxl did not show significant differences in expression in females until dpf, mainly because its expression decreased slowly in males. tdrd pattern of expression resembled that of vasa, being more highly expressed in females at dpf and onwards. other genes increased their expression at dpf in females, while their expression remained low or even decreased in males. fold change values of females vs. males at dpf were above for sox, dnmt, dact, rdh, ctnnb, sf and piwi and these differences increased from to dpf, when the fcs were between and for the seven genes. fxr was also up-regulated in females at dpf and showed a 47 % increase at dpf, but presented higher expression levels in undifferentiated individuals, irrespective of sex. a similar pattern was also observed for wnt dmrt and zar, genes which at some point during sex differentiation, 105– dpf, showed a higher expression in females, but its expression decreased from to dpf. in contrast to that observed in females, in our study there were very few genes in the subset assayed whose expression were higher in males and, even in these cases, the differences between males and females were low. one of them was sox9a for which sex differences increased from dpf up to dpf. another sox family gene, sox highly expressed in undifferentiated individuals in both sexes, was up-regulated in males at dpf, but the difference remained constant at dpf. two additional genes, fshb and cyp11a, showed mean expression values slightly higher in males, but not significant. some genes like ptges hh hsp or t4_ did not show sex differences. among these, a gene of the sox family, sox showed some groups of outliers whose expression was not explained either by sex or by length/age. other two genes of this family, sox and sox did not present any clear expression pattern along development or by sex, and androgen receptor ar did not show sex dimorphic expression either. discriminant analysis a discriminant analysis considering the earliest dimorphic expressed genes -cyp19a1a, amh and vasa– correctly classified 100 % of the genetic males and 82 % of the genetic females at dpf, representing as a whole 91 % of individuals correctly classified. the remaining 18 % of the genetic females were grouped as males. these three genes can constitute a fairly efficient genetic tool for early sexing of turbot. furthermore, from dpf onwards, the expression of cyp19a1a alone is capable of perfectly discriminating males and females without error.fig. discriminant analysis 3d plot. ninety days post fertilization samples were plotted in a three dimensions graph according to their fold change values for cyp19a1a, amh and vasa. female samples are colored in magenta and males in blue network analysis to further understand the relationships between genes we performed a network analysis based on gene-to-gene correlations. a tight cluster with several female up-regulated genes was found with all their genes inter-connected. also, sox and amh constituted a small male-like cluster together with fshb and cyp11a. the two clusters were connected through two genes: fxr gene located at the main sex determining region of turbot, and gsdf. the absence of some genes in the network suggests that they did not show significant relationships with any other gene in our analysis and for the chosen correlation threshold. this does not mean that they do not have any role in sex differentiation, since our study analyzed the expression of a limited number of genes. if more genes were added, it is possible that these genes could have shown connection to the network through them.fig. network representation. weighed correlation network performed with the fold change expression values of the genes is shown. genes are represented as blue circles if they are up-regulated in males at any age, magenta if they are up-regulated in females, or dark green if no differences were found. lines connecting genes indicate significant correlations, red lines are positive correlations and blue lines are negative correlations temperature effects on gene expression we found a higher proportion of females at both 15 °c and 18 °c than at 23 °c, where the male:female proportion was close to 1:. however, these differences were not significant. genotyping of the sd marker in the main sd qtl in lg strongly suggested that some genetic males developed as females: a total of genetic males were classified as females by histology. this discrepancy can either be caused by 1) temperature effects on sd or 2) interaction of a secondary qtl with the main sd qtl in lg further, we analyzed the effects of temperature on gene expression in males and females separately and, since turbot length was different between temperatures in almost every development stage, we checked if the detected temperature differences in gene expression were independent of length or not. among those genes with length-independent temperature effects on expression, only sox showed temperature effects which are not sex dependent. this gene showed higher expression at 15 and 23 °c, although the difference between and 23 °c was not significant in females. among those genes showing sex-specific temperature effects amh, sox9a and cyp11a were more highly expressed at low temperatures in females, while sox and dmrt showed the opposite pattern with higher expression at 23 °c in females. on the contrary, ctnnb piwi sf and sox were up-regulated at low temperatures in males, and the four genes showed a very similar pattern.fig. sex proportions and temperature. percentage of male and female turbot, histologically sexed, at and 23 °c. also, the percentage of phenotypic females which are genetic males is shown. no genetic females developed as males. sample size =  fish per temperature fig. temperature effects on gene expression. mean fold change gene expression values at and 23 °c in the whole dataset are shown for males and females. error bars represent standard deviation. significate differences between temperatures are indicated by *, **, *** or ****. black asterisks indicate that the detected difference is independent of fish length, while red asterisks indicate that fish length has an influence on the temperature differences finally, since dpf is the first age where we find sex dimorphic expression and can be considered the onset of sex differentiation, we decided to test for temperature differences at this age which were independent of both length and sex, since most genes did not show dimorphic gene expression at this age. a single gene showed significant differences between temperatures at dpf: ctnnb. ctnnb highest expression was observed at 18 °c, while its lowest expression was at 23 °c. discussion recently, sex determination and differentiation have begun to be perceived as a single modular process rather than a cascade, where sex behaves in many instances as a threshold-like character. under this view, gonad fate depends on several factors acting coordinately, including many genes and possibly also environmental variables. in this sense, and given the interest in obtaining monosex stocks in many finfish aquaculture species, understanding the different gene patterns during early sex differentiation, supposedly the moment where gonad fate can be more easily controlled or altered, is of great importance to manipulate sex determination in order to control sex ratios. in this study we have analyzed the expression of genes previously connected to sex differentiation in other species, some of them studied in fish regarding sex differentiation for the first time. a total of genes were found to show dimorphic expression at some point during early sex differentiation in turbot. the influence of temperature was also assessed, finding differences between temperatures for genes. this study has broadened our knowledge of gene expression patterns during early sex determination in turbot in particular and in fish in general. early sex differentiation and primodial germ cells although morphological gonad differences between sexes were not detected before dpf, the first molecular signs of sex differentiation were observed between and dpf, characterized by an expression increase of gsdf and tdrd but also by the expression increase and differential expression of cyp19a1a, amh and vasa, which allowed discriminating males and females at dpf. vasa, tdrd and gsdf are genes related to primordial germ cell development. tdrd proteins were detected in the primordial germ cells of zebrafish at dpf and were involved in both oocyte and sperm development. tdrd was found to interact with vasa, which is also a highly specific marker of germ cells required for their development and conserved along several invertebrate and vertebrate taxa. while vasa and tdrd are specific germ cell markers, gsdf is a teleost-specific member of the tgf-β superfamily which has been reported to be expressed in the somatic cells surrounding the primordial germ cells in rainbow trout and promoting their proliferation. gsdf has shown higher levels of expression in testis in zebrafish, three-spot wrasse and coelacanth. furthermore, gsdf has been proposed as the male sex determining gene in anoploma fimbria, and a copy named gsdfy has been found to be the sex determinant gene in oryzias luzonensis, with a higher expression in males 10 days after hatching. gsdf has also been found to be directly up-regulated by dmy, male sdg of oryzias latipes, and sox male sdg of oryzias dancena during the first stages of gd. although gsdf does not seem to have such a male-like function in turbot, since it did not present sex dimorphic expression, it seems to be important for gonad development in both sexes since its expression greatly increased between and dpf. even more, this is the only gene which showed a significant expression increase from to dpf in our study, consistent with a function as germ cell inductor since its expression increase precedes that of vasa or tdrd the expression patterns of vasa, tdrd and gsdf in turbot suggested that primordial germ cells start proliferating between and dpf in both sexes, either faster in females or suffering a certain delay in males, as suggested by vasa and tdrd expression levels. another germ-line specific gene, piwi was investigated in our study, but it showed a different expression pattern and its induction was limited to females and delayed until dpf. in zebrafish, piwi has been found to play a crucial role in meiosis and perhaps its different pattern of expression in turbot may be related to the onset of meiosis in female germ cells. the amount of primordial cells is recognized as one of the initial differences between male and female gonads in some fish species like zebrafish, medaka and stickleback, although not in others such as loach and goldfish. further, germ cell proliferation has been found to be associated with sd in several fish species. in medaka, germ cell proliferation is inhibited in males when the sex determining gene, dmy, is expressed at the sd stage before testis differentiation. when dmy is not active in xy embryos, germ cells proliferate and enter meiosis as in xx embryos. surprisingly, amh does not present a sex dimorphic expression in medaka during gd. although the müllerian ducts are not present in modern teleosts, amh orthologs have been described in fish and characterized as gonad specific key factors for male sex differentiation. yet, even if amh does not present sex dimorphic expression in medaka, it has been found to control germ cell proliferation in this species, and mutations on its receptor lead to excessive proliferation of germ cells which caused male-to-female sex reversal, although phenotypic female xy gonads still expressed dmy. furthermore, if these amhrii mutants are depleted of germ cells, testis development takes place. thus, amh seems to be a repressor of germ cell proliferation in medaka and to play a major role in sex determination. this is also the case in fugu, where a single snp in the coding region of amhrii is likely responsible for sd. this snp encodes a protein with a reduced function and is fixed in females, which are not sensitive to amh. also, in the patagonian pejerrey, a copy of amh, amhy, has been found to be the sex determining gene and its action has been suggested to regulate germ cell proliferation, being upstream to the autosomal amh and relegating the former to a function in testicular maturation and/or spermatogenesis. interestingly, in the female gonad of zebrafish the maintenance of cyp19a1a expression, but not its activation, has been related to the presence of the primordial germ cells, pointing towards a model where amh is responsible for the control of germ cell proliferation while germ cells aid to maintain cyp19a1a expression. amh and cyp19a1a are genes involved in sex differentiation across all vertebrate taxa, and have been reported as male and female like genes, respectively, in several fish species. in our experiment, cyp19a1a showed higher expression in females and amh in males at dpf. a threshold expression of amh could be controlling sexual fate. if amh does not reach the required expression level, primordial germ cells will proliferate and maintain cyp19a1a levels while amh levels decrease. on the contrary, if amh expression reaches a certain threshold, germ cells stop proliferating and cyp19a1a expression decreases. although further data is required, this hypothesis seems to be consistent with the findings in turbot and other fish species, and could be a common mechanism controlling the balance between male and female gonad differentiation in fish. female sex differentiation turbot gd is in progress at 105– dpf, either towards males or females, and sex was easily identified by cyp19a1a expression alone at these developmental stages. several female up-regulated genes were detected during this period of gd, while male gene expression pattern was very similar to undifferentiated fish for the assayed genes, excluding a few classical male-like genes. foxl is a transcription factor that activates cyp19a1a transcription and both genes are strictly co-expressed in mammals. the expression pattern of foxl was similar to that previously described for cyp19a1a, vasa, tdrd and gsdf. its expression in turbot is consistent with an activation of cyp19a1a; however, the later decrease of cyp19a1a while foxl expression is still high in males suggests other roles for foxl not related to cyp19a1a activation at early stages of development. foxl expression has also been described in the male gonad of tilapia, southern catfish and goldfish. also, cells with cyp19a1a expression without foxl expression have been reported in medaka. although the authors did not exclude that those cells had earlier foxl expression, at least foxl did not seem to be essential for cyp19a1a expression maintenance. foxl seems to be involved in the initiation of the female gonad development cascade through the activation of cyp19a1a expression in fish, as previously suggested, however foxl up-regulation alone is not sufficient for the activation of cyp19a1a and it might have additional roles important for gonad development both in males and females. ctnnb is another female-like gene with a critical function in female differentiation in mammals, antagonizing sox and blocking testis development, thus promoting ovarian development. ctnnb is the key downstream component of the canonical wnt signaling pathway and our results are in agreement with a conserved female function of this pathway, which has also been shown in zebrafish and rainbow trout. however, in turbot, the wnt pathway seems to be at least partially independent of wnt the highest expression of this gene was found at undifferentiated stages, consistent with a role in early gonad development in both sexes also observed in mammals, although with a different function since in mammals it is responsible for the development of müllerian ducts. wnt expression dropped at dpf irrespective of sex, and later was found scarcely up-regulated in females. wnt is a key female gene in mammals which antagonizes fgf and down-regulates sox expression. wnt not so clear female-pattern during gd in turbot is likely related to the absence of fgf in teleosts. no dimorphic wnt expression has been observed in rainbow trout, zebrafish or the more distant rana rugosa. wnt does not show a conserved function in female sd throughout evolution, and the results in our study suggest that it is not involved in the expression of ctnnb in the female gonad development and therefore other wnt proteins could be responsible for activating the wnt pathway, which seems to have a conserved female prominent function. six genes, dnmt rdh sox dact ctnnb and sf showed a similar increasing expression pattern in females from to dpf. these genes showed high pair-wise correlation values and also the highest negative correlations with amh amongst all the assayed genes. hence, amh down-regulation may be required for the activation of genes involved in female gonad development. the coactivation of these genes in female differentiation does not seem to be related to a specific pathway considering their functional diversity, but with the activation of several concomitant pathways at the beginning of ovarian development, and their up-regulation is possibly connected to the down-regulation of amh. male-like genes besides amh, we found two up-regulated genes in males from dpf onwards: sox9a and sox furthermore, the heatmap and network analyses suggested that cyp11a and fshb are related to male development as well, but they did not display a dimorphic expression pattern at any stage. sox9a is an essential player in sex differentiation and its male-like nature seems to be rather conserved throughout evolution, but its relevance seems to be variable. in mammals, this gene is directly activated by sry and is responsible and sufficient for fating the male gonad, also activating amh transcription. in birds, sox is co-expressed with amh and induced by the sd gene dmrt. however, in medaka, sox is not required for testis development. in turbot, sox9a dimorphic expression is found later than that of amh, and its expression is more stable along the assayed stages, suggesting a less important role in gd in this species. genes in the main sd region fxr is located in the turbot genome very close to sma-usce, the highest sex-associated marker in this species, within the main sd region at lg, thus representing a potential sd candidate gene. in this experiment, fxr is highly expressed before the first gd signs and precedes the expression of cyp19a1a and amh, although at this time its expression is not sexually dimorphic. fxr is a rna-binding protein and is an autosomal paralogue of fmrp, important for normal female reproductive function and cognition development in humans. further, fxr it has been related to female gametogenesis in pigs and xenopus laevis. in turbot, fxr position in the functional network, connecting male and female clusters, and its higher expression at low temperatures, suggest that further analysis on this gene in younger individuals looking for a putative dimorphic pattern between sexes or other type of dimorphic cue are desirable. on the contrary sox also located in the main sd region of turbot, did not present a dimorphic expression pattern and also showed a relatively steady expression along all the assayed stages, suggesting no role of this gene in turbot sd. temperature effects a higher proportion of phenotypic females were obtained at and 18 °c, however sex ratio differences between temperatures were not significant. this effect was not due to genetic segregation distortion, since genotyping of the sd marker strongly suggested that some genetic males developed as females. although proper male-to-female sex reversal at low temperatures has not been described in fish, it has been described in reptiles, with genetic males developing as females. the influence of temperature on sex ratio in turbot was previously reported, but it was family-dependent and not always in the same way: in this study, two families presented a higher proportion of females than expected at 23 °c while another family presented more females at 15 °c. in any case, if temperature effects exist in turbot, they seem to be limited and family dependent, however further work should be done to evaluate this issue. another possible explanation for genetic males developing as phenotypic females is the existence of a family-specific genetic interaction between a secondary sd qtl and the main sd qtl in lg which could be responsible for genetic males developing as phenotypic females. independently of the mechanism, families with higher proportions of females could be very interesting for aquaculture, since it would lead to higher growth rates given the important sex size dimorphism of turbot. nevertheless, temperature effects on gene expression are very important to explain possible sex ratio shifts and other effects on sex differentiation. we found temperature effects on the expression of some of the genes assayed: amh, ctnnb cyp11a, dmrt piwi sf sox sox sox9a and sox among these, only sox temperature effects were sex independent. this gene showed higher expression both at and 23 °c, which may indicate some kind of stress response. sox is a transcription factor regulating several genes and it is also involved in the maintenance of stem-cell identity. among the genes which presented sex dependent temperature effects, ctnnb piwi sf and sox showed higher expression at low temperatures in males. among these, ctnnb is remarkable because it not only occupies a downstream position in the wnt pathway, which needs to be up-regulated for the development of an ovary in zebrafish, but also because it showed temperature-dependent expression differences in our study at dpf, the onset of sd. genes found upstream in the wnt signalling pathway are likely regulated by temperature and responsible for this increase in ctnnb expression in turbot. consistently with our results, elevated ctnnb expression has been reported connected to low temperatures in rats and tilapia, albeit in other tissues. a recent study in oyster found a biased sex ratio towards females related to higher ctnnb expression at lower rearing temperatures. thus, ctnnb and the wnt pathway are good candidates for future studies aimed at investigating temperature effects on sex ratios in turbot and possibly other fish. other two genes showed the opposite pattern in females, with higher expression at higher temperatures: dmrt and sox dmrt was upregulated during gonad development and also expressed in germ cells in the swamp eel, while in rana rugosa it was expressed in the developing gonad during sd without any dimorphic pattern, suggesting a function both in testicular and ovarian differentiation. sox has been associated with ovarian development in dicentrarchus labrax, although it did not present dimorphic expression in turbot and so, apparently, it is not related to female differentiation in this species. yet, this gene seemed to have a peak of expression at dpf in both sexes, so it could have some function in early gonad development in turbot for both males and females. in the swamp eel, sox was also expressed both in testis and ovary. temperature effects are gene and sex specific. given the labile nature of sd in fish related to specific morphogenetic thresholds, several genes could be responsible for sex ratio shifts. as seen in this study, several genes involved in sex show expression differences due to temperature, and so these genes are potential candidates for sex ratio alterations.  turbot sex differentiation is ongoing at dpf and sex can be distinguished by the expression levels of three genes when fish are 5–6 cm length: cyp19a1a, amh and vasa; while later females are easily discriminated by the expression of cyp19a1a. the first molecular signs of sex differentiation are the dimorphic expression of these three genes and an increase in the expression of vasa, gsdf and tdrd connected with primordial germ cells, suggesting their proliferation from to dpf and an important role in sex differentiation. the primary sex determining gene of turbot remains unknown, since neither sox nor fxr genes located in the main sd region of turbot, showed expression patterns clearly consistent with this role. our data suggest that female development has more complex machinery and is strongly regulated, suggesting the involvement of both methylation and splicing mechanisms. furthermore, we have observed that temperature affects the expression of several genes and suggest that the wnt/β-catenin pathway could be a likely candidate to explain possible temperature-induced sex ratio shifts. the present study, analyzing the expression pattern of many genes related to sex differentiation, has revealed that turbot sex differentiation is a complex process with many factors involved. these results are more compatible with a view of sex determination as a network where the activation or repression of several genes can affect gonad fate. this view of sex determination as a threshold character could help us to understand temperature effects during sex differentiation. 
although important knowledge has been achieved on the genetic basis of sex determination in this species, the master sd gene remains unknown and precise information on gene expression at the critical stage of sex differentiation is lacking. in the present work, we examined the expression profiles of relevant genes related to sex differentiation, from up to 135 days post fertilization, when gonads are differentiating. we also considered the influence of three temperature regimes on sex differentiation. the first sex-related differences in molecular markers could be observed at 90 days post fertilization and so we have called that time the onset of sex differentiation. three genes were the first to show differential expression between males and females and also allowed us to sex turbot accurately at the onset of sex differentiation: cyp19a1a, amh and vasa. furthermore, we profiled several genes involved in sex differentiation and found specific temperature effects on their expression. the expression of genes related to primordial germ cells started to increase between 75– dpf and vasa and tdrd later presented higher expression in females. our results enabled us to identify the first sex-associated genetic cues at the initial stages of gonad development in turbot and to accurately sex turbot at this age, establishing the correspondence between gene expression profiles and histological sex. we also detected changes in the expression levels of several genes depending on culture temperature. controlling sex ratios is essential for the aquaculture industry, especially in those species with sex dimorphism for relevant productive traits, hence the importance of knowing how the sexual phenotype is established in fish.
 controlling sex ratios is essential for the aquaculture industry, especially in those species with sex dimorphism for relevant productive traits, hence the importance of knowing how the sexual phenotype is established in fish. turbot, a very important fish for the aquaculture industry in europe, shows one of the largest sexual growth dimorphisms amongst marine cultured species, being all-female stocks a desirable goal for the industry. although important knowledge has been achieved on the genetic basis of sex determination in this species, the master sd gene remains unknown and precise information on gene expression at the critical stage of sex differentiation is lacking. in the present work, we examined the expression profiles of relevant genes related to sex differentiation, from up to 135 days post fertilization, when gonads are differentiating. we also considered the influence of three temperature regimes on sex differentiation.  the first sex-related differences in molecular markers could be observed at 90 days post fertilization and so we have called that time the onset of sex differentiation. three genes were the first to show differential expression between males and females and also allowed us to sex turbot accurately at the onset of sex differentiation: cyp19a1a, amh and vasa. the expression of genes related to primordial germ cells started to increase between 75– dpf and vasa and tdrd later presented higher expression in females. two genes placed on the sd region of turbot did not show any expression pattern suggestive of a sex determining function. we also detected changes in the expression levels of several genes depending on culture temperature.  our results enabled us to identify the first sex-associated genetic cues at the initial stages of gonad development in turbot and to accurately sex turbot at this age, establishing the correspondence between gene expression profiles and histological sex. furthermore, we profiled several genes involved in sex differentiation and found specific temperature effects on their expression. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords fishgonaddevelopmentqpcrgenessex ratioaromatasemalefemaletemperaturespanish ministerio de ciencia e innovaciónsánchez laura spanish governement consolider ingenio aquagenomicscsd2007-00002martínez paulino issue-copyright-statement© the author 2015 