Gene expression microarray analysis of early oxygen-induced retinopathy in the rat.

Different inbred strains of rat differ in their susceptibility to oxygen-induced retinopathy (OIR), an animal model of human retinopathy of prematurity. We examined gene expression in Sprague-Dawley (susceptible) and Fischer 344 (resistant) neonatal rats after 3 days exposure to cyclic hyperoxia or room air, using Affymetrix rat Genearrays. False discovery rate analysis was used to identify differentially regulated genes. Such genes were then ranked by fold change and submitted to the online database, DAVID. The Sprague-Dawley list returned the term "response to hypoxia," absent from the Fischer 344 output. Manual analysis indicated that many genes known to be upregulated by hypoxia-inducible factor-1alpha were downregulated by cyclic hyperoxia. Quantitative real-time RT-PCR analysis of Egln3, Bnip3, Slc16a3, and Hk2 confirmed the microarray results. We conclude that combined methodologies are required for adequate dissection of the pathophysiology of strain susceptibility to OIR in the rat. ELECTRONIC SUPPLEMENTARY MATERIAL: The online version of this article (doi:10.1007/s12177-009-9041-7) contains supplementary material, which is available to authorized users.

Introduction

Retinopathy of prematurity (ROP) is a blinding condition of some premature infants who require neonatal intensive care and supplemental oxygen therapy [1]. Severe retinal disease affects approximately 5% of these neonates, and despite best available care, 10–15% will suffer significant visual impairment [2]. ROP shows variation in incidence, progression, and disease severity among premature neonates, and strong evidence exists to indicate a heritable predisposition [2-4]. Thus, in a retrospective study of monozygotic and dizygotic twins, 70% of the variance in susceptibility to ROP was found to be the result of genetic factors alone, even after accounting for the effects of gestational age and duration of supplemental oxygen [5]. Furthermore, many studies have identified an influence of ethnicity in susceptibility to ROP [2, 6–8].

The retinae of small rodents and some other mammals are largely avascular at birth, and retinal vascularization occurs in the early neonatal period. Neonatal rodents exposed to hyperoxia and that are then returned to the relative hypoxia of room air develop rapid, reliable neovascularization and the development of microvessel abnormalities that closely resemble human ROP. Models of oxygen-induced retinopathy (OIR) have long informed our understanding of the pathophysiology of ROP [9]. A heritable component to ocular angiogenesis has been well documented [10], and we have shown that inbred strains of rat exhibit differential susceptibility to OIR [11]. Three albino strains, Lewis, Fischer 344 (F344), and Wistar–Furth, were found to be resistant to the effects of hypoxia following hyperoxic exposure, whereas the albino Sprague–Dawley (SD) strain was relatively susceptible, and the pigmented Dark Agouti (DA) strain was very sensitive. The last two strains develop florid microvascular abnormalities and neovascular tufts. These differences in susceptibility are independent of litter size, body mass, arterial oxygen tension, and major histocompatibility complex haplotype. Genetic backcross analysis of the susceptibility trait suggested an autosomal dominant pattern of inheritance in the F344 and DA strains [12], confirming that susceptibility to OIR in the rat is a heritable trait. Further, we and others have demonstrated that disease susceptibility correlates with differential expression of several angiogenesis-related genes very early in retinal microvascular development [13, 14].

To identify the pathways responsible for the difference in retinal vascular phenotype exhibited by different albino rat strains following exposure to cyclic hyperoxia and subsequent hypoxia, we performed an initial screen of retinal mRNAs from F344 (resistant to OIR) and SD (moderately susceptible to OIR) at an early time-point (day 3), using the Affymetrix GeneChip rat arrays covering over 27,000 genes. Quantitative real-time reverse-transcription polymerase chain reaction (qRT-PCR) was used to confirm differential expression of some candidates, and the online Database for Annotation, Visualization, and Integrated Discovery (DAVID) was used for pathway analysis [15, 16].

Materials and methods

Experimental animals

Inbred F344 and SD rats were obtained from the institutional animal facility. Rats were allowed unlimited access to rat chow and water and were exposed to a 12-h light–dark cycle. Room temperature was maintained at 24°C. All animal experimentation was approved by the institutional Animal Welfare Committee.

Exposure of neonatal rats to cyclic hyperoxia

The oxygen exposure protocol used in this study is a modification of those used in other studies of OIR in the rat [17-21]. Within 12 h of birth, female rats and their litters were placed in a humidified chamber and exposed to alternating 24-h cycles of hyperoxia (80% oxygen in air) and normoxia (21% oxygen in air) until the rats were euthanized [11]. This protocol was chosen because it reliably induces retinopathy in susceptible strains, without evidence of systemic toxicity [11]. Age and strain-matched room air (normoxia)-exposed rats were used as controls.

Retinal dissection and histochemistry

Rats were killed with an inhaled overdose of isoflurane anesthesia on post-natal day 3 (P3), after the second 24-h exposure to hyperoxia. Retinal dissections and histochemistry labeling were performed as previously described [11]. Briefly, eyes were enucleated and fixed in 2% weight/volume paraformaldehyde in phosphate-buffered saline, pH 7.4 for 90 min on ice. The retinae were dissected, and four incisions were made to flatten the tissue. Within 6 h of enucleation, the dissected retinae were stained for 18 h with fluorochrome-conjugated Griffonia simplicifolia isolectin B4 (GS-IB4; Alexa Fluor 488 conjugate; Molecular Probes, Eugene, OR, USA) to highlight the vasculature.

