Transcriptome Analysis of the Preterm Rabbit Lung after Seven Days of Hyperoxic Exposure

The neonatal management of preterm born infants often results in damage to the developing lung and subsequent morbidity, referred to as bronchopulmonary dysplasia (BPD). Animal models may help in understanding the molecular processes involved in this condition and define therapeutic targets. Our goal was to identify molecular pathways using the earlier described preterm rabbit model of hyperoxia induced lung-injury. Transcriptome analysis by mRNA-sequencing was performed on lungs from preterm rabbit pups born at day 28 of gestation (term: 31 days) and kept in hyperoxia (95% O2) for 7 days. Controls were preterm pups kept in normoxia. Transcriptomic data were analyzed using Array Studio and Ingenuity Pathway Analysis (IPA), in order to identify the central molecules responsible for the observed transcriptional changes. We detected 2217 significantly dysregulated transcripts following hyperoxia, of which 90% could be identified. Major pathophysiological dysregulations were found in inflammation, lung development, vascular development and reactive oxygen species (ROS) metabolism. To conclude, amongst the many dysregulated transcripts, major changes were found in the inflammatory, oxidative stress and lung developmental pathways. This information may be used for the generation of new treatment hypotheses for hyperoxia-induced lung injury and BPD.


Introduction
Preterm birth leads to a dysregulated development in many organs that are not yet adapted to postnatal life. In the pulmonary system, this often results in bronchopulmonary dysplasia (BPD), a complex disease in which multiple factors interact. Premature lungs (most frequently in the saccular stage of lung development) are exposed to hyperoxic and hyperbaric conditions during ventilation and administration of supplementary oxygen. This process is often amplified by pre-or postnatal infections, fluid imbalance, malnutrition, genetic predispositions, etc. Persistent inflammation overwhelms natural tissue repair and leads to an arrest in alveolar development and vasculogenesis. [1] The resulting lung parenchyma is composed of rudimentary alveoli, interstitial thickening and a dysmorphic capillary configuration. [2] These morphological changes can be categorized as a developmental arrest.
The aforementioned risk factors are well-known and the application of less invasive ventilation strategies and permissive hypoxemia have proven their efficacy. They became the cornerstones of current neonatal management. [3,4] Despite this, BPD continues to be a frequent complication of premature birth. About 15-25% of very low birth weight (VLBW) infants develop BPD, and rates in extremely low birth weight (ELBW) infants are even higher. [5] Moreover, BPD remains an important risk factor for lung disease later in life. [6] At the molecular level BPD is still poorly studied compared to other pathologies. Several individual molecules assessed in bronchoalveolar fluid have been proposed to play a key role in BPD like interleukin-8 or matrix metalloproteinase-3. [7,8] Most studies however where hypothesis driven or they examined certain pathways in lung disease and repair.
Animal models are needed to study disease mechanisms and to evaluate new preventive or therapeutic strategies for BPD. Most research into BPD has been performed in hyperoxiaexposed rodent models. [9][10][11] Unfortunately, their lung development differs from humans as birth occurs in the early saccular stage of lung development. In rodents, alveolization only starts several days postnatally while humans start alveolizing in utero. This entails that not all the relevant findings can be extrapolated to the human context. As such, there is an advantage to study the hyperoxia induced lung injury in animal models that mimic human development more closely. The rabbit is considered a large animal model that has favorable characteristics for this research question. In contrast to rodents, rabbits indeed start alveolizing prior to birth, as do pigs, sheep, primates and humans. [12] Moreover, rabbits are easy to handle and house, have a large litter size and allow technical manipulation of the fetus at relevant developmental stages which makes them ideal for the study of effects of perinatal interventions. Therefore we used the rabbit as a model for the study of hyperoxia-induced lung injury in the preterm born pup. [13] Rather than focus on single putative molecule or pathways, we herein used a more complete 'systems biology'-approach, [14] that allows the exploration of larger patterns and networks. Herein we analyze transcriptome data using software (IPA) that combines expression data with current generic knowledge of molecular interactions. Our goal was to obtain more complete insights in the pathophysiology and generate new theoretical therapeutic strategies for preterm hyperoxia induced lung injury.

