A meta-analysis study of gene expression datasets in mouse liver under PPARα knockout

Summary

Gene expression profiling of peroxisome-proliferator-activated receptor α (PPARα) has been used in several studies, but there were no consistent results on gene expression patterns involved in PPARα activation in genome-wide due to different sample sizes or platforms. Here, we employed two published microarray datasets both PPARα dependent in mouse liver and applied meta-analysis on them to increase the power of the identification of differentially expressed genes and significantly enriched pathways. As a result, we have improved the concordance in identifying many biological mechanisms involved in PPARα activation. We suggest that our analysis not only leads to more identified genes by combining datasets from different resources together, but also provides some novel hepatic tissue-specific marker genes related to PPARα according to our re-analysis.

1. Introduction

The transcription factor peroxisome-proliferator-activated receptor α (PPARα) is one of the nuclear-hormone receptors (NRs), which need to be ligand-activated (Issemann & Green, 1990; Bookout et al., 2006; Michalik et al., 2006). In mammals, PPARα is thought to participate in several biological processes, especially in the nutrient metabolism and inflammation. The expression levels of PPARα are high in most of tissues with active fatty acid catabolism, such as liver, heart, kidney, brown adipose tissue, muscle, small and large intestines (Issemann & Green, 1990; Braissant et al., 1996; Devchand et al., 1996). Particularly, the liver is a central player in the whole body energy homeostasis by its ability to orchestrate fatty acid and glucose metabolism, and it plays an important role in fatty acid oxidation and this catabolic energy burning is regulated by PPARα (Reddy & Hashimoto, 2001). Nowadays, the PPARα-null mouse model is frequently used according to the study of the critical role for PPARα in the cellular fasting response (Leone et al., 1999), since we could get the significantly PPARα-dependent genes by the comparison between wild-type (WT) and knockout (KO) mice. These identified genes may be associated with metabolism or diseases, such as obesity and diabetes. During the past 15 years, numerous studies with PPARα null mice have demonstrated the critical roles played by this receptor in energy metabolism, hepatic steatosis, inflammation, cardiac pathophysiology, cell-cycle alterations and hepatocarcinogenesis (Lee et al., 2002; Lefebvre et al., 2006). Particularly in recent years, microarray technology has been highly used to map PPARα-dependent genes and further characterizes PPARα function in different tissues in genome wide, especially for liver. We could establish gene expression patterns that would facilitate understanding of almost every aspect of cellular and molecular biology related to PPARα. Since mass databases have been created using a variety of gene expression data by different researchers, many studies have considered gene expression determination by PPARα activation at the transcriptome level. In these liver studies, PPARα was shown to be critical for the coordinated transcriptional activation of genes involved in lipid catabolism, including cellular fatty acid uptake and activation, mitochondrial β-oxidation, peroxisomal fatty acid oxidation, ketone body synthesis, fatty acid elongation and desaturation and apolipoprotein synthesis (Leone et al., 1999; Mandard et al., 2004). However, the reproducibility of datasets and the credibility of results are both questionable, especially for the fact that the transcription patterns of genes revealed by different databases frequently contradict each other. Therefore, combining the gene-expression datasets generated at different laboratories or conditions in mouse hepatic tissue is needed to identify tissue-specific pathways or genes involved in PPARα activation. It may be helpful for us to get more powerful and comprehensive knowledge on the biological functions of PPARα in mice liver.

A non-parametric method named rank product (RP) is one of the most useful statistical models to identify the significant genes in the analysis of microarray experiments, which is based on the calculation of RP from replicate datasets (Breitling et al., 2004; Nemhauser et al., 2004). The RP analysis can also be applied to datasets from multiple origins. Owing to the lack of experimental standards for microarray experiments that leads to heterogeneous datasets, direct comparison is not possible. Instead of using actual expression data, this approach can combine the gene rank from different origins together to select differentially expressed genes. The smaller the RP value, the more important the observed gene is. The significance of the list of up- or down-regulated genes is based on the estimated percentage of false positive (pfp) predictions, which is also known as false discovery rate (FDR). It is true that the meta-analysis of related microarray datasets based on this approach would be powerful to identify regulatory factors required for response to PPARα KO.

2. Materials and methods

(i) Microarray datasets collection and preprocessing

