Pulmonary transcriptome analysis in the surgically induced rabbit model of diaphragmatic hernia treated with fetal tracheal occlusion

The prevalence of congenital diaphragmatic hernia (CDH) ranges between 1-4/10,000 births, which translates to 542 to 2168 children in EU-27 member states every year (2008) (Kotecha et al., 2012). It can occur either as an isolated condition, associated with other anomalies or as part of a genetic syndrome. In CDH, lung development is disturbed already from the embryonic period, but progresses because of herniation of the viscera through the defect into the chest, competing for space with the developing lungs. As a consequence, lungs of babies with CDH have fewer and less mature airway branches, a smaller cross-sectional area of pulmonary vessels, remodeled vascular architecture, and an altered vasoreactivity (Kinsella et al., 2005). At birth, this leads to respiratory insufficiency and persistent pulmonary hypertension (PPHT). This is lethal in up to 30% of cases, despite prenatal referral to high-volume centers that offer standardized neonatal care (Grushka et al., 2009; Hayakawa et al., 2013; van den Hout et al., 2011).

In countries with prenatal ultrasound (US) screening practices, >60% of cases are diagnosed at the latest by the second trimester. In isolated cases, an individualized prognosis can be made based on the degree of liver herniation and parenchymal lung size, predicting mortality and early neonatal morbidity (Cannie et al., 2008; Claus et al., 2011; Cruz-Martinez et al., 2011; Jani et al., 2007,, 2009; Mayer et al., 2011; Ruano et al., 2012). The ability to prenatally identify a future non-survivor enables clinicians to propose a prenatal intervention that avoids this outcome. Today, the clinical method to stimulate lung development is fetal tracheal occlusion (TO). TO prevents egress of lung liquid, causing increased pulmonary stretch, which accelerates lung growth (DiFiore et al., 1994). This operation is currently evaluated within the framework of a randomized clinical trial, comparing fetoscopic endoluminal tracheal occlusion (FETO) to expectant management during pregnancy (Deprest et al., 2009,, 2014).

Because TO still fails to save half of the fetuses with severe CDH and carries a risk for preterm delivery, there is a need for alternative prenatal strategies. Preclinical evaluation of these therapies involves the use of animal models, i.e. rodents, rabbits, and eventually lambs. Although, in all mammalians, lung development progresses through five stages (embryonic, pseudoglandular, canalicular, saccular, alveolar), this is not at the same pace (Pringle, 1986). This needs to be kept in mind when choosing an animal model in studies on lung development. For instance, rodents are ‘early’ models because CDH is induced in the embryologic phase by the administration of a teratogen (Clugston et al., 2006; Keijzer et al., 2000). In these models, nitrofen (NF) is prenatally administered at embryonic day (E)11.5 (mice) or E13.5 (rats), and effects are studied at term, when lungs are in the canalicular to saccular phase (Beurskens et al., 2007; Correia-Pinto et al., 2010). All NF-exposed pups have lung hypoplasia, half of them also the diaphragmatic defect (Keijzer et al., 2000; Rottier and Tibboel, 2005; Schittny et al., 2000). These effects are believed to be caused by interference of NF with the retinoic acid (RA) signaling pathway, which is key in lung development (Malpel et al., 2000). This pathway is also disturbed in CDH in humans (Beurskens et al., 2009).

In larger animals, hypoplasia is induced by surgical creation of a diaphragmatic defect, typically in the pseudoglandular phase. Rabbits are widely available, have low housing needs, a timed gestation and large litter size. Rabbits alveolize in utero as in humans (Roubliova et al., 2010). Pups with diaphragmatic hernia (DH) display both histological and functional changes, such as reduced airway and vascular development, and pathologic compliance, airway resistance, tissue damping and elastance – mimicking the clinical phenotype (Flemmer et al., 2007; Roubliova et al., 2004; Wu et al., 2000). Gene expression of a number of critical signaling molecules relevant to alveolarization, angiogenesis and regulation of vascular tone, but not to surfactant production, have been shown to be disrupted just as in humans (Vuckovic et al., 2013,, 2012). However, a broader study on gene expression levels in this model has not been done so far. The use of RNA-sequencing (RNA-seq) for transcriptome analysis has become increasingly widespread with the advent of massively parallel sequencing technologies, in part owing to reductions in costs and increased throughput, and improved knowledge of non-model-organism reference genomes. Therefore, we wanted to investigate the pulmonary transcriptome after surgically induced DH creation and subsequent TO in the rabbit model. The provided gene expression database can be used to develop further treatment strategies for CDH.

