Detection of transcriptional difference of porcine imprinted genes using different microarray platforms

Technical Reproducibility

All arrays were normalized (treatment of each array described in methods) and compared by pairwise correlations between technical replicates, with the average Pearson correlation coefficient given below. Affymetrix Human and Porcine arrays were both highly correlated between technical replicates indicating high technical reproducibility. The lower correlation between replicates of the Affymetrix Human array is likely due to the greater percentage of inherently more variable non-hybridizing probes caused by sequence divergence mismatches. Also, it was necessary to remove approximately 10% of the probes on the long glass oligonucleotide array due to printing defects prior to normalization [see Additional File 1]. After this procedure, it was found that although these spotted arrays performed relatively well, they nonetheless showed significantly more variability than the Affymetrix arrays. Decreased error variance due to high technical reproducibility is one of the key contributing factors to a platform's ability to identify differentially expressed genes (Figure 1).

Cross-species hybridization onto Affymetrix Human Genechips

Utilizing the design described above, porcine cRNA was hybridized to human Affymetrix arrays and data analyzed by filtering and subsequent analysis via a linear mixed model as described in the methods. Due to the effects of cross-species hybridization, the complexity of downstream analysis for Affymetrix Human U133+2.0 arrays in the context of cross-species hybridization was significantly greater than for the remaining arrays. Specifically, a high degree of variability within probe sets, likely due to probes with low sequence identity between human and pig was noted. The plots of probe expression profiles highlights the difficulty of this problem (Figure 2a,2b).

In the Affymetrix Human probe set 212092_at, targeting the gene PEG10 (Figure 2b), the small black arrows indicate that only 1^st, 2^nd, and 11^th probes appear to be differentially expressed in control versus parthenogenetic samples. The remaining probes show intensities that are randomly distributed around the median of the array. In contrast, Ssc.13476.1.A1_at, a porcine-specific probe for the same transcript (Figure 2a), showed clear evidence of differential expression across all probes. This inconsistent hybridization within probe sets is representative of the probes on the Affymetrix Human array under these cross-species hybridization conditions. Figure 2c further illustrates difficulties with the filtering process. In this figure, the dotted red line in the expression profile on the left is a filtering threshold. The probe set in the lower left shows results of filtering at a particular threshold (> one standard deviation from the mean intensity of the array). This illustrates how this type of filtering is imperfect as at a fairly stringent threshold, only two of three probe sets that exhibit evidence for differential expression are retained along with two probe sets that do not.

Estimates of Differential Expression

Differential expression for all three platforms was determined by estimating treatment effects after fitting a linear mixed model using SAS and JMP/Genomics via the method of Wolfinger et al.. [9]. The volcano plots (plotting estimate of treatment effects against the negative log of the p-value) (Figure 3) demonstrate that the number of significant differentially expressed genes identified varies greatly between these three platforms. All known imprinted transcripts have been highlighted in red; the dotted red lines correspond to a Bonferroni adjusted (p < 0.05) on the vertical axis and > 2-fold change on the horizontal. The Affymetrix Porcine array identified both the greatest number of differentially expressed genes as well as the greatest number of differentially expressed known imprinted genes. This same trend was exhibited when the threshold for differential expression was set at q<0.05 and q<0.20. Table 2 summarizes the transcriptional differences identified.

Sensitivity

A cumulative distribution plot shows the proportion of genes that are at or below a full range of p-value thresholds (Figure 4). It corroborates that at all significance thresholds, the Affymetrix Porcine array detects genes as being differentially expressed with higher frequency.

Identification of Sequence-Matched Probes

In order to assess the performance of the microarrays independently of coverage, we identified a set of 333 probe clusters common to all platforms by sequence-based mapping on the probe level. Briefly, we selected short oligonucleotide probes that mapped to long oligonucleotide probes with complete sequence identity and included them in a matching probe cluster where there were mappings for both short oligonucleotide Affymetrix microarrays to the single porcine glass long oligonucleotide microarray. 10,212 Affymetrix Porcine short oligonucleotides mapped to 5,452 unique Porcine Glass long oligonucleotides; 727 Affymetrix Human short oligonucleotides map to 520 unique Porcine Glass long oligonucleotides. The intersection of these two mappings results in the 333 probe clusters which are used to assess intra-platform reproducibility. 82 of these probe clusters contain short oligonucleotide sequences that match exactly between the two Affymetrix platforms; correlation coefficients are also calculated for this subset. (Figure 5)

