Transcriptome Profiling of Rabbit Parthenogenetic Blastocysts Developed under In Vivo Conditions

Parthenogenetic embryos are one attractive alternative as a source of embryonic stem cells, although many aspects related to the biology of parthenogenetic embryos and parthenogenetically derived cell lines still need to be elucidated. The present work was conducted to investigate the gene expression profile of rabbit parthenote embryos cultured under in vivo conditions using microarray analysis. Transcriptomic profiles indicate 2541 differentially expressed genes between parthenotes and normal in vivo fertilised blastocysts, of which 76 genes were upregulated and 16 genes downregulated in in vivo cultured parthenote blastocyst, using 3 fold-changes as a cut-off. While differentially upregulated expressed genes are related to transport and protein metabolic process, downregulated expressed genes are related to DNA and RNA binding. Using microarray data, 6 imprinted genes were identified as conserved among rabbits, humans and mice: GRB10, ATP10A, ZNF215, NDN, IMPACT and SFMBT2. We also found that 26 putative genes have at least one member of that gene family imprinted in other species. These data strengthen the view that a large fraction of genes is differentially expressed between parthenogenetic and normal embryos cultured under the same conditions and offer a new approach to the identification of imprinted genes in rabbit.


Introduction
Embryonic stem cells (ESCs) have enormous potential in biomedicine for cell replacement, drug screening, predictive toxicology and developmental studies [1] and are envisaged as a powerful source of pluripotent cells for differentiation into desirable tissue for regenerative medicine and cell therapy [2,3]. Despite the tremendous potential of ESCs, their handicap is the isolation method, as they are obtained from the inner cell mass of a blastocyst, making the embryo unviable [4].
Parthenogenetic embryos are being studied as an alternative source of ESCs, which would avoid ethical concerns related to destruction of the embryo [4,5]. ESCs derived from parthenogenetic embryos (pESCs) have been shown to differentiate into all cell types and functional organs in the body [6]. However, several studies have evaluated similarities and differences between parthenogenetic and conventional ESCs in pluripotency, karyotype, in vivo and in vitro differentiation ability and RNA expression levels in human, nonhuman primates and rabbit [1,2,3,5,7,8]. Generally, they present normal karyotypes and are similar in their undifferentiated state, expressing normal pluripotency markers, but present different transcriptomes, with different expression patterns of extracellular matrix proteins and methylation.
In rabbit, ESCs lines from different origin have been derived and characterised [8,9]. Fang et al. [8] showed that ESCs derived from fertilised, parthenogenetic and nuclear transfer embryos seem to be similar, in that all three types were able to give rise to cells and tissue types of the three primary germ layers when ESCs are cultured in vivo and in vitro. In this case, ESCs of parthenogenetic and nuclear transfer embryos were derived using the same protocol. However, the origin of the source of the cell line has important consequences [1]. Piedrahita et al. [10] showed that ESCs lines from mice and pigs derived with the same protocol have some similar characteristics, but not all. Under in vitro culture, parthenote embryos present altered mRNA expression patterns, while in vivo developed parthenotes seem to be similar to normal embryos for the expression of factor OCT-4, Vascular Endothelial Growth Factor, Epidermal Growth Factor Receptor 3 and Transforming Growth Factor b2 genes [11]. In fact, in parthenote embryos the maximum development reached in all mammals species has been reported when embryos were transferred to subrogate females in early stages of development, providing a large in vivo culture.
In the present work, we employed a microarray to characterise transcriptome differences between 6-day parthenote embryos and 6-day fertilised blastocysts developed in vivo. In addition, based on the list of candidate genes identified by microarray, we studied the expression levels of selected transcripts in the parthenotes and fertilised blastocyst derived in vivo and checked this list with a database of genes previously listed as imprinted, while also reporting the identification of putative imprinted genes in rabbit blastocysts.

