Sexually dimorphic gene expression emerges with embryonic genome activation and is dynamic throughout development

Sexually dimorphic gene expression emerges very early in development

Eight-cell embryos provide the first snap-shot in mice after the initiation of embryonic genome activation at the 2 cell stage, when maternal transcripts have been depleted [13,14]. We isolated 12 eight-cell embryos for a discovery set and a further 12 for validation. Importantly, the embryos were produced by natural mating of inbred C57BL/6 J mice, without super-ovulation or in vitro culture, thereby providing a unique opportunity to examine sexually dimorphic expression in vivo. We generated transcriptome profiles for individual embryos using an adapted single-cell RNA-seq protocol (sequencing statistics are given in Additional file 1: Table S1). To our knowledge, this is the largest number of genome-wide transcriptome profiles from individual mouse embryos to date. We detected 13 469 transcripts expressed across all eight cell embryos, which is greater than previous microarrays of blastocysts, but less than other RNA-seq data from eight cell embryos that was sequenced at greater depth [11,14]. The sex of the embryos was determined by plotting expression of the Y chromosome gene, Eif2s3y and the X chromosome gene, Xist, as both are expressed in a sex-specific manner in murine blastocysts [10]. A clear separation of 2 distinct groups consisting of 6 female and 6 male embryos for the discovery set was found (Additional file 2: Figure S1). We then sought to identify genes that were differentially expressed between male and female embryos, revealing 69 genes with a genome-wide corrected p-value < 0.1, (Benjamini-Hochberg correction) and −0.5 < log₂ (male/female) > 0.5 (Table 1). The majority of sex dimorphic differences (78%) originated from the X and Y chromosomes (51/69 and 3/69, respectively; Figure 1A), whilst ~22% (15/69) were autosomal in origin (Figure 1B). All X-linked genes were more expressed in female eight-cell embryos. Sex-dimorphic autosomal genes were equally distributed amongst the sexes.

Substantial variation between individual embryos is apparent (Figures 1A,B, Additional file 2: Figure S2), and may reflect small differences in the timing of embryonic genome activation. However, we validated sex-biased expression of these genes in another 12 eight-cell embryos isolated in an independent experiment (8 female, 4 male) (Figure 1C). Biological replicates showed a strong directional correlation (ρ = 0.95, p-value < 2.2x10⁻¹⁶), with all but one of the 69 genes showing the same sex-based expression bias, even though transcripts show average log(male/female) values small in magnitude. Three non-coding transcripts were included in the eight-cell signature, consistent with a role for non-coding RNAs in early embryo development [15]. Two of these were involved in X inactivation, Xist and B230206F22Rik (also known as Ftx) [16,17]. The third, D7Ertd715e is located immediately 3′ to the Snrpn/Snurf imprinted cluster on chromosome 7, but its function is unknown. We performed KEGG analysis of X-linked genes, using detected X chromosome genes as background and separately for the autosomal genes, using all detected autosomal genes as background. Neither analysis revealed enrichment after Bonferroni correction (p < 0.05).

The dynamics of sexual dimorphism at different developmental stages

Having defined a high confidence, replicated gene expression signature from eight-cell embryos, we explored the temporal dynamics of these genes with respect to sex in pre- and post-implantation development. We obtained publically available transcriptome data derived from multiple stages (late two-, four- and sixteen-cell) of mouse pre-implantation development [14]. Not all embryonic stages profiled by Deng et al., were included because sex could not be determined (earlier than late two-cell), only one sex was sampled, or the sex ratio was highly skewed (e.g. blastocysts). Several factors make using this data for defining stage-specific signatures problematic. The embryos profiled by Deng et al. are F1 hybrids of two genetically diverse inbred mouse strains (CAST/EiJ females × C57BL/6 J males). As such, embryo sex is directly confounded by genotype. They were also produced using super-ovulation, which may influence embryo development [18]. Furthermore, only small numbers of embryos were profiled at each stage (<5 total, Table 2), reducing the confidence of signatures defined de novo. Nonetheless, the data is useful for comparative analyses using the high confidence signature we defined from eight-cell embryos. To ensure that the coverage and dynamic range between these datasets is comparable, as we used different RNA-seq protocols, we first compared our eight-cell male embryos (n = 6) to eight-cell male embryos profiled by Deng et al., (n = 4), establishing a high degree of correlation between these profiles (Additional file 2: Figure S3, ρ =0.84, p-value < 2.2 × 10⁻¹⁶). To get a more expansive view of the developmental dynamics of dimorphic expression, we additionally generated embryonic stem (ES) cells derived from eight-cell embryos (2 male, 2 female) using the 2i method [19]. The 2i methodology results in ES cell lines that are transcriptionally and epigenetically similar to pre-implantation epiblast [20-22]. This was confirmed by a strong correlation of the 2i ESC transcriptome to recently published inner cell mass single-cell profiles (Additional file 2: Figure S4, ρ = 0.70, p-value < 2.2 × 10⁻¹⁶) [23]. To profile post-implantation development we generated RNA-seq profiles from 17.5 dpc (days post-coital) fetal and adult liver (12 male, 12 female at each stage). We selected liver, because sex-based differences are well characterised and prolific in adults [24]. The ES cell, fetal and adult transcriptomes were all from a C57BL/6 J background.