Within the constraints of the technology (70-mer probes versus 25-mer probes), this sequence-based probe-to-probe mapping is the most rigorous method of identifying comparable matching probes to assess inter-platform reproducibility.

Inter/Intra-Platform Reproducibility

These correlation values compare favorably with the best correlations reported by Pylatuik et al. on biological replicates between platforms and are comparable to the correlations for A values obtained by Barczak et al. where technical replicates were used to compare the same sample between short and long oligonucleotide platforms (Table 3; [10, 11]).

Validation of known imprinted genes by qRT-PCR

Experimental Design

Gene expression profiles of fibroblast cell lines derived from day 27 control and parthenogenetic embryos were compared in each of the three platforms. For each platform, three biological replicates were used. Each biological replicate consisted of fibroblasts derived from a randomly selected fetus and cultured for two passages. One of the biological replicate was further split into three technical replicates. For biparental controls, sex of fetuses was determined by PCR and only female fetuses were used to avoid sex-related gene expression differences. For the technical replicates, one of the biological replicates was split into three identical pools of RNA and hybridized independently. For cross-platform comparisons, the same starting pool of total RNA was used to generate labeled targets for each of the three individual experiments. A balanced dye swap design was employed for the two-channel glass oligonucleotide microarray and one control and one parthenogenetic biological replicate were each divided into three technical replicates (Figure 6).

Generation of control pregnancies

Control crossbred gilts were mated by artificial insemination with boars to produce the biparental control fetuses for this study. Gilts were mated at 12 and 24 hr after their natural standing heat.

Generation of parthenogenetic pregnancies

Oocyte collection and maturation: Porcine ovaries were collected from sows at a local slaughterhouse and transported in 0.9% saline solution at 30–35°C. At the lab, ovaries were washed four times with warmed saline solution. Cumulus oocyte complexes (COCs) were aspirated from ovarian follicles 3–8 mm in diameter using a 5 ml syringe fitted with a short bevel 18-gauge needle. Follicular fluid was collected in 50 ml centrifuge conical tubes at room temperature. Collected COCs were washed three times in TLH-PVA medium. Oocytes with uniform cytoplasm and at least two layers of compacted cumulus cells were used for maturation. COCs were matured in TC199-Hepes supplemented with 10% porcine follicular fluid (pFF), 5 μg/ml insulin, 10 ng/ml EGF, 0.6 mM cysteine, 0.2 mM pyruvate, 25 μg/ml kanamycin and 5 μg/ml of each eCG and hCG. Fifty COCs were cultured in 500 μl medium in a 4-well Nunc dish at 38.5°C, 5% CO₂ in a humidified atmosphere. COCs were cultured for 22 hr before being changed to the same but eCG- and hCG-free culture medium and cultured for additional 15 hr [18].

Electrical activation of pig oocytes

After 40 hr of maturation in vitro, cumulus cells of IVM pig oocytes were removed by repeated pipetting in 0.1% hyaluronidase. Denuded oocytes were washed three times in Ca²⁺ free-NCSU23 medium [19] and then exposed for 5 min to activation medium consisting of 0.3 M mannitol, 0.05 mM MgSO₄, and 0.1 mM CaCl₂. Oocytes were then transferred between electrodes (1 mm apart) covered by 3 ml of the activation medium in a chamber connected to an electrical pulsing machine (BTX ECM 2001). Oocytes were stimulated by a single DC pulse of 150 V/mm for 100 μsecs. After activation, oocytes were washed twice in NCSU-23 medium supplemented with 0.4% BSA (IVC medium) and moved into IVC medium containing 10 μg/ml of cyclohexamide and incubated for 6 hr in this medium. Then, oocytes were washed three times in NCSU-23 medium supplemented with 0.05% BSA for transfer.

Embryo transfer into recipient

Activated oocytes were transferred into naturally cycling gilts on the first day of the standing estrous. Ventral laparotomy was performed and oocytes were transferred into the oviduct [18].

Collection of biparental and parthenogenetic fetuses

