Integrated Analysis of Quantitative Proteome and Transcriptional Profiles Reveals the Dynamic Function of Maternally Expressed Proteins After Parthenogenetic Activation of Buffalo Oocyte*

Maternal-effect genes are especially critical for early embryonic development after fertilization and until massive activation of the embryonic genome occurs. By applying a tandem mass tag (TMT)-labeled quantitative proteomics combined with RNA sequencing approach, the proteome of the buffalo was quantitatively analyzed during parthenogenesis of mature oocytes and the two-cell stage embryo. Of 1908 quantified proteins, 123 differed significantly. The transcriptome was analyzed eight stages (GV, MII, 2-cell, 4-cell, 8-cell, 16-cell, morula, blastocyst) of Buffalo using the RNA sequencing approach, and a total of 3567 unique genes were identified to be differently expressed between all consecutive stages of pre-implantation development. Validation of proteomics results (TUBB3, CTNNA1, CDH3, MAP2K1), which are involved in tight junction and gap junction, revealing that the maternal expression of the proteins possibly plays a role in the formation of cellular junctions firstly after parthenogenetic activation. Correlation and hierarchical analyses of transcriptional profiles and the expression of NPM2 and NLRP5 mRNA of buffalo in vitro developed oocytes and parthenogenetic embryos indicated that the “maternal-to-zygotic transition” (MZT) process might exist in the model of parthenogenesis, which is similar to a normally fertilized embryo, and may occur between the 8-cell to 16-cell stage. These data provide a rich resource for further studies on maternal proteins and genes and are conducive to improving nuclear transfer technology.

Oocytes of all animals contain mRNAs and proteins that are deposited during oogenesis. The oocyte provides the specialized environment within which the zygotic nucleus initiates and executes its developmental program and as a result, also plays an obligatory role in its regulation. Maternal effects are pronounced in mammals especially. The uterine environment has a major effect on offspring development in the early stage, regulated by the mechanisms including hormone transmission (1) across the placenta and nutrient transfer (2). On the other hand, the postnatal environment is structured by the joint influence of milk quality and maternal behaviors, including time spent nursing (3,4), nest building, maintenance (5), and offspring grooming, and the quality of milk varies in its composition, such as fats, carbohydrates, and protein (6), and in biologically active proteins that can have effects on growth and immunological (7). To fully describe maternal effects in mammals, therefore, the identification of qualitative differences in mRNA transcript and protein content of maternal genome will enhance our understanding of the molecular genetic factors essential for maternally directed control of early embryonic development.
And maternal-effect genes were first defined in Drosophila 30 years ago, and in mammals in 2000 (8,9). Maternal genes are those genes whose products, RNA or protein, are produced or deposited in the oocyte or are present in the fertilized egg or embryo before expression of zygotic genes is initiated (10). Microarray expression studies have revealed numerous genes that are maternally expressed in vertebrates (11)(12)(13)(14)(15)(16)(17)(18)(19)(20). Maternal-gene products are particularly important in early embryos after fertilization occurs and up to the point when vast activation of the embryonic genome appears. The maternal factors that are stored in the oocyte greatly affect the maternal-to-embryonic transition, including the subcellular organelles and macromolecules, and they are acquired during folliculogenesis and oogenesis. Various technological approaches have been used to identify maternal-effect genes in mammals since 2000. Several maternal-effect genes have been identified in mice using typical knockout strategies (8,9,21,22). The conventional knockout method can lead to embryonic lethality, but conducted with conditional knockdown via oocyte injection and RNAi (the RNA interference pathway) can make precise exploration of maternal genes. In recent years, an increasing list of maternal factors involved in the regulation of MET has been identified. These factors include NLRP5 (official name for MATER) (9), HSF1 (8), NPM2 (21), ZARL (22), ZFP36l2 (23), BRG1 (24), PADI6 (25), CTCF (26), ATG5 (27), and AGO2 (official name for EIF2C2) (28).
Also, many maternal genes have been described whose products are involved in protecting imprinted methylation sites during preimplantation development, including DPPA3 (29), ZFP57 (30), DNMT1 (31), and TRIM28 (32). Although the importance of maternal effects (MEs) is becoming widely appreciated, we know little about dynamic function and interaction between maternally expressed proteins or genes at the level of protein and transcriptome.
DNA microarray research has provided valuable information for the dynamic changes of important factors at the transcription level. However, there is a large amount of material that is transcribed and stored in oocytes that has no apparent effect on oogenesis, but that requires the regulation of embryonic development. Thus, a precise picture of the functional status of oocytes is needed and can be achieved through transcriptomic and proteomic analyses. Also, the regulation of maternally expressed proteins from maternal genome remains poorly understood. How to obtain the maternal information only? To this end, we used parthenogenetic activation of oocyte as our research model because it does not require any paternal factors, which makes it an ideal model for the study of maternal genome. Also, parthenogenetic activation of mammalian oocytes is a very important research ingredient in reproductive biology. And it is not only provides alternative embryonic resources of embryos, but also is the most important technical condition for activation of nuclear transfer recombinant embryos in somatic cell nuclear technology; Parthenogenetic activation of mammalian oocytes also are used to establish embryonic stem cell lines.
Proteomic methods have been used in the exploration of mammalian oocytes and embryos, including mouse (33)(34)(35)(36)(37), bovine (38 -40), and pig (41)(42)(43)(44), and the amount of proteomics data has increased. Also, transcriptomic methods are also widely used in this area of research (45)(46)(47)(48). But these related studies were mainly focused on identifying protein expression profiles of oocytes and embryos at different developmental stages, however, to date, there is no report on the dynamic changes of maternally expressed proteins and genes from maternal genome (without paternal genome) that occur during the early stages of the embryo using the parthenogenesis as an experimental model at the level of transcription and translation for mammal.
Buffalo (Bubalus bubalis) is an important economic livestock species in Southern Asia. Buffalo was used as a tool for farming in the past, but with the rapid development of agricultural machinery, its advantages in farming have disappeared, so this species faces the risk of being eliminated. However, we are pleased to discover that buffalo have many good characteristics, such as resistance to roughage, the high value of dairy products, and good quality meat. However, its economic value is limited by its low reproductive efficiency. There are many factors which result in the low reproductive efficiency of Buffalo, for example, the beginning of suitable breeding age for buffalo is latter than cattle and cow; and the level of progesterone of female buffalo during its period of pregnancy is much lower than the cattle; and the characteristics of estrus are inconspicuous in female buffalo, and its interval of calving is longer than cattle and cow, also the female buffalo is not sensitive to the treatment of reproductive hormone. Notably, under the laboratory conditions, we found that the quantity of vesicular follicle and the oocyte suited to the IVM (in vitro maturation) collected from the buffalo ovary is less than the cattle and cow, and we found that the rate of cleavage and blastocyst in buffalo is lower than the cattle during its period of pre-implantation, which indicated that the different regulatory mechanism may exist in the development of pre-implantation, and before the activation of zygotic genome, the early development of embryo is regulated by maternal genome mainly.
The application of somatic cell nuclear transfer technology in buffalo reproduction has been reported (49 -52). However, there are still many problems in nuclear transfer technology, which limits the promotion of nuclear transfer technology, such as low cloning efficiency, low embryo implantation rate, abnormal fetal development, and many abnormal gene expressions during genome development of cloned embryos, including abnormalities of DNA methylation, acetylation, abnormal expression of imprinted genes, and so on. It is well known that the improvement of reproductive technology for the domestic animal must be based on the understanding of its reproductive mechanism and characteristics. However, the current theoretical knowledge about the reproductive mechanism of this species is still not enough, and the improvement of our breeding technology is limited by this. This study made the research of maternal genome as the entry point to study the reproductive mechanism of buffalo, and this knowledge will enhance the efficacy of in vitro maturation protocols utilized in either an agricultural or clinical setting and attempts to clone domestic animals by nuclear transfer. In addition, the strategy of this research may provide a valuable reference to explore the function of the maternal genome for other species.
In this study, we used the tandem mass tag (TMT) 1 -labeled quantitative proteomic combined with transcriptome of oocytes/embryos to research the dynamic function of the maternal genome during the early embryonic development. And mature oocytes (metaphase II [MII]) and two-cell stage embryos (parthenogenetic cleavage [PC]) were used as a quantitative proteomic object and the oocytes and embryos from 8 stages of buffalo were used as RNA sequencing-based transcriptomic object. The aim of this study was to show the dynamic function of maternally expressed proteins or genes using the research model of parthenogenetic activation of oocyte and to identify the important role of maternally expressed proteins in early embryo development. We also applied bioinformatic methods to an analysis of the protein expression profiles by Gene Ontology (GO) classification and a mapping of the interaction control network for partial proteins. Finally, we verified the expression of differential proteins at transcription and translation levels using real-time polymerase chain reaction (RT-PCR) and Western blotting, and we analyzed protein migration and expression using changes in immunofluorescence.