Taking the two-, four- and sixteen-cell embryo data generated by Deng et al., together with our ES cell, fetal and adult liver data, we next examined whether the dimorphic expression of the 69 genes defined in eight-cell embryos is conserved throughout development. Unfortunately, Deng et al., only sampled males at the eight-cell embryo stage, precluding a direct comparison at this developmental time-point. Of the 69 signature genes, only one gene located on the Y chromosome (Eif2s3y) was expressed in a sexually dimorphic manner (logFC > 0.5) across all profiles, from 2 cell embryos onwards. From the four-cell stage onwards, a substantial proportion of transcripts are of embryonic origin [13,14] and the sexually dimorphic expression of the other Y chromosome genes (Ddx3y, Ube1y1) is established. The sex-specific expression of the majority of X chromosome genes is established from the two- to four-cell stage and maintained throughout pre-implantation development and in the ES cells (Figure 2A). However, the majority of the sex-chromosome genes within this signature do not show dimorphic expression in fetal or adult liver.

Interestingly, the 15 autosomal genes that are sexually dimorphic in eight-cell embryos did not show consistent directional changes at other time-points in pre-implantation development. As these genes were expressed at comparable levels to the sex-chromosome encoded transcripts, it seems unlikely that an increased signal to noise ratio is responsible for this. Despite not showing dimorphism at other pre-implantation stages, 14/15 of these genes show biological validation at the eight-cell stage (Figure 1C), suggesting that this stage-specific autosomal sexual dimorphism is bona fide. The expression level of these genes seems to be very dynamic throughout pre-implantation development (Additional file 2: Figure S5). Similar pre-implantation stage-specific expression has also been shown for genes identified as dimorphic in bovine blastocysts [11]. Indeed, global autosomal gene expression reaches a nadir at the eight-cell stage (Additional file 2: Figure S6). Collectively, we show that sex-chromosome genes show largely consistent sexual dimorphism throughout pre-implantation development and in ground state pluripotent stem cells, regardless of genetic background, but that dimorphism of autosomal genes is more temporally dynamic.

X-linked dimorphic genes escape paternal X inactivation in pre-implantation development

The emergence of sexually dimorphic expression of Xist at the four-cell stage coincides with the initiation of transcriptional inactivation of the paternal X chromosome in mice [25]. Higher X-linked expression in females must therefore reflect either an up-regulation of these genes from the maternal X chromosome in females, or escape from paternal X chromosome inactivation in the pre-implantation stages. Consistent with the escape from silencing of the paternal X we find that the majority of X-linked sex dimorphic genes are located distally from the X inactivation centre (Additional file 2: Figure S7), confirming previous observations [26]. Indeed, allelic data generated by Deng et al., confirms that 35/38 (~92%) genes that had adequate allele-specific information were expressed from both maternal and paternal X chromosomes at one or more developmental time-points (four-, sixteen- or early blastocyst) [14]. Although there is some suggestion that inter-specific crosses may have altered X inactivation, reassuringly, two of the three transcripts solely of paternal origin are known to be involved in silencing the chromosome from which they are transcribed, Xist and B230206F22Rik (also known as Ftx) [16,17]. The other transcript (Gla) showed paternal-specific expression at the four-cell stage, with allelic data for later stages not available. Notably, silencing of paternal expression of this transcript has previously been shown to be initiated after the eight-cell and completed only by the blastocyst stage [27].