Pregnancies were confirmed by ultrasound two days before the collection on day 27. Fetuses were collected following euthanasia of the gilt and dissection of the reproductive tract. Fetuses were removed from their placenta, weighed and placed into 50 ml conical tubes containing DMEM supplemented with 10% fetal bovine serum (FBS). Tubes were kept on ice for transportation to the laboratory. Placentas were taken separately, weighed and placed in liquid nitrogen for later studies.

Isolation of fibroblasts from biparental and parthenogenetic fetuses

The head and viscera of fetuses were removed and the remaining tissue was minced with a sterile razor blade. The tissue was added to 10 ml of 0.05% trypsin (Gibco) supplemented with 0.9 mM potassium chloride, 0.9 mM dextrose, 0.7 mM sodium bicarbonate, 0.1 mM EDTA (all from Sigma), and 20 mM sodium chloride (EMD Bioscience). The tissue/trypsin solution was shaken at 37°C for 15 min a total of three times. After incubation, the supernatant was collected, pooled, and pelletted. The cell pellet was resuspended in DMEM/F12 media (Gibco) supplemented with 10% FBS and 5% calf serum (CS) (both from Hyclone), 30 mM sodium bicarbonate, 0.5 mM pyruvic acid, and 2 mM N-acetyl-L-cysteine (all from Sigma). In addition, 100 units penicillin and 100 ug streptomycin,(Gibco) were added per 100 ml media to inhibit microbial growth. The cells were placed in the appropriate number of 10 cm tissue culture plates (Corning), incubated in a 5% CO₂ incubator at 37°C, expanded once and frozen in 50% FBS, 40% media, and 10% DMSO (Sigma) for long time storage and future use.

Determination of sex of fetuses by PCR

The sex of the fetuses from which each of the biparental control fibroblast cell lines was derived was determined by SRY genotyping using the following primers: 5'-TGAACGCTTTCATTGTGTGGTC-3', 5'-TCCTCCGTGTCTCTGATGACCG-3' [20]. The PCR thermocycling conditions were 95°C for 2 min, 35 cycles of 95°C for 20 s, 55°C for 30 s, 72°C for 1 min, followed by 72°C for 7 min.

RNA Isolation

Cells derived from biparental female and parthogenetic fetuses were grown in 10 cm tissue culture plates in DMEM/F12 media (Gibco) supplemented with 10% FBS and 5% calf serum CS (both from Hyclone), 30 mM sodium bicarbonate, 0.5 mM pyruvic acid, and 2 mM N-acetyl-L-cysteine (all from Sigma). At 90% confluency, the RNA was extracted using RNAqueous Kit (Ambion) as per the instructions of the manufacturer. Briefly, media was removed from the plates and cells were lysed in 1 ml lysis buffer. To this was added 1 ml of 64% alcohol. The contents were mixed and passed through the column. The RNA bound to the column was washed once with wash solution 1 and twice with wash solution 2/3. Finally RNA was eluted in 40 μ l of hot elution buffer. RNA was quantified by spectrophotometry and quality verified by running 5 μ g of RNA on 1% agarose gel. The resulting RNA was used for microarray analyses.

Target Production and Hybridization: Affymetrix Human and Porcine arrays

Before target production, the quality and quantity of each RNA sample was assessed using a 2100 BioAnalyzer (Agilent). Target was prepared and hybridized according to the Affymetrix Technical Manual. Total RNA (10 ug) was converted into cDNA using Reverse Transcriptase (Invitrogen) and a modified oligo(dT)₂₄ primer that contains T7 promoter sequences (GenSet). After first strand synthesis, residual RNA was degraded by the addition of RNaseH and a double-stranded cDNA molecule was generated using DNA Polymerase I and DNA Ligase. The cDNA was then purified and concentrated using phenol:chloroform extraction followed by ethanol precipitation. The cDNA products were incubated with T7 RNA Polymerase and biotinylated ribonucleotides using an In Vitro Transcription kit (Enzo Diagnostics). One-half of the cRNA product was purified using an RNeasy column (Qiagen) and quantified with a spectrophotometer. The cRNA target (20 ug) was incubated at 94°C for 35 min in fragmentation buffer (Tris, Magnesium Acetate, Potassium Acetate). The fragmented cRNA was diluted in hybridization buffer (MES, NaCl, EDTA, Tween 20, Herring Sperm DNA, Acetylated BSA) containing biotin-labeled OligoB2 and Eukaryotic Hybridization Controls (Affymetrix). The hybridization cocktail was denatured at 99°C for 5 min, incubated at 45°C for 5 min and then injected into a GeneChip cartridge. The GeneChip array was incubated at 42°C for at least 16 hr in a rotating oven at 60 rpm. GeneChips were washed with a series of nonstringent (25°C) and stringent (50°C) solutions containing variable amounts of MES, Tween20 and SSPE. The microarrays were then stained with Streptavidin Phycoerythrin and the fluorescent signal was amplified using a biotinylated antibody solution. Fluorescent images were detected in a GeneChip^® Scanner 3000 and expression data were extracted using the MicroArray Suite 5.0 software (Affymetrix). All GeneChips were scaled to a median intensity setting of 500.

