The Integrative Analysis of microRNA and mRNA Expression in Mouse Uterus under Delayed Implantation and Activation

Background Delayed implantation is a developmental arrest at the blastocyst stage and a good model for embryo implantation. MicroRNAs (miRNAs) have been shown to be involved in mouse embryo implantation through regulating uterine gene expression. This study was to have an integrative analysis on global miRNA and mRNA expression in mouse uterus under delayed implantation and activation through Illumina sequencing. Methodology/Principal Findings By deep sequencing and analysis, we found that there are 20 miRNAs up-regulated and 42 miRNAs down-regulated at least 1.2 folds, and 268 genes up-regulated and 295 genes down-regulated at least 2 folds under activation compared to delayed implantation, respectively. Many different forms of editing in mature miRNAs are detected. The percentage of editing at positions 4 and 5 of mature miRNAs is significantly higher under delayed implantation than under activation. Although the number of miR-21 reference sequence under activation is slightly lower than that under delayed implantation, the total level of miR-21 under activation is higher than that under delayed implantation. Six novel miRNAs are predicted and confirmed. The target genes of significantly up-regulated miRNAs under activation are significantly enriched. Conclusions miRNA and mRNA expression patterns are closely related. The target genes of up-regulated miRNAs are significantly enriched. A high level of editing at positions 4 and 5 of mature miRNAs is detected under delayed implantation than under activation. Our data should be valuable for future study on delayed implantation.


Introduction
Embryo implantation is a mutual interaction between blastocyst and uterus. The successful implantation of an embryo is dependent on both proper preparation of active blastocyst and receptive endometrium [1]. Delayed implantation is a developmental arrest at the blastocyst stage and a good model for deciphering the molecular interaction between embryo and uterus. There are around 100 species of mammals undergoing delayed implantation [2]. Because estrogen is essential for on-time uterine receptivity and blastocyst activation in mice [3], ovariectomy on day 4 of pregnancy will lead to blastocyst dormancy [4]. Many specific factors have been identified to be essential for embryo implantation through large-throughput analysis [5,6], and global gene expression in mouse uterus during delayed implantation and activation was also examined by Reese et al [5]. The global gene expression in mouse blastocysts during delayed implantation and activation was also reported [4]. However, the mechanism underlying delayed implantation and activation is still unclear.
Except for protein-coding RNAs, microRNAs (miRNAs) have been shown to be involved in mouse embryo implantation through regulating uterine gene expression [7,8]. Extensive sequence variations (isomiRs) for almost all miRNA and miRNA* species add additional complexity to the miRNA transcriptome [9]. RNA editing from A to I is widely present in human [10,11]. Additionally, this kind of editing was also detected in the seed sequences of miRNAs and may have effects on the recognition of target genes [11]. Illumina sequencing has opened the door for detecting and profiling known and novel miRNAs and mRNAs at unprecedented sensitivity. These latest high-throughput strategies permit high-resolution views of expressed miRNAs over a wide dynamic range of expression levels [9]. Direct sequencing also offers the potential to detect variations in mature miRNA length, as well as enzymatic modifications of miRNAs [12].
The large-scale proteomic analysis in mouse uterus during embryo implantation is still lacking. Because miRNAs can downregulate some of their targets not only at the translational but also at the transcriptional level [13], and the expression profiles of intragenic miRNAs and of their corresponding host genes are very similar both at the tissue and cellular level [14,15], it is therefore possible to use the paired expression analysis of miRNAs and mRNAs to identify mRNA targets of miRNAs. Serial analysis of gene expression (SAGE) is a high-throughput method for global gene expression analysis that allows the quantitative and simultaneous analysis of a large number of transcripts [16]. Therefore, the combination of SAGE and Illumina sequencing seems to be perfectly suited for deep transcriptome analysis [17].
This study was to have an integrative analysis on global miRNA and mRNA expression in mouse uterus under delayed implantation and activation through Illumina sequencing. We found that miRNA and mRNA expression patterns are closely related. A higher level of editing at positions 4 and 5 of mature miRNAs is detected under delayed implantation than under activation. The data from this study would provide a combined and comprehensive tissue-specific analysis of diverse miRNAs and transcriptional activity and also shed new light into the fine-tuning process of implantation.

Illumina sequencing of small RNAs
Total RNAs from mouse uteri under delayed implantation and activation were used to construct small RNA libraries for sequencing. The raw data from Illumina sequencing is available at Gene Expression Omnibus (GEO: GSE19473). In our two libraries, there were 5,334,521 reads for delayed implantation and 5,618,688 reads for activation, respectively. The read size was mainly ranged from 21 to 23 nt. The percentage of the 22 nt reads in total reads was 53.17% for delayed implantation and 54.10% for activation, respectively.
Read counts from delayed and activated uterus were normalized to tags per million (TPM) for each library. Differentiallyexpressed miRNAs were selected according to their fold changes (.1.2 fold), TPM of either library (.100) and p-values (,0.001). Based on the above-mentioned standards, there were 20 miRNAs up-regulated (Table 2) and 42 miRNAs down-regulated (Table 3). There were three up-regulated miRNA* sequences detected, including miR-17*, miR-145* and miR-21*.