Sex-biased expression of most X chromosome genes that are dimorphic in eight-cell embryos, with the exception of Xist, and Eif2s3x is erased (−0.5 < logFC > 0.5) in fetal and adult liver (Figure 2A), but female-biased expression was retained for 5/51 X-linked genes in publically available gene expression array data from late gestation placenta [28], at log(male/female) < −0.5. As extra-embryonic tissues in the mouse demonstrate imprinted inactivation of the paternal X chromosome, we asked if dimorphic X-linked expression is more similar in general between these two tissues. Using a cut-off of log(male/female) < −0.5, we identified 79 X-linked genes showing female biased expression in placenta and found that there was a trend towards female-biased expression of these genes in the eight-cell embryos (t test, mean = −0.14, p-value = 0.046). Consistent with re-activation of the inactive paternal X chromosome in the epiblast-like 2i ES cells, the female bias in expression of the dimorphic X-linked genes identified from eight-cell embryos is maintained and expanded upon (Additional file 2: Figure S8). In fetal and adult somatic tissues, a different set of X chromosome genes are dimorphic.

Sex-biased expression is correlated with post-translational histone modifications

Specific post-translational histone modifications are associated with particular transcriptional states and genomic features. We sought to explore the relationship of sexually dimorphic expression in pre-implantation development with chromatin features. To address this, we generated genome-wide profiles for three post-translational histone modifications in our 2i ES cells using chromatin immunoprecipitation combined with deep sequencing (ChIP-seq). 12–20 million mapped reads were generated for each mark (Additional file 1: Table S1). Determining average profiles for the 2000 most high- or low-expressed genes confirmed the expected enrichment of H3K4me3 at active gene promoters, whilst H3K27me3 was depleted (Additional file 2: Figure S9). Both H3K27me3 and H3K9me3 were enriched at repetitive elements and H3K9me3 was associated with imprinted genes, confirming the distribution of these modifications is typical of ground state pluripotency [20]. Genes located on the X chromosome conform to genome-wide distribution patterns.

Given the association of H3K4me3 with transcriptional activity, we next asked whether the X-linked genes identified as sex-dimorphic in the eight-cell embryos are more enriched for H3K4me3 in female 2i ES cells. Dimorphic expression of 46/51 of these genes is conserved in the ES cells. Before assessing relative enrichment for X-linked dimorphic genes, we first needed to account for the difference in chromosome dosage between the sexes. When ChIP-seq reads were counted in 3 kb windows surrounding the TSS, we were able to show that all three histone marks show relative enrichment in females, as expected. We developed a model allowing us to correct for chromosome dosage in subsequent analyses (see Methods). Interestingly, after applying the model for dosage correction, we found that sex-dimorphic genes show female-biased enrichment of H3K4me3 surrounding the TSS (FE_female = 1.47; p-value = 0.002), whilst H3K9me3 (FE_female = 0.72; p-value = 0.028) and H3K27me3 (FE_female = 0.63; p-value = 0.005) show greater enrichment in the males (Figure 3). This is consistent with enrichment of H3K4me3 in female compared to male ES cells at the X-linked Rhox6/9 genes which is lost upon differentiation and loss of expression [29].

Although power for de novo calling is limited by the number of biological replicates for the ES cells, we were able to identify sex-biased expression (genome-wide corrected p-value < 0.1, −0.5 < logFC > 0.5) for use in correlative analyses. Whilst the X-linked genes identified as sex-dimorphic in eight-cell embryos show evidence for escape from paternal X inactivation, they represent only a subset of the X-linked genes identified as sex-dimorphic in ES cells, where the paternally silenced X has been reactivated (Additional file 2: Figure S8). Consistent with what was found for the transcripts identified from eight-cell embryos, we observe that transcriptional activity of X-linked genes defined from ES cells is associated with enrichment of H3K4me3 in the 3 kb surrounding the TSS (Figure 4A; t-test p-value < 2.2x10⁻¹⁶ (mean = −0.136)). However, the ES cell defined X-linked dimorphic genes differ from the eight-cell embryo defined subset in that whilst male cells show H3K9me3 enrichment (Figure 4B; t-test p-value < 2.2x10⁻¹⁶ (mean = 0.619)), male specific H3K27me3 enrichment is very weak (Figure 4C; t-test p-value = 0.04 (mean = 0.037)), suggesting differential regulation of these gene subsets at the chromatin level, despite both maternal and paternally inherited alleles in both subsets being transcriptionally active in the ES cells.

