BACKGROUND
hypothesis testing is the foundation of the scientific methodology  <cit> .
in biomedical research, focused and hypothesis-driven projects are especially encouraged due to the enormous complexity of the field.
but little attention has been paid to how such hypotheses are generated.
beyond intuition, a quantitative, evidence-based methodology for generating research hypotheses will increase the tempo of biomedical research.
the idea of leveraging the “big data” of tens of thousands of large genome-wide datasets has been discussed  <cit> .
such data-driven, inductive methods can complement hypothesis-driven approaches to form an iterative process of ongoing research  <cit> .
our knowledge about many fundamental biological processes remains limited and fragmented.
we can take advantage of the massive genome-wide datasets to learn about biological systems, akin to the study of an unknown planet in space exploration.
first proposed by tukey, the goal of exploratory data analysis   <cit>  is to explore the data and find what it can tell us.
it is open-ended and much broader than testing pre-defined hypotheses or building desired models.
eda is a philosophically different approach to statistical analysis, not a new set of tools.
in addition to data visualization, all statistical techniques can be used to ask various questions  <cit> .
eda can assess distributions, structures and dependencies that are useful for modeling.
more importantly, eda can reveal important characteristics and trends, which can help formulate new hypotheses for further investigation  <cit> .
indeed, many big scientific discoveries are made accidentally  <cit> .
open-minded exploration of scientific data should be an essential part of a research project.
eda can serve as a strategy to generate new hypotheses from biomedical data.
the infrastructure is largely in place due to many open, collective efforts for sharing of data, annotation and software.
the main challenge is to organically combine multiple datasets to gain a holistic understanding of complex biological systems.
this requires interpretation of results within biological contexts and the ability to survey the literature while analyzing massive data.
such broad, data-driven effort aimed at general understanding of a biological process can be referred to as exploratory bioinformatics investigation .
this study represents an attempt of ebi on gene regulation in pre-implantation mouse embryos.
in the initial stage of mammalian development, zygotes undergo maternal-to-zygotic transition   <cit> , during which maternal factors are eliminated while mrna and protein synthesis using zygotic genome are initiated to take control of embryo development.
a highly coordinated cascade of regulatory mechanisms unfolds rapidly to give rise to several cell lineages and the formation of blastocysts  <cit> .
understanding the complex process of pre-implantation development  is important for both fertility related interventions as well as manipulation of embryonic stem cells .
many gene expression studies of the early mammalian embryos has been carried out using expressed sequence tags  , dna microarrays , rna-sequencing   <cit> , and, more recently, single-cell rna-seq .
these studies documented the dynamic waves of gene expression at different stages.
using the powerful single-cell rna-seq technique  <cit> , deng et al.
analyzed hundreds of cells from embryos of mixed background mice and found strong evidence for random and widespread monoallelic expression  <cit> .
with sequence-level detail at single-cell resolution, this and other similar datasets  <cit>  can be used for in-depth study of gene regulation during pd.
the epigenetic remodeling of maternal and paternal genomes was also revealed by dna methylation profiling .
in zebra fish, transcription factors  such as pou5f <dig> , nanog , and soxb <dig>  were found to activate zygotic gene expression  <cit> .
but the molecular mechanisms of pd remain poorly understood.
retrotransposon expression is a defining event in genome reprogramming during pd  <cit> .
transposable elements  cover 30–50% of mammalian genomes  <cit> .
some are actively transcribed  <cit> , even retrotransposed, as much of the genome is briefly hypomethylated in early embryos.
the expression of retrotransposons is dynamic and stage-specific .
evsikov et al.
showed that retrotransposons, especially long terminal repeats , are abundantly represented in transcripts derived from mouse embryo at 2-cell stage  <cit> .
class iii ltrs such as merv-l family ltrs are transcribed at extremely high levels, accounting for about 3% of total transcriptional output at this stage  <cit> .
human endogenous retroviruses  hervk ltr are expressed at zygotic genome activation  at eight-cell stage  <cit> .
long interspersed elements , another major type of retrotransposon, are also expressed during pd  <cit> , similar to short interspersed nuclear elements .
interestingly, some endogenous retroviral activities have been found to be associated with and can serve as markers of pluripotency and totipotency  <cit>  in stem cell populations, underlying the importance of studying retrotransposon expression.
transcription of retrotransposons can directly influence the expression of neighboring genes.
by analyzing est libraries, peaston et al.
analyzed the expression of different types of tes during pd and found that some tes, especially merv-l family ltrs, provide alternative 5′ first exons to  <dig> chimeric transcripts in 2-cell mouse embryos  <cit> .
besides ltrs, lines and sines also contributed a small number of chimeric transcripts.
lines have also been found to initiate fusion transcripts in other studies  <cit> .
more broadly, by analyzing cap-selected 5′ end of mouse and human transcripts from various embryonic and adult tissues, faulkner et al.
<cit>  found that 6–30% of all transcripts initiate from tes.
these transcripts are often tissue-specific.
studies on long noncoding rnas  also identified ~ <dig>  tes critical for the biogenesis of about 30% of total lncrnas sequences  <cit> .
tes can also influence gene expression indirectly by shaping the epigenetics landscape , as gc-rich sines tend to be found near genes while lines has the opposite distribution.
some transcripts from tes give rise to small rnas that involve in post-transcriptional regulation during pd .
line- <dig> rna can regulate its own expression  <cit> .
chip-seq data in human cells shows that tes contribute 25% of binding sites for tfs critical for embryonic stem cells such as pou5f <dig>  nanog and ctcf  <cit> .
xie et al.
<cit>  compared expression dynamics of pd in human, mouse and bovine and found substantial difference in co-expression network, which could be partly due to the cis-regulatory modules provided by species-specific transposons.
recently, tӧhӧnen et al.
used single-cell rna-seq to analyze  <dig> single-cells from early human development  <cit> .
they found that the promoters of  <dig> genes upregulated early at 4-cell stage are enriched with alu elements that harbor motifs, including a “taatcc” core motif bound by pitx and otx homeobox family  <cit> .
these studies provide substantial evidence for tes’ broad role in regulating gene expression during pd.
in this study, we use ebi approach to study the dynamic expression of both regular genes and tes in mouse pd.
starting from the large single-cell rna-seq data of deng et al.
<cit> , our approach is to systematically observe the gene expression dynamics to help develop a broad understanding of the regulatory mechanisms of this complex process.
our goal is to produce insights that can lead to novel, testable hypotheses on early embryo development.
the data was analyzed alongside other expression and epigenetic studies of pd, using various tools and annotation databases .
our analyses provide evidence for co-regulation of transposons and regular genes in early embryo development in mouse, followed up with similar observation in human, bovine and zebrafish.
this adds to existing evidence that the expansion of species-specific tes may help rewire developmental pathways  <cit> .
motif analysis of promoter sequences of stage-specific genes identified many homeobox domain tfs as potential regulators of pd, which could be experimentally tested.
we also examined the non-random distribution of tes in the mouse genome and the enrichment of sines in the promoters of genes specifically expressed in escs.
finally, we discussed the evolutionary benefits of transposons in promoting genetic diversity, especially in slow-reproducing animals.fig.
<dig>
a exploratory bioinformatics investigation on gene regulation in mouse pre-implantation development.
b hierarchical clustering of gene expression during pd.
each of the  <dig>  rows represents a gene.
columns correspond to samples labeled by developmental stages .
red indicates expression levels higher than average for the row.
expression lower than average is shown in green
RESULTS
to study allelic-specific transcription, deng et al.
<cit>  analyzed individual cells from mixed lineage mouse embryos  using the single-cell rna-seq technique  <cit> .
excluding technical control samples,  <dig> rna-seq libraries were generated, representing cells from zygote to blastocyst stages.
the 2-cell  stage is divided into early, middle  and late phases.
each library contains about  <dig> million reads, mostly  <dig> base pair  long.
with a total of about  <dig> billion reads, this massive data enables in-depth ebi on the temporal and spatial regulation of transcription during pd.
similar expression patterns for regular genes and retrotransposons
the raw sequence reads were re-analyzed to quantify gene expression using tophat and cufflinks programs  <cit>  with updated genome annotation from ensembl  <cit> .
gene expression data is available in additional file 1: table s <dig>  hierarchical clustering was used to analyze the expression pattern of  <dig>  genes with sufficient change in expression across  <dig> individual cells at various developmental stages.
in addition to temporal dynamics, this data can also reveal variability among cells at the same stage.
see additional file  <dig> for more details.
ten gene clusters were defined according to the similarity in temporal expression patterns , similar to previous efforts based on dna microarray  <cit> .
see additional file 1: table s <dig> for gene lists.
cluster a includes  <dig> genes highly expressed in oocytes and quickly reduced at 2-cell  and 4-cell  stages.
these are mostly maternal mrnas undergoing degradation, which is evident from allele-specific mapping .
zga occurs during the 2c stage, evidenced by marked changes in gene expression between early and mid 2c stages .
the  <dig> transcripts in cluster b are exclusively expressed at 2c and 4c stages.
cluster d genes are induced at mid 2c, but their transient expression lasts until the 16-cell  stage.
cluster f genes are gradually upregulated at the 2c stage and downregulated in the blastocyst.
genes in clusters g to j are activated at various stages and remain highly expressed.
our goal is to find the regulatory mechanism behind these different patterns of expression.
during pd, much of the genome is de-methylated and transposons are actively transcribed.
previous studies have shown that their expression patterns are dynamic and stage-specific  <cit> .
to estimate their abundance, we used tetranscripts software  <cit>  with a special index file derived from repeatmasker data available at ucsc genome browser website  <cit> .
additional file 2: figure s <dig> shows the expression pattern of some of the highly transcribed transposons.
the additional file 1: table s <dig> contains the detailed expression levels of all tes.
the most highly expressed is a retrovirus-like element mt-int , a ltr element of the ervl-malr family.
it is expressed from the zygote to the 4c stage, similar to cluster a genes.
merv-l  is an ervl family ltr that is sharply induced by more than 500-fold at mid 2c before decreasing to low levels at 8c, showing an expression pattern similar to cluster b genes.
this is in agreement with a previous report that merv-l accounts for about 3% of total transcriptional output at this stage  <cit> .
there are several other ltr elements with this type of expression, including mt2c_mm, orr1a <dig>  orr1a <dig>  and mt2b <dig>  intracisternal a particle  elements are also transcribed between 2c and 16c as expected  <cit> , similar to cluster d genes.
besides ltrs, lines and sines are also expressed as expected.
sine transcripts increase modestly at the 2c stage and remain at that level through the blastocyst stage, similar to cluster f and g genes.
retrotransposons are transcribed in a highly regulated manner during pd, with expression patterns mirroring those of regular genes.
some merv-l elements seem to give rise to a microrna  from several genomic loci .
low level expression of mir- <dig> in escs has been detected  <cit> .
further study is needed to determine whether this microrna aids the clearance of maternal mrnas or transcripts from transposons.
figure 1b also indicates the heterogeneity in expression profiles among cells of the same stage.
this could be due to technical variations.
however, there is a 16c sample identified as 16cell_4- <dig> with expression profile similar to early 2c stage or zygote, characterized by higher expression of oocyte-specific markers such as oog <dig> , oog <dig> .
few reads  of this library map to paternal allele.
there is another cell  at mid-blastocyst with the expression of such genes.
reads from this cell contains substantial amount of paternal reads , but these oocyte-specific genes seems to mostly originate from maternal alleles.
some variations are also observed in the expression patterns of tes .
these may also be early signs of cell fate leading to different cell lineages  <cit> .
enrichment of transposable elements in promoters
we systematically studied the distribution of all repetitive elements  in the mouse genome  based on repeatmasker data available at ucsc genome browser website  <cit> .
covering 44% of the genome, the  <dig> , <dig> res are classified into  <dig> different types with unique names and consensus sequences in the repbase database  <cit> .
they belong to  <dig> repeat families which are grouped into  <dig> repeat classes.
figure 2a shows highly enriched res in the 2 kb promoters of  <dig> gene clusters.
b <dig> and b <dig> sine elements are overrepresented in the promoters of cluster c, d, f and g genes.
most of these genes are upregulated at the 2c stage .
this is congruent with the fact that transcription of b <dig> increases at the 2c stage.
for example,  <dig> % of genes in cluster c contain at least one b1_mus <dig> element in their promoter, which is much higher than the percentage  in cluster a. note that b1_mus <dig> is just one of many forms of b <dig> elements.fig.
<dig>
a enrichment of repeat elements in the promoters of genes.
original data matrix represents percentage of genes in a cluster containing a repeat in their promoter.
red indicates a certain repeat is enriched in the promoters of the gene clusters.
b
zfp <dig> gene contains a ltr element  around the translation starting site.
another ltr element  is in the intron.
these two positions correspond to known promoter regions marked as p <dig> and p <dig>  expression levels of zfp <dig>  obox <dig>  and retroelement mt2_mm are shown in , respectively
promoters of cluster d genes are  <dig> -fold enriched with mt2c_mm ltr compared with other genes with false discovery rate   <cit>  of  <dig>  × 10− <dig> .
the mt2c_mm retrotransposon itself is transcribed between 2c and 16c, similar to cluster d genes.
mt2b <dig> and mltr31d are also overrepresented in the promoters of cluster d genes.
although relatively enriched, these ltrs are only found near a small proportion of group d genes.
as shown in additional file 1: table s <dig>  about 5% of cluster d genes contain mt2c_mm elements, which is a  <dig> -fold enrichment compared to all other clusters combined.
promoters of cluster b genes are enriched with other ltr elements, namely mervl-int, orr1a2-int, orr1a3-int, mt2_mm, and mt2b <dig>  the most significant is a  <dig> -fold enrichment of merv-l  with fdr <  <dig>  × 10− <dig>  compared with genes in other groups combined.
this retrotransposon itself is transcribed at very high levels only in 2c-4c stages , similar to cluster b genes, suggesting a shared mechanism of regulation.
this agrees with previous reports  <cit> .
additional file 1: table s <dig> includes  <dig> genes that contain ervl elements in promoters and show an expression pattern similar to these ltrs.
as an example, transcription of zfp <dig>  starts at the middle of the ltr element mt2b <dig> .
this gene is sharply upregulated at mid 2c stage, before quickly decreasing to very low levels at the 4c stage .
this gene has been studied experimentally by liu et al.
<cit> , who found one major promoter p <dig> and an alternative, weaker promoter p <dig> in the intron.
p <dig> actually lies within the mt2b <dig> element and p <dig> is in another ltr in the intron .
zfp <dig> is likely generated by retrotransposition, similar to homologous pseudogene zfp353-ps  <cit> , where mrna sequences are reverse-transcribed and inserted into the genome.
retrocopies of mrna typically have no intron and are not expressed due to the lack of a functional promoter.
but a subsequent or preceding insertion of ltr elements upstream of the retrocopy can provide a promoter.
the expression of  <dig> such genes  was observed during pd when ltrs are actively transcribed.
these expressed retrogenes can be further studied.
early genes share motifs for obox homeobox factors
we scanned the proximal promoter region  of all the genes for transcription factor binding sites  using the comprehensive cis-bp database  <cit> , which includes thousands of binding motifs determined using protein binding microarrays or inferred based on shared protein domains.
out of the  <dig> binding motifs scanned, four were found enriched in the most highly induced gene group at mid 2c .
sharing the same “taatc” core sequence, these four motifs can be bound by obox <dig>  obox <dig>  otx <dig>  pitx <dig>  and other factors.
noticeably, the obox <dig> gene is upregulated by more than 100-fold at mid 2c, and quickly downregulated at 4c .
this is similar to the expression pattern of its potential target genes such as zfp <dig> , suggesting a potential regulatory role of obox factors.fig.
<dig> enrichment of transcription factor binding motifs in genes highly induced at mid 2-cell stage compared with early 2-cell stage.
a enriched motifs in the promoters of genes with highest fold-change at mid 2c.
g1-g <dig> represents  <dig> groups of  <dig> genes sorted by fold-change.
red indicates enrichment of motifs.
b the expression patterns of tfs that bind to the corresponding motifs shown in .
c binding motifs by obox families in ltrs belong to the ervl family.
d expression patterns of obox family homeobox factors are highly regulated during pd
oocyte specific homeobox  family tfs have not been studied in detail because they seem to be rodent-specific and are only expressed during pd  <cit> .
they are highly regulated in oocytes and cleavage stage embryos .
they are among the top  <dig> tfs when ranked by variance in gene expressions during pd.
since these tfs bind to similar motifs , it is difficult to distinguish the true driver of gene expression.
three members  are similar in sequence, and are treated as one group in our rna-seq mapping; they are highly expressed in oocytes and are maternally derived.
due to their higher expression at 2c and shared binding motifs, obox1/2/ <dig> could work together with obox <dig> to regulate gene expression.
experimental studies are needed to confirm whether these factors are redundant or have differences.
obox <dig> expression is elevated from 2c to 16c, and low in the blastocyst stage , in agreement with a previous study  <cit> .
obox binding sites are enriched in ervl family ltr elements that are overrepresented in the promoters of cluster b genes .
for example, the 493 bp consensus sequence of merv-l contains three obox binding sites, while the 521 bp mt2b <dig> contains four.
although not as highly transcribed as merv-l, mt2b <dig> is one of reliable predictors of nearby gene expression among ltrs, as the corresponding transcripts almost always start right from the repeat, similar to zfp <dig>  interestingly, the mt2b <dig> element in the main promoter of zfp <dig> is interrupted by a  tandem repeat of 21 bp sequence, which contains the obox binding motif .
these tandem repeats may be selectively retained during evolution.
furthermore, genes with multiple obox binding motifs are more highly induced at mid 2c when compared to genes with one or no such motifs .
the retrovirus-like ltrs contain their own promoters and enhancers.
these promoters may be bound by obox tf families to drive the expression of both retrotransposons and nearby genes.
among the genes with transient 2c expression is the zscan <dig> family .
zscan <dig> proteins are involved in telomere elongation and genomic stability in escs  <cit> , and have been reported to restore developmental potency in escs  <cit>  and to transiently activate embryonic genes in induced pluripotent cells   <cit> .
among the genes transiently expressed at 2c, zscan <dig> genes ranked highest in terms of number of obox <dig> binding sites, with three binding sites at 90% similarity level in the 350 bp around the tss.
further study is needed to verify if obox family tfs are upstream regulators of zscan <dig> genes, and whether induction of obox gene expression in adult cells could promote pluripotency.
thus, there is evidence that the poorly-characterized obox family induces gene expression at mid 2c to jump start zga.
even though obox <dig> mutants develop normally  <cit> , loss of function of other family members may be lethal to the embryo, due to potential effects on hundreds of downstream genes and even ltr transposons.
we also examined the expression of these genes in another single-cell rna-seq dataset  <cit> , and found that the patterns are similar, except that obox <dig> is expressed lower .
using this this dataset, we also identified the enrichment of the same “taatc” motif among genes upregulated in 2c stage .
the overrepresented “taatc” motif can also bound by other homeobox tfs like otx <dig> and pitx <dig> .
these tfs may also regulate gene expression, as their transcripts are also present in oocyte, even though at a much lower level .
recently, a similar motif was discovered by töhönen et al.
<cit>  in the promoters of  <dig> human genes upregulated in early zga using single-cell rna-seq.
this agreement may suggest a conserved mechanism of zga, and warrants experimental confirmation.
sines associated with broad genome activation
to further study expression change in late 2c, we ranked genes by their fold-changes from late 2c compared to mid 2c, and then divided them into  <dig> groups of  <dig> genes.
as shown in fig.
4a–c, the res differentially distributed in the promoters of these gene groups are mainly sine elements, including the alu family  and b <dig> family.
the most significantly associated elements are b1_mm, b1_mus <dig>  and b1_mus <dig> .
these b <dig> elements are highly similar in their sequences, with only a few base-pair differences in most cases .
with over  <dig>  copies, b <dig> is one of the most prevalent retrotransposons in the mouse genome.
similar to human alu elements, b <dig> originated from initial duplication of the 7sl rna  <cit> .
the distribution of sine, but not line, is conserved across species  <cit> .
b <dig> elements originated from trnas  <cit> .
it is possible that b <dig> and b <dig> elements play a role in gene regulation during pd, and their fast expansion may have benefited mammalian development.fig.
<dig> enrichment of sine elements in genes upregulated in the late 2-cell stage.
a average fold-change for  <dig> groups of  <dig> genes.
b relative enrichment  or depletion  of repeat elements in the promoters of genes in these groups.
c fdr values derived from chi-square tests of independency of repeat element frequency and gene groups.
d correlation of fold-change in 2c with the presence of b <dig> elements in gene groups.
each point in the plot represents one group of  <dig> genes.
the vertical axis represents the percentage of genes in the group with b1_mus <dig> in their promoters, while the horizontal axis represents average fold-change.
e association of repeats with fold-change during the late 2-cell stage.
the average fold change of genes with one or more repeats within 2 kb on either sense or antisense strand
figure 4b also shows that genes downregulated at the 2c stage are depleted in these elements.
the average fold-change observed in these groups is highly associated with the percentage of genes containing alu family elements in their promoters.
as shown in fig.
4d, the pearson’s correlation coefficient  is  <dig>   for one of such element b1_mm.
about 36% of the  <dig> genes that are upregulated by more than  <dig> fold contain b1_mm elements, which is much higher than the 18% observed in gene downregulated by  <dig> fold.
in addition to b1_mm elements, other b <dig> and b <dig> elements also show this trend.
using a different dataset of single-cell rna-seq data  <cit> , we were able to confirm this remarkably linear and consistent correlation .
multiple b <dig> elements in promoters are associated with stronger upregulation in a dosage-dependent manner .
as shown in fig.
5c, the  <dig> genes with five b <dig> elements or more  in promoters are upregulated by  <dig> -fold on average, which is significantly higher than the  <dig> -fold observed in  <dig> genes with  <dig> b <dig> elements .
this is in turn higher than genes with three b <dig> elements .
the effect of b <dig> elements is surprisingly linear.
each additional b <dig> element is associated with an approximately 20–40% upregulation .
among genes that are upregulated more than  <dig> fold at late 2c stage, 59% contain alu family repeats in promoter, which is much higher than the 41% observed other genes.fig.
<dig> presence of some repeats strongly correlated with activation of gene expression during 2c stage in a dosage- and distance -dependent manner.
the heights of the bars indicate average fold-change in late 2-cell stage compared to early 2-cell.
significant deviations from zero are indicated by stars.
the numbers on the bar represent the number of genes affected.
the mt <dig> repeats in  and  can be mt2b, mt2b <dig>  mt2b <dig>  mt2_mm, or mt2c_mm.
in  and , alu family repeats mostly represent b <dig> elements.
e result from multiple linear regressions shows that sine and ltr elements are correlated with gene expression change.
f correlation of the number of alu family repeats in promoters and broad gene expression in various cells and tissues.
error bar shows standard error calculated from replicates.
r <dig> and j <dig> are mouse embryonic stem cell lines
the effects of the b <dig> elements are also dependent on the distance to tss .
the correlation is stronger when b <dig> elements are located closer to tss.
further upstream, the effect of alu elements is weaker, diminishing after 5–7 kb.
this is similar to what was observed in the mt <dig> elements .
b1_mm and other mouse-specific b <dig> elements are more efficient than other alu family repeats common in muridae  or rodentia .
but it is difficult to isolate a specific type of element, as different subtypes of alu family repeats often co-appear.
to systematically investigate the effects of various repeats on zga, we used multiple linear regression to model 2c fold change as a function of the numbers of various kinds of repeats within the 2 kb promoter.
the model also included cpg islands and tfbs of  <dig> tfs determined by chip-seq  <cit> .
the results  confirm the effect of sine and ltr elements.
while the effects of ervl family ltrs are strong and only seen on the same strand, b <dig> elements have weaker effects on both strands.
genes with multiple c-myc sites show a bigger fold-change.
e2f <dig> binding sites are associated with weaker but significant upregulation in a larger number of genes.
the effect of alu elements and c-myc and e2f <dig> binding sites can be confirmed using independent single-cell rna-seq data  <cit>  .
sines are a major source of cpg dinucleotides in mammalian genomes.
it is possible that the effect of alu elements in promoters is through the contribution of cpg sites that affect epigenetic modifications.
linear regression analysis shows b <dig> elements have a much more significant  correlation with zga fold-change than those of cpg dinucleotides , or cpg islands .
also, gene clusters defined by methylation of promoter regions during pd  <cit>   have little in common with the gene clusters by expression .
therefore, the effect of b <dig> repeats cannot be fully explained by cpg sites.
sines correlate with gene expression in adult tissues and stem cells
it has been reported that alu family repeats are enriched near housekeeping genes  <cit> , which are both broadly and highly expressed across tissue types.
we calculated the correlation coefficient between expression levels and the number of alu family repeats in promoters among genes.
as shown in fig.
5e, there is a significant positive correlation for all cells and tissues.
the pcc dramatically increases from  <dig>  at early 2c to  <dig>  at late 2c, and gradually decreases to  <dig>  at the late blastocyst stage.
pcc is higher in embryonic cells than in adult tissues.
this association is independent of cpg islands .
occurrences of alu family repeats are a stronger predictor of gene expression in undifferentiated ipscs compared with day  <dig> definitive endoderm  <cit> .
a similar pattern of correlation is observed with b <dig> family repeats .
sines may play a role in regulating gene expression in both pre-implantation embryos and escs.
urrutia et al.
<cit>  found no association between alu content and peak  gene expression level across tissues.
however, we found a significant correlation  with average gene expression.
in addition, even among genes expressed in all tissues, more alu elements in the promoter are associated with higher average expression .
the same is true for tissue-specific genes .
promoters of esc–specific genes are rich in sines and low in cpg_island
to further delineate the role of alu family repeats in gene regulation, we re-analyzed rna-seq data  <cit>  of normal tissues from both fetal and adult mice, as well as escs.
we divided  <dig>  protein-coding genes into  <dig> groups using k-means clustering based on their expression pattern in various tissues/cells .
in addition to ubiquitously-expressed housekeeping genes, we identified many clusters of tissue-specific genes.
cluster  <dig> contains  <dig> genes specifically expressed in mesc cells.
as shown in additional file 2: figure s <dig>  this includes known tfs , as well as many other factors such as zscan <dig>  esrrb foxn <dig>  and sox <dig>  these genes contain many alu family repeats in their promoters .
as shown in fig.
6b, their average alu/b <dig> coverages of the promoters of these genes are as high as housekeeping genes, much higher than other tissue-specific gene clusters.
unlike the promoters of housekeeping genes, the promoters of esc-specific genes are less likely to have cpg islands.
for example, the alu-rich promoter regions of nanog and oct <dig> are shown in fig.
6c.
human orthologs of these two genes are also enriched with alu elements.
sine elements might be important for the expression of pluripotency related genes.fig.
<dig> genomic content of promoters of mesc-specific genes.
a results of k-means clustering of normal tissue gene expression show housekeeping genes  and tissue-specific genes.
b gene clusters are plotted by average coverage of alu/b <dig> repeats and cpg island coverage.
housekeeping genes  are high in both alu/b <dig> coverage and cpg island.
mesc specific genes  are high in alu/b <dig>  but low in cpg island.
gene clusters specifically expressed in testis, intestine, and placenta are lower in both.
c devoid of cpg islands, the promoters of pou5f <dig> and nanog are enriched with alu/b <dig> elements
sines correlate with zga in other species
single-cell rna-seq data of early human embryogenesis  <cit>  were used to investigate the correlation between transposons and zga, which occurs between 4- and 8-cell stages  <cit> .
regression analysis at the repeat family level shows a significant association between the presence of alu family repeats  and gene expression .
the most prevalent alu element, alujb, is significant on both sense and antisense strands.
presence of an alu element is associated with 15–48% upregulation.
figure 7a shows that genes with multiple alu elements in promoters are upregulated in a dosage-dependent manner, similar to what was observed in mouse.
as suggested by fig.
7a–b, genes with fewer than two alu elements in the promoter are downregulated at 8c, similar to what was observed in mouse.
this is in agreement with another single-cell rna-seq analysis of human preimplantation embryos, which shows that alu elements are enriched in upstream of tss of the  <dig> genes upregulated during pd  <cit> .
alu elements were found to contain binding motifs for pitx <dig> and tbx <dig>  <cit> .fig.
<dig>
a genes with more alus in promoter are upregulated in a dosage-dependent manner.
stars indicate significant difference from zero based on t-test.
b among gene groups defined by fold-change, highly expressed ones tend to contain more alus in promoters in humans.
c association of gene expression with trna family of sine repeats in the bovine genome.
stars indicate significant difference from zero.
d dna transposon, dna11ta1_dr, is associated with gene upregulation during zga in zebrafish in a dosage-dependent manner.
e mouse genes with multiple alu family repeats, mostly b <dig> elements, are associated with go: <dig>  intracellular part.
genes with l <dig> elements in their promoter, on the other hand, are depleted in genes related to intracellular part.
f genes with l <dig> elements are enriched with gpcr activity, while in genes containing alu elements, such genes are depleted.
g mouse genes with sine element in promoters are enriched in genes with yeast orthologs , compared to genes that does not contain such elements in promoters.
h among sine elements, alu family, mostly b <dig> elements, are enriched with genes with yeast orthologs
we also compared the occurrences of human alu and mouse b <dig> elements in promoters of orthologous gene pairs.
the number of alu elements in human gene promoters is highly correlated with the number of b <dig> elements in orthologous mouse gene promoters .
this has been noted as surprising, as b <dig> and primate alu elements replicated in these genomes independently  <cit> .
the differences in zga fold-changes between orthologous gene pairs in human and mouse are significantly associated with the differences in the number of alu family repeats in their promoters .
the rapid expansion of b <dig> in the mouse and alu in the human genome may contribute to gene expression divergence.
bovine zga is associated with trna sine repeats.
we used rna-seq data based on pooled embryos  <cit> .
the repeat families associated with expression change during zga between 8c and 16c are erv <dig>  simple repeat, and trna.
the trna family sine repeats in bovine have a weaker but significant association  based linear regression.
it is also dosage-dependent , as genes with three or more trna family repeats are more highly upregulated than genes with one or two such repeats .
the trna family repeats are associated with higher gene expression in both sense and antisense strands.
the most significant repeats are sine2-1_bt and sine2-2_bt, which are 120 bp bovine-specific sine repeats derived from trna.
the association with trna family repeats can also be confirmed using dna microarray data  <cit> .
zga in zebrafish is associated with at-rich dna transposons .
the zebrafish  genome  <cit>  is dominated by more than  <dig> million dna transposons.
there are fewer retrotransposons compared with mammals.
using rna-seq, harvey et al.
studied the zebrafish zga  <cit> , which happens at about  <dig>  h post-fertilization.
results from regression analysis  show that some dna transposons are highly associated with gene upregulation at zga.
the most significant is dna11ta1_dr, a non-autonomous dna transposon.
there are  <dig> genes containing this repeat in the 5 kb promoter region, and their expression is significantly higher than other genes .
the  <dig> genes with two or more dna11ta1_dr elements are induced at significantly  higher levels than the  <dig> genes with one element, which is in turn higher  than genes without such an element .
the accumulation of sines near genes is likely to result from a positive selection process  <cit> .
since there is no known mechanism to remove sines from genomes, it is difficult to explain the lack of sines in gene-poor regions  <cit> .
the transposons significantly correlated with zga are often prevalent in and specific to the host species.
specific transposons seem to be encouraged to expand during the course of evolution.
non-random distribution of transposons in the genome
we studied the distribution of all repeats across the mouse genome using the ebi approach.
while some res like b <dig> are prevalent, others are only observed dozens of times.
we found that the frequency follows lognormal distribution .
lognormal distribution implies that growth rate is independent of existing occurrence  <cit> .
the distribution of the distances between repeats is power-law like  <cit> , which could be expected as transposons often “copy-and-paste” to nearby loci and form clusters on the genome.
some repeats show enrichment and strand-preference near genes .
we found that sines are enriched in introns, promoters, and downstream regions, suggesting that sines are located near genes.
on the contrary, lines are depleted from these regions and are away from genes.
there are  <dig> types of ltrs that are specifically enriched in promoter regions.
some repeats demonstrate strand-preference.
for example, rltr10-int repeats align  <dig> -times more frequently on the same strand relative to the nearby gene than on the opposite strand, which is highly unusual .
there are a total of  <dig> repeats that are enriched in a strand-specific manner, including several prevalent ervl-malr family members , ervk family members , the ervl family , and the erv <dig> family .
this could be explained by new ltr elements generating new genes by activating retrogenes, as discussed.
we have shown that mt2b <dig> contains obox <dig> binding sites and is strongly associated with gene expression during pd.
other elements may regulate gene expression in other situations  <cit> .
surprisingly, most intronic retrotransposons are more likely found on the opposite strand of the host gene .
it is possible that intronic sequences, once spliced off transcripts in the nucleus, are spliced, reverse-transcribed, and inserted back into the genome, resulting in intronic retrotransposons on the antisense strand.
dna transposons do not show such a strong strand-specificity.
more investigation on the intronic strand-specificity is needed to verify this possible mechanism.
genes with multiple b <dig> in promoters form core cellular machinery
based on gene ontology   <cit> , we found that genes with multiple b <dig> elements in promoters are more likely to code for proteins that constitute intracellular parts  and less likely to be related to g-protein coupled receptor  activity, an extracellular signaling process .
on the contrary, genes with l <dig> line elements in their promoter are enriched in gpcr related genes and depleted of intracellular parts.
low complexity repeats are found to be enriched in promoters of genes related to the rna metabolic process .
thus, the distribution of repeats in the genome is related to the function of genes.
it is possible that the mouse embryo utilizes b <dig> elements to quickly establish core proteins during zga.
mouse genes with sines in promoters are more evolutionarily conserved.
figure 7g–h suggests that mouse genes with sines in the promoter are much more likely to have yeast orthologs .
this is similar to, but independent of, the effect of cpg islands .
a similar trend is observed in human genes .
thus, the distribution of transposons is correlated with the function of nearby genes.
the distribution of transposons may be regulated through selection.
key transcription factors in early development
taking advantage of high-resolution expression data  <cit> , we also systematically analyzed tfbs in the promoters of genes co-regulated at other stages beyond 2c.
enriched tfbs and expression patterns of corresponding genes are shown in additional file 2: figures s24–s <dig>  which are summarized in fig.
8a.fig.
<dig>
a transcription factors  identified during pd.
upregulated and downregulated factors are shown in red and purple, respectively.
homeobox tfs are underlined.
genes coding for the tfs that are upregulated or downregulated at the same stages provide more confidence and are shown in bold.
b a possible correlation between genome size and generation time across selected model organisms.
c a possible correlation between genome size and longevity among animals.
taxonomical classes are represented by different shapes and colors.
n =  <dig>  r =  <dig> , p =  <dig>  × 10−14
we identified several tfs known to regulate embryo development, such as sox <dig>  oct <dig> , and klf <dig>  <cit> .
in addition, many tfs in fig.
8a are upregulated at the same developmental stage as their potential target genes, thus giving more support to their involvement in gene regulation.
this includes obox <dig> and klf <dig>  which are upregulated at mid 2c, as well as nr2c <dig>  zscan <dig>  and elk <dig> at late 2c.
zscan <dig> is known to be expressed during pd  <cit>  and is involved in maintaining pluripotency in escs  <cit> .
similarly, many motifs enriched in the promoters of downregulated genes are bound by tfs that are downregulated at the same time.
some of the tfs are likely maternally derived with high expression in the zygote and reduced at the 2c stage: lbx <dig>  lbx <dig>  sebox, zfp <dig>  pou2f <dig>  stat <dig>  grhl <dig>  emx <dig>  prdm <dig>  hoxa <dig>  lhx <dig>  ebf <dig>  and e2f <dig>  sebox is known as one of the maternal effect genes , and its rna products are carried over from oocyte to regulation of gene expression in early pd  <cit> .
further study should verify whether other tfs in this list are megs.
the ggaa motif bound by the e <dig> transformation-specific  domain tfs is repeatedly identified as enriched at several stages.
one of the tfs, gabpa, has been shown to be involved in early embryogenesis  <cit> .
the most highly expressed ets domain tfs in oocyte and 2c is ehf.
elf <dig> is highly expressed at the blastocyst stage.
ets domain tfs are a large and conserved family of tfs involved in a variety of developmental processes in animals  <cit> .
similarly, the g-rich motifs are repeatedly identified at various stages.
this motif can be bound by klf <dig> or other tfs, such as klf <dig>  klf <dig>  sp <dig>  sp <dig>  zbtb7b, etc.
the most highly expressed is klf <dig>  which was reported to regulate lineage formation in the pre-implantation embryo  <cit> , alongside klf <dig>  <cit> .
to account for the heterogeneity among cells in the blastocyst stage, we used hierarchical clustering to divide them into three types/lineages, which were then analyzed separately.
sox <dig> and oct <dig> binding sites are enriched in promoters of genes upregulated in type  <dig> cells at the late blastocyst stage.
these markers indicate that type  <dig> cells likely correspond to the epiblast, which is derived from the inner cell mass  and leads to the embryo proper  <cit> .
figure 8a shows that two homeobox tfs, hesx <dig>  and six <dig>  are both upregulated together with their target genes in this type of cell.
hesx <dig> is believed to be downstream of multiple pluripotency related pathways  <cit> .
six <dig> is considered to be an oncogene and is known to be involved in embryonic muscle formation  <cit> .
in type  <dig> cells, sox <dig> binding sites are overrepresented in promoters of genes upregulated at the late blastocyst stage, when the sox <dig> gene itself is also upregulated.
sox <dig> directly promotes differentiation towards extraembryonic cells, which leads to the primitive endoderm  <cit> , which is also derived from icm but develops into the yolk sac.
therefore, type  <dig> cells are likely committed to the primitive endoderm.
msx <dig> and isx are two other homeobox tfs identified as inducing gene expression in these cells.
little is known about isx in embryogenesis, but it is highly induced at the late blastocyst stage in type  <dig> cells.
interestingly, we found msx <dig> and msx <dig> are associated with primitive endoderm and trophoblasts, respectively.
their differential role should be further studied.
type  <dig> cells may represent trophoblasts, which make up the outer layer of the blastocyst.
in the promoter of genes highly induced in this type of cell, we failed to detect binding sites for cdx <dig>  a key regulator for trophoblasts  <cit> .
but sites bound by msx <dig> are enriched.
msx <dig> and msx <dig> were shown to be critical for the interaction between the blastocyst and the uterus  <cit> .
these two proteins are highly conserved in mammals  <cit> , and mutations of msx <dig> and msx <dig> lead to failures in implantation  <cit> .
our analysis identified many other tfs that could potentially contribute to the complex gene regulatory network during pd.
for example, mecom is upregulated at early 2c and downregulated at 4c, along with its potential target genes.
mecom is highly expressed in the embryo, and mutant is embryonic lethal  <cit> .
figure 8a also shows that we identified many homeobox tfs, which are believed to be regulators of morphogenesis and development  <cit> .
in addition, many tfs in fig.
8a have been studied in relation to embryogenesis and cancer.
further study of these tfs will elucidate their role in gene regulation in pd.
discussion
transposons, mutagenesis, and the c-value paradox
some organisms in our ecosystem can finish a reproduction cycle in 20 min, while others require more than 10 years to reach sexual maturity.
if genetic mutations happen stochastically at similar rates, a huge imbalance in how quickly organisms evolve and adapt would result.
slow-reproducing organisms, therefore, are under pressure to find ways to dramatically promote mutagenesis and genotype diversity.
transposition of mobile dna elements may be a necessary “copy-and-paste” mechanism that promotes not only insertional mutations but also homologous recombination.
we showed here that they may also regulate many coding genes during early development and thus will have substantial influence of morphogenesis.
some retrotransposons are even active in somatic cells and lead to genetic mosaicism within individuals  <cit> .
there is additional evidence that tes are drivers of genome evolution  <cit> , rather than just “junk” dna.
following this argument, we would expect slow-reproducing organisms to have more tes in their genome, leading to larger genomes.
across model organisms, generation time seems to be proportional to genome size on a log-scale .
using larger datasets  <cit> , longevity and genome size are weakly but significantly correlated  across  <dig> animal species .
further study is needed to confirm the correlation in more organisms.
but this may shed some light on the c-value  paradox  <cit> : eukaryote haploid dna contents vary greatly, but are unrelated to organismic complexity.
even though tes can be disruptive for individuals, they might be necessary for the adaptation and survival of the species, especially in slow-reproducing organisms.
otherwise, it is hard to imagine how the accumulation of tens of thousands of transposons near essential genes could be tolerated over millions of years.
our analysis show that these tes may influence development from an early formative stage.
thus expansion of different tes that contribute to the rewiring of developmental pathways may facilitate speciation and adaptation.
gene regulation by transposons
tes can be involved in epigenetic regulation, as they can recruit the silencing machinery  <cit> .
many examples have been reported that retrotransposons can serve as tf binding sites to promote nearby gene expression.
in addition to ltrs, which contain promoters that can be used to drive expression of nearby genes  <cit> , lines elements contain an antisense promoter that can be used by nearby genes  <cit> .
mer <dig>  a dna transposon, was found to contribute to pregnancy-related gene network and its evolution  <cit> .
in escs, kunarso et al.
<cit>  reported that about 25% of pou5f <dig>  nanog and ctcf binding sites are provided by tes.
thus tes is important in the regulatory network of escs  <cit> .
more importantly, retroviral activity was found to be a hallmark of pluripotency  <cit> .
erv-derived ltr elements may be contributing to the gene regulatory network of innate immunity  <cit> .
adding to these results, this study systematically investigated the correlation of tes and the genomics reprogramming in pd.
we show that tes maybe play a more profound role than previously thought, affecting thousands of genes.
we also provide some evidence for the potential role of sines in activating housekeeping and esc-specific genes.
possible mechanisms of b <dig> and alu in gene regulation
mouse b <dig> and human alu elements originated from 7sl rna  <cit>  and contain rna polymerase iii promoters.
the a-box and b-box included in alu sequences are bound by a multi-subunit transcription factor tfiiic, to form the pol iii complex.
although some micrornas are shown to be transcribed by pol iii using upstream alus  <cit> , it is unlikely that pol iii would produce thousands of essential genes.
this would predict a strand-specific correlation and alternative tss’, similar to ltrs.
alu family repeats also contain binding sites for many factors associated with rna pol ii  <cit> , including p <dig>  <cit> , retinoic acid receptors  <cit> , yy <dig>  <cit> , pit <dig>  <cit> , etc.
our analysis shows that b <dig> elements also contain tfbs’ for obox family proteins, especially obox <dig>  these tfs may act upon alu elements to drive gene expression.
in addition, it is well documented that tfiiic, without the rest of the pol iii apparatus, has the so-called extra-transcriptional effects  ranging from nucleosome positioning, genome organization, and direct effect on pol ii transcription  <cit> .
alu provides the majority of the tfiiic binding sites in humans and mice  <cit> .
b <dig> elements are also enriched with ctcf binding motifs defined by chip-seq  <cit> .
some human alus serve as estrogen receptor -dependent enhancers for brca <dig>  <cit> .
su et al.
found that human alu elements in the proximal upstream region are more conserved and show many properties of enhancers  <cit> .
indeed, similar to enhancer rnas   <cit> , transcription of alu transposons may boost the expression of downstream genes through an enhancer mechanism specific to embryonic cells.
further study is needed to investigate these possible mechanisms.
transposons near genes should be treated more like regulatory elements.
in order for a single tf to regulate a large number of genes, it can evolve to take advantage of existing mobile elements.
this is more likely than the scenario where hundreds of genomic loci converge to the binding motif of an existing tf, which is especially true for highly specific motifs such as that of ctcf  <cit> .
CONCLUSIONS
guided by biological curiosity, exploratory bioinformatics analysis of the single-cell rna-seq and related data yields many actionable insights.
one of the surprising observations is that genes with similar expression patterns in early embryogenesis share specific transposons in their promoters.
during zga, while ltrs are linked to transient, forceful and early induction of several hundred genes, sine elements are associated with the upregulation of thousands of essential genes.
the machinery that transcribes retro-transposons may also be used to establish the expression landscape of early embryos.
this study also demonstrates the power of single-cell rna-seq, especially when applied to the study of normal developmental processes.

 hypothesis testing is the foundation of the scientific methodology. in biomedical research, focused and hypothesis-driven projects are especially encouraged due to the enormous complexity of the field. but little attention has been paid to how such hypotheses are generated. beyond intuition, a quantitative, evidence-based methodology for generating research hypotheses will increase the tempo of biomedical research. the idea of leveraging the “big data” of tens of thousands of large genome-wide datasets has been discussed. such data-driven, inductive methods can complement hypothesis-driven approaches to form an iterative process of ongoing research. our knowledge about many fundamental biological processes remains limited and fragmented. we can take advantage of the massive genome-wide datasets to learn about biological systems, akin to the study of an unknown planet in space exploration. first proposed by tukey, the goal of exploratory data analysis is to explore the data and find what it can tell us. it is open-ended and much broader than testing pre-defined hypotheses or building desired models. eda is a philosophically different approach to statistical analysis, not a new set of tools. in addition to data visualization, all statistical techniques can be used to ask various questions. eda can assess distributions, structures and dependencies that are useful for modeling. more importantly, eda can reveal important characteristics and trends, which can help formulate new hypotheses for further investigation. indeed, many big scientific discoveries are made accidentally. open-minded exploration of scientific data should be an essential part of a research project. eda can serve as a strategy to generate new hypotheses from biomedical data. the infrastructure is largely in place due to many open, collective efforts for sharing of data, annotation and software. the main challenge is to organically combine multiple datasets to gain a holistic understanding of complex biological systems. this requires interpretation of results within biological contexts and the ability to survey the literature while analyzing massive data. such broad, data-driven effort aimed at general understanding of a biological process can be referred to as exploratory bioinformatics investigation. this study represents an attempt of ebi on gene regulation in pre-implantation mouse embryos. in the initial stage of mammalian development, zygotes undergo maternal-to-zygotic transition, during which maternal factors are eliminated while mrna and protein synthesis using zygotic genome are initiated to take control of embryo development. a highly coordinated cascade of regulatory mechanisms unfolds rapidly to give rise to several cell lineages and the formation of blastocysts. understanding the complex process of pre-implantation development is important for both fertility related interventions as well as manipulation of embryonic stem cells. many gene expression studies of the early mammalian embryos has been carried out using expressed sequence tags, dna microarrays, rna-sequencing, and, more recently, single-cell rna-seq. these studies documented the dynamic waves of gene expression at different stages. using the powerful single-cell rna-seq technique, deng et al. analyzed hundreds of cells from embryos of mixed background mice and found strong evidence for random and widespread monoallelic expression. with sequence-level detail at single-cell resolution, this and other similar datasets can be used for in-depth study of gene regulation during pd. the epigenetic remodeling of maternal and paternal genomes was also revealed by dna methylation profiling. in zebra fish, transcription factors such as pou5f, nanog, and soxb were found to activate zygotic gene expression. but the molecular mechanisms of pd remain poorly understood. retrotransposon expression is a defining event in genome reprogramming during pd. transposable elements cover 30–50% of mammalian genomes. some are actively transcribed, even retrotransposed, as much of the genome is briefly hypomethylated in early embryos. the expression of retrotransposons is dynamic and stage-specific. evsikov et al. showed that retrotransposons, especially long terminal repeats, are abundantly represented in transcripts derived from mouse embryo at 2-cell stage. class iii ltrs such as merv-l family ltrs are transcribed at extremely high levels, accounting for about 3% of total transcriptional output at this stage. human endogenous retroviruses hervk ltr are expressed at zygotic genome activation at eight-cell stage. long interspersed elements, another major type of retrotransposon, are also expressed during pd, similar to short interspersed nuclear elements. interestingly, some endogenous retroviral activities have been found to be associated with and can serve as markers of pluripotency and totipotency in stem cell populations, underlying the importance of studying retrotransposon expression. transcription of retrotransposons can directly influence the expression of neighboring genes. by analyzing est libraries, peaston et al. analyzed the expression of different types of tes during pd and found that some tes, especially merv-l family ltrs, provide alternative 5′ first exons to chimeric transcripts in 2-cell mouse embryos. besides ltrs, lines and sines also contributed a small number of chimeric transcripts. lines have also been found to initiate fusion transcripts in other studies. more broadly, by analyzing cap-selected 5′ end of mouse and human transcripts from various embryonic and adult tissues, faulkner et al. found that 6–30% of all transcripts initiate from tes. these transcripts are often tissue-specific. studies on long noncoding rnas also identified ~ tes critical for the biogenesis of about 30% of total lncrnas sequences. tes can also influence gene expression indirectly by shaping the epigenetics landscape, as gc-rich sines tend to be found near genes while lines has the opposite distribution. some transcripts from tes give rise to small rnas that involve in post-transcriptional regulation during pd. line- rna can regulate its own expression. chip-seq data in human cells shows that tes contribute 25% of binding sites for tfs critical for embryonic stem cells such as pou5f nanog and ctcf. xie et al. compared expression dynamics of pd in human, mouse and bovine and found substantial difference in co-expression network, which could be partly due to the cis-regulatory modules provided by species-specific transposons. recently, tӧhӧnen et al. used single-cell rna-seq to analyze single-cells from early human development. they found that the promoters of genes upregulated early at 4-cell stage are enriched with alu elements that harbor motifs, including a “taatcc” core motif bound by pitx and otx homeobox family. these studies provide substantial evidence for tes’ broad role in regulating gene expression during pd. in this study, we use ebi approach to study the dynamic expression of both regular genes and tes in mouse pd. starting from the large single-cell rna-seq data of deng et al., our approach is to systematically observe the gene expression dynamics to help develop a broad understanding of the regulatory mechanisms of this complex process. our goal is to produce insights that can lead to novel, testable hypotheses on early embryo development. the data was analyzed alongside other expression and epigenetic studies of pd, using various tools and annotation databases. our analyses provide evidence for co-regulation of transposons and regular genes in early embryo development in mouse, followed up with similar observation in human, bovine and zebrafish. this adds to existing evidence that the expansion of species-specific tes may help rewire developmental pathways. motif analysis of promoter sequences of stage-specific genes identified many homeobox domain tfs as potential regulators of pd, which could be experimentally tested. we also examined the non-random distribution of tes in the mouse genome and the enrichment of sines in the promoters of genes specifically expressed in escs. finally, we discussed the evolutionary benefits of transposons in promoting genetic diversity, especially in slow-reproducing animals.fig. a exploratory bioinformatics investigation on gene regulation in mouse pre-implantation development. b hierarchical clustering of gene expression during pd. each of the rows represents a gene. columns correspond to samples labeled by developmental stages. red indicates expression levels higher than average for the row. expression lower than average is shown in green  to study allelic-specific transcription, deng et al. analyzed individual cells from mixed lineage mouse embryos using the single-cell rna-seq technique. excluding technical control samples, rna-seq libraries were generated, representing cells from zygote to blastocyst stages. the 2-cell stage is divided into early, middle and late phases. each library contains about million reads, mostly base pair long. with a total of about billion reads, this massive data enables in-depth ebi on the temporal and spatial regulation of transcription during pd. similar expression patterns for regular genes and retrotransposons the raw sequence reads were re-analyzed to quantify gene expression using tophat and cufflinks programs with updated genome annotation from ensembl. gene expression data is available in additional file 1: table s hierarchical clustering was used to analyze the expression pattern of genes with sufficient change in expression across individual cells at various developmental stages. in addition to temporal dynamics, this data can also reveal variability among cells at the same stage. see additional file for more details. ten gene clusters were defined according to the similarity in temporal expression patterns, similar to previous efforts based on dna microarray. see additional file 1: table s for gene lists. cluster a includes genes highly expressed in oocytes and quickly reduced at 2-cell and 4-cell stages. these are mostly maternal mrnas undergoing degradation, which is evident from allele-specific mapping. zga occurs during the 2c stage, evidenced by marked changes in gene expression between early and mid 2c stages. the transcripts in cluster b are exclusively expressed at 2c and 4c stages. cluster d genes are induced at mid 2c, but their transient expression lasts until the 16-cell stage. cluster f genes are gradually upregulated at the 2c stage and downregulated in the blastocyst. genes in clusters g to j are activated at various stages and remain highly expressed. our goal is to find the regulatory mechanism behind these different patterns of expression. during pd, much of the genome is de-methylated and transposons are actively transcribed. previous studies have shown that their expression patterns are dynamic and stage-specific. to estimate their abundance, we used tetranscripts software with a special index file derived from repeatmasker data available at ucsc genome browser website. additional file 2: figure s shows the expression pattern of some of the highly transcribed transposons. the additional file 1: table s contains the detailed expression levels of all tes. the most highly expressed is a retrovirus-like element mt-int, a ltr element of the ervl-malr family. it is expressed from the zygote to the 4c stage, similar to cluster a genes. merv-l is an ervl family ltr that is sharply induced by more than 500-fold at mid 2c before decreasing to low levels at 8c, showing an expression pattern similar to cluster b genes. this is in agreement with a previous report that merv-l accounts for about 3% of total transcriptional output at this stage. there are several other ltr elements with this type of expression, including mt2c_mm, orr1a orr1a and mt2b intracisternal a particle elements are also transcribed between 2c and 16c as expected, similar to cluster d genes. besides ltrs, lines and sines are also expressed as expected. sine transcripts increase modestly at the 2c stage and remain at that level through the blastocyst stage, similar to cluster f and g genes. retrotransposons are transcribed in a highly regulated manner during pd, with expression patterns mirroring those of regular genes. some merv-l elements seem to give rise to a microrna from several genomic loci. low level expression of mir- in escs has been detected. further study is needed to determine whether this microrna aids the clearance of maternal mrnas or transcripts from transposons. figure 1b also indicates the heterogeneity in expression profiles among cells of the same stage. this could be due to technical variations. however, there is a 16c sample identified as 16cell_4- with expression profile similar to early 2c stage or zygote, characterized by higher expression of oocyte-specific markers such as oog, oog. few reads of this library map to paternal allele. there is another cell at mid-blastocyst with the expression of such genes. reads from this cell contains substantial amount of paternal reads, but these oocyte-specific genes seems to mostly originate from maternal alleles. some variations are also observed in the expression patterns of tes. these may also be early signs of cell fate leading to different cell lineages. enrichment of transposable elements in promoters we systematically studied the distribution of all repetitive elements in the mouse genome based on repeatmasker data available at ucsc genome browser website. covering 44% of the genome, the, res are classified into different types with unique names and consensus sequences in the repbase database. they belong to repeat families which are grouped into repeat classes. figure 2a shows highly enriched res in the 2 kb promoters of gene clusters. b and b sine elements are overrepresented in the promoters of cluster c, d, f and g genes. most of these genes are upregulated at the 2c stage. this is congruent with the fact that transcription of b increases at the 2c stage. for example, % of genes in cluster c contain at least one b1_mus element in their promoter, which is much higher than the percentage in cluster a. note that b1_mus is just one of many forms of b elements.fig. a enrichment of repeat elements in the promoters of genes. original data matrix represents percentage of genes in a cluster containing a repeat in their promoter. red indicates a certain repeat is enriched in the promoters of the gene clusters. b zfp gene contains a ltr element around the translation starting site. another ltr element is in the intron. these two positions correspond to known promoter regions marked as p and p expression levels of zfp obox and retroelement mt2_mm are shown in, respectively promoters of cluster d genes are -fold enriched with mt2c_mm ltr compared with other genes with false discovery rate of  × 10−. the mt2c_mm retrotransposon itself is transcribed between 2c and 16c, similar to cluster d genes. mt2b and mltr31d are also overrepresented in the promoters of cluster d genes. although relatively enriched, these ltrs are only found near a small proportion of group d genes. as shown in additional file 1: table s about 5% of cluster d genes contain mt2c_mm elements, which is a -fold enrichment compared to all other clusters combined. promoters of cluster b genes are enriched with other ltr elements, namely mervl-int, orr1a2-int, orr1a3-int, mt2_mm, and mt2b the most significant is a -fold enrichment of merv-l with fdr <   × 10− compared with genes in other groups combined. this retrotransposon itself is transcribed at very high levels only in 2c-4c stages, similar to cluster b genes, suggesting a shared mechanism of regulation. this agrees with previous reports. additional file 1: table s includes genes that contain ervl elements in promoters and show an expression pattern similar to these ltrs. as an example, transcription of zfp starts at the middle of the ltr element mt2b. this gene is sharply upregulated at mid 2c stage, before quickly decreasing to very low levels at the 4c stage. this gene has been studied experimentally by liu et al., who found one major promoter p and an alternative, weaker promoter p in the intron. p actually lies within the mt2b element and p is in another ltr in the intron. zfp is likely generated by retrotransposition, similar to homologous pseudogene zfp353-ps, where mrna sequences are reverse-transcribed and inserted into the genome. retrocopies of mrna typically have no intron and are not expressed due to the lack of a functional promoter. but a subsequent or preceding insertion of ltr elements upstream of the retrocopy can provide a promoter. the expression of such genes was observed during pd when ltrs are actively transcribed. these expressed retrogenes can be further studied. early genes share motifs for obox homeobox factors we scanned the proximal promoter region of all the genes for transcription factor binding sites using the comprehensive cis-bp database, which includes thousands of binding motifs determined using protein binding microarrays or inferred based on shared protein domains. out of the binding motifs scanned, four were found enriched in the most highly induced gene group at mid 2c. sharing the same “taatc” core sequence, these four motifs can be bound by obox obox otx pitx and other factors. noticeably, the obox gene is upregulated by more than 100-fold at mid 2c, and quickly downregulated at 4c. this is similar to the expression pattern of its potential target genes such as zfp, suggesting a potential regulatory role of obox factors.fig. enrichment of transcription factor binding motifs in genes highly induced at mid 2-cell stage compared with early 2-cell stage. a enriched motifs in the promoters of genes with highest fold-change at mid 2c. g1-g represents groups of genes sorted by fold-change. red indicates enrichment of motifs. b the expression patterns of tfs that bind to the corresponding motifs shown in. c binding motifs by obox families in ltrs belong to the ervl family. d expression patterns of obox family homeobox factors are highly regulated during pd oocyte specific homeobox family tfs have not been studied in detail because they seem to be rodent-specific and are only expressed during pd. they are highly regulated in oocytes and cleavage stage embryos. they are among the top tfs when ranked by variance in gene expressions during pd. since these tfs bind to similar motifs, it is difficult to distinguish the true driver of gene expression. three members are similar in sequence, and are treated as one group in our rna-seq mapping; they are highly expressed in oocytes and are maternally derived. due to their higher expression at 2c and shared binding motifs, obox1/2/ could work together with obox to regulate gene expression. experimental studies are needed to confirm whether these factors are redundant or have differences. obox expression is elevated from 2c to 16c, and low in the blastocyst stage, in agreement with a previous study. obox binding sites are enriched in ervl family ltr elements that are overrepresented in the promoters of cluster b genes. for example, the 493 bp consensus sequence of merv-l contains three obox binding sites, while the 521 bp mt2b contains four. although not as highly transcribed as merv-l, mt2b is one of reliable predictors of nearby gene expression among ltrs, as the corresponding transcripts almost always start right from the repeat, similar to zfp interestingly, the mt2b element in the main promoter of zfp is interrupted by a tandem repeat of 21 bp sequence, which contains the obox binding motif. these tandem repeats may be selectively retained during evolution. furthermore, genes with multiple obox binding motifs are more highly induced at mid 2c when compared to genes with one or no such motifs. the retrovirus-like ltrs contain their own promoters and enhancers. these promoters may be bound by obox tf families to drive the expression of both retrotransposons and nearby genes. among the genes with transient 2c expression is the zscan family. zscan proteins are involved in telomere elongation and genomic stability in escs, and have been reported to restore developmental potency in escs and to transiently activate embryonic genes in induced pluripotent cells. among the genes transiently expressed at 2c, zscan genes ranked highest in terms of number of obox binding sites, with three binding sites at 90% similarity level in the 350 bp around the tss. further study is needed to verify if obox family tfs are upstream regulators of zscan genes, and whether induction of obox gene expression in adult cells could promote pluripotency. thus, there is evidence that the poorly-characterized obox family induces gene expression at mid 2c to jump start zga. even though obox mutants develop normally, loss of function of other family members may be lethal to the embryo, due to potential effects on hundreds of downstream genes and even ltr transposons. we also examined the expression of these genes in another single-cell rna-seq dataset, and found that the patterns are similar, except that obox is expressed lower. using this this dataset, we also identified the enrichment of the same “taatc” motif among genes upregulated in 2c stage. the overrepresented “taatc” motif can also bound by other homeobox tfs like otx and pitx. these tfs may also regulate gene expression, as their transcripts are also present in oocyte, even though at a much lower level. recently, a similar motif was discovered by töhönen et al. in the promoters of human genes upregulated in early zga using single-cell rna-seq. this agreement may suggest a conserved mechanism of zga, and warrants experimental confirmation. sines associated with broad genome activation to further study expression change in late 2c, we ranked genes by their fold-changes from late 2c compared to mid 2c, and then divided them into groups of genes. as shown in fig. 4a–c, the res differentially distributed in the promoters of these gene groups are mainly sine elements, including the alu family and b family. the most significantly associated elements are b1_mm, b1_mus and b1_mus. these b elements are highly similar in their sequences, with only a few base-pair differences in most cases. with over copies, b is one of the most prevalent retrotransposons in the mouse genome. similar to human alu elements, b originated from initial duplication of the 7sl rna. the distribution of sine, but not line, is conserved across species. b elements originated from trnas. it is possible that b and b elements play a role in gene regulation during pd, and their fast expansion may have benefited mammalian development.fig. enrichment of sine elements in genes upregulated in the late 2-cell stage. a average fold-change for groups of genes. b relative enrichment or depletion of repeat elements in the promoters of genes in these groups. c fdr values derived from chi-square tests of independency of repeat element frequency and gene groups. d correlation of fold-change in 2c with the presence of b elements in gene groups. each point in the plot represents one group of genes. the vertical axis represents the percentage of genes in the group with b1_mus in their promoters, while the horizontal axis represents average fold-change. e association of repeats with fold-change during the late 2-cell stage. the average fold change of genes with one or more repeats within 2 kb on either sense or antisense strand figure 4b also shows that genes downregulated at the 2c stage are depleted in these elements. the average fold-change observed in these groups is highly associated with the percentage of genes containing alu family elements in their promoters. as shown in fig. 4d, the pearson’s correlation coefficient is for one of such element b1_mm. about 36% of the genes that are upregulated by more than fold contain b1_mm elements, which is much higher than the 18% observed in gene downregulated by fold. in addition to b1_mm elements, other b and b elements also show this trend. using a different dataset of single-cell rna-seq data, we were able to confirm this remarkably linear and consistent correlation. multiple b elements in promoters are associated with stronger upregulation in a dosage-dependent manner. as shown in fig. 5c, the genes with five b elements or more in promoters are upregulated by -fold on average, which is significantly higher than the -fold observed in genes with b elements. this is in turn higher than genes with three b elements. the effect of b elements is surprisingly linear. each additional b element is associated with an approximately 20–40% upregulation. among genes that are upregulated more than fold at late 2c stage, 59% contain alu family repeats in promoter, which is much higher than the 41% observed other genes.fig. presence of some repeats strongly correlated with activation of gene expression during 2c stage in a dosage- and distance -dependent manner. the heights of the bars indicate average fold-change in late 2-cell stage compared to early 2-cell. significant deviations from zero are indicated by stars. the numbers on the bar represent the number of genes affected. the mt repeats in and can be mt2b, mt2b mt2b mt2_mm, or mt2c_mm. in and, alu family repeats mostly represent b elements. e result from multiple linear regressions shows that sine and ltr elements are correlated with gene expression change. f correlation of the number of alu family repeats in promoters and broad gene expression in various cells and tissues. error bar shows standard error calculated from replicates. r and j are mouse embryonic stem cell lines the effects of the b elements are also dependent on the distance to tss. the correlation is stronger when b elements are located closer to tss. further upstream, the effect of alu elements is weaker, diminishing after 5–7 kb. this is similar to what was observed in the mt elements. b1_mm and other mouse-specific b elements are more efficient than other alu family repeats common in muridae or rodentia. but it is difficult to isolate a specific type of element, as different subtypes of alu family repeats often co-appear. to systematically investigate the effects of various repeats on zga, we used multiple linear regression to model 2c fold change as a function of the numbers of various kinds of repeats within the 2 kb promoter. the model also included cpg islands and tfbs of tfs determined by chip-seq. the results confirm the effect of sine and ltr elements. while the effects of ervl family ltrs are strong and only seen on the same strand, b elements have weaker effects on both strands. genes with multiple c-myc sites show a bigger fold-change. e2f binding sites are associated with weaker but significant upregulation in a larger number of genes. the effect of alu elements and c-myc and e2f binding sites can be confirmed using independent single-cell rna-seq data. sines are a major source of cpg dinucleotides in mammalian genomes. it is possible that the effect of alu elements in promoters is through the contribution of cpg sites that affect epigenetic modifications. linear regression analysis shows b elements have a much more significant correlation with zga fold-change than those of cpg dinucleotides, or cpg islands. also, gene clusters defined by methylation of promoter regions during pd have little in common with the gene clusters by expression. therefore, the effect of b repeats cannot be fully explained by cpg sites. sines correlate with gene expression in adult tissues and stem cells it has been reported that alu family repeats are enriched near housekeeping genes, which are both broadly and highly expressed across tissue types. we calculated the correlation coefficient between expression levels and the number of alu family repeats in promoters among genes. as shown in fig. 5e, there is a significant positive correlation for all cells and tissues. the pcc dramatically increases from at early 2c to at late 2c, and gradually decreases to at the late blastocyst stage. pcc is higher in embryonic cells than in adult tissues. this association is independent of cpg islands. occurrences of alu family repeats are a stronger predictor of gene expression in undifferentiated ipscs compared with day definitive endoderm. a similar pattern of correlation is observed with b family repeats. sines may play a role in regulating gene expression in both pre-implantation embryos and escs. urrutia et al. found no association between alu content and peak gene expression level across tissues. however, we found a significant correlation with average gene expression. in addition, even among genes expressed in all tissues, more alu elements in the promoter are associated with higher average expression. the same is true for tissue-specific genes. promoters of esc–specific genes are rich in sines and low in cpg_island to further delineate the role of alu family repeats in gene regulation, we re-analyzed rna-seq data of normal tissues from both fetal and adult mice, as well as escs. we divided protein-coding genes into groups using k-means clustering based on their expression pattern in various tissues/cells. in addition to ubiquitously-expressed housekeeping genes, we identified many clusters of tissue-specific genes. cluster contains genes specifically expressed in mesc cells. as shown in additional file 2: figure s this includes known tfs, as well as many other factors such as zscan esrrb foxn and sox these genes contain many alu family repeats in their promoters. as shown in fig. 6b, their average alu/b coverages of the promoters of these genes are as high as housekeeping genes, much higher than other tissue-specific gene clusters. unlike the promoters of housekeeping genes, the promoters of esc-specific genes are less likely to have cpg islands. for example, the alu-rich promoter regions of nanog and oct are shown in fig. 6c. human orthologs of these two genes are also enriched with alu elements. sine elements might be important for the expression of pluripotency related genes.fig. genomic content of promoters of mesc-specific genes. a results of k-means clustering of normal tissue gene expression show housekeeping genes and tissue-specific genes. b gene clusters are plotted by average coverage of alu/b repeats and cpg island coverage. housekeeping genes are high in both alu/b coverage and cpg island. mesc specific genes are high in alu/b but low in cpg island. gene clusters specifically expressed in testis, intestine, and placenta are lower in both. c devoid of cpg islands, the promoters of pou5f and nanog are enriched with alu/b elements sines correlate with zga in other species single-cell rna-seq data of early human embryogenesis were used to investigate the correlation between transposons and zga, which occurs between 4- and 8-cell stages. regression analysis at the repeat family level shows a significant association between the presence of alu family repeats and gene expression. the most prevalent alu element, alujb, is significant on both sense and antisense strands. presence of an alu element is associated with 15–48% upregulation. figure 7a shows that genes with multiple alu elements in promoters are upregulated in a dosage-dependent manner, similar to what was observed in mouse. as suggested by fig. 7a–b, genes with fewer than two alu elements in the promoter are downregulated at 8c, similar to what was observed in mouse. this is in agreement with another single-cell rna-seq analysis of human preimplantation embryos, which shows that alu elements are enriched in upstream of tss of the genes upregulated during pd. alu elements were found to contain binding motifs for pitx and tbx.fig. a genes with more alus in promoter are upregulated in a dosage-dependent manner. stars indicate significant difference from zero based on t-test. b among gene groups defined by fold-change, highly expressed ones tend to contain more alus in promoters in humans. c association of gene expression with trna family of sine repeats in the bovine genome. stars indicate significant difference from zero. d dna transposon, dna11ta1_dr, is associated with gene upregulation during zga in zebrafish in a dosage-dependent manner. e mouse genes with multiple alu family repeats, mostly b elements, are associated with go: intracellular part. genes with l elements in their promoter, on the other hand, are depleted in genes related to intracellular part. f genes with l elements are enriched with gpcr activity, while in genes containing alu elements, such genes are depleted. g mouse genes with sine element in promoters are enriched in genes with yeast orthologs, compared to genes that does not contain such elements in promoters. h among sine elements, alu family, mostly b elements, are enriched with genes with yeast orthologs we also compared the occurrences of human alu and mouse b elements in promoters of orthologous gene pairs. the number of alu elements in human gene promoters is highly correlated with the number of b elements in orthologous mouse gene promoters. this has been noted as surprising, as b and primate alu elements replicated in these genomes independently. the differences in zga fold-changes between orthologous gene pairs in human and mouse are significantly associated with the differences in the number of alu family repeats in their promoters. the rapid expansion of b in the mouse and alu in the human genome may contribute to gene expression divergence. bovine zga is associated with trna sine repeats. we used rna-seq data based on pooled embryos. the repeat families associated with expression change during zga between 8c and 16c are erv simple repeat, and trna. the trna family sine repeats in bovine have a weaker but significant association based linear regression. it is also dosage-dependent, as genes with three or more trna family repeats are more highly upregulated than genes with one or two such repeats. the trna family repeats are associated with higher gene expression in both sense and antisense strands. the most significant repeats are sine2-1_bt and sine2-2_bt, which are 120 bp bovine-specific sine repeats derived from trna. the association with trna family repeats can also be confirmed using dna microarray data. zga in zebrafish is associated with at-rich dna transposons. the zebrafish genome is dominated by more than million dna transposons. there are fewer retrotransposons compared with mammals. using rna-seq, harvey et al. studied the zebrafish zga, which happens at about  h post-fertilization. results from regression analysis show that some dna transposons are highly associated with gene upregulation at zga. the most significant is dna11ta1_dr, a non-autonomous dna transposon. there are genes containing this repeat in the 5 kb promoter region, and their expression is significantly higher than other genes. the genes with two or more dna11ta1_dr elements are induced at significantly higher levels than the genes with one element, which is in turn higher than genes without such an element. the accumulation of sines near genes is likely to result from a positive selection process. since there is no known mechanism to remove sines from genomes, it is difficult to explain the lack of sines in gene-poor regions. the transposons significantly correlated with zga are often prevalent in and specific to the host species. specific transposons seem to be encouraged to expand during the course of evolution. non-random distribution of transposons in the genome we studied the distribution of all repeats across the mouse genome using the ebi approach. while some res like b are prevalent, others are only observed dozens of times. we found that the frequency follows lognormal distribution. lognormal distribution implies that growth rate is independent of existing occurrence. the distribution of the distances between repeats is power-law like, which could be expected as transposons often “copy-and-paste” to nearby loci and form clusters on the genome. some repeats show enrichment and strand-preference near genes. we found that sines are enriched in introns, promoters, and downstream regions, suggesting that sines are located near genes. on the contrary, lines are depleted from these regions and are away from genes. there are types of ltrs that are specifically enriched in promoter regions. some repeats demonstrate strand-preference. for example, rltr10-int repeats align -times more frequently on the same strand relative to the nearby gene than on the opposite strand, which is highly unusual. there are a total of repeats that are enriched in a strand-specific manner, including several prevalent ervl-malr family members, ervk family members, the ervl family, and the erv family. this could be explained by new ltr elements generating new genes by activating retrogenes, as discussed. we have shown that mt2b contains obox binding sites and is strongly associated with gene expression during pd. other elements may regulate gene expression in other situations. surprisingly, most intronic retrotransposons are more likely found on the opposite strand of the host gene. it is possible that intronic sequences, once spliced off transcripts in the nucleus, are spliced, reverse-transcribed, and inserted back into the genome, resulting in intronic retrotransposons on the antisense strand. dna transposons do not show such a strong strand-specificity. more investigation on the intronic strand-specificity is needed to verify this possible mechanism. genes with multiple b in promoters form core cellular machinery based on gene ontology, we found that genes with multiple b elements in promoters are more likely to code for proteins that constitute intracellular parts and less likely to be related to g-protein coupled receptor activity, an extracellular signaling process. on the contrary, genes with l line elements in their promoter are enriched in gpcr related genes and depleted of intracellular parts. low complexity repeats are found to be enriched in promoters of genes related to the rna metabolic process. thus, the distribution of repeats in the genome is related to the function of genes. it is possible that the mouse embryo utilizes b elements to quickly establish core proteins during zga. mouse genes with sines in promoters are more evolutionarily conserved. figure 7g–h suggests that mouse genes with sines in the promoter are much more likely to have yeast orthologs. this is similar to, but independent of, the effect of cpg islands. a similar trend is observed in human genes. thus, the distribution of transposons is correlated with the function of nearby genes. the distribution of transposons may be regulated through selection. key transcription factors in early development taking advantage of high-resolution expression data, we also systematically analyzed tfbs in the promoters of genes co-regulated at other stages beyond 2c. enriched tfbs and expression patterns of corresponding genes are shown in additional file 2: figures s24–s which are summarized in fig. 8a.fig. a transcription factors identified during pd. upregulated and downregulated factors are shown in red and purple, respectively. homeobox tfs are underlined. genes coding for the tfs that are upregulated or downregulated at the same stages provide more confidence and are shown in bold. b a possible correlation between genome size and generation time across selected model organisms. c a possible correlation between genome size and longevity among animals. taxonomical classes are represented by different shapes and colors. n =  r = , p =   × 10−14 we identified several tfs known to regulate embryo development, such as sox oct, and klf. in addition, many tfs in fig. 8a are upregulated at the same developmental stage as their potential target genes, thus giving more support to their involvement in gene regulation. this includes obox and klf which are upregulated at mid 2c, as well as nr2c zscan and elk at late 2c. zscan is known to be expressed during pd and is involved in maintaining pluripotency in escs. similarly, many motifs enriched in the promoters of downregulated genes are bound by tfs that are downregulated at the same time. some of the tfs are likely maternally derived with high expression in the zygote and reduced at the 2c stage: lbx lbx sebox, zfp pou2f stat grhl emx prdm hoxa lhx ebf and e2f sebox is known as one of the maternal effect genes, and its rna products are carried over from oocyte to regulation of gene expression in early pd. further study should verify whether other tfs in this list are megs. the ggaa motif bound by the e transformation-specific domain tfs is repeatedly identified as enriched at several stages. one of the tfs, gabpa, has been shown to be involved in early embryogenesis. the most highly expressed ets domain tfs in oocyte and 2c is ehf. elf is highly expressed at the blastocyst stage. ets domain tfs are a large and conserved family of tfs involved in a variety of developmental processes in animals. similarly, the g-rich motifs are repeatedly identified at various stages. this motif can be bound by klf or other tfs, such as klf klf sp sp zbtb7b, etc. the most highly expressed is klf which was reported to regulate lineage formation in the pre-implantation embryo, alongside klf. to account for the heterogeneity among cells in the blastocyst stage, we used hierarchical clustering to divide them into three types/lineages, which were then analyzed separately. sox and oct binding sites are enriched in promoters of genes upregulated in type cells at the late blastocyst stage. these markers indicate that type cells likely correspond to the epiblast, which is derived from the inner cell mass and leads to the embryo proper. figure 8a shows that two homeobox tfs, hesx and six are both upregulated together with their target genes in this type of cell. hesx is believed to be downstream of multiple pluripotency related pathways. six is considered to be an oncogene and is known to be involved in embryonic muscle formation. in type cells, sox binding sites are overrepresented in promoters of genes upregulated at the late blastocyst stage, when the sox gene itself is also upregulated. sox directly promotes differentiation towards extraembryonic cells, which leads to the primitive endoderm, which is also derived from icm but develops into the yolk sac. therefore, type cells are likely committed to the primitive endoderm. msx and isx are two other homeobox tfs identified as inducing gene expression in these cells. little is known about isx in embryogenesis, but it is highly induced at the late blastocyst stage in type cells. interestingly, we found msx and msx are associated with primitive endoderm and trophoblasts, respectively. their differential role should be further studied. type cells may represent trophoblasts, which make up the outer layer of the blastocyst. in the promoter of genes highly induced in this type of cell, we failed to detect binding sites for cdx a key regulator for trophoblasts. but sites bound by msx are enriched. msx and msx were shown to be critical for the interaction between the blastocyst and the uterus. these two proteins are highly conserved in mammals, and mutations of msx and msx lead to failures in implantation. our analysis identified many other tfs that could potentially contribute to the complex gene regulatory network during pd. for example, mecom is upregulated at early 2c and downregulated at 4c, along with its potential target genes. mecom is highly expressed in the embryo, and mutant is embryonic lethal. figure 8a also shows that we identified many homeobox tfs, which are believed to be regulators of morphogenesis and development. in addition, many tfs in fig. 8a have been studied in relation to embryogenesis and cancer. further study of these tfs will elucidate their role in gene regulation in pd. discussion transposons, mutagenesis, and the c-value paradox some organisms in our ecosystem can finish a reproduction cycle in 20 min, while others require more than 10 years to reach sexual maturity. if genetic mutations happen stochastically at similar rates, a huge imbalance in how quickly organisms evolve and adapt would result. slow-reproducing organisms, therefore, are under pressure to find ways to dramatically promote mutagenesis and genotype diversity. transposition of mobile dna elements may be a necessary “copy-and-paste” mechanism that promotes not only insertional mutations but also homologous recombination. we showed here that they may also regulate many coding genes during early development and thus will have substantial influence of morphogenesis. some retrotransposons are even active in somatic cells and lead to genetic mosaicism within individuals. there is additional evidence that tes are drivers of genome evolution, rather than just “junk” dna. following this argument, we would expect slow-reproducing organisms to have more tes in their genome, leading to larger genomes. across model organisms, generation time seems to be proportional to genome size on a log-scale. using larger datasets, longevity and genome size are weakly but significantly correlated across animal species. further study is needed to confirm the correlation in more organisms. but this may shed some light on the c-value paradox: eukaryote haploid dna contents vary greatly, but are unrelated to organismic complexity. even though tes can be disruptive for individuals, they might be necessary for the adaptation and survival of the species, especially in slow-reproducing organisms. otherwise, it is hard to imagine how the accumulation of tens of thousands of transposons near essential genes could be tolerated over millions of years. our analysis show that these tes may influence development from an early formative stage. thus expansion of different tes that contribute to the rewiring of developmental pathways may facilitate speciation and adaptation. gene regulation by transposons tes can be involved in epigenetic regulation, as they can recruit the silencing machinery. many examples have been reported that retrotransposons can serve as tf binding sites to promote nearby gene expression. in addition to ltrs, which contain promoters that can be used to drive expression of nearby genes, lines elements contain an antisense promoter that can be used by nearby genes. mer a dna transposon, was found to contribute to pregnancy-related gene network and its evolution. in escs, kunarso et al. reported that about 25% of pou5f nanog and ctcf binding sites are provided by tes. thus tes is important in the regulatory network of escs. more importantly, retroviral activity was found to be a hallmark of pluripotency. erv-derived ltr elements may be contributing to the gene regulatory network of innate immunity. adding to these results, this study systematically investigated the correlation of tes and the genomics reprogramming in pd. we show that tes maybe play a more profound role than previously thought, affecting thousands of genes. we also provide some evidence for the potential role of sines in activating housekeeping and esc-specific genes. possible mechanisms of b and alu in gene regulation mouse b and human alu elements originated from 7sl rna and contain rna polymerase iii promoters. the a-box and b-box included in alu sequences are bound by a multi-subunit transcription factor tfiiic, to form the pol iii complex. although some micrornas are shown to be transcribed by pol iii using upstream alus, it is unlikely that pol iii would produce thousands of essential genes. this would predict a strand-specific correlation and alternative tss’, similar to ltrs. alu family repeats also contain binding sites for many factors associated with rna pol ii, including p, retinoic acid receptors, yy, pit, etc. our analysis shows that b elements also contain tfbs’ for obox family proteins, especially obox these tfs may act upon alu elements to drive gene expression. in addition, it is well documented that tfiiic, without the rest of the pol iii apparatus, has the so-called extra-transcriptional effects ranging from nucleosome positioning, genome organization, and direct effect on pol ii transcription. alu provides the majority of the tfiiic binding sites in humans and mice. b elements are also enriched with ctcf binding motifs defined by chip-seq. some human alus serve as estrogen receptor -dependent enhancers for brca. su et al. found that human alu elements in the proximal upstream region are more conserved and show many properties of enhancers. indeed, similar to enhancer rnas, transcription of alu transposons may boost the expression of downstream genes through an enhancer mechanism specific to embryonic cells. further study is needed to investigate these possible mechanisms. transposons near genes should be treated more like regulatory elements. in order for a single tf to regulate a large number of genes, it can evolve to take advantage of existing mobile elements. this is more likely than the scenario where hundreds of genomic loci converge to the binding motif of an existing tf, which is especially true for highly specific motifs such as that of ctcf.  guided by biological curiosity, exploratory bioinformatics analysis of the single-cell rna-seq and related data yields many actionable insights. one of the surprising observations is that genes with similar expression patterns in early embryogenesis share specific transposons in their promoters. during zga, while ltrs are linked to transient, forceful and early induction of several hundred genes, sine elements are associated with the upregulation of thousands of essential genes. the machinery that transcribes retro-transposons may also be used to establish the expression landscape of early embryos. this study also demonstrates the power of single-cell rna-seq, especially when applied to the study of normal developmental processes. 