Image analysis of isolectin-labeled retinae

Imaging of labeled retinae was performed within 12 h of staining. Images were captured using a fluorescence microscope (Olympus Optical Co. Ltd., Tokyo, Japan) attached to a CCD-digital camera (Roper Scientific, Trenton, NJ, USA) and image acquisition software (CoolSNAP, Roper Scientific). Serial, overlapping high-resolution images of the entire retina were taken using a ×4 objective lens. Images were then merged to create a montage (Adobe Photoshop CS3, Adobe systems Inc., San Jose, CA, USA).

Preparation of RNA

At P3, rats were removed from the oxygen chamber and euthanized with an inhaled overdose of isoflurane anesthetic. The eyes were immediately enucleated and placed into chilled diethylpyrocarbonate (DEPC, Sigma, St. Louis, MO, USA)-treated saline. Total RNA preparations were made from whole retinae with the retinal pigment epithelium attached and snap-frozen in liquid nitrogen. The total time in normoxia from removal of rat pups from the oxygen chamber to snap-freezing of the retinae was about 20–30 min. Tissue handling, thereafter, was kept to a minimum to prevent RNA degradation. Total RNA was isolated from pooled retinae of individual rats using TRIzol (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s method. Total RNA was quantified using the Nanodrop 8000 (Thermo Scientific, Wilmington, DE, USA). The absorbance ratio of the RNA at 260 and 280 nm was used as a guide in determining RNA purity. Samples with a ratio of 1.9 or more were considered to be acceptable. RNA integrity and residual DNA contamination was assessed using agarose gel electrophoresis. RNA with a 28S:18S ratio of approximately 2:1 and minimal DNA contamination was considered acceptable for use. RNA was further purified using a Qiagen RNeasy Micro column purification kit (Qiagen, Hilden, Germany) according to the manufacturer’s protocol.

Preparation of test RNA samples for mRNA microarray analysis

Retinae from post-natal day-3 rats (P3) were collected for each experimental condition as shown in Supplemental Material Fig. S1. Briefly, four experimental conditions were tested: F344 room air-exposed (RA), F344 cyclic hyperoxia-exposed (O2), SD RA, and SD O2. For each of these experimental conditions, total RNA from a minimum of nine rats from up to three different litters was collected. A proportion of several RNA samples were used to create a test pool consisting of three randomly selected rats. RNA from the remaining rats was used to make a second and third test sample as required, using RNA from a maximum of three rats per sample. These acted as biological replicates for the same experimental conditions.

Affymetrix mRNA microarrays

Samples were analyzed at the Adelaide Microarray Centre (University of Adelaide and Hanson Institute, South Australia, Australia). Prior to arraying, samples were analyzed to determine RNA integrity using the Agilent Bioanalyser 2100 (Agilent Technologies, Santa Clara, CA, USA). Samples with a RNA integrity number of seven or above were deemed suitable for further use in mRNA microarrays.

Two samples were provided for each of the four experimental strain/treatment conditions. Sample grouping was initially assessed using a principal component analysis (PCA) plot (Partek Genomics Suite) of the corrected data (Supplemental Material Fig. S2). As a result of the original PCA, additional samples for F344 O2 and SD RA were submitted for analysis, and the enhanced dataset was reanalyzed by PCA (Supplemental Material Fig. S3). A maximum of three samples were used for each experimental condition. Single-color arrays allowed for comparisons to be made among all experimental conditions. A total of six comparisons were made: F344 versus SD (indifferent to treatment), O2 versus RA (indifferent to strain), F344 RA versus SD RA, F344 O2 versus SD O2, F344 O2 versus F344 RA, and SD O2 versus SD RA.

Affymetrix GeneChip Gene 1.0 ST Rat Arrays (Affymetrix Inc., Santa Clara, CA, USA) covering 27,342 genes using 722,254 distinct 25-mer oligonucleotide probes were used for gene expression analysis.

Data analysis of microarrays

Affymetrix raw gene array data were processed using the Partek Genomics Suite software (Partek Inc., St. Louis, MO, USA). Briefly, .cel files were imported using RMA background correction, Partek’s own GC content correction, and mean probe summarization. Differential gene expression was assessed by ANOVA with the p value adjusted using step-up [22] multiple test correction to control the false discovery rate (FDR). Adjusted p values < 0.05 were considered to be significant.

A bioinformatics approach was used to determine the biological context of the large amounts of gene expression data generated by the microarray screen. Gene lists from comparisons showing significant differences in gene expression were submitted to the freely available database DAVID [15, 16]. Pathways showing enrichment for specific genes were found using the functional annotation chart pathway-mining tool.

Preparation of cDNA for RT-PCR

One microgram of RNA pools previously prepared for microarray analysis was DNase-treated with Turbo DNA-free according to the manufacturer’s instructions (Ambion, Austin, TX, USA) and reverse-transcribed in 20-μl reactions using a first-strand cDNA synthesis kit (SuperScript III First-Strand Synthesis System, Invitrogen). An RT-free control sample was prepared in parallel with each cDNA sample, with substitution of DEPC-treated water for reverse transcriptase. For purposes of standardization, a standard cDNA pool was prepared from pooled retinal RNA extracted from 10 F344, SD, and DA rats from a range of developmental stages that had been exposed to room air or cyclic hyperoxia. Samples of cDNA were frozen at −20°C.

Primers for RT-PCR

Primers for candidate genes were designed to flank an intron using Primer3 software [23] and tested in silico for specificity against sequences for Rattus norvegicus using BLAST software (NCBI, Bethesda, MD, USA) (Table 1). Primers were synthesized by Geneworks Ltd (Thebarton, South Australia, Australia). Primer pairs were aligned with the rat genome (NCBI BLAST software) to identify the potential for genomic DNA amplification and were considered acceptable if the genomic amplicon was more than twice the length of the target cDNA amplicon. Acidic ribosomal phosphoprotein (Arbp) and hypoxanthine guanine phosphoribosyl transferase (Hprt) were used as housekeeping reference genes.

Table 1

Primer sequences for qRT-PCR

Gene (accession number)	Primer sequence (5′–3′)	Nucleotide position	Amplicon size (bp)
EGL nine homolog 3	TTGGGACGCCAAGTTACATG (for)	1003–1078	76
(Egln3) (NM 019371.1)	GGGCTCCACGTCTGCTACAA (rev)
BCL2/adenovirus E1B 19 kDa-interacting protein 3	GCGCACAGCTACTCTCAGCA (for)	478–627	150
(Bnip3) (NM 053420.2)	GTCAGACGCCTTCCAATGTAGA (rev)
Solute carrier family 16, member 3	ATGTGTGTGAACCGCTTTGG (for)	327–449	123
(Slc16a3) (NM 030834.1)	GACCCCTGTGGTGAGGTAGATC (rev)
Hexokinase 2	CAACATTCTCATCGATTTCACGAA (for)	2463–2607	145
(Hk2) (NM 012735.1)	GATGGCACGAACCTGTAGCA (rev)
Acidic ribosomal phosphoprotein	AAAGGGTCCTGGCTTTGTCT (for)	766–856	91
(Arbp) (NM 022402.1)	GCAAATGCAGATGGATCG (rev)
Hypoxanthine guanine phosphoribosyl transferase	TTGTTGGATATGCCCTTGACT (for)	629–733	105
(Hprt) (NM 012583.2)	CCGCTGTCTTTTAGGCTTTG (rev)

Quantitative RT-PCR

qRT-PCR was performed with reference to the MIQE guidelines [24]. Each 20-µl reaction mixture contained 10 µl of SYBR Green master-mix (Power SYBR Green PCR Master Mix; Applied Biosystems), 2 µl each forward and reverse primers (0.5 µM final concentration), and 6 µl cDNA sample diluted 1/100 with Ultra Pure water (Fisher Biotech, Western Australia, Australia). Reactions were run on a StepOnePlus Real-Time PCR System (Applied Biosystems). Cycling conditions were: initial denaturation (95°C, 10 min) followed by 40 cycles of denaturation (95°C, 15 s), annealing, and extension (60°C, 1 min). Candidate genes were tested in triplicate wells and in duplicate experiments, giving a total of six replicates for each sample pool per gene. The standard cDNA pool was included in triplicate in each PCR run. A single RT-free negative control for each sample and two water controls (no-template control) were included in each experiment. Melt-curve analysis was used to confirm amplicon specificity. The melt curve of each PCR product was compared with that of the corresponding sequenced product. Representative PCR products of all genes were separated on agarose gels, purified, sequenced, and compared with the predicted amplicon sequence to confirm identity. DNA was extracted from agarose gels using the Qiaquick column purification system (Qiagen) and eluted in 10 mM Tris–Cl pH 8.5 for quantification and sequencing. Purified DNA was labeled using the BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems) and resolved using the ABI 3100 Genetic Analyser (Applied Biosystems).

