Gene Expression of Chicken Gonads Is Sex- and Side-Specific

Many mechanisms of avian sex determination are still unknown. Although sex in birds and mammals is genetically determined, the mammalian and avian sex chromosomes are not homologous [Fridolfsson et al., 1998]. In contrast to mammals, female birds are heterogametic (ZW), while males are homogametic (ZZ). It is not clear whether it is the presence of the W chromosome in females, the double representation of the Z chromosome in males vis-à-vis females, or both of these characteristics that are crucial for the determination of sex in birds [Ellegren, 2000].

In every gonochoristic species, genetic pathways are differentially activated in males and females in order to initiate testis or ovary development [Merchant-Larios and Moreno-Mendoza, 2001; Ayers et al., 2013]. These sex-biased genes are often linked to reproduction, rendering many of them subject not only to natural selection but as well to the pressures of sexual [Civetta and Singh, 1999; Swanson and Vacquier, 2002] and anthropogenic selection (e.g. hormone-mimicking chemicals in food and environment). Sex-biased gene expression is relatively common in metazoans [Ranz et al., 2003; Marinotti et al., 2006; Yang et al., 2006; Eads et al., 2007] and has important evolutionary [Meiklejohn et al., 2003; Zhang et al., 2004; Ellegren and Parsch, 2007; Mank et al., 2007], medical [Ivakine et al., 2005; Lu et al., 2007; Xu et al., 2008], and genomic implications [Connallon and Knowles, 2005, 2007; Hambuch and Parsch, 2005]. Strong right-left differences in females, however, are restricted to birds - except for birds of prey [Jacob and Bakst, 2007] - and characterized by a regressing right ovary with a medulla only, while the left one develops ongoing from day 7 of incubation to a functional ovary with cortex and medulla [Ukeshima and Fujimoto, 1991; Ukeshima, 1996; Smith, 2007; González-Morán, 2011]. Previous studies on avian gonad development mainly focused on early embryonic events based on histological descriptions [González-Morán, 2011] or on mRNA expression profiles of a distinct subset of genes [Hoshino et al., 2005; Hudson et al., 2005a, b; Rodríguez-Léon et al., 2008; Carré et al., 2011] as e.g. those involved in steroidogenesis, paracrine signaling, transcription, and homeostasis [Yamamoto et al., 2003; Diaz et al., 2011]. PITX2 and RALDH2 were shown to be important for ovarian dimorphism in the chicken embryo [Rodríguez-Léon et al., 2008; Smith et al., 2008a]. However, the molecular consequences of this differential expression are far from being completely understood. At present, the publication of the chicken genome [Wicker et al., 2005] along with the advent of next-generation sequencing (NGS)-based genome-wide gene expression profiling techniques provides the opportunity to monitor the differential gene expression of developing chicken gonads and to infer common and diverted mechanisms of gonadal differentiation from these data.

DeepSuperSAGE, a genome-wide, open-architecture transcriptome profiling technique, is the adaptation of SuperSAGE [Matsumura et al., 2003] to NGS [Matsumura et al., 2010, 2011; Zawada et al., 2011]. For the generation of comprehensive transcription profiles, so-called SuperTAGs, originating from the 3′ end of a given cDNA, are sequenced in their millions and annotated to their respective genes in an according database [Matsumura et al., 2003, 2005]. These SuperTAGs contain sufficient information to reliably characterize the cDNA and thus the transcript they are derived from [Velculescu et al., 1995; Saha et al., 2002]. Contrary to microarrays, deepSuperSAGE due to its open-architecture is able to quantify the expression of virtually all polyadenylated RNAs, including previously unknown and lowly expressed ones as well as natural antisense transcripts and polyadenylated non-coding RNAs without prior knowledge of the respective sequences.

There is only rudimentary research regarding the female right gonad, since most of the research was done focusing on the left gonads of both sexes. The present study provides a comprehensive description not only of the sex but especially of the side-specific differential sexual development of male and female chicken gonads on the level of gene expression at embryonic day 19 (E19). We wanted to see if there is any difference in male right-left gene expression in E19 and if there is any parallelism to the female gonads. Our genome-wide characterization of differences in the gene expression of mature gonads in an unprecedented depth represents an important step to a molecular characterization of side-specific gonadal development in birds and has implications beyond this taxon for other vertebrate classes, including mammals.