it also demonstrates that exploratory bioinformatics investigation can pinpoint developmental pathways for further study, and serve as a strategy to generate novel insights from big genomic data. the dynamics of gene expression can be partially explained by the enrichment of transposable elements in gene promoters and the similarity of expression profiles with those of corresponding transposons. starting with a single-cell rna-seq dataset consisting of mouse embryonic cells derived from zygote to blastocyst stages, we reconstructed the temporal and spatial gene expression pattern during pd. following this strategy, we re-analyzed a large body of genomic data to study the complex gene regulation in mouse pre-implantation development. together with other recent studies, our results provide further evidence that many transposable elements may play a role in establishing the expression landscape in early embryos. we also identified many potential transcription factors underlying pd and discussed the evolutionary necessity of transposons in enhancing genetic diversity, especially for species with longer generation time. such enhancer-like effects are also found for human alu and bovine trna sines. sines also seem to be predictive of gene expression in embryonic stem cells, raising the possibility that they may also be involved in regulating pluripotency. long terminal repeats are associated with transient, strong induction of many nearby genes at the 2- cell stages, probably by providing binding sites for obox and other homeobox factors. b and b sines are correlated with the upregulation of thousands of nearby genes during zygotic genome activation.
 instead of testing predefined hypotheses, the goal of exploratory data analysis is to find what data can tell us. following this strategy, we re-analyzed a large body of genomic data to study the complex gene regulation in mouse pre-implantation development.  starting with a single-cell rna-seq dataset consisting of mouse embryonic cells derived from zygote to blastocyst stages, we reconstructed the temporal and spatial gene expression pattern during pd. the dynamics of gene expression can be partially explained by the enrichment of transposable elements in gene promoters and the similarity of expression profiles with those of corresponding transposons. long terminal repeats are associated with transient, strong induction of many nearby genes at the 2- cell stages, probably by providing binding sites for obox and other homeobox factors. b and b sines are correlated with the upregulation of thousands of nearby genes during zygotic genome activation. such enhancer-like effects are also found for human alu and bovine trna sines. sines also seem to be predictive of gene expression in embryonic stem cells, raising the possibility that they may also be involved in regulating pluripotency. we also identified many potential transcription factors underlying pd and discussed the evolutionary necessity of transposons in enhancing genetic diversity, especially for species with longer generation time.  together with other recent studies, our results provide further evidence that many transposable elements may play a role in establishing the expression landscape in early embryos. it also demonstrates that exploratory bioinformatics investigation can pinpoint developmental pathways for further study, and serve as a strategy to generate novel insights from big genomic data. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords single-cell rna-seqexploratory data analysispre-implantation developmentearly embryogenesistransposonsrepetitive dnahttp://dx.doi.org/ /100000002national institutes of healthgm083226ge steven xijin http://dx.doi.org/http://dx.doi.org/ /100000001national science foundationiia-1355423issue-copyright-statement© the author 2017 