At harvest, there were seven surviving DH+TO fetuses [mean lung-to-body weight ratio (LBWR): 0.017; standard deviation (s.d.): 0.002; confidence interval (CI) 95%: 0.013-0.022) and seven DH fetuses (mean LBWR: 0.011; s.d.: 0.003; CI 95%: 0.003-0.018). We also took one random control for every third litter (n=4) (mean LBWR: 0.016; s.d.: 0.001; CI 95%: 0.013-0.018). RNA integrity (RIN) values were determined, which led to the exclusion of one DH+TO fetus. All six remaining fetuses from the DH+TO group were analyzed. We selected four (DH) fetuses displaying the best RIN values (data not shown). Principal component analysis (PCA) and unsupervised hierarchical clustering revealed that three DH+TO samples clustered together as a single homogenous group, and two DH+TO samples clustered with the DH group. Those two samples correlate with a clinical sub-group of non-responders to the TO treatment, analysis of which was not the aim of our study. The remaining DH+TO sample was considered as a statistical outlier based upon quality control (QC) correlation analysis and calculation of modified allele dosage (MAD) scores using the Array Studio software (OmicSoft, Cary, USA); hence, it was excluded from further analysis. This distinct gene expression pattern of DH+TO fetuses correlated with two patterns of lung growth, as evidenced by the LBWR, which is a pathologic measure of the degree of lung development (Table 1).

Statistical inference analysis was performed between the groups defined above [DH (n=4), DH+TO (n=3) and control (n=4)] using the general linear model function in Array Studio to calculate fold change (FC) values and false-discovery rate (FDR) values (applying the Benjamini–Hochberg procedure).

PCA and unsupervised hierarchical clustering of those dysregulated genes in any group comparison (described above) demonstrated clustering of homogenous groups for DH, DH+TO and control as expected, and revealed that DH+TO and control groups were more comparable than the DH group (Figs S1, S2).

The heat map generated by the unsupervised hierarchical clustering for the total of 641 unique genes found to be dysregulated (FC±2; FDR<0.1) in any group comparison is shown in Fig. 1. This reveals three large gene clusters whose expression values define the changes due to DH and after treatment with TO. The first gene cluster, containing 157 genes (Fig. 1, top) demonstrates low expression in both DH and control, and high expression in DH+TO only. This cluster represents genes whose expression is induced by TO, but which are unchanged in DH. The second, and largest, gene cluster, containing 378 genes (Fig. 1, mid) demonstrates high expression in DH only and low expression in control and DH+TO. This cluster represents genes whose expression is induced by DH, and which are reversed to lower levels by subsequent TO. The third gene cluster, containing 106 genes (Fig. 1, bot), demonstrates high expression in controls, and low expression in both DH and DH+TO. This cluster represents genes whose expression is reduced following DH, and which remain expressed at low levels even after TO.

To summarize, from the hierarchical cluster heat map shown in Fig. 1 it is apparent that the major change in gene expression (Fig. 1, gene cluster 2, mid) is upregulation of genes in DH, which are downregulated by TO towards, or beyond, levels seen in the controls (normal state) at the same developmental time point. One subset of genes is unchanged by TO and remains expressed at low levels, and one subset is highly upregulated by TO in comparison to both control and DH without TO.

The database of all dysregulated genes found in our study can be found in Table S1.