Newly laid fertile chicken eggs (White Leghorn, deriving from the same genetic background) were obtained from a commercial local supplier (LSL Rhein-Main, Schaafheim, Germany) and bred in a ThermoStar 100 egg incubator (J. Hemel Brutgeräte, Verl, Germany). Incubation was performed at 37.6°C and 60% humidity and egg turning once every 2 h. On E19, 2 days before anticipated hatching, embryos were decapitated immediately after removal from the egg, dissected under a microscope, and sexed by gonad morphology.

For genetic sex confirmation, blood samples from each embryo were collected in ethanol and stored at -20°C till isolation of the DNA with the DNeasy (Qiagen, Hilden, Germany) isolation kit. Genetic sexing was carried out by a standard PCR protocol using the primers 2550F (5′-GTT ACT GAT TCG TCT ACG AGA-3′) and 2718R (5′-ATT GAA ATG ATC CAG TGC TTG-3′) previously described by Fridolfsson and Ellegren [1999]. These primers are focused on the CHD1 introns, placed on the Z (CHD1Z, 600 bp) and W chromosome (CHD1W, 450 bp). Thermal cycling comprised DNA polymerase activation at 95°C for 1 min followed by 30 cycles of denaturation at 95°C for 30 s, annealing at 52°C for 30 s, elongation at 72°C for 1 min, and a final extension step at 72°C for 3 min. All amplifications were performed on an advanced primus 96 thermocycler (Peqlab, Erlangen, Germany). The amplicons of this modified protocol were separated into 1 band (Z) in the case of male or 2 bands (Z + W) in the case of female animals on a 1.4% agarose gel.

Gonads were separated, cleaned from adhesive tissue, and immediately frozen individually in 200 µl RNA lysis buffer (Promega, Mannheim, Germany) for later RNA isolation. Total RNA was extracted using the SV Total RNA Isolation System Kit according to manual 048 (Promega). Deviating from the protocol, on-column DNaseI digestion of genomic DNA was elongated from original 15 min to 30 min. Additionally, DNaseI digestion was carried out with Baseline-Zero™ DNase (Epicentre; provided by Biozym Scientific GmbH, Hessisch Oldendorf, Germany) in solution to ensure that the samples were completely free of DNA. Total RNA concentration was estimated in a dilution series with the LabelGuard NanoPhotometer (Implen, München, Germany). RNA quality and quantity was further determined using a Caliper lab-on-a-chip system (Agilent, Hopkinton, Mass., USA). All isolated total RNA samples had an RNA Integrity Number (RIN) ranging from 8.5-10.0 (highest quality).

DeepSuperSAGE libraries were constructed from total RNA of gonads from 10 pooled individuals by GenXPro GmbH essentially as described previously [Matsumura et al., 2010, 2011; Zawada et al., 2012] with slight modifications. Sequencing was performed on Illumina's Genome Analyzer IIx, and subsequent base calling was carried out by Illumina's GAPipeline v1.0. Distinct libraries were sorted out from the bulked sequencing data according to their respective indices, followed by elimination of PCR-derived tags identified by TrueQuant technology (GenXPro). The 26-bp SuperTAGs were extracted from the remaining sequences, counted, and subsequently combined to TAG groups of common origin (UniTAGs) according to Akmaev and Wang [2004]. For each library, the UniTAG read numbers were normalized to a million sequenced reads in total (tags per million; TPM) to establish comparability. The UniTAG numbers of the pooled datasets from left and right gonads were generated by in silico calculation of the arithmetic mean of the normalized read counts from the left and right gonad in the respective libraries. Logarithmic fold changes (base 2) were determined by pair-wise comparison of the normalized UniTAG numbers in the particular libraries, and statistical significance was assessed by χ² tests according to Man et al. [2000]. If a given UniTAG was only present in 1 of the 2 libraries, the respective TPM count was adjusted from 0 to 0.05 to allow for calculation of fold changes.

