Global DNA methylation profiles of buffalo (Bubalus bubalis) preimplantation embryos produced by handmade cloning and in vitro fertilization

Somatic cell nuclear transfer technique (SCNT) has proved to be an outstanding method of multiplication of elite animals but accompanied with low efficiency and live birth rate of cloned animals. Epigenetic alterations of DNA has been one of the culprits behind this issue. Cloned embryos are found to deviate slightly from regular pattern of demethylation and re-methylation at the time of nuclear reprogramming and embryonic development when compared with embryos produced by in vitro fertilization (IVF). Thus, the present study was aimed at evaluating global DNA methylation profiles of cloned embryos at 2-cell, 8-cell and blastocyst stages and compare it with corresponding stages of embryos produced by IVF by using MeDIP-Sequencing on Illumina-based platform. We found out that cloned embryos exhibited significantly different DNA methylation pattern as compared to IVF embryos with respect to distribution of differentially methylated regions in different components of genome, CpG islands distribution and methylation status, gene ontological profiles and pathways affected throughout the developmental stages. The data generated from MeDIP-Seq was validated at blastocyst stage cloned and IVF embryos by bisulfite-sequencing PCR on five randomly selected gene regions.

www.nature.com/scientificreports/ This involves changes in the expression profile of 10,000 to 12,000 genes which drive the embryonic and fetal development 24 . During nuclear reprogramming, the epigenetic marks specific to the differentiated somatic cell must be fully erased, and those associated with the totipotent stage of the embryo must be re-established. The epigenetic reprogramming after SCNT has been often found to be abnormal or aberrant 23 . Therefore, epigenetics plays a crucial role in determining the cloning efficiency and the fate of cloned offspring born. DNA methylation is one of the major epigenetic marks that plays crucial role at the time of nuclear reprogramming of SCNT embryos. DNA methylation occurs predominantly at CpG dinucleotides and is involved in several key genome functions such as imprinting, X-chromosome inactivation, genome stability, silencing of retrotransposons and inactivation of cancer-related genes 25,26 . Basically, DNA methylation suppresses gene expression by recruiting methyl-CpG-binding proteins such as MECP2 (methyl-CpG-binding protein 2), MBD1 (methyl-CpG-binding domain protein 1), MBD2, and MBD3, as well as associated histone deacetylases, corepressor proteins and chromatin remodelling machineries at the promoter regions of specific genes 27 . DNA methylation is necessary for transcription of octamer-binding transcription factor 4 (OCT4) which is critical factor for pluripotency thus, certifies initiation and maintenance of early embryonic development after successful cloning 28,29 . Therefore, aberrant DNA methylation may cause abnormal gene expression during initial stages of development, ultimately leading to embryo/fetal developmental failure 23,[30][31][32] .
Several studies have been made to investigate DNA methylation profile in SCNT embryos and compare it with their IVF counterparts. Umpteen of studies has been done with immunofluorescence technique to show aberrant pattern of DNA methylation in cloned embryos in many species. Dean et al. 33 used indirect immunofluorescence method and observed that cloned mouse embryos showed reduction in methylation only up to one-cell stage. Also, pattern of DNA methylation of many cloned embryos was found to be similar to donor fibroblast cells. Santos et al. 34 reported genome-wide failure of nuclear reprogramming in cloned bovine pre-implantation embryos at blastocyst stage using immunofluorescence. H3-K9 and DNA hypermethylation was found to be the major reason behind this failure. Likewise, aberrant DNA methylation pattern was observed along with variable methylation levels among nuclei within swamp buffalo SCNT embryos 35 . Histone hyperacetylation was also observed at 4-and 8-cell stage of SCNT embryos compared with IVF counterparts. Bisulfite sequencing was further comprehensively used to establish genome-wide methylation in cloned embryos. DNA methylation of bovine alpha satellite I repeat sequence was found to hypermethylated in blastocyst as well as in SCNT foetuses in comparison with in-vivo produced embryos 36 . Couldrey et al. 37 reported repetitive satellite I sequence was found to differ significantly with regards DNA methylation and cell type (inner cell mass or trophectoderm) between cloned and in-vivo produced embryos. Likewise, comprehensive genome-wide methylome was also studied by applying several other Next Generation Sequencing (NGS) strategies such as Whole Genome Bisulfite Sequencing (WGBS) [38][39][40] and Reduced Representation Bisulfite Sequencing (RRBS) 41,42 . To our information, there is only one report available on genome-wide DNA methylation profile of SCNT embryos in farm animals 43 . However, there is no report till date on the global DNA methylation profiling examined using a high throughput NGS technique in embryos produced by SCNT and IVF in any species. Thus, in the present study, we carried out genome-wide DNA methylation profiling in buffalo embryos produced by handmade cloning (HMC) and in-vitro fertilization (IVF) by using methylated DNA immunoprecipitation combined with high-throughput sequencing (MeDIP-Seq). To our knowledge, this is the first report on global DNA methylation pattern of cloned embryos relative to IVF counterparts in buffalo.