Animal model
We earlier described in detail the preterm rabbit model for hyperoxia-induced lung injury. [13] Briefly, time mated pregnant does underwent cesarean section at 28 days (term = 31 days) of gestation (early saccular lung developmental phase). The pups were randomly divided into two groups: (1) the normoxia group, where pups were housed in 21% oxygen and (2) the hyperoxia group, where pups were nursed in hyperoxia ( 95% oxygen); both for seven days. The model has been earlier described in detail elsewhere [15][16][17]. Briefly, immediately after delivery, pups were placed in an incubator at 32°C, fed twice daily via an orogastric tube and received prophylactic antibiotics and vitamin K. All animals were treated according to current guidelines of animal well-being, and the experiments were approved by the Ethics Committee for Animal Experimentation of the Faculty of Medicine of the Katholieke Universiteit Leuven (Leuven, Belgium).

Harvesting of specimens
At day 7 of life, pups were harvested for histological and transcriptome analysis. After euthanizing the pups with a mixture of embutramide 200mg, mebezonium 50mg and tetracain hydrochloride 5mg (intracardiac injection of 0.1mL T61, Intervet Belgium NV, Mechelen, Belgium), thoracotomy was performed and both lungs and trachea were removed "en bloc". The left bronchus was ligated and the left lung was snap-frozen, while the right lung was processed for histological analysis (see below). Six pups per group were used for histological analysis. Out of the snap frozen samples, four left lungs from the normoxia group, and four from the hyperoxia group were randomly selected for transcriptome analysis.

Lung injury score
A 20G catheter was inserted in the trachea where after the right lung was fixed with 4% paraformaldehyde by immersion and perfusion under a constant hydrostatic pressure of 25cm H 2 O for 24 hours before embedding. Paraffin sections were stained with hematoxylin and eosin (HE) for further histological analysis. Lung injury score was assessed on 20 random highpower fields (400x total magnification) for each pup. The selection of random fields was obtained by successive random displacements (each at least one high power field in length) from the initial position provided that at least 50% of each field was occupied by lung alveoli. A three-tiered scheme was used to quantify each of five histological parameters: presence of neutrophils in alveolar and interstitial space, hyaline membranes, debris filling the airspaces and septal thickening. [18] The data were analyzed statistically using student-T test, significance was set with a p-value < 0.05.

RNA isolation and sequencing
RNA isolation was performed on snap frozen lung tissue using the RNeasy mini kit (Qiagen). Tissue lysis and homogenization was performed in 1200μL Buffer RLT using the TissueLyser system (Qiagen). Following tissue disruption and homogenization, samples were centrifuged for 3 min at 14000rpm in a benchtop micro centrifuge. Lysate was transferred to fresh tubes and an equal volume of 70% ethanol was added. 600μL of sample was added twice to a spin column, with 2 RNeasy spin columns used per sample. Following wash steps RNA was eluted in 50μL RNAse-free H2O. Total RNA quantification was performed using the Nanodrop 1000 spectrophotometer (Thermo Scientific). RNA integrity was assessed using the RNA 6000 Nano Kit and the Bioanalyser (Agilent Technologies) according to the manufacturer's recommendations. Between 1 and 2μg of total RNA was used as input material for sequencing library preparation which was performed with the TruSeq RNA Library Preparation Kit (Illumina) according to the manufacturers protocol. Fragmentation was performed for 6 minutes. 8 PCR cycles were used for PCR enrichment step. Samples were indexed to allow for multiplexing. Sequencing libraries were quantified using the Qubit fluorometer (Life Technologies). Library quality and size range was assessed using the Bioanalyser (Agilent Technologies) with the DNA 1000 Kit (Agilent Technologies) according to the manufacturer's recommendations. Individual libraries were diluted to a final concentration of 2nM and pooled for sequencing. Pooled libraries were then sequenced in a single lane of an Illumina HiSeq2000 flow cell generating single end 50bp reads.

Data analysis
Expression values were calculated per gene and normalized to reads per kilobase per million reads (RPKM) values as described by Mortazavi. [19] An RPKM value of >1 was applied as cut-off to distinguish between 'noisy' or 'leaky' transcription and 'true' mRNA expression levels for comparison between samples. For comparison between both groups, fold changes (FC) were calculated, indicating the ratio of expression. We applied a cut-off of >2 or <-2, in order to decrease the number of transcripts in our analysis. False discovery rates (FDR) were calculated as a measure for statistical significance of the difference in expression of a gene between the two groups. We considered a transcript change significant if FDR was <0.01.
Further analysis was performed using the Ingenuity Pathway Analysis (IPA) software. The IPA 'Upstream Regulator Analysis' predicts upstream regulators by combining the directional expression changes from our mRNA-sequencing, and knowledge from prior experimental reports on causal effects between molecules (endogenous and exogenous), compiled in the IPA Knowledge Base. Upstream Regulator Analysis calculates a z-score based on the edge of dysregulation of all the downstream molecules and the uniformity of the existing evidence about the upstream-downstream relation, for every upstream regulator known to have a causal effect on at least 4 dysregulated transcripts. Z-scores of <-2 and >2 respectively indicate a significant inhibition and activation state of the upstream regulator, regardless of the actual expression level of these molecules. [20] The Network Generation Algorithm links molecules based on experimentally observed interactions, and orders these molecule based on their interconnectedness. In general, the more interactions with other network members, the more central a molecule will be in a network.

