Genome-Wide DNA Methylation Patterns of Bovine Blastocysts Developed In Vivo from Embryos Completed Different Stages of Development In Vitro

Early embryonic loss and altered gene expression in in vitro produced blastocysts are believed to be partly caused by aberrant DNA methylation. However, specific embryonic stage which is sensitive to in vitro culture conditions to alter the DNA methylation profile of the resulting blastocysts remained unclear. Therefore, the aim of this study was to investigate the stage specific effect of in vitro culture environment on the DNA methylation response of the resulting blastocysts. For this, embryos cultured in vitro until zygote (ZY), 4-cell (4C) or 16-cell (16C) were transferred to recipients and the blastocysts were recovery at day 7 of the estrous cycle. Another embryo group was cultured in vitro until blastocyst stage (IVP). Genome-wide DNA methylation profiles of ZY, 4C, 16C and IVP blastocyst groups were then determined with reference to blastocysts developed completely under in vivo condition (VO) using EmbryoGENE DNA Methylation Array. To assess the contribution of methylation changes on gene expression patterns, the DNA methylation data was superimposed to the transcriptome profile data. The degree of DNA methylation dysregulation in the promoter and/or gene body regions of the resulting blastocysts was correlated with successive stages of development the embryos advanced under in vitro culture before transfer to the in vivo condition. Genomic enrichment analysis revealed that in 4C and 16C blastocyst groups, hypermethylated loci were outpacing the hypomethylated ones in intronic, exonic, promoter and proximal promoter regions, whereas the reverse was observed in ZY blastocyst group. However, in the IVP group, as much hypermethylated as hypomethylated probes were detected in gene body and promoter regions. In addition, gene ontology analysis indicated that differentially methylated regions were found to affected several biological functions including ATP binding in the ZY group, programmed cell death in the 4C, glycolysis in 16C and genetic imprinting and chromosome segregation in IVP blastocyst groups. Furthermore, 1.6, 3.4, 3.9 and 9.4% of the differentially methylated regions that were overlapped to the transcriptome profile data were negatively correlated with the gene expression patterns in ZY, 4C, 16C and IVP blastocyst groups, respectively. Therefore, this finding indicated that suboptimal culture condition during preimplantation embryo development induced changes in the DNA methylation landscape of the resulting blastocysts in a stage dependent manner and the altered DNA methylation pattern was only partly explained the observed aberrant gene expression patterns of the blastocysts.


Introduction
In vitro embryo production (IVP) using oocytes matured and fertilized under in vitro culture condition has been a common practice for commercial and research purposes. Despite years of optimization, early embryonic losses, placental dysfunction, fetal death and over sized fetuses are still observed in embryos produced under in vitro culture conditions [1,2]. In addition, in vitro originated embryos are marked by alterations in their transcriptome profile compared to their in vivo counterparts [3][4][5][6][7][8] and this altered gene expression could be partly due to epigenetic reprogramming errors caused by aberrant DNA methylation [2].
DNA methylation is believed to be one of the mechanisms involved in regulating the gene expression profile of the embryo during its subsequent development. DNA modifications resulted from the addition of a methyl (CH 3 ) group to the cytosine residue in CpG dinucleotides within the DNA sequence by the catalytic activity of DNMT enzymes [9,10] make the DNA to be less accessible to the transcription machinery, consequently hindering the gene expression [11].
Unlike differentiated somatic cells whose genomic methylation patterns seem to be stable, the DNA methylation pattern in germ cells and preimplantation embryos is changing dynamically to maintain cell reprograming during development [10]. For instance, prior to fertilization, the sperm genome is relatively methylated while the oocyte genome is hypomethylated [10,12], but the methylation level in the CpG islands may be greater in oocyte compared to sperm. In mouse, the paternal genome undergoes active genome-wide demethylation immediately after fertilization, while demethylation of the maternal genome happens in a sequential fashion [13]. On the other hand, another study indicated the presence of active methylation in the paternal and maternal genomes during embryonic development [14], while others reported stable methylation levels in the early cleavage stages of embryos [12].
In bovine, the DNA methylation pattern changes dynamically in a sex and stage dependent manner during preimplantation embryo development [15]. For instance, while a dramatic loss of DNA methylation occurs in the male pronucleus immediately after the union of gametes [16], higher DNA methylation appears at the 8-cell stage in the female embryos and at the blastocyst stage in male embryos [15].
Apart from sex and stage of development, the DNA methylation pattern of the embryo is also modified by the culture environment [17]. Previous reports indicated altered placental DNA methylation patterns in human pregnancy established after transfer of in vitro produced embryos [18]. Similarly, a study in rabbit embryos has also shown faster demethylation in in vitro embryos compared to their in vivo counterparts [19]. Altered DNA methylation patterns in the fetal liver due to changes in oocyte maturation media and embryo culture have also been reported [20]. Moreover, in vitro conditions during preimplantation period could cause epigenetic abnormalities in the placenta during the later stages of pregnancy [21]. All in all, previous data indicated that suboptimal culture conditions may induce aberrant DNA methylation patterns during the preimplantation stage and fetal development. However, the preimplantation embryonic stages which are sensitive to culture conditions to cause significant modifications of the DNA methylation patterns in the resulting blastocysts remained unclear. In our previous study, blastocysts derived from in vitro embryos transferred to the recipient animals before or after embryonic genome activation displayed altered transcriptome profiles compared to the blastocysts developed completely in the vivo conditions [22]. Therefore, we hypothesized that stage specific suboptimal culture condition during preimplantation embryo development could cause significant modifications of the DNA methylation in the resulting blastocysts. Therefore, the main objective of this study was designed to explore the global DNA methylation patterns of the blastocysts derived from embryos cultured under in vitro culture conditions in a stage specific manner before being transferred to the in vivo condition.

Animal handling, management and estrous synchronization
Simmental healthy heifers were selected for this experiment. All experimental heifers received similar total mixed ration, and they were kept in the same farm and housing conditions. Handling and management of experimental animals were adhered to the rules and regulations of the German law of animal protection. Moreover, the experiment was approved by the Animal Welfare committee of the University of Bonn with proposition number 84-02.05.20.12.075.