Results
Under this study, global DNA methylation profile of pre-implantation buffalo embryos at 2-cell, 8-cell and blastocyst stages, produced by HMC and IVF, was generated by MeDIP-Seq for comparing the differentially methylated regions (DMRs) between cloned and IVF embryos at each developmental stage.
Production of cloned and IVF embryos. In this experiment, somatic cells from an adult bull  were used for the production of cloned embryos at 2-cell, 8-cell and blastocyst stages by HMC. Semen obtained from the same bull was used for the production of IVF embryos at each of these stages (Fig. 1A). This was done to minimize the genetic variability between the cloned and IVF embryos since the embryos of the two groups produced in this manner were genetically half-identical. Blastocyst rate for cloned embryos varied from 30 to 35% for 15 trials whereas cleavage rate for 2-cell and 8-cell developmental stage was recorded to be 95-98% and 70-80% respectively. Likewise, for IVF embryos, blastocyst rate was observed to vary from 10-15% rate for 10 trials with cleavage rate of 15-22% was observed for 2-and 8-cell stages.
It has been reported that poor quality and low developmental competence of oocytes lead to lower cloning and IVF efficiency in terms of blastocyst rate and quality of embryos 44 . Therefore, it is necessary that competent and good quality oocytes are used for the production of cloned and IVF embryos that can be subjected to MeDIP-Seq. We achieved this by two means. Firstly, we selected only those embryos which appeared to be morphologically normal in appearance. Second, we selected oocytes of high developmental competence for cloning and IVF by staining them with Brilliant Cresyl Blue (BCB). It has been shown in a previous study in our laboratory that the blastocyst rate was higher for BCB + (oocytes with high developmental competence) than for BCB-(oocytes with low developmental competence) oocytes. Moreover, the blastocysts produced from BCB + oocytes had inner cell mass (ICM) cell number, ICM/trophectoderm (TE) cell number ratio, global level of H3K18ac, apoptotic index and expression level of BCL-XL similar to that of blastocysts produced through IVF 45 . Therefore, BCB + blastocysts had better developmental competence and were closer to IVF blastocysts in terms of quality, epigenetic status and gene expression than BCB-blastocysts, In view of this, only BCB + oocytes were used for HMC and IVF in the present study.   Supplementary Fig. 1). The alignment statistics of the three biological replicates each of cloned and IVF embryos are given in (Supplementary Table 1).
Chromosomal distribution of MeDIP-Seq reads. The chromosomal distribution of the three biological replicates each of cloned and IVF 2-cell, 8-cell and blastocyst stage embryos is given in (Fig. 1B). MeDIP-Seq raw reads were detected in most of chromosomal regions (chromosome 1-29 and chromosome X) at the 2-cell, 8-cell and blastocyst stage in both cloned and IVF embryos. Maximum number of reads was found to map on chromosome number 1 as chromosome 1 is the largest chromosome in Bos taurus genome.
MA and volcano plot of cloned vs IVF embryos at different developmental stages. MA plot is commonly used to visualize the log fold change (FC) versus mean expression between two samples. It is a type of 2D scatter plot which allows visual display of base-2 log FC along the Y-axis and normalized mean expression along the X-axis, where each data point represents single CpG site. In the present study, this plot has been used to visualize log FC versus mean expression between cloned relative to IVF embryos at different developmental stages. Each data point with extreme values on the Y-axis represents the CpGs with significant differential methylation. Data points falling above the 1 threshold on the Y-axis indicate significant number of CpGs which were hypermethylated whereas, those below − 1, indicate CpGs which were hypomethylated. The data points falling close to 0 on the Y-axis indicates that both the samples share highly similar methylation pattern. Volcano plot is another type of scatter plot which allows quick visual identification of those data points (CpGs, genes etc.) that display changes of high magnitude which are statistically significant on plotting the measure of statistical significance (p-value) against the magnitude of the change (FC). In the present study, this scatter plot has been used to predict the status of CpGs which are significantly differentially methylated, either hyper-or hypomethylated in cloned relative to IVF embryos at given stages of cloned and IVF embryos (Fig. 1C). The CpGs sites seen dispersed towards the right side represent the hypermethylated sites, and those seen dispersed towards the left side represent the hypomethylated sites. The CpGs seen towards the top are those which are statistically most significantly differentially methylated.
Overall distribution of differentially methylated regions (DMRs). The percentage of distribution of DMRs throughout the three developmental stages examined i.e., 2-cell, 8-cell, and blastocyst stages in cloned and IVF embryos is given in (Fig. 1D). In the cloned embryos, the number of highly methylated regions (HMRs) which was found to be higher at the 2-cell stage, decreased sharply at the 8-cell stage and then increased to the highest level at the blastocyst stage. In IVF embryos, number of HMRs increased from the 2-cell stage to the 8-cell stage and then declined sharply at the blastocyst stage. Also, distribution of HMRs in cloned 2-cell and 8-cell stage embryos was significantly (padj ≤ 0.05) less in comparison with IVF counterparts. The blastocyst stage of cloned embryos showed significantly (padj ≤ 0.05) higher numbers of HMRs. These results suggest that DNA methylation is altered dynamically throughout the developmental stages in cloned embryos in comparison with their IVF counterparts.
Determination of DMRs at each developmental stages. Total number of DMRs were detected at each developmental stage of cloned and IVF embryos. The number of uniquely present DMRs in cloned and IVF embryos at 2-cell stage was found to be 90,501 (11.2%) and 1,49,925 (18.6%), respectively. From among the commonly present DMRs i.e., 5,66,734, a total of 2,81,276 were found to be hypermethylated in cloned relative to IVF embryos. The number of uniquely present DMRs in cloned and IVF embryos at 8-cell stage was found to be 68,826 (14.7%) and 1,78,171 (38%), respectively. From among the commonly present DMRs i.e., 2,21,383, a total of 1,04,174 were found to be hypermethylated in cloned relative to IVF embryos. The number of uniquely present DMRs in cloned and IVF embryos was 1,49,797 (18.4%) and 1,02,839 (12.6%) respectively, whereas, 5,62,134 were commonly expressed in the two groups. From among the commonly present DMRs 2,82,582 were hypermethylated in cloned relative to IVF embryos. The Venn diagrams showing the DMRs commonly present in cloned and IVF 2-cell, 8-cell and blastocyst stage embryos and the DMRs uniquely present in the two groups are presented in (Fig. 2A).

