Genome-wide changes accompanying knockdown of fatty acid synthase in breast cancer

Changes in gene expression resulting from knockdown of FAS

MDA-MB-435 mammary carcinoma cells were selected as the model for defining the FAS knockdown signature. For our experiment, four independent siRNA duplexes targeting FAS (FAS #1-#4) were chosen based on there ability to knockdown the enzyme and to induce tumor cell apoptosis after 72 h [17]. Inhibition of FAS by each duplex was verified using 25 nM of siRNA as demonstrated by a decrease in FAS mRNA, protein, and fatty acid biosynthesis after 48 h relative to non-silencing control siRNA (Figure 1a–c). The non-silencing control siRNA was selected based on minimal cross reactivity with known targets and had no impact on FAS expression or activity when compared against Lipofectamine 2000 transfection alone (data not shown). Gene expression profiles were examined on two separate occasions following transfection with FAS siRNA. Treatment times (12, 24, 36 and 48 h) were chosen to capture early and late gene changes associated with the block in cell cycle progression and the advent of apoptosis occurring in response to knockdown of FAS. Abrogation of FAS (>70%) was verified within 12 h of transfection and persisted throughout the 48 h experiment. Significant changes in cell viability were not observed 48 h post-transfection indicating cells were viable during the time of gene analysis (data not shown).

The first approach we used in analyzing the effect of knockdown of FAS on gene expression was to identify target genes modified by at least 3 of the siRNA duplexes. The expression signature of each FAS siRNA was determined by identifying genes significantly changed over the time course using MB statistics and two way ANOVA analysis. For each siRNA duplex only genes with a p-value ≤ 0.05 and at least 1.2 fold changes in both biological replicates at a given time point were maintained for expression analysis, and are published as supporting information online (see Additional Files 1, 2, 3, and 4). The FAS knockdown signature was defined as the overlap of significant gene changes occurring in response to at least 3 of the siRNA duplexes with an average gene expression change 1.5 fold that observed in controls. This allowed us to identify gene changes specifically associated with knockdown of FAS, while at the same time eliminate potential off-target effects unique to the individual siRNA sequences. Using this approach, we identified 279 genes whose expression changed in response to knockdown of FAS (169 genes up-regulated and 110 down-regulated; see Additional File 5). Alterations in gene expression occurred as early as 12 h after knockdown of FAS with the majority of genes affected by 24 h (Figure 1d). These results show that inactivation of FAS has a profound affect on gene transcription.

Several of the genes significantly altered by knockdown of FAS are known to play a role in the regulation of lipid metabolism. Inactivation of FAS increased the expression of INSIG1, and the leptin receptors, LEPR and LEPROTL1, which block transcription, proteolytic cleavage and transcriptional activation of the sterol regulatory element-binding protein (SREBP) family of transcription factors that promote lipid biosynthesis (Table 1) [19, 20]. Additionally, loss of FAS led to down-regulation of ID2 and ID3, which are dominant negative helix-loop-helix (HLH) proteins that bind SREBP 1c and functionally repress FAS promoter activity [21]. ID2 and ID3 also repress expression of the cyclin-dependent kinase inhibitor p21^cip1 [22, 23], thus, as expected the suppression of these genes in response to knockdown of FAS was coupled with up-regulation of the p21^{cip 1}gene (Table 1). These data provide evidence of a direct link between regulation of FAS and control of cell cycle progression.

Another important biological category affected by knockdown of FAS is regulation of protein ubiquitination (Table 1). We found that knockdown of FAS altered the expression of several E2 ubiquitin conjugation enzymes (2 up-regulated) and E3 ubiquitin ligases (4 down-regulated and 6 up-regulated) which function to target proteins for degradation by the proteosome [24]. In addition, we identified up-regulation of 2 splice variants of UBE2V1, which are catalytically inactive E2 enzymes [25]. UBE2V1 variants result from the co-transcription of UBE2V1 with the neighboring upstream gene, Kua, which shares sequence consensus to fatty acid hydrolases (variant 2, Kua-UEV fusion gene [GI_40806191-I]; variant 4, Kua gene [GI_40806192-I]). This targeting of ubiquitination enzymes indicates that the knockdown of FAS affects the tumor cell proteome not only through changes in transcription, but also on the post-translational level.

Knockdown of FAS leads to widespread changes in pathways involved in metabolism

The second approach we used to analyze the gene array results was to identify gene pathways coordinately up-regulated or down-regulated by knockdown of FAS using Gene Set Enrichment Analysis (GSEA). GSEA is a unique application that coordinates small changes in gene expression from a large number of functionally related genes in order to identify pathways that are significantly changed [26]. As expected, fatty acid metabolism (NES = -1.808, NOM p-val = 0.002, FDR q-val = 0.085) was significantly suppressed by knockdown of FAS (Figure 2). Fatty acid elongation (NES = -1.537, NOM p-val = 0.037, FDR q-val = 0.394) was also reduced, which is consistent with the loss of available palmitate.

The effects of knockdown of FAS on the expression of metabolic genes extended to pathways other than fatty acid metabolism (Figure 2). Inactivation of FAS down-regulated both the glycolysis/gluconeogenesis (NES = -1.560, NOM p-val = 0.029, FDR q-val = 0.377) and krebs-TCA cycle (NES = -1.737, NOM p-val = 0.014, FDR q-val = 0.094) pathways. Down-regulation of the krebs-TCA cycle was consistent with an overall suppression in mitochondrial genes involved in energy metabolism and oxidative phosphorylation (GO 005739; NES = -1.515, NOM p-val = 0.039, FDR q-val = 0.365). Furthermore, this down-regulation in genes regulating glucose utilization was coupled with a reduction in insulin signaling (NES = -1.635, NOM p-val = 0.005, FDR q-val = 0.280). This suggests that inhibition of FAS may restrict the ability of tumor cells to produce the energy necessary to thrive.

Knockdown of FAS up-regulates pathways involved in cell cycle arrest and apoptosis

The anti-tumorigenic effects of knockdown of FAS were traced to an up-regulation in genes that modify tumor proliferation (Figure 3). GSEA analysis showed that knockdown of FAS up-regulated genes involved in DNA damage signaling (NES = 1.809, NOM p-val = 0.000, FDR q-val = 0.068) and cell cycle arrest (NES = 1.504, NOM p-val = 0.036, FDR q-val = 0.163). Induction of cell cycle arrest was coupled with up-regulation of the p27 (NES = 1.582, NOM p-val = 0.009, FDR q-val = 0.128) and Rb (NES = 1.581, NOM p-val = 0.010, FDR q-val = 0.126) pathways consistent with our previous finding that FAS regulates the G₁ checkpoint through these pathways [9]. Additionally, knockdown of FAS elevated the expression of several cyclin-dependent kinase inhibitors (p16INK4 (CDKN2A), p15INK4 (CDKN2B) and p21^{Cip 1}) which function as negative regulators of the Rb pathway [27], demonstrating FAS induces a high degree of regulation over the G₁/S transition (Figure 3). Results also showed that knockdown of FAS elevated genes involved in the G₂ pathway including the checkpoint regulators CHEK2, WEE1 and PLK1 (Table 1; Figure 3). These data indicate that FAS controls tumor proliferation by regulating several aspects of the cell division cycle.

The apoptotic effects of knockdown of FAS were recently shown to be blocked by co-treatment with the pan-caspase inhibitor z-VAD-fmk [28], suggesting a role for caspases in mediating FAS induced tumor cell death. Here, we found up-regulation of caspase 7 and caspase 8 (receptor-mediated apoptosis) in response to knockdown of FAS (Figure 3). This was accompanied by enhanced expression of pathways and genes that activate caspase 8 such as tumor necrosis factor (tnfr1, NES = 1.408, NOM p-val = 0.038, FDR q-val = 0.230; tnfr2, NES = 1.788, NOM p-val = 0.002, FDR q-val = 0.070; Figure 3). GSEA results also demonstrated compensatory up-regulation of NF-κB (NF-κB Induced, NES = 1.5, NOM p-val = 0.022, FDR q-val = 0.166) and the MAP kinase (NES = 1.539, NOM p-val = 0.016, FDR q-val = 0.133) survival pathways in line with previous reports (see Additional File 7; [14]).

Cell line and culture conditions

The MDA-MB-435 mammary tumor cell line was obtained from Janet Price at the University of Texas Southwestern. MDA-MB-435 cells were maintained in Minimum Essential Medium Eagle (MEM) with Earl's Salts (Mediatech, Inc., Herndon, VA, USA) supplemented with 10% fetal bovine serum (Irvine Scientific, Santa Ana, CA, USA), 2 mM L-glutamine (Invitrogen Life Technologies, Inc., Carlsbad, CA, USA), MEM vitamins (Invitrogen Life Technologies, Inc.), nonessential amino acids (Mediatech, Inc.) and antibiotics (Omega Scientific, Inc., Tarzana, CA, USA). Cells were grown at 37°C under a humidified, 5% CO₂ atmosphere.

FAS gene silencing using siRNA

Four individual FAS siRNA sequences corresponding to 5'-GAGCGUAUCUGUGAGAAACUU-3' (nucleotides 6241–6259; FAS#1), 5'-GACGAGAGCACCUUUGAUGUU-3' (nucleotides 1741–1749; FAS#2), 5'-UGACAUCGUCCAUUCGUUUUU-3' (nucleotides 1758–1776; FAS#3), 5'-UGACAUCGUCCAUUCGUUUUU-3' (nucleotides 6236–6254; FAS#4) were designed and synthesized by Dharmacon (Lafayette, CO, USA). MDA-MB-435 cells were plated at 7.81 × 10³ cells/cm² in 10-cm² plates and grown for 24 h prior to transfection with 25 nM FAS#1, FAS#2, FAS#3, FAS#4 or non-silencing control duplex #2 (D-001210-02) siRNA in Opti-MEM medium (Invitrogen) using Lipofectamine 2000 reagent (Invitrogen). After 5 h, transfection medium was replaced with normal culture medium and cells were grown for 12–48 h post-transfection.

Real-time quantitative PCR

Total RNA was isolated from adherent cells using the Qiagen RNeasy kit (Qiagen, Inc, Valencia, CA, USA). Genomic DNA contamination was removed using RNase-Free DNase (Qiagen). RNA purity was assessed by A₂₆₀/A₂₈₀ absorption and RNA integrity was verified by agarose gel electrophoresis. cDNA was synthesized from 4 μg RNA using the Superscript III First-Strand Synthesis System (Invitrogen). Real-time quantitative PCR reactions were performed using Power SYBR^® Green PCR Master Mix (Applied Biosystems, Foster City, CA, USA) on a Stratagene Mx3000p QPCR System (La Jolla, CA, USA). cDNA was amplified using forward and reverse primers for FAS (5'-AACTCCATGTTTGGTGTTTG-3' and 5'-CACATGCGGTTTAATTGTG-3') and normalized to the housekeeping gene P0 (5'CAAGACTGGAGACAAAGTGG-3' and 5'AATCTGCAGACAGACACTGG-3'). The formation of a single PCR product was verified using automated melting curve analysis.

Western blot analysis

Adherent cell populations were lysed in 2× SDS sample buffer, separated by SDS-PAGE and transferred to nitrocellulose as previously described [9]. Membranes were probed for FAS and β-tubulin protein expressions using anti-FAS (BD Transduction Laboratories, San Jose, CA, USA) and anti-β-tubulin (Upstate Biotechnology, Lake Placid, NY, USA) antibodies. Immunoreactivity was detected using anti-mouse IgG conjugated peroxidase and visualized by enhanced chemiluminescence.

FAS activity

FAS activity was determined by measuring the incorporation of [¹⁴C]malonyl-CoA into cellular fatty acids as previously described [9]. Briefly, MDA-MB-435 cells were lysed in 20 mM Tris (pH 7.5), 1 mM dithiothreitol, and 1 mM EDTA by sonication. Insoluble material was removed by centrifugation (14,000 rpm) for 15 min at 4°C. The resulting lysates (70 ug) were incubated in reaction buffer containing 500 μM NADPH, 166.6 μM acetyl CoA, 100 mM KCl and 0.4 μCi [¹⁴C]malonyl-CoA (GE Healthcare, Piscataway, NJ, USA) for 25 min at 37°C. Reactions were chased with 25 nM cold malonyl-CoA for 15 min and terminated by the addition of chloroform:methanol (1:1). The chloroform extract was dried under N₂ and extracted with water-saturated butanol. The butanol extract was evaporated under N₂, and labeled fatty acids were detected using a Beckman model LS 5000CE Scinillation Counter (Beckman Coulter, Inc., Fullerton, CA, USA).

BeadArray gene expression

MDA-MB-435 cells were transfected with siRNA on two separate occasions to generate independent measurements for each treatment, and monitored for whole genome expression at 12 h intervals over a 48 h time period. Total RNA was isolated and examined for purity/integrity as described above. Labeled cRNA was prepared from 500 ng RNA using the Illumina^® RNA Amplification Kit from Ambion (Austin, TX, USA). The labeled cRNA (1500 ng) was hybridized overnight at 55°C to the Sentrix^® HumanRef-6 Expression BeadChip (>46,000 gene transcripts; Illumina, San Diego, CA, USA) according to the manufacturer's instructions. BeadChips were subsequently washed and developed with fluorolink streptavidin-Cy3 (GE Healthcare). BeadChips were scanned with an Illumina BeadArray Reader and hybridization efficiency was monitored using BeadStudio software (Illumina) to measure internal controls built into the Illumina system.

Expression data analysis

Expression data was filtered to identify genes expressed on the BeadChip at >0.99 confidence (16,585 genes), and normalized per chip using non-linear Normalize.quantiles or VSN normalization to remove array to array variability (Bioconductor Project [43]). The efficiency of both normalization methods were validated by measuring gene to gene correlations across BeadChips according to the methods of Ploner et al. [44] using R 2.0.1 software (R Development Core Team, Vienna, Austria).

Differentially expressed genes that changed over time in response to each of the four FAS siRNA duplex were identified using two independent methods. In the first method, per chip normalized Normalize.quantiles data was imported into the GeneSpring GX 7.3.1 software package (Agilent Technologies, Palo Alto, CA, USA), and per gene normalized to the appropriate biological control for each time point. Genes that significantly changed over time in response to each FAS siRNA duplex were compared against control using a two-way ANOVA with a Benjamini & Hochberg FDR of 0.05. In the second method, VSN normalized data was imported into the Bioconductor Timecourse package [43]. Genes regulated over time by each FAS siRNA duplex in both biological replicates were determined using mb. long() statistics, with an mb>0 used as criterion for the identification of significantly regulated genes. Genes identified as significant using these two complementary approaches were combined, and genes with a p-value ≤ 0.05 and with 1.2 fold changes in both the biological replicates of each duplex at a given time point, were maintained for expression analysis. Genes identified as significantly modified by each of the 4 FAS siRNA duplexes were then overlapped using venn diagrams to identify a set of genes modified by 3 or more of the siRNA treatments. Expression profiles of the core gene list were averaged, and genes modified by 1.5 fold were visually examined and ordered by hierarchical clustering using the Pearson Correlation similarity measurement (GeneSpring).

Gene set enrichment analysis

For pathway analysis, data from the 16,585 genes expressed in MDA-MB-435 cells was exported from GeneSpring, filtered for duplicate symbols and analyzed using GSEA software (Broad Institute) according to published methods [26]. Briefly, data was overlapped on 522 gene sets (s2.symbols.gmt) downloaded with the GSEA package and measured for the enrichment of genes at the top or bottom of the gene list to determine their correlation with the gene set's phenotype. The GSEA parameters used included: metric for ranking genes, signal2noise; enrichment statistic, weighted; permutation type, phenotype; permutation number, 1000; and gene set size restrictions, 4 minimum, 500 maximum. Gene sets significantly modified by FAS siRNA treatment were identified using a nominal p-value < 0.05 and a multiple hypothesis testing FDR < 0.25. NES represents the enrichment of genes in the designated GSEA gene set, ranked according to the overrepresentation of genes at the top or bottom of the list, normalized to gene set size.