Lung injury score
Histological assessment of both normoxia as well as hyperoxia exposed lungs reveals a clear increase in neutrophil infiltration and alveolar wall thickening (Fig 1A and 1B). [13] Generating a combined lung injury score revealed a significant difference between both groups (Fig 1C) with an increase in lung injury observed in animals held in hyperoxia (p<0.01).

Hyperoxia and gene expression profiles
We identified 2217 transcripts that were dysregulated (FDR < 0.01) in hyperoxia by at least a factor 2 (FC>2 or FC<-2). A heat map of these genes visualizes a clear distinction between hyperoxia and normoxia. (Fig 2). We were able to link IPA-gene names to 1989 of the 2217 dysregulated transcripts. This group was used for further analysis in IPA. 1023 of these transcripts are upregulated in hyperoxia, with fold changes up to +235.786. The other 966 are downregulated, as far down as -38.476. The ten most upregulated molecules in hyperoxia are shown in Table 1, the ten most downregulated in Table 2. A complete list of these 1989 dysregulated molecules can be found in S1 Table. Pathway analysis of differentially expressed genes Running the Upstream Regulator Analysis, we identified 58 inhibited and 160 activated significant upstream regulators. Both endogenous and exogenous molecules are included in this analysis. All significant upstream regulators which are classified as transcription factors are displayed in Table 3. All other upstream regulator types can be found in S2 Table. Significantly dysregulated molecules were sorted by their function or involvement in pathophysiological processes. An overview of the most relevant dysregulated molecules is provided in Table 4.

Discussion
Under the assumption that a rabbit model has several advantages over rodent models for studying BPD with regard to lung development and the potential of experimental therapeutic interventions, we used this preterm rabbit model [13] to comprehensively describe possible pathways involved in the pathophysiology of BPD using modern techniques. Herein, we used whole transcriptome analysis as this offers an comprehensive amount of expression information, that allows to document the intricate and complex processes that occur in the damaged developing lung. Lung injury score demonstrated results that parallel lung damage seen in clinical BPD-samples. [2] This confirms the relevance of insights in the molecular mechanisms of BPD obtained in our model. We focused on four relevant biological networks in this context: inflammation, lung development, vasculogenesis and ROS metabolism. Several inflammatory mediators were significantly dysregulated in our model, confirming an important role for inflammation in the pathogenesis of hyperoxia-induced lung injury. The pattern of dysregulation of these molecules is consistent with a pro-inflammatory state. This is in line with data from other animal models [21] and clinical studies. [22] The increased TNFαactivity induces many of the transcriptional changes seen in our model, some of these downstream mediators have been described in hyperoxia exposed term rabbit lungs. [23] It remains unclear if this TNFα-activity is increased as a protective pathway or merely secondary to the inflammation in the preterm lung. Experimental therapy both with TNFα and anti-TNFα Table 4. Dysregulated genes of special interest. Dysregulated genes of special interest for hyperoxia-induced lung injury; in order of appearance; FC: fold change, FDR: false discovery rate.