Analysis of qRT-PCR data

Gene expression in each sample was determined relative to the standard cDNA pool and normalized to the reference genes Arbp and Hprt using the method of Vandesompele et al [25]. Arbp and Hprt are stably expressed under conditions of cyclic hyperoxia and during early development in a number of rat strains [26].

Results

Retinal vascularization in post-natal day-3 rats exposed to cyclic hyperoxia

Montages of representative sectors from the retinae of 3-day-old F344 and SD rats exposed to cyclic hyperoxia from birth or raised in room air are shown in Fig. 1. During normal retinal vascular development, vessels develop from the optic nerve head (indicated by arrow) towards the peripheral retina. At post-natal day 3, vessels in the room air-exposed and cyclic hyperoxia-exposed F344 rats, and in the room air-exposed SD rats, had begun to develop and showed normal morphology (Fig. 1a–c). SD rats exposed to cyclic hyperoxia showed aberrant vascular development, with fewer vessels radiating from the optic nerve head and reduced density of the capillary network (Fig. 1d).

Fig. 1

Montages of representative sectors of retinae of room air-exposed and cyclic hyperoxia-exposed rats collected at post-natal day 3 and stained with fluorochrome-conjugated GS-IB4 to highlight the vasculature. a F344 rat raised in room air; b F344 rat exposed to cyclic hyperoxia from birth; c SD rat raised in room air; d SD rat exposed to cyclic hyperoxia from birth

Changes in gene expression in response to cyclic hyperoxia exposure and/or strain

The Affymetrix microarrays were used to investigate changes in gene expression in response to cyclic hyperoxia exposure and/or strain. The FDR represents the number of false positives that may be present within the dataset. The number of genes that were significantly overexpressed or underexpressed in each comparison is shown in Table 2 at a FDR of 5%. The majority of significant differences in gene expression were observed between strains. When the false discovery rate was set at 10% or 20%, no differences were found in gene expression relating to both strain and treatment combined (data not shown).