Distribution of DMRs at gene level.
Genes which overlap with HMRs in the promoter or gene-body regions are considered as methylated genes. If for HMRs coverage, edge R. logFC = + 1 then the gene is designated as hypermethylated and if edge R. logFC = − 1, it is designated as unmethylated. In the present study, overall 27,518 genes were identified at the 2-cell stage, out of which 12,392 genes were found to be hypermethylated in cloned relative to IVF embryos. In case of 8-cell stage embryos, the total number of genes was found to be 23,621 of which, 14,594 genes were hypermethylated. Similarly, at the blastocyst stage, 24,556 genes were identified out of which, 9111 genes were hypermethylated in cloned relative to IVF embryos (Fig. 2B), Supplementary Table 2). The maximum number of hypermethylated genes was found at the 8-cell stage whereas, the minimum number was observed at the blastocyst stage. The number of genes showing hypermethylation at 2-cell, 8-cell and blastocyst stage in cloned and IVF embryos is shown in (Supplementary Table 3  www.nature.com/scientificreports/ Distribution of DMRs in different components of the genome. Genome-wide DNA methylation profile in cloned and IVF embryos was observed with the distribution of MeDIP-Seq reads in different component of genome including upstream and downstream 1 kb and 2 kb promotor, exons, introns, intergenic region, CpG islands and repetitive elements. Out of the 8,07,160 DMRs evaluated (P < 0.05), 3,71,777 DMRs were found to be hypermethylated in cloned relative to IVF 2-cell stage embryos. Of 4,68,379 DMRs evaluated (P < 0.05), 1,73,000 DMRs were found to be hypermethylated in cloned relative to IVF 8-cell stage embryos. Out of the 8,14,770 DMRs evaluated (P < 0.05), 4,32,379 DMRs were found to be hypermethylated in cloned relative to IVF blastocyst stage embryos. Major part of hypermethylated DMRs found at 2-, 8-, and blastocyst stages belong to intergenic region at upstream and downstream 1 kb and 2 kb, respectively. Venn diagram showing the distribution of hypermethylated DMRs in different genomic regions (promoter, exon, intron and intergenic) in cloned relative to IVF 2-, 8-cell and blastocyst stage embryos is presented in (Fig. 2C). The number and percentage of hypermethylated DMRs in upstream and downstream 1 kb and 2 kb, upstream 2 kb and downstream 2 kb in cloned relative to IVF embryos at each developmental stage is presented in (Supplementary Table 4) respectively.
CpG islands in cloned and IVF embryos. Analysis of the percentage of reads which mapped distinctly with CpGs islands (CGIs) in all the biological replicates of cloned and IVF embryos at 2-cell, 8-cell and blastocyst stage showed a range of 10.95% to 22.91% (Fig. 3A). Maximum percentage of reads which mapped with CGIs for three biological replicates of cloned and IVF embryos at 2-cell stage was 11% and 16.62%, respectively. Similarly, for 8-cell stage cloned and IVF embryos, it was found to be 22.91% and 16.37%, whereas, in case of blastocysts, the corresponding values were 16.39% and 17.87%, respectively. In the present study, the total number of CpG islands distributed in buffalo embryo genome were observed. For the first time, we identified total 2,17,246 CpG islands, out of which 1,04,188 (42.55%) and 1,13,058 (48.07%) islands were present in cloned and IVF embryos, respectively. In the cloned embryos, 52,168 (23.72%) were found to be hypermethylated, whereas, the corresponding figures in IVF embryos were 41,365 (30.25%). The number and percentage of CpG islands hypermethylated in cloned versus IVF embryos at each developmental stage (i.e., 2CC vs 2CI, 8CC vs 8CI, CBL vs BLI) is given in (Table 1).
CpGs enrichment in cloned and IVF embryos. The assessment of CpGs enrichment value of methylated regions (frequency of CpGs wise) was done by dividing the relative frequency of CpGs (region.relH) of the regions by the relative frequency of CpGs (genome.relH) of the reference genome. Similarly, the final value of CpGs enrichment of methylated regions (ratio CpG wise) was done by dividing the observed/expected ratio of CpGs within the region (region.goge) by the observed/expected ratio of CpGs within the reference genome (genome.goge). The enrichment score for both frequency CpGs and ratio CpGs in the genomic regions sequenced compared to the reference genome was above 1 in all the cases (Table 2), indicating a successful enrichment of methylated fragments in the data sets.
Distribution of CpG islands in repetitive elements of genome. Distribution of CpG islands with hypermethylation in the repetitive elements of genome i.e., LINE, SINE, repetitive DNA, LTR, simple repeats and low complexity repeats etc. was examined between cloned and IVF embryos at different developmental stages and between different developmental stages (Supplementary Table 5). Relative percentage of CpG islands found in particular repeats in comparison with reference genome (UMD 3.1.1) at different developmental stages in cloned and IVF embryos is presented in (Fig. 3B). It was observed that at 2-, 8-cell and blastocyst stage cloned embryos, maximum distribution of CpGis was found in LTRs, simple repeats and LINE respectively.
Gene ontology (GO) analysis. GO analysis was done in order to explore the biological significance, detailed annotation of gene function, biological process and cellular distribution of differentially methylated genes between cloned and IVF embryos. GO analysis is based on a controlled vocabulary of terms of three domains-"Biological Process" (BP), "Molecular Function" (MF) and "Cellular Compartment" (CC) 46 . Using GO terms, MeDIP-Seq results can be summarized to provide insights into the methylation status of genes between cloned and IVF embryos. In the present study, genes with hypermethylation (Edge R. Log fold change ≥ + 1) in cloned relative to IVF embryos were used for GO analysis.
A total of 483 Biological Process, 85 Molecular functions and 121 cellular components were found to be affected at 2-cell stage cloned embryos relative to IVF counterparts. GO analysis of differentially methylated genes in cloned and IVF 8-cell stage embryos revealed a total of 526 Biological Processes, 97 Molecular functions and 115 cellular components were found to be affected with DNA hypermethylation. Likewise, a total of 594 Biological Process, 131 Molecular functions and 176 Cellular components were affected in blastocyst stage cloned embryos relative to IVF ones. Among all three stages of cloned embryos relative to IVF counterparts, molecular function, binding, and protein binding were majorly affected in case of Molecular function GO term category. Cellular Component GO term including cell part, cellular component, and intracellular anatomical structure being the most enriched ones (Fig. 3C). The functional annotations, which were found to be hypermethylated are given in (supplementary Tables 6, 7, 8).
Pathways analysis. At 2-cell stage cloned embryos, total 30 pathways were detected, out of which 15 pathways were affected by hypermethylation of genes related to embryonic development relative to IVF counterparts. Out of 28 pathways detected, 15 pathways were affected by hypermethylation of genes related to embryonic development in 8-cell stage cloned relative to IVF embryos. At blastocyst stage cloned embryos, 37 pathways were found, of which 27 pathways were affected by hypermethylation of genes related to embryonic development  www.nature.com/scientificreports/ relative to IVF embryo (Table 3). Wnt signaling pathway, inflammation mediated by chemokine and integrin signaling pathway and gonadotropin-releasing hormone receptor pathway were majorly affected by hypermethylation of DNA throughout the developmental stages.
Methylation status of important genes involved in embryonic development in cloned relative to IVF embryos. Methylation status of some important genes involved in embryonic development viz. pluripotency-related genes, imprinted genes, apoptosis-related genes, cell cycle-related genes, methylationspecific genes and some functionally important genes is given in (Table 4). A particular gene is depicted as hypermethylated if it overlaps with highly methylated region (HMR) even in a single gene region. However, the actual genes expression level of a particular gene can vary irrespective of presence of methylation in any region of that gene. There is a complex relationship between DNA methylation and gene expression as higher gene expression levels are generally linked with low methylation at promotor region but higher gene body methylation 47 . The exact expression status of a gene can be known only after examining the correlation of DNA methylation with mRNA expression profile of that particular gene.

Results validation.
To confirm the results of MeDIP-Seq data, five genes (IGF2R, PEG10, MDM2, TCG7 and COL4A1) were selected randomly to perform bisulfite sequencing PCR (BSP). Primers were designed using MethPrimer software. The total of number of blastocysts used for bisulfite sequencing PCR in cloned and IVF groups was 25 for each group. Each sample was sequenced in triplicate. Data generated from Sanger sequencing was analysed by BioEdit version 7.2.5. All the five gene regions which had showed hypermethylation in cloned relative to IVF embryos in MeDIP-Seq data, showed hypermethylation in bisulfite sequencing PCR also which validated and confirmed the reliability of our data obtained from MeDIP-Seq (Fig. 4).

Discussion
Although, some reports are available in bovine and other species on the DNA methylation in cloned compared to IVF embryos at different developmental stages, these reports provide a limited insight on the aberrations in DNA methylation in cloned embryos since low throughput techniques like immunostaining and RT-PCR were used for examining the methylation status in these studies [33][34][35] . To our information, there is no report in any species on the comparative global DNA methylation profile of cloned and IVF embryos at different developmental stages using a high throughput technique like MeDIP-Seq. Over the years, many methods have been developed for analyzing genome-wide DNA methylome. These include whole genome bisulfite sequencing (WGBS), reduced representation bisulfite sequencing (RRBS), TET-assisted bisulfite sequencing (TAB-seq), comprehensive high-throughput arrays for relative methylation (CHARM), methylation array and MeDIP-Seq. However, each method has its own merits and demerits when it comes to sequence whole genome. WGBS is considered to be the best method for examination of DNA www.nature.com/scientificreports/ methylation profile of mammalian genome due to its high resolution. But it is costly, time consuming and generates complicated data that further needs comprehensive bioinformatics analysis which limits the popularity of this technique. MeDIP-Seq is a widely used affinity-based method for examining DNA methylation. Global DNA methylation analysis is feasible with MeDIP-Seq even with very low amount of starting (as low as 1 ng) DNA sample 48,49 . Many recent studies have used MeDIP-Seq to generate relative genome-wide DNA methylation profile [50][51][52] . Therefore, we selected this method to study genome-wide DNA methylation in cloned embryos in comparison with IVF embryos in buffalo.
In the present study, DNA methylation profile of cloned and IVF embryos generated by MeDIP-Sequencing showed significant difference between both kind of embryos in reference to number and distribution of DMRs and CpGis in different parts of genome. Also, cloned embryos were found to be hypermethylated with respect to IVF counterparts which affected various important pathways involved in embryonic development. Hypermethylation of DNA throughout all the developmental stages studied in cloned embryos affected the methylation status of many functional, structural and imprinted genes that may lead to abnormal gene expression.
An unusual pattern of de-methylation and re-methylation in cloned embryos has been reported in several 5-methylcytosine immunofluorescence-based studies in past years [33][34][35]43 . Our study was not immunofluorescence-based but our MeDIP-Seq data is in agreement with the results of these low-throughput studies as cloned embryos showed hypermethylation of DNA at blastocyst stage in comparison to IVF counterparts.
Dean et al. 33 reported that cloned bovine embryos lack passive demethylation mechanism and have high expression of DNMT3A and DNMT3B which is responsible for de novo methylation at 4-cell and 8-cell stages. Likewise, in the present study, 23,621 genes were identified at the 8-cell stage in cloned embryos out of which 14,594 showed hypermethylation. This scenario reflects re-methylation of genes due to re-establishment of methylation marks that coincides with time of embryonic genome activation (EGA). Memili et al. 53 suggested the same possible reason for this pattern of methylation at 8-cell stage as this stage involves major change in transcriptional activation of the embryonic genome which might activate Dnmt3a and Dnmt3b enzymes.
Dean et al. 33 reported increased methylation level in the bovine blastocyst stage embryos specifically in the trophectoderm cells. Contrary to this, we found that the less number of genes were showing hypermethylation (9111) in cloned embryos relative to IVF embryos at the blastocyst stage.
Using single-cell RRBS to generate DNA methylomes for SCNT embryos in mice, Liu et al. 54 reported abnormally high levels of DNA methylation in 2-cell and 4-cell stage SCNT embryos. It was found that, SCNT embryos at 4-cell stage were relatively hypermethylated compared with 2-cell stage embryos which opposed the demethylation pattern of normal embryos. Also, the aberrant DNA methylation in arrested SCNT embryos was attributed to differences in expression of Dnmt1 and Tet1 enzymes. However, our study showed a slightly different trend as we found that 27,517 genes were identified in 2-cell stage cloned embryo, with 12,392 genes showing hypermethylation at 2-cell stage cloned embryos in comparison to IVF counterparts. This suggests that the majority of genes undergo demethylation at early developmental stages during nuclear reprogramming. In our data, DNMT1 and TET1 showed hypermethylation which could be the reason for lower methylation level of genes in 2-cell stage cloned embryos relative to their IVF counterparts. It is suggested that cytoplasmic factors of oocytes might help in triggering demethylation with remodelling of chromatin.
To our information, there is only one study in which genome-wide DNA methylation has been studied in SCNT preimplantation embryos but it is focused on specific parts of genome. Zhang et al. 43 compared DNA methylation reprogramming in bovine SCNT and IVF preimplantation embryos and analyzed the influence of vitamin C (VC) on the reprogramming of DNA methylation by using bisulfite sequencing. The results showed that global DNA methylation followed a typical pattern of demethylation and remethylation in IVF preimplantation embryos however, the global genome remained hypermethylated in SCNT preimplantation embryos. Compared with the IVF group, pluripotency genes POU5F1 and NANOG showed insufficient demethylation and hypermethylation in the SCNT group. Similarly, in our study, we found that pluripotency genes (POU5F1, TCF3, LIN28b, and DUSP1) showed high number of regions with hypermethylation which suggests lower expression of these genes.
Gene ontology analysis (GO) of differentially methylated genes in the present study showed 483 biological process, 85 molecular function and 121 cellular components to be affected by hypermethylation of genes in cloned 2-cell stage embryos relative to their IVF counterparts. Similarly, in the 8-cell stage cloned embryos, 526 biological process, 97 molecular function and 115 cellular components were found to be affected by hypermethylation of genes with respect to IVF embryos. GO analysis for cloned blastocyst-stage embryos relative to IVF embryos revealed 594 biological process, 131 molecular function and 176 cellular component to be affected by hypermethylation.
In the present study, 15 pathways were found to be affected at the 2-cell stage, 15 at the 8-cell stage and 27 at the blastocyst stage in cloned embryos relative to IVF counterparts. For all the developmental stages covered  www.nature.com/scientificreports/ under present study, the major pathways affecting embryonic development in cloned embryos relative to IVF counterparts were found to be Wnt signaling pathway, inflammation mediated by chemokine and cytokine signalling pathways, gonadotropin-releasing hormone receptor pathway and Integrin signalling pathway. These pathways were affected by hypermethylation of DNA. It was found in the present study that imprinting genes such as IGF2R, PEG10, GRB10, MEST, and MAOA showed presence of hypermethylated regions that suggest lower gene expression at the blastocyst stage in cloned relative to IVF embryos. Some important apoptosis-related genes such as CASP9 and DFFA, FAS, APAF1, and PTEN were found to have hypermethylation in gene region in blastocyst stage cloned embryos compared to their IVF counterparts.
Similarly, cell cycle-related genes such as MDM2, and NPAT showed presence of hypermethylation in gene regions in blastocyst stage cloned embryos relative to IVF counterparts.
We found that the methylation status of methylation-specific genes viz. DNMT1, DNMT3A, DNMT3B, TET1, EZH2, and MBD3 showed no methylation in gene regions with whereas TET2 and TET3 gene found to have hypermethylation in gene regions in blastocyst stage cloned embryos relative to their IVF counterparts. This pattern of methylation reflects high expression of abovementioned genes in the cloned blastocyst-stage embryos which might be responsible for the abnormal hypermethylation observed in cloned blastocysts in comparison with IVF counterparts.
There are some functionally important genes which contribute to embryonic development in bovine embryos. We found that YY1, SOD1, IL6, G6PD, PCGF2, MCL1, and IFN-τ represented no methylation in the gene regions whereas TIMP2, MTPN, COL4A1, TCF7, AQP9, BMP7, FGF7, WNT2, DKK3, CRABP2, PPARA, DYSF, LEF1 and HDAC8 showed high number of regions with hypermethylation in blastocyst stage cloned embryos compared to their IVF counterparts. Shyam et al. 55 also observed a high expression of WNT signalling pathway genes i.e. TCF7 and LEF1 in cloned buffalo blastocyst-stage embryos. They also reported lower gene expression level of IFN-τ (primary maternal signal for recognition of pregnancy in ruminants) in cloned blastocyst was lower than in IVF counterparts.
Similarly, Sood et al. 56 also reported down-regulated expression of IFN-τ in cloned embryos in comparison with IVF counterparts. In contrast, present study indicates hypermethylation of TCF7 and LEF1 at gene level which suggests lower expression of these genes in cloned embryos. Also, we found the gene region for IFN-τ to be unmethylated which points to a high expression of this gene in cloned compared to IVF embryos. Thus, our results regarding methylation status of TCF7, LEF1 and IFN-τ genes were found to be not in agreement with above mentioned studies.
Kim et al. 57 reported abnormally up-regulated expression of tissue inhibitor of metalloproteinase-2 (TIMP-2) and superoxide dismutase (SOD) in the placentae of SCNT cloned Korean native cattle that died immediately after birth and in normal placentae obtained by artificial insemination. Abnormal expression of these two genes reported to be partly responsible to for abnormal placental function and low survivability of cloned claves. Results from our study regarding these two gene showed that showed hypermethylation in gene region of TIMP-2 gene and no methylation in gene regions of SOD1 gene in cloned blastocysts relative to IVF counterparts. For TIMP2 gene, our results are not in agreement with those of above study but our results for SOD1 gene, are in agreement with above study as its unmethylated status can suggesting high level of gene expression.
Gene such as AQP9 encode a water channels protein which is involved in coordinating urea transport. Gao et al. 58 reported abnormally lower expression level of AQP9 gene in placenta of SCNT cattle by both RNA-seq and q-PCR. Similarly, in present study, methylation level of AQP9 gene was found to be very high which suggests lower expression level of this gene in cloned embryos relative to IVF counterparts. Sood et al. 56 also reported down-regulated expression of DKK3 gene in cloned embryos in comparison with IVF embryos. Likewise, our study showed high number of hypermethylated regions of DKK3 gene in cloned blastocysts relative to IVF   www.nature.com/scientificreports/ counterparts that might indicate towards low expression level of this gene. Together, this study suggests that cloned embryos and IVF embryos present a profound differences in global DNA methylation profile which is consistent throughout all developmental stages. These differences are related to distinct DMRs and CpGis distribution pattern, gene ontology terms and pathways in cloned embryos relative to IVF counterparts. A large number of pathways are found to be affected by hypermethylation in cloned relative to IVF embryos. Among these, the major pathways related to embryonic development are Integrin signalling pathway, Wnt signaling pathway, apoptosis signaling pathway, inflammation mediated by chemokine & cytokine signaling pathway, gonadotropin-releasing hormone receptor pathway, CCKR signaling map, angiogenesis. Further investigation is needed to understand the crosstalk among other epigenetic mechanisms that all together regulate these pathways.
There is a need of further more targeted studies stimulated by this study to evaluate hypermethylation in exact genomic locations, promoters and gene body and role of regulatory proteins in establishing and eliminating methylation in the genome. It is essential to clarify whether this DNA hypermethylation alone is responsible for causing low efficiency of blastocyst rate and live birth rate of cloned offsprings and how approved epigenetic modifiers drugs can be used to prevent developmental failure of cloned embryos.

Methods
Ethics statement. All   In vitro embryo production by HMC and IVF. Buffalo ovaries from the slaughterhouse were aspirated for Cumulus-oocyte complexes (COCs). COCs were classified into useable and unusable categories based on their morphological appearance as ones with evenly granular homogenous ooplasm and unexpanded cumulus mass with ≥ 2 layers of cumulus cells are considered as useable. The oocytes which were wholly or partially denuded of the cumulus mass or which had expanded cumulus mass or irregular ooplasm were classified as non-useable. For selection of COCs of high developmental competence, those of useable category were stained with Brilliant Cresyl Blue (BCB), as described previously 45 . COCs with distinct blue colour (BCBþ) were washed several times with the IVM medium, which consisted of TCM-199 + 10% FBS þ 5 mg/mL porcine follicle stimulating hormone (pFSH) + l mg/mL estradiol-17β + 0.81 mM sodium pyruvate + 50 mg/mL gentamicin sulfate. For IVM, the COCs were placed in 100 mL droplets of the IVM medium (15-20 COCs per droplet), overlaid with sterile mineral oil in 35 mm Petri dishes and cultured in a CO 2 incubator at 38.5 °C for 21 h as previously described 60,61 . Culture of donor fibroblast cells to full confluence to synchronize them in G1 stage of cell cycle and all other procedures of HMC were carried out as described previously 62 . For production of IVF embryos, oocytes were subjected to IVM for 24 h to achieve metaphase II stage with production of first polar body. IVF A portion of the sonicated DNA was reserved as Input (control) and was not subjected to immunoprecipitation. The denatured DNA was placed on ice for 10 min and was then incubated with 5mC monoclonal antibody for 2 h at 4 °C on an orbital shaker (50-100 rpm). After immunoprecipitation, methylated DNA fragments were column purified and eluted in 20 µL of elution buffer and stored at − 20 °C till further used for library preparation.

MeDIP-sequencing.
Library preparation and quality control. Global DNA methylation analysis of Handmade cloned and IVF embryos at 2-cell, 8-cell and blastocyst stages was done using MeDIP-Seq by a commercial service provider, M/S DNA Xpert, New Delhi, India. Illumina library preparation was performed by using the NEBNext Ultra DNA Library Prep Kit for Illumina (NEB, New England Biolabs Inc., USA #E7370S). Sonicated DNA (25 ng from each pool) was end repaired into dA-tailed fragments. Then the NEBNext-indexed adaptor was ligated into sonicated DNA fragments. Adaptorligated DNA was further cleaned up by AMPure XP beads (Beckman Coulter, Inc., Life Sciences, United States #A63881) and then subjected to the MeDIP procedure. The Input DNA was linked to different index primers. After immunoprecipitation, adaptor-ligated DNA was extracted by phenol/chloroform and precipitated by ethanol. The antibody-enriched DNA was subsequently amplified by using an index primer and universal PCR primer provided in NEB Multiplex Oligos for Illumina. PCR was performed with 20 cycles for each MeDIP sample. The PCR products were then purified by AMPure XP beads and eluted with TE buffer. The library concentration and size distribution was analyzed by Agilent 4200 Tapestation.The library was run on Agilent 4200 Tape Station and the resulting profile was evaluated. QC-passed libraries (peak size varies from 28 to 400 bp) were used to generate high-quality reads for further analysis.
Sequencing and quality control. The QC-passed were then set up on Illumina HiSeq 2500 platform for sequencing. Quality check of raw reads was carried out using FastQC (v0.11.5). Based upon Phred score value, all the reads were found to be of good quality, therefore, no filtering was required for any of the samples. The QC-passed high-quality reads (Phred Score Cut off of Q 20) were used for further analysis.
MeDIP-Seq data analysis. R package MEDIPS. R (v3.6.3) and MEDIPS (v1.36.0) softwares was used for analysis of raw data. The MeDIP-Seq data analysis workflow employed in the present study consisted of following major steps viz., loading of raw data into MEDIPS software, quality check of data, identification of differentially methylated regions (DMRs) based upon fold change criteria (with length of DMR around 1999 bp in the program window) and classification/clustering of genes. The gene ontology (GO) enrichment analysis and pathway analysis was done by using R package MEDIPS and Panther Classification System, respectively. The reads generated were aligned using bowtie2 (v2.2.8) software to Bos taurus reference genome, UMD 3.1.1. CpG islands of genome were identified using cpgiscan (v1.0). CpGi annotation for each sample was done by intersecting bam file against the cpgiscan bed output using bedtools. Distribution of CpGs and feature annotation at different components of genome (promoter,exon, intron and intergenic region) was done by using R package genomation. genomation (v1.16.0). Distribution of CpG islands and their methylation status in different repetitive elements of genome (LINE, SINE, repetitive DNA, LTR, simple repeats, low complexity repeats) was done by using RepeatMasker.
Bisulfite sequencing PCR analysis. Five gene regions were selected randomly to validate MeDIP-Seq data. For that, 25 blastocyst-stage for each cloned and IVF embryos were used for validation of MeDIP-Seq data by bisulfite sequencing PCR analysis. Genomic DNA was isolated from respective pools of cloned and IVF blastocyst-stage embryos using the DNeasy ® Blood & Tissue Kit (Qiagen, catalogue no. 69504) according to the manufacturer's instructions. The purified DNA thus obtained is ready to be used for further processing or can be stored at − 20 °C. DNA quality control analysis was carried on NanoDrop™ 2000 (Thermo Scientific, USA). Bisulfite conversion of DNA samples isolated from cloned and IVF blastocyst-stage embryos was done using EZ DNA Methylation-Direct™ Kit (Cat. No. D-5020, Zymo Research, Orange, CA, USA) according to the instruc- www.nature.com/scientificreports/ tion manual, with minor modifications as detailed; Genomic DNA isolated using above protocol was used for bisulfite conversion. DNA (20 μL) was added to 130 μL of CT Conversion Reagent solution and was centrifuged briefly to ensure that there were droplets in the cap or on sides of the tube. The mixture was then heated in a thermal cycler at 98 °C for 8 min and 64 °C for 3.5 h and was then stored at 4 °C for up to 20 h. Bisulfite converted DNA was desalted, purified and eluted with 40 μL of elution buffer. The O.D. of the bisulfite-converted DNA was recorded. Consequently, bisulfite sequencing PCR (BSP) was carried out with 1 μL of modified DNA per PCR reaction Meth-primers for selected five genes were designed from The Li Lab (MethPrimer) software at 2 kb upstream regions covering at least one CpG island. For designing bisulfite PCR primers, the sequences for each gene were taken from NCBI (http:// www. ncbi. nlm. nih. gov). The sequences of the BSP primers used to amplify the targeted products are shown in (Supplementary Table 9). Hotstart PCR was carried out for bisulfiteconverted DNA pairs using Hotstart DNA polymerase (Zymo Taq™ premix, Zymo Research). The amplified product was taken as template DNA (with 1:5 dilution) for second PCR amplification with the same primers. The PCR cycle included denaturation for 10 min at 95 °C followed by 35 repeated cycles of 95 °C for 30 s, annealing at variable temperatures for 30 s and extension at 72 °C for 30 s followed by a final extension at 72 °C for 10 min. Resulting PCR products were run on 2% agarose gel to confirm the amplicon size. Purification of the bisulfite-treated amplified product was done with GeneJET PCR Purification Kit (Thermo Fisher Scientific, Inc. USA #K0701). Purified PCR products were sent for custom sequencing to M/S DNA Xpert Pvt. Ltd., (New Delhi, India). The PCR products were sequenced in triplicates with both forward and reverse primers of selected region of interest. The data files generated with Sanger sequencing were processed using the Sequencing Analysis V5.3 software. Sequences with non-overlapping peaks are considered to be good for further analysis. Peaks with high background of < 50% of the co-efficient of variation (CV) are also considered to be good for further analysis. By using Seq Scanner, PDF and Fasta files were created from ABI files. Analysis of data generated for calculating methylation level (%) and creating lollipop diagrams for each sample was done by using BioEdit version 7.2.5 software.

Data availability
The high-throughput sequencing data (Raw data/ FastQ) have been submitted in NCBI's Sequence Read Archive (SRA) under BioProject ID: PRJNA723972.