Analysis of the dysregulated gene sets for each comparison detailed above and for each of the gene clusters defined above was performed using the IPA web application (Ingenuity Pathway Analysis, Ingenuity Systems, Redwood City, CA, USA). Of the total 641 unique genes, 560 genes were mapped within IPA to known functional genes with human, mouse or rat homologs. This gene set was used to create a protein-protein interaction (PPI) network where experimental evidence supports direct interactions (i.e. predictions and indirect interactions were not allowed), and orphan molecules were subsequently removed. The resulting network is displayed in Fig. S3, where the molecules are colored according to the gene clusters identified in Fig. 1 (blue=genes with significantly increased expression in DH compared to DH+TO and/or control; yellow=genes with significantly increased expression in DH+TO compared to DH and/or control; and green=genes with significantly decreased expression in DH and/or DH+TO in comparison to control).

The Cytoscape plugin iRegulon (Janky et al., 2014) was used to identify potential transcription factors (TFs) of interest, using the 641 genes dysregulated in any comparison as the input. iRegulon searches for enrichment of cis-regulatory TF motifs in the target gene set, and predicted FOXJ1 as a key regulatory TF in our study. The input list contains 139/641 (22%) genes that are ‘target’ genes of FOXJ1.

In addition, we used the STEM application (Ernst and Bar-Joseph, 2006) to identify clusters of genes with similar changes in expression over time or, in our case, the changes observed from control to disease (DH) to treatment (DH+TO). The input was the FC values for the 641 genes significantly dysregulated in any comparison, where time point 1 is control data (i.e. FC=1), time point 2 is the CDH vs control data, and time point 3 is the TO vs control data. This identified four significant clusters containing genes with similar changes in expression profile, of which Profile 9 contained the largest number of genes and was of highest significance (Fig. 2, and Table S2).

PCR analysis of ten of the 12 analyzed genes yielded results that were in agreement with RNA-seq data (Fig. S4). For eight genes, the differential expression between groups was concordant with the RNA-seq data (P<0.05, Mann–Whitney U test). However, qPCR did not show a significant increase of RFX3 and LRRQI1 in the DH group compared to control (Fig. S5). For BMPR2 and PDE5A, qPCR failed to show any difference between DH and wild type. This distinct gene expression pattern of DH+T fetuses correlated with two patterns of lung growth, as evidenced by the LBWR, which is a pathologic measure of the degree of lung development. However, PDE5A was significantly downregulated in the TO group.

In this study, we describe for the first time the pulmonary transcriptome analysis of specimens obtained in a rabbit model for pulmonary hypoplasia. The latter was induced by creating a diaphragmatic defect during the pseudoglandular phase. Conversely, forced lung growth was induced by fetal TO at the transition of the canalicular to saccular phase. We found that the largest group of genes that were significantly dysregulated were 378 genes that were both upregulated by DH creation and downregulated by TO to a level similar to that of controls. Furthermore, this study gives a database of genes that are significantly influenced by DH creation and consecutive TO (Table S1). This database can be used for further understanding of the disease process and development of treatment modalities for CDH. Below, we discuss some of the most relevant genes that we found were dysregulated.

Many studies have documented expression changes for numerous genes in association with CDH and, to a lesser extent, also the effects of TO, all of this in various animal models of CDH. This is typically done by using PCR for selected genes, or using broader arrays, at least for experiments done in (NF-exposed) rodents, a species in which molecular tools are abundantly available. Using a more modern technique such as RNA-seq, one can now also document and screen for changes in gene expression in relevant animal models for pulmonary hypoplasia and induced lung growth, even if the genome has not been completely identified. We herein used this technology in rabbits, and studied dysregulations associated with pulmonary hypoplasia and TO. The latter is done to force lung growth, which has been abundantly demonstrated in several animal models of CDH using other tools (Khan et al., 2007). In this experiment, we confirmed the effects of TO on cell proliferation. MKI67 is a key cell proliferation marker and was significantly upregulated by TO. The effect of TO was in proportion to the pre-existing lung size (FC 3.6056, FDR 0.0874; TO vs control), a phenomenon already clinically demonstrated (Peralta et al., 2008). Other genes related to the mechanisms of TO in previous studies were identified. To name only one, AQP4 was upregulated in DH (FC 2.3562, FDR 0.018) yet went down to a normal level after TO (FC −2.7286, FDR 0.0106). AQP4 is a transport molecule that has previously been shown to be downregulated by stretch force on fetal lung AT2 cells in vitro (Wang et al., 2006).

