|Home | About | Journals | Submit | Contact Us | Français|
This study aimed to identify key genes contributing to pathological complete response (pCR) to chemotherapy by mRNA sequencing (RNA-seq).
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.
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.
We found key genes that may contribute to pCR to chemotherapy, such as PLCB4, ADCY6, and CNR1, as well as some transcription factors.
Breast cancer is a leading gynecologic cancer. Currently, neoadjuvant chemotherapy (NCT) is widely used for the treatment of early breast cancer . 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  and loss of bcl-2  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.
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.
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.
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).
Reads count data were first preprocessed using the edgeR package of R (version 3.1.0) , 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 .
DEGs were submitted to the online tool STRING  (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 .
Key proteins in the PPI network were identified by comparing the following parameters: degree of centrality , betweenness centrality , subgraph centrality , and closeness Centrality) . The parameters were calculated using a cytoscape plug-in CytoNCA .
Next, module analysis was performed to indent the subnetwork from the PPI network using ClusterONE . Pathway enrichment analysis of the genes in the module was also performed using DAVID; p<0.05 was used as a cutoff.
Transcriptional regulators are responsible for the transcriptional regulation of gene expression . Transcriptional factors that may regulate the DEGs in a module were searched by use of the anplug-in iRegulon  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.
Quality statistics of reads are shown in Table 1.
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).
The upregulated DEGs were significantly enriched in 12 KEGG pathways, and the downregulated ones were significantly enriched in 2 pathways.
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).
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.
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  and ovarian carcinoma by microarray analysis . 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 , 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 . 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  and has been suggested to have a critical role in breast cancer [25–27]. 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) , 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 . Remarkably, decrease in its expression has also been reported in chemosensitive serous epithelial ovarian cancer tumors , 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 . Sellar et al. reported that this gene can modulate cisplatin sensitivity in human epithelial ovarian cancer . 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 . Some of the HOX members (Hoxa9, Hoxa10, and Hoxa11) have been suggested to contribute to epithelial ovarian cancer .
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.
Disclosure of conflict of interest
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]