A positive correlation between expression and H3K4me3 enrichment in the 3 kb surrounding the TSS of sex-dimorphic autosomal genes expressed in ES cells was found (Figure 4D; ρ = 0.57, p-value < 2.2x10⁻¹⁶). TSS H3K9me3 enrichment does not correlate with gene expression, although is enriched in males (Figure 4E; ρ = 0.0052, p-value = 0.89), whilst H3K27me3 is very weakly anti-correlated (Figure 4F; ρ = −0.077, p-value = 0.036). The autosomal genes identified in the eight-cell embryos were not called as significantly dimorphic in the ES cells, although they were still expressed. Furthermore, as there were so few sexually dimorphic autosomal genes in the eight-cell signature, we could not make correlations with the ESC histone marks. Our findings suggest that sex-biased gene expression associates with enrichment of H3K4me3 at the TSS in the sex that has higher expression. This relationship is true for both X-linked and autosomal genes.

Sexual dimorphism is more pronounced in adult compared to fetal liver

Sexually dimorphic expression of the eight-cell signature genes in liver at adult and fetal stages did not recapitulate that found in pre-implantation development (Figure 2). To investigate the inter-relationship of sex-biased gene expression at later developmental stages, we identified a signature of 1488 dimorphic genes within adult liver (genome-wide corrected P < 0.1, −0.5 < logFC > 0.5; Additional file 3: Table S2). Autosomal genes account for ~96% of the dimorphic genes, in contrast to what is observed at pre-implantation stages (Additional file 2: Figure S8). KEGG analyses identified many enriched pathways, with drug metabolism as most significant, consistent with previous findings (Additional file 4: Table S3) [30]. Reassuringly, our data strongly correlated with previously published adult liver dimorphic differences, despite the different profiling platforms and genetic backgrounds between the two datasets (ρ = 0.82; p-value < 2.2x10⁻¹⁶) [24]. On the premise of this, we established a cut off (−0.5 < logFC > 0.5) by which we could define a signature for the other tissues; adipose, muscle and brain, profiled by Yang et al. Consistent with previous analyses we found that there was not a significant correlation for sex-dimorphic expression across tissues (data not shown), whilst X-linked genes showed modest female-biased expression differences [31]. Furthermore, the signature defined from eight-cell embryos did not show directional consistency in any of the adult tissues (data not shown).

We also defined a sexually dimorphic signature of 394 genes from late gestation fetal liver (genome-wide corrected p-value < 0.1, −0.5 < logFC > 0.5; Additional file 3 Table S2). As with adult liver, the majority (~94%) of sex-biased differences were autosomal in origin (Additional file 2: Figure S8). KEGG analysis identified 8 signaling and disease-associated pathways that were enriched in fetal liver after Bonferroni correction p-value < 0.05; (Additional file 4: Table S3). Using the same criteria, we identified 3.8X as many sex-dimorphic transcripts in adult compared to fetal liver. There was substantial overlap between the two developmental stages, with 72% (264/369) of autosomal, 67% (14/21) of X-linked and 100% (4/4) of Y-linked genes identified as dimorphic in fetal liver also being identified as dimorphic in adult liver (Additional file 3: Table S2), collectively representing a 11.95 fold enrichment (p-value < 0.001), when all expressed genes are considered as background. Yet, this enriched subset of dimorphic transcripts that is common across both fetal and adult liver do not show directional consistency across datasets (Figure 5) (X-squared = 2.07, p-value = 0.15). Our findings suggest that sexually dimorphic gene expression is present in late gestation liver, but is further expanded upon in the equivalent adult tissues, with some sexual-dimorphism being developmental stage specific. Consistent with previous findings, we did not find any up-regulation of X chromosome expression to equate with the level of autosomal gene expression, thereby refuting Ohno’s hypothesis [32]. The expression level of the X chromosome relative to autosomal genes was remarkably similar in males and females (1.84X and 1.82X, respectively for adult liver), but showed some variation according to developmental stage (1.43X and 1.43X, respectively for fetal liver).

Sex-based expression differences in adult liver are known to be largely dependent on differences in growth hormone regulation in response to sex-specific gonadal hormones [33,34]. To test if growth hormone responsiveness underlies the differences in dimorphic expression between fetal and adult liver, we asked if the dimorphic genes present in adult, but not fetal liver were enriched for previously identified growth hormone responsive transcripts [35]. Indeed, adult-specific sex-dimorphic genes had a 5.4X enrichment (p < 0.001). Collectively, our results suggest that sex-biased expression of some genes in liver is intrinsic, but that sex differences in growth hormone induces dimorphism of additional genes. Our data suggests that this latter mechanism is not yet operational in late gestation.

Animal maintenance and tissue isolation