In vitro embryo production, embryo transfer and flushing
Oocytes were collected from bovine ovaries obtained from local abattoirs. For this, cumulusoocyte complexes (COCs) were aspirated from follicle size of 2 to 8 mm diameter. Good quality oocytes were matured and fertilized in vitro as described previously [7,8,23]. After 18 hrs of oocytes and sperm co-incubation, the cumulus cells were removed and cumulus-free presumptive zygotes were washed three times in CR1aa supplemented with 10% heat inactivated estrous cow serum. The embryos were then cultured in 400 μl of culture medium in four well dishes (Nunc, Roskilde, Denmark) covered with mineral oil at 38.7°C in 5% CO 2 in humidified air. Some of the zygotes were then transferred to the recipient animals and others were transferred at 4-cell or 16-cell and other group of zygotes in vitro cultured until blastocysts stage. Therefore, for this study five blastocyst groups were generated. The first (ZY), second (4C) and the third (16C) blastocyst groups were obtained from embryos developed in vivo after transfer of in vitro developed zygotes, 4-cell and 16-cell stage embryos, respectively to estrous synchronized recipients using transvaginal endoscopic tubal transfer technique and the blastocysts were collected at day 7 of the estrous cycle by uterine flushing. Estrous synchronization was performed as indicated previously in our lab [22,24]. The fourth (IVP) blastocyst group was obtained from zygotes culture under in vitro condition until the blastocysts stage. The number of zygotes, 4-cell or 16-cell stage embryos transferred and the blastocysts recovered from each embryo group is indicated in S1 Table. The fifth (VO) blastocyst group was obtained from embryos developed completely under in vivo conditions. For this, Simmental heifers were estrous synchronized, super ovulated and artificially inseminated using the semen collected from the bull whose semen was used for in vitro fertilization as described previously [8,22,24] and the blastocyst stage embryos were collected 7 days post insemination as described above by uterine flushing. A total of four pools of blastocysts each consisting of 10 blastocysts from ZY, 4C, 16C, IVP and VO groups were used for global DNA methylation analysis and additional three pools each consisting of 10 blastocysts were generated to validate the array result using bisulfite sequencing.

Genomic DNA isolation and DNA fragmentation
Genomic DNA was isolated from all blastocyst groups using the Allprep DNA/RNA micro kit (Qiagen) according to the manufacturer's protocol. The genome DNA was then fragmented using the MseI enzyme in the presence of 0.5 μl bovine serum albumin (100x), 10 μl NEB buffer 4 (10x) and spike-in controls at 37°C for 16 hrs followed by further incubation at 65°C for 20 min. After the end of incubation period, gDNA was precipitated using linear acrylamide (Ambion), sodium acetate (3M, pH 5.2) and ethanol. After repeated washing in 70% ethanol, the gDNA pellets were dissolved in 5 μl nuclease free water.