Materials and Methods
All chemicals in this study were purchased from Sigma-Aldrich Química S.A. (Madrid, Spain) unless stated otherwise.

Parthenogenetic oocyte activation
To obtain oocytes for parthenogenetic activation, 32 receptive does were induced to ovulate with an intramuscular dose of 1 mg of Buserelin acetate. Does were slaughtered 16-18 h postinduction of ovulation and the reproductive tract was immediately removed. Oocytes were recovered by perfusion of each oviduct with 5 mL of pre-warmed Phosphate Buffered Saline without calcium chloride (PBS) and supplemented with 0.1% of Bovine Serum Albumin (BSA). Recovered oocytes were submitted to two sets 1 h apart of two DC electrical pulses of 3.2 kv/cm for 20 ms at 1 sec apart in an activation medium (0.3 M mannitol supplemented with 100 mM MgSO 4 and 100 mM CaCl 2 ), followed by 1 h exposure in TCM199 medium supplemented with 5 mg/mL of cycloheximide and 2 mM of 6-DMAP. A total of 369 oocytes were activated.

Oviductal transfer by laparoscopy
Presumptive parthenotes were transferred by laparoscopy into oviducts of 13 synchronised receptive does just after activation, whose ovulation was induced as previously described [12,13]. About 28 activated oocytes per doe were transferred. Receptive does were anaesthetised by an intramuscular injection of 16 mg xylazine (Rompun; Bayern AG, Leverkusen, Germany), followed by an intravenous injection of ketamine hydrochloride at the rate of 25 mg/kg body weight (Imalgene 1000; Merial S.A, Lyon, France) to keep does under anaesthesia during laparoscopy. Females were slaughtered 6 days later and parthenote blastocysts were recovered by uterine horns perfusion with 20 mL of Dulbecco Phosphate Buffered Saline (DPBS) supplemented with 0.1% of BSA.
Control embryo recovery at day 6 of development Six receptive does were artificially inseminated with pooled sperm from fertile males [14] and induced to ovulate as previously described. In vivo fertilised embryos were collected from does slaughtered at 6 days of pregnancy by flushing uterine horns as previously described.

RNA extraction, amplification and sample labelling
As the amount of RNA present in a single embryo is rather limited [15], for each experimental group (parthenotes and in vivo fertilised embryos) four independent pools consisting of seven embryos were produced. Total RNA was isolated using traditional phenol/chloroform extraction by sonication in the Trizol reagent (Invitrogen). Concentration, quality and integrity of RNA were evaluated by Bioanalyzer 2100 (Agilent Technologies). Afterwards, 150 ng of Total RNA were amplified and labelled using QuickAmp Labelling Kit (Agilent Technologies, Madrid, Spain), following the manufacturer's instructions, which employs a linear amplification method with T7 polimerase. Control embryo samples were labelled with Cyanine 5 dye (Cy5) and parthenote

Hybridisation, washing and scanning of Microarrays
Equal amounts of Cy3 and Cy5 labelled samples (825 ng) were mixed with 106 Blocking Agent and Fragmentation Buffer, and then 55 mL of the mixture were hybridised into the commercial microarray specific for rabbit (Rabbit 446 oligonucleotide array; cat: G2519F -020908, Agilent Technologies, Madrid, Spain). This microarray was manufactured using the Agilent 60-mer SurePrint technology, which represented sequences of Refseq, Unigene and Ensembl databases (specifically 12083 identifiers of genes corresponding to the ENSEMBL database). After 17 hours at 65uC, hybridised slides were washed and scanned using the Agilent DNA Microarray Scanner G2565B (Agilent Technologies, Madrid, Spain). The resulting images were processed using the Feature Extraction v.10 Software (Agilent Technologies, Madrid, Spain) with default parameters. Only microarrays which passed control quality tests of Feature Extraction Software were used in posterior analysis.