Table 2

False discovery rate report with significance levels set at 5%

Experimental condition	Number of gene with significant p values
p Value (Strain; F344 versus SD)	251
p Value (Treatment; O2 versus RA)	0
p Value (Strain * Treatment)	0
p Value (F344 O2 versus SD O2)	10
p Value (F344 RA versus SD RA)	11
p Value (F344 O2 versus F344 RA)	0
p Value (SD O2 versus SD RA)	0

Total number of p values = 27,342

O exposed to cyclic hyperoxia; RA exposed to room air

Comparisons of gene expression between F344 and SD rats showed 251 genes to be significantly regulated by strain at p < 0.05 (adjusted for the number of comparisons; see Supplemental Material Table S1 for entire list). The top 20 differentially regulated genes are shown in Table 3. Using DAVID, genes related to strain were found to be involved in processes such as peptide transport, cell localization, immune response, and nervous system development (Table 4). When room air exposure was taken into consideration, 11 genes were found to be significantly differentially regulated in their fold changes (Table 5).

Table 3

Top 20 of 251 genes significantly regulated by strain

Gene	Gene description	Fold change	Fold change description	Adj. p value
Mospd1a	Motile sperm domain containing 1	22.11	F344 up vs. SD	0.002
RGD1564739	Akin to spermatogenesis-associated glutamate (E)-r	−3.62	F344 down vs. SD	0.002
Ly6g6e	Lymphocyte antigen 6 complex, locus G6E	3.84	F344 up vs. SD	0.002
Slc26a7	Solute carrier family 26, member 7	−4.42	F344 down vs. SD	0.002
Hey2	Hairy/enhancer-of-split related to YRPW motif 2	1.27	F344 up vs. SD	0.003
Tap2	Transporter 2, ATP-binding cassette, sub-family B	5.60	F344 up vs. SD	0.003
10722435b	Unknown	−19.64	F344 down vs. SD	0.005
Stk32a	Serine/threonine kinase 32A	−5.33	F344 down vs. SD	0.005
RGD1562660	RGD1562660	−1.30	F344 down vs. SD	0.006
RGD1560289	Similar to chromosome 3 open reading frame 20	−2.38	F344 down vs. SD	0.008
Arsb	Arylsulfatase B	−2.50	F344 down vs. SD	0.008
RT1-Bb	RT1 class II, locus Bb	−3.59	F344 down vs. SD	0.008
Ly6g6d	Lymphocyte antigen 6 complex, locus G6D	3.50	F344 up vs. SD	0.010
10804316b	Unknown	−1.72	F344 down vs. SD	0.010
Oprk1	Opioid receptor, kappa 1	2.00	F344 up vs. SD	0.010
10824257b	Unknown	2.69	F344 up vs. SD	0.010
Ank2	Ankyrin 2, neuronal	−1.27	F344 down vs. SD	0.011
10724713b	Unknown	−1.76	F344 down vs. SD	0.011
Fgd3	FYVE, RhoGEF and PH domain containing 3	2.32	F344 up vs. SD	0.011
10900112b	Unknown	−3.56	F344 down vs. SD	0.011

Adjusted p value < 0.05; ± indicates overexpression or underexpression of the gene, respectively. For the entire gene list, see Supplemental Material Table S1

aGene is represented by a second exon in the array with a fold change of 17.45 and an adjusted p value of 0.002

bAffymetrix probeset ID of a gene which has not been identified

Table 4

DAVID analysis for Gene Ontology terms enriched in the strain comparison ranked by fold enrichment

GO term	Number of genes	Fold enrichmenta	False discovery rateb
Peptide transport	6	14.2	0.1107
Cell localization	10	3.4	4.5728
Immune response	11	2.9	8.7035
Nervous system development	16	2.8	0.9572

aThe ratio between the frequency of term-annotated genes in the query gene list compared to the reference gene list, in this case the rat genome

bThe expected percentage of false positives in a given proportion of genes considered to be significant

Table 5

Genes significantly regulated (adjusted p value < 0.05) by strain in room air-exposed F344 and SD rats

Gene	Gene description	Fold change	Fold change description	Adj. p value
Mospd1a	Motile sperm domain containing 1	22.43	F344 RA up vs. SD RA	0.015
RGD1564739	Akin to spermatogenesis associated glutamate (E)-r	−3.98	F344 RA down vs. SD RA	0.015
Slc26a7	Solute carrier family 26, member 7	−4.10	F344 RA down vs. SD RA	0.022
Ly6g6e	Lymphocyte antigen 6 complex, locus G6E	3.25	F344 RA up vs. SD RA	0.026
Tap2	Transporter 2, ATP-binding cassette, sub-family B	5.56	F344 RA up vs. SD RA	0.029
Hey2	Hairy/enhancer-of-split related with YRPW motif 2	1.23	F344 RA up vs. SD RA	0.035
10722435b	Unknown	−20.69	F344 RA down vs. SD RA	0.035
Stk32a	Serine/threonine kinase 32A	−4.99	F344 RA down vs. SD RA	0.041
Oprk1	Opioid receptor, kappa 1	2.34	F344 RA up vs. SD RA	0.041
Arsb	Arylsulfatase B	−2.66	F344 RA down vs. SD RA	0.048
RGD1310110	Similar to 3632451O06 Rik protein	−1.44	F344 RA down vs. SD RA	0.048