We searched Gene Expression Omnibus (GEO, www.ncbi.nlm.nih.gov/geo/) for the gene expression profiling studies related to PPARα KO. Data were included in our re-analysis if they met the following conditions: (1) the data are genome-wide; (2) comparison was conducted between PPARα KO and WT; (3) complete microarray raw or normalized data are available. Finally, we chose the datasets GSE9786 and GSE8295 for our re-analysis studies, which were respectively contributed by Rosen and Hooiveld in different periods, both for mouse liver (Rakhshandehroo et al., 2007; Rosen et al., 2008). Both laboratories measured gene expression under two classes (WT and KO) with similar condition. The details of them were shown in Table 1.

Table 1. The details of individual study and meta-analysis of the chosen two datasets

GSE9786 and GSE8295 were both related to PPARα in mouse liver using the platform of Affymetrix Mouse 430 2.0 GeneChips, the sample sizes used were both 16 and the numbers of replicates were both 4. However, GSE9786 was contributed by Rosen on 4 February 2008 for 7 days treatment and GSE8295 was contributed by Hooiveld on 24 July 2007 for 5 days treatment. According to individual analysis on the above two with comparison of WT to KO, there were 289 and 1143 significantly down-regulated and 198 and 624 up-regulated genes, respectively, in the study of GSE9786 and GSE8295. The numbers of significantly enriched down- or up-regulated pathways were 14 and 9 for GSE8295 and 13 and 9 for GSE9786. Based on meta-analysis (RP) of the chosen two datasets, there appeared to 1695 down-regulated and 1708 up-regulated genes as well as 13 and 6 enriched down- and up-regulated pathways.

In GSE9786, 129S1/SvlmJ WT and PPARα-null mice were exposed to 3 mg/kg/day PFOA or water for 7 days. Total RNA was isolated from liver samples and gene expression analysed using Affymetrix Mouse 430 2.0 GeneChips. Data from 16 samples, with four mice in each of the four treatment groups, were analysed. In GSE8295, 3–5-month-old male pure bred WT (129S1/SvImJ) and PPARα-null (129S4/SvJae) mice were used. WT and PPARα-null mice were treated with the synthetic PPARα ligand Wy14643 (0·1% w/w) mixed in the food or normal food (control) for 5 days (n=4 per group). Liver total RNA from biological replicates was hybridized onto Affymetrix mouse genome 430 2.0 GeneChip arrays.

For the assessment of the influence of preprocessing on the comparison, the data preprocessing was performed using software packages developed in version 2.6.0 of Bioconductor and R version 2.10 (Gentleman et al., 2004). Each Affymetrix dataset was background adjusted, normalized and log 2 probe-set intensities calculated using the Robust Multichip Averaging (RMA) algorithm in Affy package (Gautier et al., 2004).

(ii) Individual study of each dataset

In our individual study of each dataset, statistical analysis was performed by two-way ANOVA with a Benjamini and Hochberg FDR (BH-FDR=0·05) for multiple testing correction followed by Tukey's post-hoc tests (Park et al., 2003). Differentially expressed genes between two neighbour stages were identified by 2-fold changes. Clustering on groups and genes was performed based on the identified genes' expression using the method of Hierarchical clustering. All of these processes were performed using software packages developed in version 2.6.0 of Bioconductor and R version 2.10.0.

(iii) Meta-analysis on two datasets

In order to identify regulatory elements required for response to PPARα KO, we divided the samples in each dataset into two classes. Concretely, for GSE9786, eight samples from GSM247116 to GSM247123 were set to be the first class (WT) and the other eight samples from GSM247124 to GSM247131 were set to be the second class (KO); for GSE8295, eight samples from GSM205784 to GSM205788 and GSM205804 to GSM205816 were set to be the first class (WT) and the other eight samples from GSM205789 to GSM205803 and GSM205817 to GSM205820 were set to be the second class (KO). Our meta-analysis on the chosen two datasets was performed by the RankProd package to identify up- or down-regulated genes under WT against the treatment of PPARα KO. The cut-off of pfp value was chosen as 0·05.

(iv) Pathway analysis