Microarray data analysis
Filtering of problematic probes identified as flag outliers and identification of differentially expressed genes between both experimental groups were performed using the software Gene-Spring v.11.5 (Agilent Technologies, Madrid, Spain). A nonsupervised analysis of global gene expression was performed using the principal components analysis (PCA). To identify differentially expressed genes, we used the T-test with Benjamini and Hochberg multiple test correction implemented in the GeneSpring (Agilent Technologies). Probe sets were considered differentially expressed between two conditions if they had a false discovery rate (FDR) of p-value,0.05. Gene Ontology analysis and functional annotation of differentially expressed genes were performed by Blast2GO software v.2.5.1 with default parameters [16]. All data sets related to this study were deposited in NCBI's Gene Expression Omnibus [17] and are accessible through GEO Series accession number GSE41043.

Real-time qPCR
To validate the microarray results obtained, six genes (IMPACT; SMARCA2: SWI/SNF related matrix associated actin dependent regulator of chromatin subfamily A member 2; EMP1: Epithelial membrane protein 1; DPY30; CALC: calcitonin gene-related peptide variant 1; SCGB1A1: secretoglobin family 1A member 1) that showed a significant difference between experimental groups were selected and analysed in twelve independent pool samples (microarray samples plus additional pools). To prevent DNA contamination, one deoxyribonuclease treatment step (gDNA Wipeout Buffer, Qiagen Iberia S.L, Madrid, Spain) was performed from total RNA (1000 ng). Reverse transcription was then carried out using the Reverse Transcriptase Quantitect kit (Qiagen Iberia     S.L, Madrid, Spain) according to the manufacturer's instructions. Real-time qPCR (RT-qPCR) reactions were conducted in an Applied Biosystems 7500 (Applied Biosystems, Foster City, CA). Every PCR was performed with 5 mL of 1/10 diluted cDNA of each sample used in each reaction in a final volume of 20 mL of 10 mL of SYBR Green Master Mix (Applied Biosystems) and 200 nM of forward and reverse primers (list of RT-qPCR primers is shown in Table 1). The PCR protocol included an initial step of 50uC (2 min), followed by 95uC (10 min) and 40 cycles of 95uC (15 sec) and 60uC (1 min). After RT-qPCR, a melting curve analysis was performed by slowly increasing the temperature from 65uC to 95uC, with continuous recording of changes in fluorescent emission intensity. Serial dilutions of cDNA pool made from several samples were run in triplicate to assess PCR efficiency and decide which dilution to use for unknown samples. Target and reference genes in unknown samples were run in duplicate. Nontemplate controls (cDNA was replaced by water) for each primer pair were run in all plates. A DDC t method adjusted for PCR efficiency was used [18], employing the geometric average of H2AFZ and GAPDH as normalisation factor [19] and relative expression of cDNA pooled from various samples was used as a calibrator. The products of RT-qPCR were confirmed by

Statistical Analysis
Data were analysed using the Statgraphics version Plus 5.1 (Statistical Graphics Co., Rockville, MD, USA,) software package. The relative expression data were analysed using General Linear Model (GLM). For SMARCA2 a Neperian logarithmic transformation was done before analysis for data normalisation. Differences in mean values were tested using ANOVA followed by a multiple pair wise comparison using t-test. Differences of p,0.05 were considered to be significant.

Parthenote embryo production and blastocyst recovery
From the total of 369 oocytes activated and transferred to recipient does, 49 blastocysts properly developed were recovered at day 6 post-activation (13.3%). Sixty-four in vivo fertilised blastocysts were recovered at day 6 post-insemination (88.9% related to ovulation rate, estimated as the number forming corpora lutea).

