Although mammals have apparently symmetric gonads, most avian species including ducks [1
] and chickens [2
] exhibit an unusual asymmetry in gonadal development. In female birds (genetically ZW), only the left gonad becomes a functional ovary, whereas the right ovary completes regression at the adult stage. In male birds (genetically ZZ), testicular development occurs bilaterally [3
In chick embryos, gonadogenesis begins approximately on day 3 of development during Hamburger Hamilton (HH) stage 23 (E3) [5
]. At this stage (the “indifferent stage”), there is no detectable morphological asymmetry between the left and right embryonic gonads of either sex. However, a greater proportion of the circulating primordial germ cells (PGCs) colonize the left than the right gonad [6
]. After sexual differentiation, the morphological appearance of the gonads is initially similar in males and females (HH28 E5.5), but by HH36 (E10.5), the gonads in males and females are distinctly different. The morphological differences between the left and right gonads in the female become more pronounced, while the asymmetry between male gonads diminishes. By HH44 (E18.5), only a few germ cells are randomly found throughout the core region in the right ovary, whereas the germ cells are distributed predominately in the cortex in the left ovary [8
In addition to the morphological asymmetry between the left and right gonads, several genes are related to the asymmetrical development and the right ovary degeneration [9
]. For example, Bmp7
is expressed asymmetrically in the indifferent ridges of both sexes [10
], whereas the estrogen receptor (ER) is expressed in the left but not the right cortex of both sexes [12
has been shown to control asymmetric gonadal development in both sexes of the chick and can rescue the degenerative fate of the right ovary [11
The aforementioned studies focused on early embryonic events mainly based on histological descriptions or on mRNA expression profiles of a distinct gene subset. However, the comprehensive molecular descriptions of this differential expression are far from being completely understood. Many questions have still been left unanswered; for example, is there a global gene expressional regression in the development of the right ovary along with the morphological regression process? What are the molecular factors underlying side-specific development in female chicken gonads? With the advent of next-generation sequencing (NGS), RNA-seq has been used to assess chicken gonadal sex- and side-biased genes at different developmental stages [13
Here, we provide a comprehensive description of the side-specific differential sexual development of female chicken gonads at the level of gene expression on E6 (the onset of morphological differentiation), E12, and post-hatching day 1 (D1, the end of embryonic development), and we discuss the molecular basis of this unusual asymmetry. There was no global regression in the development of right gonad from E6 to D1. Several new genes, gene ontology (GO) terms, and pathways were identified as important for the left-right asymmetry of the gonad. A total of 111 genes, five GO terms, and three pathways were significantly differences between the left and right gonads among all the development stages. We also present the number and the percentage of genes within eight development-dependent expression patterns of differently expressed genes (DEGs) in the left and right gonads of female chicken.
The chicken embryo represents a suitable model for studying sex determination and gonadal asymmetry. Many sex-biased genes have been identified, such as SOX9
, and Lhx9
]. However, only a small number of genes have been identified as expressed in a side-biased manner. Here, we provide the gene expression profiles of female chicken gonads and identify novel candidate side-biased genes by RNA-seq. In addition, a detailed view of the female chicken gonads transcriptome has been revealed by comparing three development stages.
After data filtering, 16.6 to 18.0 M reads per sample were acquired. Based on the results of saturation curves analysis, the gene quantification was reliable for genes with medium or high expression levels. However, some of rare mRNA will be missed.
As gonadal development proceeds, the morphological differences between the female left and right gonads become more pronounced [8
]. However, our results showed that the number of DEGs between the left and right gonads increased from E6 to E12, and then decreased from E12 to D1. We found that the major DEGs were located on the autosome in the third development stage, which is consistent with a previous report [13
]. There were no obvious differences in the percentage of DEGs annotated to the autosome in all development stages. The hierarchical clustering analysis and PCA revealed that E6FL and E6FR were clustered together, and E12FR and D1FR were more similar to the gonad at E6. These results show that there is not a global regression of gene expression profiles in the development of the right gonad.
In addition to morphological asymmetry, asymmetric germ cell distribution was also reported. Previous studies have suggested that the number of PGCs and germ cells in the left gonad was greater than that in the right gonad [7
]. Accordingly, we especially focused on specific germ cell markers including VASA (CVH) [21
], PIWIL1 [22
], DAZL [23
], and TDRD [24
]. PIWIL1 plays a central role during gametogenesis by repressing transposable elements and preventing their mobilization, which mediates the repression of transposable elements during meiosis by forming complexes composed of piRNAs and Piwi proteins and governs the methylation and subsequent repression of transposons [25
]. Recent studies identified the TDRD members, including TDRD1, TDRD5, TDRD7, and TDRD9, as participating in the Piwi pathway and/or retrotransposon silencing [26
, encoding germ plasm components, were required for germ cell formation in many metazoan species [24
]. Piwi was co-localized with VASA
mRNA in germ cells [30
]. We found that the expression levels of CVH, PIWIL1, and TDRD5 were significantly higher in the left gonad than the right gonad in female chicken according to both RNA-seq and qRT-PCR. In addition, the RNA-seq results also showed that MAELSTROM
and the other TDRD members, including TDRD1
, and TDRD9
, were expressed asymmetrically among the three development stages. MAELSTROM
encodes a protein that co-localizes with VASA and a RDE1/AGO1 homolog and plays a crucial role in the piRNA pathway [31
Our GO analyses further revealed that the strongest differences occur in the piRNA metabolic process (Piwi-associated RNA metabolic process), pole plasm, chromatoid body, P granule, pi-body, and germ plasm. The chromatoid body, a germ-cell specific RNA-processing center, has been suggested to be the mammalian counterpart of germ plasm [32
mRNA co-localizes with both Drosophila melanogaster
germ plasm and the mouse chromatoid body. Furthermore, TDRD5 is required for retrotransposon silencing, chromatoid body assembly, and spermiogenesis in mice [29
]. The P granule, a small cytoplasmic, non-membranous RNA/protein complex, aggregates in the primordial germ cells of nematodes [33
Pluripotency genes including cPouV
, and ERNI
were also found to be expressed asymmetrically in embryonic gonads [34
]. Scheider et al. reported that pluripotent markers such as cNanog and cPouV did not display asymmetric expression in male or female gonads at higher development stages [13
]. However, in our study, RNA-seq analysis revealed that cPouV
were expressed at higher levels in the left gonad than the right gonad on both E6 and D1. On E12, the expression levels of cPouV
between the left and right gonads were not significantly different. The expression levels of cSox2
between the left and right gonads were not different among all development stages.
The asymmetric distribution of germ cells was also related to cell proliferation. Our RNA-seq and qRT-PCR analyses have shown that G0S2
2 were preferentially expressed in the left gonad among all the development stages. G0S2 is a multifaceted protein involved in proliferation, apoptosis, and metabolism [35
]. PITX2 was reported to regulate gonadal cell proliferation and morphogenesis. Misexpression of PITX2
in the right gonad is sufficient to induce symmetric development of the gonads and rescue the degeneration of the right gonad [11
To evaluate the development-dependent transcriptomic activities in the left and right gonads, we performed a time course differential gene expression analysis by comparing any two adjacent developmental stages. Here, we present the number and percentage of genes within each expression pattern. Unexpectedly, we found that the percentage of genes within the patterns of UM and MD was markedly different between the left and right gonads. Future studies will focus on the molecular mechanism of these differences in more detail.
4. Materials and Methods
4.1. Embryo Incubation and Tissue Collection
Fertilized White Leghorn chicken eggs (Gallus gallus) were obtained from the Experimental Station of China Agricultural University (Beijing, China) and bred at 37.5 °C under a relative humidity of 55–65% (P-008B Biotype, Showa Furanki, Saitama, Japan). All animals received humane care as outlined in the Institutional Guidelines of the Care and Use of Laboratory Animals at China Agricultural University (Permit Number: SKLAB-2015-06-06, Beijing, China).
The embryonic gonads were collected on embryonic day 6 (E6), E12, and post-hatching day 1 (D1). At each stage, the left and right gonads were dissected separately, placed individually in TRIzol reagent (Invitrogen, Carlsbad, CA, USA) and stored at 4 °C. Embryonic blood cells (1 μL) were used for sex determination as previously described [36
]. After sexing, 6–10 gonads, depending on the stage, were pooled from each replicate according to sex, stage, and side. Each set included two biological replicates. Pooling has been documented as an appropriate approach to prepare samples for expression analysis [37
4.2. Library Construction and Sequencing
Total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s instructions and was treated with DNase I. RNA integrity and concentration were assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). The sequencing libraries were generated using the NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (#E7530L, NEB, Ipswich, MA, USA) following manufacturer’s recommendations. Library concentration was measured using the Qubit dsDNA BR assay kit (Life technologies, Waltham, MA, USA). Insert size was assessed using the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA), and qualified insert size was accurately quantified using StepOnePlus™ Real-Time PCR System (Thermo Fisher Scientific, Waltham, MA, USA). The libraries were sequenced by Annoroad (Beijing, China) on an Illumina NextSeq 500 platform (Illumina, Santiago, CA, USA), and 50 bp single-read reads were generated. The RNA-seq data from the 12 samples have been submitted to the National Center for Biotechnology Information Sequence Read Archive with accession number SRP081829.
4.3. Data Filtering and Alignment
Raw reads generated by the Illumina NextSeq 500 were filtered to remove low quality reads (reads containing more than 15 % bases with a q
-value ≤ 19), adaptor-containing reads, and reads containing more than 5% ambiguous nucleotides. After pre-processing, clean reads were obtained and aligned to the chicken genome (galGal4) using TopHat v2.0.12 (Baltimore, MD, USA) [38
]. Reads mapped to multiple locations and unmapped reads were excluded from gene expression analysis.
4.4. Gene Expression Analysis
The gene expression level was estimated by the RPKM (reads per kilobase million mapped reads) method using HTSeq v0.6.0 (California Institute of Technology, Pasadena, CA, USA) [39
]. Pearson’s correlation coefficient between biological replicates was calculated using genes expressed in at least one of the samples. Hierarchical clustering was performed using Pearson’s correlation distance. Analysis of differentially expressed genes was calculated by DESeq (v1.16, Boston, MA, USA) with a false discovery rate (FDR) <0.05 and an absolute value of fold change (FC) >2 [40
4.5. Functional Annotation and Enrichment Analysis
The GO (Gene Ontology, http://geneontology.org/
) enrichment of DEGs was implemented using the hypergeometric test, in which the p
-value was calculated and adjusted as a q
-value, and data background was genes in the whole genome. GO terms with q
< 0.05 were considered to be significantly enriched. Pathway analysis of DEGs was performed using the KEGG PATHWAY database (http://www.kegg.jp
). The KEGG enrichment of DEGs was also implemented by the hypergeometric test. KEGG terms with q
< 0.05 were considered to be significantly enriched.
4.6. Quantitative Real-Time PCR (qRT-PCR) Analysis
qRT-PCR analysis was performed by using the LightCycler 480 system and the LightCycler 480 SYBR Green Master kit (Roche, Mannheim, Germany) according to the manufacturer’s instructions. The primers sequences used in this study are listed in Table S8
. The cycling parameters were as follows: 95 °C for 10 min, 40 cycles of 95 °C for 15 s and 60 °C for 1 min, followed by one cycle of 95 °C for 15 s, 60 °C for 15 s and 95 °C for 15 s. A final step was performed to obtain a melting curve for each PCR product to determine the specificity of amplification. All samples were analyzed in triplicate on the same plate. The expression levels of genes were calculated relative to the expression of the GAPDH
using the 2-ΔΔCt
4.7. Analysis of Development-Dependent Gene Expression Patterns
Development-dependent gene expression patterns were analyzed by comparing the union of the DEGs (the RPKMs of the DEGs in all stages were >0) between any two adjacent developmental stages, using the younger developmental stage group as the denominator. Each gene showed at least one significant difference between the developmental stages. By indicating a significant upregulation with “up”, a significant downregulation with “down” and an insignificant difference with “maintain”, there were eight possible patterns, including up-up (UU), up-maintain (UM), up-down (UD), maintain-up (MU), maintain-down (MD), down-up (DU), down-maintain (DM), and down-down (DD).
The results are given as the mean ± standard deviation (SD). Statistically significant differences were computed using Student’s t test with the statistical software SPSS (Version 20.0, IBM Corp., Armonk, NY, USA). Probability (p) values of less than 0.05 were considered to be statistically significant.