All animal procedures were conducted in accordance with the Home Office Animals (Scientific Procedures) Act 1986 (Project License number - 70/6693). C57BL/6 J mice were obtained from Charles River Laboratories, UK. All animals were maintained on standard laboratory chow and a 12 hr light/dark cycle. Male mice were housed with virgin females overnight. The detection of a vaginal plug the following morning was considered 0.5 dpc. Females were killed at 2.5 dpc and eight-cell embryos isolated from the fallopian tubes by flushing. Alternatively, females were killed at 17.5 dpc and embryonic liver was collected and snap frozen with liquid nitrogen. Adult animals were killed at 15–20 weeks of age and tissues collected and snap frozen. All animals were killed between 10 am and 12 pm.

Single cell isolation from 8 cell embryos

Single embryos were transferred to acidic Tyrode’s solution to remove the zona pellucida, then washed in PBS-BSA (1 mg/mL) and dissociated into single cells as described previously [50]. All cells from a single embryo were used to generate a single RNA-seq library. For the discovery set, embryos of both sexes were derived from 4 independent litters, embryos used for biological replication were derived from an Additional 5 independent litters. The sires of the replicate litters were exposed to in utero protein restriction.

Derivation of embryonic stem (ES) cell lines

ES cells were derived from eight-cell embryos as previously described [51]. Briefly, isolated eight-cell embryos were cultured in KSOM medium supplemented with (2i): mitogen-activated protein kinase inhibitor (PD0325901, 1 μM) and glycogen synthase kinase-3 inhibitor (CHIR99021, 3 μM). After two days, embryos were transferred to NDiff227 medium with 2i and LIF and allowed to develop into blastocysts. The trophectoderm was lysed by immunosurgery, and the ICM (inner cell mass) of each embryo plated in NDiff227 supplemented with 2i and LIF and for expansion to generate 2i ES cell lines.

Generation and sequencing of RNA-seq libraries

The eight-cell embryo RNA-seq libraries were generated using a single-cell protocol adapted from [50]. Further details are provided in Supplementary Information.

Generation and sequencing of ChIP-seq libraries

The chromatin immunoprecipitation (ChIP) assay was performed according to previously published protocols with minor modifications [52]. Chromatin was sonicated to get fragments of 100 to 700 bp and immunoprecipitated with the following antibodies anti-H3K27me3 (07–449, Millipore), anti-H3K4me3 (39159, Active Motif) and anti-H3K9me3 (ab8898, Abcam). ChIP-seq libraries were prepared using the Illumina ChIP-seq librarary prep kit, according to the manufacturers’ instructions.

Analysis

All analysis was performed on UCSC reference genome mm10 and the gene annotation file downloaded from Tophat website (downloaded on September 25^th 2013).

RNA-Seq analysis

All RNA-Seq data was mapped using Tophat v2.0.4 with Bowtie 2 v2.1.0 and samtools v0.1.18 using default settings. The mean insert sizes and standard deviations were calculated in silico using Picard Tools v1.98 from 1,000,000 reads. Duplicates were filtered using Picard Tools and reads were assigned to a gene using HTSeq. Differential analysis was performed using DESeq2 v1.4.5 which uses a generalised linear model in which counts are modelled using a negative binomial distribution. Genome wide corrected p-values were calculated using Benjamini-Hochberg multiple testing adjustment procedure. Normalised variance stabilizing transformed counts were used for all plots and further analysis. Existing public data was extracted as raw counts and analysed with DESeq2 in the same manner as our data.

ChiP-Seq analysis

All ChiP-Seq data was mapped using Bowtie 2 v2.1.0 with default settings. Coverage across the genome was calculated using genomeCoverageBed from Bedtools v2.17.0 and converted to BigWig format using bedGraphToBigWig downloaded from the UCSC website. For correlation with RNA-Seq data the average coverage over +/−1.5 kbp of the TSS of each gene was calculated using bigWigAverageOverBed.

Chromosome dosage model

To account for the difference in chromosome dosage between the sexes we defined a model based on the distribution of reads for the 2 male samples. We first removed the Y chromosome and then counted the (average) total number of reads (across the 2 samples) within a 3 kb window around the TSS of each of the genes across all remaining chromosomes including chrX (N_total). We then calculated the average total number of reads on the X chromosome (N_x). A scaled read density for chromosome X genes was then calculated:

Where N_i represents the number of reads for gene i on chromosome X. We can then compare this to the actual (average) read density for the X chromosome of the female samples calculated as:

where N_ftotal is the total number of reads (excluding the y-chromosome) for the female samples and N_i represents the number of reads for gene i on chromosome X.

Data access

ChIP-seq data and RNA-seq data are available in the NCBI Gene Expression Omnibus (GSE59222), and will be made public upon acceptance of the manuscript for publication.