PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of medscimonISLHome PageSearch ArticleTable of ContentsSubmit ManuscriptSubscribe
 
Med Sci Monit. 2017; 23: 4321–4327.
Published online 2017 September 7. doi:  10.12659/MSM.903272
PMCID: PMC5600194

RNA Sequencing Uncovers Molecular Mechanisms Underlying Pathological Complete Response to Chemotherapy in Patients with Operable Breast Cancer

Abstract

Background

This study aimed to identify key genes contributing to pathological complete response (pCR) to chemotherapy by mRNA sequencing (RNA-seq).

Material/Methods

RNA was extracted from the frozen biopsy tissue of patients with pathological complete response and patients with non-pathological complete response. Sequencing was performed on the HiSeq2000 platform. Differentially expressed genes (DEGs) were identified between the pCR group and non-pCR (NpCR) group. Pathway enrichment analysis of DEGs was performed. A protein-protein interaction network was constructed, then module analysis was performed to identify a subnetwork. Finally, transcription factors were predicted.

Results

A total of 673 DEGs were identified, including 419 upregulated ones and 254 downregulated ones. The PPI network constructed consisted of 276 proteins forming 471 PPI pairs, and a subnetwork containing 18 protein nodes was obtained. Pathway enrichment analysis revealed that PLCB4 and ADCY6 were enriched in pathways renin secretion, gastric acid secretion, gap junction, inflammatory mediator regulation of TRP channels, retrograde endocannabinoid signaling, melanogenesis, cGMP-PKG signaling pathway, calcium signaling pathway, chemokine signaling pathway, cAMP signaling pathway, and rap1 signaling pathway. CNR1 was enriched in the neuroactive ligand-receptor interaction pathway, retrograde endocannabinoid signaling pathway, and rap1 signaling pathway. The transcription factor-gene network consists of 15 transcription factors and 16 targeted genes, of which 5 were downregulated and 10 were upregulated.

Conclusions

We found key genes that may contribute to pCR to chemotherapy, such as PLCB4, ADCY6, and CNR1, as well as some transcription factors.

MeSH Keywords: Achaete-Scute Complex Genome Region, Breast Neoplasms, Male, Clinical Conference

Background

Breast cancer is a leading gynecologic cancer. Currently, neoadjuvant chemotherapy (NCT) is widely used for the treatment of early breast cancer [1]. However, response to chemotherapy varies among patients, even when they have identical histological and molecular subtypes [2,3]. Pathological complete response (pCR) has been implicated as a surrogate indicator of long-term clinical benefit, such as disease-free survival and overall survival (OS), by many studies [4,5]. To date, only negative estrogen receptor (ER) status [6] and loss of bcl-2 [7] have been suggested as predictors of pCR in patients receiving neoadjuvant chemotherapy for operable and locally advanced breast cancer.

However, the mechanisms underlying pCR are largely unknown. Here, we attempted to use RNA sequencing (RNA-seq) technique to seek key genes and transcription factors that may contribute to pCR by comparing the transcriptome data between breast cancer patients with pCR to chemotherapy and those with non-pCR.

Material and Methods

Patient enrollment and sampling

Patients diagnosed with invasive breast cancer were included in this study. The inclusion criteria were: (1) patients with stage II or IIIA, operable breast cancer; (2) without distant metastasis; (3) without other serious disease, who can complete the chemotherapy as planned; (4) at least 10 ug total RNA can be extracted from the frozen biopsy tissue.

Finally, 6 patients, including 3 showing pCR (MXL-3-2, MXL-5-1 and MXL-P2) and 3 showing non-pCR (MXL-N1, MXL-N2 and MXL-1-3), were included. Written informed consent was provided by each patient and volunteer before sampling.

Total RNA extraction, library construction, and sequencing