± indicates over- or underexpression of the gene, respectively

aGene is represented by a second exon in the array with a fold change of 17.45 and an adjusted p value of 0.015

bAffymetrix probeset ID of a gene which has not been identified

A comparison of treatment groups (cyclic hyperoxia and room air-exposed rats), irrespective of strain, showed minimal fold changes in gene expression, and adjusted p values indicated the fold changes were insignificant (Supplemental Material Table S2). A comparison of cyclic hyperoxia-exposed F344 and SD rats showed that ten genes were differentially expressed to a significant extent (Table 6). The majority of these genes were represented in the previous comparisons.

Table 6

Genes significantly regulated (adjusted p value < 0.05) by strain in cyclic hyperoxia-exposed (O2) F344 and SD rats

Gene	Gene description	Fold change	Fold change description	Adj. p value
Mospd1a	Motile sperm domain containing 1	21.80	F344 O2 up vs. SD O2	0.012
Ly6g6e	Lymphocyte antigen 6 complex, locus G6E	4.55	F344 O2 up vs. SD O2	0.012
Slc26a7	Solute carrier family 26, member 7	−4.77	F344 O2 down vs. SD O2	0.012
Hey2	Hairy/enhancer-of-split related with YRPW motif 2	1.32	F344 O2 up vs. SD O2	0.012
RGD1564739	Akin to spermatogenesis associated glutamate (E)-r	−3.30	F344 O2 down vs. SD O2	0.016
Tap2	Transporter 2, ATP-binding cassette, sub-family B	5.64	F344 O2 up vs. SD O2	0.024
10722435b	Unknown	−18.64	F344 O2 down vs. SD O2	0.033
Stk32a	Serine/threonine kinase 32A	−5.69	F344 O2 down vs. SD O2	0.033
RGD1562660	RGD1562660	−1.33	F344 O2 down vs. SD O2	0.033
RGD1560289	Similar to chromosome 3 open reading frame 20	−2.50	F344 O2 down vs. SD O2	0.042

± indicates over- or underexpression of the gene, respectively

aGene is represented by a second exon in the array with a fold change of 17.46 and an adjusted p value of 0.012

bAffymetrix probeset ID of a gene which has not been identified

No significant changes in gene expression were seen in F344 rats exposed to cyclic hyperoxia for 3 days compared to room air-exposed rats (Supplemental Material Table S3). Similarly, no significant changes were observed in cyclic hyperoxia-exposed SD rats compared with room air-exposed controls (Supplemental Material Table S4). Adjusted p values were greater than 0.05. The entire Affymetrix microarray dataset can be found in Supplemental Material Table S5. The dataset has also been deposited in NCBI's Gene Expression Omnibus and is accessible through GEO Series accession number GSE18998 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE18998).

Fold change analysis of differentially regulated genes

In each of the six comparisons, the number of genes that showed a change of at least 1.5-fold accounted for 1.3 % or less of the total number of genes represented on the microarrays (Table 7). Because no hyperoxia-regulated genes reached statistical significance in the initial microarray analysis, an alternative approach was used. For each strain, genes that differed between cyclic hyperoxia and room air-treated groups were ranked by fold change. The top 50 upregulated and downregulated genes were submitted to DAVID to identify enriched gene groups (Table 8) [15, 16]. Genes from both strains were enriched for Gene Ontology (GO) terms related to carbohydrate metabolism. Many more of these terms appeared in the SD list compared to the F344 list. The SD list showed strong enrichment for genes involved in glycolysis, a process that is known to be affected by oxygen levels. In addition, the SD list returned the GO term “Response to Hypoxia” which was absent from the F344 output.

Table 7

Number of genes differentially expressed in the retina, based on fold changes irrespective of up- or downregulation of the gene

Comparison	Total number of genes	No. genes changed <1.5-fold	No. genes changed >1.5-fold	No. genes changed >2-fold	No. Genes changed >3-fold
F344 vs. SD (strain)	27,342	27,145 (99.3%)	141 (0.5%)	41 (0.1%)	15 (0.1 %)
O2 vs. RA (treatment)	27,342	27,309 (99.9%)	31 (0.1%)	2 (0.0%)	0 (0.0%)
F344 RA vs. SD RA	27,342	27,016 (98.8%)	269 (1.0%)	42 (0.2%)	15 (0.1%)
F344 O2 vs. SD O2	27,342	27,105 (99.1%)	168 (0.6%)	51 (0.2%)	18 (0.1%)
F344 O2 vs. F344 RA	27,342	27,293 (99.8%)	47 (0.2%)	2 (0.0%)	0 (0.0%)
SD O2 vs. SD RA	27,342	27,259 (99.7%)	70 (0.3%)	11 (0.0%)	2 (0.0%)

Percentages of genes changed are shown in brackets

Table 8

DAVID analysis for Gene Ontology terms enriched in the top 50 genes (ranked by fold change) that were up- or downregulated by cyclic hyperoxia in each strain