TO is, however, in essence performed to enhance alveologenesis and reverse vascular changes that are typical of this disease. We observed dysregulation of connective tissue growth factor (CTGF). CTGF indeed enhances alveologenesis and microvascularization during late lung development (Hall-Glenn and Lyons, 2011). CTGF was earlier shown to be upregulated by TO in the saccular phase of E21 rats (Burgos et al., 2010; Mesas-Burgos et al., 2009). In line with that, CTGF was, in our rabbits, also firmly upregulated by TO (FC 2.8135, FDR 0.0151). The same dysregulation was observed for a number of other genes that have been previously associated with DH and TO. BMPR2 expression was increased in DH (FC 1.9265, FDR 0.0228) and went down after TO (FC −2.7956, FDR 0.0097). This gene encodes a member of the bone morphogenetic protein (BMP) receptor family of transmembrane serine/threonine kinases and the ligands of this receptor are the BMPs, which are members of the TGF-β superfamily and are involved in embryogenesis. Disruption of BMPR2 activity and downstream signaling has been demonstrated in the NF rodent model of CDH (Gosemann et al., 2013; Makanga et al., 2013). BMPR2 plays a key role in pulmonary vasculogenesis during the late stages of fetal lung development and is essential for control of endothelial and smooth-muscle cell proliferation. Dysfunction of BMPR2 and downstream signaling have been shown to disturb the crucial balance of proliferation of smooth muscle cells, contributing to the pathogenesis of human and experimental pulmonary hypoplasia (Gosemann et al., 2013; Makanga et al., 2013). In line with this was downregulation of BMP7 in DH, yet significant upregulation by TO (FC 2.311, FDR 0.0163). BMP7 was reportedly downregulated in the lungs of the NF rodent (Makanga et al., 2013).

We were, however, mainly interested in the RNA-seq technology because it has the potential to identify new targets that are involved in the condition and that can respond to prenatal intervention. One example is that DH in our rabbits was associated with a significant increase in PDE5A gene expression. The latter has been well documented in rodents and sheep, and it has been tied to the postnatal problem of pulmonary hypertension (Noori et al., 2007). The interesting observation here is that TO dramatically reduced PDE5A gene expression (FC −2.49; FDR 0.0097). PDE5A gene expression might also be altered by an alternative prenatal intervention, i.e. the transplacental administration of drugs. Sildenafil, which is already used to treat PPHT of the newborn (Noori et al., 2007), was recently shown to effectively inhibit cyclic guanosine monophosphate (cGMP)-specific phosphodiesterase type 5 in utero (Luong et al., 2011). In the latter experiment, NF rodents were used, and we confirmed a similar effect in rabbits, and others in the fetal lamb (Luong et al., 2011; Russo et al., 2015; Shue et al., 2014).

When starting from our RNA-seq data, the Cytoscape plugin iRegulon predicted that over 20% of dysregulated genes were target genes of FOXJ1. FOXJ1 is a marker for ciliated cells; it is actually involved as a transcription factor in ciliogenesis (Yu et al., 2008). To our knowledge, it has so far not been named in DH. Our experiments demonstrate that DH rabbit pups display upregulation of FOXJ1, whereas it was downregulated after TO. Other ciliated-cell marker genes (CCDC19, SPAG17, CCDC39, LRRIQ1, EFHC1, LRRC23, TTC18, WDR16, FANK1, ENKUR, CCDC113) were also upregulated in DH lungs compared to controls, and downregulated by TO to levels lower than in controls. These findings point to an increased number of ciliated cells in DH lungs, which is reversed by TO. Loss of ciliated cells after TO has been previously described in the upper conducting airways in fetal lambs with CDH, although that study did not quantify changes in the more distant airways (Deprest et al., 2000). The study of effects of DH and TO on the variety of individual lung cell types can benefit from modern technology such as single-cell RNA-seq (Treutlein et al., 2014). Those investigators applied this method to the developing distal lung epithelium in the mouse, which revealed a number of putative markers of ciliated cells, Clara cells, alveolar type 1 (AT1) cells and AT2 cells. A number of genes in our data are also present within the top 30 putative markers for each lung cell type (both known and novel) previously identified in the Treutlein et al. (2014) study (Table 2). For instance, the Clara cell markers SCGB1A1, PPAP2B and KDR (VEGFR2) were also upregulated in DH, and subsequently downregulated after TO. These findings are consistent with findings in rodents (Asabe et al., 1998; Santos et al., 2007).