EXPERIMENTAL PROCEDURES
Collection of Buffalo In Vitro Developed Oocytes and Parthenogenetic Embryos-The collection and maturation of oocytes were performed as previously reported (53). In brief, we collected buffalo ovaries from the local slaughter house and transported them to the laboratory at 25°C in physiological saline. We washed the ovaries three times in saline solution, extracted cumulus-oocyte complexes (COCs) from 2 to 6 mm follicles in the ovaries, and selected COCs with compact cumulus-cell layers for maturation in vitro (IVM). Before mature culture, we collected the immature oocyte of germinal vesicle (GV stage).Then we washed the immature oocytes in TCM199 (Gibco, Grand Island, NY) complemented with 0.3% bovine serum albumin, and transported them to the maturation medium (TCM199 complemented with 5% estrous bovine serum and 10 g/ml follicle-stimulating hormone (FSH). The oocytes were cultured in an incubator at 38.5°C under 5% CO 2 for 22 h-24 h. Matured oocytes (MII stage) were collected under a microscope based on a morphological evaluation; completely matured oocytes included with homogenous cytoplasm that had at least three cumulus layers and were accompanied with the first polar body. Parthenogenetic embryos were collected by chemical activation, as previously described (54). Briefly, maturated oocytes were induced by exposure to 5 mM ionomycin in embryo culture medium for 5 min and subsequent incubation in 2 mM 6-dimethylaminopurine for 3 h at 39°C and 5% CO 2 . The oocytes were activated to new embryo culture medium, and two-cell stage parthenogenetic embryos were collected after 22 h-24 h, and after 2-4 days to collect embryos at the 2-to 16-cell stages, morula and blastocysts were collected at days 5 and 7. All oocytes and embryos were examined and staged under light microscopy.
RNA Sequencing and Identification of DEGs-The oocytes and embryos from 8 stages of buffalo were conducted the experimental process of transcriptional sequencing includes RNA extraction, sample detection, library construction, quality control and sequencing. We isolated total RNA from individual oocytes/embryos using Array-Pure™ Nano-scale.
RNA Purification Kit (Epicenter). We conducted Smart-Seq (55), and firstly we lysed each cell in hypotonic solution and converted poly(A)ϩ RNA to full-length cDNA using oligo (dT) priming and SMART template switching technology, followed by 12-18 cycles of PCR preamplification of cDNA, and the amplified cDNA was used to construct standard Illumina sequencing libraries. After that, Using Qubit2.0 for preliminary quantification and Agilent 2100 to test the insert size of the library, then insert size is expected and proceed with the next experiment. The qPCR method is used to accurately quantify the effective concentration of the library (library effective concentration Ͼ2 nm) and complete the checking. After the library is checked, and the Illumina Hiseq X ten platform is used for sequencing. The adaptor sequences and low-quality sequence reads were removed from the data sets. Raw sequences were transformed into clean reads after data processing. These clean reads were then mapped to the reference genome sequence. (ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/ 000/471/725/GCA_000471725.1_UMD_CASPUR_WB_2.0/). Only reads with a perfect match or one mismatch were further analyzed and annotated based on the reference genome. Tophat2 tools (56) soft were used to map with the reference genome. Cuffquant and Cuffnorm, which are the components of the Cufflinks software (57), were used to perform gene expression analysis. During the construction of the database, it was performed by mixing RNA from three oocytes or embryos of the same stage, to reduce the influence of cellular heterogeneity and differences between individual cell. Therefore, the sample at each stage resulted in a data set, so EBseq (58) was used to make a differently expressed analysis of genes.
Protein Extraction-Pooled groups of 1,000 buffalo maturated oocytes (MII oocytes) and 1,000 two-cell stage parthenogenetic embryos (PC embryos) were dissolved using 30 l of lysate (a mixture of 8 M urea, 2 mM ethylenediaminetetraacetic acid, 10m M dithiothreitol, and 1% protease inhibitor) for 30 min on ice. The samples were centrifuged (12,000 ϫ g at 4°C for 10 min), and the supernatant was transferred to a new 1.5 ml tube for further experiment.
Digestion and Peptide Labeling Using a Tandem Mass Tag (TMT) Reagent-The extracted protein was dissolved in 45 l of 100 mM triethyl ammonium bicarbonate, 5 l of 2% sodium dodecyl sulfate (SDS), and 55 l of ultrapure water. Next, we used 5 l of 200 mM trichloroethyl phosphate for 1 h at 55°C to reduce the sample, followed by 5 l of 375 mM iodoacetamide for 30 min at room temperature in an environment protected by light. Each sample was then precipitated in cold acetone (1:6) at Ϫ20°C overnight, centrifuged (12,000 ϫ g; 4°C; 10 min), and the supernatant was transferred to a new 1.5 ml tube and dissolved in 100 l of 100 mM triethyl ammonium bicarbonate. Finally, 2.5 l of Trypsin (1 g/l) was added to each sample for digestion at 37°C overnight. Peptides were labeled with TMT reagents according to the manufacturer's instructions. Briefly, 41 l of anhydrous acetonitrile (ACN) was added to the TMT reagent in each tube (0.8 mg) and dissolved for 5 min at room temperature. Next, 20 l of each tag were added to each sample for labeling and the samples were incubated for 1 h at room temperature. We then added 8 l of 5% hydroxylamine into each sample and incubated for 15 min at room temperature to terminate the reaction. Samples were pooled and evaporated in a vacuum; samples of MII oocytes and PC embryos were labeled with TMT-128 and TMT-129, respectively.
Peptide Prefractionation Using High-pH Reversed-phase High-performance Liquid Chromatography-Each sample was evaporated in a vacuum and re-suspended in 50 l of buffer A (98% ddH 2 O, 2% ACN; pH 10.0) and separated by high-pH reversed-phase liquid chroma-tography (RP-HPLC) column (2.1 ϫ 100 mm, 3 m, 150 Å, C18). For separation, we used buffer A and buffer B (98% ACN, 2% ddH 2 0; pH 10.0) for a 60-min linear gradient (4 -20% buffer B for 30 min, 20 -95% buffer B for 25 min, and 95% buffer B for 5 min. The fractions that were eluted from the column every 1.5 min were collected, and 20 fractions were collected according to the peak intensity. The collected fractions were desalted using a ZipTips C18 column (Millipore, Billerica, MA). The processed samples were stored at Ϫ80°C until being used for mass spectrometry analyses.
Peptide Identification by Liquid Chromatography-Tandem Mass Spectrometry Analysis-After being desalted, the fractions were evaporated with a vacuum, and 10 l of solvent A (2% ACN and 0.1% formic acid) were used for resuspend the samples. A volume of 2 l from each sample was loaded onto a trap column (PepMap RSLC C18 column, 50 um ϫ 15 cm, 2 m nanoViper, Thermo Fisher Scientific, Bremen, Germany), a maximum pressure (600 bar) was applied, and an analytical column (Acclaim® PepMap100C18 column, 0.075 ϫ 150 mm, 3 m, 100Å, Thermo Fisher Scientific, Bremen, Germany) was used for elution at a rate of 300 nL/min. Next, we used buffer A (2% ACN and 0.1% formic acid) and buffer B (98% ACN and 0.1% formic acid), followed by a 60-min gradient (5-40% buffer B for 45 min, 40 -100% buffer B for 10 min, 100% buffer B for 5 min) to separate the peptides. An LTQ-Orbitrap Elite hybrid mass spectrometer (Thermo Fisher Scientific, Bremen, Germany) connected to an Easy-nLC 1000 nano-liquid chromatography system (Thermo Fisher Scientific, Odense, Denmark) was used for analyses of all peptide fractions online. A data-dependent mode in the scan range of 350 -1800 m/z was carried out for the mass spectrometry (MS) analyses, and the survey scans were captured at a mass resolution of 60,000 at 400 m/z by the Orbitrap analyzer. In the linear ion trap, ten of the most intense precursor ions were selected for secondary mass spectrometry analysis (MS2) under the mode of high-energy collision dissociation. The dynamic exclusion parameter included an exclusion count of 2 and an exclusion time of 40 s. Siloxane ions were used for internal calibration (m/z ϭ 445.1200).
Raw Data Processing, Protein Quantification, and Bioinformatic Analyses-The raw data files were processed and quantified by Proteome Discoverer software v1.3.0.339 (Thermo Fisher Scientific, Massachusetts), and searched against the UniProt Bos Taurus protein database (UP000009136, 24078 sequences, release 2017_03) using the SEQUEST algorithm. The search parameters used were as follows: fixed modifications, including carbamidomethylation of cysteine (ϩ57.02146 Da); TMT reagent adducts (ϩ229.162932 Da) on the lysine and peptide N termini; and variable modification, including oxidation of methionine (ϩ15.99492 Da). Mass tolerance of precursor ions was 20 ppm, and fragment ions were set to Ϯ 0.5 Da, Enzyme specificity, trypsin. The false discovery rate (FDR) was calculated using the peptide validator and based on a decoy database search. A p value Ͻ 0.05 and a fold change Ͼ1.2 were significant.
The identified proteins were analyzed by gene ontology annotation, including by biological processes, cellular components, and molecular functions, which were based on the UniProt-GOA database (http:// www.ebi.ac.uk/GOA/). The pathway analyses of differently expressed proteins were based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. The online KEGG Automatic Annotation Server (KAAS) service tools were used to annotate each protein's KEGG database description, and then to map the annotation result using other KEGG online service tools and the KEGG mapper. The interaction between differential expression protein (DEPs) groups was derived from the Search Tools for the Retrieval of Interacting Genes/ Proteins (STRING) database. A high confidence (0.7) for the required interaction score and the active interaction sources including text mining, experiments, and databases were chosen to draw the protein-protein interaction map using Cytoscape 3.2.1.

Real-time Quantitative Reverse Transcription Polymerase Chain
Reaction-Total RNA was extracted from oocytes or embryos at different stages (MII and PC) using an RNeasy® Micro Kit (Qiagen Co.). A Takara RNA polymerase chain reaction (PCR) kit was used to reverse-transcribe the cDNAs. Quantitative reverse transcription PCR (qRT-PCR) analyses were performed on a LightCycler 96 system (Roche, Switzerland). The reaction solution (20 l) consisted of 1.5 l of cDNA, 10 l of SYBR® Premix Ex TaqTM II (Takara), 2 l of primer mix (2 nM), and 7 l of ddH 2 O. The 2 Ϫ⌬⌬CT method was used to calculate the relative expression levels of targeted genes. Each gene was analyzed using different sets of oocytes, and experiments were performed in three replicates. The primer sequences used for the qRT-PCR analyses are presented in supplemental Table S8.
Immunocytochemical Analysis-The oocytes were removed from the granulosa cells and placed in 4% polyoxymethylene for 30 min at room temperature, the oocytes were transferred to a 1% Triton X-100 phosphate-buffered saline (PBS) solution at room temperature for 20min, and the oocytes were blocked with a PBS solution containing 1% BSA (albumin from bovine serum). The oocytes were mixed with a primary antibody at an effective concentration and incubated overnight at 4°C. The oocytes were washed three times for 2 min each time with a PBS solution containing 0.1% Tween-20 and 0.01% Triton X-100. The oocytes were transferred into a diluted solution containing a secondary antibody labeled with fluorescein isothiocyanate (FITC) for 1 h at room temperature. The oocytes were then washed three times for 2 min each time with a PBS solution containing 0.1% Tween-20 and 0.01% Triton, covered with a coverslip coated in a ProLong Gold antifade reagent with DAPI (Life Technologies), and kept in the dark until an inverted fluorescence microscope was used to assess fluorescence (Olympus, IX73, Japan). As a negative control, the primary antibody was replaced with normal rabbit IgG (Immunoglobulin G) or mouse lgG.
Experimental Design and Statistical Rationale-For mRNA expression profiling, the transcriptomic analysis was conducted for eight stages for the oocytes and early parthenogenetic activated embryos of buffalo in vitro(GV, MII, 2cells, 4cells, 8cells, 16cells, morula, blastocyst) using paired-end Sequencing in Illumina Hiseq X ten platform. Quantification of gene expression levels was estimated by fragments per kilobase of transcript per million fragments mapped (59). The formula is shown as follow: FPKM ϭ cDNA Fragments Mapped Fragments (Millions) ϫ Transcript Length (kb) Before differential gene expression analysis, for each sequenced library, the read counts were adjusted by edgeR program package through one scaling normalized factor. Differential expression analysis of two samples was performed using the EBSeq R package. The resulting FDR (false discovery rate) were adjusted using the PPDE-(posterior probability of being DE) .The FDR Ͻ 0.05 and log2(fold change) Ն1 was set as the threshold for significantly differential expression. A matrix of Pearson correlation coefficient was created using R, which was in turn used to create the heatmap. Principle component analysis (PCA) was analyzed by using R.A matrix of Pearson correlation coefficient was created using R, which was in turn used to create the heatmap. Principle component analysis (PCA) was also analyzed by using R. Genes were deemed differentially expressed between hierarchical subsequent developmental stages if they showed a p value of less than 0.05 (Negative Binomial Distribution). For protein expression profiling, a TMT-based quantitative proteomic approach was employed to screen DEPs in a total of 2000 MII stage oocytes and 2-cell stage embryos. Two biological replicates are included in the proteomics experiments. The raw data files were processed and quantified by Proteome Discoverer software v1.3.0.339 (Thermo Fisher Scientific, MA), and searched against the UniProt Bos taurus protein database (UP000009136, 24078 sequences, release 2017_03) using the SEQUEST algorithm. The false discovery rate (FDR) was calculated using the peptide validator and based on a decoy database search. A p value Ͻ 0.05 and a fold change Ͼ1.2 were significant. All quantitative data are presented as the mean Ϯ S.E. from three independent experiments. A p value Ͻ0.05 was statistically significant. Statistical analyses were done using IBM SPSS Statistics version 17.0 and GraphPad Prism version 5.0. One-way Analysis of Variance (ANOVA) combined with Tukey's multiple comparison Test was used to calculate the significance of differences between different groups.

Quantitative Proteomic Analyses of Early Parthenogenetic
Embryo Stages Using Tandem Mass Tags-A total of 2000 oocytes or embryos from MII oocytes (matured oocyte) and the two-cell stage of the early embryo after parthenogenesis were treated with 5 mM ionomycin and 2 mM 6-dimethylaminopurine. Next, TMT-6-plex reagents were combined with liquid chromatography-tandem mass spectrometry (LC-MS/ MS) analyses. Proteins isolated from the MII stage and the two-cell stage were labeled with mass tags TMT-128 and TMT-129, respectively. Finally, a total of 1908 maternally expressed proteins were identified (FDR Ͻ 1%) in two biological replicates (supplemental Table S1 and S2). We used the Pearson correlation coefficient and the coefficient of variation to analyze the correlation of protein expression levels for each sample, which demonstrated that the experiment had high reproducibility ( Fig.1A and 1B). A total of 1224 maternally expressed proteins were found in both experiments. A complete list of all identified proteins with annotation from every repeat are shown in supplemental Table S3. Compared with the data from the MII stage, 123 proteins were expressed differently (fold change Ͼ 1.2; p Ͻ 0.05); of the123 proteins, 67 proteins were up regulated and 56 proteins were downregulated (Fig. 1C). We also analyzed the subcellular structure distribution for these proteins and found that they were mainly located in the following areas of the cell: cytosol (34.15%), nucleus (18.70%), mitochondria (14.63%), extracellular matrix (14.63%), and plasma membrane (9.76%) (Fig. 1D). These findings demonstrate that the differentially expressed proteins are widely distributed in the cell. The top 10 differently maternally expressed proteins are listed in Table I, and they are sorted by fold change. A complete list of all differently expressed proteins with annotation are shown in supplemental Table S4.
We also used the Cytoscape 3.2.1 mapping software to visualize the interaction between the 1,908 identified proteins ( Fig. 1E and supplemental Table S5). It can be seen from the figure that the relationship between these proteins is complicated. Using the BinGO package for gene ontology (GO) analysis, we found that the proteins were mainly involved in the regulation of biological processes, metabolism, transport, cellular processes, and signaling processes.
Gene Ontology (GO) Analysis of Differentially Expressed Proteins-The differentially expressed proteins (upregulated and downregulated) were classified by gene ontology annotation based on three categories: biological process, cellular component, and molecular function. The GO annotation proteome was derived from the UniProt-GOA database (http:// www.ebi.ac.uk/GOA/) (supplemental Table S6). After GO analysis of up-regulated proteins and an analysis of the biological processes, we found that the biological processes include plasma membrane organization, actin cytoskeleton organization, processes associated with the cell cycle, positive regulation of cell junction assembly, processes associated with the cell junctions, oocyte differentiation, and regulation of cell communication, signal transduction, and cell proliferation ( Fig. 2A). The molecular function of these proteins is important in actin-dependent ATPase activity (Fig. 2B). The results of the cellular component analysis showed that these proteins were mainly as the component of caveola, lamellipodium and actin cytoskeleton (Fig. 2C).
A GO analysis was conducted on the downregulated proteins, and we found that the results of this group (41 GO terms) was much less than the upregulated proteins (273 GO terms), which may indicate that upregulated proteins are more important than the downregulated proteins during parthenogenetic development in early embryos. In an analysis of biological processes, we identified some processes such as intermediate filament cytoskeleton organization and regulation of hormone secretion (Fig. 2a). The molecular function of these proteins plays a major role in the binding of vitamins and carbohydrates, and in the activity of oxidoreductase, transaminase, and methyltransferase (Fig. 2b). Furthermore, the downregulated proteins are mainly a component of keratin filament, intermediate filament cytoskeleton, and the extracellular matrix (Fig. 2c).
KEGG Pathway Analysis of Differentially Expressed Proteins-We analyzed differentially expressed proteins based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (which is shown in Fig. 3A-3B and supplemental Table S7). We found that the upregulated proteins were mainly enriched in tight junctions, phosphatidylinositol 3-kinase (P13k-Akt) signaling pathways, the regulation of actin cytoskeleton, Rap1 signaling pathways, gap junctions, and cellular adhesion molecules (CAMs); these signaling pathways may be associated with early embryonic development. Of note, we found that the number of proteins enriched in the tight junctions is the largest; the gap junctions are also significantly enriched, which is like the tight junction pathway. These revealed that these proteins might play an important role in the formation of intercellular junctions during early embryonic development. The connective complex of cells plays an important role in the polarization of cells and the formation of the blastocele, and these connections include the tight junction at the top of the cell, the gap junction on the side, and the desmosomes. We then picked out the identified proteins that were involved in the two signaling pathways and used the information to create a visualization of the interactions ( Fig. 3C and 3D) and to conduct KEGG mapping ( Fig. 3E and 3F). We found nine proteins (MYH8, GNAI1, MYH4, NRAS, CTNNA1, MYL12B, CTTN, MYH10, KRAS) that were upregulated in tight junctions ( Fig. 3C and 3E) and five proteins (GNAI1, NRAS, TUBB3, MAP2K1, KRAS) that were upregulated in gap junctions, and we also found that MAP2K5 was downregulated in gap junctions ( Fig. 3D and 3F). The time of compaction of the blastocyst was at the 8-to 16-cell stage during embryonic development in bovine (60 -63). However, in this study we found that many proteins were involved in the formation of tight junctions and gap junctions at the two-cell stage of early parthenogenetic embryo development in the buffalo. We also found that the development rate of parthenogenetic oocyte cleavage was advanced by at least 6 h compared with the embryo that underwent in vitro fertilization in the actual experiment, which indicates that the time of compaction during early development of parthenogenetic embryo was earlier than the fertilized eggs. Maternal proteins may play an important role in the compaction of preimplantation embryos, and these proteins may participate in the formation of cellular junctions shortly after parthenogenetic activation. To verify our ideas, we selected four upregulated proteins in the tight junction and gap junction for subsequent validation analysis (Table II).
Transcriptomic Analysis-Using linearly amplified RNA from oocytes or embryos, we obtained more than 570 million sequencing reads from 8 stages of buffalo immature and matured oocytes and pre-implantation embryos (supplemental Table S9). The total numbers of detectable genes ranged from immature oocytes to blastocysts of buffalo in our research with the use of RNA-seq are shown in supplemental Table  S10.This experiment identified 9112 differential expressed genes(p Ͻ 0.05) (supplemental Table S11). A total of 3567 unique genes were identified to be differently expressed between all consecutive stages of development (FPKMϾ0) (supplemental Table S12).

Correlation and Hierarchical Analyses of Transcriptional Profiles and the Expression of NPM2 and NLRP5 mRNA of Buffalo In Vitro Developed Oocytes and Parthenogenetic
Embryos-In oocytes/embryos RNA sequencing research, Pearson correlation coefficients ( Fig. 4A and supplemental Table  S10) and principal component analyses (PCA) (Fig. 4C) of all detected genes revealed the two distinct segmentations of buffalo pre-implantation development: the first from GV (immatured oocyte) to 8-cell,and second from 16-cell to blastocyst stage (Fig. 4A and 4C), which is consisted with the result of Hierarchical clustering of differentially expressed genes in developed oocytes and parthenogenetic embryos of buffalo (Fig. 4B). This segmentation demonstrated that parthenogenetic activation of embryonic development exists an embryonic genome activation process which is the same as the normal fertilization of embryos. Also, the EGA in buffalo occurs between 8-cells-16-cells, consistent with the time of the development of the normal fertilized embryonic development (60 -63).
In the quantitative proteomic study, we found two maternal proteins-NPM2 and NLRP5 (official name for MATER)-that have been previously reported (Table III). NPM2 and NLRP5 FIG. 1. A, The correlations between the protein expression levels of each sample were calculated by using the Pearson correlation coefficient. The red represents the MII stage and the blue represents the two-cell stage after parthenogenetic activation. B, The coefficient of variation (CV) is used to reflect the degree of discretization of each sample data, which is the ratio of the standard deviation of the data to the mean, expressed as a percentage (CV ϭ [standard deviation/average value] ϫ 100%). The horizontal axis is the CV value, and the vertical axis represents the proportion of total protein. The red line represents the MII stage and the green line represents the two-cell (that is, the parthenogenetic cleavage [PC]) stage. C, Each point in the volcanic map represents a protein; the abscissa represents the logarithmic value of the expression of a certain protein in the two samples. The green point represents the downregulation of the expressed protein, the red point represents the upregulated differential expression protein, and the black dot represents the non-differentially expressed protein. D, An analysis of the distribution of proteins in the subcellular structure found that the were mainly located in the cytosol (34.15%), nucleus (18.70%), mitochondria (14.63%), extracellular matrix (14.63%), and plasma membrane (9.76%). E, The interaction between the identified proteins visualized using Cytoscape 3.2.1 software. The color of the node is represented by the betweenness, which means the influence in the network; if the color is redder, the influence is greater. The size of the node is represented by the degree, which means the stability in the network; if the shape is bigger, the stability is stronger. Using the BinGO package, we found that the identified proteins were mainly involved in the following processes: regulation of biological processes, metabolism, transport, cellular process, and signaling process.
are maternal genes that are both involved in the regulation of embryonic gene activation (9,21). Notably, expression of NPM2 and NLRP5 is upregulated at the two-cell stage after parthenogenetic activation. We used the RT-qPCR method to detect the expression of the two genes during the early stages of buffalo embryonic development, and GV (germinal vesicle) oocyte was used as control. We found that there was not a significant difference between the two genes before the 16cell stage. Expression of the GV stage was highest compared with all other stages, and the expression during the two-cell stage was higher compared with during the MII stage; these findings were consistent with proteome data. After the eightcell stage, the expression of the two genes was significantly decreased, and they were little expressed during the blastocyst stage (Fig. 4D-I and II). The results of the analysis of NPM2 and NLRP5 were consistent with a previous report (64,65). In addition, we analyzed the tissue expression specificity of NPM2 and Nlrp5, and we found that they only expressed in the ovary, MII oocyte, the two-cell stage embryo, and granular cells (Fig. 4D-III and V). NLRP5 and NPM2 showed similar expression changes, indicating that the two genes may be an important maternal protein during early embryonic development in the buffalo. The analysis of quantitative proteomic study combined with research of RNA sequencing demonstrates that the "maternal-to-zygotic transition"(MZT) possibly exists in parthenogenesis similar to the fertilized embryo(maternal genome and paternal genome) and occurs between the 8-cell stage to 16-cell stage.
Proteomic Profile Validation-The four upregulated proteins involved in the formation of gap junctions and tight junctions were chosen for Western blotting, RT-qPCR, and cellular immunofluorescence localization analysis (Fig. 5A-5C). The results of Western blot analysis of the four proteins (CDH3, CTNNA1, MAP2K1, TUBB3) are presented in Fig. 5B. We analyzed the gray value of each lane of Western blotting using Image J software, and the results were consistent with the analyses of TMT quantification and transcription levels (Fig.  5C). All results used MII stage as a control. To further explore the potential functions of the differentially expressed proteins, we used cellular immunofluorescence localization analysis to investigate these proteins because there is little-known information on their cellular localization in the MII oocyte and the two-cell embryo. As shown in Fig.  5A, CTNNA1 was mainly expressed in the cell membrane, CDH3 was mainly expressed in cytoplasm and nucleus, and the other two proteins (MAP2K1 and TUBB3) were both mainly located in the cytoplasm. The fluorescence intensity of the four proteins in the two-cell embryo was higher compared with the MII oocyte, which indicates that the gap junctions and tight junctions in parthenogenetic embryos may regulate embryo densification processes earlier than fertilized egg embryos. DISCUSSION Proteomic and transcriptomic methods have been applied to the research of oocytes and embryos in mammals, and proteomic data have been obtained that include protein expression profiles of oocytes and embryos at different developmental stages. A growing list of maternal proteins encoded by maternal-effect genes have been characterized using proteomic approaches, such as CDH1 (66 -68), Dnmt1 (31,69), DPPA3 (29,70), H3.3 (71)(72)(73), NLRP5 (official name for Mater) (9,65), NPM2 (21,74), and NOBOX (75,76). Various aspects of early embryonic development are regulated by these factors, such as epigenetic modifications and chromatin remod-eling, onset of embryonic genome activation, maternal mRNA degradation, cell compaction and signal transduction protein translation. However, there are no reports to date that focus on expression profiles of maternally expressed proteins using parthenogenetic activation as an experimental model during early embryo development in the buffalo. Parthenogenetic embryos only contain maternal genetic information that can be used to facilitate the study of the role of the maternal genome in early embryonic development and to discover the maternally expressed genes and proteins that are involved in cleavage and cell differentiation. In addition, transplantation studies allow for an investigation of the role that the oocyte plays in nuclear donor replication; the start of this development is completely dependent on the activation of oocytes. Therefore, the study of the parthenogenesis process is also conducive to improving nuclear transfer technology.  General Aspects-In this study, we undertook a panoramic analysis of the proteomes of two-cell embryos, which is the first developmental stage after parthenogenetic activation, and made a comparison to MII oocytes using TMT-labeled quantitative proteomics. We generated two biological repli-cates for each stage, and we were able to identify 1908 proteins in 1908 protein clusters with high confidence (FDR Ͻ 1%). We used a Pearson correlation coefficient and a coefficient of variation for our analysis and found that the two experiments had high reproducibility (Fig. 1A and 1B). We

FIG. 4. Correlation and hierarchical analyses of transcriptional profiles and the expression of NPM2 and NLRP5 mRNA of buffalo in vitro developed oocytes and parthenogenetic embryos.
A, Heatmap of all stages from buffalo oocytes and embryonic development of parthenogenesis. The color spectrum, ranging from red through white to blue, represents Pearson correlation coefficients ranging from 1 to 0.2, indicating high to low correlations, and the correlation between the eight stages of expression can be divided into two stages: the first stage from GV to 8cells which is high correlation between this five stages, and the second stage of higher correlation is from 16cells to blastocyst. B, Hierarchical clustering of differentially expressed genes in developed oocytes and parthenogenetic embryos of buffalo. Two major clusters are shown, one consisted of the GV (immatured oocytes), MII (matured oocytes) and embryos at the 2-cell, 4-cell and 8cell stages. The second cluster is consisted of embryos at the 16cell, morula and blastocyst stages. The separation of embryos into two groups demonstrated that parthenogenetic activation of embryonic development exist an embryonic genome activation process which is the same as the normal fertilization of embryos, also the EGA in buffalo occurs between 8 cells-16 cells. C, Principal component analysis (PCA) of the transcriptomes for eight stages of buffalo in vitro developed oocytes and parthenogenetic embryos. Embryos from the same stage are shown by symbols of the same shape. PCA1, PCA2 and PCA3 represent the top three dimensions of the differentially expressed genes among the parthenogenetic embryos. D-I and II, Real-time quantitative polymerase chain reaction analysis of buffalo NPM2 and NLRP5.The mRNA expression during in vitro oocyte maturation and early embryogenesis after parthenogenetic activation (mean Ϯ standard deviation; **, p Ͻ 0.01; ***, p Ͻ 0.001). (D-III and V). The tissue distribution of NPM2 and NLRP5 in buffalo.GAPDH (glyceraldehyde-3-phosphate dehydrogenase) was used as a positive control for RNA quality. then used Cytoscape 3.2.1 software to visualize the interactions between all of the identified proteins and, using the BinGO package, we found that the biological processes of those proteins were mainly involved in the regulation of metabolism; transport; cellular processes; signaling processes ( Fig. 1E), which revealed that embryo DNA replication is rapid in the cleavage stage; and the mechanism of cleavage, which involves changes to mitochondria and surface cytoplasm, including the karyokinesis and cytokinesis. Many proteins, RNAs, and other substances that are required in the early stages of embryonic development are mainly derived from oocytes; maternal RNAs and proteins that control early embryonic development begin to synthesize and accumulate during the oocyte development, growth, and maturation in mammals. We analyzed the distribution of subcellular structures for all identified proteins and found that they were mainly located in the cytosol (34.15%), nucleus (18.70%), mitochondria (14.63%), extracellular matrix (14.63%), and plasma membrane (9.76%) (Fig. 1D). The distribution of the structures indicates that metabolism, signal transduction, and transcriptional translation are active after the activation of oocytes, which is consistent with the results of the BinGO package analysis. We successfully identified 123 differentially expressed proteins (fold change Ͼ 1.2; p Ͻ 0.05) at the two-cell stage compared with the data from the MII stage-67 of the proteins were upregulated and 56 of the proteins were downregulated (Fig. 1C). The 123 proteins were classified by gene ontology annotation, and we found that the results of the downregulated proteins (42 GO terms) were much less than the upregulated proteins (273 GO terms). This finding may indicate that upregulated proteins are more active than the downregulated proteins during the first developmental stage following parthenogenetic activation, and that two maternal proteins that have been previously reported (that is, NPM2 and NLRP5) have both identified in the present study as being upregulated (Table III). After an analysis using the KEGG database to investigate the differentially expressed proteins, we found that the upregulated proteins were mainly enriched in the tight junctions, P13k-Akt signaling pathway, regulation of actin cytoskeleton, Rap1 signaling pathway, gap junctions, and CAMs (Fig. 3A), and the downregulated proteins are related to the metabolism of arginine, proline, glyoxylate, and dicarboxylate (Fig. 3B). Notably, we found that the largest number of enriched proteins was in the tight junctions, followed by the gap junctions, which led us to infer that these proteins may play an important role in the formation of intercellular junctions during early embryonic development. We chose four upregulated proteins (CDH3, CTNNA1, MAP2K1, TUBB3) involved in tight junctions and gap junctions for validation analyses using RT-qPCR (Fig. 5C), Western blotting (Fig. 5B), and cell immunofluorescence localization (Fig. 5A), and we found that these results were consistent.
Also, we conducted the RNA sequencing-based transcriptomic analysis of the oocytes and embryos from eight stages (GV, MII, 2-cell, 4-cell, 8-cell, 16-cell, morula, blastocyst) of Buffalo, and we obtained more than 570 million sequencing reads from these stages; for the first time, it is revealed that the buffalo oocytes and early embryos expressed more than 60% of the total estimated 29,424 genes in the buffalo genome, which is a dynamic expression profile of maternal genome during pre-implantation development in Buffalo. We conducted correlation and hierarchical analysis and PCA (NPM2 and Nlrp5) between 8 and 16 cells, which found in proteomic data. In addition, to confirm the results from the transcriptomic data, we conducted quantitative real-time PCR (qRT-PCR) on 10 genes (SNAI1, LFNG, POLRMT, RBP4, RPRD1B, FOSL1, F2, CDK5R1, WWTR1, and SETD2) that upregulated in 16-cell stage using parthenogenetic embryos at the 8-cell, 16-cell, and blastocyst (n ϭ 2) stages. We found that all these genes are upregulated in 16-cell embryos, whose change trends are consistent with those from the RNA-seq (supplemental Fig. S2).
Maternally Expressed Proteins Are Possibly Involved in the Formation of Cellular Junctions First After Parthenogenetic Activation-In the tight junctions, there were nine proteins (MYH8, GNAI1, MYH4, NRAS, CTNNA1, MYL12B, CTTN, MYH10, KRAS) that were upregulated ( Fig. 3C and 3E). In the gap junctions, there were five proteins (GNAI1, NRAS, TUBB3, MAP2K1, KRAS) that were upregulated and one protein (MAP2K5) that was downregulated ( Fig. 3D and 3F). Four proteins (CDH3, CTNNA1, MAP2K1, TUBB3) were validated by a combined analysis (Fig. 5C); CTNNA1 was mainly expressed in the cell membrane, CDH3 was mainly expressed in the cytoplasm and nucleus, and other two proteins (MAP2K1, TUBB3) were both located mainly in the cytoplasm (Fig. 5A). Those findings indicate that the four proteins may be involved in the biological processes of metabolism, transport, cellular processes, and signaling processes (Fig. 1E). The cellular junctions exist in the same tissue or different organizations and are formed by the adjacent cell surface and cell gap, and then organizational structure can maintain a relatively stable structure between cells.
Tight junction biogenesis regulates pattern formation during embryogenesis (77). The tight junctions begin to form after compaction occurs at the late eight-cell stage in the mouse embryo. In the present study, CDH3 and CTNNA1 are involved in the formation of tight junctions and they are upregulated during the two-cell stage (Fig. 5B-5C). CDH3 encodes a classic cadherin from the cadherin superfamily. Alternative splicing results in multiple transcript variants, at least one of which encodes a preproprotein that is proteolytically processed to generate the mature glycoprotein. CDH3 is a calcium-dependent cell-cell adhesion protein that is comprised of five extracellular cadherin repeats that make up a transmembrane region and a highly conserved cytoplasmic tail. GO annotations related to this protein include calcium ion binding, which is involved in cell junction organization and adhesion. In the present research, the result of cell immunofluorescence localization analysis revealed that CDH3 is in the cytoplasm and nucleus (Fig. 5A), which is consistent with its function. The process of tight junction membrane assembly is dependent on prior activation of the E-cadherin adhesion system. When cultured with neutralizing anti-E-cadherin antibody, embryo compaction is inhibited and membrane assembly of tight junctions is seriously disrupted and cannot target the apicolateral site of cell contact (78). Similarly, the normal assembly of tight junctions cannot develop in E-cadherin null mutant embryos (79). In addition, there is evidence that Ecadherin adhesion at compaction acts to balance tight junction proteins from turnover in case they are assembled at the membrane by their anchorage to the cytoskeleton (80). CTNNA1 encodes a member of the catenin protein family that is critical to cell adhesion, and it is located on the plasma membrane with conjoining cadherins; this finding is consistent with immunofluorescence localization (Fig. 5A). The association between catenins and cadherins forms a solid complex that is linked to the actin filament network. It seems to be critical for cadherin cell-adhesion properties, and it is associated with both E-and N-cadherins. The solid complex of E-cadherin-catenin adhesion mediates the connection between cadherins and the actin cytoskeleton at adherence junctions. We also found many upregulated proteins that are primarily involved in Golgi localization ( Fig. 2A), which is like the occludin protein (an important protein for the formation of tight junctions). The ZO-1␣ϩ isoform is transcribed and translated for the first time in the embryo; it associates with occludin at the Golgi complex, which occurs concurrently with maturation steps in occludin processing (81).
Gap junctions are essentially composed of only one protein, connexin, which is an integral membrane protein made up of four membrane-spanning domains (82). In the mouse embryo, the dominant proteins contributing to gap-junction formation are Cx32 and Cx43. Cx32 is maternally expressed in the mouse, whereas Cx43 expression appears in the embryonic genome and increases during later stages of cleavage (83)(84)(85). The onset of gap-junction coupling at the eight-cell stage is regulated developmentally, and this evidence suggests that DNA replication during the two-cell stage may be required for coordination of this process (86). In the present study, we found many maternally expressed proteins that are upregulated at the two-cell stage. At the time of assembly, Cx43 becomes phosphorylated, and stimulation of protein kinase A enhances both phosphorylation and assembly events (87). In the present study, we also found that MAP2K1 and TUBB3 were involved in gap junctions and were upregulated at the two-cell stage. MAP2K1 belongs to the dual specificity protein kinase family, which acts as a mitogen-activated protein (MAP) kinase. MAP kinases, also known as extracellular signal-regulated kinases (ERKs), act as an integration point for multiple biochemical signals; its related pathways include Trk receptor signaling, which is mediated by the MAPK pathway, and development of VEGF signaling via the VEGFR2-generic cascades. GO annotations for this gene include transferase activity, transferring phosphorus-containing groups and protein tyrosine kinase activity, and MAP2K1 may have the function of making proteins that are associated with gap junction phosphorylation, such as Cx43 and Cx32. TUBB3 encodes a class III member of the beta tubulin protein family. Beta tubulins are one of two core protein families (alpha and beta tubulins) that heterodimerize and assemble to form microtu-bules. Tubulin is the major constituent of microtubules, and it binds two moles of GTP, one at an exchangeable site on the beta chain and one at a non-exchangeable site on the alpha chain. TUBB3 plays a critical role in proper axon guidance and maintenance. In the present study, we found that TUBB3 was in the cytoplasm, which was consistent with its role in cytoskeleton formation. Gap-junction formation between mouse blastomeres occurs for the first time during the eight-cell stage (88,89) following synthesis of constituents during the preceding cell cycle (86,90).
The "Maternal-to-Zygotic" Process Possibly Exists in Parthenogenesis and Occurs Between the 8-cell Stage to 16-cell Stage Like Normal Embryo-In general, the first cleavages are primarily supported by maternal factors, whereas blastocyst formation requires a combination of maternal and embryonic factors (65). NPM2 protein, also known as nucleoplasmin, is one of the crucially maternal proteins in the Xenopus oocyte nucleus (91)(92)(93). NPM2 is a nuclear chaperone protein that is expressed in oocytes. It is an acidic and highly phosphorylated protein that joins histones to form a stable pentamer complex, and it participates in nucleosome assembly (94,95). Also, nucleoplasmin is found natively bound to histones H2A and H2B, facilitating the storage of these proteins in the oocytes (96). NPM2 may also be involved in the formation of nucleolus-like bodies during oocyte growth (97), and it is necessary for the composition of oocyte nuclear and nucleolar domains and the compaction of oocyte chromatin during the final stages of oocyte development (21). NPM2 knockout mice (NPM2-null) are infertile or subfertile, and embryos that lack NPM2 exhibit a reduced initial cleavage rate, obstruct the oocyte development to the two-cell stage, and lead to cell death 50 h after fertilization (21).
In the present experiment, TMT quantitative methods were used to determine that the expression of NPM2 is upregulated at the two-cell stage following parthenogenetic activation (Table III). We used an RT-qPCR method to detect the expression of the two genes(NPM2 and NLRP5) during the early stages of buffalo embryonic development. When compared with the GV stage, we found that there was no significant difference in the presence of NPM2 before the eight-cell stage. We also found that the expression of the GV stage was highest compared with all other stages, and the expression of the two-cell stage was higher than the MII stage, which was consistent with proteome data. Also, despite the existence of species specific differences in the precise timing of the onset of zygotic transcription in mammals, commonalties exist including evidence for a minor degree of transcriptional activity preceding the major activation of the genome in bovine, and patterns of gene expression may also be conserved, at least between the cow and the mouse, including a transient expression pattern which may represent a hallmark of chromatin remodeling during ZGA (98). After the development of the eight-cell stage, the expression of NPM2 sharply decreased, and it was not expressed in the blastocyst stage ( Fig. 4D-I), which may indicate that the NPM2 protein is maternally inherited in buffalo eggs. The expression pattern of buffalo NPM2 mRNA during early embryogenesis was very similar to some known maternal effect genes that are essential for early embryonic development, such as MATER (65), ZAR1 (99,100), and KPNA7 (101). We analyzed the tissue expression specificity of NPM2 and found that it only expressed in the ovary, MII oocyte, the two-cell stage embryo, and granular cells (Fig.  4D-III), which revealed that NPM2 may be reproduction-related maternal proteins in Buffalo.
MATER (a maternal antigen that embryos require; also known as NLRP5) is an oocyte-specific maternal effect gene that is necessary for early embryonic development beyond the two-cell stage in mice (65). In the present experiment, NLRP5 mRNA decreased strongly after the eight-cell stage, it progressively decreased during the embryo cleavage stages, and it was hardly detected in blastocysts ( Fig. 4D-II) , which is similar to the behavior of NPM2 and consistent with previous research (65). The mRNA expression of NPM2 and NLRP5 sharply decreased after the eight-cell stage in buffalo. We also observed a consistent decrease in NPM2 mRNA levels during early development and, more specifically, around the time of MZT, which is a period during which developmental control is transferred from maternally provided gene products to those synthesized by the zygotic genome (102). We also analyzed the tissue expression specificity of MATER and found that it only expressed in the ovary, MII oocyte, the two-cell stage embryo, and granular cells (Fig. 4D-V), which revealed that MATER might be reproduction-related maternal proteins in buffalo, like NPM2.
Also, combing with the analysis of RNA sequencing, Pearson correlation coefficients ( Fig. 4A and supplemental Table  S13) and principal component analyses (PCA) (Fig. 4C) of all detected genes revealed the two distinct segmentations of buffalo pre-implantation development: the first from GV (immatured oocyte) to 8-cell,and second from 16-cell to blastocyst stage ( Fig. 4A and 4C), which is consisted with the result of Hierarchical clustering of differentially expressed genes in developed oocytes and parthenogenetic embryos of buffalo (Fig. 4B).
All these results demonstrate that parthenogenetic embryonic development and normal fertilization embryo development may have a similar genome activation transcription regulation process and the MZT process also possibly exists during parthenogenesis, and it may occur between the 8-cell stage to 16-cell stage in buffalo.
In conclusion, this study provided the first insights into the numerous maternally expressed proteins and genes. Our research demonstrates that the maternally expressed proteins after parthenogenetic activation are possibly involved in the formation of cell junctions firstly. Also, the two maternal proteins (NPM2, NLRP5) that have been reported in previous studies were confirmed again in the present study. Also, correlation and hierarchical analyses of transcriptome from oocytes/embryos and the expression of NPM2 and NLRP5 mRNA of buffalo in vitro developed oocytes and parthenogenetic embryos indicate that parthenogenetic embryos may exist an "MZT" like normally fertilized embryos, and that they may occur between the 8-cell stage to 16-cell stage. These data provide a rich resource for further studies on maternal proteins and genes, and they improve our understanding of the mechanisms of maternal genomics regulation and are conducive to improving nuclear transfer technology, also the strategy of this research may provide a valuable reference to explore the function of maternal genome for other species.
Acknowledgments-We thank Dr. Bernard A Goodman for providing language assistance and editorial help. We also appreciate Yizhi Zhong (YanLeShang Biotech, Guangxi, China) for providing the technical support of Western blotting, RT-qPCR, and cellular immunofluorescence localization analysis.

DATA AVAILABILITY
The MS raw data for TMT-LC-MS/MS have been deposited to the Proteome Xchange Consortium via the PRIDE partner repository (103) with the dataset identifier PXD008549.The RNA-Seq data have been deposited in NCBI's Gene Expression Omnibus (104) and are accessible through GEO Series accession number GSE108578 (https://www.ncbi.nlm.nih. gov/geo/query/acc.cgi?accϭGSE108578).