Breast Cancer Res Treat (2015) 150:321–334 DOI 10.1007/s10549-015-3327-1
PRECLINICAL STUDY
Integrative transcriptome-wide analyses reveal critical HER2regulated mRNAs and lincRNAs in HER2+ breast cancer Callie R. Merry • Sarah McMahon • Cheryl L. Thompson • Kristy L. S. Miskimen Lyndsay N. Harris • Ahmad M. Khalil
•
Received: 7 October 2014 / Accepted: 28 February 2015 / Published online: 8 March 2015 Ó Springer Science+Business Media New York 2015
Abstract Breast cancer is a major health problem affecting millions of women worldwide. Over 200,000 new cases are diagnosed annually in the USA, with approximately 40,000 of these cases resulting in death. HER2-positive (HER2?) breast tumors, representing 20–30 % of early-stage breast cancer diagnoses, are characterized by the amplification of the HER2 gene. However, the critical genes and pathways that become affected by HER2 amplification in humans are yet to be specifically identified. Furthermore, it is yet to be determined if HER2 amplification also affects the expression of long intervening non-coding (linc)RNAs, which are involved in the epigenetic regulation of gene expression. We examined changes in gene expression by next generation RNA sequencing in human tumors pre- and post- HER2 inhibition by trastuzumab in vivo, and changes in gene expression in response to HER2 knock down in cell culture models. We integrated our results with gene expression analysis of
HER2? tumors vs matched normal tissue from The Cancer Genome Atlas. The integrative analyses of these datasets led to the identification of a small set of mRNAs, and the associated biological pathways that become deregulated by HER2 amplification. Furthermore, our analyses identified three lincRNAs that become deregulated in response to HER2 amplification both in vitro and in vivo. Our results should provide the foundation for functional studies of these candidate mRNAs and lincRNAs to further our understanding of how HER2 amplification results in tumorigenesis. Also, the identified lincRNAs could potentially open the door for future RNA-based biomarkers and therapeutics in HER2? breast cancer. Keywords Breast cancer HER2 Gene expression lincRNAs RNA-seq Genome-wide
Introduction Electronic supplementary material The online version of this article (doi:10.1007/s10549-015-3327-1) contains supplementary material, which is available to authorized users. C. R. Merry S. McMahon A. M. Khalil (&) Departments of Genetics and Genome Sciences, Case Western Reserve University, Case Western Reserve University, Cleveland, OH 44106, USA e-mail: [email protected] C. R. Merry A. M. Khalil Departments of Biochemistry, Case Western Reserve University, Case Western Reserve University, Cleveland, OH 44106, USA
Breast cancer is the leading cause of cancer among women, with over 200,000 new cases diagnosed in the United C. L. Thompson Departments of Family Medicine and Community Health, Case Western Reserve University, Case Western Reserve University, Cleveland, OH 44106, USA L. N. Harris Departments of Medicine, Case Western Reserve University, Case Western Reserve University, Cleveland, OH 44106, USA
C. L. Thompson K. L. S. Miskimen L. N. Harris A. M. Khalil Case Comprehensive Cancer Center, Case Western Reserve University, Case Western Reserve University, Cleveland, OH 44106, USA
123
322
States each year [1]. Breast cancer is composed of several subtypes [2], with 20–30 % of early-stage breast cancers characterized by the amplification of the human epidermal growth factor receptor 2 (HER2) gene, also known as HER2/neu and ERBB2 [3]. In normal cells, HER2 signaling is initiated through heterodimerization of HER2 with other ERBB family members in response to extracellular signals. Dimerization and subsequent autophosphorylation through the tyrosine kinase domain initiate a signaling cascade that ultimately leads to changes in gene expression patterns that regulate cell proliferation and growth. However, how HER2 gene amplification in breast cancer leads to uncontrolled cell proliferation is yet to be fully understood [3]. We postulated that identifying the key genes (both coding and non-coding) that become deregulated in response to HER2 amplification could provide important insights into how HER2 amplification affects cell proliferation. Furthermore, identifying novel noncoding genes, such as lincRNAs, that become deregulated when HER2 becomes amplified may potentially provide clues into global changes in gene expression patterns observed in HER2? breast cancer (see below). The human genome encodes, in addition to mRNAs and small non-coding RNAs, over 8,300 long intervening noncoding RNAs (lincRNAs) [4–6]. Recent studies have demonstrated important roles for lincRNAs in embryonic and organ development in vivo [7–9]. Additionally, gene expression studies revealed that each lincRNA can regulate numerous mRNA genes, mostly at genomic sites far away from a lincRNA site of transcription, and thus, lincRNAs are thought to exert their effects in trans [5, 10–14]. The impact of lincRNAs on gene expression patterns has implicated lincRNAs in a wide range of biological functions including dosage compensation, genomic imprinting, alternative splicing, nuclear organization, and regulation of mRNA translation [15–20]. Intriguingly, lincRNAs utilize various mechanisms to exert their effects in the cell. These mechanisms include acting as scaffolds and guides for chromatin-modifying complexes to specific regions of the genome, serving as decoys that regulate transcription factors binding to specific DNA sequences, and as microRNA ‘‘sponges’’ that regulate miRNA:mRNA interactions [5, 11–14, 21–24]. Recently, there has been growing interest in unraveling the roles of lincRNAs in human disease with the hope that these novel transcripts can be used as diagnostic biomarkers and/or therapeutic targets [18, 25–27]. In addition, there is a significant interest in studying the functional roles of lincRNAs in cancer since numerous lincRNAs have been shown to be dysregulated across multiple cancer types [8, 28, 29]. These studies have also shown that the expression of some lincRNAs correlates with clinical parameters, such as overall patient prognosis
123
Breast Cancer Res Treat (2015) 150:321–334
and metastasis [14, 25, 30, 31]. For example, the lincRNA HOTAIR is highly up-regulated in breast cancer, hepatocellular carcinoma, and colorectal cancers [14, 30, 31]. Additionally, in vitro and in vivo functional studies have implicated HOTAIR in promoting cancer metastasis [14]. The lincRNA linc-p21 is transcriptionally activated by the tumor suppressor p53 and subsequently represses a subset of p53 gene targets in trans through binding of hnRNP-K, thereby facilitating the p53 response in cells [13]. A final example focuses on the X inactive specific transcript (Xist), which has been known for decades to be the key transcript for X chromosome inactivation (Xi) in mammalian females. Recently, it has been shown that genetic deletion of Xist in hematopoietic stem cells results in a variety of hematological cancers and premature death in mice [20]. Taken together, these selected examples, as well as many others, demonstrate important roles of lincRNAs in tumorigenesis and metastasis, and their potential utilization as biomarkers and/or therapeutic targets. To further our understanding of the effects of HER2 signaling on gene expression of both mRNAs as well as novel non-coding genes such as lincRNAs, we utilized a combined in vitro and in vivo transcriptomic approach to pinpoint critical downstream genes. The expression of these identified mRNAs and lincRNAs was subsequently examined in RNA-seq datasets obtained from The Cancer Genome Atlas (TCGA) project to further refine our top candidates. Our results led to the identification of potentially key mRNAs and lincRNAs that may contribute to HER2? breast cancer.
Results Inhibition of HER2 by trastuzumab in vivo affects the expression of both mRNAs and lincRNAs Previous studies have shown that HER2 inhibition in cell culture models of HER2? breast cancer leads to dramatic changes in the expression of mRNAs. However, it is not currently known which of these mRNAs are affected as a result of HER2 inhibition in tumors in vivo, as well as the effects of HER2 inhibition on regulatory non-coding RNAs such as lincRNAs. To identify the mRNAs and lincRNAs that become affected by HER2 inhibition in vivo, we analyzed RNA-seq data from a clinical trial that we originally designed to predict benefit from the HER2 inhibitor, trastuzumab (HerceptinÒ) (Clinical Trials.gov ID NCT00617942, also see ‘‘Materials and methods’’ section). Trastuzumab is a monoclonal antibody drug designed to target the HER2 receptor by potentially inhibiting its heterodimerization [32]. The clinical trial accrued 80 patients of which 50 pairs of tumors could be biopsied pre-
Breast Cancer Res Treat (2015) 150:321–334
trastuzumab and post-trastuzumab (post-one dose, *10–14 days). All patients continued receiving a combination of chemotherapy and trastuzumab for 4 months, and a subset of thirteen tumor pairs representing the extremes of response to treatment was subjected to total RNA isolation and RNA-seq analysis (see ‘‘Materials and methods’’ section). Eleven out of these thirteen patients were identified as responders to trastuzumab as measured by pathological complete response (pCR) (Supporting file 1). We reasoned that since trastuzumab was an effective therapeutic drug in these eleven HER2? breast cancer patients, changes in gene expression pre- and post-one dose of trastuzumab would reflect mRNAs and lincRNAs regulated through HER2. Therefore, we identified differentially expressed mRNAs and lincRNAs pre- versus postone dose of trastuzumab treatment in our RNA-seq dataset of these eleven patients that were responsive to treatment as measured by pCR. We calculated FPKM (fragments per kilobase of exon per million fragments mapped) values of each known mRNA and lincRNA in the human genome in each RNAseq sample (22 total: 11 pre- vs. 11 post-trastuzumab). We calculated an average FPKM of each mRNA and lincRNA in pre- versus post-trastuzumab, and subsequently calculated a fold change of Post/Pre for each mRNA and lincRNA. We identified 228 mRNAs and 28 lincRNAs to be differentially expressed (Fig. 1a, b, and Supporting Files 2 and 3). We also performed pathway analysis on differentially expressed mRNAs and identified five cancer-related pathways to be significantly affected (p \ 0.05, p values ranged from 0.006 to 0.00004) (Fig. 1c). In summary, we identified a set of mRNAs and lincRNAs that become affected when HER2 signaling is inhibited in vivo. However, because of genetic and environmental differences of these 11 patients, we observed heterogeneity in the level of differential expression of both mRNAs and lincRNAs. Thus, to further refine our list of key mRNAs and lincRNAs that are critical components of the HER2 pathway, we decided to modulate HER2 in a cell culture model (see below). HER2 signaling affects the expression of mRNAs as well as lincRNAs in cell culture model of HER2? breast cancer We began our studies in BT474 cells, which are HER2? breast cancer cells that have been utilized extensively as a cell culture model to study HER2? breast cancer. To identify mRNAs and lincRNAs that are downstream of HER2 signaling in BT474 cells, we utilized a loss-offunction approach in which we depleted HER2 using validated siRNAs and examined changes in gene expression using RNA-seq. We transfected BT474 cells with
323
siRNAs targeting HER2, and simultaneously transfected the same number of cells with negative control siRNAs. We initially examined HER2 mRNA levels using RTqPCR at 48 and 72 h. At these time points, we found that the siRNAs were effective at knocking down HER2 by 80 and 87 %, respectively (Fig. 2a). Subsequently, we performed a new round of transfections at 24, 48, and 72 h, and examined HER2 protein levels by Western blot analysis. At 24 h, there was a modest reduction in HER2 protein levels; however, a significant reduction in HER2 protein levels in comparison to cells transfected with negative control siRNAs was achieved 48 h post-transfection (Fig. 2b). Based on these gene expression analyses of HER2 knock down with siRNAs, the 48 h post-transfection time point was chosen as optimal to sufficiently deplete HER2 at both the mRNA and protein levels, while still capturing some of the early changes in mRNAs and lincRNAs expression due to loss of HER2. To that end, equal numbers of BT474 cells were transfected with either HER2 siRNAs or negative control siRNAs in three biological replicates, and at 48 h post-transfection total RNA was isolated, quantified, and subjected to RNA-seq (see ‘‘Materials and methods’’ section). We found 1,015 mRNAs (303 up-regulated and 712 down-regulated) and 167 lincRNAs (139 up-regulated and 28 down-regulated) to be differentially expressed by Ctwofold in the HER2-depleted BT474 cells in comparison to cells transfected with negative control siRNAs (Fig. 2c, d, and Supporting Files 4 and 5). To place our findings into context, we identified a previous study that utilized mRNA microarrays to examine the effects of HER2 inhibition by trastuzumab on mRNAs in BT474 and SKBR3 cells [33]. This study reported sixteen mRNA genes that are significantly down-regulated in response to HER2 inhibition by trastuzumab. We found that eight of these sixteen mRNA genes were also significantly down-regulated in our current study. We designed primers for those mRNA genes, and confirmed in an independent set of knockdown experiments that their expression is altered in response to HER2 depletion in BT474 breast tumor cells by RT-qPCR (Fig. 2e). Furthermore, pathway analysis of differentially expressed mRNAs revealed numerous affected pathways post-HER2 depletion in BT474 cells (Supporting File 6). Since no previous studies of lincRNAs modulated in response to HER2 depletion or inhibition are available, we selected five lincRNAs from our RNA-seq data for validation by qPCR in an independent knockdown experiment in BT474 cells. We selected these lincRNAs based on fold changes and p values and we included both up- and down-regulated lincRNAs in response to HER2 depletion. Four out of the five lincRNAs showed statistically significant up- or down-regulation in response to HER2 depletion by RT-qPCR similar to what we observed
123
324
Breast Cancer Res Treat (2015) 150:321–334
A
B
C
Fig. 1 Inhibition of HER2 in vivo by trastuzumab results in significant changes in the expression of mRNAs and lincRNAs in HER2? breast cancer patients. Gene expression analysis by next generation RNA sequencing (RNA-seq) was performed on RNA isolated from HER2? tumors pre- and post-one dose of trastuzumab. Expression of genes was calculated as FPKM values and differentially expressed genes between pre- and post-treatment conditions were
identified. Differentially expressed mRNAs (a) and lincRNAs (b) in all 22 RNA-seq samples (11 pre- vs 11 post-trastuzumab) are represented by heatmaps; c Over enrichment analysis of differentially expressed mRNAs pre- versus post-trastuzumab reveals several key pathways affected in vivo when HER2 signaling is inhibited. A number of genes and p values are shown for each pathway
by RNA-seq analysis (Fig. 2f). To further confirm that these lincRNAs are downstream of the HER2 pathway, and that changes in their expression are not due to off-target effects of siRNAs, we utilized a second independent siRNA against HER2. First, we confirmed that this second siRNA is effective at knocking down HER2 protein levels by Western blot analysis (Fig. 3a). Next, we performed qPCR analysis on the same 5 lincRNAs (see Fig. 2f), and found that four out of the five lincRNAs show a similar response to knocking down HER2 with the second siRNA,
but only three lincRNAs pass a p value of \0.05 (Fig. 3b). Also, the expression of lincMCl1-1, which was not affected with first siRNA, was responsive to the second siRNA, similar to what we observed by RNA-seq. These experiments demonstrate that HER2 depletion in BT474 cells affects the expression of both mRNAs and lincRNAs; however, some variability is observed due to off-target effects of siRNA-mediated depletion. Thus, to overcome these limitations, we identified mRNAs and lincRNAs that are affected by both HER2 inhibition by trastuzumab
123
Breast Cancer Res Treat (2015) 150:321–334 Fig. 2 HER2 affects the expression of mRNAs and lincRNAs in BT474 cells. a Real-time qPCR analysis of HER2 mRNA levels in BT474 cells transfected with either control siRNAs or HER2 siRNAs (15 nM final concentration). We observed a significant knockdown of HER2 mRNA levels at 48 h postsiRNA transfections. b Western blot analysis of HER2 protein levels in BT474 cells transfected with control siRNAs versus cells transfected with HER2 siRNAs at 48 h posttransfections. A significant reduction in HER2 protein levels is observed. Actin was used as a loading control; c, d. Gene expression analysis by RNA-seq led to the identification of differentially expressed mRNAs and lincRNAs in BT474 cells treated with either control or HER2 siRNAs. Heatmaps of differentially expressed mRNAs and lincRNAs between BT474 cells transfected with control siRNAs vs HER2 siRNAs are shown. e, f. Validation of RNAseq data by RT-qPCR of eight mRNA genes and five lincRNAs in cells transfected with control siRNAs versus HER2 siRNAs. For both mRNAs and lincRNAs, we observed similar patterns in our original knockdown experiments analyzed by RNA-seq and subsequent knockdown experiments analyzed by RTqPCR
325
A
C
E
in vivo and HER2 depletion in BT474 cells by siRNAs (see below). To identify mRNAs and lincRNAs that are affected by both HER2 inhibition by trastuzumab in tumors and HER2 depletion in cell culture, we intersected differentially expressed mRNAs and lincRNAs that we identified postHER2 knockdown in BT474 cells with differentially
B
D
F
expressed mRNAs and lincRNAs that we identified in response to HER2 inhibition by trastuzumab in vivo. We found 44 mRNAs and 3 lincRNAs to be common between the two datasets (Fig. 4a, b, and Supporting file 7 and 8). Some of the 44 mRNA genes identified are key components of the PLK1 signaling pathway and E2F transcription factor network (Fig. 4c), which are key pathways affected
123
326
A
Breast Cancer Res Treat (2015) 150:321–334
role of these genes in HER2? breast cancer, we next examined their expression in a cohort of HER2? breast tumors and their matched normal control tissue in the Cancer Genome Atlas (TCGA) database. Thousands of mRNAs and hundreds of lincRNAs are differentially expressed in HER2? breast cancer tumors
B
Fig. 3 Validation of HER2-regulated lincRNAs in BT474 cells using a second independent siRNA. a We designed and confirmed the knockdown of HER2 protein using a second independent siRNA (distinct from first siRNA used for knock down experiments and subsequent RNA-seq analysis). Western blot analysis of HER2 protein levels in BT474 cells treated with siRNA#2 alone or in combination with first siRNA (siRNA pool) demonstrates the effectiveness of this second siRNA in knocking down HER2 protein in comparison to negative control siRNAs, b We examined the expression of the same 5 lincRNAs in Fig. 2f post-HER2 knock down in BT474 cells with a second independent siRNA by qPCR analysis. We found that HER2 knock down with a second independent siRNA also affects the expression of some of these lincRNAs. These results demonstrate that some of the lincRNAs identified in our cell culture studies are downstream of the HER2 signaling pathway and some are affected due to off-target effects of siRNAs
in HER2? tumors (see next ‘‘Result’’ section). In summary, we have identified a small set of mRNAs and lincRNAs that are affected in response to HER2 inhibition/ depletion in vivo and in cell culture suggesting an important role of not only mRNAs but also lincRNAs in HER2? breast cancer. To gain further insights into the potential
123
Using transcriptomic analyses pre- versus post-HER2 knockdown in cell culture and pre- versus post-HER2 inhibition in tumors, we are able to identify a small set of lincRNAs and mRNAs that are putative members of the HER2 regulatory network. To further validate our observations, we turned to RNA-seq data from The Cancer Genome Atlas (TCGA) project to determine the expression patterns of these mRNAs and lincRNAs in HER2 tumors versus matched normal tissue. The TCGA represents one of the most comprehensive studies of RNA expression in thousands of cancers including breast cancer [34]. We mined the TCGA RNA-seq data for HER2? breast tumors with matched normal tissue. In total, we identified 12 patients (24 RNA-seq samples: from 12 tumors and 12 matched normal tissue) (Supporting file 9). Within this TCGA cohort, we identified 2,521 mRNAs and 283 lincRNAs to be differentially expressed between these tumors and their matched normal pairs (Fig. 5a, b, and Supporting files 10 and 11). We performed pathway analysis of differentially expressed mRNAs in the TCGA dataset, and as expected, cancer-related pathways were highly enriched (Fig. 5c). To determine the expression patterns of the 44 mRNAs and 3 lincRNAs that are affected in response to both HER2 knockdown (in BT474 cells) and HER2 inhibition (in tumors in vivo), we examined their expression in the TCGA RNA-seq data. Of the 44 mRNAs identified, 35 mRNAs are also dysregulated in the TCGA cohort. We graphed the expression values of each of these 35 mRNAs in all three datasets: HER2 inhibition in vivo (Post-trastuzumab/Pretrastuzumab), HER2 knockdown in BT474 cells (siHER2/ siControl), and TCGA cohort (Tumor/Normal). Strikingly, each mRNA shows a similar directionality in HER2 knockdown in BT474 cells and in HER2 inhibition in tumors in vivo, and as expected, negatively correlated in the TCGA cohort (tumor/normal) (Fig. 6). Gene names are shown in Supporting File 12. We also examined the expression of the three intersected lincRNAs in our analysis of TCGA RNA-seq data. In this analysis, it was expected that a lincRNA that is up-regulated in response to HER2 depletion and inhibition to be downregulated in tumor samples in comparison to matched normal tissue, and vice versa. Linc-STARD6-2, which we found to be down-regulated when HER2 was knocked-down or inhibited, is up-regulated in 8/12 tumors in comparison to their
Breast Cancer Res Treat (2015) 150:321–334
327
A
B
C
Fig. 4 Identification of commonly affected mRNAs and lincRNAs post-HER2 inhibition in tumors in vivo and post-HER2 knockdown in cell culture by siRNAs. To identify mRNAs and lincRNAs that are likely critical targets of HER2, we intersected differentially expressed mRNAs and lincRNAs that we have identified in pre- vs post-HER2 inhibition by trastuzumab in tumors in vivo and differentially expressed mRNAs and lincRNAs in BT474 cells treated with either control or HER2 siRNAs. We identified 44 mRNAs (a) and 3
lincRNAs (b) to overlap between the two datasets.c Of the 44 common mRNAs identified in our aforementioned analysis, several genes are known components of the E2F transcription factor network and PLK1 signaling pathway. These two pathways were also identified in our analysis of gene expression of HER2? tumors and matched normal breast tissues from TCGA, suggesting an important role for these two pathways in HER2? breast cancer
matched normal tissue (Fig. 7a). Linc-GJA1-2 and lincSLC39A10-10, which become up-regulated in response to HER2 knockdown or inhibition, are down-regulated in 10/11 and 9/12 tumors in comparison to their matched normal tissue, respectively (Fig. 7b, c).
possible to examine all differentially expressed genes (coding and non-coding) in greater depth and accuracy [35]. In addition to mRNAs, we can now analyze the entire transcriptome including novel classes of non-coding RNAs such as lincRNAs. Many of these non-coding RNAs are regulatory in nature, and can greatly influence mRNA gene expression patterns suggesting that their expression may also affect human health and disease [8, 25, 36, 37]. Thus, in our current study, we have utilized several RNA-seq datasets that take advantage of both cell culture models of HER2? breast cancer and clinically relevant tumor data in vivo to identify both lincRNAs and mRNAs that are affected in HER2? breast cancer. These RNA-seq datasets were generated from RNA isolated from tumors pre- and post-inhibition of HER2 by trastuzumab in vivo, RNA isolated from BT474 cells treated with either negative control or HER2 siRNAs and RNA isolated from HER2? breast tumors and their matched normal controls. By interrogating these datasets and integrating the results, we are able to identify a small set of mRNAs and lincRNAs that are affected by HER2 amplification during tumorigenesis. We believe that these mRNA and lincRNA candidates are the most promising for follow-up functional studies.
Discussion The amplification of the HER2 gene in 20–30 % of earlystage breast cancer patients demonstrates the important role of this genetic event in breast cancer tumorigenesis. To determine the effects of HER2 amplification on gene expression, a number of previous studies have examined changes in mRNA expression in HER2? breast cancer cell lines in response to HER2 inhibition using mRNA microarrays. However, these microarrays are limited in both the dynamic range as well as in the number of gene targets that can be interrogated. Also, these microarrays examined a very small number of long noncoding RNAs. Furthermore, since these studies were all performed using cell lines, it was not possible to determine their relevance to what takes place in tumors in vivo. The recent advances in high-throughput sequencing technologies made it
123
328
Breast Cancer Res Treat (2015) 150:321–334
A
B
C
Fig. 5 Dysregulation of mRNAs and lincRNAs in HER2? tumors. We utilized public RNA-seq of 12 HER2? tumors and 12 matched normal breast tissues to identify differentially expressed mRNAs and lincRNAs. Differentially expressed mRNAs (a) and lincRNAs (b) in HER2? breast cancer patient tumors versus adjacent matched normal tissue obtained from The Cancer Genome Atlas (TCGA) are represented by heatmaps. In total, 2,521 mRNAs and 283 lincRNAs were identified as differentially expressed (see ‘‘Materials and
methods’’ section). c Cancer-associated pathways are affected when HER2 is amplified. Over enrichment analysis of differentially expressed mRNAs between tumors and matched control samples from 12 HER2? breast cancer patients (TCGA cohort) reveals several biological pathways that become altered due to increased HER2 expression in these tumors. A number of genes and p values associated with each pathway are also shown
In our current studies, we have identified a small number of mRNAs that are clearly dysregulated in HER2? tumors and become significantly affected by HER2 depletion or inhibition. Importantly, several of these genes are known to affect cell growth and proliferation and could potentially impact tumorigenesis as a result of HER2 amplification. For example, we identified the RRM2 gene, which has
been previously shown to be associated with decreased survival in breast cancer [38], to be highly up-regulated in HER2? tumors and is significantly down-regulated by both HER2 inhibition and depletion. Another key gene that we identified is TOP2A, which is required for gene transcription and has been a target of several anti-cancer drugs, is also affected by HER2 inhibition/depletion. A final
123
Breast Cancer Res Treat (2015) 150:321–334
329
Fig. 6 Identification of 35 mRNAs that are affected in all three datasets with the expected directionality of expression. Of the 44 mRNAs that are affected by both HER2 inhibition in tumors by trastuzumab and HER2 knock down in BT474 cells, 35 of these mRNAs were also deregulated in HER2? tumors vs normal tissues (TCGA). These 35 mRNA transcripts are found to be dysregulated in all three datasets: TCGA HER2? (Tumor/Normal), Trastuzumab
clinical trial samples (Post-/Pre-trastuzumab treatment), and in BT474 cells (HER2 siRNA/control siRNA). We graphed the expression values of these 35 mRNAs in all three datasets, and strikingly each mRNA shows similar directionality in HER2 knockdown in BT474 cells and in HER2 inhibition in vivo, and as expected, negatively correlated in the TCGA cohort (tumor/normal)
example is PRC1, which is highly up-regulated during S and G2/M phases of mitosis is also down-regulated by both HER2 knock down and trastuzumab-mediated inhibition (See supporting file #12 for complete list of genes). These few examples provide insights into how HER2 amplification leads to increased cell proliferation during tumorigenesis, and potentially opens the door for designing combined therapies targeting both HER2 by trastuzumab and other genes by either currently approved or new drugs. In future studies, we will carry out high-throughput drug screens to identify new drugs that target specific proteins identified in our current study. Although the overexpression of HER2 in breast cancer has been known and studied for over 20 years, most studies have mainly focused on identifying proteins as therapeutic targets. However, recent advances in developing epigenetic-based therapies suggest that it is also critical to
explore non-coding RNAs as drug targets. Since many lincRNAs are also known to interact with epigenetic complexes and in some cases guide them to their genomic targets [5, 39], this opens up the opportunity for RNAbased therapeutic approaches. There are a few advantages to targeting lincRNAs instead of epigenetic protein complexes directly. First, because many lincRNAs have tissuespecific expression, therapeutically targeting lincRNAs in various cancers may reduce negative side effects in comparison to current treatments. Additionally, since some lincRNAs direct chromatin-modifying complexes to specific genomic loci [23, 39], this opens up the possibility of lincRNA-mediated epigenetic reprograming of genomic regions responsible for dysregulation in cancer [40]. While the developments of RNA-based therapies are still in their infancy, some successes have been achieved. For example, a recent study has targeted the extra chromosome 21 in
123
330 Fig. 7 Validation of lincRNAs expression in 12 tumors vs matched normal tissue (TCGA cohort). We have identified three lincRNAs that are affected by both HER2 inhibition and HER2 knockdown in tumors and BT474 cells, respectively. To determine the potential role of these lincRNAs in HER2? cancer, we examined their expression in HER2? tumors and matched normal tissues from TCGA RNA-seq datasets. The fold change of a lincSTARD6-2, b linc-GJA1-2, and c linc-SLC39A10-10 was graphed for each TCGA HER2? patient sample (Tumor/ Normal). Linc-STARD6-2 shows up-regulation in 8/12 TCGA Tumor/Normal samples. linc-GJA1-2 and LincSLC39A10-10 show downregulation in 10/11 and 9/12 tumor samples compared to matched normal controls, respectively. One patient (A1LB) did not show expression of linc-GJA1-2 in either normal or tumor tissue. Three samples, A0DZ, A1EN, and A18U, have no linc-GJA1-2 expression in the tumor sample preventing calculation of the exact fold change, and thus, bars were extended to the edge of the graph for these samples to indicate that there is an observable difference between tumor and normal samples
Breast Cancer Res Treat (2015) 150:321–334
A
B
C
Down Syndrome for heterochromatization and silencing by the lincRNA Xist using genome editing technologies [41]. Therapeutic approaches that target lincRNAs could include small molecules or antagonistic oligonucleotides that would prevent lincRNA–protein interactions, since most lincRNAs function within ribonucleoprotein (RNP)
123
complexes [42]. A direct inhibition or down-regulation of lincRNA expression could also be utilized; however, a major challenge in targeting lincRNAs is their low expression levels in human cells in comparison to mRNAs [8, 28, 43–46]. This is due to the fact that many lincRNAs are regulatory molecules, and in many cases a few copies per
Breast Cancer Res Treat (2015) 150:321–334
331
cell are sufficient for these lincRNAs to exert their effects [21]. Thus, in our analysis of RNA-seq datasets, we have set a 0.25 fpkm or higher cutoff for a lincRNA and 1 fpkm or higher for an mRNA to be called expressed and to be able to capture lowly expressed lincRNAs that can still exert biological activity in cells. In summary, our current study led to the identification of both lincRNAs and mRNAs that appear to be involved in HER2 regulatory network. Our future work will focus on elucidating the functional significance of these transcripts in breast cancer pathology.
deemed differentially expressed if the fold change was C2.0 or B0.5. Heatmaps were generated in R using heatmap.2, and Z scores were scaled by row using standard Z score calculation of log fold change.
Materials and methods
RNA isolation from cell lines
Next generation RNA-sequencing files
RNA from BT474 cells was isolated using RNeasyÒ Mini Kit (Qiagen) according to the manufacturer’s protocol. An added DNase (Qiagen) treatment step was included after the first wash to remove DNA contamination of RNA preps.
All raw files have been deposited in GEO: GSE60182. Clinical trial information Institute: Brown University Oncology Group (Brown University, Yale University, Cedar-Sinae Center), PI: William Sikov MD; Correlative Science PI: Lyndsay Harris MD; BrUOG Study ID: BR-211B; Clinical Trials.gov ID NCT00617942. RNA isolation and next generation RNA sequencing (RNA-seq) from tumor samples (clinical trial) Frozen biopsied cores were processed for RNA isolation using AllPrep (Qiagen), and the Ovation RNA-seq System (NuGen) was used for RNA amplification. Library preparation was performed using TruSeq v3 (Illumina) and then sequenced on an Illumina HiSeq 2500. The sequenced reads were aligned to the human genome version hg19 using GSNAP (Wu and Nacu, 2010). The uniquely aligned reads were further analyzed using Cufflinks V2.0.2 and aligned to human mRNAs and lincRNA databases (Cabili et al., 2011). Expression values were calculated as FPKM (fragment per kilobase of exon per million of mapped fragments) and were used to determine expression of lincRNAs and mRNAs in pre- and post-treatment samples. Transcripts were called expressed if FPKM values across either all pre-treatment samples or all post-treatment samples were C1.0 or C0.25, for mRNAs and lincRNAs, respectively. The mean expression level was calculated, and differences in expression between pre- and post-treatment samples were assessed for the 11 patients that achieved a pathological complete response (pCR) using the non-parametric Wilcoxon test for paired samples. Additionally, fold change (post/pre) was calculated to identify differentially expressed transcripts. Transcripts were
Cell culture of breast cancer cell lines Human breast cancer cell lines were grown in Hybri-Care Medium (ATCCÒ 46-XTM) supplemented with 10 % fetal bovine serum (FBS)(Bioexpress) and 100 units/ml of Penicillin and 100 lg/ml Streptomycin (Life Technologies) at 37 °C with 5 % CO2.
siRNA transfections The knockdown of HER2 was achieved through the transfection of HER2-specific siRNAs (Life Technologies, Catalogue #4390824, s611 and s613) at a final concentration of 15 nM with Lipofectamine RNAiMax (Life Technologies) at 7.5 ll/well of 6-well plate. Negative control siRNA #1 and #2 (Ambion Cat. #M4611 and AM4613) were used as negative controls. Western blot analyses Protein lysates were prepared with Laemmli Sample Buffer (BioRad) and separated on a gradient 4–20 % SDS-PAGE Mini-ProteanÒ TGXTM Gels (BioRad). The gel was transferred to a nitrocellulose membrane (Thermo Scientific) and probed with primary antibodies overnight. Antib-actin (Ambion, AM4302, 3.1 mg/ml) was used as a loading control at a dilution of 1:1000, and anti-HER2/ Erb2 (Cell Signaling, 2242S) was used to detect HER2 at a dilution of 1:500. Anti-mouse HRP (Thermo Scientific, 32230) and anti-rabbit HRP (Abcam, ab6721) were used as secondary antibodies. HRP was activated using SuperSignalÒ West Pico Chemiluminescent Substrate (Thermo Scientific) for autoradiography. Next generation RNA sequencing of BT474 cells After RNA was isolated from BT474 samples, we assessed the quality of RNA using BioRad Experion. RNA samples with RNA integrity number (RIN) larger than 8 (max is 10) were considered high quality and suitable for RNA-seq.
123
332
Library preparation was performed using ScriptseqTM Complete Gold (Human/Mouse/Rat) (Illumina) and sequenced on Illumina HiSeq 2500. Six samples were run on a single flow cell. We generated 100 bp paired-end strandspecific sequences, which were mapped to human genome release hg19 using TopHat with 2 mismatches allowed for full-length reads. The raw reads were mapped to human genes annotated in RefSeq database and lincRNAs annotated in Cabili et al. [4] using Cufflinks V2.0.2, and subsequently used for differential gene expression analysis after normalizing the values to the total mapped reads in each sample, see supporting file 14. Expression values were calculated as FPKM (fragment per kilobase of exon per million of mapped fragments) and were used to determine expression of both lincRNAs and mRNAs in HER2 knock down (KD) and negative control (NC) siRNA samples. Transcripts were considered expressed if FPKM values across either all BT474 HER2 KD samples or all NC samples were C1.0 for mRNAs and C0.25 for lincRNAs. The mean expression level was calculated, and from this statistically significant differences in expression between HER2 KD samples and NC samples was determined using a paired t test in R. Additionally, fold changes (HER2 KD/ NC) were calculated to identify differentially expressed transcripts. Transcripts were deemed differentially expressed if the fold change was C2.0 or B0.5. Heatmaps were generated in R using heatmap.2, and Z scores were scaled by row using standard Z score calculation of log fold change. Real-time quantitative PCR (RT-qPCR) Real-time quantitative PCR (RT-qPCR): RNA was converted to cDNA using RNA to cDNA EcoDryTM Premix Random Hexamers (Clontech). Primers pairs were designed using primer3 software, and most primers used were designed to span exon–exon boundaries. A complete list of all primers is included (Supporting file 13). Maxima SyBr Green/ROX qPCR Master Mix (Thermo Scientific) was used for qRT-PCR. A comparative CT quantitation was performed with a hold stage of 50 °C for 2 min and 95 °C for 10 min followed by 409 cycle of 95 °C for 15 s and 60 °C for 1 min and finally melt curve at 95 °C for 15 s, 60 °C for 1 min, and a ramp to 95 °C at 0.3 °C increments. Analysis was done using the 2 DDCT method with GAPDH as the reference gene [47]. Analysis of next generation RNA sequencing from The Cancer Genome Atlas (TCGA) RNA sequencing (RNA-seq) fastq files from 12 tumor and 12 adjacent matched normal breast cancer pairs were
123
Breast Cancer Res Treat (2015) 150:321–334
obtained through The Cancer Genome Atlas (TCGA) consortium via the Cancer Genomics Hub (https://cghub. ucsc.edu) [2] (Supporting file 9). Sequences were aligned to UCSC hg19 with default parameters in TopHat v2.0.11 specifying an unstranded library strategy [48]. Aligned sequences were then assembled into both mRNA and lincRNA transcripts using default parameters in Cufflinks v2.1.1 [48]. Relative transcript abundance was reported from Cufflinks as Fragments Per Kilobase of exon per Million fragments mapped (FPKM). Transcripts were deemed as expressed if FPKM values across either all tumor samples or all matched normal samples were C1.0 for mRNAs and C0.25 for lincRNAs. The non-parametric Wilcoxon test for paired samples was used to test for statistical significance (p \ 0.05) using the statistical software R [R Development Core Team (2011) R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria]. Fold changes were calculated as tumor/normal, and transcripts with greater than a twofold change were selected as differentially expressed. Over enrichment analysis of differentially expressed mRNAs using NCBI cancer pathways was also performed in R using fisher’s exact test. P values were corrected for multiple testing using the Bonferroni correction. Heatmaps were generated in R using heatmap.2, and Z-scores were scaled by row using standard Z-score calculation of log fold change. Acknowledgments We would like to thank Maria Sandoval, Jennifer Yori, and Ruth Keri at CWRU Pharmacology for BT474 cells and other reagents; Vinay Varadan at Case Comprehensive Cancer Center and Sheldon Bai at the RNA Center for discussion of bioinformatic analyses; Megan E Forrest, Jessica Sabers, Thomas LaFramboise, and Anthony Wynshaw-Boris for discussion of results. Conflict of interest Callie R. Merry, Sarah McMahon, Cheryl L. Thompson, Kristy L.S. Miskimen, Lyndsay Harris, and Ahmad M Khalil declare that there is no conflict of interest to be disclosed.
References 1. Siegel R, Naishadham D, Jemal A (2012) Cancer statistics, 2012. CA Cancer J Clin 62(1):10–29 2. Cancer Genome Atlas Network (2012) Comprehensive molecular portraits of human breast tumours. Nature 490(7418):61–70 3. Mitri Z, Constantine T, O’Regan R (2012) The HER2 receptor in breast cancer: pathophysiology, clinical use, and new advances in therapy. Chemother Res Pract 2012:743193 4. Cabili MN, Trapnell C, Goff L, Koziol M, Tazon-Vega B, Regev A, Rinn JL (2011) Integrative annotation of human large intergenic noncoding RNAs reveals global properties and specific subclasses. Genes Dev 25(18):1915–1927 5. Khalil AM, Guttman M, Huarte M, Garber M, Raj A, Rivea Morales D, Thomas K, Presser A, Bernstein BE, van Oudenaarden A et al (2009) Many human large intergenic noncoding
Breast Cancer Res Treat (2015) 150:321–334
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16. 17. 18.
19.
20.
21.
22.
23.
RNAs associate with chromatin-modifying complexes and affect gene expression. Proc Natl Acad Sci USA 106(28):11667–11672 Guttman M, Amit I, Garber M, French C, Lin MF, Feldser D, Huarte M, Zuk O, Carey BW, Cassady JP et al (2009) Chromatin signature reveals over a thousand highly conserved large noncoding RNAs in mammals. Nature 458(7235):223–227 Sauvageau M, Goff LA, Lodato S, Bonev B, Groff AF, Gerhardinger C, Sanchez-Gomez DB, Hacisuleyman E, Li E, Spence M et al (2013) Multiple knockout mouse models reveal lincRNAs are required for life and brain development. eLife 2(0):e01749 Fatica A, Bozzoni I (2014) Long non-coding RNAs: new players in cell differentiation and development. Nat Rev Genet 15(1):7–21 Grote P, Wittler L, Hendrix D, Koch F, Wahrisch S, Beisaw A, Macura K, Blass G, Kellis M, Werber M et al (2013) The tissuespecific lncRNA Fendrr is an essential regulator of heart and body wall development in the mouse. Dev Cell 24(2):206–214 Kretz M, Webster DE, Flockhart RJ, Lee CS, Zehnder A, LopezPajares V, Qu K, Zheng GX, Chow J, Kim GE et al (2012) Suppression of progenitor differentiation requires the long noncoding RNA ANCR. Genes Dev 26:338–343 Guttman M, Donaghey J, Carey BW, Garber M, Grenier JK, Munson G, Young G, Lucas AB, Ach R, Bruhn L et al (2011) lincRNAs act in the circuitry controlling pluripotency and differentiation. Nature 477(7364):295–300 Loewer S, Cabili MN, Guttman M, Loh YH, Thomas K, Park IH, Garber M, Curran M, Onder T, Agarwal S et al (2010) Large intergenic non-coding RNA-RoR modulates reprogramming of human induced pluripotent stem cells. Nat Genet 42(12):1113–1117 Huarte M, Guttman M, Feldser D, Garber M, Koziol MJ, Kenzelmann-Broz D, Khalil AM, Zuk O, Amit I, Rabani M et al (2010) A large intergenic noncoding RNA induced by p53 mediates global gene repression in the p53 response. Cell 142(3):409–419 Gupta RA, Shah N, Wang KC, Kim J, Horlings HM, Wong DJ, Tsai MC, Hung T, Argani P, Rinn JL et al (2010) Long noncoding RNA HOTAIR reprograms chromatin state to promote cancer metastasis. Nature 464(7291):1071–1076 Moran VA, Perera RJ, Khalil AM (2012) Emerging functional and mechanistic paradigms of mammalian long non-coding RNAs. Nucleic Acids Res 40:6391–6400 Wang KC, Chang HY (2011) Molecular mechanisms of long noncoding RNAs. Mol Cell 43(6):904–914 Clark MB, Mattick JS (2011) Long noncoding RNAs in cell biology. Semin Cell Dev Biol 22(4):366–376 Qureshi IA, Mattick JS, Mehler MF (2010) Long non-coding RNAs in nervous system function and disease. Brain Res 1338:20–35 Carrieri C, Cimatti L, Biagioli M, Beugnet A, Zucchelli S, Fedele S, Pesce E, Ferrer I, Collavin L, Santoro C et al (2012) Long noncoding antisense RNA controls Uchl1 translation through an embedded SINEB2 repeat. Nature 491(7424):454–457 Yildirim E, Kirby JE, Brown DE, Mercier FE, Sadreyev RI, Scadden DT, Lee JT (2013) Xist RNA is a potent suppressor of hematologic cancer in mice. Cell 152(4):727–742 Wang KC, Yang YW, Liu B, Sanyal A, Corces-Zimmerman R, Chen Y, Lajoie BR, Protacio A, Flynn RA, Gupta RA et al (2011) A long noncoding RNA maintains active chromatin to coordinate homeotic gene expression. Nature 472(7341):120–124 Ulitsky I, Shkumatava A, Jan CH, Sive H, Bartel DP (2011) Conserved function of lincRNAs in vertebrate embryonic development despite rapid sequence evolution. Cell 147(7):1537–1550 Tsai MC, Manor O, Wan Y, Mosammaparast N, Wang JK, Lan F, Shi Y, Segal E, Chang HY (2010) Long noncoding RNA as
333
24.
25.
26. 27. 28.
29. 30.
31.
32.
33.
34.
35. 36. 37.
38.
39.
40.
41.
modular scaffold of histone modification complexes. Science 329(5992):689–693 Cesana M, Cacchiarelli D, Legnini I, Santini T, Sthandier O, Chinappi M, Tramontano A, Bozzoni I (2011) A long noncoding RNA controls muscle differentiation by functioning as a competing endogenous RNA. Cell 147(2):358–369 Niland CN, Merry CR, Khalil AM (2012) Emerging roles for long non-coding RNAs in cancer and neurological disorders. Front Genet 3:25 Taft RJ, Pang KC, Mercer TR, Dinger M, Mattick JS (2010) Noncoding RNAs: regulators of disease. J Pathol 220(2):126–139 Mattick JS (2009) The genetic signatures of noncoding RNAs. PLoS Genet 5(4):e1000459 Iyer MK, Niknafs YS, Malik R, Singhal U, Sahu A, Hosono Y, Barrette TR, Prensner JR, Evans JR, Zhao S et al (2015) The landscape of long noncoding RNAs in the human transcriptome. Nat Genet Morris KV, Mattick JS (2014) The rise of regulatory RNA. Nat Rev Genet 15(6):423–437 Yang Z, Zhou L, Wu LM, Lai MC, Xie HY, Zhang F, Zheng SS (2011) Overexpression of long non-coding RNA HOTAIR predicts tumor recurrence in hepatocellular carcinoma patients following liver transplantation. Ann Surg Oncol 18(5):1243–1250 Kogo R, Shimamura T, Mimori K, Kawahara K, Imoto S, Sudo T, Tanaka F, Shibata K, Suzuki A, Komune S et al (2011) Long noncoding RNA HOTAIR regulates polycomb-dependent chromatin modification and is associated with poor prognosis in colorectal cancers. Cancer Res 71(20):6320–6326 Sendur MA, Aksoy S, Ozdemir NY, Zengin N, Altundag K (2012) What is the mechanism of progression with trastuzumab treatment—escape or resistance? Asian Pacific J Cancer Prevent APJCP 13(11):5929–5930 Le XF, Lammayot A, Gold D, Lu Y, Mao W, Chang T, Patel A, Mills GB, Bast RC Jr (2005) Genes affecting the cell cycle, growth, maintenance, and drug sensitivity are preferentially regulated by anti-HER2 antibody through phosphatidylinositol 3-kinase-AKT signaling. J Biol Chem 280(3):2092–2104 Cancer Genome Atlas Research N, Weinstein JN, Collisson EA, Mills GB, Shaw KR, Ozenberger BA, Ellrott K, Shmulevich I, Sander C, Stuart JM (2013) The cancer genome atlas pan-cancer analysis project. Nat Genet 45(10):1113–1120 Wang Z, Gerstein M, Snyder M (2009) RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet 10(1):57–63 Ulitsky I, Bartel DP (2013) lincRNAs: genomics, evolution, and mechanisms. Cell 154(1):26–46 Moran VA, Perera RJ, Khalil AM (2012) Emerging functional and mechanistic paradigms of mammalian long non-coding RNAs. Nucleic Acids Res 40(14):6391–6400 Putluri N, Maity S, Kommagani R, Creighton CJ, Putluri V, Chen F, Nanda S, Bhowmik SK, Terunuma A, Dorsey T et al (2014) Pathway-centric integrative analysis identifies RRM2 as a prognostic marker in breast cancer associated with poor survival and tamoxifen resistance. Neoplasia 16(5):390–402 Kaneko S, Bonasio R, Saldana-Meyer R, Yoshida T, Son J, Nishino K, Umezawa A, Reinberg D (2014) Interactions between JARID2 and noncoding RNAs regulate PRC2 recruitment to chromatin. Mol Cell 53(2):290–300 Moskalev EA, Schubert M, Hoheisel JD (2012) RNA-directed epigenomic reprogramming: an emerging principle of a more targeted cancer therapy? Genes Chromosomes Cancer 51(2):105–110 Jiang J, Jing Y, Cost GJ, Chiang JC, Kolpa HJ, Cotton AM, Carone DM, Carone BR, Shivak DA, Guschin DY et al (2013) Translating dosage compensation to trisomy 21. Nature 500(7462):296–300
123
334 42. Khalil AM, Rinn JL (2011) RNA-protein interactions in human health and disease. Semin Cell Dev Biol 22(4):359–365 43. Mustafi D, Kevany BM, Bai X, Maeda T, Sears JE, Khalil AM, Palczewski K (2013) Evolutionarily conserved long intergenic non-coding RNAs in the eye. Hum Mol Genet 22:2992–3002 44. Geisler S, Coller J (2013) RNA in unexpected places: long noncoding RNA functions in diverse cellular contexts. Nat Rev Mol Cell Biol 14(11):699–712 45. Wapinski O, Chang HY (2011) Long noncoding RNAs and human disease. Trends Cell Biol 21(6):354–361
123
Breast Cancer Res Treat (2015) 150:321–334 46. Dinger ME, Pang KC, Mercer TR, Crowe ML, Grimmond SM, Mattick JS (2009) NRED: a database of long noncoding RNA expression. Nucleic Acids Res 37(Database issue):D122–D126 47. Livak KJ, Schmittgen TD (2001) Analysis of relative gene expression data using real-time quantitative PCR and the 2(-delta delta C(T)) Method. Methods 25(4):402–408 48. Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L (2012) Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc 7(3):562–578
PRECLINICAL STUDY
Integrative transcriptome-wide analyses reveal critical HER2regulated mRNAs and lincRNAs in HER2+ breast cancer Callie R. Merry • Sarah McMahon • Cheryl L. Thompson • Kristy L. S. Miskimen Lyndsay N. Harris • Ahmad M. Khalil
•
Received: 7 October 2014 / Accepted: 28 February 2015 / Published online: 8 March 2015 Ó Springer Science+Business Media New York 2015
Abstract Breast cancer is a major health problem affecting millions of women worldwide. Over 200,000 new cases are diagnosed annually in the USA, with approximately 40,000 of these cases resulting in death. HER2-positive (HER2?) breast tumors, representing 20–30 % of early-stage breast cancer diagnoses, are characterized by the amplification of the HER2 gene. However, the critical genes and pathways that become affected by HER2 amplification in humans are yet to be specifically identified. Furthermore, it is yet to be determined if HER2 amplification also affects the expression of long intervening non-coding (linc)RNAs, which are involved in the epigenetic regulation of gene expression. We examined changes in gene expression by next generation RNA sequencing in human tumors pre- and post- HER2 inhibition by trastuzumab in vivo, and changes in gene expression in response to HER2 knock down in cell culture models. We integrated our results with gene expression analysis of
HER2? tumors vs matched normal tissue from The Cancer Genome Atlas. The integrative analyses of these datasets led to the identification of a small set of mRNAs, and the associated biological pathways that become deregulated by HER2 amplification. Furthermore, our analyses identified three lincRNAs that become deregulated in response to HER2 amplification both in vitro and in vivo. Our results should provide the foundation for functional studies of these candidate mRNAs and lincRNAs to further our understanding of how HER2 amplification results in tumorigenesis. Also, the identified lincRNAs could potentially open the door for future RNA-based biomarkers and therapeutics in HER2? breast cancer. Keywords Breast cancer HER2 Gene expression lincRNAs RNA-seq Genome-wide
Introduction Electronic supplementary material The online version of this article (doi:10.1007/s10549-015-3327-1) contains supplementary material, which is available to authorized users. C. R. Merry S. McMahon A. M. Khalil (&) Departments of Genetics and Genome Sciences, Case Western Reserve University, Case Western Reserve University, Cleveland, OH 44106, USA e-mail: [email protected] C. R. Merry A. M. Khalil Departments of Biochemistry, Case Western Reserve University, Case Western Reserve University, Cleveland, OH 44106, USA
Breast cancer is the leading cause of cancer among women, with over 200,000 new cases diagnosed in the United C. L. Thompson Departments of Family Medicine and Community Health, Case Western Reserve University, Case Western Reserve University, Cleveland, OH 44106, USA L. N. Harris Departments of Medicine, Case Western Reserve University, Case Western Reserve University, Cleveland, OH 44106, USA
C. L. Thompson K. L. S. Miskimen L. N. Harris A. M. Khalil Case Comprehensive Cancer Center, Case Western Reserve University, Case Western Reserve University, Cleveland, OH 44106, USA
123
322
States each year [1]. Breast cancer is composed of several subtypes [2], with 20–30 % of early-stage breast cancers characterized by the amplification of the human epidermal growth factor receptor 2 (HER2) gene, also known as HER2/neu and ERBB2 [3]. In normal cells, HER2 signaling is initiated through heterodimerization of HER2 with other ERBB family members in response to extracellular signals. Dimerization and subsequent autophosphorylation through the tyrosine kinase domain initiate a signaling cascade that ultimately leads to changes in gene expression patterns that regulate cell proliferation and growth. However, how HER2 gene amplification in breast cancer leads to uncontrolled cell proliferation is yet to be fully understood [3]. We postulated that identifying the key genes (both coding and non-coding) that become deregulated in response to HER2 amplification could provide important insights into how HER2 amplification affects cell proliferation. Furthermore, identifying novel noncoding genes, such as lincRNAs, that become deregulated when HER2 becomes amplified may potentially provide clues into global changes in gene expression patterns observed in HER2? breast cancer (see below). The human genome encodes, in addition to mRNAs and small non-coding RNAs, over 8,300 long intervening noncoding RNAs (lincRNAs) [4–6]. Recent studies have demonstrated important roles for lincRNAs in embryonic and organ development in vivo [7–9]. Additionally, gene expression studies revealed that each lincRNA can regulate numerous mRNA genes, mostly at genomic sites far away from a lincRNA site of transcription, and thus, lincRNAs are thought to exert their effects in trans [5, 10–14]. The impact of lincRNAs on gene expression patterns has implicated lincRNAs in a wide range of biological functions including dosage compensation, genomic imprinting, alternative splicing, nuclear organization, and regulation of mRNA translation [15–20]. Intriguingly, lincRNAs utilize various mechanisms to exert their effects in the cell. These mechanisms include acting as scaffolds and guides for chromatin-modifying complexes to specific regions of the genome, serving as decoys that regulate transcription factors binding to specific DNA sequences, and as microRNA ‘‘sponges’’ that regulate miRNA:mRNA interactions [5, 11–14, 21–24]. Recently, there has been growing interest in unraveling the roles of lincRNAs in human disease with the hope that these novel transcripts can be used as diagnostic biomarkers and/or therapeutic targets [18, 25–27]. In addition, there is a significant interest in studying the functional roles of lincRNAs in cancer since numerous lincRNAs have been shown to be dysregulated across multiple cancer types [8, 28, 29]. These studies have also shown that the expression of some lincRNAs correlates with clinical parameters, such as overall patient prognosis
123
Breast Cancer Res Treat (2015) 150:321–334
and metastasis [14, 25, 30, 31]. For example, the lincRNA HOTAIR is highly up-regulated in breast cancer, hepatocellular carcinoma, and colorectal cancers [14, 30, 31]. Additionally, in vitro and in vivo functional studies have implicated HOTAIR in promoting cancer metastasis [14]. The lincRNA linc-p21 is transcriptionally activated by the tumor suppressor p53 and subsequently represses a subset of p53 gene targets in trans through binding of hnRNP-K, thereby facilitating the p53 response in cells [13]. A final example focuses on the X inactive specific transcript (Xist), which has been known for decades to be the key transcript for X chromosome inactivation (Xi) in mammalian females. Recently, it has been shown that genetic deletion of Xist in hematopoietic stem cells results in a variety of hematological cancers and premature death in mice [20]. Taken together, these selected examples, as well as many others, demonstrate important roles of lincRNAs in tumorigenesis and metastasis, and their potential utilization as biomarkers and/or therapeutic targets. To further our understanding of the effects of HER2 signaling on gene expression of both mRNAs as well as novel non-coding genes such as lincRNAs, we utilized a combined in vitro and in vivo transcriptomic approach to pinpoint critical downstream genes. The expression of these identified mRNAs and lincRNAs was subsequently examined in RNA-seq datasets obtained from The Cancer Genome Atlas (TCGA) project to further refine our top candidates. Our results led to the identification of potentially key mRNAs and lincRNAs that may contribute to HER2? breast cancer.
Results Inhibition of HER2 by trastuzumab in vivo affects the expression of both mRNAs and lincRNAs Previous studies have shown that HER2 inhibition in cell culture models of HER2? breast cancer leads to dramatic changes in the expression of mRNAs. However, it is not currently known which of these mRNAs are affected as a result of HER2 inhibition in tumors in vivo, as well as the effects of HER2 inhibition on regulatory non-coding RNAs such as lincRNAs. To identify the mRNAs and lincRNAs that become affected by HER2 inhibition in vivo, we analyzed RNA-seq data from a clinical trial that we originally designed to predict benefit from the HER2 inhibitor, trastuzumab (HerceptinÒ) (Clinical Trials.gov ID NCT00617942, also see ‘‘Materials and methods’’ section). Trastuzumab is a monoclonal antibody drug designed to target the HER2 receptor by potentially inhibiting its heterodimerization [32]. The clinical trial accrued 80 patients of which 50 pairs of tumors could be biopsied pre-
Breast Cancer Res Treat (2015) 150:321–334
trastuzumab and post-trastuzumab (post-one dose, *10–14 days). All patients continued receiving a combination of chemotherapy and trastuzumab for 4 months, and a subset of thirteen tumor pairs representing the extremes of response to treatment was subjected to total RNA isolation and RNA-seq analysis (see ‘‘Materials and methods’’ section). Eleven out of these thirteen patients were identified as responders to trastuzumab as measured by pathological complete response (pCR) (Supporting file 1). We reasoned that since trastuzumab was an effective therapeutic drug in these eleven HER2? breast cancer patients, changes in gene expression pre- and post-one dose of trastuzumab would reflect mRNAs and lincRNAs regulated through HER2. Therefore, we identified differentially expressed mRNAs and lincRNAs pre- versus postone dose of trastuzumab treatment in our RNA-seq dataset of these eleven patients that were responsive to treatment as measured by pCR. We calculated FPKM (fragments per kilobase of exon per million fragments mapped) values of each known mRNA and lincRNA in the human genome in each RNAseq sample (22 total: 11 pre- vs. 11 post-trastuzumab). We calculated an average FPKM of each mRNA and lincRNA in pre- versus post-trastuzumab, and subsequently calculated a fold change of Post/Pre for each mRNA and lincRNA. We identified 228 mRNAs and 28 lincRNAs to be differentially expressed (Fig. 1a, b, and Supporting Files 2 and 3). We also performed pathway analysis on differentially expressed mRNAs and identified five cancer-related pathways to be significantly affected (p \ 0.05, p values ranged from 0.006 to 0.00004) (Fig. 1c). In summary, we identified a set of mRNAs and lincRNAs that become affected when HER2 signaling is inhibited in vivo. However, because of genetic and environmental differences of these 11 patients, we observed heterogeneity in the level of differential expression of both mRNAs and lincRNAs. Thus, to further refine our list of key mRNAs and lincRNAs that are critical components of the HER2 pathway, we decided to modulate HER2 in a cell culture model (see below). HER2 signaling affects the expression of mRNAs as well as lincRNAs in cell culture model of HER2? breast cancer We began our studies in BT474 cells, which are HER2? breast cancer cells that have been utilized extensively as a cell culture model to study HER2? breast cancer. To identify mRNAs and lincRNAs that are downstream of HER2 signaling in BT474 cells, we utilized a loss-offunction approach in which we depleted HER2 using validated siRNAs and examined changes in gene expression using RNA-seq. We transfected BT474 cells with
323
siRNAs targeting HER2, and simultaneously transfected the same number of cells with negative control siRNAs. We initially examined HER2 mRNA levels using RTqPCR at 48 and 72 h. At these time points, we found that the siRNAs were effective at knocking down HER2 by 80 and 87 %, respectively (Fig. 2a). Subsequently, we performed a new round of transfections at 24, 48, and 72 h, and examined HER2 protein levels by Western blot analysis. At 24 h, there was a modest reduction in HER2 protein levels; however, a significant reduction in HER2 protein levels in comparison to cells transfected with negative control siRNAs was achieved 48 h post-transfection (Fig. 2b). Based on these gene expression analyses of HER2 knock down with siRNAs, the 48 h post-transfection time point was chosen as optimal to sufficiently deplete HER2 at both the mRNA and protein levels, while still capturing some of the early changes in mRNAs and lincRNAs expression due to loss of HER2. To that end, equal numbers of BT474 cells were transfected with either HER2 siRNAs or negative control siRNAs in three biological replicates, and at 48 h post-transfection total RNA was isolated, quantified, and subjected to RNA-seq (see ‘‘Materials and methods’’ section). We found 1,015 mRNAs (303 up-regulated and 712 down-regulated) and 167 lincRNAs (139 up-regulated and 28 down-regulated) to be differentially expressed by Ctwofold in the HER2-depleted BT474 cells in comparison to cells transfected with negative control siRNAs (Fig. 2c, d, and Supporting Files 4 and 5). To place our findings into context, we identified a previous study that utilized mRNA microarrays to examine the effects of HER2 inhibition by trastuzumab on mRNAs in BT474 and SKBR3 cells [33]. This study reported sixteen mRNA genes that are significantly down-regulated in response to HER2 inhibition by trastuzumab. We found that eight of these sixteen mRNA genes were also significantly down-regulated in our current study. We designed primers for those mRNA genes, and confirmed in an independent set of knockdown experiments that their expression is altered in response to HER2 depletion in BT474 breast tumor cells by RT-qPCR (Fig. 2e). Furthermore, pathway analysis of differentially expressed mRNAs revealed numerous affected pathways post-HER2 depletion in BT474 cells (Supporting File 6). Since no previous studies of lincRNAs modulated in response to HER2 depletion or inhibition are available, we selected five lincRNAs from our RNA-seq data for validation by qPCR in an independent knockdown experiment in BT474 cells. We selected these lincRNAs based on fold changes and p values and we included both up- and down-regulated lincRNAs in response to HER2 depletion. Four out of the five lincRNAs showed statistically significant up- or down-regulation in response to HER2 depletion by RT-qPCR similar to what we observed
123
324
Breast Cancer Res Treat (2015) 150:321–334
A
B
C
Fig. 1 Inhibition of HER2 in vivo by trastuzumab results in significant changes in the expression of mRNAs and lincRNAs in HER2? breast cancer patients. Gene expression analysis by next generation RNA sequencing (RNA-seq) was performed on RNA isolated from HER2? tumors pre- and post-one dose of trastuzumab. Expression of genes was calculated as FPKM values and differentially expressed genes between pre- and post-treatment conditions were
identified. Differentially expressed mRNAs (a) and lincRNAs (b) in all 22 RNA-seq samples (11 pre- vs 11 post-trastuzumab) are represented by heatmaps; c Over enrichment analysis of differentially expressed mRNAs pre- versus post-trastuzumab reveals several key pathways affected in vivo when HER2 signaling is inhibited. A number of genes and p values are shown for each pathway
by RNA-seq analysis (Fig. 2f). To further confirm that these lincRNAs are downstream of the HER2 pathway, and that changes in their expression are not due to off-target effects of siRNAs, we utilized a second independent siRNA against HER2. First, we confirmed that this second siRNA is effective at knocking down HER2 protein levels by Western blot analysis (Fig. 3a). Next, we performed qPCR analysis on the same 5 lincRNAs (see Fig. 2f), and found that four out of the five lincRNAs show a similar response to knocking down HER2 with the second siRNA,
but only three lincRNAs pass a p value of \0.05 (Fig. 3b). Also, the expression of lincMCl1-1, which was not affected with first siRNA, was responsive to the second siRNA, similar to what we observed by RNA-seq. These experiments demonstrate that HER2 depletion in BT474 cells affects the expression of both mRNAs and lincRNAs; however, some variability is observed due to off-target effects of siRNA-mediated depletion. Thus, to overcome these limitations, we identified mRNAs and lincRNAs that are affected by both HER2 inhibition by trastuzumab
123
Breast Cancer Res Treat (2015) 150:321–334 Fig. 2 HER2 affects the expression of mRNAs and lincRNAs in BT474 cells. a Real-time qPCR analysis of HER2 mRNA levels in BT474 cells transfected with either control siRNAs or HER2 siRNAs (15 nM final concentration). We observed a significant knockdown of HER2 mRNA levels at 48 h postsiRNA transfections. b Western blot analysis of HER2 protein levels in BT474 cells transfected with control siRNAs versus cells transfected with HER2 siRNAs at 48 h posttransfections. A significant reduction in HER2 protein levels is observed. Actin was used as a loading control; c, d. Gene expression analysis by RNA-seq led to the identification of differentially expressed mRNAs and lincRNAs in BT474 cells treated with either control or HER2 siRNAs. Heatmaps of differentially expressed mRNAs and lincRNAs between BT474 cells transfected with control siRNAs vs HER2 siRNAs are shown. e, f. Validation of RNAseq data by RT-qPCR of eight mRNA genes and five lincRNAs in cells transfected with control siRNAs versus HER2 siRNAs. For both mRNAs and lincRNAs, we observed similar patterns in our original knockdown experiments analyzed by RNA-seq and subsequent knockdown experiments analyzed by RTqPCR
325
A
C
E
in vivo and HER2 depletion in BT474 cells by siRNAs (see below). To identify mRNAs and lincRNAs that are affected by both HER2 inhibition by trastuzumab in tumors and HER2 depletion in cell culture, we intersected differentially expressed mRNAs and lincRNAs that we identified postHER2 knockdown in BT474 cells with differentially
B
D
F
expressed mRNAs and lincRNAs that we identified in response to HER2 inhibition by trastuzumab in vivo. We found 44 mRNAs and 3 lincRNAs to be common between the two datasets (Fig. 4a, b, and Supporting file 7 and 8). Some of the 44 mRNA genes identified are key components of the PLK1 signaling pathway and E2F transcription factor network (Fig. 4c), which are key pathways affected
123
326
A
Breast Cancer Res Treat (2015) 150:321–334
role of these genes in HER2? breast cancer, we next examined their expression in a cohort of HER2? breast tumors and their matched normal control tissue in the Cancer Genome Atlas (TCGA) database. Thousands of mRNAs and hundreds of lincRNAs are differentially expressed in HER2? breast cancer tumors
B
Fig. 3 Validation of HER2-regulated lincRNAs in BT474 cells using a second independent siRNA. a We designed and confirmed the knockdown of HER2 protein using a second independent siRNA (distinct from first siRNA used for knock down experiments and subsequent RNA-seq analysis). Western blot analysis of HER2 protein levels in BT474 cells treated with siRNA#2 alone or in combination with first siRNA (siRNA pool) demonstrates the effectiveness of this second siRNA in knocking down HER2 protein in comparison to negative control siRNAs, b We examined the expression of the same 5 lincRNAs in Fig. 2f post-HER2 knock down in BT474 cells with a second independent siRNA by qPCR analysis. We found that HER2 knock down with a second independent siRNA also affects the expression of some of these lincRNAs. These results demonstrate that some of the lincRNAs identified in our cell culture studies are downstream of the HER2 signaling pathway and some are affected due to off-target effects of siRNAs
in HER2? tumors (see next ‘‘Result’’ section). In summary, we have identified a small set of mRNAs and lincRNAs that are affected in response to HER2 inhibition/ depletion in vivo and in cell culture suggesting an important role of not only mRNAs but also lincRNAs in HER2? breast cancer. To gain further insights into the potential
123
Using transcriptomic analyses pre- versus post-HER2 knockdown in cell culture and pre- versus post-HER2 inhibition in tumors, we are able to identify a small set of lincRNAs and mRNAs that are putative members of the HER2 regulatory network. To further validate our observations, we turned to RNA-seq data from The Cancer Genome Atlas (TCGA) project to determine the expression patterns of these mRNAs and lincRNAs in HER2 tumors versus matched normal tissue. The TCGA represents one of the most comprehensive studies of RNA expression in thousands of cancers including breast cancer [34]. We mined the TCGA RNA-seq data for HER2? breast tumors with matched normal tissue. In total, we identified 12 patients (24 RNA-seq samples: from 12 tumors and 12 matched normal tissue) (Supporting file 9). Within this TCGA cohort, we identified 2,521 mRNAs and 283 lincRNAs to be differentially expressed between these tumors and their matched normal pairs (Fig. 5a, b, and Supporting files 10 and 11). We performed pathway analysis of differentially expressed mRNAs in the TCGA dataset, and as expected, cancer-related pathways were highly enriched (Fig. 5c). To determine the expression patterns of the 44 mRNAs and 3 lincRNAs that are affected in response to both HER2 knockdown (in BT474 cells) and HER2 inhibition (in tumors in vivo), we examined their expression in the TCGA RNA-seq data. Of the 44 mRNAs identified, 35 mRNAs are also dysregulated in the TCGA cohort. We graphed the expression values of each of these 35 mRNAs in all three datasets: HER2 inhibition in vivo (Post-trastuzumab/Pretrastuzumab), HER2 knockdown in BT474 cells (siHER2/ siControl), and TCGA cohort (Tumor/Normal). Strikingly, each mRNA shows a similar directionality in HER2 knockdown in BT474 cells and in HER2 inhibition in tumors in vivo, and as expected, negatively correlated in the TCGA cohort (tumor/normal) (Fig. 6). Gene names are shown in Supporting File 12. We also examined the expression of the three intersected lincRNAs in our analysis of TCGA RNA-seq data. In this analysis, it was expected that a lincRNA that is up-regulated in response to HER2 depletion and inhibition to be downregulated in tumor samples in comparison to matched normal tissue, and vice versa. Linc-STARD6-2, which we found to be down-regulated when HER2 was knocked-down or inhibited, is up-regulated in 8/12 tumors in comparison to their
Breast Cancer Res Treat (2015) 150:321–334
327
A
B
C
Fig. 4 Identification of commonly affected mRNAs and lincRNAs post-HER2 inhibition in tumors in vivo and post-HER2 knockdown in cell culture by siRNAs. To identify mRNAs and lincRNAs that are likely critical targets of HER2, we intersected differentially expressed mRNAs and lincRNAs that we have identified in pre- vs post-HER2 inhibition by trastuzumab in tumors in vivo and differentially expressed mRNAs and lincRNAs in BT474 cells treated with either control or HER2 siRNAs. We identified 44 mRNAs (a) and 3
lincRNAs (b) to overlap between the two datasets.c Of the 44 common mRNAs identified in our aforementioned analysis, several genes are known components of the E2F transcription factor network and PLK1 signaling pathway. These two pathways were also identified in our analysis of gene expression of HER2? tumors and matched normal breast tissues from TCGA, suggesting an important role for these two pathways in HER2? breast cancer
matched normal tissue (Fig. 7a). Linc-GJA1-2 and lincSLC39A10-10, which become up-regulated in response to HER2 knockdown or inhibition, are down-regulated in 10/11 and 9/12 tumors in comparison to their matched normal tissue, respectively (Fig. 7b, c).
possible to examine all differentially expressed genes (coding and non-coding) in greater depth and accuracy [35]. In addition to mRNAs, we can now analyze the entire transcriptome including novel classes of non-coding RNAs such as lincRNAs. Many of these non-coding RNAs are regulatory in nature, and can greatly influence mRNA gene expression patterns suggesting that their expression may also affect human health and disease [8, 25, 36, 37]. Thus, in our current study, we have utilized several RNA-seq datasets that take advantage of both cell culture models of HER2? breast cancer and clinically relevant tumor data in vivo to identify both lincRNAs and mRNAs that are affected in HER2? breast cancer. These RNA-seq datasets were generated from RNA isolated from tumors pre- and post-inhibition of HER2 by trastuzumab in vivo, RNA isolated from BT474 cells treated with either negative control or HER2 siRNAs and RNA isolated from HER2? breast tumors and their matched normal controls. By interrogating these datasets and integrating the results, we are able to identify a small set of mRNAs and lincRNAs that are affected by HER2 amplification during tumorigenesis. We believe that these mRNA and lincRNA candidates are the most promising for follow-up functional studies.
Discussion The amplification of the HER2 gene in 20–30 % of earlystage breast cancer patients demonstrates the important role of this genetic event in breast cancer tumorigenesis. To determine the effects of HER2 amplification on gene expression, a number of previous studies have examined changes in mRNA expression in HER2? breast cancer cell lines in response to HER2 inhibition using mRNA microarrays. However, these microarrays are limited in both the dynamic range as well as in the number of gene targets that can be interrogated. Also, these microarrays examined a very small number of long noncoding RNAs. Furthermore, since these studies were all performed using cell lines, it was not possible to determine their relevance to what takes place in tumors in vivo. The recent advances in high-throughput sequencing technologies made it
123
328
Breast Cancer Res Treat (2015) 150:321–334
A
B
C
Fig. 5 Dysregulation of mRNAs and lincRNAs in HER2? tumors. We utilized public RNA-seq of 12 HER2? tumors and 12 matched normal breast tissues to identify differentially expressed mRNAs and lincRNAs. Differentially expressed mRNAs (a) and lincRNAs (b) in HER2? breast cancer patient tumors versus adjacent matched normal tissue obtained from The Cancer Genome Atlas (TCGA) are represented by heatmaps. In total, 2,521 mRNAs and 283 lincRNAs were identified as differentially expressed (see ‘‘Materials and
methods’’ section). c Cancer-associated pathways are affected when HER2 is amplified. Over enrichment analysis of differentially expressed mRNAs between tumors and matched control samples from 12 HER2? breast cancer patients (TCGA cohort) reveals several biological pathways that become altered due to increased HER2 expression in these tumors. A number of genes and p values associated with each pathway are also shown
In our current studies, we have identified a small number of mRNAs that are clearly dysregulated in HER2? tumors and become significantly affected by HER2 depletion or inhibition. Importantly, several of these genes are known to affect cell growth and proliferation and could potentially impact tumorigenesis as a result of HER2 amplification. For example, we identified the RRM2 gene, which has
been previously shown to be associated with decreased survival in breast cancer [38], to be highly up-regulated in HER2? tumors and is significantly down-regulated by both HER2 inhibition and depletion. Another key gene that we identified is TOP2A, which is required for gene transcription and has been a target of several anti-cancer drugs, is also affected by HER2 inhibition/depletion. A final
123
Breast Cancer Res Treat (2015) 150:321–334
329
Fig. 6 Identification of 35 mRNAs that are affected in all three datasets with the expected directionality of expression. Of the 44 mRNAs that are affected by both HER2 inhibition in tumors by trastuzumab and HER2 knock down in BT474 cells, 35 of these mRNAs were also deregulated in HER2? tumors vs normal tissues (TCGA). These 35 mRNA transcripts are found to be dysregulated in all three datasets: TCGA HER2? (Tumor/Normal), Trastuzumab
clinical trial samples (Post-/Pre-trastuzumab treatment), and in BT474 cells (HER2 siRNA/control siRNA). We graphed the expression values of these 35 mRNAs in all three datasets, and strikingly each mRNA shows similar directionality in HER2 knockdown in BT474 cells and in HER2 inhibition in vivo, and as expected, negatively correlated in the TCGA cohort (tumor/normal)
example is PRC1, which is highly up-regulated during S and G2/M phases of mitosis is also down-regulated by both HER2 knock down and trastuzumab-mediated inhibition (See supporting file #12 for complete list of genes). These few examples provide insights into how HER2 amplification leads to increased cell proliferation during tumorigenesis, and potentially opens the door for designing combined therapies targeting both HER2 by trastuzumab and other genes by either currently approved or new drugs. In future studies, we will carry out high-throughput drug screens to identify new drugs that target specific proteins identified in our current study. Although the overexpression of HER2 in breast cancer has been known and studied for over 20 years, most studies have mainly focused on identifying proteins as therapeutic targets. However, recent advances in developing epigenetic-based therapies suggest that it is also critical to
explore non-coding RNAs as drug targets. Since many lincRNAs are also known to interact with epigenetic complexes and in some cases guide them to their genomic targets [5, 39], this opens up the opportunity for RNAbased therapeutic approaches. There are a few advantages to targeting lincRNAs instead of epigenetic protein complexes directly. First, because many lincRNAs have tissuespecific expression, therapeutically targeting lincRNAs in various cancers may reduce negative side effects in comparison to current treatments. Additionally, since some lincRNAs direct chromatin-modifying complexes to specific genomic loci [23, 39], this opens up the possibility of lincRNA-mediated epigenetic reprograming of genomic regions responsible for dysregulation in cancer [40]. While the developments of RNA-based therapies are still in their infancy, some successes have been achieved. For example, a recent study has targeted the extra chromosome 21 in
123
330 Fig. 7 Validation of lincRNAs expression in 12 tumors vs matched normal tissue (TCGA cohort). We have identified three lincRNAs that are affected by both HER2 inhibition and HER2 knockdown in tumors and BT474 cells, respectively. To determine the potential role of these lincRNAs in HER2? cancer, we examined their expression in HER2? tumors and matched normal tissues from TCGA RNA-seq datasets. The fold change of a lincSTARD6-2, b linc-GJA1-2, and c linc-SLC39A10-10 was graphed for each TCGA HER2? patient sample (Tumor/ Normal). Linc-STARD6-2 shows up-regulation in 8/12 TCGA Tumor/Normal samples. linc-GJA1-2 and LincSLC39A10-10 show downregulation in 10/11 and 9/12 tumor samples compared to matched normal controls, respectively. One patient (A1LB) did not show expression of linc-GJA1-2 in either normal or tumor tissue. Three samples, A0DZ, A1EN, and A18U, have no linc-GJA1-2 expression in the tumor sample preventing calculation of the exact fold change, and thus, bars were extended to the edge of the graph for these samples to indicate that there is an observable difference between tumor and normal samples
Breast Cancer Res Treat (2015) 150:321–334
A
B
C
Down Syndrome for heterochromatization and silencing by the lincRNA Xist using genome editing technologies [41]. Therapeutic approaches that target lincRNAs could include small molecules or antagonistic oligonucleotides that would prevent lincRNA–protein interactions, since most lincRNAs function within ribonucleoprotein (RNP)
123
complexes [42]. A direct inhibition or down-regulation of lincRNA expression could also be utilized; however, a major challenge in targeting lincRNAs is their low expression levels in human cells in comparison to mRNAs [8, 28, 43–46]. This is due to the fact that many lincRNAs are regulatory molecules, and in many cases a few copies per
Breast Cancer Res Treat (2015) 150:321–334
331
cell are sufficient for these lincRNAs to exert their effects [21]. Thus, in our analysis of RNA-seq datasets, we have set a 0.25 fpkm or higher cutoff for a lincRNA and 1 fpkm or higher for an mRNA to be called expressed and to be able to capture lowly expressed lincRNAs that can still exert biological activity in cells. In summary, our current study led to the identification of both lincRNAs and mRNAs that appear to be involved in HER2 regulatory network. Our future work will focus on elucidating the functional significance of these transcripts in breast cancer pathology.
deemed differentially expressed if the fold change was C2.0 or B0.5. Heatmaps were generated in R using heatmap.2, and Z scores were scaled by row using standard Z score calculation of log fold change.
Materials and methods
RNA isolation from cell lines
Next generation RNA-sequencing files
RNA from BT474 cells was isolated using RNeasyÒ Mini Kit (Qiagen) according to the manufacturer’s protocol. An added DNase (Qiagen) treatment step was included after the first wash to remove DNA contamination of RNA preps.
All raw files have been deposited in GEO: GSE60182. Clinical trial information Institute: Brown University Oncology Group (Brown University, Yale University, Cedar-Sinae Center), PI: William Sikov MD; Correlative Science PI: Lyndsay Harris MD; BrUOG Study ID: BR-211B; Clinical Trials.gov ID NCT00617942. RNA isolation and next generation RNA sequencing (RNA-seq) from tumor samples (clinical trial) Frozen biopsied cores were processed for RNA isolation using AllPrep (Qiagen), and the Ovation RNA-seq System (NuGen) was used for RNA amplification. Library preparation was performed using TruSeq v3 (Illumina) and then sequenced on an Illumina HiSeq 2500. The sequenced reads were aligned to the human genome version hg19 using GSNAP (Wu and Nacu, 2010). The uniquely aligned reads were further analyzed using Cufflinks V2.0.2 and aligned to human mRNAs and lincRNA databases (Cabili et al., 2011). Expression values were calculated as FPKM (fragment per kilobase of exon per million of mapped fragments) and were used to determine expression of lincRNAs and mRNAs in pre- and post-treatment samples. Transcripts were called expressed if FPKM values across either all pre-treatment samples or all post-treatment samples were C1.0 or C0.25, for mRNAs and lincRNAs, respectively. The mean expression level was calculated, and differences in expression between pre- and post-treatment samples were assessed for the 11 patients that achieved a pathological complete response (pCR) using the non-parametric Wilcoxon test for paired samples. Additionally, fold change (post/pre) was calculated to identify differentially expressed transcripts. Transcripts were
Cell culture of breast cancer cell lines Human breast cancer cell lines were grown in Hybri-Care Medium (ATCCÒ 46-XTM) supplemented with 10 % fetal bovine serum (FBS)(Bioexpress) and 100 units/ml of Penicillin and 100 lg/ml Streptomycin (Life Technologies) at 37 °C with 5 % CO2.
siRNA transfections The knockdown of HER2 was achieved through the transfection of HER2-specific siRNAs (Life Technologies, Catalogue #4390824, s611 and s613) at a final concentration of 15 nM with Lipofectamine RNAiMax (Life Technologies) at 7.5 ll/well of 6-well plate. Negative control siRNA #1 and #2 (Ambion Cat. #M4611 and AM4613) were used as negative controls. Western blot analyses Protein lysates were prepared with Laemmli Sample Buffer (BioRad) and separated on a gradient 4–20 % SDS-PAGE Mini-ProteanÒ TGXTM Gels (BioRad). The gel was transferred to a nitrocellulose membrane (Thermo Scientific) and probed with primary antibodies overnight. Antib-actin (Ambion, AM4302, 3.1 mg/ml) was used as a loading control at a dilution of 1:1000, and anti-HER2/ Erb2 (Cell Signaling, 2242S) was used to detect HER2 at a dilution of 1:500. Anti-mouse HRP (Thermo Scientific, 32230) and anti-rabbit HRP (Abcam, ab6721) were used as secondary antibodies. HRP was activated using SuperSignalÒ West Pico Chemiluminescent Substrate (Thermo Scientific) for autoradiography. Next generation RNA sequencing of BT474 cells After RNA was isolated from BT474 samples, we assessed the quality of RNA using BioRad Experion. RNA samples with RNA integrity number (RIN) larger than 8 (max is 10) were considered high quality and suitable for RNA-seq.
123
332
Library preparation was performed using ScriptseqTM Complete Gold (Human/Mouse/Rat) (Illumina) and sequenced on Illumina HiSeq 2500. Six samples were run on a single flow cell. We generated 100 bp paired-end strandspecific sequences, which were mapped to human genome release hg19 using TopHat with 2 mismatches allowed for full-length reads. The raw reads were mapped to human genes annotated in RefSeq database and lincRNAs annotated in Cabili et al. [4] using Cufflinks V2.0.2, and subsequently used for differential gene expression analysis after normalizing the values to the total mapped reads in each sample, see supporting file 14. Expression values were calculated as FPKM (fragment per kilobase of exon per million of mapped fragments) and were used to determine expression of both lincRNAs and mRNAs in HER2 knock down (KD) and negative control (NC) siRNA samples. Transcripts were considered expressed if FPKM values across either all BT474 HER2 KD samples or all NC samples were C1.0 for mRNAs and C0.25 for lincRNAs. The mean expression level was calculated, and from this statistically significant differences in expression between HER2 KD samples and NC samples was determined using a paired t test in R. Additionally, fold changes (HER2 KD/ NC) were calculated to identify differentially expressed transcripts. Transcripts were deemed differentially expressed if the fold change was C2.0 or B0.5. Heatmaps were generated in R using heatmap.2, and Z scores were scaled by row using standard Z score calculation of log fold change. Real-time quantitative PCR (RT-qPCR) Real-time quantitative PCR (RT-qPCR): RNA was converted to cDNA using RNA to cDNA EcoDryTM Premix Random Hexamers (Clontech). Primers pairs were designed using primer3 software, and most primers used were designed to span exon–exon boundaries. A complete list of all primers is included (Supporting file 13). Maxima SyBr Green/ROX qPCR Master Mix (Thermo Scientific) was used for qRT-PCR. A comparative CT quantitation was performed with a hold stage of 50 °C for 2 min and 95 °C for 10 min followed by 409 cycle of 95 °C for 15 s and 60 °C for 1 min and finally melt curve at 95 °C for 15 s, 60 °C for 1 min, and a ramp to 95 °C at 0.3 °C increments. Analysis was done using the 2 DDCT method with GAPDH as the reference gene [47]. Analysis of next generation RNA sequencing from The Cancer Genome Atlas (TCGA) RNA sequencing (RNA-seq) fastq files from 12 tumor and 12 adjacent matched normal breast cancer pairs were
123
Breast Cancer Res Treat (2015) 150:321–334
obtained through The Cancer Genome Atlas (TCGA) consortium via the Cancer Genomics Hub (https://cghub. ucsc.edu) [2] (Supporting file 9). Sequences were aligned to UCSC hg19 with default parameters in TopHat v2.0.11 specifying an unstranded library strategy [48]. Aligned sequences were then assembled into both mRNA and lincRNA transcripts using default parameters in Cufflinks v2.1.1 [48]. Relative transcript abundance was reported from Cufflinks as Fragments Per Kilobase of exon per Million fragments mapped (FPKM). Transcripts were deemed as expressed if FPKM values across either all tumor samples or all matched normal samples were C1.0 for mRNAs and C0.25 for lincRNAs. The non-parametric Wilcoxon test for paired samples was used to test for statistical significance (p \ 0.05) using the statistical software R [R Development Core Team (2011) R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria]. Fold changes were calculated as tumor/normal, and transcripts with greater than a twofold change were selected as differentially expressed. Over enrichment analysis of differentially expressed mRNAs using NCBI cancer pathways was also performed in R using fisher’s exact test. P values were corrected for multiple testing using the Bonferroni correction. Heatmaps were generated in R using heatmap.2, and Z-scores were scaled by row using standard Z-score calculation of log fold change. Acknowledgments We would like to thank Maria Sandoval, Jennifer Yori, and Ruth Keri at CWRU Pharmacology for BT474 cells and other reagents; Vinay Varadan at Case Comprehensive Cancer Center and Sheldon Bai at the RNA Center for discussion of bioinformatic analyses; Megan E Forrest, Jessica Sabers, Thomas LaFramboise, and Anthony Wynshaw-Boris for discussion of results. Conflict of interest Callie R. Merry, Sarah McMahon, Cheryl L. Thompson, Kristy L.S. Miskimen, Lyndsay Harris, and Ahmad M Khalil declare that there is no conflict of interest to be disclosed.
References 1. Siegel R, Naishadham D, Jemal A (2012) Cancer statistics, 2012. CA Cancer J Clin 62(1):10–29 2. Cancer Genome Atlas Network (2012) Comprehensive molecular portraits of human breast tumours. Nature 490(7418):61–70 3. Mitri Z, Constantine T, O’Regan R (2012) The HER2 receptor in breast cancer: pathophysiology, clinical use, and new advances in therapy. Chemother Res Pract 2012:743193 4. Cabili MN, Trapnell C, Goff L, Koziol M, Tazon-Vega B, Regev A, Rinn JL (2011) Integrative annotation of human large intergenic noncoding RNAs reveals global properties and specific subclasses. Genes Dev 25(18):1915–1927 5. Khalil AM, Guttman M, Huarte M, Garber M, Raj A, Rivea Morales D, Thomas K, Presser A, Bernstein BE, van Oudenaarden A et al (2009) Many human large intergenic noncoding
Breast Cancer Res Treat (2015) 150:321–334
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16. 17. 18.
19.
20.
21.
22.
23.
RNAs associate with chromatin-modifying complexes and affect gene expression. Proc Natl Acad Sci USA 106(28):11667–11672 Guttman M, Amit I, Garber M, French C, Lin MF, Feldser D, Huarte M, Zuk O, Carey BW, Cassady JP et al (2009) Chromatin signature reveals over a thousand highly conserved large noncoding RNAs in mammals. Nature 458(7235):223–227 Sauvageau M, Goff LA, Lodato S, Bonev B, Groff AF, Gerhardinger C, Sanchez-Gomez DB, Hacisuleyman E, Li E, Spence M et al (2013) Multiple knockout mouse models reveal lincRNAs are required for life and brain development. eLife 2(0):e01749 Fatica A, Bozzoni I (2014) Long non-coding RNAs: new players in cell differentiation and development. Nat Rev Genet 15(1):7–21 Grote P, Wittler L, Hendrix D, Koch F, Wahrisch S, Beisaw A, Macura K, Blass G, Kellis M, Werber M et al (2013) The tissuespecific lncRNA Fendrr is an essential regulator of heart and body wall development in the mouse. Dev Cell 24(2):206–214 Kretz M, Webster DE, Flockhart RJ, Lee CS, Zehnder A, LopezPajares V, Qu K, Zheng GX, Chow J, Kim GE et al (2012) Suppression of progenitor differentiation requires the long noncoding RNA ANCR. Genes Dev 26:338–343 Guttman M, Donaghey J, Carey BW, Garber M, Grenier JK, Munson G, Young G, Lucas AB, Ach R, Bruhn L et al (2011) lincRNAs act in the circuitry controlling pluripotency and differentiation. Nature 477(7364):295–300 Loewer S, Cabili MN, Guttman M, Loh YH, Thomas K, Park IH, Garber M, Curran M, Onder T, Agarwal S et al (2010) Large intergenic non-coding RNA-RoR modulates reprogramming of human induced pluripotent stem cells. Nat Genet 42(12):1113–1117 Huarte M, Guttman M, Feldser D, Garber M, Koziol MJ, Kenzelmann-Broz D, Khalil AM, Zuk O, Amit I, Rabani M et al (2010) A large intergenic noncoding RNA induced by p53 mediates global gene repression in the p53 response. Cell 142(3):409–419 Gupta RA, Shah N, Wang KC, Kim J, Horlings HM, Wong DJ, Tsai MC, Hung T, Argani P, Rinn JL et al (2010) Long noncoding RNA HOTAIR reprograms chromatin state to promote cancer metastasis. Nature 464(7291):1071–1076 Moran VA, Perera RJ, Khalil AM (2012) Emerging functional and mechanistic paradigms of mammalian long non-coding RNAs. Nucleic Acids Res 40:6391–6400 Wang KC, Chang HY (2011) Molecular mechanisms of long noncoding RNAs. Mol Cell 43(6):904–914 Clark MB, Mattick JS (2011) Long noncoding RNAs in cell biology. Semin Cell Dev Biol 22(4):366–376 Qureshi IA, Mattick JS, Mehler MF (2010) Long non-coding RNAs in nervous system function and disease. Brain Res 1338:20–35 Carrieri C, Cimatti L, Biagioli M, Beugnet A, Zucchelli S, Fedele S, Pesce E, Ferrer I, Collavin L, Santoro C et al (2012) Long noncoding antisense RNA controls Uchl1 translation through an embedded SINEB2 repeat. Nature 491(7424):454–457 Yildirim E, Kirby JE, Brown DE, Mercier FE, Sadreyev RI, Scadden DT, Lee JT (2013) Xist RNA is a potent suppressor of hematologic cancer in mice. Cell 152(4):727–742 Wang KC, Yang YW, Liu B, Sanyal A, Corces-Zimmerman R, Chen Y, Lajoie BR, Protacio A, Flynn RA, Gupta RA et al (2011) A long noncoding RNA maintains active chromatin to coordinate homeotic gene expression. Nature 472(7341):120–124 Ulitsky I, Shkumatava A, Jan CH, Sive H, Bartel DP (2011) Conserved function of lincRNAs in vertebrate embryonic development despite rapid sequence evolution. Cell 147(7):1537–1550 Tsai MC, Manor O, Wan Y, Mosammaparast N, Wang JK, Lan F, Shi Y, Segal E, Chang HY (2010) Long noncoding RNA as
333
24.
25.
26. 27. 28.
29. 30.
31.
32.
33.
34.
35. 36. 37.
38.
39.
40.
41.
modular scaffold of histone modification complexes. Science 329(5992):689–693 Cesana M, Cacchiarelli D, Legnini I, Santini T, Sthandier O, Chinappi M, Tramontano A, Bozzoni I (2011) A long noncoding RNA controls muscle differentiation by functioning as a competing endogenous RNA. Cell 147(2):358–369 Niland CN, Merry CR, Khalil AM (2012) Emerging roles for long non-coding RNAs in cancer and neurological disorders. Front Genet 3:25 Taft RJ, Pang KC, Mercer TR, Dinger M, Mattick JS (2010) Noncoding RNAs: regulators of disease. J Pathol 220(2):126–139 Mattick JS (2009) The genetic signatures of noncoding RNAs. PLoS Genet 5(4):e1000459 Iyer MK, Niknafs YS, Malik R, Singhal U, Sahu A, Hosono Y, Barrette TR, Prensner JR, Evans JR, Zhao S et al (2015) The landscape of long noncoding RNAs in the human transcriptome. Nat Genet Morris KV, Mattick JS (2014) The rise of regulatory RNA. Nat Rev Genet 15(6):423–437 Yang Z, Zhou L, Wu LM, Lai MC, Xie HY, Zhang F, Zheng SS (2011) Overexpression of long non-coding RNA HOTAIR predicts tumor recurrence in hepatocellular carcinoma patients following liver transplantation. Ann Surg Oncol 18(5):1243–1250 Kogo R, Shimamura T, Mimori K, Kawahara K, Imoto S, Sudo T, Tanaka F, Shibata K, Suzuki A, Komune S et al (2011) Long noncoding RNA HOTAIR regulates polycomb-dependent chromatin modification and is associated with poor prognosis in colorectal cancers. Cancer Res 71(20):6320–6326 Sendur MA, Aksoy S, Ozdemir NY, Zengin N, Altundag K (2012) What is the mechanism of progression with trastuzumab treatment—escape or resistance? Asian Pacific J Cancer Prevent APJCP 13(11):5929–5930 Le XF, Lammayot A, Gold D, Lu Y, Mao W, Chang T, Patel A, Mills GB, Bast RC Jr (2005) Genes affecting the cell cycle, growth, maintenance, and drug sensitivity are preferentially regulated by anti-HER2 antibody through phosphatidylinositol 3-kinase-AKT signaling. J Biol Chem 280(3):2092–2104 Cancer Genome Atlas Research N, Weinstein JN, Collisson EA, Mills GB, Shaw KR, Ozenberger BA, Ellrott K, Shmulevich I, Sander C, Stuart JM (2013) The cancer genome atlas pan-cancer analysis project. Nat Genet 45(10):1113–1120 Wang Z, Gerstein M, Snyder M (2009) RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet 10(1):57–63 Ulitsky I, Bartel DP (2013) lincRNAs: genomics, evolution, and mechanisms. Cell 154(1):26–46 Moran VA, Perera RJ, Khalil AM (2012) Emerging functional and mechanistic paradigms of mammalian long non-coding RNAs. Nucleic Acids Res 40(14):6391–6400 Putluri N, Maity S, Kommagani R, Creighton CJ, Putluri V, Chen F, Nanda S, Bhowmik SK, Terunuma A, Dorsey T et al (2014) Pathway-centric integrative analysis identifies RRM2 as a prognostic marker in breast cancer associated with poor survival and tamoxifen resistance. Neoplasia 16(5):390–402 Kaneko S, Bonasio R, Saldana-Meyer R, Yoshida T, Son J, Nishino K, Umezawa A, Reinberg D (2014) Interactions between JARID2 and noncoding RNAs regulate PRC2 recruitment to chromatin. Mol Cell 53(2):290–300 Moskalev EA, Schubert M, Hoheisel JD (2012) RNA-directed epigenomic reprogramming: an emerging principle of a more targeted cancer therapy? Genes Chromosomes Cancer 51(2):105–110 Jiang J, Jing Y, Cost GJ, Chiang JC, Kolpa HJ, Cotton AM, Carone DM, Carone BR, Shivak DA, Guschin DY et al (2013) Translating dosage compensation to trisomy 21. Nature 500(7462):296–300
123
334 42. Khalil AM, Rinn JL (2011) RNA-protein interactions in human health and disease. Semin Cell Dev Biol 22(4):359–365 43. Mustafi D, Kevany BM, Bai X, Maeda T, Sears JE, Khalil AM, Palczewski K (2013) Evolutionarily conserved long intergenic non-coding RNAs in the eye. Hum Mol Genet 22:2992–3002 44. Geisler S, Coller J (2013) RNA in unexpected places: long noncoding RNA functions in diverse cellular contexts. Nat Rev Mol Cell Biol 14(11):699–712 45. Wapinski O, Chang HY (2011) Long noncoding RNAs and human disease. Trends Cell Biol 21(6):354–361
123
Breast Cancer Res Treat (2015) 150:321–334 46. Dinger ME, Pang KC, Mercer TR, Crowe ML, Grimmond SM, Mattick JS (2009) NRED: a database of long noncoding RNA expression. Nucleic Acids Res 37(Database issue):D122–D126 47. Livak KJ, Schmittgen TD (2001) Analysis of relative gene expression data using real-time quantitative PCR and the 2(-delta delta C(T)) Method. Methods 25(4):402–408 48. Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L (2012) Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc 7(3):562–578
