Rewirable gene regulatory networks in the preimplantation embryonic development of three mammalian species

Abstract

Mammalian preimplantation embryonic development (PED) is thought to be governed by highly conserved processes. While it had been suggested that some plasticity of conserved signaling networks exists among different mammalian species, it was not known to what extent modulation of the genomes and the regulatory proteins could "rewire" the gene regulatory networks (GRN) that control PED. We therefore generated global transcriptional profiles from three mammalian species (human, mouse, and bovine) at representative stages of PED, including: zygote, two-cell, four-cell, eight-cell, 16-cell, morula and blastocyst. Coexpression network analysis suggested that 40.2% orthologous gene triplets exhibited different expression patterns among these species. Combining the expression data with genomic sequences and the ChIP-seq data of 16 transcription regulators, we observed two classes of genomic changes that contributed to interspecies expression difference, including single nucleotide mutations leading to turnover of transcription factor binding sites, and insertion of cis-regulatory modules (CRMs) by transposons. About 10% of transposons are estimated to carry CRMs, which may drive species-specific gene expression. The two classes of genomic changes act in concert to drive mouse-specific expression of MTF2, which links POU5F1/NANOG to NOTCH signaling. We reconstructed the transition of the GRN structures as a function of time during PED. A comparison of the GRN transition processes among the three species suggested that in the bovine system, POU5F1's interacting partner SOX2 may be replaced by HMGB1 (a TF sharing the same DNA binding domain with SOX2), resulting in rewiring of GRN by a trans change.

Figure 1.

Summary of preimplantation embryonic development. (A) Photomicrographs of human embryos from representative stages of preimplantation development. (B) Table of all the transcriptional arrays performed with pertinent collection parameters. The average Pearson correlations are listed for every experimental group. (C) Hierarchical clustering of gene expression in human, mouse, and bovine embryos during representative stages of PED. This clustering is based on 8479 (human), 7093 (mouse), and 9474 (bovine) informative probe sets. Adjacent stages tend to exhibit similar transcriptional profiles, except for those stages near the time of zygote genome activation. Later PED stages, characterized by expression of genes important for cellular differentiation, segregate into clusters distinct from earlier stages in all three species.

Figure 2.

Cross-species comparison of gene expression. (A) Pictorial description of trends in transcription of a subset of genes with distinct roles during PED. Human, mouse, and bovine probe sets are depicted by red, green, and blue curves, respectively. When a gene is targeted by multiple probe sets, all probe sets are shown. The maternally deposited transcript for TCF7 is degraded at two-, four-, and eight-cell stages in mouse, human, and bovine embryos, respectively. SIN3A is highly abundant in zygotes, and exhibits consistent decreasing expression patterns in all three species. POU5F1 exhibits a consistent up-regulation in all three species, however its expression peaked at the eight-cell stage in mouse and at the morula stage in human and bovine. MTF2 transcript exhibits a reversed expression pattern in mouse and in the other two species, as reproduced by two probe sets in each species. (B) Conserved and species-specific MT and ZGAT. A percentage of MT, 33.7%, is shared in all three species, whereas only 3.2% of the ZGAT is shared in all three species. MT has a 10.7-fold larger chance of being shared across species than ZGAT (P-value < 10⁻³⁰).

Figure 3.

CRM-containing transposons. (A) Proportions of mouse ChIP-seq detected binding regions that overlap with transposons (blue) and a specific class of transposons (red). An asterisk (*) denotes a specific class of transposons that is significantly overrepresented in the ChIP-seq detected regions. (B) Numbers of CRM-containing transposons detected by ChIP-seq. (C) The percentage of transposon-flanked genes with mouse-specific expression (solid lines) and the overall percentage of genes with mouse-specific expression (dashed lines).

Figure 4.

Murine-specific expression of the transcription factor Csda is induced by a murine-specific transposome carrying POU5F1 binding site. (A) CSDA is expressed during mouse PED, but not in human and bovine PED. CSDA is analyzed by two mouse probe sets (red), four human probe sets (green), and one bovine probe set (blue). (B) ChIP-seq data in mouse ES cells. The region of the mouse genome upstream of the Csda gene contains a functional POU5F1 binding site, which would appear to be competent to bind the POU5F1 protein. Two CTCF-associated TFBSs attract the CTCF insulator and define the boundaries of a genomic regulatory region (yellow) between them. The functional POU5F1 TFBS is carried by a murine-specific transposon Lx8, which is conserved in the mouse and the rat, but not in other mammals. (C) Expression and genome sequence data jointly suggest that the murine lineage gained a regulatory relationship between Pou5f1 and Csda.

Figure 5.

Mouse-specific expression of Mtf2 is due to the creation of transcription factor binding sites produced through the combined effects of a single nucleotide mutation and transposon-introduced sequences. (A) Binding sites for POU5F1, SOX2, NANOG, and RNA Polymerase II lie in close proximity to the mouse Mtf2 gene, between two CTCF binding sites (yellow region). (B) A single nucleotide mutation generates a NANOG binding site in the mouse, but not in human or cow. ChIP-seq data accumulated over 62 sequence reads, overlapping on a single nucleotide in the mouse genome. (C) Electrophoretic mobility shift assay results confirm that mutation of the mouse binding site back to the version in other mammals (A-to-T) eliminates NANOG binding in vitro. (A>T) A-to-T mutation; (A>C) A-to-C mutation. (D) A mouse-specific transposable element inserts a CRM for POU5F1 and NANOG into the region upstream of Mtf2: This CRM is bound by POU5F1 and NANOG in mouse ES cells. (E) Expression levels of Mtf2, Pou5f1, and Nanog in mouse ES and trophoblastic stem (TS) cells. (F) Mtf2 knockdown in mouse ES cells affects the expression of other transcription factors and signaling proteins, including the down-regulation of another ICM-specific gene, Notch3. pSUPER and Luci are control vectors that carry an empty vector (pSUPER) and shRNA against a luciferase gene (Luci). RNAi-1 and RNAi-2 are knockdowns mediated by two different shRNAs.

Figure 6.

Comparison of the state transition of GRN across species. For each time point in each species, the active state of a GRN module is drawn. The nodes represent the genes with detectable transcripts; the edges represent the predicted interactions of the connecting nodes. Interactions were predicted by the MATISSE program, using correlation of expression changes in a window of three time points centered at the current time point, and protein–protein interaction data.