Target Production and Hybridization: U.S. Pig Genome Coordination Program Glass Arrays

RNA was extracted from primary cultures of control and gynogenote fibroblasts using RNAqueous^® (Ambion) following the manufacturer's suggested protocol and stored at -80°C. One microgram of purified RNA was converted to aminoallyl-coupled RNA (aRNA) and coupled with Cy3 or Cy5 using Amino Allyl Message Amp II aRNA Kit (Ambion), again suggested protocols were followed. Specific activity and aRNA concentration of the purified labeled aRNA was determined by assaying one μl of sample on a NanoDrop^® ND 1000. Specific activities (pmol dye/pmol aRNA) of all probes were between 25 and 40. Control and parthenogenetic probes were pooled so that equal molar amounts of each dye were used per array. Pooled probes were dried to completion using an Eppendorf Vacufuge then fragmented using Fragmentation Reagent (Ambion). Fragmented probes were dried to a 10 μl volume and immediately used for hybridization. Glass arrays were generated at the University of Minnesota Microarray Printing Facility and obtained through the U.S. Pig Genome Coordination Program. Arrays were used within two weeks of receipt. Slides were pre-hybridized, hybridized and washed according to GAPS II Coated Slides instruction manual (Corning) with the exception that 300–400 picomoles of dye were used per slide and 0.1 μg/μl of Thymus DNA (Sigma) was used in the hybridization buffer. The arrays were scanned with ScanArray Express (Packard Bioscience). Acquired images were analyzed using QuantArray software version 3.0 (Packard Bioscience). During the quantification process, approximately 10% of probes were discarded due to visually identified printing defects on the arrays.

Normalization and Filtering

For Affymetrix arrays, probe intensity values were log2 transformed and quantile normalization was applied [21]. The average of the three technical replicates was taken to determine the probe intensities for the corresponding biological replicate. The Affymetrix Human arrays were treated as a special case due to the effects of sequence divergence on cross-species hybridization. For these arrays, we employed quantile normalization, and then subsequently corrected for the increased variability of probe expression profiles by filtering out non-hybridizing probes from within probe sets. We tried two approaches to filtering, one solely based on the intensity of the perfect match probe, and the second based on both the difference and the ratio between the perfect match and the mismatch probe as described by Ji et al. [4]. In the perfect match only approach, we filtered out non-hybridizing probes which did not exceed arbitrary filtering thresholds in any of the samples. We tried four filtering thresholds: 0, the median array intensity, the third quartile, and one standard deviation above the mean. In the perfect match and mismatch approach, we implemented a filter at PM - MM > 200, PM/MM > 2.

The spotted glass oligonucleotide arrays were normalized with a lowess normalization with a smoothing parameter of 0.2 to broadly correct for dye effects.

Statistical Analysis

The Affymetrix arrays were fitted to the following gene by gene linear mixed model using SAS and JMP/Genomics (Cary, NC) [9].

y _ijk = μ + T _i + P _j + A _k + ε _ijkl

For each probe set, y is the log2 transformed intensity of the i^th treatment, j^th probe, and k^th array. This model included fixed effects for treatment (control or parthenote, T) and probe (P) and random effects for array (A).

The glass spotted oligonucleotide array was fitted to the following gene by gene linear mixed model.

y _ijk = μ + T _i + D _j + A _k + D _j *A _k + ε _ijkl