A total of 10 g RNA was extracted from the frozen biopsy tissue of each patient, then poly(A)+ mRNA was isolated using Sera-mag Magnetic Oligo(dT) Beads (Thermo, part # 1004815). Afterwards, the purified RNA was broken into short segments in Fragment buffer. Next, reverse transcription was performed to construct a cDNA library. RNA concentration was measured with a Qubit® 2.0 Fluorometer and RNA integrity number (RIN) was measured by use of a Bioanalyzer 2100 (Agilent, CA, USA). Sequencing was performed on the HiSeq2000 platform (Illumina Inc., San Diego, CA) and single-end reads (50 bp in length) were generated.

Sequence quality control and alignment

First, reads were filtered by removing reads containing more than 40% bases with quality value <Q20 and reads containing more than 4% N bases. Sequence quality control was done with the fastx_toolkit. Next, clean reads were mapped to the human reference genome hg19 using TopHat 2.1.1 software (download site: http://ccb.jhu.edu/software/tophat/index.shtml).

Identification of differentially expressed genes correlated with mild and severe pneumonia

Reads count data were first preprocessed using the edgeR package of R (version 3.1.0) [8], which were then processed into log2-counts per million (logCPM) matrix by voom method using the R/limma package. Then, differential expression analysis between the pCR and NpCR groups was performed using the R/limma package, and | log2FC | >1 and p<0.05 were used as cutoffs for a differentially expressed gene (DEG). A heat map was also used to display the gene expression profile of differentially expressed genes, based on the hierarchical clustering results using R package gplots [9].

Construction of a protein-protein interaction network

DEGs were submitted to the online tool STRING [10] (version: 10.0; http://string-db.org/) to identify protein-protein interaction (PPI) pairs. Only pairs with PPI score of 0.7 or higher were retained to construct a PPI network using Cytoscape [11].

Key proteins in the PPI network were identified by comparing the following parameters: degree of centrality [12], betweenness centrality [13], subgraph centrality [14], and closeness Centrality) [15]. The parameters were calculated using a cytoscape plug-in CytoNCA [16].

Next, module analysis was performed to indent the subnetwork from the PPI network using ClusterONE [17]. Pathway enrichment analysis of the genes in the module was also performed using DAVID; p<0.05 was used as a cutoff.

Prediction of transcription factors and construction of a transcription factor-gene network

Transcriptional regulators are responsible for the transcriptional regulation of gene expression [18]. Transcriptional factors that may regulate the DEGs in a module were searched by use of the anplug-in iRegulon [19] of Cytoscape, which integrates transcription factor information from the databases Transfac, Jaspar, Encode, Swissregulon, and Homer. The parameters were set as follows: Minimum identity between orthologous genes: 0.05, and Maximum false discovery rate on motif similarity: 0.001. A gene-transcription factor pair with normalized enrichment score (NES) >3 was retained.

Result

Quality control of reads and statistics

Quality statistics of reads are shown in Table 1.

Table 1
Quality statistics of reads.

Identification of DEGs

A total of 673 DEGs were identified in patients with pCR to chemotherapy, including 419 upregulated ones and 254 downregulated ones. Overall, the identified DEGs can distinguish the pneumonia samples from the control sample (Figure 1).

Figure 1
Heat map showing the gene expression profile of differentially expressed genes based on their hierarchical clustering. The upregulated genes and downregulated genes are indicated by red and green, respectively.

The upregulated DEGs were significantly enriched in 12 KEGG pathways, and the downregulated ones were significantly enriched in 2 pathways.

PPI network construction and module analysis

The constructed PPI network consisted of 276 proteins forming 471 PPI pairs (Figure 2). The top 10 proteins for each centrality parameter are listed in Table 2.

Figure 2
The protein-to-protein interaction (PPI) network of differentially expressed genes. The upregulated proteins and downregulated proteins are indicated by red and green, respectively.
Table 2
The top 10 proteins for each centrality parameter in the protein-protein interaction network.

One subnetwork was obtained, which included 18 protein nodes forming 79 protein-protein interaction pairs (Figure 3). More importantly, the heatmap revealed that these 18 genes can accurately distinguish the pCR samples and NpCR samples (Figure 4). Pathway enrichment analysis revealed that PLCB4 and ADCY6 were enriched in the renin secretion pathway, gastric acid secretion pathway, gap junction pathway, inflammatory mediator regulation of TRP channels pathway, retrograde endocannabinoid signaling pathway, melanogenesis pathway, cGMP-PKG signaling pathway, calcium signaling pathway, chemokine signaling pathway, cAMP signaling pathway, and rap1 signaling pathway. CNR1 was enriched in the neuroactive ligand-receptor interaction pathway, the retrograde endocannabinoid signaling pathway, and the rap1 signaling pathway (Table 3).

Figure 3
Subnetwork of protein-to-protein interaction (PPI). The upregulated proteins and downregulated proteins are indicated by red and green, respectively.
Figure 4
Heatmap based on the 18 protein nodes. The upregulated proteins and downregulated proteins are indicated by red and green, respectively.
Table 3
KEGG pathway enrichemnt of genes in the protein-protein interaction subnetwork.

Construction of transcription factor-gene network

The transcription factor-gene network consisted of 15 transcription factors and 16 targeted genes (Figure 5). These transcription factors were also differentially expressed: 5 downregulated and 10 upregulated.

Figure 5
Transcription factor-gene network for subnetwork module. The upregulated proteins and downregulated proteins are indicated by red and green, respectively. Triangles represent the transcription factor and ellipses represents the target protein.

Discussions

Because pCR to chemotherapy is recognized as an independent prognostic and predictive factor for survival of BC patients, we attempted to learn more about what contributes to the difference between patients with pCR and those with NpCR. Based on the transcriptome data from RNA-seq, we first identified DEGs that were associated pCR, and then tried to find key proteins (such as PLCB4, ADCY6, CNR1, and MAPK14) to from a protein-protein interaction perspective.

Among the upregulated genes, PLCB4 and ADCY6 were speculated to contribute to pCR to chemotherapy via multiple pathways, such as via renin secretion, gastric acid secretion, gap junction, inflammatory mediator regulation of TRP channels, retrograde endocannabinoid signaling, melanogenesis, cGMP-PKG signaling pathway, calcium signaling pathway, chemokine signaling pathway, cAMP signaling pathway, and rap1 signaling pathway. CNR1 was a downregulated gene, which was enriched in pathways of neuroactive ligand-receptor interaction, retrograde endocannabinoid signaling, and rap1 signaling pathway.

ADCY6 encodes adenylyl cyclase 6, which is a membrane-associated enzyme catalyzing the formation of the secondary messenger cyclic adenosine monophosphate (cAMP). The differential expression of this gene has also been reported in breast cancer [20] and ovarian carcinoma by microarray analysis [21]. By microarray analysis, ADCY6 expression was found to be decreased by 2-fold more in the 81-fold gemcitabine-resistant variant MiaPaCa2-RG than in the gemcitabine-sensitive cell line MiaPaCa2 [22], suggesting that high ADCY6 expression is associated with gemcitabine sensitivity. This seems to agree with its upregulation in patients with pCR to chemotherapy.

PLCB4 encodes an enzyme 1-phosphatidylinositol-4,5-bisphosphate phosphodiesterase beta-4, which catalyzes the formation of inositol 1,4,5-trisphosphate (IP3) and diacylglycerol from phosphatidylinositol 4,5-bisphosphate (PIP2). Woodfield et al. identified the primary gene targets of TFAP2C in hormone responsive breast carcinoma cells and found that PLCB4 was induced by TFAP2C [23]. TFAP2C is a member of the retinoic acid-inducible developmentally regulated family of AP-2 factors, which regulates its target genes via a helix-loop-helix motif [24] and has been suggested to have a critical role in breast cancer [2527]. Thus, we supposed that PLCB4 may also have a role in the pathogenesis of breast cancer. Zheng et al. reported that PLCB4 was among the 12 genes showing significant upregulated expression (ratio >5) in acute myeloid leukemia cells (HL-60/MDR) with resistance to 6 ordinary antitumor drugs (Vincristine, Adriamycin, Etoposide, Daunorubicin, Mitomycin C, and Cyclophosphamide) [28], suggesting that PLCB4 contributes to resistance to antitumor drugs. However, this seems contradictory to our finding here that PLCB4 was upregulated in patients with complete response to chemotherapy. Thus, whether this helps to induce pCR to chemotherapy needs to be further studied. Further pathway enrichment analysis suggests that this gene may contribute to pCR to chemotherapy.

CNR1 was downregulated and encodes a G protein-coupled cannabinoid receptor (also known as cannabinoid CB1), located primarily in the central nervous system, and its activation is necessary for the antiemetic action of cannabinoids used for preventing emesis in cancer patients receiving chemotherapy [29]. Remarkably, decrease in its expression has also been reported in chemosensitive serous epithelial ovarian cancer tumors [30], suggesting that this gene may also have a role in breast cancer.

Interestingly, these 3 genes were predicted to be regulated by transcription factors such as BARX2, CRX, TOX, POU2F1, VAX2, HOXA6, HOXB7, MYBL1, and NFATC2. Some of these transcription factors were upregulated and some were downregulated, indicating that their regulation upon the expression of a gene is complicated. Among them, BARX2 was found to be downregulated in patients showing pCR to chemotherapy. BARX2 is located within a minimal region at 11q25, which has been reported to promote breast cancer progression [31]. Sellar et al. reported that this gene can modulate cisplatin sensitivity in human epithelial ovarian cancer [32]. In contrast, HOXA6 and HOXB7 were upregulated in patients showing pCR to chemotherapy. These 2 transcription factors are encoded by the homeobox (HOX) gene family, and were found to have an important role in embryogenesis and differentiation of adult cells [33]. Some of the HOX members (Hoxa9, Hoxa10, and Hoxa11) have been suggested to contribute to epithelial ovarian cancer [34].

Conclusions

By RNA-seq, we found some key genes that may contribute to pCR to chemotherapy, such as PLCB4, ADCY6, CNR1, as well as some transcription factors. Their roles need to be further validated by experimental evidence.

Footnotes

Disclosure of conflict of interest

None.

Source of support: This work was supported by the Natural Science Foundation of Zhejiang Province [No. LY13H160029]; the Key Science and Technology Project of Zhejiang Province [No. 2014C03004]; and the Natural Science Foundation of Zhejiang Province [No. LQ17H160013]

References

1. Kaufmann M, Von Minckwitz G, Bear H, et al. Recommendations from an international expert panel on the use of neoadjuvant (primary) systemic treatment of operable breast cancer: New perspectives 2006. Ann Oncol. 2007;18(12):1927–34. [PubMed]
2. Gajdos C, Tartter PI, Estabrook A, et al. Relationship of clinical and pathologic response to neoadjuvant chemotherapy and outcome of locally advanced breast cancer. Eur J Surg Oncol. 2002;80(1):4–11. [PubMed]
3. Rouzier R, Perou CM, Symmans WF, et al. Breast cancer molecular subtypes respond differently to preoperative chemotherapy. Clin Cancer Res. 2005;11(16):5678–85. [PubMed]
4. von Minckwitz G, Untch M, Blohmer J-U, et al. Definition and impact of pathologic complete response on prognosis after neoadjuvant chemotherapy in various intrinsic breast cancer subtypes. J Clin Oncol. 2012;30(15):1796–804. [PubMed]
5. Maas M, Nelemans PJ, Valentini V, et al. Long-term outcome in patients with a pathological complete response after chemoradiation for rectal cancer: A pooled analysis of individual patient data. Lancet Oncol. 2010;11(9):835–44. [PubMed]
6. Buzdar A, Valero V, Theriault R, et al. Pathological complete response to chemotherapy is related to hormone receptor status. Breast Cancer Res Treat. 2003;82:S69.
7. von Minckwitz G, Sinn H-P, Raab G, Loibl S, et al. Clinical response after two cycles compared to HER2, Ki-67, p53, and bcl-2 in independently predicting a pathological complete response after preoperative chemotherapy in patients with operable carcinoma of the breast. Breast Cancer Res. 2008;10(2):R30. [PMC free article] [PubMed]
8. Robinson MD, Mccarthy DJ, Smyth GK. edgeR: A Bioconductor package for differential expression analysis of digital gene expression dataSMotn. Bioinformatics. 2009;26(1):139–40. [PMC free article] [PubMed]
9. Warnes GR, Bolker B, Bonebakker L, et al. gplots: Various R programming tools for plotting data. 2009. https://rdrr.io/cran/gplots/
10. Szklarczyk D, Franceschini A, Wyder S, et al. STRING v10: Protein–protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 2015;43(Database issue):D447–52. [PMC free article] [PubMed]
11. Shannon P, Markiel A, Ozier O, et al. Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504. [PubMed]
12. Rito T, Deane CM, Reinert G. The importance of age and high degree, in protein-protein interaction networks. J Comput Biol. 2012;19(6):785–95. [PubMed]
13. Wang H, Hernandez JM, Van Mieghem P. Betweenness centrality in a weighted network. Phys Rev E Stat Nonlin Soft Matter Phys. 2008;77(4 Pt 2):046105. [PubMed]
14. Estrada E, Rodriguez-Velazquez JA. Subgraph centrality in complex networks. Phys Rev E Stat Nonlin Soft Matter Phys. 2005;71(5 Pt 2):056103. [PubMed]
15. Du Y, Gao C, Chen X, et al. A new closeness centrality measure via effective distance in complex networks. Chaos. 2015;25(3):033112. [PubMed]
16. Tang Y, Li M, Wang J, et al. CytoNCA: A cytoscape plugin for centrality analysis and evaluation of protein interaction networks. Biosystems. 2015;127:67–72. [PubMed]
17. Nepusz T, Yu H, Paccanaro A. Detecting overlapping protein complexes in protein-protein interaction networks. Nat Met. 2012;9(5):471–72. [PMC free article] [PubMed]
18. Blancafort P, Segal DJ. Designing transcription factor architectures for drug discovery. Mol Pharmacol. 2004;66(6):1361–71. [PubMed]
19. Verfaillie A, Imrichová H, Van de Sande B, et al. iRegulon: From a gene list to a gene regulatory network using large motif and track collections. PLoS Comput Biol. 2014;10(7):e1003731. [PMC free article] [PubMed]
20. Kumar R, Selth LA, Schulz RB, et al. Genome-wide mapping of ZNF652 promoter binding sites in breast cancer cells. J Cell Biochem. 2011;112(10):2742–47. [PubMed]
21. Ko H-H, Tseng K-J, Wei L-M, Tsai M-H. Possible disease-link genetic pathways constructed by hierarchical clustering and conditional probabilities of ovarian carcinoma microarray data. IEEE; 2010. pp. 1–4.
22. Nakahira S, Nakamori S, Tsujie M, et al. Involvement of ribonucleotide reductase M1 subunit overexpression in gemcitabine resistance of human pancreatic cancer. Int J Cancer. 2007;120(6):1355–63. [PubMed]
23. Woodfield GW, Chen Y, Bair TB, et al. Identification of primary gene targets of TFAP2C in hormone responsive breast carcinoma cells. Genes Chromosomes Cancer. 2010;49(10):948–62. [PMC free article] [PubMed]
24. Feng W, Williams T. Cloning and characterization of the mouse AP-2epsilon gene: A novel family member expressed in the developing olfactory bulb. Mol Cell Neurosci. 2003;24(2):460–75. [PubMed]
25. Allouche A, Nolens G, Tancredi A, et al. The combined immunodetection of AP-2α and YY1 transcription factors is associated with ERBB2 gene overexpression in primary breast tumors. Breast Cancer Res. 2008;10(1):R9. [PMC free article] [PubMed]
26. Turner BC, Zhang J, Gumbs AA, et al. Expression of AP-2 transcription factors in human breast cancer correlates with the regulation of multiple growth factor signalling pathways. Cancer Res. 1998;58(23):5466–72. [PubMed]
27. Pellikainen J, Naukkarinen A, Ropponen K, et al. Expression of HER2 and its association with AP-2 in breast cancer. Eur J Cancer. 2004;40(10):1485–95. [PubMed]
28. Zheng G-H, Fu J-R, Xu Y-H, et al. Screening and cloning of multi-drug resistant genes in HL-60/MDR cells. Leuk Res. 2009;33(8):1120–23. [PubMed]
29. Darmani NA. ΔD9-Tetrahydrocannabinol and synthetic cannabinoids prevent emesis produced by the cannabinoid CB1 receptor antagonist/inverse agonist SR 141716A. Neuropsychopharmacology. 2001;24(2):198–203. [PubMed]
30. Bachvarov D, L’Esperance S, Popa I, et al. Gene expression patterns of chemoresistant and chemosensitive serous epithelial ovarian tumors with possible predictive value in response to initial chemotherapy. Int J Oncol. 2006;29(4):919–34. [PubMed]
31. Stevens T, Meech R. BARX2 and estrogen receptor-α (ESR1) coordinately regulate the production of alternatively spliced ESR1 isoforms and control breast cancer cell growth and invasion. Oncogene. 2006;25(39):5426–35. [PubMed]
32. Sellar GC, Watt KP, Li L, et al. The homeobox gene BARX2 can modulate cisplatin sensitivity in human epithelial ovarian cancer. Int J Oncol. 2002;21(5):929–34. [PubMed]
33. Garcia-Fernàndez J. The genesis and evolution of homeobox gene clusters. Nat Rev Genet. 2005;6(12):881–92. [PubMed]
34. Cheng W, Liu J, Yoshida H, et al. Lineage infidelity of epithelial ovarian cancers is controlled by HOX genes that specify regional identity in the reproductive tract. Nat Med. 2005;11(5):531–37. [PubMed]

Articles from Medical Science Monitor : International Medical Journal of Experimental and Clinical Research are provided here courtesy of International Scientific Literature, Inc.