GO term	Number of genes	Fold enrichment a	False Discovery Rate b
F344
Cellular carbohydrate catabolic process	8	9.6	0.033
Carbohydrate catabolic process	8	9.1	0.049
SD
Glycolysis	8	18.6	3.72E-04
Glucose catabolic process	8	15.9	0.0011
Response to hypoxia	8	15.6	0.0013
Hexose catabolic process	8	15.2	0.0015
Monosaccharide catabolic process	8	15.2	0.0015
Alcohol catabolic process	8	14.7	0.0019
Cellular carbohydrate catabolic process	8	11.8	0.0083
Glucose metabolic process	9	11.4	0.0019
Carbohydrate catabolic process	8	11.1	0.0122
Monosaccharide metabolic process	11	11.0	8.62E-05
Hexose metabolic process	10	10.1	9.22E-04
Soluble fraction	9	7.3	0.0378
Alcohol metabolic process	12	7.1	0.0013
Cellular carbohydrate metabolic process	12	6.9	0.0018
Carbohydrate metabolic process	13	5.3	0.0072

GO terms are ranked by fold enrichment

aThe ratio between the frequency of term-annotated genes in the query gene list compared to the reference gene list, in this case the rat genome

bThe expected percentage of false positives in a given proportion of genes considered to be significant

Manual analysis of the gene lists showed that many genes known to be upregulated by hypoxia-inducible factor-1α (Hif-1α) were present (Table 9) and were almost exclusively downregulated by cyclic hyperoxia, the inverse response of these genes to hypoxia. The SD list contained almost twice as many Hif-1α-regulated genes as the F344 list. Those in the SD list were more strongly regulated by cyclic hyperoxia and represented 18 of the top 20 downregulated genes, compared with seven of the top 20 downregulated genes in the F344 list. qRT-PCR analysis of expression levels for Egln3, Bnip3, Slc16a3, and Hk2 (Fig. 2) using RNA pools identical to those used in the microarray analysis confirmed the strain-dependent difference in regulation of these genes by hyperoxia. Fold changes measured by qPCR for both strains were greater than those seen in the microarray results and all four candidate genes were downregulated in both strains in response to cyclic hyperoxia exposure. A substantially greater inhibition of gene expression was seen in SD rats compared to F344 rats.

Table 9

Known hypoxia-regulated genes identified from the top fifty genes (ranked by fold change) that were downregulated by cyclic hyperoxia in each strain compared to room air exposure

Gene	Gene description	Fold change
Genes downregulated by hyperoxia in F344 rats
Adra1b	Adrenergic, alpha-1B-, receptor	−2.12
Igfbp2	Insulin-like growth factor binding protein 2	−1.86
Slc16a3	Solute carrier family 16, member 3	−1.86
Cox4i2	Cytochrome c oxidase subunit IV isoform 2	−1.58
Genes downregulated by hyperoxia in SD rats
Slc16a3	Solute carrier family 16, member 3	−4.10
Igfbp2	Insulin-like growth factor binding protein 2	−3.42
Hk2	Hexokinase 2	−2.95
Pfkfb3	6-phosphofructo-2-kinase/fructose-2,6-biphosphatase 3	−2.80
Egln3	EGL nine homolog 3	−2.79
Bnip3	BCL2/adenovirus E1B 19 kDa-interacting protein 3	−2.63
Cox4i2	Cytochrome c oxidase subunit IV isoform 2	−2.45
Pdk1	Phosphoglycerate kinase 1	−2.31
Vegfa	Vascular endothelial growth factor A	−2.23
P4ha1	Procollagen-proline, 2-oxoglutarate 4-dioxygenase	−2.20
Mif	Macrophage migration inhibitory factor	−2.01
Ak3l1	Adenylate kinase 3-like 1	−1.97
Pfkl	Phosphofructokinase, liver	−1.93
Tpi1	Triosephosphate isomerase 1	−1.89
Egln1	EGL nine homolog 1	−1.72
Bhlhb2	Basic helix-loop-helix family, member e40	−1.67
Ldha	Lactate dehydrogenase A	−1.62
Pgk1	Phosphoglycerate kinase 1	−1.61
Aldoa	Aldolase A, fructose-bisphosphate	−1.57
Gpi	Glucose phosphate isomerase	−1.56

Only genes downregulated by 1.5-fold or greater are shown

Fig. 2

RT-PCR quantification of retinal gene expression in cyclic hyperoxia-treated and untreated F344 and SD rats at post-natal day 3. a Egln3; b Bnip3; c Slc16a3; d Hk2. Quantification was performed on RNA pools identical to those used in the microarray analysis. Expression levels were normalized to the reference genes Arbp and Hprt. Numbers above the bars are fold inhibition by hyperoxia. Error bars ± standard deviation. *p < 0.05

Discussion

Animal models of ROP [19] aim to recreate the biphasic progression of the disease seen in humans [27]. In our rat model, phenotypic changes representative of OIR are maximal among strains at post-natal day 14 [11], at which time SD rats show large retinal avascular areas indicative of susceptibility to OIR, whereas retinal avascular areas in F344 rats are small, indicating resistance to the disease [13]. At post-natal day 3, the retinae of room air- and cyclic hyperoxia-exposed F344 rats and room air-exposed SD rats showed normal morphology. However, cyclic hyperoxia-exposed SD rats showed aberrant vascular development, suggesting that the phenotypic changes reflecting differential susceptibility to OIR begin to manifest as early as post-natal day 3. We would, thus, expect the phenotypic changes to be driven by changes in gene expression.