The above is a simplified interpretation and biased presentation of our results. First, we identified a number of discrepancies in results between rabbits and other models. For example, AQP4, TGFBR3, EPAS1 and TGFB-I were upregulated in rabbits with CDH, whereas they were downregulated in the NF rat. In rodents, DH is induced early in gestation by a teratogen interfering directly with the RA pathway, whereas, in rabbits, a surgical defect is created beyond mid-gestation. This means that, in both models, different mechanisms or pathways are involved. Therefore, prenatal interventions might have different effects. For instance, maternal retinoic acid administration in rabbits does not affect the gross anatomical, morphological or proliferative characteristics (Gallot et al., 2008), whereas it rescues lungs in NF rats (Montedonico et al., 2006). Next to that there are also large differences between the time points in gestation and/or lung development at which prenatal interventions are done.

Actually, the largest effect measured was that for WNT inhibitor factor 1 (WIF-1) gene expression [up in DH (FC 8.8107, FDR 0.0153) and down after TO (FC −4.1432, FDR 0.0272)]. However, WIF-1 downregulation was previously reported in the NF rodent lungs during the saccular stage of lung development (Fujiwara et al., 2012). WNT signaling plays an essential role in embryonic development. WIF1 is a target gene of SMAD1 in the developing lung epithelial cells and acts to inhibit WNT proteins (Xu et al., 2011). SMAD1 plays a key role in organogenesis, including in lung development and maturation, and SMAD1-knockout mice display reduced sacculation, which is an important feature of pulmonary hypoplasia (References?). How this needs to be tied together with our rabbit data remains unclear.

In addition, we observed changes that were different from what was previously described in clinical specimens. One example is that of the endothelin receptors EDNRA and EDNRB, which play a complex role in vascular tone. They are overexpressed in the thickened media of the pulmonary arteries of newborns with CDH (de Lagausie et al., 2005), as well as in the NF rodent model (Dingemann et al., 2010). Binding of endothelin 1 to its receptor EDNRA on pulmonary artery smooth-muscle cells results in vasoconstriction, whereas binding to the receptor EDNRB, present on endothelial cells, results in vasodilation mediated by endogenous nitric oxide. In our rabbit model, however, EDNRB was not differently expressed in DH as compared to controls.

We also observed discrepancies between the effects of DH and TO within the rabbit model. For instance c-fos induced growth factor [FIGF; also known as vascular endothelial growth factor D (VEGFD)] was downregulated in CDH as well as following TO. This is a member of the PDGF/VEGF family, and activates VEGFR2 and VEGFR3 receptors. Because it was not differentially expressed in the controls, we assume it was not involved. Conversely, VEGFR2 was upregulated in DH (FC 1.7819, FDR 0.0281) and down to normal levels after TO (FC −2.321, FDR 0.0121), which demonstrates the complex interactions of these factors.

We further had a number of discrepancies between RNA-seq and qPCR, which we touched upon already when reporting on results for RFX3, LRRQL1, BMPR2 and PDE5A. Such discrepancies have been observed in previous reports (Duressa et al., 2013; Shi and He, 2014). This lack of concordance might be due to the poor knowledge of the rabbit genome and the absence of any single nucleotide polymorphism (SNP) database, which could result in poor primer design for certain genes.