Further classifications of pathway analysis were performed by using the Database for Annotation, Visualization and Integrated Discovery (DAVID) (Dennis et al., 2003; Huang et al., 2009) following the above significance analysis, revealing overrepresented pathways of identified genes associated with PPARα KO. The pathway with P<0·01 was identified as the significantly associated one.

3. Results and discussion

(i) Individual significance analysis results

Through our individual analysis on the above two datasets by two-ANOVA–FDR test with FDR ⩽5% and Fold Change (FC) ⩾2 with comparison of WT to KO, in total 1767 genes were significantly identified in GSE8295 including 1143 down-regulated and 624 up-regulated. In contrast, the total number in GSE9786 was relatively lower as 487 including 289 down-regulated and 198 up-regulated (Table 1). The details of significantly expressed genes were shown in Additional Files 1 and 3 for GSE8295 and GSE9786, respectively. To explore the gene expression patterns of these two datasets, we performed hierarchical clustering on the identified genes across all samples and made a comparison between GSE8295 and GSE9786. As shown in Fig. 1, there were remarkable differences not only in the number of significant genes, but also in the expression patterns. We next sought to determine the pathways associated with the regulation of PPARα KO in individual dataset, analysis of overrepresented pathways were performed by using DAVID with P-value cutoff as 0·01. The numbers of significantly enriched down- or up-regulated pathways were 14 and 9 for GSE8295 and 13 and 9 for GSE9786 (Table 1). The full lists of significantly related pathways were shown in Additional Files 2 and 4, respectively, for GSE8295 and GSE9786.

Fig. 1. Heatmaps of the selected two datasets. The heatmaps were showing hierarchical clusterings on the filtered entities using all samples included in GSE8295 (Figure 1A) and GSE9786 (Figure 1B). X-axis represents all samples and their groups included in each dataset; Y-axis represents all the filtered probe sets. Low gene expression was shown in green and high gene expression was in red.

In summary, based on our individual study on the chosen two datasets, there were differences between GSE8295 and GSE9786 in many aspects, such as the number of PPARα-dependent genes and the significance of related pathways. These differences were possibly due to the quality and quantity of samples chosen in each microarray experiment and the treatment condition may also play essential roles in effecting the accuracy of results. Therefore, a gene expression meta-analysis on these two datasets is necessary to get more reliable results.

(ii) Meta-analysis results and comparison with individual study

Through our gene expression meta-analysis of the chosen two datasets based on the RP method, there appeared to be totally 3403 significantly identified genes including 1695 down-regulated genes and 1708 up-regulated genes with the pfp value <0·05. The details of all of them were shown in Additional file 5 and the top 10 ranked down- or up-regulated genes were showing in Table 2. On comparison with the results of individual studies, it was noted that the test of meta-analysis lead to more identified genes by combining datasets from two different origins together. In addition, according to the pfp values and P values of identified genes by our meta-analysis, the significances of down- or up-regulated genes were improved and the results become more powerful.

Table 2. The top 10 ranked significant genes in meta-analysis