Gene expression profiling and validation by real-time qPCR
PCA showed that samples from the same group clustered together (Figure 1). Analysis of expression data identified a total of 2541 differentially expressed transcripts between 6-day-old parthenotes and in vivo fertilised embryos. Among these, 1185 were upregulated whereas the 1356 remaining transcripts were downregulated. Table 2 shows a classification of differentially expressed transcript probes based on fold-changes. Specifically, parthenogenetic blastocysts exhibited changes in the expression of 92 genes, of which 16 had lower expression and 76 showed higher expression than in vivo fertilised embryos using a minimal 3-fold change as a cut-off. The lists of the upregulated and downregulated genes in the parthenogenetic blastocysts are shown in Table 3 and 4, respectively.
All genes selected to validate the microarray analysis exhibited expression patterns in line with previous results. Similarly, the three genes that exhibited lower expression in parthenotes in the microarray experiment (MPACT, DPY30 and CALC) also showed decreased expression by RT-qPCR (Table 5), while three genes showing higher expression in parthenogenetic blastocysts by the microarray analysis (SCGB1A1, EMP1 and SMARCA2) also  (Table 5). Comparisons between fold-change of results for RT-qPCR and microarray are shown in Table 5. The PCR experiments reproduced the microarray profiling for selected genes, although fold changes differed between RT-qPCR and microarray, which can be explained by different probes used for RT-qPCR and microarray [20].
Biological process, molecular function and cellular component vocabulary items assigned to upregulated and downregulated genes in parthenote embryos are shown in Figures 2, 3, and 4 respectively. For Biological Process, the most represented categories of altered genes were those related to cellular macromolecule process, transport, regulation of cellular process, protein metabolic process, nucleic acid metabolic process and macromolecule modifications ( Figure 2). As far as molecular function is concerned, the most represented GO terms were DNA and RNA binding, receptor binding and transferase activity (Figure 3). Finally, main annotations for cellular components are those related to mitochondrion, nuclear lumen, nucleus and cytoskeleton ( Figure 4).