The UniTAG reads were annotated in a multi-step procedure using the BLAST software [Altschul et al., 1997] to ensure an unambiguous assignment to their corresponding transcripts and to eliminate any remaining adaptor sequences. Reference datasets were extracted from the publicly accessible databases at the Harvard University (HU; ftp://occams.dfci.harvard.edu/pub/bio/tgi/data/Gallus_gallus) and the National Center for Biotechnology Information (NCBI; ftp://ftp.ncbi.nih.gov/refseq/release/vertebrate_other). A third reference comprising only sequences from Gallus gallus was generated based on the above-mentioned NCBI dataset, and finally, all 3 databases were trimmed to generate subsets (26-bp references) that only comprise sequences from the last possible anchoring site for a given SuperTAG. Beginning with a minimum required BLAST score of 52, UniTAG reads were successively aligned against these reference datasets in the following order of precedence: (1) 26-bp HU dataset, (2) 26-bp NCBI subset of G. gallus, (3) full-length HU, and (4) full-length NCBI subset. Reads which did not attain the specified BLAST score were re-annotated in the same fashion with a lowered required score of 48, and in a third round of 44. Still not annotated reads were finally also aligned against the complete NCBI dataset with a required BLAST score of 42 or above. In this final annotation step, the complete 26-bp NCBI dataset was used as the third reference, followed by the full-length HU and G. gallus NCBI reference and the complete full-length NCBI dataset as final reference. UniTAG reads, which could not be annotated with this multi-step procedure, were excluded from further analysis.

For comparison of the tendencies in gene expression of both datasets (deepSuperSAGE vs. qRT-PCR), log₂ fold changes were calculated with the TPM values of the respective UniTAGs and compared to the qPCR-based fold change (fig. 1). The quantitative real-time PCR assays (TaqMan/dual labeled probe qPCR assays) were performed on transcripts previously selected according to their differential expression profile in the sequencing data. Total RNA isolations of additional tissues were performed as described above. The primer assays for this study were designed and provided by GenXPro GmbH. Quantitative real-time PCR reactions were run as a one-step qPCR (Applied BioSystems StepOne RT-PCR systems) with dual labeled probes (FAM-BHQ) and ROX reference dye (Invitrogen, Darmstadt, Germany) as passive reference. An amount of 5 ng of total RNA served as template for each reaction of 15 µl. The applied One-Step RT-PCR Mastermix (Clontech-Takara QTaq-Mastermix) contained hot start Taq DNA polymerase, an optimized reaction buffer, 5 mM MgCl₂ (final concentration), 2.5-3.5 mM nucleotides (including 200 μM dUTP), and a reverse transcriptase combined with an RNase inhibitor. Specific primers were applied with a final concentration of 0.2-1.0 µM, and the dual labeled probes of 0.016-0.08 µM. The PCR regime consisted of the following steps: (1) reverse transcription: 48°C for 20 min; (2) activation of the hot start Taq DNA polymerase (inhibited by antibody during reverse transcription) at 95°C for 10 min, and (3) 45 cycles with denaturing at 95°C for 15 s, and an annealing/ extension step at 60°C for 1 min.

Amplification of the target genes was monitored by qPCR probe-released fluorescence (FAM dye) at each cycle. The threshold cycle (C_t), defined as the PCR cycle at which a statistically significant increase of reporter fluorescence is first detected (10× above background), was used as a measure for the starting copy numbers of the respective gene. Relative quantification of the amplified targets was performed via the comparative ΔΔC_t method [Pfaffl, 2001]. The amount of target, normalized to an endogenous reference (AURKAIP1) and relative to the passive reference, is given by 2^-ΔΔC_t.

Though deepSuperSAGE and other modern transcriptome analysis technologies nowadays provide in-depth descriptions of gene expression, the assignment of a function to differentially expressed genes is not trivial in organisms that are not as comprehensively characterized as human and mice or rats. When there is no direct experimental evidence available, the most efficient way of function assigning is based on strict orthology, one of the central concepts of comparative genome analysis. By definition, orthologs are genes or proteins in 2 or more species that share significant similarity and that are thought to have diverged from a common ancestral gene that existed in their last common ancestor [Remm et al., 2001; Li et al., 2003; O'Brien et al., 2004; Chen et al., 2006; Hulsen et al., 2006]. Determination of orthology relations facilitates the knowledge transfer between species and can be used to improve both structural and functional annotation in organisms that are less well annotated. At present, there are several ortholog prediction methods and search tools available [O'Brien et al., 2005; Wright et al., 2005; Li et al., 2006; Hubbard et al., 2007]. However, the number of proteins from one species that is considered to be part of the same orthologous group in another species varies from one method to another due to differently employed algorithms and the diversity of species included in these methods [Hulsen et al., 2006]. Since most gene ontology (GO) annotations for newly sequenced species provided by the EBIGOA Project [Camon et al., 2004] are based on Electric Annotation (IEA) that is not assigned by a curator, these annotations do not necessarily include the gene products determined by deepSuperSAGE. Here, we used IEA-based GO (http://www.geneontology.org/) to categorize differentially expressed genes and to assign a potential function to them. The final data set containing information from all the libraries was subjected to GO enrichment analysis using the GenXPro in- house pipeline system.

The original SuperTAG sequences from this project will be available at Genome Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo).

For the description of gene expression profiles in the gonads, we created 4 deepSuperSAGE datasets with about 2 million 26-bp SuperTAGs each, based on the sequences from distinct pools of 10 individual right and left ovaries and testes, respectively. Table 1 summarizes all tags and transcripts in the 4 libraries used for analysis.

From all annotated SuperTAGs, a total of 67,848 tags revealed high homology hits (BLAST score = 52) in 1 or more of the respective databases; 3,339 annotated tags were considered to be predicted by the respective databases and therefore excluded from further analysis.

The reliability and reproducibility of deepSuperSAGE data was tested via qRT-PCR analysis of a subset of 21 transcripts with similar and dissimilar expression between male and female gonads. In most cases, SuperTAG counts from the pooled gonads represented the mean of the ΔΔC_t values obtained for a given RNA from individual specimens. The overall similarity and the high correlation (R² = 0.77) between the tested individual qRT-PCR and pooled deepSuperSAGE expression ratios (fig. 1) demonstrates the exact quantification of transcript levels by deepSuperSAGE profiling.

The cumulated and normalized UniTAGs could be classified into 61,312 G. gallus-specific UniTAGs representing 48,079 UniTAGs with low abundance (0.13-10 tags per million), 6,075 UniTAGs of intermediate abundance (10-100 tags per million), and 1,216 highly abundant UniTAGs (>100 tags per million). These abundance categories contain on average 78, 10, and 2% of the total number of UniTAGs, respectively (upon request, data can be provided by the author). A total of 836 UniTAGs was found to be differentially expressed (M-F, p < 10^-5) between male and female gonads, and a considerable number of UniTAGs were exclusively detected in either the male or female gonads; 114 significantly differentially expressed UniTAGs (p < 10^-5) between the sexes showed more than a 2-fold difference in gene expression. Of these, 40 were up- and 74 down-regulated in testes compared to ovaries. The UniTags found laterally differently expressed are shown in table 2. The 20 most significantly regulated mRNAs between male and female gonads are summarized in table 3 which is an extract of table 2.

The list of transcripts up-regulated in male gonads mainly consists of ribosomal proteins and subunits, inhibin, CDH1-D, DMRT1, and SPARC as the most up-regulated mRNA in testes as compared to ovaries. SPARC (secreted protein acidic and rich in cysteine, originally described as osteonectin) is a highly conserved and abundant component of the extracellular matrix in bone tissues in vertebrate species [Termine et al., 1981; Bassuk et al., 1993]. SPARC binds to collagen type I-V in a calcium-dependent manner and is prominent in regions of tissue morphogenesis which correlates with its proposed function as a mediator of tissue remodeling [Sage et al., 1989; Wu et al., 1996; Delostrinos et al., 2006]. CDH1-D and its homologs are suggested to play a role in temporal and spatial regulation of the anaphase-promoting complex in and outside of the cell cycle in chicken embryos of stage E15-20 [Wan and Kirschner, 2001].

The female up-regulated transcripts comprise phosphatidylcholine 2-acylhydrolase precursor (PLA2G2E precursor), WPKCI, tetraspanin-8 (TSPAN8), gallinacin (GAL-9, GAL-10), and gastrotropin.

Further differentially expressed transcripts between sexes include galectin-2 (LGALS2), galectin-3TM1 (LGALS3-TM1), thrombomucin and its precursor, aromatase, GATA binding factor 2 and 5 (GATA2, GATA5), hRSPO1, tetraspanin-1 (TSPAN1), prolactin receptor precursor, truncated testis-specific box 1-B or box 1-less prolactin receptor, SOX9, AMH and its precursor (fig. 2b, c).

General right-left comparison revealed 1,056 and 822 right-left differentially (Fch > |1|; p < 10^-5) expressed UniTAGs for male and female, respectively. Interestingly, 72 of these UniTAGs are right-left specific in both sexes, with relatively low differences or significance between the sexes. The ovaries display a comprehensive side-specific expression of genes, while differential expression in right and left testes is limited to genes that are found to be generally right-left differentially expressed in both sexes (fig. 2a; table 2, 3). Of these, chromogranin A, specifically androgen-regulated isoform 1, androgen receptor-associated co, somatostatin precursor, followed by surfactant protein A precursor, and cytidine deaminase display the strongest fold changes, and might represent some of the key player genes involved not only in sexual differentiation but rather in lateral development in birds (fig. 2a).

At day E19, differential gene expression in left and right ovaries clearly reflects the tendency of female gonads for growth or regressing in size. Growth-promoting expression patterns comprise the transcripts of e.g. TSPAN1, TSPAN8, GIIEp, galectin-2, gastrotropin, and SOX9. Cytidine deaminase and carboxypeptidase O are almost exclusively expressed in the growing left ovary. The most up-regulated mRNAs in the regressed right ovary (compared to female left ovary and male testes) are thrombomucin/thrombomucin precursor and aromatase. These genes are not located on the W or Z chromosome, but since the gonads are already well developed at day E19, this is not to be expected. These genes seem to be involved in the growth of the left or regression of the right female gonad, rather than in the differentiation of new tissues (fig. 2; table 3).

The subfamily 3 of nuclear receptors includes steroid receptors as e.g. estrogen, glucocorticoid, progesterone, and androgen receptors as well as the estrogen receptor-related receptors [Nuclear Receptors Nomenclature Committee, 1999; Germain et al., 2006]. Neither the estrogen, androgen, nor glucocorticoid receptors are differentially expressed between the sexes and between left and right gonads. Interestingly, the progesterone receptor is slightly up-regulated in the female right gonad, while there is no significant difference between male gonads.

The mRNAs encoding steroidogenic enzymes required to convert cholesterol to androgens are already present in the avian embryo before gonadal differentiation [Nomura et al., 1999]. Besides significantly differential expression of P450arom and W-chromosome-based WPKCI transcripts as key players in the differentiation of gonads, we found the chicken orthologs of estrogen-related transcripts minimally differentially expressed between male and females, as e.g. estrogen receptor binding fragment associated gene 9 protein (EBAG9). The EBAG9 gene is conserved in human, chimpanzee, dog, cow, mouse, rat, and zebrafish. For most of the human EBAG9counterparts, roles in ovarian or breast cancers have been described [Ikeda et al., 2000]. Their expression in chicken provides the opportunity to study their role in the development of sex-specific cancers in an easily accessible experimental system.

The classification of transcripts from male and female gonads (M-F) into GO categories showed that 25,771 UniTAGs grouped to ‘molecular function', 25,864 to ‘cellular component', and 25,798 to ‘biological process'. Differences in the abundance of tags in male and female gonads were considered significant at an enriched p < 10^-5. The most regulated GO:terms are listed in table 4. The GO classification shows that the strongest (p < 10^-15) differences occur in larval development(GO:0002164, up-regulated in left gonads of both sexes), nematode larval development(GO:0002119, up-regulated in left gonads of both sexes), and structural constituent of ribosome(GO:0003735, up-regulated in males). Additionally, the GO:terms post-embryonic development (GO:0009791) and embryo development ending in birth or egg hatching (GO:00009792), both up-regulated in left gonads, are slightly lower significantly regulated than the development GO:terms.

Additionally, the GO:term growth (GO:0040007) is highly affected: it is up-regulated in females, especially in the growing left gonad. The good reliability of the presented data is proven by the fact that the affected GO:terms comprise different developmental stages from pre-hatching to post-hatching. This shows the validity of the identified transcription profiles: the embryonic development is nearly finished and the individual is shortly before hatching - 2 days left at E19.

Since mammalian meiotic sex chromosome silencing [Turner, 2007] makes it difficult to investigate sexually antagonistic genes, birds are the model of choice for studying sex- and right-left biased gene expressions. The presented dataset gives the opportunity to analyze not only differential gene expression in females and males but above all lateral differences within a gender.

The fact that structural gonadal differentiation at day E19 is already very advanced in comparison to other surveys, though the left and right ovaries are still of almost similar size compared to mature gonads, has to be taken into account for a correct interpretation of the presented data. The maturation of gonads at E19 mainly involves growth or regression in size and not the development of new tissues. The chicken Z sex chromosome comprises over 700 known protein-coding and at least 45 non-coding genes (www.ensembl.org/Gallus_gallus). Any of these genes could play a role in sex determination and/or downstream of gonadal sex differentiation. The undifferentiated expression of these sex-determining genes at the already advanced stage of gonadal differentiation at E19 indicates that these mechanisms are not active anymore and explains the discrepancies between sex-chromosome-specific gene expression described in earlier developmental stages [McNagny et al., 1997; Carré et al., 2011] and our findings.

Pluripotent markers like PITX2,Nanog,or PouV, found in other surveys especially in earlier stages [Ryan et al., 1998; Lavial et al., 2007; Intarapat and Stern, 2013] did not display an asymmetric expression in male or female gonads at E19. Additionally, neither Sox2 nor ERNI were found to be differentially expressed in the left or right gonads. All of them had been associated with asymmetric expression patterns in the gonads of both sexes by other authors in Hamburger-Hamilton (HH) stages 33 and 35 [Hamburger and Hamilton; 1951]. Especially PouV, ERNI, and Nanog were reported to be less abundant in higher developmental stages [Lavial et al., 2007; Intarapat and Stern, 2013]. Our data confirm these findings for the above-mentioned genes for a speculation: the pluripotent genes lose importance in higher development stages, and other genes have taken the wheel, reacting to the status the pluripotent markers have predetermined.

In contrast to this, the germ cell marker Cvh as well as TP53 are both slightly up-regulated in the male left gonad. Cvh-expressing cells are present in the central zone of the area pellucida in HH stage 10, Cvh-positive cells were found in the hypoblast layer and anteriorly in the germinal crescent due to morphogenetic movement [Tsunekawa et al., 2000]. Cvh plays an essential role in germline formation, and Cvh-positive primordial germ cells have been located in the left gonads of both sexes in HH stage 35 [Intarapat and Stern, 2013]. TP53 belongs to a group of transport proteins carrying specific substances in the blood or across cell membranes. It is involved in cell cycle regulation as a trans-activator that negatively regulates cell division by controlling a set of genes required for this process.

Another mRNA implicated in left-right asymmetry is BMP4 which is slightly up-regulated in the left gonad of both sexes and significantly up-regulated in the female compared to male gonads. BMP4 is a bone morphogenetic protein that is a potent inducer of bone formation. BMPs belong to the TGF-beta family and negatively regulate the structure and function of the limb apical ectodermal ridge [Pizette and Niswander, 1999]. It is conceivable that this asymmetric expression relates to a similar function in gonadal development [Hoshino et al., 2005].

There are several other genes that are strongly up-regulated both in male and female left gonads (fig. 2a) and thus seem to be involved in later stages of side-specific development than the previously reported genes. Overexpression of genes in the right gonads is less prominent. The mRNAs encoding ovoinhibitor (serine peptidase inhibitor, Kazal type 5, SPINK5) and surfactant protein A, for instance, are only slightly up-regulated in male and female right gonads. These genes are not located on either of the sex chromosomes. Neither are the mRNAs almost exclusively expressed in the female left gonad as e.g. the tetraspanins TSPAN1 and 8, PLA2G2E, LGALS2, and LGALS3-TM1 as well as somatostatin and myosin-11 (MYO11). TSPANs act as molecular facilitators and thus are involved in diverse processes. For example, proteins from the tetraspanin subfamily associate with integrin and support cell motility. Also, they have been ascribed a role in the development of the nervous system [Perron and Bixby, 1999]. PLA2G2E (GO:0004623) catalyzes the calcium-dependent hydrolysis of the 2 acyl groups in 3-sn-phosphoglycerides. LGALS2 (GO:0030246) and LGALS3-TM1are interacting selectively and non-covalently with lactose and β-galactose, respectively. These genes are implicated in cell cycle regulation, and LGALS3 is furthermore involved in the mammalian endometrical system with a substantial role in uteral embryo implantation [Yang et al., 2012]. LGALS3 is also used as a tumor marker in mammals [Chiu et al., 2010; Righi et al., 2010]. Somatostatin acts as an inhibitor of growth hormone secretion [Florio et al., 1994]. Its functions include the modulation of neurotransmission, cell secretion, and cell proliferation [Patel, 1999]. None of the above mentioned genes - as many others differentially expressed between the sexes - is located on the W or Z chromosome which is to be expected in the advanced developmental stage of E19. Whether these differentially expressed genes have been (indirectly) activated by pluripotent markers or not warrants further research.

We found a significantly higher expression of the W-linked WPKCIgene (also known as HINTW or ASW) in embryonic ovaries than in testes. It is already expressed at E4.5 in the developing female gonads and has been suggested to be involved in avian sex determination [Hori et al., 2000; O'Neill et al., 2000; Smith and Sinclair, 2004]. Early WPKCIexpression (E5) has also been detected in the avian spinal cord, spinal ganglion, and myotomes [Hori et al., 2000; Scholz et al., 2006]. The conserved vertebrate embryo gene DMRT1 is one of the mostly favored Z-linked candidate genes controlling gonadal sex differentiation [Nanda et al., 1999; Raymond et al., 1999; Smith et al., 2009; Chue and Smith, 2011], and our data confirms former observations [Shimada, 2002; Alam et al., 2008; Koba et al., 2008] that it is more expressed in testes versus ovaries. This is a further confirmation of the low effectiveness of dosage compensation of birds in higher developmental stages found in other surveys [Ellegren et al., 2007; Itoh et al., 2007]. The transcription profile generated via deepSuperSAGE allows for discrimination between DMRT1 and its isoform e. Interestingly, expression levels in the female left and right gonad differ severely, although both are generally more abundant in testes (fig. 2c). DMRT1 is about 1-fold up-regulated in the female right gonad compared to the left one, while its isoform e is more abundant in the female left gonad. The isoform e lacks the complete 3′ untranslated region from DMRT1 which comprises several microRNA binding sites. The fact that DMRT1 isoform e is far more abundant in all tissues but lacking these microRNA binding sites leads to the speculation of a post-transcriptional regulation of the DMRT1 mRNA via microRNA-mediated degradation aside from alternative splicing mechanisms.

SOX9 is also found to be up-regulated in males, while RSPO1 as well as aromatase are up-regulated in females, concretizing the findings of Ayers et al. [2013] as well as Chue and Smith [2011]. DMRT1 has been suggested to activate the testicular marker SOX9 indirectly because of the time lag between DMRT1 and SOX9expression (day 4.5 and day 6.5 of incubation, respectively), and the chicken homolog of hemogen (cHEMGN) was shown to link the expression of both [Nakata et al., 2013]. In line with this, we observed an exclusive expression of cHEMGN in male gonads.

The secretion of 17β-estradiol in chicken E8 embryos is higher in the left ovary than in the right [Pedernera et al., 1999], and the higher amount of estrogen produced in the left ovary is suggested to protect the left Müllerian duct from the effects of anti-Müllerian hormone (AMH) [Villalpando et al., 2000; Smith and Sinclair, 2001]. This is responsible for the regression of the right Müllerian duct in earlier stages [Oréal et al., 2002]. In our embryonic gonads, the AMH transcript is up-regulated in both male gonads compared to the females. Since males usually do not develop Müllerian ducts because of the elevated AMH level in early development stages, the up-regulation of AMH in male gonads at E19 seems to have another - or additional - function. EBAG9 is the only mRNA related to estrogen that is slightly down-regulated in the male right gonad. The other estrogen receptors, in general, are not differentially expressed between the sexes. Additionally, we did not find any differential expression of WNT4or β-catenin in our dataset which are both thought to be especially involved in ovary development [Chue and Smith, 2011; Ayers et al., 2013]. The most likely reason for the discrepancy between our and previous findings is that the ovaries are nearly finished in their structural differentiation at E19 compared to days E4.5-12.5 [Smith et al., 2008b].

Table 2 as well as the heat map in figure 3 both illustrate the differential expression of a selection of male and female as well as right-left differentially expressed genes. The male gonads in general do not display big differences in right-left comparisons. On the contrary, the regressing female right gonad exhibits strong differences in gene expression of transcripts involved in almost all metabolic pathways not only compared to left and right testes but also compared to the left ovary. Apart from genes involved in sexual development, also genes coding immune reactivity/defense/stress, blood building, as well as transcription factors and kinase transcripts show strong differences in male-female and, especially, right-left (female) gene expression.

Many strongly up-regulated genes in ovaries encode structural proteins (e.g. elastin, collagen, fibrilles), hemoglobin and - especially in the female right gonad - thrombomucin. This cell surface protein is distantly related to CD34 and besides haemoglobin one of the most used human hematopoietic stem cell markers defining thrombocytes and multipotent hematopoietic progenitors [McNagny et al., 1997]. Surprisingly, our results seem to broaden the findings of McNagny et al. [1997] that this gene is highly expressed in female gonads. McNagny and co-workers stained early embryonic stages (E3-E5) and found that thrombomucin is mainly expressed on the surface of intra- and extraembryonic hematopoietic cells. Additionally, immune globulins, cytidine deaminase followed by carboxypeptidase O and defensins (Gal-10), as well as fibrinogen are also up-regulated in the ovaries. Our findings lead to the speculation that, within the scope of evolution, the rudimental right ovary might have gained new function(s) like e.g. hematopoiesis, because mucins have been shown to play both positive and negative roles in adhesion processes and thrombomucin may prevent the differentiation of hematopoietic cells [McNagny et al., 1997].

By application of deepSuperSAGE as a potent technique for the identification of differential gene expression in chicken E19 gonads, we identified several genes that seem to be more important for sexual and lateral development than previously thought. Some of these might be specific for the highly developed gonads in E19 and consequently may exert important functions in developed gonads. The presented dataset is an outstanding possibility for research in the field of right-left development, especially in view of gynandromorph animals [Zhao et al., 2010; Chue and Smith, 2011], and is going to assist studies on gene expression, especially in the context of comparative genomics. Our results will facilitate the study of evolutionarily important traits at the molecular level, deepening our knowledge of avian genomics, avian proteomics, endocrinology, as well as developmental biology, including sex differentiation.

The authors thank Ruth Jüngling and Anne Plöttner, University of Frankfurt, for technical advice, and Sabrina Giebner, University of Frankfurt, for assistance in preparation of the tissues as well as Nico Krezdorn, GenXPro, for handling of big datasets. We additionally thank 2 anonymous reviewers for very valuable advices.

This research was financially supported in the framework of Hessen ModellProjekte (HA project no. 155/08-11), financed with funds of LOEWE - Landes-Offensive zur Entwicklung Wissenschaftlich-ökonomischer Exzellenz, Förderlinie 3: KMU-Verbundvorhaben (State Offensive for the Development of Scientific and Economic Excellence).