Quantification of small non-coding RNAs across human preimplantation developmentTo comprehensively evaluate the sncRNAome during preimplantation development and lineage segregation, we profiled single cells from human preimplantation embryos spanning embryonic days (E) 3 to 7 (E3-E7) (Fig. 1a). Female patient age ranged from 26-45, with an average of 36.6 years. Following quality control, we analysed 972 cells from 69 embryos with an average of 65.2k unique molecular identifiers (UMIs) annotated to known sncRNAs from available databases (see methods) at E3, decreasing to 29.3k at E6/7 (Fig. 1b, c, Supplementary Data 1), which was largely attributed to a decrease in piRNAs and tRNAs (Fig. 1c). The main biotypes of sncRNAs and their proportion were identified including miRNAs (9.62%), snoRNAs (23.8%), rRNA (18.6%), piRNAs (17.2%) and tRNAs (29.3%) (Fig. 1c). We observed that while piRNAs and tRNAs decreased with development, suggesting the majority are parentally inherited, miRNAs, snoRNAs, and rRNA fragments gradually increased (Fig. 1c) – a trend that has been reported to begin in human zygotes16. The length distributions from each class show typical patterns, with primary peaks at 22 nt for miRNA, 74 nt for tRNA, and 30 nt for piRNA (Fig. 1d), further verifying the quality of data obtained.Fig. 1: Single-cell profiling of sncRNAs in human embryos.a Overview of the experimental workflow. Donated embryos were cultured from E3-E7, dissociated into single cells, and manually picked for Small-seq or Co-seq. The dotted arrow indicates where a subset of cells were used. Number of cells and embryos remaining following QC filtering were attached below (created with BioRender). b Total UMIs and (c) proportion of UMIs annotated to microRNA (miRNA), small nucleolar RNA (snoRNA), small nuclear RNA (snRNA), transfer RNA (tRNA), ribosomal RNA (rRNA), and piwi-interacting RNA (piRNA) by day of embryonic development. The number of cells for each day is as follows: E3: 20, E4: 106, E5: 304, E6: 394, E7: 148. The boxplot rectangles represent the first and third quartiles, a horizontal line inside the box indicates the median value. d Percentage of total annotated reads by sncRNA class of various sequence lengths. e Gene expression of sncRNA processing genes across embryonic development and lineage (inner cell mass (ICM) and trophectoderm (TE)). Error bars represent standard error of the mean (SEM) values. The number of cells for each day is as follows: E3: 61, E4: 144, E5: 245, E6: 278, E7: 361. Created in BioRender. Kwok, Y. (2023) BioRender.com/v41c600.We next aimed to better understand the change in expression patterns observed for the sncRNA biotypes by determining the expression of key components involved in the biogenesis of sncRNAs, as these have not been previously reported across preimplantation development. To do this, we leveraged our published embryonic data which mainly focused on protein-coding RNAs17. We observed an overall increase in the pri-miRNA processing machinery (DGCR8, DROSHA), with a relative decrease in piRNA-interacting proteins and maturation factors (PIWIL2, PIWIL3, HENMT1) with development (Fig. 1e). The general pattern in expression of these components fits with the decrease in piRNA observed with development and further supports their lack of de novo generation. Further, the gradual increase in miRNA micro-processing complex and DICER supports the increase in miRNA molecules observed with development in the embryo. The Argonaute proteins AGO1-4 form the base for the RNA-induced silencing complex (RISC) and appear to fluctuate by developmental time and lineage, but are expressed quite abundantly, suggesting a functional role of miRNAs associated with blastocyst formation.We also examined the genomic location of reads that do not overlap with the main types of small RNAs, by comparing the mapping positions with annotations for protein-coding regions, lncRNA, and repeat regions from Repbase18, as well as other annotations from GENCODE19 (predominantly pseudogenes) (Supplementary Fig. 1). The majority of these reads are primarily from the sense strand of annotated protein-coding regions (10.8%) and lncRNAs (9.8%), with read length peaks at 18 and 21 nt. Reads overlapping with repeat regions are mainly from annotated LTR retroposons (1.5% on the sense strand and 0.8% on the antisense strand), with a read length peak at 18 nt. It remains unclear whether these reads have a biological or technical origin.Identification of cell-type specific miRNA signaturesOur knowledge surrounding the role of sncRNAs is limited during mammalian preimplantation development, particularly in the human. We reasoned that if we could identify biotype signatures which correspond to cell type and determine their developmental dynamics, we may begin to unravel some additional mechanism(s) underlying blastocyst formation. Given the lack of known sncRNA markers for lineages in the human embryo, we used our previously developed method, Co-seq, to first accurately classify our cell types14. Co-seq relies on a split-cell experiment, followed by library amplification of sncRNA using Small-seq and RNA using Smart-seq2. Leveraging our previous single-cell RNA-sequencing (scRNA-seq) data, we integrated it with the current scRNA-seq from Co-seq to identify cell types observed in our Co-seq data (Fig. 2a, b)17. This allowed us to accurately classify the ICM and TE, in addition to the EPI, PE, polar and mural lineages (Fig. 2b, c). Using the Small-seq portion of each split cell, we were then able to overlay cell type with clusters observed and generate unique miRNAs for each of these lineages. We observed the greatest enrichment of miR-376c-3p and miR-376a-3p in ICM cells and miR-525-5p and miR-518b in TE cells (Fig. 2d). Although relatively little is known about the role of the miR-376 cluster in embryo development, the functions of the chromosome 19 miRNA cluster (C19MC) which include miR-525 and miR-518b, have been well established as markers of TE differentiation11.Fig. 2: Lineage identification and miRNA dynamics across embryo development.a Integration of Co-Seq transcriptome data with5 to establish cell identities. b Violin plots of marker gene expression in Co-Seq data. c UMAP of miRNA expression data from Co-Seq cells, coloured by cell identity and lineage. d Expression of marker miRNA miR-376c-3p (ICM), miR-376a-3p (ICM), miR−525 − 5p (TE), miR−518b (TE) in miRNA portion of Co-Seq data. e UMAP of miRNA expression in all cells passing QC, f coloured by lineage and (g) developmental time. h Heatmap of top significantly enriched miRNA expression in each lineage. i The proportion of miRNA markers in each miRNA family, split by lineage. EPI epiblast, ICM inner cell mass, PE primitive endoderm, TE trophectoderm, E3-E7 embryonic days 3–7, EB early blastocyst.To verify the cell types present in our comprehensive dataset and identify sncRNA signatures for each cell type, we integrated the lineage-specific miRNA signatures obtained from the Co-seq experiment (Fig. 2e). Taking into consideration embryonic day, we identified miRNA expression dynamics from E3 to E7, which we found corresponded to pre-lineage (8-cell, morula, and early blastocyst (EB; just prior to ICM and TE specification)), followed by transition into ICM and TE (Fig. 2f, g). A cluster of unclassified cells was found to consist of cells from all developmental timepoints (Supplementary Fig. 2a–c), which were excluded from downstream analysis as they did not appear to be developmentally related. Developmental time contributed to the largest changes in miRNA expression; with lineage specific miRNA emerging at E5, just after cavitation, corresponding to the first lineage split of ICM and TE, and increasing to E7 (Fig. 2f, g). In addition, we investigated the relationship between embryo morphology (stage) and miRNA expression, and observed that miRNAs clustered more congruently with embryonic day (in which we also considered morphological parameters) (Supplementary Fig. 2d, e). We next identified those unique miRNA markers that corresponded to each lineage, with the top 10 displayed (Fig. 2h, Supplementary Data 2). We observed an enrichment of miR-506 family members in pre-lineage cells, including: hsa-mir-508-3p, hsa-mir-509-3p, hsa-mir-513a-5p, and hsa-mir-514a-3p (Fig. 2h, i, Supplementary Data 2). This cluster is a developmentally-related locus with high expression in spermatogenesis20,21, however its function(s) in early embryogenesis are unknown. High levels of expression of the pluripotency-associated miR-182-5p and miR-183-5p were also detected in pre-lineage cells22. All of the top 10 enriched miRNAs in ICM cells, including miR-381-3p, miR-376c-3p and miR-299-3p, have been implicated in the maintenance of ground-state pluripotency in mouse ESCs and ICM (Fig. 2h)22,23. The majority of the TE-enriched miRNA we identified, including miR-518b, miR-519c-3p, and miR-27b-3p, are well-known trophoblast miRNA markers24,25. We now speculate that these miRNAs may be of particular importance given that they are now observed in the preimplantation TE and maintained during postimplantation. Perhaps, early dysregulation of these miRNAs during the preimplantation window may result in suboptimal placental formation and pathological conditions, providing possible targets for predicting embryo competence in clinic. We further investigated the marker miRNAs to determine their genomic origin. We observed that approximately ~60% of those identified in the TE belong to the miR-515 gene family, embedded in the C19MC (Fig. 2i). Similarly, the miR-154 family originating from a miRNA cluster on chromosome 14, is highly represented in the ICM lineage and contributes to ~20% of enriched miRNA.Given that human ESCs are derived from the preimplantation embryo, we next asked whether the expression of embryonic miRNA markers observed in the pre-lineage, ICM, and TE were similar to naïve and primed ESCs, as well as the cancer cell line, HEK293T, as a control (Supplementary Fig. 3a). We observed no obvious enrichment of the pre-lineage markers in any of the three cell lines. However, a clear enrichment of almost all of the ICM and TE markers was observed in naïve ESCs, many of which have been shown to be required for naive pluripotency in ESCs11. There were three notable exceptions in the ICM markers, miR-302a-5p, miR-302b-5p, and miR-367-3p, all of which were higher in primed ESCs. We also did the reverse analysis and examined the expression of miRNAs determined to be specific to naïve and primed hESCs (Supplementary Fig. 3b), which confirmed an enrichment of naïve pluripotency miRNAs in the ICM. Naïve hESCs are known to be representative of the E6/E7 human preimplantation embryo and these data further corroborate the close molecular relationship between naïve hESCs and blastocyst derived cells.The most striking finding from our enrichment analysis is the identification of two genomic hotspots associated with the lineage clusters: upregulation of a 123 kb region in TE cells located at Chr.19q13 (C19MC), and of a 180 kb region in ICM cells positioned at Chr.14q32 (C14MC; Fig. 3a, b). In the ICM, we observed 67% (14/21) of the marker miRNA originating from the sense strand of the C14MC cluster (Fig. 3a, b). We first detected C14MC expression in the early blastocyst, just prior to distinct ICM-TE lineage specification (Fig. 3c), and then remained highly enriched only in the ICM (Fig. 3c). The enrichment of the C14MC cluster exclusively in the ICM was surprising as it is reported in post-implantation trophoblast cells26, similar to the C12MC cluster which is the mouse counterpart (reviewed in ref. 27). Mir-381-3p, originating from the C14MC cluster, was the most highly expressed miRNA in our ICM cells and is one of the major drivers of ground-state pluripotency in mouse ESCs22. Perhaps, there is a switch in the functional role of C14MC that is dependent on the window of development and a ramping up of expression is observed in the trophoblast postimplantation. Of note, although highly ICM-enriched, the average C14MC miR expression was 19-fold lower compared to the miRNAs of the C19MC (Supplementary Data 3).Fig. 3: Enrichment and expression of the C19MC and C14MC in early embryos.a Circular expression plot of ICM-enriched miRNA and their genomic origin, with expanded visualisation of C14MC, and similarly for (b) TE-enriched miRNA. IGV expression plots of relative sncRNA expression in our dataset from the (c) C14MC and (d) C19MC loci between lineages. Additional tracks subtracting the TE from ICM expression on E5 and E6/7 demonstrates enrichment of C19 in TE, over basal expression levels in the ICM. ICM inner cell mass, TE trophectoderm, E3-E7 embryonic days 3–7, EB early blastocyst.C19MC is a maternally imprinted locus, which is upregulated in trophoblast and cell line derivatives, ESCs and germ cell maintenance11,28. In our TE cells, 72% (65/90) of the marker miRNA originated from the C19MC, all of which are located on the sense strand. Members of the C19MC were expressed at all developmental time points and lineages, though more highly expressed (22%) in the TE upon formation of the lineages, as evidenced in our enrichment analysis and TE-normalised ICM expression (Fig. 3d, Supplementary Data 3). We also observed that the magnitude of expression in the TE compared to ICM continues to diverge with time. Here, we demonstrate enrichment of C19MC in the preimplantation TE lineage, suggesting an early expression and perhaps a role in the regulation of trophoblast differentiation (Fig. 3d). Further studies are warranted to improve our understanding of the roles these clusters play in preimplantation development.Non-canonical sncRNAs in preimplantation developmentNext, we examined whether we could identify cell type-specific signatures for the other classes of sncRNAs. We observed minimal cell clustering when analysing expression of tRNA fragments, piRNA or snRNA by embryonic day or lineage (Supplementary Fig. 4a). This suggests that these biotypes undergo cell type-independent changes in expression level during the transition from 8-cell to differentiated blastocyst. In contrast, we were able to obtain distinct cell type signatures for snoRNAs expression (Fig. 3c, Supplementary Fig. 4a). SnoRNAs SNORD115-116 from the paternally expressed long non-coding RNA host gene SNHG14 were enriched at E3-E4, prior to lineage specification, while the snoRNAs from the paternally imprinted MEG8 gene, SNORD112-114, were enriched in the ICM (Fig. 3c, Supplementary Fig. 4b,c). Although limited data implicate snoRNAs in lineage formation, some previous data suggest a dynamic profile in human ESCs29 and during the morula to blastocyst transition in the bovine30. Interestingly, spliceosomal repression can effectively reprogramme mouse ESCs to a totipotent state, and SNORD116 has unique effects on the spliceosome31. Previous profiling of 14 different tissues and cell lines suggest SNORD114 has the highest expression in the ovary and is predicted to interact with and post-transcriptionally modify 18S rRNA32. The genomic co-localization of SNORD112-114 and the ICM-enriched miRNA suggest that C14MC may be generating both snoRNA and miRNA precursors upon emergence of the ICM lineage. This is the first report of C14MC expression in the ICM, pointing to possible important role(s) of this locus in the maintenance of pluripotency and ICM lineage formation.We next investigated tRNA expression dynamics as they have been previously shown to be developmentally regulated29,30. Although no obvious clustering by embryonic day or lineage was observed with UMAP analysis of aggregated tRNA count data (Supplementary Fig. 4a), we further separated tRNAs into 5’ halves (5’-tRFs), 3’ halves (3’-tRFs), and full length sequences and observed an increase in proportion of full length tRNAs from E3 to E7 with a corresponding decrease in 5’-tRFs and 3’-tRFs (Supplementary Fig. 5a and b). Surprisingly, the overall distribution of Arginine and Asparagine codon usage represented by 3’-tRFs drastically decreased and increased over developmental time, respectively (Supplementary Fig. 5c). A similar expression pattern of 5’-tRFs were observed, in which Asparagine codons increased and Glutamyl codons decreased across development. The fragment data is in contrast to the codon usage in full length tRNAs which appeared consistent from E3 to E7, and in other codons which have limited differences in 5’-tRFs and 3’-tRFs (Supplementary Fig. 5c). While little is known about the roles of tRNA fragments in development, tRFArg is able to donate its arginyl to Arginyltransferase, a highly conserved enzyme that causes mid-gestational embryonic lethality when knocked out in mice33,34. New mechanisms of gene regulation and post-translational modifications by tRFs and their binding partners are emerging, and additional work is needed to delineate their function in preimplantation embryogenesis (reviewed in ref. 7).Trajectory and pseudo-time analysis of miRNAs throughout preimplantationNext, we sought to evaluate the changes in miRNA expression at a higher resolution throughout E3 to E7 development. The largest changes in miRNA differential expression (DE) were observed between 8-cell to morula (96 DE miRNAs) and morula to EB (163 DE miRNAs), while the initial formation of the ICM and TE was accompanied by modest changes (5 and 27 significantly differentially expressed miRNAs respectively, with log2 fold change more than 0.1 and fdr less than 0.05; Fig. 4a; Supplementary Data 4). During the maturation of TE from E5 (early blastocyst) to E6/7 (late blastocyst), we observed 74 DE miRNA (Fig. 4a), which included upregulation of the majority of the C19MC cluster miRNAs (e.g. miR-522-3p, miR-516b-5p, miR-512-3p) and downregulation of many pluripotency-related miRNA (e.g. miR-372-5p, miR-373-5p, miR-302b/d-5p). In contrast, no DE miRNA were observed in the ICM between E5 to E6/7, suggesting that once expressed, ICM associated miRNAs are maintained. These data demonstrate the divergence in miRNA expression patterns whereby in the TE, increased expression of trophectoderm-related miRNA markers are observed while concordantly decreasing their pluripotency, meanwhile ICM cells maintain expression of a set of miRNA during lineage maturation to E6/7.Fig. 4: Pseudotime analysis of miRNA expression in early embryonic development.a Differential expression between developmental stages and lineages. b Developmental- and lineage-specific pseudotime trajectories with heatmaps of lineage trajectory-associated miRNA expression. Expression of several key lineage-associated miRNAs across pseudotime for (c) pre specification, (d) ICM, and (e) TE lineages. ICM: inner cell mass; TE trophectoderm, E3-E7 embryonic days 3–7, EB early blastocyst. The confidence interval (error bands, 95%) is indicated by bandwidth. The measure of centre and confidence intervals were calculated using the “loess” function with default parameters in R software.We then inferred the pseudotime of our cells to address how miRNA expression patterns dynamically progress during human embryo development. We were particularly interested in identifying miRNAs which could be potentially driving lineage specification or maintenance. We observed a clear trajectory from 8-cell to EB, at which point two branches formed, one for the ICM and one for the TE (Fig. 4b). Generally, miRNA expression dynamics across pseudotime could be grouped into three categories. 1. Highly expressed in the pre-lineage, with a decrease in expression across time, irrespective of lineage. 2. Drastic increase in expression following ICM formation and 3. An enrichment in TE following specification. Several miRNA that were found to be highly abundant in the pre-lineages included miR-184, miR-204-5p, miR-508-3p, and miR-513a-5p (Fig. 4c). Given the relative quiescence of the human genome prior to embryonic genome activation (EGA), these miRNAs may be maternally inherited. For example, miR-184 is known to be required for drosophila germline development35, and has been shown to be present in human ova along with miR-204-5p and miR-508-3p16,36. In contrast, development of the ICM trajectory is associated with an increase in several key miRNA (Fig. 4d). In mouse embryonic stem cells, miR-182 is required for self-renewal and maintenance of pluripotency37 and has been shown to be among the top 20 miRNA secreted from human blastocysts, independent of lineage38. miR-376 cluster members, miR-376a-3p and miR-376c-3p, are significantly increased along the ICM trajectory, and have previously been shown to target IGFR1, a key regulator of trophectoderm development39,40. Along the trophoblast trajectory, miR-27b-3p shows an exceptional rise, possibly participating in trophoblast differentiation through the suppression of pluripotency markers (Fig. 4e)41. Similarly, miR-519c-3p and miR-1323 are key members of C19MC which are highly expressed in trophoblast lineages42.With the abundance of C19MC and C14MC members among DE and trajectory-associated miRNAs, we aggregated expression across the entire cluster to better understand temporal changes in expression (Supplementary Fig. 6). miRNAs from the C14MC are relatively low until the EB window, during which their expression increases dramatically along the ICM trajectory but not along that of the TE (Supplementary Fig. 6a). Conversely, C19MC expression is already present at the 8-cell stage, rising to EB, and then diverges with higher expression in the TE vs ICM trajectory. These patterns are supported in ESCs, in which C14MC expression is drastically, and C19MC expression modestly, higher in naïve ESCs (Supplementary Fig. 6b)29. These data support a model in which the C14MC locus expression is activated in ICM, conferring a high-level of pluripotency, while expression of the C19MC locus is relatively high in both lineages but increases in the TE and decreases in the ICM with differentiation.Furthermore, to verify the embryo and lineage specific expression of selected miRNAs using a different approach, we performed RT-qPCR on embryos biopsies (Supplementary Fig. 7). While expression of all miRNAs was confirmed in the human embryo, only miRNAs determined to be significantly enriched in the ICM were enriched in the RT-qPCR results (miR-146b-5p and miR-182-5p). RT-qPCR lacked the sensitivity to validate miRNAs determined to be significantly enriched in the TE by Small-seq, possibly due to the low number of embryos. This is likely due to the minimal difference in expression of ‘TE enriched’ miRNAs compared to those identified in the ICM (Supplementary Fig. 7).Developmental and lineage-specific miRNA targeting analysisCanonically, miRNAs exert their influence through binding to the RNA-induced silencing complex (RISC), targeting the 3’-untranslated region of target gene RNA through base-pairing complementarity, and then either preventing translation (partial complementarity) or by endonucleolytic cleavage43. Next, we wanted to determine how developmentally-regulated miRNA could direct differentiation. We clustered our miRNAs into the 3 common expression dynamics across the ICM or TE trajectories, starting at 8-cell: C1=Low → High → Low, C2=Low → High, C3=High → Low (Fig. 5a, b). We restricted our targeting analysis to miRNAs following the C2 and C3 patterns, as we reasoned that they are the most likely to exert consistent directional effects on differentiation. By filtering for target genes with significant negative correlations with previous embryonic gene expression data17, we identified 646 and 1098 miRNA-target genes in ICM, and 1230 and 1031 in TE from the C2 and C3 cluster types, respectively (Fig. 5c–f, Supplementary Data 5). We were specifically interested in miRNA-target-gene pairs that overlapped in both cluster direction and trajectory, such that the expression patterns of the miR-target pair were similar during progression from 8-cell to establishment of ICM or TE. This resulted in 1898 miR-target pairs common between the TE and ICM trajectories in the C2 pattern, and 1845 in the C3 pattern (Fig. 5g). Several developmentally-related predicted interactions were observed including the targeting and potential downregulation of the totipotency-related, zygotic genome activation (ZGA) markers PRAMEF1 by miR-367-3p, DPPA3 by miR-18a-5a, and SLC34A2 by miR-373-5p in both trajectories17,44. These data further support the known role for a miRNA-mediated loss of totipotency during the transition from 8-cell to blastocyst. In zebrafish and xenopus embryos, zygotic miRNAs have been shown to participate first in translation repression followed by degradation of maternal RNAs during the ZGA45,46, in contrast to mouse models where miRNA function is dispensable prior to the blastocyst stage47. Our analysis revealed 188 genes potentially targeted and therefore downregulated by miR-373-5p, which follows a C2 pattern (steady increase) along both the ICM and TE trajectories (Supplementary Data 5). miR-373-5p, along with a host of other miRNAs which followed the C2 pattern (miR-302a/b-3p, miR-520b/c-3p, miR-372-3p, miR-367-3p and miR-373-3p), may drive differentiation and an exit from 8-cell totipotency through coordinated targeting and downregulation of the early embryonic gene programme (Fig. 5h).Fig. 5: Trajectory-related miRNA patterns and potential embryonic gene regulation.K-means clustering of trajectory-related miRNA expression patterns in (a) ICM and (b) TE. C1=Low → High → Low, C2=Low → High, C3=High → Low. Error bars represent standard deviation values calculated among the miRNAs belonging to each cluster, the number of miRNAs of each cluster is indicated above. c–f Heatmaps of miRNA and target gene expressions in TE and ICM trajectories in C2 and C3 patterns. g Overlap between miRNA-target pairs between ICM and TE trajectories within C2 (top) and C3 (bottom) patterns of expression. h Number of target genes for C2 miRNAs with significant target enrichment over background. ICM inner cell mass, TE trophectoderm, E3-E7 embryonic days 3–7, EB early blastocyst.To determine the potential influence of miRNAs on lineage formation, we identified ICM- and TE-enriched miRNA signatures (n = 31 and 19, respectively; Fig. 6a; Supplementary Data 6). Given the redundancy of miRNA targets, we ‘pre-filtered’ potential gene targets by focusing on the biologically relevant genes in human embryos17. As such, we leveraged our previously published single-cell RNA-seq dataset in the human and identified ICM- and TE- enriched genes (n = 843 and 659, respectively; Fig. 6b). Among the significant miRNAs, 21 and 14 were significantly negatively correlated (p value < 0.05) with the gene targets identified (n = 197 and 71, respectively; Fig. 6c, d; Supplementary Data 6)5. Based on the repressive nature of miRNAs, we speculated that highly enriched miRNAs in the ICM may participate in lineage specification through targeting of TE-related genes, and vice versa. Through screening of both predicted and validated gene targets, we generated miRNA-target networks for both lineages48. ICM-enriched miR-369-3p, miR-409-3p, miR-146b-3p, and miR-381-3p are predicted to target well known TE transcription factors, GATA3, GATA2, TEAD1, and TET2 (Fig. 6e; Supplementary Data 6)49,50,51. Other well-characterised TE markers potentially targeted by enriched miR-409-3p, miR-374a-3p, and miR-378a-3p in the ICM include, NR2F2, PEG10, and LRP2, respectively52,53. Although none of these predicted miRNA-target interactions have been validated in embryos or hESCs, the 3’UTR of NR2F2 has been shown to be targeted by miR-409-3p in both HEK293 and TZM-bl cell lines54,55. The enrichment of ICM miRNAs suggest a role in maintaining self-renewal and preventing differentiation to TE. Similarly, we observed TE-enriched miRNAs that have corresponding gene targets enriched in the ICM. These include miR-524-5p, targeting the well-known naïve pluripotency gene KLF4, and a particular enrichment of miRNAs targeting PE related genes such as miR27b-3p, miR-151a-5p, miR-524-5p, and miR-27b-3p which target PDGFRA, SOX17, HNF1B, and IL6ST, respectively17,56,57,58 (Fig. 6f). Of these, miR-151a-5p has been experimentally validated to reduce the expression of endoderm marker SOX17 by luciferase reporter assays in cancer cell lines59. Further, ICM-enriched lysine demethylase KDM7A is targeted in the TE by four different TE-enriched miRNA: miR-27b-3p, miR-23b-3p, miR-23a-3p, and miR-524-5p. Zygotic knockdown of KDM7A in porcine embryos reduced blastocyst formation by 69%, dysregulated H3K9 and H3K27 methylation, and resulted in a reduced numbers of ICM cells, indicating the importance of KDM7A on development of the ICM60. Collectively, these miRNAs may enforce TE identity.Fig. 6: miRNA targeting analysis and potential influence on lineage segregation.a Volcano plots of miRNA and (b) their target genes in the ICM vs TE in E5-E7 cells. c Heatmaps of ICM-enriched miRNA expression (left) between lineages from E5-E7 and the corresponding TE enriched target genes (right). d Similarly for TE-enriched miRNA and their ICM-enriched target genes. P values of correlation were calculated by one-sided Student’s t test. e miRNA-target gene network plot for ICM-enriched miRNAs and (f) TE-enriched miRNAs. For visualisation, the top 35 target genes of each miRNA (ordered by fold change) with at least 0.3 log2 Fold change are presented. Gene expression data reanalysed from5. ICM inner cell mass, TE trophectoderm, E3-E7 embryonic days 3–7, EB early blastocyst.Building on our findings of miRNA-mediated repression in lineage specification, we expanded our targeting analysis to explore the potential for miRNAs to also play a role in the upregulation of gene expression. Though we acknowledge that this phenomenon may be rare, studies suggest that miRNAs can also enhance gene expression under certain conditions61,62. With this in mind, we examined the possibility that lineage-specific miRNAs might not only suppress but also positively regulate target genes, thereby reinforcing lineage identity. Our analysis revealed 255 miRNA-target pairs within the ICM showing significant positive correlations, such as miR-493-5p with SOX2, miR-182-5p with PDGFRA, and miR-493-5p with MAPK1 (Supplementary Data 6). In the TE, we identified 163 positively correlated pairs, including miR-27b-3p with GATA3 and GATA2, as well as miR-524-5p with NR2F2 (Supplementary Data 6). Therefore, it’s possible that miRNAs play a dual role in regulating lineage-specific gene expression, highlighting pairs for future functional validation to elucidate their roles in reinforcing lineage specification or maintenance.Investigating miRNA function in mammalian preimplantation embryosWe next wanted to determine the functional role of miR-381-3p, as it was highly expressed in the ICM cells at the blastocyst stage (Figs. 2h, 4d), suggesting it may have an important function in ICM-TE specification or pluripotency. In accordance with this, miR-381-3p has been shown to be highly expressed and have important roles in both mouse and human naive pluripotency23,63. Given the limited access to 8-cell human embryos and the conserved expression of miR-381-3p in mouse and human23,63, we chose to use the mouse model.First, we confirmed that mouse embryos passively internalised miRNA mimics (Invitrogen™) (ThermoFisher Scientific), by spiking in a scrambled negative control mimic conjugated to Alexa Fluor 488, as others have also demonstrated64 (Supplementary Fig. 8a). We next cultured mouse embryos with 10 µM miR-381-3p mimic, which prevented developmental progression to blastocyst stage, as expected given its enriched expression in the ICM and predicted targets related to TE formation (Fig. 7a). As we were interested in observing the role specifically in the ICM, we wanted to see if we could obtain some embryos that developed to the blastocyst stage by lowering the abundance of miRNA mimic. Indeed, with 5 µM miR-381-3p mimic, 45.8% of embryos were able to reach the blastocyst stage by E3.5 compared to controls (Fig. 7a).Fig. 7: The role of miR-381-3p during mouse preimplantation development.a Percentage of total embryos treated with scrambled negative control or miR-381-3p mimic that developed into blastocysts at embryonic day (E)3.5. Statistical analysis was conducted using the Chi-Square Test. Values on the bar graphs correspond to: n embryos that developed to blastocyst/total n embryos in the group. b Principal component analysis (PCA) plot based on scRNA-seq gene expression data from mouse embryos in control and miR-381-3p mimic-treated cells, where colour represents cell identity and shape represents the different treatments. c Ternary plot showing the average expression of DEGs in 8-cell, inner cell mass (ICM), and trophectoderm (TE) cells in control and mimic-treated embryos. DEGs were detected in control embryonic 8-cell, ICM, and TE cells, and further classified by the lineage with the highest gene expression. The vertices of the grey triangle inside the plot indicate the mean expression proportion in the three lineages for each type of DEG. d Violin plot shown selective marker expression in control and mimic-treated embryos. e Gene set enrichment analysis (GSEA) plot showing enrichment scores for genes with at least a 0.1 log2(fold change) difference related to DEGs determined between ICM vs. TE control, mimic-treated TE vs. control TE, and mimic-treated ICM vs. control ICM. P values from the one-sided Kolmogorov–Smirnov test are indicated in the figure heading. f Analysis of cell lineage allocation in early (32-cell) blastocysts treated with scrambled control versus 5 μM miR-381-3p mimic, where each data point represents the proportion of total embryonic cells belonging to the ICM or TE in a single embryo (n = 20 control embryos, pink and n = 19 mimic embryos, blue). Lineage was assigned based on the exclusive expression of known lineage markers Sox2 (ICM) and Cdx2 (TE) using immunofluorescence and confocal microscopy. Statistical analysis was conducted using a two-sided Student’s t-test, where values are reported as means ± standard deviation. g Average nuclear fluorescence intensity (protein expression, measured in arbitrary fluorescence units (AFUs)) of Sox2 (in ICM cells), Cdx2 (in TE cells), and Tead1 (in all cells) in early blastocysts treated with scrambled control versus 5 μM miR-381-3p mimic. Each data point represents the average intensity of the cells in one embryo. Statistical analysis was conducted using the Mann-Whitney or a two sided Student’s t test where appropriate. Values are reported as means ± standard deviation. h Representative immunofluorescence and brightfield (BF) images of control and mimic-treated blastocysts analysed in (f, g). Staining/imaging was performed in 5 separate batches on n = 11–17 control and n = 9–19 mimic-treated embryos. DAPI = nuclear stain. Comp = composite image. BF Brightfield. ns = p > 0.05, *** = p < 0.001, **** = p < 0.0001.To obtain a comprehensive overview of the impact of miR-381-3p mimic treatment on mRNA content and gene expression, we sequenced the transcriptomes of mouse embryos at the 8-cell and 32-cell stages from both control and miR-381-3p mimic-treated groups for 10 µM and 5 µM, respectively. After dimensional reduction using the top 2000 variable genes, we successfully identified the lineages for mouse embryo cells, confirmed by the expression of lineage markers such as Spp1, Upp1, Cdx2, and Dppa1 (Supplementary Fig. 8d, e). Although cells from control and treated groups were clustered together in UMAP space (calculated using the top 25 principal components), we found that PC_2 and PC_3 represent gene expression variance from lineages and treatment, respectively (Fig. 7b). Further, we observed that ICM and TE lineage-specific genes were less distinguishable in treated groups (Fig. 7c, d), suggesting that the miR-381-3p mimic (5 µM) was modifying the transcriptional profiles of cells at the 32-cell blastocyst stage to resemble an earlier developmental time point. Considering the generally lower expression of ICM and TE related genes in 32-cell embryos, we then performed GSEA for the DEGs between ICM and TE in 32-cell embryos to more broadly examine the role of miR-381-3p at this specific stage (Fig. 7e and Supplementary Data 7). We observed that the expression level of genes enriched in the ICM (determined by comparing ICM vs TE control), were decreased in the ICM cells derived from miR-381-3p treated embryos. In the TE cells derived from miR-381-3p treated embryos, the expression level of genes enriched in the TE (determined by comparing ICM vs TE control) decreased while the ICM enriched genes increased (Fig. 7e). This suggests that miR-381-3p is modifying the global gene signatures of the embryonic cells in both lineages making them less distinguishable. The regulation direction of DEGs in treated TE cells aligned with miR-381’s role as an ICM-specific miRNA. However, it is puzzling to observe that DEGs highly expressed in ICM lineage were also underexpressed in the miR-381-3p mimic-treated group.To examine the regulatory role of miR-381-3p at the protein level, including the potential impact on cell lineage allocation and the expression of ICM-TE lineage markers, we analysed surviving blastocysts treated with the 5 µM dose of the mimic using immunofluorescence. Treatment with 5 μM miR-381-3p mimic did not shift the proportion of ICM (Sox2 + ) and TE (Cdx2 + ) cells (Fig. 7f) nor the ICM:TE Ratio (Supplementary Fig. 8c) compared to controls. However, the expression of the lineage markers Sox2, Cdx2, and the computationally predicted target using TargetScan and miRDB, Tead1, significantly decreased in mimic-treated blastocysts (Fig. 7g, h and Supplementary Data 7). Notably, this was not due to a decrease in the total cell count of these early blastocysts (Supplementary Fig. 8b). Given that Sox2 and Cdx2 protein expression is lower in mouse embryos at earlier developmental timepoints, this suggests that although these embryos have reached cavitation and formed blastocysts, they are more representative of earlier timepoints in terms of molecular expression, in corroboration with our scRNA-seq analysis above. Together, when considering both the low and high dose of miR-381-3p mimic, our experiments suggest that this miRNA may play an important role in maintenance of ICM-like phenotype and/or prevent the TE programme.Novel embryonic miRNAsGiven that sncRNA expression in human preimplantation is relatively unexplored, we speculated that our dataset may provide insight into the presence of novel, highly-context specific miRNA. We identified six novel miRNAs, two of which passed strict abundance, sequence, and loci characteristics utilising miRDeep2 (Supplementary Data 8)65. These novel miRNA share a seed sequence with miRNAs from mouse and C. elegans66, suggesting they are evolutionarily conserved. Furthermore, these two miRNA candidates met additional criteria for annotating new miRNA67, and have been reported previously16,68, providing further confidence in our detection. The predicted structures and expression levels are presented in Supplementary Figs. 9–13.The previous identification of novel_miRNA_4_chr4_29031 in pre-EGA human oocytes and embryos16, and our observation that it decreases from 8-cell to blastocyst (Supplementary Fig. 9), suggests that it may be maternally inherited. Sequence abundance and structural analysis predicted the 3’ end of this novel miRNA to be the mature form (Supplementary Fig. 9). Next, we applied the same method as with the DE miRNA targeting and identified potential gene targets of this novel miRNA which display significantly inverse expression patterns (Supplementary Data 8). Significantly downregulated targets of novel_miRNA_4_chr4_29031 included KMT2A, a histone methyltransferase and SMARCD2, a chromatin remodeler. Previous knockdown experiments of KMT2A by miRNA targeting in mouse zygotes reduced H3K4 trimethylation and blastocyst development69, leading to the hypothesis that the observed decrease in novel_miRNA_4_chr4_29031 from E3 to E7 may permit translation of KMT2A. Similarly, SMARCD2 is a chromatin remodeler that is upregulated during preimplantation development70. We now speculate that novel_miRNA_4_chr4_29031 contributes to epigenetic reprogramming from zygote to blastocyst. Future functional studies of these novel miRNAs and their embryonic gene targets are warranted.miRNA isoform analysisSeveral studies have described the presence of miRNA variations in length and/or sequence due to physiological events occurring in vivo, termed isomiRs (reviewed in71). These variations include elongations/trimming of either end of the mature sequence and non-templated (post-transcriptional) additions to the 3’ end72. Therefore, we asked if embryonic cells exhibited any consistent modifications across development or within specific isomiRs. First, we observed that the overall proportion of miRNAs without length variations was 39.2% at E3, consistently rising to 43.5% at E7 (Supplementary Fig. 14a). The most abundant length variations included 3’ elongations, 3’ trimming, and 5’ trimming. The 21 isomiRs with the highest abundance were clustered by expression from E3 to E7 (Supplementary Fig. 14b). Several of the isomiRs displayed interesting dynamics across development. The 3’ trimmed isoform of miR-23a-3p is highly expressed at E3, but decreased drastically across development, meanwhile its 3’ elongated isoform consistently increases starting at E5. miR-23a-3p has been shown to display tissue-specific isoform expression and strong evolutionary selection73, while isomiRs as a whole have distinct regulatory activities from their canonical counterparts74. The dynamic regulation of miR-23a-3p isomiRs, among others, points to potential developmental-specific roles which are adjacent or complementary to their canonical sequences.We also asked if isomiRs with non-templated additions (NTAs) were present and variable in our dataset. NTAs were less abundant than length variations (average of 13.9% of UMI), but still followed a decreasing trend from E3 to E7 (Supplementary Fig. 14c). The most common NTA was the addition of one to three 3’ adenosines, which has been shown to have a stabilising or destabilising effect on miRNA, depending on the cellular context9,72. In our data, 8 of the top 10 isomiRs with the NTA of adenosine decreased from E3 to E7 (Supplementary Fig. 14d). In contrast, two of the most influential miRNA based on our trajectory analysis, miR-367-3p and miR-302c-5p, were the two miRNA with an increasing abundance of 3’ adenosine additions from E3 to E7. Our current knowledge pertaining to the role of isomiRs during embryo development, irrespective of species, is limited and warrants further exploration.
