Distinct gene regulatory programs define the inhibitory effects of liver X receptors and PPARG on cancer cell proliferation
© The Author(s). 2016
Received: 18 April 2016
Accepted: 14 June 2016
Published: 11 July 2016
The liver X receptors (LXRs, NR1H2 and NR1H3) and peroxisome proliferator-activated receptor gamma (PPARG, NR1C3) nuclear receptor transcription factors (TFs) are master regulators of energy homeostasis. Intriguingly, recent studies suggest that these metabolic regulators also impact tumor cell proliferation. However, a comprehensive temporal molecular characterization of the LXR and PPARG gene regulatory responses in tumor cells is still lacking.
To better define the underlying molecular processes governing the genetic control of cellular growth in response to extracellular metabolic signals, we performed a comprehensive, genome-wide characterization of the temporal regulatory cascades mediated by LXR and PPARG signaling in HT29 colorectal cancer cells. For this analysis, we applied a multi-tiered approach that incorporated cellular phenotypic assays, gene expression profiles, chromatin state dynamics, and nuclear receptor binding patterns.
Our results illustrate that the activation of both nuclear receptors inhibited cell proliferation and further decreased glutathione levels, consistent with increased cellular oxidative stress. Despite a common metabolic reprogramming, the gene regulatory network programs initiated by these nuclear receptors were widely distinct. PPARG generated a rapid and short-term response while maintaining a gene activator role. By contrast, LXR signaling was prolonged, with initial, predominantly activating functions that transitioned to repressive gene regulatory activities at late time points.
Through the use of a multi-tiered strategy that integrated various genomic datasets, our data illustrate that distinct gene regulatory programs elicit common phenotypic effects, highlighting the complexity of the genome. These results further provide a detailed molecular map of metabolic reprogramming in cancer cells through LXR and PPARG activation. As ligand-inducible TFs, these nuclear receptors can potentially serve as attractive therapeutic targets for the treatment of various cancers.
KeywordsNuclear receptors LXR PPARG Cell proliferation Metabolism Energy homeostasis Transcription Chromatin state dynamics RNA-seq ChIP-seq
Metazoan systems are exposed to a multitude of extracellular stimuli, including both compounds that are endogenously synthesized and those obtained through diet [1–3]. In response to these environmental signals, cells activate intricate transcriptional programs [1–3]. The careful interplay between extracellular cues and intracellular gene expression responses maintains normal cellular equilibrium and is key for cellular differentiation and development [4–7]. Despite this importance, more work is needed for a thorough, molecular characterization of the temporal gene regulatory responses to diverse external signals.
The liver X receptors (LXRs), also referred to as Nuclear Receptor Subfamily 1, Group H, Member 2 (NR1H2) and Nuclear Receptor Subfamily 1, Group H, Member 3 (NR1H3), and peroxisome proliferator-activated receptor gamma (PPARG), also referred to as Nuclear Receptor Subfamily 1 Group C Member 3 (NR1C3), are ligand-activated nuclear receptor transcription factors (TFs) that respond to oxysterols  and fatty acids , respectively. LXRs are comprised of two isoforms (alpha and beta; LXRA and LXRB) that have related functions and tissue expression profiles, including a common affinity for oxysterols , while PPARG is activated by distinct fatty acids . LXRs and PPARG maintain energy homeostasis by regulating lipid [10, 11] and glucose [12, 13] metabolism. Intriguingly, they are also key modulators of inflammatory responses . These metabolic TFs are also involved in several diseases, including diabetes [12, 15–19], obesity [20–24], atherosclerosis [25–29], and cancer [30–33].
Protective functions have consistently been reported for LXRs in a variety of diverse cancers. For instance, diets rich in plant-derived phytosterols, a putative LXR agonist , were shown to decrease risk of breast , gastric , and lung  cancers. Importantly, a protective effect for phytosterols in colorectal cancer has also been demonstrated using cell lines and animal models . Administration of an LXR agonist also suppressed the growth of LNCaP prostate cancer cell tumor xenografts in mice  while genetic ablation of LXRB results in gallbladder carcinogenesis . Moreover, specific LXR activation leads to anti-proliferative effects in breast , prostate , and colorectal  cancer cells. Several mechanisms have been proposed for these LXR-mediated effects. For instance, cell cycle inhibition [32, 33], induction of apoptosis , and ligand deprivation  have all been described. The role of PPARG in tumorigenesis is complex and its effects appear to be tissue-specific. Several studies have suggested protective effects in breast , hepatic , and lung  cancers. A recent study characterized a PPARG anti-proliferative effect in lung cancer that was proposed to be due to the regulation of cellular reactive oxygen species via activation of pyruvate dehydrogenase kinase 4 (PDK4) and beta-oxidation of fatty acids . However, in colorectal cancer, both tumor-suppressive  and tumor-facilitating [30, 47] functions have been proposed for PPARG.
To better understand the underlying gene regulatory functions of LXRs and PPARG that lead to metabolic reprogramming in cancer cells, we performed a comprehensive genome-wide analysis of LXR and PPARG activity in HT29 colorectal cancer cells. Despite common physiological effects on energy homeostasis, paradoxically, the genome-wide molecular programs elicited by LXRs and PPARG appear to be temporally distinct. Our results highlight the advantages of integrating multiple levels of transcriptional regulation rather than relying on a single feature. From a broader perspective, this work supports the notion that common phenotypic effects can be mediated through distinct gene regulatory mechanisms.
HT29 colorectal cancer cells were obtained from ATCC and grown under recommended cell culture conditions using McCoy’s 5A media containing 10 % fetal bovine serum (FBS) and 1 % penicillin/streptomycin. GW3965 (Sigma), T0901317 (Santa Cruz), and rosiglitazone (Sigma) powder were diluted in DMSO and added to media at a final concentration of 10 μM. Prior to drug treatments, HT29 cells were grown in phenol-red free McCoy’s 5A media containing 10 % charcoal/dextran-treated FBS and 1 % penicillin/streptomycin.
Cell phenotype assays
Cell-based phenotypic assays were performed in 96-well cell culture plates. Cell proliferation measurements were determined using the CyQuant assay (Invitrogen) for DNA content and CellTiter-Glo Luminescent Cell Viability Assay (Promega) for ATP levels. Glutathione levels for oxidative stress were determined using the GSH-Glo Glutathione Assay (Promega).
For analyses of metabolites, cell cultures were washed with phosphate-buffered saline, molecular biology-grade water and then subsequently pelleted for mass spectrometry. Cell pellets were resuspended in 50 % methanol, shaken at 4 °C for 30 minutes and centrifuged at 12,000 rpm for 10 minutes. Duplicate aliquots of supernatant were dried at 55 °C for 60 minutes using a vacuum concentrator system (Labconco). Derivatization by methoximation and trimethylsilyation was performed as previously described . All derivatized samples were analyzed on a Leco Pegasus 4D system (GCxGC-TOFMS), controlled by the ChromaTof software (Leco, St. Joseph, MI, USA). Peaks were identified by spectral match using the NIST, GOLM, and Fiehn libraries (Leco) and confirmed by running derivatized standards (Sigma). Peaks present in less than two-thirds of samples were excluded from further analysis. Sample replicates were averaged and peak areas were sum normalized prior to comparisons.
Transfection of LXRA transgene plasmid constructs
Flag-tagged (3-prime 3X-Flag tags separated by glycine spacers) LXRA cDNA plasmid constructs were transfected (Fugene) into HT29 cells. HT29 cells were subsequently selected using G418 (Invitrogen) and expanded prior to GW3965 drug treatment and ChIP-seq experimentation.
ChIP-seq assays and analysis
ChIP-seq was performed as previously outlined . Antibodies for H3K27ac (Abcam, ab4729), H3K36me3 (Abcam, ab9050), RNA polymerase II (RNAP2; Abcam, ab5408), LXRB (Active Motif, 61177), PPARG (Santa Cruz, sc-7273), and Flag (Sigma, F1804) were used.
RNA-seq experimentation and analysis
RNA-seq was performed as previously outlined  using Nextera transposases (Illumina). Briefly, cells were lysed with Buffer RLT (Qiagen) containing 10 % beta-mercaptoethanol. The Norgen Animal Tissue RNA Purification Kit (Norgen Biotek) was used to isolate RNA from cells. The Dynabead mRNA Purification Kit (Invitrogen) beads were used to isolate mRNA transcripts and cDNA synthesis was performed using Superscript reverse transcriptase (Invitrogen). Nextera transposases (Illumina) were used to fragment cDNA prior to next-generation sequencing. All experimental treatments at all time points were performed in quadruplicate.
Two-sided Student’s t-tests were performed to assign significance for cellular phenotypic assays and metabolite levels. Only reproducible binding sites identified by replicate ChIP-seq experiments were reported and used for downstream analyses. ChIP-seq peaks were identified using MACS peak caller , while motif analysis was performed using MEME  on the top 500 most enriched binding sites. For normalized ChIP-seq read-depth analyses, the number of reads were tabulated across binding site coordinates and normalized to the total number of aligned reads obtained for each ChIP-seq experiment. For mapping reads, we utilized a 100-bp sequence centered on the peak summit of each binding site. For ranking temporal changes in read depth, we generated a linear model comparing replicate ChIP-seq read enrichment values between two time points. Sites with altered enrichment were subsequently divided by slope and p value ranked. The top 50 % of p value-ranked sites from each category were used for cofactor motif analyses. For assigning TF motif fold enrichments, we compared the total number of motifs found across a set of binding sites with the number of motifs identified after scrambling binding site sequences. For RNAP2 read depth promoter analyses, sequences within 500 bp of transcription start sites were used as promoters. Connectivity maps were generated using Cytoscape (http://www.cytoscape.org/) and Enrichment map . Kyoto Encyclopedia of Genes and Genomes (KEGG) gene pathway enrichments were generated using Gene Set Enrichment Analysis (GSEA; http://software.broadinstitute.org/gsea/msigdb/annotate.jsp). Wilcoxon rank sum tests were used to assign significance for H3K36me3 enrichment analyses. Differential gene expression was determined using DESeq . High-confidence GW3965 + T0901317 responsive genes were assigned for all targets that passed significance and fold change cutoffs across both drug treatments. For comparison of gene expression with RNAP2 and H3K36me3 occupancy, if a gene harbored multiple promoters and/or transcripts, the promoter and/or gene body coordinate exhibiting the strongest read enrichment was used.
LXRs and PPARG generate common phenotypic effects
To better characterize cellular metabolic alterations, we next assessed metabolite levels through mass spectrometry of HT29 cell lysate after drug treatment (Fig. 1d–f). We focused our analyses on carbohydrate metabolism and evaluated key metabolites involved in oxidative phosphorylation and aerobic glycolysis (Fig. 1d). We found that LXR and PPARG activation led to significant decreases in glucose and citrate levels, as well as to significant increases in pyruvate and lactate levels after drug treatment (Fig. 1e, f), supporting an increase in cellular aerobic glycolysis at the expense of oxidative phosphorylation. Although both drug treatments generated the same cellular response, the temporal dynamics of these metabolic effects were distinct; rosiglitazone produced more substantial effects at 24 h, while several of the GW3965-mediated changes were evident only after 48 h of drug treatment. Apart from glucose pathway metabolites, we also identified a variety of compounds that were significantly altered in response to drug treatment (Additional file 1: Figure S1). In all these cases, LXR and PPARG activation generated identical alterations, highlighting a common mechanism of metabolic reprogramming.
Distinct transcriptional profiles from LXR and PPARG signaling
Total number of differentially regulated genes
Fold change ≥1.5
Fold change ≤ −1.5
Fold change ≥2
Fold change ≤ −2
Fold change ≥1.5
Fold change ≤ −1.5
Fold change ≥2
Fold change ≤ −2
GW3965 + T0901317
Analysis of overlapping LXR- and PPARG-responsive genes
Overlap 24 h (%)
Concordant 24 h
Discordant 24 h
Percentage concordant 24 h
Overlap 48 h (%)
Concordant 24 h
Discordant 24 h
Percentage concordant 48 h
GW3965 (fc ±2)
32 (25.8 %)
203 (18.2 %)
T0901317 (fc ±2)
9 (15.5 %)
189 (16.9 %)
GW3965 (fc ±1.5)
119 (38.9 %)
790 (35.1 %)
T0901317 (fc ±1.5)
60 (35.1 %)
775 (32.1 %)
We performed genetic pathway enrichment analyses based on this list of differentially expressed genes. Gene Ontology enrichment mapping  identified wide alterations in pathways affected by GW3965 treatment between 24 and 48 h (Fig. 2b, c). For instance, the 24-h network consisted predominantly of up-regulated pathways involved in lipid metabolism, signal transduction, small molecule biosynthesis, and transmembrane transport. By contrast, at 48 h, we identified a large set of down-regulated pathways involved in cell cycle regulation, transcription, and tissue morphogenesis. This LXR network topology was distinct from the PPARG response (Additional file 1: Figure S4), where 24- and 48-h rosiglitazone inductions produced concordant pathway enrichments predominantly involved in cell metabolism (up-regulated) and membrane transport (down-regulated).
We also performed enrichment analyses using the KEGG pathway dataset from GSEA . At 24 h, the GW3965 and T0901317 drug treatment produced significantly up-regulated pathways involved in lipid and energy metabolism (PPAR signaling, fatty acid metabolism, pyruvate metabolism, etc.), including ABC transporters. Although these patterns were maintained at 48 h, several key pathways and processes, including MAPK, TGF-beta, cell cycle, and Wnt, were significantly down-regulated by both drugs. For rosiglitazone treatment, the up-regulated and down-regulated KEGG pathways were more consistent between 24 and 48 h. A detailed list of all KEGG pathway enrichments can be found in Additional file 2: Table S1. Pathway enrichments for genes commonly induced by both LXR and PPARG activation are also tabulated in Additional file 3: Table S2.
Chromatin state dynamics support cellular transcriptional effects
We further determined whether the more gradual GW3965 gene expression coincides with slower chromatin state dynamics. We performed replicate ChIP-seq experiments for the histone-3 lysine-27 acetylation (H3K27ac) chromatin modification under control (i.e., DMSO carrier only) culture conditions and after 24- and 48-h drug treatments. Because H3K27ac is a modification that appears at active regulatory elements , this temporal characterization provides a genome-wide picture illustrating changes to the cis-regulatory architecture upon LXR and PPARG activation. Using data from control culture conditions as a reference, we identified H3K27ac sites that were gained or lost after drug treatment at both time points (Fig. 2d; Additional file 1: Figure S5). We next compared the chromatin responses between drug treatments and identified a substantially larger number of H3K27ac sites that were gained at 24 h in response to rosiglitazone (6500 versus 3162 loci; 51 % more), mirroring transcriptional effects. By 48 h, the difference between both drugs substantially diminished (7179 versus 6273 loci; 13 % more). Interestingly, the degree of H3K27ac site loss was similar between drug treatments (Fig. 2d; Additional file 1: Figure S5). We further identified a substantial amount of overlap in H3K27ac binding events that were gained between drug treatments and their associated genetic pathway enrichments (Additional file 1: Figure S6).
We next compiled all time points for each drug treatment separately to identify loci showing dynamic changes in H3K27ac enrichment across the entire duration of treatment (Additional file 1: Figure S7). We identified all GW3965- and rosiglitazone-induced regulatory sites showing differences in normalized read depth from 0 to 24 h and from 24 to 48 h of drug treatment. The vast majority of dynamic H3K27ac sites exhibited a gradual increase or decrease in H3K27ac occupancy over time, suggesting these analyses are identifying nuclear receptor-mediated regulatory events. Further validating our approach and conclusions above, regulatory sites gaining H3K27ac occupancy throughout the time course were enriched near the promoters of genes up-regulated at 48 h after GW3965 and rosiglitazone drug treatments (Additional file 1: Figure S8). This enrichment was also observed using GW3965 + T0901317 target genes (Additional file 1: Figure S9). We further performed Genomic Regions Enrichment of Annotations Tool (GREAT) pathway analyses , which confirmed that sites with temporally increasing H3K27ac occupancy were enriched near relevant genes involved in cell metabolism and/or implicated in colorectal cancer (Additional file 1: Figures S10 and S11). Taken together, these additional data further support the notion that LXR and PPARG direct temporally distinct gene regulatory responses through alterations in chromatin state.
LXRB and PPARG have temporally opposing binding patterns
We next integrated the 2-h and 48-h binding information for LXRB and PPARG proteins, as these datasets illustrate temporal occupancies of TFs under endogenous promoter control. Intriguingly, LXRB and PPARG generated opposing profiles (Fig. 3c, d); the number of LXRB binding events increased significantly with time (514 and 5928 binding sites specific to 2 and 48 h, respectively), while PPARG binding exhibited strongest occupancy at 2 h (6454 and 575 binding sites specific to 2 and 48 h, respectively). As these results relied on qualitative metrics, we next applied quantitative approaches to validate these observations. We pooled all LXRB and PPARG binding events for both time points and calculated the normalized sequencing read depth for all ChIP-seq experiments at this complete set of sites. In agreement with our data above, the resulting histograms for all LXRB and PPARG binding sites identified global shifts in distribution of binding enrichment (Fig. 3e, f); LXRB exhibited stronger occupancy at 48 h for the majority of all binding events, while PPARG binding showed the opposing pattern. Taken together, these data highlight temporally distinct, genome-wide occupancy for LXRB and PPARG proteins and further support the more immediate cellular response from rosiglitazone compared with GW3965 treatment.
AP1 TFs are key cofactors of LXRs and PPARG
We next performed ChIP-seq to analyze the binding of JunD, a known heterodimer of the AP1 TF complex , to validate these computational observations. We performed replicate JunD ChIP-seq experiments under control culture conditions (DMSO carrier only). As AP1 complexes act as pioneer factors that sit at regulatory elements to open up chromatin , their binding is believed to be precede and be independent of subsequent nuclear receptor binding. Replicate experiments were highly concordant (Fig. 4b). As expected, the canonical AP1 motif was highly enriched in these datasets (Fig. 4c). We next compared JunD binding events to LXRB and PPARG cistromes (Fig. 4d, e). Supporting our motif enrichment analysis, 54 % of LXRB sites overlapped with JunD at 48 h, compared with 41 % at 2 h (Fig. 4d). PPARG generated more pronounced changes, with 37 % of PPARG sites co-occurring with JunD at 2 h and with only 12 % of sites at 48 h (Fig. 4d). We also performed the inverse comparison by assessing the percentage of all JunD binding events that co-occurred with LXRB and PPARG sites and this analysis generated the same results (Fig. 4e); a 1.9-fold stronger concordance was identified at 48 h for LXRB, while a 1.4-fold higher overlap was identified at 2 h for PPARG. Using the LXRA ChIP-seq data, we identified similar enrichments patterns for AP1 motif enrichment and JunD co-occupancy as the endogenous LXRB binding profile (Additional file 1: Figure S12).
LXRs and PPARG maintain temporally distinct gene regulatory activities
Diverse genomic and functional studies support the role of RNAP2 as a marker of active, promoter-distal enhancer elements [60, 61]. To determine if these gene promoter enrichments were driven by nearby, functional nuclear receptor binding events, we performed replicate RNA polymerase II (RNAP2) ChIP-seq experiments under control culture conditions (DMSO) and after 24 and 48 h of GW3965 and rosiglitazone treatment. We first compared RNAP2 binding events to the LXR ChIP-seq datasets and divided LXR sites into those that were co-occupied by RNAP2 and sites that were devoid of RNAP2 occupancy. By integrating these distinct loci with GW3965 + T0901317 target genes that were repressed at 48 h, we identified RNAP2 co-occupied LXR sites as the principle drivers of the enrichment with down-regulated gene promoters, pointing to the functionality of these binding events (Additional file 1: Figure S18). We obtained comparable results when we used only GW3965-responsive genes (Fig. 5b, c) and lower fold change cutoffs (Additional file 1: Figure S19). Moreover, the majority of LXRA (56 %) and LXRB (79 %) binding events at 48 h were coincident with RNAP2. Identical analyses using 2-h rosiglitazone PPARG binding results and 24-h rosiglitazone RNAP2 occupancy identified a similar enrichment with RNAP2 co-occupancy at genes up-regulated by rosiglitazone at 24 h (Additional file 1: Figure S20).
To validate the functional likelihood of these binding events further, we compared our nuclear receptor ChIP-seq data with H3K27ac binding data. For all analyses, and using GW3965 + T0901317 repressed genes at 48 h, as well as 24-h rosiglitazone-activated genes, we confirmed that PPARG, LXRA, and LXRB binding sites overlapping with H3K27ac exhibited more pronounced enrichment near these drug responsive genes compared with sites devoid of H3K27ac occupancy (Additional file 1: Figures S21–S23). This observation held true using both GW3965 and high-confidence GW3965 + T0901317 gene sets for LXR comparisons and by controlling for confounding effects from negative feedback loops. Similar enrichments were observed for 24-h up-regulated genes for both LXR nuclear receptors and 48-h up-regulated genes for PPARG (Additional file 1: Figure S24). Notably, LXRA and LXRB binding events coincident with the JunD binding were also preferentially enriched near repressed genes (Additional file 1: Figure S25).
Finally, we performed quantitative analyses at LXRA and LXRB sites that co-occurred with RNAP2 at 48 h to ascertain temporal changes in RNAP2 occupancy, which may further provide mechanistic insight (Fig. 5d, e). Intriguingly, an analysis of normalized RNAP2 read depth at both LXRA and LXRB sites revealed a bimodal distribution, with a subset of sites exhibiting stronger RNAP2 enrichment prior to drug treatment and a fraction of binding events displaying greater occupancy after 48 h of GW3965 treatment. Related assessments of H3K27ac temporal patterns generated more stable distributions (Additional file 1: Figure S26), limiting this bimodal pattern to RNAP2 occupancy. Collectively, these data point to gene-activating roles for PPARG, while LXR proteins seem to harbor more dynamic regulatory activities characterized by early activating effects and late repressing functions.
LXR repression is consistent with RNAP2 promoter-proximal pausing
We next divided the LXR sites overlapping with RNAP2 sites (LXR + RNAP2) into promoter binding events and LXR + RNAP2 binding events distal to promoters (±500 bp) to determine if promoter association can explain the observed bimodal distribution of RNAP2 read enrichment (Fig. 5d, e). Notably, by separating LXR + RNAP2 co-occupied sites based on overlap with promoter regions (Fig. 6c), we were able to explain the temporal RNAP2 occupancy distribution; overlapping sites at promoters exhibited stronger RNAP2 enrichment at 48 h, while LXR binding events outside of promoters displayed diminished RNAP2 occupancy at 48 h. The University of California Santa Cruz (UCSC) Genome browser images (https://genome.ucsc.edu/) of raw ChIP-seq signal at several loci spanning 48-h, GW3965-repressed genes are provided in Additional file 1: Figures S32–S35.
The observation of a late, stronger RNAP2 promoter occupancy at LXR responsive genes is difficult to reconcile with the decrease in RNA expression during the late drug response. These effects may be mediated by post-transcriptional regulation, for example, by enhancing transcript degradation or, alternatively, by RNAP2 promoter-proximal pausing, with increased enrichment stemming from longer RNAP2 promoter residency times. To potentially distinguish between these two possibilities, we performed ChIP-seq experiments using an antibody targeting histone-3 lysine-36 trimethylation (H3K36me3), a histone modification found in gene bodies that is correlated with expression [62, 63]. We ran H3K36me3 ChIP-seq experiments under control culture conditions (DMSO) and after 48 h of GW3965 treatment. We confirmed the predictive value of H3K36me3 gene body occupancy for measuring global gene expression (rho > 0.77; Additional file 1: Figure S36) and subsequently compared H3K36me3 gene body read depth with expression of GW3965-responsive genes (Fig. 6d; Additional file 1: Figures S37–S39). For 78 % of repressed genes at 48 h, H3K36me3 displayed a concordant decrease in gene body occupancy (Fig. 6d; Additional file 1: Figure S37). This pattern was further maintained using lower fold change cutoffs (Additional file 1: Figures S38 and S39). To investigate whether H3K36me3 changes were directly related to LXRs, we evaluated dynamic H3K36me3 occupancy at all 48-h repressed genes that were enriched for nearby LXR binding events (<10 kb from promoter) and for genes that were devoid of nearby LXR binding (>100 kb from promoters). Intriguingly, genes with nearby LXRA and LXRB binding exhibited a significantly stronger decrease (Wilcoxon; LXRA p = 1.265e-8, LXRB p = 1.394e-9) in H3K36me3 gene body occupancy (Fig. 6e, f). Incorporating GW3965-responsive genes with lower fold change cutoffs led to similarly significant patterns (Wilcoxon; LXRA and LXRB p < 2.2e-16) for both proteins (Additional file 1: Figure S40). Taken together, these data are consistent with a RNAP2 promoter-pausing mechanism mediated by LXRA and LXRB nuclear receptors.
Although gene regulatory networks govern genome function and underlie diverse physiological and developmental processes, our overall understanding of how these transcriptional programs are spatiotemporally regulated remains rudimentary. Here we sought to understand the gene regulatory networks that control tumor cell metabolism and proliferation in a temporal manner by generating a thorough characterization of the gene regulatory effects of LXR and PPARG signaling. Collectively, our work illustrates how two distinct gene regulatory cascades lead to a common phenotypic outcome and further provides a better molecular understanding of how extracellular metabolic cues impact tumor cell physiology.
Multiple studies have suggested a role for LXR and PPARG metabolic nuclear receptors in tumor cell biology [64–67]. Our identification of a pronounced inhibitory effect on tumor cell proliferation and decreased glutathione levels that is consistent with increased oxidative stress agrees with previous reports [31, 32]. A recent study suggested increased reactive oxygen species via PDK4 expression and beta-oxidation of fatty acids as the likely anti-proliferative mechanism in lung cancer cells . Indeed, rosiglitazone treatment in our study also led to a significant up-regulation of PDK4 expression (fold change = 6, adjusted p value 2.09e-151). Our metabolomic interrogation further identified an increase in aerobic glycolysis, demonstrating a role for this regulatory network in the Warburg effect. Although this observation is contrary to the observed increase in oxidative stress, the assessment of metabolite concentrations cannot distinguish cause from effect and, therefore, the increase in aerobic glycolysis may reflect an indirect cellular response to oxidative stress rather than a direct effect stemming from LXR and PPARG signaling.
Through an investigation that interrogated multiple tiers of transcriptional regulation and genome organization, we demonstrated that these overlapping cell phenotypes are mediated though widely distinct signaling cascades at the molecular level. The genome-wide cellular response to rosiglitazone was rapid and supported by quick changes to metabolite levels, linear transcriptional profiles, abrupt chromatin alterations, and a short-lived burst of PPARG binding. By contrast, the more gradual changes to cellular metabolites, a stalled transcriptional response, gradual increases in the chromatin landscape, and dramatic increase in genome-wide LXRB occupancy with time highlight a more sustained GW3965 response.
Despite these differing molecular networks, both sets of nuclear receptors were strongly associated with AP1 transcriptional regulators, with PPARG displaying an early preference for AP1 motifs and co-occupancy, while LXRs exhibited a more pronounced late association with AP1 motifs and proteins. In light of the consistent association with upregulated genes, PPARG appears to behave as a transcriptional activator. This is consistent with the canonical model of type II nuclear receptor function, wherein ligand binding leads to the displacement of co-repressor complexes from nuclear receptors and subsequent recruitment of co-activator proteins . The regulatory functions of LXRA and LXRB were more complex; our genomic data points to both activating and repressive activities at early and late time points, respectively. Taken together, our results further point to LXR-mediated genetic repression through RNAP2 promoter-proximal pausing. Unlike the canonical model for PPARG activity, our results suggest the intriguing possibility that these transcriptional profiles reflect a LXR transrepression of the AP1 signaling program through protein–protein interactions at AP1 binding sites that lead to RNAP2 promoter-proximal pausing. Similar mechanisms have been reported for LXRs in macrophages on inflammatory Toll-like receptor genes [69, 70], including transrepression of the AP1 machinery . Further experimentation will be required to fully define the molecular mechanism of this putative LXR-mediated repression.
By relying on a unique, multi-tiered approach, our detailed results provide a molecular understanding of how extracellular nutrients impact cancer cell physiology at the genome level. This comprehensive analysis illustrates the complexity of genome function and structure by elucidating how common phenotypic outcomes are genetically encoded through diverse transcriptional programs. Our results finally allude to the tantalizing possibility that tumor cell growth can be altered, or even fully inhibited, through metabolic reprogramming.
FBS, fetal bovine serum; KEGG, Kyoto Encyclopedia of Genes and Genomes; LXR, liver X receptor; LXRA, liver X receptor alpha; LXRB, liver X receptor beta; PPARG, peroxisome proliferator-activated receptor gamma; RNAP2, RNA polymerase II; TF, transcription factor.
We thank Gregory Cooper and Gregory Barsh for helpful suggestions and advice, as well as Tatsuya Uechi and Amy Ridgeway for key technical support. We also thank the HudsonAlpha Genomics Services Laboratory, especially Dr. Shawn Levy and Ms. Angela Jones, for their contributions to the large-scale DNA sequencing performed for the ChIP-seq and RNA-seq experiments in this study.
This work was supported by the State Cancer Fund of Alabama (to RMM), the NIH grant R01HL117334 (DS, RMM, and MJG) and by the UAB MSTP (NIH-NIGMS 5T32GM008361-21 to RCR).
Availability of data and materials
All ChIP-seq and RNA-seq data from this study have been submitted to the NCBI Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) under the accession number GSE77039.
DS participated in conceiving the study design, coordinated and performed experiments, analyzed the results, and wrote the manuscript. RCR performed experiments, analyzed data, and edited the manuscript. BSR and TCB analyzed the results. ECD and SKM performed experiments. SJC, MJG, and JG contributed to the study design and edited the manuscript. RMM participated in conceiving the study design, coordinated the experiments, and edited the manuscript. All authors read and approved the final manuscript.
The authors declare that they have no competing interests.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Macneil LT, Walhout AJ. Gene regulatory networks and the role of robustness and stochasticity in the control of gene expression. Genome Res. 2011;21(5):645–57.View ArticlePubMedPubMed CentralGoogle Scholar
- Herschman HR. Extracellular signals, transcriptional responses and cellular specificity. Trends Biochem Sci. 1989;14(11):455–8.View ArticlePubMedGoogle Scholar
- Lopez-Maury L, Marguerat S, Bahler J. Tuning gene expression to changing environments: from rapid responses to evolutionary adaptation. Nat Rev Genet. 2008;9(8):583–93.View ArticlePubMedGoogle Scholar
- Bulger M, Groudine M. Functional and mechanistic diversity of distal transcription enhancers. Cell. 2011;144(3):327–39.View ArticlePubMedPubMed CentralGoogle Scholar
- Sakabe NJ, Savic D, Nobrega MA. Transcriptional enhancers in development and disease. Genome Biol. 2012;13(1):238.View ArticlePubMedPubMed CentralGoogle Scholar
- Levine M. Transcriptional enhancers in animal development and evolution. Curr Biol. 2010;20(17):R754–63.View ArticlePubMedPubMed CentralGoogle Scholar
- Tjian R, Maniatis T. Transcriptional activation: a complex puzzle with few easy pieces. Cell. 1994;77(1):5–8.View ArticlePubMedGoogle Scholar
- Calkin AC, Tontonoz P. Transcriptional integration of metabolism by the nuclear sterol-activated receptors LXR and FXR. Nat Rev Mol Cell Biol. 2012;13(4):213–24.PubMedPubMed CentralGoogle Scholar
- Berger J, Moller DE. The mechanisms of action of PPARs. Annu Rev Med. 2002;53:409–35.View ArticlePubMedGoogle Scholar
- Joseph SB, Laffitte BA, Patel PH, Watson MA, Matsukuma KE, Walczak R, Collins JL, Osborne TF, Tontonoz P. Direct and indirect mechanisms for regulation of fatty acid synthase gene expression by liver X receptors. J Biol Chem. 2002;277(13):11019–25.Google Scholar
- Li AC, Glass CK. PPAR- and LXR-dependent pathways controlling lipid metabolism and the development of atherosclerosis. J Lipid Res. 2004;45(12):2161–73.View ArticlePubMedGoogle Scholar
- Lehmann JM, Moore LB, Smith-Oliver TA, Wilkison WO, Willson TM, Kliewer SA. An antidiabetic thiazolidinedione is a high affinity ligand for peroxisome proliferator-activated receptor gamma (PPAR gamma). J Biol Chem. 1995;270(22):12953–6.View ArticlePubMedGoogle Scholar
- Mitro N, Mak PA, Vargas L, Godio C, Hampton E, Molteni V, Kreusch A, Saez E. The nuclear receptor LXR is a glucose sensor. Nature. 2007;445(7124):219–23.Google Scholar
- Hong C, Tontonoz P. Coordination of inflammation and metabolism by PPAR and LXR nuclear receptors. Curr Opin Genet Dev. 2008;18(5):461–7.View ArticlePubMedPubMed CentralGoogle Scholar
- Jaziri R, Lobbens S, Aubert R, Pean F, Lahmidi S, Vaxillaire M, Porchay I, Bellili N, Tichet J, Balkau B et al. The PPARG Pro12Ala polymorphism is associated with a decreased risk of developing hyperglycemia over 6 years and combines with the effect of the APM1 G-11391A single nucleotide polymorphism: the Data From an Epidemiological Study on the Insulin Resistance Syndrome (DESIR) study. Diabetes. 2006;55(4):1157–62.Google Scholar
- Rangwala SM, Lazar MA. Peroxisome proliferator-activated receptor gamma in diabetes and metabolism. Trends Pharmacol Sci. 2004;25(6):331–6.View ArticlePubMedGoogle Scholar
- Cao G, Liang Y, Broderick CL, Oldham BA, Beyer TP, Schmidt RJ, Zhang Y, Stayrook KR, Suen C, Otto KA et al. Antidiabetic action of a liver x receptor agonist mediated by inhibition of hepatic gluconeogenesis. J Biol Chem. 2003;278(2):1131–6.Google Scholar
- Steffensen KR, Gustafsson JA. Putative metabolic effects of the liver X receptor (LXR). Diabetes. 2004;53 Suppl 1:S36–42.View ArticlePubMedGoogle Scholar
- Laffitte BA, Chao LC, Li J, Walczak R, Hummasti S, Joseph SB, Castrillo A, Wilpitz DC, Mangelsdorf DJ, Collins JL et al. Activation of liver X receptor improves glucose tolerance through coordinate regulation of glucose metabolism in liver and adipose tissue. Proc Natl Acad Sci U S A. 2003;100(9):5419–24.Google Scholar
- Kliewer SA, Lenhard JM, Willson TM, Patel I, Morris DC, Lehmann JM. A prostaglandin J2 metabolite binds peroxisome proliferator-activated receptor gamma and promotes adipocyte differentiation. Cell. 1995;83(5):813–9.View ArticlePubMedGoogle Scholar
- Hu E, Kim JB, Sarraf P, Spiegelman BM. Inhibition of adipogenesis through MAP kinase-mediated phosphorylation of PPARgamma. Science. 1996;274(5295):2100–3.View ArticlePubMedGoogle Scholar
- Walczak R, Tontonoz P. PPARadigms and PPARadoxes: expanding roles for PPARgamma in the control of lipid metabolism. J Lipid Res. 2002;43(2):177–86.PubMedGoogle Scholar
- Gao M, Liu D. The liver X receptor agonist T0901317 protects mice from high fat diet-induced obesity and insulin resistance. AAPS J. 2013;15(1):258–66.View ArticlePubMedGoogle Scholar
- Kalaany NY, Gauthier KC, Zavacki AM, Mammen PP, Kitazume T, Peterson JA, Horton JD, Garry DJ, Bianco AC, Mangelsdorf DJ. LXRs regulate the balance between fat storage and oxidation. Cell Metab. 2005;1(4):231–44.Google Scholar
- Feng J, Han J, Pearce SF, Silverstein RL, Gotto Jr AM, Hajjar DP, Nicholson AC. Induction of CD36 expression by oxidized LDL and IL-4 by a common signaling pathway dependent on protein kinase C and PPAR-gamma. J Lipid Res. 2000;41(5):688–96.Google Scholar
- Collins AR, Meehan WP, Kintscher U, Jackson S, Wakino S, Noh G, Palinski W, Hsueh WA, Law RE. Troglitazone inhibits formation of early atherosclerotic lesions in diabetic and nondiabetic low density lipoprotein receptor-deficient mice. Arterioscler Thromb Vasc Biol. 2001;21(3):365–71.Google Scholar
- Chen Z, Ishibashi S, Perrey S, Osuga J, Gotoda T, Kitamine T, Tamura Y, Okazaki H, Yahagi N, Iizuka Y et al. Troglitazone inhibits atherosclerosis in apolipoprotein E-knockout mice: pleiotropic effects on CD36 expression and HDL. Arterioscler Thromb Vasc Biol. 2001;21(3):372–7.Google Scholar
- Joseph SB, McKilligin E, Pei L, Watson MA, Collins AR, Laffitte BA, Chen M, Noh G, Goodman J, Hagger GN et al. Synthetic LXR ligand inhibits the development of atherosclerosis in mice. Proc Natl Acad Sci U S A. 2002;99(11):7604–9.Google Scholar
- Bradley MN, Hong C, Chen M, Joseph SB, Wilpitz DC, Wang X, Lusis AJ, Collins A, Hseuh WA, Collins JL et al. Ligand activation of LXR beta reverses atherosclerosis and cellular cholesterol overload in mice lacking LXR alpha and apoE. J Clin Invest. 2007;117(8):2337–46.Google Scholar
- Saez E, Tontonoz P, Nelson MC, Alvarez JG, Ming UT, Baird SM, Thomazy VA, Evans RM. Activators of the nuclear receptor PPARgamma enhance colon polyp formation. Nat Med. 1998;4(9):1058–61.Google Scholar
- Srivastava N, Kollipara RK, Singh DK, Sudderth J, Hu Z, Nguyen H, Humphries CG, Carstens R, Huffman KE. Inhibition of cancer cell proliferation by PPARgamma is mediated by a metabolic switch that increases reactive oxygen species levels. Cell Metab. 2014;20(4):650–61.Google Scholar
- Vedin LL, Gustafsson JA, Steffensen KR. The oxysterol receptors LXRalpha and LXRbeta suppress proliferation in the colon. Mol Carcinog. 2013;52(11):835-44.Google Scholar
- Fukuchi J, Kokontis JM, Hiipakka RA, Chuu CP, Liao S. Antiproliferative effect of liver X receptor agonists on LNCaP human prostate cancer cells. Cancer Res. 2004;64(21):7686–9.View ArticlePubMedGoogle Scholar
- Plat J, Nichols JA, Mensink RP. Plant sterols and stanols: effects on mixed micellar composition and LXR (target gene) activation. J Lipid Res. 2005;46(11):2468–76.View ArticlePubMedGoogle Scholar
- Ronco A, De Stefani E, Boffetta P, Deneo-Pellegrini H, Mendilaharsu M, Leborgne F. Vegetables, fruits, and related nutrients and risk of breast cancer: a case-control study in Uruguay. Nutr Cancer. 1999;35(2):111–9.View ArticlePubMedGoogle Scholar
- De Stefani E, Boffetta P, Ronco AL, Brennan P, Deneo-Pellegrini H, Carzoglio JC, Mendilaharsu M. Plant sterols and risk of stomach cancer: a case-control study in Uruguay. Nutr Cancer. 2000;37(2):140–4.Google Scholar
- Mendilaharsu M, De Stefani E, Deneo-Pellegrini H, Carzoglio J, Ronco A. Phytosterols and risk of lung cancer: a case-control study in Uruguay. Lung Cancer. 1998;21(1):37–45.View ArticlePubMedGoogle Scholar
- Baskar AA, Ignacimuthu S, Paulraj GM, Al Numair KS. Chemopreventive potential of beta-Sitosterol in experimental colon cancer model--an in vitro and In vivo study. BMC Complement Altern Med. 2010;10:24.View ArticlePubMedPubMed CentralGoogle Scholar
- Chuu CP, Hiipakka RA, Kokontis JM, Fukuchi J, Chen RY, Liao S. Inhibition of tumor growth and progression of LNCaP prostate cancer cells in athymic mice by androgen and liver X receptor agonist. Cancer Res. 2006;66(13):6482–6.View ArticlePubMedGoogle Scholar
- Gabbi C, Kim HJ, Barros R, Korach-Andre M, Warner M, Gustafsson JA. Estrogen-dependent gallbladder carcinogenesis in LXRbeta-/- female mice. Proc Natl Acad Sci U S A. 2010;107(33):14763–8.View ArticlePubMedPubMed CentralGoogle Scholar
- Vedin LL, Lewandowski SA, Parini P, Gustafsson JA, Steffensen KR. The oxysterol receptor LXR inhibits proliferation of human breast cancer cells. Carcinogenesis. 2009;30(4):575–9.View ArticlePubMedGoogle Scholar
- Pommier AJ, Alves G, Viennois E, Bernard S, Communal Y, Sion B, Marceau G, Damon C, Mouzat K, Caira F et al. Liver X Receptor activation downregulates AKT survival signaling in lipid rafts and induces apoptosis of prostate cancer cells. Oncogene. 2010;29(18):2712–23.Google Scholar
- Lee JH, Gong H, Khadem S, Lu Y, Gao X, Li S, Zhang J, Xie W. Androgen deprivation by activating the liver X receptor. Endocrinology. 2008;149(8):3778–88.Google Scholar
- Kotta-Loizou I, Giaginis C, Theocharis S. The role of peroxisome proliferator-activated receptor-gamma in breast cancer. Anticancer Agents Med Chem. 2012;12(9):1025–44.View ArticlePubMedGoogle Scholar
- Pang X, Wei Y, Zhang Y, Zhang M, Lu Y, Shen P. Peroxisome proliferator-activated receptor-gamma activation inhibits hepatocellular carcinoma cell invasion by upregulating plasminogen activator inhibitor-1. Cancer Sci. 2013;104(6):672–80.View ArticlePubMedGoogle Scholar
- Girnun GD, Smith WM, Drori S, Sarraf P, Mueller E, Eng C, Nambiar P, Rosenberg DW, Bronson RT, Edelmann W et al. APC-dependent suppression of colon carcinogenesis by PPARgamma. Proc Natl Acad Sci U S A. 2002;99(21):13771–6.Google Scholar
- Jansson EA, Are A, Greicius G, Kuo IC, Kelly D, Arulampalam V, Pettersson S. The Wnt/beta-catenin signaling pathway targets PPARgamma activity in colon cancer cells. Proc Natl Acad Sci U S A. 2005;102(5):1460–5.Google Scholar
- Dunn WB, Broadhurst D, Begley P, Zelena E, Francis-McIntyre S, Anderson N, Brown M, Knowles JD, Halsall A, Haselden JN et al: Procedures for large-scale metabolic profiling of serum and plasma using gas chromatography and liquid chromatography coupled to mass spectrometry. Nat Protoc 2011, 6(7):1060–1083.View ArticlePubMedGoogle Scholar
- Gertz J, Savic D, Varley KE, Partridge EC, Safi A, Jain P, Cooper GM, Reddy TE, Crawford GE, Myers RM. Distinct properties of cell-type-specific and shared transcription factor binding sites. Mol Cell. 2013;52(1):25–36.Google Scholar
- Gertz J, Varley KE, Davis NS, Baas BJ, Goryshin IY, Vaidyanathan R, Kuersten S, Myers RM. Transposase mediated construction of RNA-seq libraries. Genome Res. 2012;22(1):134–41.Google Scholar
- Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, Nusbaum C, Myers RM, Brown M, Li W. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 2008;9(9):R137.Google Scholar
- Bailey TL, Williams N, Misleh C, Li WW. MEME: discovering and analyzing DNA and protein sequence motifs. Nucleic Acids Res. 2006;34(Web Server issue):W369–73.View ArticlePubMedPubMed CentralGoogle Scholar
- Merico D, Isserlin R, Stueker O, Emili A, Bader GD. Enrichment map: a network-based method for gene-set enrichment visualization and interpretation. PLoS One. 2010;5(11):e13984.View ArticlePubMedPubMed CentralGoogle Scholar
- Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11(10):R106.View ArticlePubMedPubMed CentralGoogle Scholar
- Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102(43):15545–50.Google Scholar
- Creyghton MP, Cheng AW, Welstead GG, Kooistra T, Carey BW, Steine EJ, Hanna J, Lodato MA, Frampton GM, Sharp PA et al. Histone H3K27ac separates active from poised enhancers and predicts developmental state. Proc Natl Acad Sci U S A. 2010;107(50):21931–6.Google Scholar
- McLean CY, Bristor D, Hiller M, Clarke SL, Schaar BT, Lowe CB, Wenger AM, Bejerano G. GREAT improves functional interpretation of cis-regulatory regions. Nat Biotechnol. 2010;28(5):495–501.Google Scholar
- Angel P, Karin M. The role of Jun, Fos and the AP-1 complex in cell-proliferation and transformation. Biochim Biophys Acta. 1991;1072(2-3):129–57.PubMedGoogle Scholar
- Biddie SC, John S, Sabo PJ, Thurman RE, Johnson TA, Schiltz RL, Miranda TB, Sung MH, Trump S, Lightman SL. Transcription factor AP1 potentiates chromatin accessibility and glucocorticoid receptor binding. Mol Cell. 2011;43(1):145–55.Google Scholar
- Savic D, Roberts BS, Carleton JB, Partridge EC, White MA, Cohen BA, Cooper GM, Gertz J, Myers RM: Promoter-distal RNA polymerase II binding discriminates active from inactive CCAAT/enhancer-binding protein beta binding sites. Genome Res. 2015;25(12):1791–800.Google Scholar
- De Santa F, Barozzi I, Mietton F, Ghisletti S, Polletti S, Tusi BK, Muller H, Ragoussis J, Wei CL, Natoli G. A large fraction of extragenic RNA pol II transcription sites overlap enhancers. PLoS Biol. 2010;8(5):e1000384.Google Scholar
- Dong X, Greven MC, Kundaje A, Djebali S, Brown JB, Cheng C, Gingeras TR, Gerstein M, Guigo R, Birney E et al. Modeling gene expression using chromatin features in various cellular contexts. Genome Biol. 2012;13(9):R53.Google Scholar
- Kolasinska-Zwierz P, Down T, Latorre I, Liu T, Liu XS, Ahringer J. Differential chromatin marking of introns and expressed exons by H3K36me3. Nat Genet. 2009;41(3):376–81.View ArticlePubMedPubMed CentralGoogle Scholar
- Bovenga F, Sabba C, Moschetta A. Uncoupling nuclear receptor LXR and cholesterol metabolism in cancer. Cell Metab. 2015;21(4):517–26.View ArticlePubMedGoogle Scholar
- Lin CY, Gustafsson JA. Targeting liver X receptors in cancer therapeutics. Nat Rev Cancer. 2015;15(4):216–24.View ArticlePubMedGoogle Scholar
- Joshi H, Pal T, Ramaa CS. A new dawn for the use of thiazolidinediones in cancer therapy. Expert Opin Investig Drugs. 2014;23(4):501–10.View ArticlePubMedGoogle Scholar
- Zhang Z, Xu Y, Xu Q, Hou Y. PPARgamma against tumors by different signaling pathways. Onkologie. 2013;36(10):598–601.View ArticlePubMedGoogle Scholar
- Sever R, Glass CK. Signaling by nuclear receptors. Cold Spring Harb Perspect Biol. 2013;5(3):a016709.View ArticlePubMedPubMed CentralGoogle Scholar
- Ghisletti S, Huang W, Ogawa S, Pascual G, Lin ME, Willson TM, Rosenfeld MG, Glass CK. Parallel SUMOylation-dependent pathways mediate gene- and signal-specific transrepression by LXRs and PPARgamma. Mol Cell. 2007;25(1):57–70.Google Scholar
- Huang W, Ghisletti S, Saijo K, Gandhi M, Aouadi M, Tesz GJ, Zhang DX, Yao J, Czech MP, Goode BL et al. Coronin 2A mediates actin-dependent de-repression of inflammatory response genes. Nature. 2011;470(7334):414–8.Google Scholar
- Ogawa S, Lozach J, Jepsen K, Sawka-Verhelle D, Perissi V, Sasik R, Rose DW, Johnson RS, Rosenfeld MG, Glass CK. A nuclear receptor corepressor transcriptional checkpoint controlling activator protein 1-dependent gene networks required for macrophage activation. Proc Natl Acad Sci U S A. 2004;101(40):14461–6.Google Scholar