Additionally, we acknowledge several other limitations to our study. The sample size is small yet it matches financial limitations and is in line with the sample size of other studies analyzing the gene expression changes by TO in DH. Another is that the study of gene expression in whole lung samples will not detect changes occurring in individual cell types or different areas along the airways. These add to the generic limitations of translational studies in animal disease models, of which the clinical relevance remains uncertain. Another limitation of the study is that the surgical model does not fully recapitulate the etiology of the disease.

In conclusion, we first describe the pattern of dysregulated pulmonary gene expression in lungs of fetal rabbits with surgically induced DH. Interestingly, TO (tracheal occlusion), which is currently investigated as a prenatal surgical method to reverse pulmonary hypoplasia, is associated with a gene expression pattern that is comparable to what was observed in littermates with normal lung development.

The animal experiments were approved by the Ethics Committee of the Faculty of Medicine of the KU Leuven and all animals were treated according to current guidelines on animal wellbeing. In total, 12 does of New Zealand white rabbits (Oryctolagus cuniculus) were operated on and premedicated, after weight estimation, with ketamine 50 mg/kg body weight (Ketamine 1000 CEVA^®; CEVA Santé Animale, Brussels, Belgium), xylazine 6 mg/kg body weight (Vexylan^®; CEVA Santé Animale) and buprenorphine 0.03 mg/kg body weight (Vetergesic^®; Reckitt Benckiser Healthcare, Brussels, Belgium), all injected intramuscularly. General anesthesia was maintained using isoflurane 1.5% (Isoba^® Vet; Abbott Laboratories Ltd, Queensborough, Kent, UK) in oxygen at 1 liter/min via a facemask. Maternal heart rate and oxygen saturation were monitored with a pulse oxymeter (Nellcor^® N-20P; Nellcor Inc., Haasrode, Belgium). Physiologic body temperature was maintained by a heating pad. The doe was placed in the supine position, the abdominal wall was shaved and disinfected with povidone iodine (Isobetadine^®; Asta Medica, Brussels, Belgium) and covered with sterile drapes. Aseptic conditions were maintained throughout all surgical procedures. We operated on two fetuses in each doe, located in the middle of the left or right uterine horn. First diaphragmatic hernia (DH) was created at 23 days gestational age (GA). A second operation was performed at 28 days GA on the previously operated fetuses: either tracheal occlusion (DH+TO group), or sham neck dissection and skin closure (DH group). The latter is performed to include the effects of the fetal surgery per se. Before and after surgery, animals were housed in separate cages at normal room temperature and daylight, with free access to food and water. At 30 days GA, the does were euthanized with an intravenous bolus of a mixture of embutramide 200 mg, mebezonium 50 mg and tetracaine hydrochloride 5 mg (0.3 ml/kg T61; Marion Roussel Hoechst, Brussels, Belgium) after previous premedication as described above. All operated fetuses as well as one control unoperated littermate (control group) were euthanized in utero and harvested by cesarean section to obtain non-ventilated lungs.

During necropsy, the fetus and its lungs were weighed on a scale measuring accurately up to 0.001 g (HF 2000; A&D Instruments, Haasrode, Belgium) and the lung-to-body weight ratio (LBWR) was calculated.

Lung development is reflected by the LBWR and pathologists have defined critical cut-offs for pulmonary hypoplasia. We minimized contamination of the pulmonary tissue with other sources of DNA/RNA by flushing the lungs with saline by vascular access through the right ventricle of the fetal heart. The parenchymal lung tissue, i.e. without the trachea, was stored in RNAlater (Qiagen Benelux B.V., Venlo NL) at −4°C immediately.

RNA isolation was performed on the entire left lung within 4 h of harvesting using the RNeasy mini kit (Qiagen Benelux B.V., Venlo NL). Tissue lysis and homogenization was performed in 1200 µl Buffer RLT using the TissueLyser system (Qiagen Benelux B.V., Venlo, The Netherlands). Following tissue disruption and homogenization, samples were centrifuged for 3 min at 8764 g (14,000 rpm) in a benchtop microcentrifuge. 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 two RNeasy spin columns used per sample. Following wash steps, RNA was eluted in 50 µl RNAse-free H₂O.