Putatively imprinted genes
In parthenote embryos expression of paternally expressed imprinted genes is not expected, since both alleles are of maternal origin. We extracted information probes from the microarray data that detected known or putative imprinted genes (Catalogue of Imprinted Genes; http://igc.otago.ac.nz/home.html). Six of the genes which appear as most specifically upregulated or downregulated in the microarray have previously been annotated as imprinted genes. GRB10 and ATP10A were upregulated in parthenotes, as expected because the maternal allele is the one expressed, while ZNF215, NDN, IMPACT and SFMBT2 were downregulated according to the paternal allele expression. Furthermore, 26 other genes of the microarray which were significantly different in parthenote embryos, also shown to have at least one member of that gene family imprinted in other species (Table 6).

Discussion
Our results demonstrated that parthenotes and in vivo fertilised rabbit blastocysts cultured under in vivo conditions differ notably in gene expression. Up till now, few works have analysed transcriptome differences between parthenotes and fertilised embryos  [20,21,22]. However, these works were carried out with parthenote embryos developed in vitro and in vitro cultured fertilised embryos. It is well documented that embryos developed under in vitro environment are still not comparable with in vivo embryos [23], as post-fertilisation culture environment is a determinant for adequate embryonic development [4,24]. For example, one of the most critical time points of preimplantation embryogenesis is the major embryonic genome activation at which the embryo switches from using the mRNA and proteins derived from the maternal genome to those resulting from de novo transcription from the embryonic genome [25]. During that time, availability of transcription factors, which are regulated by cell cycle-dependent mechanisms, is required [26]. These mechanisms are strongly influenced by a change in environmental conditions and subsequently affect the embryonic development, with potentially severe effects on foetal, prenatal and postnatal viability [27]. Corcoran et al. [20] found that a total of 384 genes were differentially expressed between in vivo and in vitro derived blastocysts, the vast majority of them (almost 85%) being downregulated in in vitro developed embryos. Likewise, the effects of developmental environment on mRNA expression in parthenogenetic embryos have also been described [11] this way. To our best knowledge, this is the first report that compared the genome-wide gene expression profiles between rabbit parthenogenetic blastocysts and fertilised blastocysts developed in vivo.
Microarray analysis of parthenotes and fertilised embryos developed in vitro indicated differences in expression of 749 genes from mouse with 1.8 fold-changes as a cut-off [20], 24 genes for early embryos and 5 for expanded embryos from bovine with 1.5 fold-changes as a cut-off [22] and 56 genes from buffalo with 1.4 fold-changes as a cut-off [21]. In this study, we observed that 1606, 557 and 199 microarray probe signals were changed in the parthenogenetic blastocyst using a minimum of 1.5, 2.0 and 3.0 fold-changes as a cut-off, respectively. The 199 probe signals represent 92 genes, of which 16 had lower expression and 76 showed higher expression in parthenotes than fertilised embryos, developed in vivo. In the present study, in terms of biological process categories, slight differences are observed between transcript percentage of up and downregulated genes. However, the main categories altered, related to transport and protein metabolic process, comprise more upregulated than downregulated genes. Genes with high fold-changes such as BZND6, ANXAL, MYL4 are involved in transport, while protein metabolic process includes genes such as ClUS, PPIL6 or CIRL. In contrast, regarding molecular function and cellular components, a higher percentage of downregulated transcripts are comprised. In this case, the main altered categories are those related to DNA and RNA binding, both located in cellular nucleus and involving genes such as GTF2B (general transcription initiation factor IIb; X), CHURC1 (Churchill domain containing 1), XRCC2 (DNA repair protein XRCC2), HNRNPD (heterogeneous nuclear ribonucleoprotein D), SAFB2 (scaffold attachment factor B2) or NEIL3 (nei endonuclease VIII-like 3) among others. So, these results suggest a great deficiency of the machinery associated with transcription and translation which might hinder basic cell functioning and thereby pre-implantatory development of parthenogenotes. Similar results of the main categories altered in biological processes have been observed before in gene expression profile studies of in vitro developed parthenotes. Processes such as proteolysis, peptidolysis, protein amino acid phosphorylation and cell transport showed to be the most representative upregulated in parthenotes, while nucleic acid binding and metabolic process were representative of the higher percentage of donwregulated transcripts in parthenotes [20,21]. To date, more than 100 imprinted genes have been identified in mice and many of them are also imprinted in humans [29]. In livestock animals, imprinted genes have also been identified [30,31,32,33]. However, to our best knowledge, few genes have been identified as subject to genomic imprinting in rabbit. All imprinted genes show either maternal-specific or paternal-specific mono-allelic expression, and their proper expression is essential for normal development, foetal growth, nutrient metabolism and adult behaviour [34]. We extracted informative probes from the microarray data that detected known or putative imprinted genes (Catalogue of Imprinted Genes; http://igc.otago.ac.nz/home. html). Of the 32 putative genes analysed in this manner (table 6), 6 were identified as conserved between rabbits, humans and mice; they included GRB10, ATP10A, ZNF215, NDN, IMPACT and SFMBT2. GRB10, SNRPN and CDKN1 were also shown to be imprinted in a previous work carried out with in vitro developed parthenotes in mouse [20]. In fact, the use of microarrays to analyse imprinted genes provided results in the same direction as quantitative allelic pyrosequencing (QUASEP) analysis [30].
In conclusion, the resulting findings of this study revealed that even under the best developmental conditions, parthenogenetic and fertilised embryos at the late blastocyst stage are different, with at least 92 genes significantly and differentially expressed. These differences have been shown to affect basic functions such as DNA and RNA binding, nucleus, mitochondrion and transport, among others. ESCs may inherit the blastocyst level of transcripts, and the alterations observed in parthenogenetic embryos could therefore be maintained in pESCs derived from them. These alterations in gene expression call for further studies to evaluate whether and to what extent these modifications are unfavourable for ESC establishment and successive transplantation therapies. Furthermore, this work represents the first approach to the study of imprinted genes in rabbit. Hence, future research into imprinted genes might also include rabbits as alternative model systems.