Adapter ligation and methyl sensitive enzymes digestion
The protocols regarding adapter/linker ligation and methyl sensitive enzyme cleavage have been previously described in Shojaei Saadi et.al [25]. Briefly, MseLig12 (100 μM) (5'-TAA CTA GCA TGC-3') and 0.5 μl MseLig21 (100 μM) (5-AGT GGG ATT CCG CAT GCT AGT-3') adapters were ligated at the fragmented gDNA sites in the presence of 0.5 μl 10× One-Phor-All Buffer PLUS (Pharmacia Biotech) and 1.5 μl nuclease free water. The reaction was incubated at 65°C for 20 min to initiate annealing and the temperature was decreased to 15°C at the rate of 1°C /min. Following this, 1 μl ATP (10 mM) and 1 μl T4 DNA Ligase (5U/μl) was added to the reaction and incubated for 16 hrs at 15°C. At the end of reaction, unmethylated genomic sites within the MseI-MseI region of the fragmented were cleaved by DNA FastDi-gest™ methyl-sensitive restriction endonucleases enzymes (HpaII, HinP1I and AciI). For this, a mix consisting of 5μl of FastDigest buffer (10X), 0.5 μl HpaII (10U/μl), 0.5 μl HinP1I and 34 μl nuclease free water was added to the reaction and incubated overnight at 37°C. Afterwards, a mix consisting of 0.5 μl AciI, 5.0 μl of FastDigest buffer 10X and 44.5 μl of water was added to the reaction and further incubated at 37°C for 4 hrs. The reaction was then terminated after inactivating the enzyme at 80°C for 10 min.

Confirmation of methyl sensitive cleavage
DNA cleavage by methyl sensitive enzyme was confirmed by analyzing the spike-in controls added during gDNA fragmentation step as previously described in the EmbryoGENE DNA Methylation Array (EDMA) Platform [25]. For this, qPCR was preformed by preparing separate mixes consisting of 0.5 μl forward and 0.5 μl reverse primers (10 μM) (targeting HpaII, HinP1I or AciI sites), 1.6 μl MgCl 2 (25 mM), 2μl LC master Fast start mix (10x) and 14.4 μl Nuclease free water. Undigested spike-in (1/1000 dilution) and non-template controls were used as positive and negative controls, respectively. The amplification plot and the dissociation curve were evaluated for each of the HpaII, AciI or HinP1I primer sets. The differences in threshold cycle (Ct) between protected (hypermethylated) and unprotected (unmethylated) control templates were determined and uniformity of the samples was evaluated based on the qPCR amplification curves. Samples with the threshold cycle (Ct) ! 5 compared to undigested controls and displaying uniformity with other samples cleaved > 97% in the same cohort were used for downstream treatment. Insufficiently cleaved samples were subjected to a second treatment with specific methyl sensitive restriction enzyme. ethidium bromide stained agarose gel electrophoresis to confirm the distribution of the digested gDNA. Once the presence of several bands were observed, further amplification was performed using 0.75 μl of amplified DNA to which a mix consisting of 0.75 μl 100 μM Mse-Lig21, 5μl PCR buffer 1 Roche (17.5mM MgCl2), 4μl dNTPs, 0.8 μl Roche DNA long template enzyme mix and 38.7μl water was added. Amplification reaction was performed in four steps, in first step samples were incubated at 94°C for 60 min, 65°C for 30 s and 72°C for 2 min followed by 14 cycles at 94°C for 40 s, 65°C for 30 s and 72°C for 90 s. The reaction continued for 9 cycles at 94°C for 40 s, 65°C for 30 s and 72°C for 2 min and finally terminated after incubating at 72°C for 5 min. The PCR products were then purified using PCR purification kit (Qiagen) and finally eluted in 30 μl of 1/10 buffer elution buffer. Following this, the gDNA was incubated in the presence of 10 μl of buffer 4 NEB (10X), 1 μl BSA (100x), 1 μl Msel (10 u/μl) and 48 μl water at 37°C followed by incubation at 65°C for 20 min. At the end of reaction, the samples were purified using PCR purification kit (Qiagen) and finally eluted in 30 μl of 1/10 elution buffer. The quantity and quality of amplified digested gDNAs were measured using NanoDrop 1000 spectrophotometer.

Dye labeling and hybridization
Two micrograms amplified gDNA from each blastocyst sample was labeled with either Cy-3 or Cy-5 dyes using Universal Linkage System Fluorescent gDNA labelling kit (Kreatech Biotechnology) at 85°C for 30 min. The samples were then purified using the Qiagen PCR purification kit (Qiagen). The dye incorporation and DNA concentration were measured using NanoDrop 1000 spectrophotometer. Hybridizations were performed by mixing 1 μg of labeled gDNA sample of ZY, 4C, 16C or IVP, 1 μg of labeled gDNA of the VO blastocyst group, 25 μl of bovine Cot-1 DNA (1.0 mg/mL, Applied Genetics Laboratories), 2.6 μl of Agilent 100x Blocking Agent and 130 μl of Agilent 2x HI-RPM Hybridization Buffer (Agilent Technologies). The samples were incubated at 95°C for 3 min, at 37°C for 30 min and then 65 μl of Agilent-CGHBlock was added to each of the samples. The samples were loaded onto the microarray and hybridization was preformed in a hybridization oven (Shel Lab) for 40 hrs at 65°C and 20 rpm. Four biological replicates of ZY, 4C, 16C or IVP blastocyst group were hybridized to the four biological replicates of VO group in EmbryoGENE DNA Methylation Array EDMA which harbors nearly 2.3 M CpG sites corresponding to 10% of the bovine CpG sites in the bovine genome [25]. At the end of hybridization, the arrays were washed and scanned with the PowerScanner (Tecan) and analyzed with Array-Pro Analyzer 6.3 software (MediaCybernetics). Array features were extracted with the Agilent's Feature Extraction software (Agilent Technologies, CA, USA). A total of 16 arrays were hybridized and biological dye-swaps were performed for each sample. The raw and normalized data of 16 arrays have been deposited on the GEO public functional genomics data repository with GEO accession number GSE69173.

Array data and downstream analysis
The array data and downstream analyses methods used for the study have been described in EmbryoGENE Microarray (EDMA) platform [25]. Briefly, the intensity cutoff (Mean + 4 Ã SD) of the negative controls was used to identify probes that displayed a signal intensity above the background. Loess normalization followed by quantile inter-array scale normalization was performed to obtain Bayesian statistics of differential methylation. Differentially methylated probes were then identified from the full set of probes using linear models for microarray data (limma) [27]. Since in EDMA platform, the measured fold changes are not quantitative but merely indicative of higher or lower odds of differentially methylation, probes which showed significant (p < 0.05) differences with absolute log 2 (fold-change) !1.5 between the treatment and the reference sample were considered as differentially methylated regions (DMRs). Therefore, probes which showed a significant increase in signal intensity by 1.5 folds and higher in ZY, 4C, 16C or IVP compared to the VO group were considered as hypermethylated probes while probes which displayed a reduced signal intensity by 1.5 folds and higher in ZY, 4C, 16C or IVP compared to the VO group were considered as hypomethylated probes. Enrichment analysis for genome-scale DNA methylation data was performed using string of integrated scripts that sorts the genome data into CpG island density, CpG island length, CpG island distance, genomic location and types of repetitive elements. Moreover, the gene ontology enrichment analysis was performed for genes that were associated with differentially methylated regions. Probes were associated with the GO terms of a gene provided the probe fragment falls within the introns, exons or promoter of a particular gene. Following this, Go enrichment was assessed using the topGO bioconductor package (http://www.bioconductor.org/packages/ release/bioc/html/topGO.html) and the weight01 algorithm defined by Alexa et al [28]. For comparative analysis of methylation profile and gene expression data, the genes absolute log 2 (fold-change) !1.5, p value < 0.05 and false discovery rate (FDR) < 0.3 were selected from our previous data [22]. The heatmaps showing the relationship between differentially methylated regions and differentially expressed genes were constructed using PermutMatrix [29].

Validation of differentially methylated genomic regions using bisulfite sequencing
Twelve differentially methylated regions (DMRs) were selected to validate the array data independently using bisulfite sequencing. For this, primers (Table 1) were designed using MethPrimer (http://www.urogene.org/cgi-bin/methprimer/methprimer.cgi), an online tool for designing primers for amplification of bisulfite converted DNA. The genomic DNA in ZY, 4C, 16C, IVP and VO blastocyst groups was bisulfite converted using EZ DNA methylation direct kit (ZymoResearch) according to the manufacturer's instruction. PCR amplification was performed in 25 μl volume containing 2 μl of bisulfite-converted DNA, 0.4 μl forward and reverses primers (10μm), 0.4 μ1 dNTP mixes (25mM), 0.2 μl Zymo Taq™ DNA polymerase (5u/μl) (ZymoResearch). The presence of PCR product was confirmed by loading 5 μl of the PCR product onto electrophoresis run on 2% w/v agarose gel stained with ethidium bromide at 120 V for 20 min. The PCR product was then purified using QIAquick PCR Purification (Qiagen) and cloned to pGEM 1 -T Easy Vector Systems (Promega, WI, USA) and transformed to E. coli competent cells. The bacterial culture was plated onto LB agar/ampicillin/IPTG/X-gal plate and incubated overnight at 37°C. Following this, 10-20 independent white colonies were selected and sequenced in GenomeLab™ GeXP Genetic Analysis System (Beckman Coulter). The bisulfite sequencing DNA methylation analysis software (BISMA) [30] was used to analyze the sequencing data. The average of methylation percentage at each CpG site in each experimental group was used for analysis.

Confirmation of candidate genes using qPCR
In the current study, the DNA methylation data and previously published transcriptome profile were merged to get insight into the relationship between methylation pattern and the gene expression pattern. Although, the transcriptome array data was validated previoulsy, here we randomly selected 7 differentially expressed genes to confirm their expression pattern using the RNA samples isolated from the same samples whose DNA was used for methylation analysis. For this, the total RNAs isolated from each group of blastocysts were reverse transcribed using Universal cDNA synthesis kit (Exiqon) according to the manufactures recommendation. Briefly, a 20 μl reaction mix consisting of total RNA, reaction buffer, enzyme mixes and nuclease-free water was incubated at 42°C for 60 min, followed by 5 min incubation at 95°C. Gene specific primers (S2 Table) were designed and their PCR product specificity was confirmed sequencing the PCR products in GenomeLab™ GeXP Genetic Analysis System (Beckman Coulter). The qPCR was then performed in 20 μl reaction volume containing iTaq SYBR Green Supermix with ROX (Bio-Rad laboratories Germany), the cDNA samples of ZY, 4C, 16C, IVP or VO blastocyst group, the specific forward and reverse primer in the StepOnePlus™ Real-Time PCR Systems (Applied Biosystems, Foster city, CA). The qPCR thermal cycling parameter was set as 95°C for 3 min followed by 40 cycles of 95°C for 15 s and 60°C for 1 min. Following this, dissociation curve was generated by starting fluorescence acquisition at 60°C and measurements were taken every 7 s interval until the temperature was reached 95°C. The qPCR was performed in four replicates and analyzed using comparative (Ct (2 −ΔΔCt )) method. The specificity of amplification was evaluated by monitoring the dissociation (melting) curve of each candidate gene. Glyceraldehyd-3-phosphate-dehydrogenase (GAPDH) was used as endogenous normalizer.

Results
Differentially methylated genomic regions in blastocysts developed from embryos exposed to in vitro culture in stage specific manner In this study, we investigated the DNA methylation landscape of blastocysts developed in vivo after embryos were cultured in vitro until zygote (ZY), 4-cell (4C) or 16-cell (16C) stage and  blastocysts developed completely under in vitro condition (IVP) with reference to the blastocysts developed completely under in vivo condition (VO). For this, the blastocysts genomic DNA was digested using the MseI restriction enzyme followed by adapter ligation, methyl-sensitive restriction enzymes and ligation mediated PCR genomic amplification. The abundance of each fragment obtained from ligation mediated amplification in each sample was then measured using EmbryoGENE DNA Methylation Array (EDMA). Probes with signal intensities higher than the background signal in each comparison (ZY vs. VO, 4C vs. VO, 16C vs. VO and IVP vs. VO) were analyzed for further analysis. Accordingly, from a total of > 400K probes, the signal intensity of 196871, 199265, 255254, 248453 and 227075 probes in ZY, 4C, 16C, IVP and VO groups, respectively exceeded the background signal and 174460 probes were commonly detected in all sample groups (S1 Fig). Differentially methylated regions were identified by fitting intensity data using a linear model [27]. Probes that exhibited fold-changes of ! log 2 1.5 at p < 0.05 in each comparison were considered differentially methylated regions (DMRs). Thus, differential methylated region analysis indicated that longer exposure of preimplantation embryos to the in vitro culture conditions increased DNA methylation profile alteration in the resulting blastocysts (Fig 1). A systematic increase in both hypermethylated and hypomethylated genomic regions were identified from the ZY to the IVP group, with hypermethylated loci outpacing that of the hypomethylated ones in 4C and 16C blastocyst groups. Furthermore, 42, 75, 69 and 51.8% of the total DMRs in the ZY, 4C, 16C and IVP blastocyst groups, respectively were hypermethylated and the rest were hypomethylated. Apart from these, the DNA methylation profiles identified in ZY, 4C, 16C and IVP blastocyst groups were evenly distributed across all chromosomes without showing any significant chromosome specific effects (

Stably hypermethylated or hypomethylated genomic loci
In order to get insight into the genomic loci that were differentially methylated in a stable or stage dependent manner, we looked into the DMRs that were specific to each blastocyst group and the DMRs that were common in two or more blastocyst groups. Hypermethylated or hypomethylated genomic loci which appeared in two or more blastocyst groups were considered as "stably" differentially methylated regions. Accordingly, a total of 137, 624, 1180, and 3086 DMRs were found to be specific to only ZY, 4C, 16C and IVP blastocyst groups, respectively, whereas 115, 350, 65 and 348 DMRs in ZY, 4C, 16C and IVP blastocyst groups, respectively were stably differentially methylated in two or more blastocyst groups (Fig 3). Moreover, differentially methylated genomic loci specific to each group were increased depending on stages of the embryos while stably differentially methylated genomic loci didn't exhibit that trend. Among stably hypermethylated or hypomethylated genomic loci, 45 DMRs were identified in ZY, 4C and IVP groups, but these genomic loci were not significantly differentially methylated in the 16C blastocyst group (Fig 3, S3 Table). Furthermore, we look into a stage wise comparisons from ZY to IVP, it was evidenced that the ZY blastocyst group had relatively higher DMRs in common with 4C group than it did with other blastocyst groups (Fig 3, S3 Table). On the other hand, the 4C blastocysts shared relatively higher number of DMRs (n = 290; 29.7% of the total) with IVP blastocyst group and thus compared to other groups, in terms of common DMRs, the 4C blastocyst group was quite closer to the IVP blastocyst group (Fig 3, S5 Table) followed by ZY group. However, the 16C blastocyst group, which the EGA is believed to occur in vitro, shared only 3% of the DMRs (45/1245) with IVP blastocyst group (Fig 3).
Differentially methylated CpG islands in blastocysts of different groups. In EDMA platform, the CpG islands covered about 13% of the total probes. We used this opportunity to analyze differentially methylated CpG islands in each blastocyst group. The result has shown that the methylation pattern of CpG islands was found to be higher in IVP (n = 357) compared to ZY (n = 14), 4C (n = 35) and the 16C (n = 22) blastocyst groups. Furthermore, the data also evidenced that 42, 86, 71, and 50% of the differential methylated CpG islands in ZY, 4C, 16C and IVP groups, respectively were found to be hypermethylated. In addition, we also looked into whether the methylation patterns of CpG islands that fall within the promoter, exon or intron genomic region could differ between the blastocyst groups. Indeed, the genomic localization of these genomic loci indicated that in ZY group, differentially methylated CpG islands were more frequently at intron ( Table 2) whereas in the 4C group, differentially methylated CpG islands were more frequently detected at exon (Table 3). Furthermore, in the 16C blastocyst group, the frequency of differentially methylated CpG islands was found to be similar in exonic and intronic regions (Table 4). Similarly, in IVP groups, the proportion of differentially methylated CpG islands detected on exonic region exceeded those located at the intron and promoter regions (S6 Table) although a total of 28 differentially methylated CpG islands were identified in promoter and proximal promoter regions ( Table 5).
Characterization of the differentially methylated CpG islands with respect to their length. In the EDMA probes, the lengths of CpG islands were partitioned as a function of percentiles in which the bottom 20 and the top 80 percentiles were classified as short and long CpG islands, respectively and the CpG islands between the bottom 20 and the top 80 percentiles were classified as intermediate islands. Thus, using this option, we determined the methylation pattern of short, long and intermediate CpG islands in blastocysts of different groups. The result showed that considerable differences in proportion differentially methylated small, intermediate and long CpG islands were observed between the blastocyst groups ( Fig 4A). The small CpG islands were hypermethylated in ZY, 4C and 16C groups, while long and intermediate-length CpG islands were hypermethylated in the 16C blastocyst group. However, in IVP group, as much hypomethylation as hypermethylation were detected in long or intermediate CpG islands (Fig 4B). Indeed, differentially methylated CpG islands covered only 2-10% of the total DMRs while 90-98% of the DMRs were located outside the CpG islands ( Fig 4C). Therefore, the non-CpG island DMRs were classified into CpG shores, Open Sea and CpG shelf when the probe is located 1-2 kbp, 2-4 kbp and > 4 kbp away from the nearest CpG island, respectively. Thus, distribution of DMRs in CpG shore, Open Sea and CpG shelf tended to be   CpG_L is the length of the island on the genome while CpG-I indicates the number of bases within the fragment which are part of the CpG island.
Numbers preceded by hyphen following the gene name indicate the exon or intron numbers depending whether the DMR is located in exon or intron of the gene. EDMA_ID refers to the identification ID of each EDMA probe. Each ID is preceded by EDMA_MET. Ppromoter = proximal promoter.
doi:10.1371/journal.pone.0140467.t002 different between the blastocyst groups ( Fig 4C). Moreover, the proportion of differentially methylated region in CpG shore, Open Sea and CpG shelf in IVP blastocyst group exhibited a similar trend to that the proportion found in the whole array (EDMA). However, the expected proportion of DMRs in CpG shore, Open Sea and CpG shelf was deviated in ZY, 4C and 16C blastocyst groups compared to the proportion found in the whole array (EDMA) indicating the presence of selective differential methylation of genomic loci adjacent to the CpG islands depending on the stage of the embryos exposed to suboptimal culture condition (Fig 4C). Moreover, unlike to ZY and IVP groups, hypermethylated genomic loci were increased in CpG shore, Open Sea and CpG shelf in 4C and 16C blastocyst groups (Fig 4D).
Differentially methylated repetitive elements. The repetitive elements, namely the longterminal-repeat (LTR) retrotransposons, short-interspersed repetitive elements (SINEs), longinterspersed repetitive elements (LINEs), low-complexity repetitive elements and simple CpG_L is the length of the island on the genome while CpG-I indicates the number of bases within the fragment which are part of the CpG Island.
Numbers preceded by hyphen following the gene name in the "Genes containing DMRs" column indicate the exon or intron number depending whether the DMR is located in exon or intron of the gene. EDMA_ID refers to the identification ID of each EDMA probe. Each ID is preceded by EDMA_MET.
PPromoter. = proximal promoter. doi:10.1371/journal.pone.0140467.t003 repeats were also differentially methylated in in ZY, 4C, 16C and IVP blastocyst groups. The number of differentially methylated SINE, LINE, LTR and low complexity repeats were relatively higher in IVP followed by 4C and 16C blastocyst groups ( Fig 5A). Moreover, the majority of these classes of bovine repetitive elements tended to be hypomethylated in ZY and hypermethylated in 4C, 16C and IVP blastocyst groups (Fig 5B).

Genomic localization of differentially methylated regions
Apart from characterizing with respect to the CpG islands, we analyzed the distribution of all the DMRs with respect to their genomic position relative to the nearby genes. This analysis indicated that 65, 73, 71 and 75% of the DMRs in ZY, 4C, 16C and IVP, respectively were located in and around the promoter or gene body regions (Fig 6A). Moreover, genomic enrichment analysis revealed that in the ZY blastocyst group, the hypermethylated probes tended to be enriched only in the intronic regions whereas in the 4C and 16C groups, hypermethylated regions were enriched in intron, exon and promoter regions. On the other hand, in the IVP group, as much hypermethylation as hypomethylation were detected in intronic, exonic and promoter regions (Fig 6B). CpG_L is the length of the island on the genome while CpG-I indicates the number of bases within the fragment which are part of the CpG Island.
Numbers preceded by hyphen following the gene name in the "Genes containing DMRs" column indicate the exon or intron number depending whether the DMR is located in exon or intron of the gene. EDMA_ID refers to the identification ID of each EDMA probe. Each ID is preceded by EDMA_MET.

Functional annotation of genes affected by differentially methylated regions
To get insight into the functional elements of the genomic region affected by the changes in methylation profile, gene ontology (GO) enrichment of the DMRs was performed using the topGO bioconductor package along with the algorithm defined by Alexa et al [28]. Probes which are located within the gene body or within 4 kbp of its transcription start site were considered for GO terms analysis. Accordingly, protein homooligomerization and ATP binding functional categories were affected in the ZY group (Fig 7A). On the other hand, at least 9 GO terms including ATP binding, programmed cell death and ATPase activity in the 4C blastocyst group (Fig 7B) and 4 gene ontology terms including glycolysis and nucleoside triphosphatase activity in the 16C blastocyst group (Fig 7C) were found to be affected by DMRs. In the IVP blastocyst group, the DMRs were associated with 83 gene ontology terms including genetic imprinting, regulation of transcription, metaphase/anaphase transition, chromosome segregation, cholesterol homeostasis and tricarboxylic acid cycle (Fig 8).

Comparative analysis of DNA methylation and transcriptome profile
To determine whether the methylation profile changes in a given loci were associaited with changes in gene expression pattern, the DNA methylation patterns of ZY, 4C, 16C and IVP blastocyst groups were overlapped to our previous transcriptome data generated from the same blastocyst groups [22]. The analysis indicated that among the DMRs that had transcriptome profile information, 4.3% (n = 11) in ZY, 3.5% (n = 35) in 4C, 4.1% (n = 52) in 16C and 5.6% (n = 191) in IVP group were overlapped to the differentially expressed genes (DEGs) while the rest were associated with genes that were not significantly differentially expressed. Considering that a higher methylation level should relate to a lower gene expression pattern [11], we opted to identify the DMRs that are negatively correlated with their corresponding On the other hand, in the 4C blastocyst group, 19 DMRs (3.4% of the DMRs overlapped to the transcriptome data) including those located at TRAPPC9, ISM1, GLDC and JAK1 genes were found to be inversely correlated with the gene expression profile (Fig 9B), while the methylation pattern of 17 DMRs was positively correlated with the gene expression (S5B Fig). Similar trend was observed in the 16C blastocyst group of which 3.9% of the total DMRs mapped to the transcriptome data were found to be inversely related to the expression pattern of their target genes (Fig 9C), but the mRNA levels of genes associated with 14 hypermethylated and 11 hypomethylated regions were positively correlated to the DNA methylation pattern (S5C Fig). In IVP blastocyst group, a total 107 DMRs (9.4% of the total DMRs mapped to the transcriptome data) comprising genes involved in cell death (ATXN3, UNC5B, NUAK2, BCL2L13, FGD3, PEG3) and cell division (STX2, ZWINT, ZNF830, BUB1B, NUP37, PPP1CC), displayed a negative correlation with the corresponding gene expression (Figs 10 and 11). However, the expression levels of genes associated with 24 hypermethylated genomic loci were increased and those associated with 60 hypomethylated DMRs were reduced compared to the in vivo blastocyst group (S6 Fig). Validation of differentially methylated region using bisulfite sequencing We selected 12 differentially methylated regions located either within or outside of CpG islands to confirm the array data using bisulfite sequencing. Three of the DMRs, namely EDMA_MET_18_16952, EDMA_MET_04_15002 and EDMA_MET_01_04641 associated with SYT5, PTPRN2 and GPR156 genes, respectively were commonly detected in two or more blastocyst groups. Among the probes analyzed, EDMA_MET_18_16952, which is a corresponding to the SYT5 gene was found to be hypomethylated both in 4C and IVP compared to the VO group and the bisulfite sequencing results confirmed the methylation level of this locus to be lower in both 4C and IVP compared to the VO group. Similarly, both bisulfite sequencing and the array results confirmed a higher methylation level of the fragment corresponding to the EDMA_MET_04_15002 (PTPRN2 gene) in 16C and IVP and lower methylation level of the genomic fragment that corresponds to the EDMA_MET_13_09053 probe (associated with CTSZ gene) in IVP compared to the VO group. Therefore, the results of the bisulfite sequencing were by and large similar to the array data (Fig 12, S7 Fig). Validation of differentially expressed genes using quantitative polymerase chain reaction (qPCR) Since comparative analysis of DNA methylation pattern and gene expression of the blastocyst groups was performed by superimposing the current DNA methylation data and our previouly published transcriptome profile data [22], we sought to confirm the expression pattern of selected candidate differentially expressed genes in RNA samples isolated from blastocysts which were used for DNA methyltion analysis. The result indicated that the expression profile of 6 of the 7 selected genes was found to be in line to the array data (Fig 13).

Discussion
Most often, in vitro produced embryos exhibited a reduced quality and viability compared to their in vivo counter parts. The molecular data also indicated compromised gene transcription, imprinting gene maintenance errors and repression of transposable elements in in vitro produced embryos due to altered epigenetic regulation such as aberrant DNA methylation [31,32]. However, it has been a challenge to identify the embryonic stages that could activate aberrant DNA methylation patterns when embryos are cultured under suboptimal culture conditions. Here, we investigated the genome wide DNA methylation patterns of blastocysts developed in vivo after preimplantation embryos were developed in vitro until zygote (ZY), 4-cell (4C) or 16-cell (16C) stages. The blastocysts developed completely under in vitro culture condition (IVP) were also included in the study to understand the effect of continuous in vitro culture on the DNA methylation pattern of the blastocyst stage bovine embryo. The DNA methylation pattern of all blastocyst groups were determined using blastocysts developed completely under in vivo condition (VO) as a common reference. Interestingly, we have identified a higher DNA methylation dysregulation in the IVP blastocysts followed by the 16C and 4C blastocyst groups, and relatively lower DNA methylation dysregulation was detected in the ZY blastocyst group. Therefore, the extent of DNA methylation dysregulations of the blastocysts seem to be strongly associated with the successive developmental stages completed under in vitro condition (Fig 1). In line to this, our previous transcriptome profile data has also indicated a higher gene expression dysregulation in IVP and 16C blastocyst groups than the 4C blastocyst group [22]. However, the critical question is: why differentially methylated genomic loci were increased in the resulting blastocysts when the preimplantation embryos advanced their development stages under in vitro culture conditions? One hypothesis is that, any stage development the embryo undergoes may be associated with the induction of a fixed level of methylation dysregulation, and early transferred embryos to the in vivo condition could have relatively more time to correct that dysregulation and revert to the normal state. Therefore, culturing embryos until blastocyst stage in vitro could induct dysregulation of DNA methylation profile in many more genomic loci than those blastocysts derived from embryos cultured in vitro only until certain specific stages of development. Despite differences in the stage of development in which the embryos spent under in vitro environment, certain genomic regions were commonly differentially methylated in two or more blastocyst groups (Fig 3). For instance, although, the ZY; 4C and 16C blastocyst groups were derived from preimplantation stages embryos exposed to in vitro culture until certain stages, 26% (66/253), 29.8% (290/ 974) and 3% (45/1245) of the DMRs in ZY, 4C and 16C blastocyst groups were found to be detected in the IVP blastocyst group indicating relatively higher communality between the 4C and IVP blastocyst groups. In fact, the 16C blastocysts were derived from embryos developed in vivo after transfer of in vitro produced 16-cell embryos, which the EGA occurred in vitro, intuitively, higher commonality would be expected between the 16C and IVP than the 4C and IVP blastocyst groups. This may in turn suggest that reprogramming errors occurred due to suboptimal culture condition before the major embryonic genome activation may prevail until the blastocyst stage although major embryonic genome activation occurs under in vivo condition. Our previous transcriptome profile data also indicated that higher communality in number of differentially expressed genes between blastocyst of the IVP and the 4C groups than IVP and 16C [22]. Then the core question is why the 4-cell transfer group tended to share relatively higher number of DMRs with IVP than the 16C blastocyst group did? Although further research is required to deduce a conclusion, this can be associated with the epigenetic memory that could be remained inscribed due to suboptimal environmental condition before embryonic genome activation. In bovine embryos, although minor embryonic genome activation (EGA) may start at 2-to 4-cell stage [33][34][35], major EGA Altered DNA Methylation Patterns in Bovine Blastocysts is believed to occur at the 8-to 16-cell stages [14,[36][37][38][39]. For instance, Kues et al. [39] indicated increased expression level of about 2400 genes in 8-cell embryos compared to 2-cell and/or 4-cell stages and Jiang et al. [38] reported also significantly higher number of differentially expressed genes as the embryo grows 4 to 8-cell stages suggesting in both cases the occurrence of major EGA at the 8-cell stage. In fact, activation of the maternal genome and transcription of the genome messages could be partly control by the pattern of DNA methylation [40]. In this regards in vitro studies indicated that the bovine embryos undergo global DNA methylation erasure soon after fertilization until 8 cell stage and an increase in methylation level occurs at 16-cell until blastocyst stage [13,15,41]. Therefore, we can speculate that some reprogramming errors that could occur before major embryonic genome activation could be displayed at the blastocyst stage although the embryos were set into the optimal culture condition.
In addition to looking into the general trends of the differentially methylated genomic regions, we have also sought to explore the genomic features whose methylation patterns were affected by the in vitro culture condition. To this end, we have also analyzed the differentially methylated CpG islands, which are believed to be the evolutionary landmarks of molecular events occurring at eukaryotic promoters [42]. Indeed, along with other genomic features, the CpG islands were also differentially methylated in all blastocyst groups although the ZY, 4C and 16C blastocyst groups were relatively less affected compared to the IVP group. Since differences in the length of the CpG islands can be associated the functional changes, we have also explored the characteristics of the differentially methylated region with respect to its CpG island length, [43]. In this regard, our result indicated differences in the proportion of long, intermediate and short CpG islands between the blastocyst groups of which the frequency of small length differentially methylated CpG islands were dominant in the 4C blastocyst group, while intermediate length CpG islands were dominant in ZY, 16C and IVP blastocyst groups (Fig 4B). This may suggest that the DNA methylation pattern of genomic regions dominated by short or long stretches of CpG sites in the blastocysts could be determined by the stage of development the embryos completed under in vitro culture condition.
Exploring the genomic region of differentially methylated CpG islands and examining the effect of these genomic features on gene expression of nearby genes may spark light of elimination to get an insight into the effect of suboptimal embryo culture condition on the embryo methylene. Indeed, our data evidenced that unlike to the ZY and IVP groups, the majority of differentially methylated CpG islands in 4C and 16C blastocyst groups were hypermethylated. However, only few hypermethylated CpG islands were located on the promoter region of these blastocyst groups. Although most CpG islands are believed to be methylation free [44], parts of the CpG islands could be methylated during development and gene silencing may occur if this happens at the promoter region [41,[45][46][47][48][49]. Nevertheless, in the current study, there was a hint that the gene regulation by CpG islands is not limited only at the promoter regions, but also this can happen at the intronic and exonic regions (S7 Table). On other hand, the hypo or hypermethylated CpG islands may not always associated with repression or induction of gene expression (S8 Table). We have also assessed the location of the non-CpG-island DMRs with respect to the nearest CpG islands to get an overview on the methylation pattern of the genomic regions which are adjacent or far from the CpG islands (Fig 4C and 4D). This type of analysis has been suggested to be useful to correlate the DNA methylation pattern with the corresponding phenotype changes [50]. In the present study, enrichment of differentially methylated regions in CpG shore, Open Sea or CpG shelf tended to be different between the ZY, 4C, 16C and IVP blastocyst groups. For instance, more than 2 fold enrichment of hypermethylated than expected was detected at the CpG shore in 4C and CpG shelf in the 16C blastocyst groups. Indeed, it is suggested that altered DNA methylation profile particularly at CpG island shores is not only regulate the gene expression patterns [51], but also it affects the cell differentiation and development [52]. Generally, characterizing differentially methylated regions with respect to the CpG islands could provide insights into the presence/absence and the level of aberrant DNA methylation at the CpG island and its vicinity. Moreover, describing differentially methylated regions in relation to proximal promoter, promoter, intron and exon is essential to better understand the embryo genome structures subjected to aberrant methylation patterns due to suboptimal culture condition. In this regard, the current study has evidenced the presence of strong difference between the ZY, 4C, 16C and IVP blastocyst groups with respect to genomic enrichment of hypermethylated genomic region. For instance, as high as four fold enrichment of hypermethylated genomic region were detected at proximal promoter and exonic region in 4C group and promoter and proximal promoters in the 16C blastocyst group. However, this was not the case in ZY blastocyst group, of which hypermethylation was favored and in the IVP blastocyst group, of which the ratios of hypermethylated and hypomethylated region in gene body categories, promoter and proximal promoters were close to 1 (which is Log 2 enrichment = 0) (Fig 6). Here the fundamental question is why blastocysts developed in vivo from in vitro produced 4-cell and 16-cell were associated with hypermethylation whereas the majority of genomic features in blastocysts developed in vivo from in vitro produced zygotes were enriched by hypomethylation and as much hypermethylation as hypomethylation were detected in blastocysts developed completely under in vitro condition in several genomic features? Although, detail and locus specific DNA methylation study may be required to uncover this puzzle, previous studies based on immunoreactive 5-methylcytosine antibody on in vitro culture bovine embryos indicated demethylation of the paternal genome in zygotes and a declining methylation pattern from cleavage until 8-cell stage and increasing methylation level in 16-cell until blastocyst stage [12,13,15].
Parallel to DNA methylation profiling, we have also compared the gene expression data [22] and the current methylation profile to get insight into the association of methylation profile and gene expression patterns in these blastocyst groups. The data revealed that 1.6%, 3.4% 3.9% and 9.4% of the total DMRs mapped to the transcriptome data in ZY, 4C, 16C and IVP, respectively displayed a negative correlation with the corresponding gene expression patterns and the rest exhibited unexpected patterns. Thus, the relationship between the DNA methylation pattern and mRNA expression level was not always as expected provided the negative correlation between gene expression and DNA methylation is the expected trend. In fact, from the DNA methylation pattern of bovine placenta, previously it has been suggested that the relationship between gene-body DNA methylation and gene expression is non-monotonic [53]. Indeed, in several occasions it has been noted that the DNA methylation is negatively correlated with gene expression at the promoter region [45,54,55] and positively correlation at the gene body region [55,56]. However, in our result, it is indicated that some of the gene body methylation were negatively and some others were positively correlated with the gene expression and the same was true for promoter regions (Figs 9, 10 and 11, S5 & S6 Figs). Thus, it is less likely to generalize that methylation in promoter regions typically leads to decreased gene expression while DNA methylation in gene bodies is often correlates with increased gene expression. In line to our study, the study by Yang et al. [57] has also indicated a negatively correlation between some gene body DNA methylations and gene expression. These authors also added that the methylation pattern of only 25 out of 105 probes at the promoter region displayed a significant inverse correlation with the gene expression stressing that even methylation in promoter region may not always result in a decrease in gene expression. Thus, the relationship between gene body/promoter methylation and gene expression may vary depending on the physiological status and the cell type under consideration. On the other hand, looking at single base resolution level may be essential to understand the association between methylation level and gene expression more precisely. For instance, it was reported that non-CG methylation level within the gene body to be positive correlated with the gene expression and however, there was no correlation between CG methylation density and gene expression [58]. However, the main issue should rather focus on extent/ proportion of DMRs that positively or negatively correlated with gene expression. Indeed in our study, it was indicated that only 4.3%, 3.5%, 4.1% and 5.6% of the total DMRs in ZY, 4C, 16C and IVP blastocyst groups were positively or negatively correlated with the corresponding gene expression and the majority of the DMRs didn't show changes in the gene expression patterns. Nevertheless, it is also known that the DNA methylation is not only regulates the mRNA expression, but also the expression of other non coding RNAs including miRNAs [59][60][61]. Thus, the DMRs which were not associated with the altered gene expression in ZY, 4C, 16C and IVP blastocyst groups could be involved in regulating small non coding RNAs molecules. Similarly, we have also identified altered gene expressions in ZY, 4C, 16C and IVP blastocyst groups which were not directly associated with DNA methylation dysregulation. In this case, other key epigenetics effectors such as histone modification, chromatin remodeling and non-coding RNAs could be involved to fine tune the regulation of gene expression during embryo development under suboptimal environmental condition.

Conclusion
In this study, we have provided detailed insights into the DNA methylation landscape of bovine blastocysts developed from embryos cultured under in vitro condition in a stage specific manner before transfer to the in vivo condition along with blastocysts developed completely developed under in vitro condition. The data has shown remarkable differences in the extent of DNA methylation profile of blastocysts depending on developmental stages completed under in vitro culture condition. The DNA methylation patterns of CpG islands and repetitive elements were also affected in the blastocysts depending on the stage of the embryo subjected to the in vitro culture. Blastocysts developed in vivo after transfer of in vitro developed 4-and 16-cell stage embryos displayed increased proportion of hypermethylated genomic loci in the gene body, promoter or proximal promoter regions. Moreover, 1.6, 3.4, 3.9 and 9.4% of the DMRs which overlapped to the transcriptome profile data in ZY, 4C, 16C and IVP blastocyst groups, respectively were negatively correlated with their corresponding gene expression while the majority exhibited unexpected patterns given that DNA methylation is usually thought to negatively regulate gene expression. Therefore, the observed systematic DNA methylation profile alteration in the blastocysts is a clear demonstration the effect of stage specific in vitro culture condition on the embryo epigenome.