Total RNA quantification was performed using the Nanodrop 1000 spectrophotometer (Thermo Scientific, Aalst, Belgium). RNA integrity was assessed using the RNA 6000 Nano Kit and the Bioanalyser (Agilent Technologies, Diegem, Belgium) according to the manufacturer's recommendations.

One µg of total RNA was used as input material for sequencing library preparation, which was performed with the TruSeq RNA Library Preparation Kit (Illumina, Eindhoven, NL) according to the manufacturer's protocol. Fragmentation was performed for 6 min. Eight PCR cycles were used for the PCR enrichment step. Samples were indexed to allow for multiplexing. Sequencing libraries were quantified using the Qubit fluorometer (Life Technologies Europe B.V., Ghent, Belgium). 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 2 nM and pooled for sequencing. Pooled libraries were sequenced in a single lane of an Illumina HiSeq2000 flow cell generating single end 50 bp reads. At least 10 million reads were obtained for all samples (range 12- to 32-million reads).

Fastq files were mapped to the rabbit genome and transcriptome (OryCun2.0 retrieved from Ensembl, release 69) with the Array Studio software (OmicSoft, Cary, NC, USA) using default parameters for single short read data. Between 72% and 78% of reads mapped uniquely to the reference genome and/or transcriptome, and a further 6% to 8% of reads mapped non-uniquely. Between 16% and 18% of reads were unmapped. Expression values were calculated per gene, normalized to ‘reads per kilobase per million reads’ (RPKM) values as previously described (Mortazavi et al., 2008). Data were filtered to remove those genes with expression of <1 RPKM in all samples. The RPKM expression values for the remaining 11,267 genes were subsequently Log2 transformed for downstream PCA, hierarchical clustering, and for the subsequent calculation of fold changes (FC) by statistical inference analysis. All data are made available via ArrayExpress (accession number: E-MTAB-3452).

Hierarchical clustering was performed using the Array Studio software application, selecting ‘Complete’ for link option, and ‘Correlation’ for distance option. The ‘Correlation’ setting is also named ‘Centered Pearson’ 1−corr(x, y), and the source algorithms are further described in http://cran.r-project.org/web/packages/amap/amap.pdf. Clustering is based upon the Log-transformed RPKM gene expression values.

Downstream pathway analysis and creation of networks was performed using the IPA software application (Ingenuity Systems, Redwood City, CA, USA). The iRegulon plugin was used with the Cytoscape analysis software (http://www.cytoscape.org/) for prediction of upstream transcription factors of interest (Janky et al., 2014). The short time-series expression miner (STEM) application (Ernst and Bar-Joseph, 2006) was used to identify gene clusters with similar changes in expression between control and CDH, and following TO.

In order to confirm RNA-seq findings, we selected for qPCR (quantitative polymerase chain reaction) 12 genes with FC>2 or FC<(−2) and the housekeeping genes SDHA, TOP1 and GAPDH, which were previously shown to be among the most stable genes in the CDH rabbit model and after TO (Vuckovic et al., 2013). PCR primers were designed using Primer 3 software and Primer-BLAST. The absence of secondary structure was checked by the Vector NTI program. All the primers were synthesized by Integrated DNA Technologies.

For a given gene, all samples were analyzed in the same qPCR run. 10-µl reactions were prepared using 2× LightCycler 480 SYBR Green I Master (Roche Applied Science) with 2 µl of a 1:10 cDNA dilution and a final concentration of 300 nM of each primer. Data were analyzed using the LightCycler 480 software (Roche Applied Science). The results were quantified using the comparative threshold cycle method as described by Livak and Schmittgen (2001). In addition, serial dilutions were used to create standard curves for relative quantification and the expression of each gene was normalized to at least two of the three housekeeping genes' expression. The sequences and the different primers and the corresponding PCR efficiencies are provided in Table S3.

We thank Rosita Kinnart and Ivan Laermans for their kind assistance in performing the animal experiments.