For each probe, y is the log2 transformed intensity of the i^th treatment, j^th dye, and k^th array. This model included fixed effects for treatment (control or parthenote, T), dye (D), the interaction between dye and treatment (D*A), and random effects for the array (A).

Least square means were estimated for the difference between treatments for each gene. In the Affymetrix Human array, corresponding p-values were converted to q-values by a method proposed by Storey that measures significance in terms of false discovery rate to optimize filtering thresholds [22]. For all arrays, p-values were adjusted with a Bonferroni correction to control the family wide error rate to <0.05.

Correlations for technical reproducibility

The control technical replicates for each of the three array platforms were compared by standard pairwise correlations. The average Pearson correlation coefficient for these three arrays is reported. The Cy5 channel of the glass array is used for comparison purposes, but both channels have similar correlation values.

Cumulative distribution of p-values

An empirical cumulative distribution function was fitted for each of the three sets of p-values. A plot of p-value by frequency was then constructed, where each point represents a gene with its corresponding p-value.

Identification of Sequence-Matched Probe Clusters

Mecham et al. have suggested that the lack of concordance in cross-platform microarray comparisons may be caused by a reliance on gene annotations without more stringent sequence-oriented matching of probes [23]. We downloaded probe sequences represented on the microarrays from Affymetrix [24, 25], and Operon [26]. Using this sequence information, probe sequences were matched at the probe level by mapping Affymetrix short oligonucleotides to porcine spotted glass long oligonucleotide sequences using the BLAST standalone program as described by Kuo et al. [27]. Using the long oligonucleotide sequences as a reference, probe clusters were identified where there is a match with complete sequence identity between the long oligonucleotide sequences and short oligonucleotide sequences for the span of the short oligonucleotide sequence for both Affymetrix microarrays (matching probe clusters). In the cases where there was more than one short oligo probe per microarray per cluster, the average of the normalized expression intensities was taken. A subset of matching probes with complete overlap were identified where the short oligonucleotide sequences on the two Affymetrix microarrays have identical sequence.

Intra/Inter-platform correlations

Average Pearson and Spearman pairwise correlation coefficients were calculated on the normalized expression intensities of the control technical replicates using JMP (Cary, NC). Additionally, Pearson and Spearman pairwise correlation coefficients were calculated for the subset of matching probes with complete overlap.

Annotation

Since the two porcine microarrays were minimally annotated, both were reannotated by BLAST against an EnsEMBL Human cDNA sequence library. This annotation was enhanced by using information from The Institute for Genome Research (TIGR) Pig Gene Index [17]. Briefly, we attempted to extend the target sequences by matching them to TIGR assembled porcine consensus sequences. These extended target sequences were compared to a library of EnsEMBL human cDNA sequences by BLAST, and the gene with the highest bit score was recorded. The same procedure was repeated for the original unextended target sequences. A subset of these original sequences with bit scores greater than 50 were evaluated for concordance with the extended target sequences and resulted in >96% agreement.

Validation by qRT-PCR

Gene transcripts were quantified by real-time reverse transcription PCR using the iCycler apparatus (BioRad Inc., Hercules, CA) and were detected with SYBR Green I as fluorochrome (Platinum SYBR Green I; Invitrogen, Carlsbad, CA). The primers used for PCR are listed in Additional File 2. The relative expression changes were determined with the $2^{-\Delta \Delta {\text{C}}_{\text{t}}}$ method, where $\Delta {\text{C}}_{\text{t}}={\text{C}}_{{\text{t}}_{\text{target gene internal reference}}}$, and $\Delta \Delta {\text{C}}_{\text{t}}=\Delta {\text{C}}_{{\text{t}}_{\text{sample}}}-\Delta {\text{C}}_{{\text{t}}_{\text{calibrator}}}$. 18S was used as internal reference gene. PCR efficiency was tested for each primer pair by 10-fold dilution series of cDNA in triplicate to make sure that efficiency is appropriate for the 2^-ΔΔCt Pfaffl et al. method [28]. To ensure the specificity and integrity of the PCR product, melt-curve analyses were performed for all PCR products. No PCR products were obtained from RNA samples when RT was omitted. Samples without template for each primer pairs were included to identify contamination. The experimental design was executed in triplicate for each control and parthenote combination.