IPA ID
Full Name FC FDR show an improvement in rats exposed to hyperoxia. [24,25] To date this has not been corroborated in larger animal models.
Many preclinical studies have reported on molecules involved in lung development [26][27][28][29] that can account for the developmental arrest in BPD. Several of these lung development molecules are also dysregulated in our model. (Fig 4) PPARγ, which was also an upregulated molecule in our study, is important for normal lung maturation. [30] Furthermore, beneficial effects  have been achieved with stimulation of PPARγ, e.g. by rosiglitazone in the rat model of hyperoxia-induced lung injury. [31][32][33] A second upregulated molecule is DKK1, which is an inhibitor of WNT-signaling. Its known involvement in branching morphogenesis and alveolization could account for a role in the developmental arrest of BPD. [28,34] Other dysregulated factors of pulmonary development include SPP1 (downregulated during secondary septation [27]), CAV1(involved in alveolar septation [35] and acute lung injury [36]) and ACE (involved in secondary septation [37]).
We also suggest a role for KLF2, which we identified as an important upstream regulator. It is highly expressed in adult mouse lung tissue and was demonstrated to be important for late stage lung development (saccular and alveolar phase). [38] Furthermore, KLF2 is involved in endothelial homeostasis and angiogenesis [39] and has negative regulatory effects on immune cell activation. [40] Normal lung development depends, certainly in the later stages, on vasculogenesis and the development of mesenchymal structures: blood vessels and respiratory epithelium are in constant interaction. Several molecular factors responsible for angiogenesis in the lung [26] are changed in our hyperoxia model. VEGFA seems a central molecule in most of our generated networks. The administration of VEGFA to attenuate lung damage has been tried in different models with varying results. [41,42] Furthermore knockdown of VEGF abolishes the protective effects of mesenchymal stem cells in a rat model of hyperoxic lung injury thereby demonstrating the critical role of VEGF. [43] The upregulation of several enzymes involved in scavenging and detoxifying of ROS can be considered as protective against the potential adverse effects of hyperoxia. NFE2L2 plays a critical role in the induction of this response. [44] Another important player we identified, is CYP1A1. Induction of CYP1A1 via the aryl hydrocarbon receptor (AHR) pathway, using βnaphthoflavone or omeprazole attenuates hyperoxic lung injury, which suggests a protective effect of CYP1A1. [45,46] We acknowledge certain limitations of our study. First, in this model we use only hyperoxia to create lung injury. This is an imperfect approximation for human 'new' BPD. [1,47] In human BPD, several other factors (barotrauma, infections, malnutrition, fluid imbalance, antenatal corticosteroids, surfactant, . . .) are involved, yet are not included in this experimental model. Secondly, there are some methodological limitations. mRNA sequencing was performed on whole, mixed, lung tissue which offers an overview of expression changes in all pulmonary cell types. This results in loss of 'spatial resolution' as specific up regulation in individual cell types constituting the airways, vasculature and parenchyma cannot be quantified in this manner. Furthermore, our samples were only examined at a single time point (7 days of hyperoxic exposure), a loss of temporal resolution has to be considered as well. Although in our previous study no functional differences were found between term and preterm animals held in normoxia [13], no transcriptome analysis has been performed on term animals. Some of the expressional differences might be explained by prematurity alone without a clear hyperoxic cause. A final limitation is the use of selective filtering of the obtained results. First, filtering for small RPKM's and fold changes may result in the loss of meaningful subtle dysregulated master switch molecules. The used software, IPA (Upstream Regulator Analysis and Network Generation Algorithm) is also dependent on prior knowledge and many expression-regulating relationships remain still to be discovered or described. This may bias our analysis, in favor of molecules which have been frequently investigated.
Despite these theoretical limitations, other results of our experiment show obvious similarities with previous studies analyzing specific molecules which may be involved in the pathophysiology of BPD. [27,48,49] The first two studies [27,48] examined the changes in expression profile during the last stages of lung development in rats without any hyperoxic lung damage. Important genes which were upregulated during the saccular stage were those involved in development. Both studies found the Wnt signaling pathway to be affected by different genes like Fzd1 (frizzled-1) or Ptn (pleiotrophin). Our data demonstrated DKK1, an inhibitor of Wnt signaling to be upregulated after exposure to hyperoxia. In the study of Bhattacharya [49] whole lung tissue was used from mice exposed to 100% hyperoxia for 10 days. They found Ahr to be a key regulator in the pathophysiology of hyperoxic lung damage. This receptor is known to activate CYP1A1, one of the main dysregulated molecules in our database. Furthermore, Cdkn1a (cyclin-dependent kinase inhibitor 1a), which is a cell cycle regulator, is upregulated both in mice as in our rabbit lung tissue. The exact mechanism how this gene is involved in the pathophysiology of hyperoxic lung injury remains to be determined.
To conclude, our study is the first to perform whole transcriptome analysis on hyperoxiainduced, preterm rabbit lung tissue identifying several central molecular mediators. The major pathophysiological features of BPD we identified are inflammation, lung developmental arrest, dysregulated vasculogenesis and increased ROS metabolism. We will use the molecular master switches revealed in this analysis for the future definition of novel preventive strategies to avoid or reduce the occurence of BPD.
Supporting Information S1