Among the top 10 ranked up-regulated genes, SULT1E1, one member of the sulphotransferase family whose function are to catalyse the sulphate conjugation of many hormones, neurotransmitters, drugs and xenobiotic compounds, was identified at the first rank with the RP value as 8·7008. The human gene of SULT1E1 was reported to encode a protein that transfers a sulpho moiety to and from oestrone and may primarily control levels of oestrogen receptors (Gong et al., 2008). Comparison of murine and human oestrogen sulphotransferase inhibition in vitro and in silico, it suggested that catalytically inactive binding modes, other than the one observed in the crystal structures, be possible in SULT1E1 (Stjernschantz et al., 2010). Moreover, the function of SULT1E1 was also found to be related to steroid sulphotransferase activity by observing the bacterial expression and characterization of a cDNA for human liver oestrogen sulphotransferase (Falany et al., 1995). These findings may be associated with the regulation of PPARα KO in lipid metabolism. Serpinb1a also named leukocyte elastase inhibitor a, one member of serine (or cysteine) peptidase inhibitors, was the second ranked one (RP=17·3453). The regulation of pulmonary innate immunity by serpinb1 was thought to be non-redundant and required to protect two key components, the neutrophil and surfactant protein-D (Benarafa et al., 2007). The most frequently identified gene in our top 10 was Meg3, which was represented by five probe sets all identified as significant. Meg3 named as maternally expressed 3, which is an imprinted gene with preferential expression from the maternal allele and closely linked to or co-expressed with reciprocally imprinted paternally expressed Dlk1 gene. This co-regulation suggests a causative role in the pathologies found in uniparental disomy animals, characterized by defects in skeletal muscle maturation, bone formation, placenta size and organization and prenatal lethality. Meg3 may play an important role in control of vascularization in the brain and may function as a tumour suppressor in part by inhibiting angiogenesis (Gordon et al., 2010). Another top gene Clta named Clathrin, light polypeptide (Lca) is a large, soluble protein composed of heavy and light chains. It functions as the main structural component of the lattice-type cytoplasmic face of coated pits and vesicles which entrap specific macromolecules during receptor-mediated endocytosis. This gene encodes one of two clathrin light chain proteins that are believed to function as regulatory elements. A B-Myb complex containing clathrin and filamin, which was previously reported to be required for normal localization of clathrin at the mitotic spindle stabilizing kinetochore fibres, may contribute to its role in chromosome stability, possibly, tumour suppression (Yamauchi et al., 2008). The seventh ranked up-regulated gene NNMT is named as nicotinamide N-methyltransferase, nicotinamide treatment may induce NNMT enzyme activity and led to an increase in adipose tissue homocysteine secretion in tissue culture (Riederer et al., 2009). The last one Lrtm1 named leucine-rich repeats and transmembrane domains 1, a gene of unknown function embedded within CACNA2D3 encoding a voltage-dependent calcium channel, which would be needed for further study in the research of our interest.

On the other hand, among the top 10 ranked down-regulated genes, the gene of Vnn1 named vanin 1 was identified twice with the represented probe sets at the first rank (RP=4·9119) and the seventh rank (RP=26·6864). Mouse vanin 1 was reported to be involved in many biological processes in mouse, from thymus homing to sexual development and indicated to be cytoprotective for islet β cells and regulates the development of type 1 diabetes (Di Leandro et al., 2008; Roisin-Bouffay et al., 2008). The second one was Cox6b2 named cytochrome c oxidase subunit 6B2 participating in some disease pathways with the function of cytochrome c oxidase activity (Pagliarini et al., 2008); however, the exact knowledge about regulation of this gene on lipid metabolism is limited. Another identified gene Cyp4a14 was a member of the family of cytochrome p450 enzymes (Cyps), which are major phase-I xenobiotic-metabolizing enzymes. The previous study on mice liver indicated that perfluorodecanoic acid (PFDA) may increase Cyp4A14 mRNA expression in WT mice, but much less in PPARα-null mice (Cheng & Klaassen, 2008). Moreover, a recent study suggests that the Cyp4a14-deficient mouse may be a useful model for evaluation of NO/20-HETE interactions, the production of which in the kidney is thought to be involved in the control of renal vascular tone and tubular sodium and chloride reabsorption (Fidelis et al., 2010). The gene Plin5 named perilipin 5 was also identified twice in our study. The perilipin is an adipose differentiation-related protein, tail-interacting protein of 47 kd (PAT) proteins that are thought to regulate intracellular lipid droplets (LDs) turnover by modulating lipolysis (Bell et al., 2008). Mouse Rad51l1 is a radiation-inducible gene that regulates cell-cycle progression. The previous study based on disruption in embryonic stem cells by homologous recombination suggested that Rad51l1 would play a role in cell proliferation and early mice embryonic development, perhaps through interaction with p53 (Shu et al., 1999). The gene of Raet1e named as retinoic acid early transcript 1E may promote tumour surveillance and autoimmunity by engaging the activating receptor NKG2D on natural killer (NK) cells and T cells (Oppenheim et al., 2005). Two hypothetical protein LOC100044218 and LOC630172 were both reported here as the top ones, whose annotation status are not on current assembly.

