miRNA profiling in bovine embryos
In the present study, we generated approximately 350GB raw data to identify the population of miRNAs present in cloned and IVF embryos at different developmental stages (2-cell, 8-cell and blastocysts) of buffalo. For all the samples, 79.67 to 85.40% of the total reads got aligned against the reference genome of Bos taurus, UMD 3.1.1 (Supplementary figure 1). The alignment statistics of the three biological replicates each of cloned and IVF embryos are given in Table 1. The number of known and novel miRNAs varied among the three biological replicates of cloned and IVF embryos at each stage (Figure 1). The raw data was normalized and the normalized signal values were used for subsequent data analysis. The Box-whisker plot showing the distribution of normalized signal values in all the replicates of cloned and IVF embryos at different stages are presented in (Supplementary figure 2). The quality of the sequencing data generated was analyzed by Principal Component Analysis (PCA), which depicted that the 2-cell 8-cell and blastocyst stage embryos were overall, grouped according to their origins (Supplementary figure 3). The cloned replicates were clustered together in one group while the IVF replicates were together in another group. The reads obtained were aligned on different chromosomes and the maximum number of reads was found to map on chromosome number 25 for all the three stages (Figure 2). Heat map consisting of all the three biological replicates of 2-cell and 8-cell and blastocysts stage of cloned and IVF were generated using Euclidean distance method of differentially expressed miRNAs with minimum fold change of 2 with p<0.05 (Figure 3).
Identification of differentially expressed miRNAs in preimplantation embryos
A total of 244, 247, 252 miRNAs were found to be differentially expressed at 2-cell, 8-cell and blastocyst in cloned relative to IVF embryos respectively. In 2-cell stage, out of 244 miRNAs, 186 were commonly expressed in both cloned and IVF 2-cell stage embryos, 35 miRNAs were expressed exclusively in cloned, whereas and 23 were expressed exclusively in IVF embryos. Similarly, at 8-cell stage 196 miRNAs were commonly expressed in both cloned and IVF 8-cell stage embryos, 27 miRNAs were expressed exclusively in cloned embryos whereas, 24 were expressed, exclusively in IVF embryos. Also, at blastocysts stage 195 were commonly expressed in both cloned and IVF blastocysts, 30 miRNAs were expressed exclusively in cloned whereas, 27 were expressed exclusively, in IVF embryos (Figure 4). The top 5 miRNAs up-regulated in cloned relative to IVF embryos at 2-cell stage were bta-mir-136, bta-mir-1187, bta-mir-431, bta-mir-494 and bta-mir-2474 whereas, the top 5 miRNAs down- regulated were, bta-mir-141, bta-mir-451, bta-mir-340, bta-mir-205 and bta-mir-760. At 8-cell stage cloned embryos, bta-mir-133a-2, bta-mir-431, bta-mir-202, bta-mir-1973 and bta-mir-282 were found to be among the top 5 upregulated miRNAs whereas, bta-mir-139, bta-mir-214, bta-mir-382, bta-mir-135 and bta-mir-2435 were the top five down-regulated miRNAs relative to IVF embryos. Similarly, in cloned blastocysts bta-mir-1949, bta-mir-1692, bta-mir-431, bta-mir-3956 and bta-mir-664a were the top 5 upregulated miRNAs whereas, were bta-mir-1973, bta-mir-1839, bta-mir-188, bta-mir-126 and bta-mir-592 were the top five down-regulated miRNAs. Top 30 differentially expressed miRNAs with fold change ≥2 at all three developmental stages of cloned embryos with respect to IVF counterparts are given in (Supplementary table 1, 2, & 3). In all the three comparisons miRNAs unique to cloned group were downregulated whereas, miRNAs unique to IVF groups were upregulated. MA and volcano plot were also used to visualize the distribution pattern of differentially expressed miRNAs in cloned embryos with respect to IVF preimplantation embryos. Different type of distribution patterns were observed in all the developmental stages of embryo (Figure 5).
In the present study, differential expression was also determined in terms of FC values for all the three stages of embryo development. The total number of differentially expressed miRNAs at different FC values (≥2 to <3-, ≥3 to <5- and ≥5-folds) in 2-cell, 8-cell and blastocysts stage is given in (Supplementary figure 4) Further, at different FC value, the number of up and down-regulated miRNAs in cloned relative to IVF embryos for all three developmental stages are shown in (Supplementary figure 5), respectively. At FC ≥2 (p<0.05) 47 miRNAs were found to differentially expressed in cloned relative to IVF 2-cell stage embryos out of which 31 were up-regulated and 16 were down-regulated. Whereas, in 8-cell stage 33 miRNAs were found to be expressed at significant level (p≤0.05) at FC ≥ 2 20 were up-regulated whereas 13 were downregulated. In case of cloned blastocyst relative to IVF 30 miRNAs were found to differentially expressed out of which 16 were upregulated and 14 were downregulated at FC ≥2 (p<0.05) (Figure 6).
Identification of commonly and uniquely expressed miRNAs
In the present study we also investigated the commonly and uniquely expressed miRNAs in cloned and IVF embryos at different FC level. Venn diagrams were used to visualize the miRNAs being expressed commonly between cloned and IVF 2-cell, 8-cell and blastocysts stage embryos and the ones uniquely expressed in the two groups respectively (Supplementary figure 6). A minimum cut off of 2-folds change revealed that 66 miRNAs were expressed differentially out of which 18 were unique to cloned embryos, 23 were unique to IVF embryos and 25 were expressed in both types of embryos. Similarly in 8-cell stage at minimum cut off of 2-folds change revealed that 76 miRNAs were expressed differentially out of which 25 were unique to the cloned embryos, 24 were unique to IVF embryos and 27 were expressed in both types of embryos. Also in case of blastocyst, at FC ≥2 73 miRNAs were expressed differentially out of which 30 were unique to the cloned blastocysts, 22 were unique to IVF blastocysts and 21 were expressed in both types of blastocysts. From among the commonly expressed miRNAs, those which were up- and down-regulated at different fold change level in cloned relative to IVF embryos are presented in (Supplementary figure 7).
At FC ≥2 and significance level of p<0.05, a total of 20 miRNAs were found to be commonly expressed between the two types of embryos out of which 4 were down- regulated and 16 were up-regulated in cloned relative to IVF 2-cell stage embryos. Whereas, in 8-cell stage out of 20 commonly expressed miRNAs in both the type of embryos 8 were down-regulated and 12 were up-regulated in cloned relative to IVF embryos. A total of 30 miRNAs were found to be commonly expressed between the two types of blastocysts at a significance level of p<0.05, out of which 14 were down-regulated and 16 were up-regulated in cloned relative to IVF blastocysts (Figure 7).
Gene ontology and KEGG analysis of cloned vs IVF preimplantation embryos
Using panther classification system, Gene ontology (GO) analysis was performed in order to inspect the biological significance, detailed annotation of gene function, biological process and cellular distribution of the targets of miRNAs expressed differentially between cloned and IVF embryos. Using GO terms, NGS results were summarized to provide insights into the changes in miRNA expression between cloned and IVF embryos at the 2-cell, 8-cell and blastocyst stage. In the present study, the targets of miRNAs expressed differentially at fold change ≥2 in cloned relative to IVF embryos were used for GO analysis. The categories most enriched under Biological Process GO term were- cellular process, metabolic process, regulation of biological process, response to stimuli. Similarly, in molecular Function GO term categories, of which the most affected ones were- binding activity and catalytic activity. In case of Cellular Component GO term, the cell part, cell and organelle being the most enriched ones (Supplementary Figure 8). In case of protein class GO term, the metabolite interconversion enzyme, protein modifying enzyme, nucleic acid binding protein and gene-specific transcriptional regulator being the most affected ones.
Pathway identification and network analysis
The KEGG pathway enrichment analysis helped us to demonstrate the relationship between the identified differentially expressing miRNA and their role in developmental processes of the embryos. At 2-cell stage, 126 pathways were detected, the major pathways related to embryonic development which were found to be affected in cloned embryos relative to IVF counterparts. Likewise, 122 pathways were detected at 8-cell stage cloned embryos relative to IVF counterparts. Also, 126 pathways were detected at blastocyst-stage cloned embryos relative to IVF counterparts (Supplementary figure 9). Through KEGG analysis we identified the most significantly enriched pathway during the course of embryo development which included apoptosis, cell cycle, MAPK signaling, mTOR, notch signaling, p53 pathway, PI3K-AKT signaling pathway, RAS signaling pathway, signaling pathways regulating pluripotency of stem cells, TGF beta signaling pathway, ERBB signaling pathway, Toll-like signaling pathway and WNT signaling pathway. Notch signaling pathway was specifically enriched pathway in 2-cell stage cloned embryo relative to IVF embryos.
KEGG analyses also revealed the significantly enriched metabolic pathway during the embryo development were pyrimidine metabolism (bta00240), purine metabolism (bta00230), fatty acid biosynthesis (bta00061), glycolysis/gluconeogenesis (bta00010).
Wnt signaling pathway is one of the important signaling pathway during embryo development. It plays a crucial role in regulating the development process of blastocyst stage. In our study, Wnt pathway was the most affected pathway during the course of embryo development. With the help KEGG Analysis we were also able to identify the target genes of upregulated miRNA which were predominantly enriched in WNT pathway (Figure 8). The total number of genes involved in Wnt signaling pathway and the number of genes affected in each stage is given in Table 2. mir-34, mir-345 and mir-409 were found to be upregulated in cloned embryos of 2-cell stage was predicted to regulate the member of Wnt signaling pathway such as WNT2, WNT9A, WNT3A, WNT7A,WNT7B. ALSO mir- 331, mir-339, mir-342 and mir-345 were predicted to target FZD2, FZD3, FZD5, and FZD7 at 2-cell stage. In 8-cell stage, mir-486, mir-487 and mir-493 were also detected to target WNT7B, WNT9B genes whereas in case of mir-338, mir-365, mir-371 and mir-378 were also found to target WNT3A, WNT9A, WNT7B and WNT2. Some more miRNAs such as mir-381, mir-455, mir-423, mir-93 and mir-574 were targeting DKK3, DKK4 and DKK1 gene respectively.
Validation of RNA seq data
In order to validate the RNA-seq data, seven miRNAs, which had been found to be expressed differentially between cloned and IVF embryos at different developmental stages (2-cell, 8-cell and blastocyst) were selected randomly and on the basis of their fold change for validation by qPCR analysis. These miRNAs were mir-218, mir-340 and mir-202 from 2-cell stage, mir-214, mir-96 and mir-139 from 8-cell stage and mir-370 from blastocyst stage. The two housekeeping miRNAs i.e., UniSp6, and miR-423 were used for the normalization of target miRNAs. The pattern and magnitude of relative expression level were found to be similar in RNA-seq and qPCR data (Figure 9).