Our rationale for examining gene changes at the early time point of day 3 was to increase the likelihood of identifying gene pathways related to the underlying genetic cause of phenotypic differences between strains. An Affymetrix microarray approach was used to screen for retinal mRNAs that are potentially involved in the development of OIR. The dataset presented herein identifies changes in gene expression as a result of strain, cyclic hyperoxia exposure, and a combination of the two. A limitation of this study is that rats were inevitably exposed to a brief period of normoxia (up to 20–30 min), following their removal from the hyperoxic environment. This may have resulted in alteration of the expression profiles of some genes; however, candidate genes were still able to be identified.

At the early time-point of P3, the majority of genes showing significant differential expression were as a result of strain alone. Of particular interest were genes which were differentially expressed as a result of strain and cyclic hyperoxia exposure because these differences may be causative of changes seen in the phenotype of rats with oxygen-induced retinopathy. Functional annotation of the top 50 genes (based on fold changes) that were upregulated or downregulated by cyclic hyperoxia in either strain identified a number of genes which are known to be regulated by Hif-1α. Differential expression of four genes: Egln3, Bnip3, Slc16a3, and Hk2 were selected for confirmation by qRT-PCR on RNA pools identical to those used in the Affymetrix microarray. Changes in gene expression presented herein are derived from qRT-PCR analysis unless otherwise stated.

Egln3 encodes PHD3, a member of the prolyl hydroxylase family of oxygen-dependent enzymes [28]. Hydroxylation of specific proline residues in HIF-1α by prolyl hydroxylases mediates the binding of the tumor suppressor protein von Hippel-Lindau to target HIF-1α for ubiquitinylation and proteasomal degradation [28, 29]. Prolyl hydroxylation is the main mechanism by which HIF-α levels are regulated in response to hypoxia [30]. There are three prolyl hydroxylases: PHD1, PHD2, and PHD3, all of which have been shown to have different functions in regulating HIF-1α and HIF-2α at the post-transcriptional level [31]. A study by Appelhoff et al. showed variation in the expression of prolyl hydroxylases in different human cell lines exposed to normoxia [31]. PHD2 was found to be the most abundantly expressed prolyl hydroxylase. Expression of PHD1 remained unchanged or significantly decreased upon exposure to hypoxia. PHD2 and, in particular, PHD3 were found to be greatly increased by hypoxia. PHD3 was shown to regulate HIF-2α to a greater extent than HIF-1α. The ability of PHD3 to remain active under hypoxic conditions suggests it may play a role in the physiological regulation of HIF, in particular HIF-2α. The activity of each prolyl hydroxylase appears to be dependent on the relative abundance of the other prolyl hydroxylases present under a specific set of conditions, adding to the complexity of HIF-α regulation [31].

The oxygen-dependent regulation of HIF-α levels by prolyl hydroxylases was identified as a pathway of interest. Egln3 contains a hypoxia response element, allowing its induction by hypoxia to occur through a HIF-dependent negative feedback loop [32]. Studies of Hif-2α knockdown mice with OIR show an absence of neovascularization, eventually resulting in the degeneration of the neural layers of the retina due to poor vascularization [33]. The importance of HIF-2α in retinal angiogenesis has been supported by another study involving mice with only one functional Hif-2α gene, in which retinal neovascularization was reduced in response to cyclic-hyperoxia exposure, due in part to an absence of pro-angiogenesis factors required to promote an angiogenic response [34]. At post-natal day 3, Egln3 was found to be decreased 5.1-fold in SD rats in response to cyclic hyperoxic exposure. Egln3 was decreased to a lesser extent (2.1-fold) in cyclic hyperoxia-exposed F344 rats compared to room air-exposed rats. In a recent study by Recchia and colleagues, both rat and mouse models of OIR were examined at early and late time points of disease progression [35]. At the early time points of P15 (rat) and P13 (mouse), an increase in Egln3 of 2.5- and 3.5-fold, respectively, was observed. At the later time points of P20 (rat) and P18 (mouse), Egln3 expression increased by 3.0- and 5.3-fold, respectively. These increases in Egln3 expression may be a result of the response to relative hypoxia, particularly at the later time point of OIR, rather than in response to hyperoxic exposure.

A decrease in Egln3, as observed in our study, might be expected to result in the stabilization of HIF-1α and increased expression of genes regulated by transcription factors, such as VEGF; however, Affymetrix data showed that VEGF was decreased 2.2-fold in SD rats in response to cyclic hyperoxia compared to room air exposure. Any stabilization of HIF-1α expression in the retina in response to hyperoxia exposure is difficult to determine because under normoxic conditions, HIF-1α is rapidly degraded by the ubiquitin proteasome pathway and has a half-life of less than 5 min [30, 36, 37], and downstream processing of retinae in normoxia following removal of neonates from a hyperoxic environment takes 20–30 min. The levels of the prolyl hydroxylases in response to hyperoxia exposure may potentially be able to be used as surrogate markers for HIF-1α in future studies.

Bnip3, expression of which was decreased 3.6-fold in cyclic hyperoxia-exposed SD rats at post-natal day 3, is a member of the Bcl-2 family of proteins that play a role in apoptosis [38]. To our knowledge, expression in the retina has not previously been reported. However, Hif-1α mediated overexpression of Bnip3 in response to hypoxia is known to exert a proapoptotic effect in human kidney epithelial cells [39]. Slc16a3 is involved in the transport of monocarboxylates such as pyruvate and lactate across cell membranes and is known to be expressed in ocular tissues including the retina [40, 41]. Slc16a3, downregulated 4.7-fold in oxygen-exposed F344 rats and 6.2-fold in oxygen-exposed SD rats, helps to export lactate, the end product of glycolysis, rapidly out of retinal cells and prevent intracellular acidosis [40]. It has been hypothesized that loss of Slc16a3 and other monocarboxylate transporters may compromise glucose metabolism and result in abnormalities in retinal function and retinal degeneration [42]. Hk2, downregulated 4.5-fold in oxygen-exposed SD rats, is involved in glycolysis and has previously been reported to be upregulated by hypoxia [43].