The subsequent pathway analysis revealed that there were 13 down-regulated and six up-related pathways participating in regulating PPARα activation. The whole information of significantly identified pathways in our meta-analysis are shown in Table 3, including the pathway names, the count of genes in each identified pathway, the significance (P-value) of identified pathways and the related genes in each pathway. According to our analysis based on KEGG pathway database (http://www.genome.jp/kegg/), we can map our identified pathways into different functional classification. For up-regulated pathways, the top significant pathway of complement and coagulation cascades (P=4·06×10⁻¹⁴) was mapped to immune system, metabolism of xenobiotics by cytochrome p450 (1·05×10⁻⁵) was mapped to xenobiotics biodegradation and metabolism, glycine, serine and threonine metabolism (2·00×10⁻⁴) was mapped to amino acid metabolism, nitrogen metabolism (8·06×10⁻³) was mapped to energy metabolism, both arachidonic acid metabolism (4·60×10⁻³) and linoleic acid metabolism (6·26×10⁻³) were mapped to lipid metabolism. For down-regulated pathways, besides five pathways including fatty acid metabolism (1·93×10⁻¹⁸), polyunsaturated fatty acid biosynthesis (3·63×10⁻⁶), fatty acid elongation in mitochondria (8·83×10⁻⁴), glycerophospholipid metabolism (5·83×10⁻³) and glycerolipid metabolism (9·20×10⁻³) in the classification of lipid metabolism, the pathway of valine, leucine and isoleucine degradation (3·80×10⁻⁷) in the classification of amino acid metabolism, the pathway of caprolactam degradation (2·39×10⁻³) in the classification of xenobiotics biodegradation and metabolism, more classifications were found including PPAR signalling pathway (2·04×10⁻¹⁶) of endocrine system, two pathways of cell cycle (2·56×10⁻⁵) and p53 signalling pathway (1·28×10⁻⁴) in the classification of cell growth and death, three pathways including butanoate metabolism (3·23×10⁻⁴), pyruvate metabolism (1·11×10⁻³), and propanoate metabolism (1·27×10⁻³) in the classification of carbohydrate metabolism. A variety of biological processes regulated by PPARα in mouse liver from our re-analysis results were summarized in Fig. 2.

Fig. 2. PPARα regulates a variety of biological processes in mouse liver. Summary of functional implication of PPARα activation as assesses by our meta-analysis based on KEGG pathways enrichment. The names in shaded boxes represent eight identified KEGG pathway maps, following with significant pathways' names, up-regulated in red and down-regulated in green under PPARα(−/−).

Table 3. The significantly identified pathways in meta-analysis

On comparison with individual results of pathway analysis, there were remarkable differences not only in the terms but also in the significance of identified pathways. The information on comparison of the significance of down- or up-regulated pathways identified among individual analysis and meta-analysis was shown in Figs 3 and 4. Actually, there appeared to be some common pathways among these studies, such as two up-regulated pathways including complement and coagulation cascades and metabolism of xenobiotics by cytochrome p450 and seven down-regulated pathways including PPAR signalling pathway, fatty acid metabolism, propanoate metabolism, the pathway of valine, leucine and isoleucine degradation, pyruvate metabolism, butanoate metabolism and glycerolipid metabolism. However, their significances were not the same. In addition, the results based on meta-analysis exhibited some novel associated pathways such as polyunsaturated fatty acid biosynthesis, fatty acid elongation in mitochondria and caprolactam degradation, related to diverse metabolisms.

Fig. 3. Comparison of the significance of down-regulated pathways identified among individual analysis and meta-analysis. The chart was showing the significance of down-regulated pathways identified under each analysis. X-axis represents −log (P values); Y-axis represents the names of each identified pathways. The result of individual analysis of GSE9786 was marked in green, GSE8295 was marked in red and the meta-analysis was marked in blue.

Fig. 4. Comparison of the significance of up-regulated pathways identified among individual analysis and meta-analysis. The chart was showing the significance of up-regulated pathways identified under each analysis. X-axis represents −log (P values); Y-axis represents the names of each identified pathways. The result of individual analysis of GSE9786 was marked in green, GSE8295 was marked in red and the meta-analysis (based on the method of RP) was marked in blue.

In order to identify the gene expression patterns of PPARα, a type of meta-analysis based on the method of RP was employed to combine datasets from different origins (datasets of GSE9786 and GSE8295) which were contradicted in some aspect of analysis results. As a result, it increases the power of the identification of differentially expressed genes and significant pathways in this study. Some novel hepatic tissue-specific marker genes as well as pathways related to PPARα were detected. These may be helpful for us to further uncover the regulation mechanisms of PPARα in mice liver.