Editing of mature miRNAs
There were different forms of editing in both libraries. T to A was the most dominant form in both libraries (2.75% for activation and 2.72% for delayed implantation), followed by A to T, C to T, G to T, T to C, T to G and G to C. Although A to G editing was a prevailing modification in some cell types [18,19], it was the eighth dominant editing form in total editing, representing 0.48% for activation and 0.55% for delayed implantation (Fig. 1A).
In both libraries, different kinds of nucleotide modification (designated as isoforms) were detected in mature miRNAs. The reference sequence in miRBase was not always the most dominant form in our study. In general, the editing in all of the miRNAs mainly occurred at positions 4, 5, 15, 19, and 21. The percentage of editing at positions 4 and 5 was significantly higher under delayed implantation than that under activation (z test, p,0.01). There was no difference between two libraries for the editing at other positions (Fig. 1B).
For let-7a, the editing rate at position 5 under delayed implantation was significantly higher than that under activation (Fig. 1C). Although the editing rate was very high at position 19, there was no difference between two libraries. Nucleotides at positions 4-5 were located right in the middle of the seed region (nucleotides 2-8) that is important for miRNA-mRNA binding. Such kind of nucleotide modification may alter genuine targets of let-7a.
For miR-143, there was a high percentage of editing at positions 4, 5, 13 and at 39-end (Fig. 1D). The editing rate at positions 4 and 5 in delayed implantation was also significantly higher than activation. However, the editing rate at position 13 and 39-end was similar between two libraries. Compared with both miR-143 and let-7a, the editing in miR-21 occurred in more positions, mainly at positions 4, 5, 9, 11, 13, 18, 19 and 21. But the percentage of editing at positions 4 and 5 under delayed implantation was also significantly higher than that under activation group, similar to let-7a and miR-143. (Fig. 1E).
Because the position 5 was located in the middle of seed sequence and the editing rate at position 5 of let-7a was significantly higher in delayed implantation than activation, TargetScan 4.2 (http://www.targetscan.org/) was used to predict target genes to see whether editing would affect the target genes. Only conserved binding sites (conserved in mouse, human, rat and dog) were considered. There were 614 target genes predicted for let-7a unedited form, while 236 target genes for edited form. There are only 27 target genes shared by both unedited and edited forms, suggesting that editing does change the genuine targets of let-7a.
Based on Gene Ontology analysis of the target genes of unedited and edited let-7a, the genes involved in ''biological process'' were compared between two groups. Compared to unedited let-7a, there were significantly more genes involved in nucleic acid metabolic process, channel activity, nucleus, biosynthetic process and transcription regulator activity in edited let-7a targets (Fig. 1F). Free energy changes between the wild type let-7a targets and edited let-7a targets were calculated by RNAeval algorithm implemented in the Vienna RNA package. When let-7a was edited from G to C at position 5, the binding between edited let-7a-5C and its target resulted in a net decrease of 24.5 kcal/mol in free energy. However, the free energy of binding between unedited let-7a and its targets was also 26.6 kcal/mol less than the binding between let-7a-5C and the targets of let-7a reference sequence.
All of the let-7a edited isoforms were listed in Table 4. The let-7a reference sequence was the dominant form, representing 56.1% of the total isoforms in delayed uterus. However, the most dominant form of miR-143 was one nucleotide deletion from 39end. miR-143 reference sequence was the second dominant form, representing 31.41% of total miR-143 isoforms (Table 5).
For miR-21, the reference sequence was the most abundant sequence. The number of tags at activation was slightly lower (6,756/7,054) than that of delayed implantation. However, the second most abundant sequence having one C addition at 39-end was significantly higher at activation than delayed implantation (8,427/2,313).The third most abundant sequence was one ''A'' deletion from 39-end, and the number of tags at delayed implantation was also slightly higher than that of activation. Because these three sequences were the same from positions 1 to 21 except for 39 deletion or addition, these sequences should be together detected by reference probe using Northern blot. If the total number of these three sequences were calculated, the tag number at activation was significantly higher than that of delayed implantation (18,964/13,581). The major difference between delayed implantation and activation was derived from the second sequence TAGCTTATCAGACTGATGTTGAC (Table 6).
In this study, the reference sequence of miR-21 at activation was slightly lower (6,756/7,054) than at delayed implantation. Because qRT-PCR was only designed for reference sequence of each miRNA, we used qRT-PCR to confirm whether miR-21 reference sequence was down-regulated in mouse uterus at activation. We found that the level of miR-21 at activation was indeed lower than that at delayed implantation. The miR-21 level in mouse uterus on day 5 of pregnancy was also checked by qRT-PCR. The level of miR-21 at implantation site was also lower than that at interimplantation sites (data not shown).

Target genes of edited let-7a
Compared to delayed implantation, let-7a was significantly upregulated under activation. Furthermore, G to C editing at position 5 was significantly higher under delayed implantation than under activation (667/6). Because G to C editing at position 5 was just located in the middle of the binding sequence 'seed sequence', further experiment was performed to examine whether this editing would shift the direction of targeting. Klf9, Gatm and Dnajb9 were predicted to be the target genes of unedited let-7a, whereas Tmem55a, Timp3 and Smad7 were predicted to be target genes of edited let-7a-5C. To confirm that these predicted target genes were indeed the target gene of their corresponding miRNAs, the 39-UTR segment of each gene was amplified by PCR from mouse cDNA and inserted into downstream of the luciferase reporter gene in the pGL3 control vector for the Dual-Luciferase assay. Compared to negative control, the luciferase activity containing 39-UTR of Klf9, Gatm and Dnajb9 was significantly inhibited by transfection with let-7a precursor, respectively ( Fig. 2A). Similarly, the luciferase activity containing 39-UTR of Tmem55a, Timp3 and Smad7 was also significantly inhibited by transfection with edited let-7a-5C precursor (Fig. 2B). Furthermore, compared to let-7a-5C precursor, the luciferase activity containing 39-UTR of Klf9, Gatm and Dnajb9 was significantly inhibited by transfection with let-7a precursor, respectively ( Fig. 2A). Conversely, compared to let-7a precursor, the luciferase activity containing 39-UTR of Tmem55a, Timp3 and Smad7 was significantly inhibited by transfection with edited let-7a-5C precursor (Fig. 2B).
In order to examine whether let-7a-5C could regulate its target genes in cultured mouse uterine cells, mouse uterine stromal cells were transfected with let-7a-5C precursor or let-7a precursor and cultured for 24 h. The expression of Tmem55a, Timp3 and Smad7 was determined by qRT-PCR. Compared to let-7a precursor, both Timp3 and Smad7 were significantly inhibited by let-7a-5C precursor (Fig. 2C). Additionally, Smad7 protein was significantly inhibited by let-7a-5C in comparison with let-7a precursor (Fig. 2D, E). There was no detectable change for Tmem55a between let-7a and let-7a-5C treatments (Fig. 2D).
Because let-7a-5C was significantly up-regulated under delayed implantation, we checked our SAGE-Illumina sequencing data to examine whether there was a reverse relation between let-7a-5C and its three target genes. In our SAGE-Illumina sequencing data, both Smad7 and Tmem55a were up-regulated under activation, which was opposite to the expression of let-7a-5C (Fig. 3A). When qRT-PCR was used to measure the expression of three target genes in mouse uterus, both Smad7 and Tmem55a were also up-regulated under activation (Fig. 3A). However, Timp3 was down-regulated under activation in our SAGE-Illumina sequencing data, which was also confirmed by qRT-PCR (Fig. 3A). Our data suggest that Timp3 may be not regulated by let-7a-5C at least at the level of mRNA expression in mouse uterus.
In order to examine whether the expression of Tmem55a, Timp3 and Smad7 was similar between activation of delayed implantation and implantation sites, the expression of Tmem55a, Timp3 and Smad7 was also examined in mouse uterus on day 5 of pregnancy. Compared to inter-implantation sites, both Smad7 and Tmem55a were also highly expressed at implantation sites. However, Timp3 expression didn't significantly change between implantation and inter-implantation sites (Fig. 3B). Because Smad7 was verified as a target gene of let-7a-5C, its protein level was also examined by Western blot. Compared to delayed implantation, Smad7 protein was up-regulated under activation (Fig. 4A, B). The level of Smad7 protein at implantation sites was also stronger than that at inter-implantation sites (Fig. 4C, D). By in situ hybridization, Smad7 mRNA expression was mainly localized in the subluminal stromal cells under activation of delayed implantation, but no Smad7 signal was detected in the uterus under delayed implantation (Fig. 4E).

Novel miRNAs
In our study, a high percentage of small RNAs were sorted as unknown RNAs, 23.78% for delayed implantation and 21.34% for activation. We would like to check whether novel miRNAs were present among unknown RNAs. Novel miRNAs were predicted by miRDeep. The miRDeep cut-off score was set at 1. Based on miRDeep software, 6 novel miRNAs were predicted. The total tag number for all of 6 novel miRNAs in both delayed implantation and activation was 538, accounting for 0.02% of the unknown RNAs and 0.005% of the total RNAs. Compared to mature miRNA strand, the percentage of miRNA* sequences was significantly less. For the miRNA* sequences of the novel miRNAs, there were 7 tags for delayed implantation and 22 tags for activation. The number of the tags for nov-miRNA-4 was 75 for delayed implantation and 76 for activation (Table 7). Since these novel miRNAs were expressed at a very low level in mouse uterus, whether these novel miRNAs are functional is still unknown.
Because novel mouse miRNAs were expressed in the delayed and activated uterus, qRT-PCR was used to verify their expression in mouse uterus. For novel miRNAs, primers were synthesized according to the paper previously published with some modifications [19]. Compared to delayed implantation, nov-miRNA-1, nov-miRNA-3, nov-miRNA-4, and nov-miRNA-6 were significantly down-regulated in mouse uterus under estrogen activation,    whereas nov-miRNA-2 was significantly up-regulated under estrogen activation. The expression of nov-miRNA-5 was not significantly different between delayed implantation and activation (Fig. 5B).
In order to examine whether these novel miRNAs were cotranscribed with their precursors, the expression of their precursors was also measured. The primers for their precursors were designed based on the pre-miRNA sequences predicted by miRDeep software. The expression trend of their precursors was very similar to that of their corresponding miRNAs, except for pre-miR-6, which was significantly up-regulated under activation (Fig. 5C).

Digital gene expression from SAGE-Illumina sequencing
To study the relationship between the miRNA and their targets, we performed SAGE-Illumina sequencing to examine the transcriptional profile of the uteri from delayed implantation and activation. The raw data from SAGE-Illumina sequencing is available at Gene Expression Omnibus (GEO: GSE19473). There were 9,912,459 and 12,869,487 reads sequenced from delayed implantation and activation, respectively. After removing the reads with 0/1 or 1/0 in both libraries, there were 8,683,980 and 11,039,265 meaningful reads remained for delayed implantation and activation, respectively. Of the meaningful reads, 356,981 and 348,124 reads were mapped to unigenes for delayed implantation and activation, respectively (Table 8). Among the total unique reads, 51,727 and 51,766 unique reads were identified for delayed implantation and activation, respectively. Of the unique reads, 45,147 (87.28%) and 45,144 (87.21%) unique reads were mapped to one gene for delayed implantation and activation, respectively. There were 10.54% tags for delayed implantation and 10.55% for activation matched to two genes. The remaining 2.18% for delayed implantation and 2.24% for activation matched to multiple genes (Table 8).
In both libraries, the top 30 most abundant tags were mainly matched to genomic repeat sequences and ribosome-related proteins (Table S2). A high percentage of the top 30 abundant tags were shared in both libraries. The counts from the two libraries (delayed and activated uterus) for each gene were compared by z-test and Bonferroni multiple test correction. Genes were designated to be significantly differentially expressed if the p value was ,0.001, and there was at least a 1.5-fold change in sequence counts between the two libraries. Under these standards, there were 3,033 genes up-regulated and 1,417 genes downregulated during activation. However, based on the standards that TPM was .100 in either library and the ratio of activation over delayed implantation was $2, there were 268 genes up-regulated and 295 genes down-regulated (Table 9, full list in Table S3).

Integrative analysis of miRNA and mRNA expression data
Because most mammalian miRNAs are intragenic and transcribed as part of their hosting transcription units [20], we hypothesized that the expression profiles of mature miRNAs and their host genes are directly correlated. miRNA expression was compared with their host mRNA expression to see whether they were co-expressed. The list of mouse intragenic miRNAs and corresponding host genes was retrieved from miRBase (Release 13.0, March 2009). Only those genes were considered as host genes if RefSeq sequences were overlapped with the miRNA either in introns, exons or UTR and were transcribed on the same strand as the miRNA. Based on our data from both miRNA Illumina sequencing and mRNA SAGE-Illumina sequencing, we found that miRNA expression was tightly correlated with the host mRNA expression (r = 0.35, p = 0.002 in delayed implantation and r = 0.37, p = 0.001 in activation), suggesting that both miRNA and the host mRNA were co-transcribed (Fig. 6A, B).
In this study, 9 miRNAs were up-regulated and 16 downregulated in the activated uterus compared to delayed implantation (fold.1.5, TPM.100 and p,0.001). There were 268 upregulated genes and 295 down-regulated genes in the activated uterus compared to delayed implantation (fold.2, TPM.100 and p,0.001). For comprehensive prediction of miRNA target genes, the results from two publically available algorithms (TargetScan and PITA) were merged. In total, 53 genes (with repeats, 49 unique genes) were predicted to be the targets of up-regulated miRNAs in activated uterus, while 196 (with repeats, 110 unique genes) genes were predicted as targets of down-regulated miRNAs. In our gene expression data, there were 52.4% down-regulated genes and 47.6% up-regulated genes. Among the target genes of the up-regulated miRNAs predicted by either TargetScan or  PITA, 63.93% of the target genes predicted was really downregulated in our study, suggesting the target genes of up-regulated miRNAs during activation were significantly enriched. However, 40.0% of the target genes of down-regulated miRNAs were upregulated in our study, which is consistent to our gene expression data (47.6%), suggesting that the target genes of down-regulated miRNAs were not enriched (Fig. 6C).
We defined a coherent target of a miRNA as a predicted target if its expression had a reverse pattern with the miRNA. Among the predicted targets of up-regulated miRNAs, down-regulated genes detected by SAGE-Illumina sequencing are considered as coherent targets and otherwise as non-coherent targets. For the down-regulated miRNA targets, up-regulated genes are coherent and otherwise non-coherent. Therefore, 44 genes (with repeat, 39 unique genes) were coherent target of up-regulated miRNAs in Figure 2. Prediction and confirmation of the target genes predicted for both let-7a (unedited) and let-7a-5C (edited). (A) The confirmation of the target genes (Klf9, Gatm and Dnajb9) of let-7a in mouse 3T3 cells using dual-luciferase assay. Cells were co-transfected with negative control, let-7a (Pre-let-7a) or let-7a-5C (Pre-let-7a-5C) precursor, respectively; (B) The confirmation of the target genes (Tmem55a, Timp3 and Smad7) of let-7a-5C in mouse 3T3 cells; (C) The relative mRNA expression level of the three target genes of edited let-7a-5C in cultured mouse stromal cells transfected with negative control, Pre-let-7a or Pre-let-7a-5C; (D) Protein level of Smad7 was detected by Western blot after mouse stromal cells were transfected with negative control, Pre-let-7a or Pre-let-7a-5C; (E) Quantification of Smad7 protein expression in (D). The difference between Pre-let-7a and Pre-let-7a-5C was compared using t-test, and the significant difference between two groups was labeled with asterisk. doi:10.1371/journal.pone.0015513.g002  Table S4.
Based on Gene Ontology analysis for the differentially expressed genes, the genes involved in cell cycle, response to stress, and metabolic process were significantly enriched in activation group compared to delayed implantation (Fig. 7A). When considering the coherent target genes of the down-regulated miRNAs, the genes involved in cell adhesion were significantly enriched among the up-regulated genes (Fig. 7B).

General comparison with other deep sequencing data
In our study, the percentage of miRNAs is 74.52% in delayed and 76.57% in activated uteri, respectively. The top 3 most abundant miRNAs are let-7c, let-7f and let-7a. These data are similar to that in mouse oocytes [21]. In amphioxus, the length distribution peaked at 22 nt and almost half of these clean reads (45.11%) are 22 nt in length [22]. In our study, 22 nt is also the dominant size among the small RNAs, consistent with the common size of miRNAs. Additionally, we found that the most abundant sequence of miRNA-143 is not the reference sequence reported. Kuchenbauer et al. [9] also found many isoforms for miR-181a and that the miRBase reference sequence was not the dominant species. However, qRT-PCR widely used for miRNA quantification is mainly based on reference sequence [19]. Therefore, the amount of the reference sequence from qRT-PCR may not reflect the amount of the dominant form of each miRNA. In our previous study, miR-21 at implantation site is up-regulated compared to interimplantation sites based on Northern blot [8]. However, by using qRT-PCR, miR-21 at implantation is down-regulated compared to  inter-implantation site in this study. Similar situation happens for miR-21 during delayed implantation and activation. Our data suggest that it should be cautious when miRNA expression level is examined by different methods. Since data from in situ hybridization and Northern blot are from the hybridization between probes and matched sequences, their data should be more closely related. The data from qRT-PCR should be matched with the reference sequence in Illumina sequencing data.
Editing and possible significance of mature miRNAs In our study, 41.5% of let-7a and 64.4% of miR-143 sequences are either edited or alternatively spliced. Among 26 cell types from human, mouse and rat, there are approximately 20% of miRNA mismatches compared with their genomic sequences, including A to I editing (identified as A to G editing), and 39 terminal A and U additions [23]. Through massively parallel sequencing, ,50% of the mature miRNAs in the let-7 family display internal insertion/ deletions and substitutions when compared to precursor miRNA and the mouse genome reference sequences [10]. For let-7a in mouse ovary, there are 35% exhibited terminal nucleotide additions and/or excisions among the sequences that aligned with pre-miR precursors [10]. A to I editing is the dominant one in several reports [10,11]. Although A to I editing was indeed found in our data, but not the dominant one. T to A editing is the most dominant form among all of the miRNAs in our study. It is shown that A to I editing in mice is often lower than that in human [23]. It is possible that the mechanism of miRNA editing may be different among mammals. The mechanism on high percentage of T to A editing is still unknown.
IsomiRs resulting from variations at the 59-end may be of particular interest as they have different seed sequences from the reference miRNA. A to I editing sites also occur within the seed region of mature miRNA sequences, showing that RNA editing can impact miRNA target recognition and function [11,24]. In our study, a significantly higher percentage of editing occurs within 59 seed region in delayed implantation group in let-7a, miR-143 and miR-21 compared to activation group. Based on our data, the target genes of unedited miRNAs are largely different from that of edited ones. Once let-7a is edited as let-7a-5C, different sets of target genes are predicted and three target genes are confirmed by luciferase analysis in our study, suggesting that a single G to C base change is sufficient to redirect silencing miRNAs to a new set of targets. For let-7a, miR-143 and miR-21, the editing within seed region in delayed implantation is significantly higher than activation group. However, the signifi-   [24]. For the physiological significance, let-7a-5C is highly expressed during delayed implantation compared to activation group. Once let-7a is edited into let-7a-5C, there are more genes involved in nucleic acid metabolic process, channel activity, nucleus, biosynthetic process and transcription regulator activity, suggesting that these functions should be suppressed during delayed implantation. This is consistent with low metabolic activity during delayed implantation [4]. Smad7 was verified as a target gene of let-7a-5C in our study. Both Smad6 and Smad7 prevent ligand-induced activation of signaltransducing Smad proteins in the transforming growth factor-b family. In cardiac myofibroblasts, ectopic Smad7 protein is associated with accelerated activation of pro-MMP-2 into MMP-2 [25]. Proper extracellular matrix degradation and blastocyst invasion are essential for embryo implantation and decidualization [26,27]. MMP-2 is a matrix metalloproteinase and strongly expressed in the stromal cells of mouse uterus during implantation period [27,28]. In our study, Smad7 expression is mainly localized in the subluminal stromal cells at implantation sites under activation, suggesting that Smad7 and MMP-2 should be colocalized at subluminal stromal cells at implantation sites. Therefore, the edited let-7a-5C may play a role for mediating a proper balance between MMP-2 level and embryo invasion for the successful pregnancy through Smad7.

Comparison with implantation-related miRNAs
In our study, there are 20 up-regulated miRNAs and 42 downregulated miRNAs at least 1.2 folds in activation. Up to date, there is no miRNA expression profile available for comparison with our data from delayed implantation. The closest data for comparison is from our previous paper on miRNA expression from implantation and inter-implantation sites in mouse uterus [8]. Compared to inter-implantation sites, there are 13 miRNAs up-regulated at least 2 folds at implantation sites. Of which, let-7f, let-7d, let-7e, let-7i, miR-20a, miR-298, let-7g, miR-21, and let-7a are up-regulated in activation group. However, miR-26a, let-7b and let-7c and miR-143 up-regulated at implantation sites either doesn't change or is down-regulated in activation group. The reason of this discrepancy is not clear. This may reflect the difference between implantation sites during early pregnancy and activation after delayed implantation. Additionally, miRNAs detected by microarray are based on reference sequences, however, the reference sequences are often not the dominant ones for some miRNAs, which might be another reason causing that discrepancy.

Comparison with implantation-related mRNA
Ultra-high-throughput sequencing is emerging as an attractive alternative to microarrays for genotyping, analysis of methylation patterns, and identification of transcription factor binding sites [29]. In 10 bp SAGE, around 53,000 tags are obtained from each library [6]. In this study, there are 356,981 tags for delayed implantation and 348,124 tags for activation. The number of tags detected is greatly increased in this study. Additionally, about 50% of the tags are single-matched to genes for 10 bp SAGE [6], while about 87% of tags are single-matched to genes in our 16 bp SAGE. Owing to its increased tag length, long SAGE tags are more reliable in direct assignment to genome sequences. In our data, only about 12% of tags are multiple-matched. Of these tags most tags are mainly B2 repeats and ribosomal proteins.
For delayed implantation, there is only one large scale study on gene expression in mouse uterus using Affymetrix murine expression arrays [5]. They reported 41 down-regulated and 21 up-regulated genes during activation compared to delayed implantation. In our study, there are much more differentially expressed genes detected, including 268 genes up-regulated and 295 genes down-regulated at activation compared with delayed implantation. All of the 21 up-regulated genes in their study are also up-regulated in our study. Among 41 down-regulated genes in their study, 37 genes are also down-regulated in our study. Ctss and Mecp2, down-regulated in their study, don't change in our study. Akr1b7 (24.41 vs 3.15) and Gzma (22.15 vs 14.82), down-regulated in their study, are up-regulated in our study, respectively. We performed qRT-PCR to solve this discrepancy(data not shown). Based on our qRT-PCR data, Akr1b1 and Gzma are indeed upregulated at activation group compared to delayed implantation. Mecp2 and Ctss are slightly down-regulated at activation group compared to delayed implantation. Additionally, the whole uterus was used in their study [5], but only the implantation sites of mouse uterus following activation were used in our study, which may also reflect the difference between these two studies. This comparison further shows that our data from SAGE-Illumina sequencing is of great value for further understanding the mechanism on delayed implantation.
Based on Gene Ontology analysis, the genes involved in metabolic process, response to stress and cell cycle are significantly enriched in the up-regulated genes of activation group compared to down-regulated genes. When considering the targeting genes of differentially expressed miRNAs, the genes involved in metabolic process, response to stress and vasculogenesis are significantly enriched among the differentially expressed genes in both groups. Among the up-regulated genes, the genes involved in cell adhesion are significantly enriched. Cell adhesion is essential for embryo implantation [44]. When examining the differentially expressed gene profile of blastocysts between dormant and activated blastocysts, the major functional categories of altered genes include the cell cycle, cell signaling, and energy metabolic pathways [4]. During delayed implantation, the uterus remains quiescent and the blastocysts become dormant [45]. Further study on these differentially expressed genes in mouse uterus should be beneficial for better understanding the mechanism on delayed implantation.
Estrogen is essential for on-time uterine receptivity and blastocyst activation in mice [3]. Ovariectomy in the morning of day 4 will lead to blastocyst dormancy. In the model of delayed implantation, estrogen is used to activate the dormant blastocyst and initiate embryo implantation [6,46]. Therefore, we compared the up-regulated genes in our study with the estrogen-stimulated genes. In mouse uterus, 102 genes are up-regulated at least 2 folds after estrogen treatment for 24 h [47]. In our study, 268 genes are up-regulated by least 2 folds in activation group compared to delayed implantation. Among 268 up-regulated genes, 16 genes are also shown to be differentially regulated in mouse uterus after estrogen treatment [47]. This suggests that at least a part of the upregulated genes in activation group is regulated by estrogen rather than embryos. The integrative analysis of mRNA and miRNAs In our study, we compared miRNA expression with their corresponding host genes. The expression trend between miRNAs and their host genes is highly correlated for both delayed implantation and activation, suggesting that host genes and miRNAs are co-transcribed. Most known miRNA genes have the same type of promoters as protein-coding genes have [48]. Perfect correlations are also found between the expression profiles of the intronic miRNAs and their host genes [49]. The strong correlation between miRNAs and their host genes indicated that they are derived from the same precursor genes and might be under the control of the same promoter. Therefore, it is possible to predict the expression of their embedded miRNAs through large scale analysis of the miRNA host genes.
Right now the information on large scale proteomic analysis of mouse uterus during embryo implantation is still lacking. However, it is shown that miRNAs could down-regulate some of their targets not only at the translational but also at the transcriptional level [13]. Therefore, it is possible to use the comparative expression analysis of miRNAs and mRNAs to identify mRNA targets of miRNAs. In our SAGE-Illumina sequencing data, there are 52.4% down-regulated genes and 47.6% up-regulated genes. Among the target genes of the upregulated miRNAs predicted by either TargetScan or PITA, 63.93% of the target genes predicted is really down-regulated, suggesting that the target genes of significantly up-regulated miRNAs during activation are significantly enriched. In our study, the enrichment level of down-regulated genes is not very high and the up-regulated genes are not significantly enriched, suggesting that some genes might not be regulated at the transcriptional, but at translational level. In MCF7 cells, the target genes of 14 downregulated miRNAs are significantly enriched, indicating the strong inverse correlation between miRNA and target gene expression might be essential in gene regulation during the acquisition of fulvestrant resistance [50]. In Drosophila, miRNAs and genes encoding predicted miRNA targets are expressed in a largely mutually exclusive manner [51]. Recent evidences suggest the involvement of miRNAs in tuning the expression of target genes to physiologically relevant levels [52]. Therefore, the differentially expressed miRNAs during embryo implantation may be essential through regulating the mRNA level of their corresponding target genes.
In conclusion, many differentially expressed miRNAs and mRNAs were identified in mouse uterus under delayed implantation and activation in our study. For miRNAs, there are 20 miRNAs up-regulated and 42 miRNAs down-regulated at least 1.2 folds under activation compared to delayed implantation. For mRNAs, there are 268 genes up-regulated and 295 genes downregulated at least 2 folds under activation compared to delayed implantation. Both miRNA and mRNA expression patterns are closely related and the target genes of up-regulated miRNAs are significantly enriched. There is a higher percentage of miRNA editing at positions 4 and 5 of mature miRNAs under delayed implantation than that under activation. Our data will shed light on further study of mouse embryo implantation.

Animal treatments
Mature mice (Kunming White outbred strain) were maintained in a controlled environment with a 14-h light/10-h dark cycle. All animal procedures were approved by the Institutional Animal Care and Use Committee of Xiamen University (XMUEA-0080). Female mice were mated with fertile males of the same strain to induce pregnancy (day 1 is the day of vaginal plug). To induce delayed implantation, pregnant mice were ovariectomized under ether anesthesia at 08:30-09:00 h on day 4 of pregnancy. Delayed implantation was maintained from days 5-7 by injecting progesterone (1 mg/mouse, Sigma) in the morning. Estradiol-17b (25 ng/mouse, Sigma) was given to progesterone-primed delayed implantation mice to initiate implantation on day 7 of pregnancy. The mice were sacrificed to collect uteri 24 h after estrogen treatment for activation group. Delayed implantation was confirmed by flushing the blastocysts from one horn of the uterus. The implantation sites of activated uterus were identified through intravenous injection of 0.1 ml of 1% Chicago blue. Uterine tissues were collected from at least 20 mice undergoing delayed implantation and activation, respectively. Equal amounts of uterine tissues from delayed implantation and activation groups were subject to the following Illumina sequencing analysis.
The 59 end of the read was treated as 59 nucleotide of the small RNA, and the 39 end of the small RNA was determined by the 39 most perfect match to the first 8 nt of the 39 adaptor. After the reads without a perfect 8-nt adaptor match were deleted, the remaining reads were retrieved from libraries of delayed and activated uteri, respectively. The adaptor-free reads were aligned to mouse genome mm9 (http://hgdownload.cse.ucsc.edu/downloads.html) by Bowtie alignment tool (http://bowtie-bio.sourceforge. net/index.shtml) [53]. The positions of all miRNA genes were also downloaded from miRBase and used for this positional annotation (http://microrna.sanger.ac.uk/sequences/ftp.shtml). Small RNA annotations other than miRNA were downloaded from piRNABank (http://pirnabank.ibab.ac.in/) [54] and fRNAdb (http://www. ncrna.org/frnadb/) [55]. Each of these reads was classified as known miRNA, piRNA, tRNA, rRNA, snRNA, snoRNA, mRNA, genomic sequence or unknown sequences.
The comparison of miRNA copies between our two deep sequencing libraries was performed according to the Z-test algorithm as described previously [56]. For a miRNA, n 1 and n 2 are the read number of this miRNA andN 1 and N 2 are the total number of reads in each library, respectively. The z-statistic is calculated according to the following formula: in which p 1~n1 =N 1 , p 2~n2 =N 2 and p 0~n1 zn 2 ð Þ = N 1 zN 2 ð Þ . Two-sided p-value was calculated from z-statistic and followed by Bonferroni multiple test correction. The same method was applied for digital gene expression data in this study.
Detection of microRNA editing miRNAs with known SNPs [18] in mature sequences were excluded for further analysis. In order to avoid cross-mapping of small RNA reads, a rough alignment was performed with Bowtie software. The potential read/known microRNA sequence pairs were then aligned by Needleman-Wunsch dynamic programming algorithm. The penalty scores for perfect match, mismatch, gap opening and gap extension were set for 1, 21, 22 and 21, respectively. Only the alignment with the highest score for each read was kept for further analysis.
Statistical analysis was performed as previous described [10]. The test is based on the null hypothesis where all positions behave the same with respect to base modification. The editing rate of each nucleotide position was calculated. Then the editing rate was transformed into standardized score (Z-score). The median Zscore across all positions was set to 0. A p-value was assigned to each Z-score according to normal distribution. A significant editing which is much higher than background noise (sequencing error) was considered if a p value is less or equal to 0.01.

SAGE-Illumina sequencing analysis of mRNAs
Total RNAs of delayed and activated uteri were extracted by TRIzol. Total RNA quality from both delayed and activation groups was comparable based on analysis with a Bioanalyzer 2100 (Agilent). The mRNA library for sequencing was prepared using Gene Expression Sample Prep Kit (Illumina) according to the manufacturer's instructions. Briefly, mRNAs were isolated through binding to a magnetic oligo(dT) beads. First strand cDNA was synthesized using the binding mRNA as a template, and followed by the synthesis of the second strand of cDNA. After cDNAs were digested with DpnII, the double stranded cDNA fragments attached to the oligo(dT) beads were collected and ligated to a GEX DpnII adapter 1 at the site of DpnII cleavage. The sequence for MmeI was included in GEX DpnII adapter 1. After purification, products coupled with oligo(dT) beads were digested with MmeI to create 16 bp tags, followed by ligating to a GEX adapter 2 at the site of MmeI cleavage. After PCR amplification, the purified DNA fragments were used directly for sequencing using the Illumina Cluster Station in Shenzhen Huada Gene Sci-Tech Company.
After sequencing, 16 bp tags were extracted from SAGE-Illumina sequencing libraries by in-house perl scripts. Then the tags were mapped to Unigene build 21. SAGEmap algorithm was used for tag-to-gene mapping with slightly modifications. Briefly, Unigene build 21 data were downloaded from NCBI. For each sequence in the Unigene database, 16 base tags adjacent to the 39most anchor enzyme DpnII site (GATC) were extracted.
The sequences were divided into 4 types: (a) mRNA: tags from GenBank submission transcripts which have poly(A) tails and/or signals; (b) High-throughput sequencing: tags from high throughput sequencing transcripts which have either poly(A) tails and/or signals; (c) ESTs with poly(A) tails: tags from EST sequences which have poly(A) tails and/or signals in the same orientation as the tags. (d) Sequences without poly(A) tails: all sequences of which no polyadenylation signal or tail was found. For each tag-UniGene pair, a reliable score was calculated as: Score~number of sequence types ð Þ |1000000 z reliable source score ð Þ |1000 z number of sequences ð Þ The number of sequence types could range from 1 to 4 depending on whether the tag-UniGene pair was defined as mRNA, mRNA sequences from high-throughput sequencing, ESTs with poly(A) tails, and/or sequences without poly(A) tails. The reliable source score = (number of mRNA + number of mRNA sequences from high-throughput sequencing) +0.56 number of ESTs with poly (A) tails. The higher the score of a tag-UniGene pair, the more reliable the mapping is considered. When more than two UniGene clusters were assigned to a certain tag, only the two UniGene clusters with the highest scores were chosen.
qRT-PCR analysis for mature miRNAs, precursor miRNAs and mRNAs Total RNAs were extracted from delayed or activation mouse uteri with TRIzol reagent, digested with DNase I and reversetranscribed into cDNA with PrimeScript TM RT reagent Kit (TaKara, Dalian, China). For mature miRNAs, there were 0.1 mg total RNAs, 2 ml PrimeScript Buffer, 25 nM stem-loop RT primers, and 0.5 ml PrimeScript RT Enzyme Mix I in 10 ml volume. Reverse transcription was performed at 16uC for 30 min followed by 60 cycles of 30uC for 30 sec, 42uC for 30 sec and 50uC for 1 sec, and ended at 85uC for 5 min. qRT-PCR was performed using a SYBR Premix Ex TaqTM kit on the Rotor-Gene 3000A system at 95uC for 10 sec, followed by 40 cycles of 95uC for 5 sec, 60uC for 5 sec and 72uC for 8 sec. For precursor miRNAs, 50 nM special reverse primers were used as RT primers. For mRNA quantification, 2.5 mM Oligo dT Primer and 5.0 mM Random 6 mers were used for reverse transcription and qRT-PCR was performed at 95uC for 10 sec, followed by 40 cycles of 95uC for 5 sec and 60uC for 34 sec. Primers used for qRT-PCR were listed in Table S5. Mouse rPL7 gene was amplified as a reference gene for normalization.

Dual-luciferase activity assay
The 39-UTR segment of each mouse gene predicted as target genes of let-7a-5C (C at the 5th position of let-7a) or let-7a (the reference sequence of let-7a, G at the 5th position) was amplified by PCR from mouse cDNAs and inserted into the downstream of luciferase reporter gene in the pGL3 control vector. Primers used for amplifying 39-UTR of each mouse gene were listed in Table  S5. The plasmid pRL-TK containing renilla luciferase was cotransfected for data normalization. Transfection and dual luciferase analysis was performed as described previously [8].

Primary culture of uterine stromal cells
Uterine stromal cells were cultured as described previously [8]. Uterine horns from mice on day 4 of pregnancy were cleaned of fat tissues, slit longitudinally, and washed thoroughly in Hanks' balanced salt solution (HBSS, Sigma) without Ca 2+ /Mg 2+ and phenol red. Tissues were then placed in 5 ml of fresh medium (HBSS with antibiotics) containing 6 mg/ml dispase (Gibco-BRL) and 10 mg/ml trypsin (Sigma), and incubated for 1 h at 4uC, 1 h at room temperature, and then 10 min at 37uC, respectively. Following the digestion, tissues were shaken several times to dislodge the sheet of luminal epithelial cells. The remaining tissues were washed three times in fresh medium and digested in HBSS containing 0.15 mg/ml collagenase (Gibco-BRL) at 37uC for 30 min. Following digestion and shaking, the digested cells were passed through a 70-mm filter to get rid of residual epithelial sheets and centrifuged. The cells were plated in 24-well culture plates and cultured in DMEM containing 10% FBS.

Western blot
Cultured cells or uterine tissues were collected in lysis buffer (50 mM Tris-HCl, pH 7.5, 150 mM NaCl, 1% Triton X-100, and 0.25% sodium deoxycholate) and briefly sonicated to shear DNA and reduce sample viscosity. Protein concentration was measured by BCA Reagent kit (Applygen, Beijing, China). Samples were run on a 10% PAGE gel and transferred onto nitrocellulose membranes. After blocked in 5% nonfat dry milk in TPBS (0.1% Tween 20 in PBS) for 1 h, membranes were incubated with monoclonal anti-human SMAD7 Antibody (R & D Systems) overnight at 4uC. After three washes in 5% milk/TPBS 10 min each, membranes were incubated in goat anti-mouse IgG conjugated with horseradish perioxidase for 1 h followed by two washes in 5% nonfat milk in TPBS, TPBS and PBS 5 min each, respectively. The signals were developed in ECL Chemiluminescent kit (Amersham Pharmacia Biotech, Arlington Heights, IL).

In situ hybridization
In situ hybridization was performed as we previously described [6]. Briefly, frozen uterine sections were hybridized with the digoxigenin-labeled Smad7 RNA probes. All of the sections were counterstained with 1% methyl green.