We have previously shown that the temporal expression of pro-angiogenic and anti-angiogenic factors contributes to susceptibility of different strains of rats to OIR [13]. A correlation between disease susceptibility and differential gene expression during early retinal vascular development as a consequence of cyclic hyperoxia exposure, rather than exposure to prolonged periods of hypoxia, was observed. Similarly in a murine model of OIR, Sato et al. showed a correlation between gene expression patterns and disease progression [44]. Analysis of gene expression at a late time point showed that genes associated with angiogenesis were upregulated in response to a sustained period of hyperoxia, coinciding with increased vascular development, consistent with our previous findings.

Little is known about the differences in gene expression at early time points in OIR. Natoli and colleagues examined the temporal gene expression patterns in mice exposed to sustained periods of hyperoxia [45]. They found that after 3 days of exposure to hyperoxia, the early response gene Oxr1 was upregulated in C57BL/6J mice susceptible to hyperoxia whereas no regulation of Oxr1 was observed in hyperoxia-resistant BALB/cJ mice. In our studies, no significant changes in Oxr1 gene expression were observed as a result of strain or cyclic hyperoxia exposure at post-natal day 3.

Using Affymetrix arrays as a screening tool to identify which genes, if any, were likely to be involved in susceptibility to OIR, we observed that most genes were not significantly differentially expressed, following correction for multiple comparisons. An alternative method for identifying genes of interest was thus required after the initial microarray screen. We chose to use fold changes as an indication of differential gene expression rather than adjusted p values and arbitrarily chose the top 50 genes found to be upregulated or downregulated in each strain following cyclic hyperoxia exposure. Most of these genes showed a fold change of 1.5-fold or more. A bioinformatics approach was then taken to examine the biological contexts of these changes in gene expression.

Analysis of gene lists using DAVID highlighted a strong enrichment for known hypoxia-responsive genes in SD rats, which was much less evident in F344 rats. This suggests that differences in one or more of these genes may contribute to the strain-specific response to cyclic hyperoxia. Almost all of the hypoxia-responsive genes identified were downregulated by cyclic hyperoxia, the inverse of their response to hypoxia. Sears et al. have shown that the vaso-obliterative effects of cyclic hyperoxia on retinal vessels can be negated by inducing a pseudo-hypoxic environment through inhibition of prolyl hydroxylases, resulting in stabilization of HIF-1α [46]. An attenuated response to cyclic hyperoxia in F344 rats may be responsible for a similar outcome. Another possibility is that alterations in oxygen sensing pathways in SD rats exposed to cyclic hyperoxia may be contributing to susceptibility to OIR. The potential role of differential ischemic preconditioning in different strains is unknown. Examining the temporal expression of the candidate genes by qRT-PCR or a proteomics approach may give a greater insight into the kinetics of OIR and identify gene expression patterns which may contribute to susceptibility to OIR in rats.

Top 251 genes regulated by Strain (XLS 67 kb)

Top 50 up- and downregulated genes for each strain in cyclic hyperoxia vs. room air comparison. Genes were ranked by fold change only (XLS 22 kb)

Top 20 genes regulated by cyclic hyperoxia exposure in F344 rats (XLS 18 kb)

Top 20 genes regulated by cyclic hyperoxia exposure in SD rats (XLS 18 kb)

Entire Affymetrix dataset (XLS 37 mb)

Summary of experimental design showing how RNA samples were pooled for microarray analysis. a RA-exposed F344 rats collected from two separate litters to create two samples; b cyclic hyperoxia-exposed F344 rats collected from three separate litters to create three samples; c RA-exposed SD rats collected from three separate litters to create three samples; d cyclic hyperoxia-exposed SD rats collected from two separate litters to create two samples (GIF 109 kb)

High resolution image file (TIFF 286 kb)

Principal component analysis of samples (n = 2). a F344 rats raised in room air, b F344 rats exposed to cyclic hyperoxia from birth, c SD rats raised in room air, d SD rats exposed to cyclic hyperoxia from birth. Samples clustered according to strain and exposure with cyclic hyperoxia. A large variation was seen in the second principal component in b. Variation was also seen in the first principal component in c. As a result, a third pool of RNA was submitted for analysis for these two experimental conditions (GIF 219 kb)

High resolution image file (TIFF 276 kb)

Principal component analysis of original and additional samples corrected for batch variation as a result of two separate sample runs. a F344 rats raised in room air (n = 2), b F344 rats exposed to cyclic hyperoxia from birth (n = 3), c SD rats raised in room air (n = 3), d SD rats exposed to cyclic hyperoxia from birth (n = 2). Samples clustered according to strain but not by exposure to cyclic hyperoxia with the additional samples. Reduced variation in principal component two was seen in b. The variation in principal component one for c also decreased. Additional samples were considered to group in an acceptable manner to be included in the overall analysis of mRNA expression (GIF 232 kb)

High resolution image file (TIFF 64 kb)