RNA-Sequence Analysis Reveals Differentially Expressed Genes (DEGs) in Patients Exhibiting Different Risks of Tumor Metastasis.

BACKGROUND Breast cancer is one of the most common malignancies in women. In a previous study, we found that for two patients who had a high risk of lymphatic metastasis, lymphatic metastasis did not occur; whereas, for two patients who had a low risk of lymphatic metastasis, lymphatic metastasis did occur. MATERIAL AND METHODS We analyzed the differential gene expressions of these four patients by RNA-sequence. The data (HRNM_T versus HRNM_N, LRYM_T versus LRYM_N, and HRNM_T versus LRYM_T) was then processed using differentially expressed genes (DEGs) analysis, functional analysis for DEGs, and PPI network construct. RESULTS For HRNM_T versus HRNM_N, there were 224 DEGs. There were 504 DEGs for LRYM_T versus LRYM_N, and 88 DEGs for LRYM_T versus LRYM_N. For HRNM_T versus HRNM_N, DEGs were up-regulated mainly in the cell cycle, the IL-17 signaling pathway, and the progesterone-mediated oocyte maturation; DEGs were down-regulated mainly in the IL-17 signaling pathway. For LRYM_T versus LRYM_N, DEGs were up-regulated mainly in protein digestion and absorption, and cytokine-cytokine receptor interaction; DEGs were down-regulated mainly in ECM-receptor interaction. For HRNM_T versus LRYM_T, DEGs were up-regulated mainly in the PPAR signaling pathway; DEGs were downregulated mainly in the adipocytokine signaling pathway. The DEGs were screened to construct PPI networks. CONCLUSIONS The GO and KEGG functional enrichments of HRNM_T versus HRNM_N, and LRYM_T versus LRYM_N were consistent with earlier studies. For HRNM_T versus LRYM_T, DEGs were up-regulated mainly in PPAR signaling; DEGs were down-regulated mainly in the adipocytokine pathway.

Background

Breast cancer is one of the most common malignancies in women, accounting for a quarter of new cases and 15% of deaths caused by cancer in females [1]. Being the leading cause of cancer death, breast cancer threatens the life and health of females and causes enormous impact on the economy, society, families, and psychological health of women [2]. In 2016, 1,677,000 women were newly diagnosed with breast cancer, representing 25.2% of all newly diagnosed cancers according to Cancer Statistics 2017 [3]. Among these new cases, more than 187,000 were from China and this number increased at twice the speed of the global growth rate [4]. Although the incidence of breast cancer in China is lower than the global rate, recently it has been rising sharply. According to the latest China Cancer Registration Report, the breast cancer incidence rate among Chinese women was 42.55/100,000 in 2012 [4]. Although this is a low incidence rate compared with the whole world, the growth rate has been twice the velocity of the global growth rate since 1990s, especially in urban areas. Currently, breast cancer is the most frequent cancer diagnosed in Chinese women and the sixth leading cause of death from cancer in Chinese women. The development of cancer involves multiply levels and stages. Metastasis plays a critical role in its progression, and ultimately leads to high fatality rates. In light of the complexity of metastasis itself, we need to more deeply explore its molecular mechanism including preclinical and clinical aspects. We need to identify early alarming signals and make timely interventions to improve therapeutic effects and quality of life, as well as relieve social and personal burdens. Tumor metastasis is very complex and involves many metastasis-associated molecules; there appears to be little effect on metastasis through single molecule-targeting interventions.

In this study, tumor tissue and the adjusted normal tissue from four patients were used for RNA-sequence analysis. Two of the patients had a high risk for lymphatic metastasis, but lymphatic metastasis had not occurred at the time of the study. The other two patients were of low risk for lymphatic metastasis, but lymphatic metastasis had occurred. We then analyzed the differential gene expression of these four patients by RNA-sequence in order to investigate elements of breast carcinoma metastasis.

Material and Methods

Ethical statement

This study was reviewed and approved by Ethics Committee of Department of Breast, Women’s Hospital, School of Medicine, Zhejiang University. Each patient involved in this study provided written informed consent.

Patient samples and clinicopathological characteristics

Samples were acquired from breast cancer patients who underwent surgical resection between May 2012 and December 2013 at the Department of Breast, Women’s Hospital, School of Medicine, Zhejiang University. Clinicopathological characteristic of each patient were recorded. Adjacent normal tissues were resected at least 3 centimeters away from the tumor, while tumor samples were carefully extracted from the center of the breast cancer carcinoma. All of the samples were frozen in liquid nitrogen and stored at −80°C. Four matched pairs of normal tissue samples and breast cancer tissues were obtained. These samples were used for RNA sequencing analysis.

RNA preparation

Total RNA from breast cancer tissue and adjacent normal tissue was extracted by TRIzol (Invitrogen, Carlsbad, CA, USA), re-dissolved in DEPC-treated water and quantified using NanoVue Plus spectrophotometry (GE Healthcare, Fairfield, CT, USA). RNA integrity was evaluated using agarose gel electrophoresis, and DNA contamination was eliminated using gDNA Eraser (Takara, Tokyo, Japan) according to the manufacturer’s guidelines.

RNA sequencing (RNA-Seq)

Breast cancer and adjacent normal tissue samples were ground in TRIzol and total RNA was isolated using the RNeasy Kit (Qiagen) according to the manufacturer’s guidelines. DNaseI was added in order to ensure that the tissues were not mixed with genomic DNA. RNA purity was detected by Agilent BioAnalyzer 2100 (Agilent, Santa Clara, CA, USA). Transcriptome sequencing libraries were prepared using the TruSeq RNA Access Library Prep Kit® (Illumina, CA, USA). Paired end sequencing (100 bp) was performed at an Illumina HiSeq 2000 instrument, which is supported by Beijing Novogene Biological Information Technology Co., LTD.

Differentially expressed genes (DEGs) analysis

The LIMMA package [5] in R language, a liner regression model, was used to identify the DEGs in breast tissue and adjacent tissue. A p value <0.05 and |log2 Fold Chance (FC)| >1.0 were set as thresholds to screen out differentially expressed genes (DEGs).

Functional analysis for DEGs

Database for Annotation, Visualization and Integrated Discovery (DAVID) [6] is considered an integrated biological knowledgebase and effective analytic tool to systematically extract biological information and explore the biological meaning behind lists of genes or proteins. GO [7] and KEGG [8] were used to perform the functional analysis of DEGs. We chose a p value less than 0.05 as cut-off criterion.

PPI network constructing

The PPI network corresponding to the screened DEGs was retrieved from the Search Tool for the Retrieval of Interacting Genes (STRING) database [9]. STRING is used for identifying known and predicted protein interactions. Proteins in the interaction network are represented with nodes, while the interactions between two proteins are represented with edges. The selection of the hub protein is based on the score of the nodes, which is calculated by the count of edges launching from a protein in the PPI network. The PPI networks in the present study were constructed using Cytoscape [10] when the confidence scores were more than 0.4.

Results

DEGs in breast cancer tissues

Patient clinicopathological characteristics are shown in Table 1, and four patients (Chenfx, Yuay, Dongyh, and Yumh) were involved in this study. In Table 2, we can see that Chenfx and Yuay showed high risk of metastasis, but there was no tumor metastasis. Dongyh and Yumh exhibited low risk of metastasis, but tumor metastasis occurred. Therefore, we analyzed the DEGs in breast cancer tissues (HRNM_T versus HRNM_N, LRYM_T versus LRYM_N, and HRNM_T versus LRYM_T). Table 3 shows that 107 genes were up-expressed and 117 genes were down-expressed in the comparison of HRNM_T versus HRNM_N. There were 220 genes up-expressed and 284 genes down-expressed in the comparison of LRYM_T versus LRYM_N. There were 61 genes up-expressed and 27 genes down-expressed in the comparison of HRNM_T versus LRYM_T.

Table 1

Patient clinicopathological characteristics.

	Chenfx	Yuay	Dongyh	Yumh
Age	43	42	46	46
Tumor size (cm)	3	4	5	3
Multifocality	0	0	0	0
Cancer embolus in vessel	1	0	0	0
ER	1	0	1	1
PR	1	0	1	1
HER2	0	1	0	0
ki67	1	1	0	0
Positive number of LN	0	0	4	8

Table 2

Abbreviation of breast tumor and adjusted tissue.

	High/low risk	Metastasis	Abbreviation
Chenfx_N	High	No	HRNM_N
Chenfx_T	High	No	HRNM_T
Yuay_N	High	No	HRNM_N
Yuay_T	High	No	HRNM_T
Dongyh_N	Low	Yes	LRYM_N
Dongyh_T	Low	Yes	LRYM_T
Yumh_N	Low	Yes	LRYM_N
Yumh_T	Low	Yes	LRYM_T

Table 3

Number of differentially expressed genes.

	Up	Down	Total
HRNM_T vs. HRNM_N	107	117	224
LRYM_T vs. LRYM_N	220	284	504
HRNM_T vs. LRYM_T	61	27	88

Significant function and pathways of DEGs

GO functional and KEGG pathway enrichment analyses were performed to both up- and down-regulated genes, and the results are shown in Figure 1 and Tables 4–6. According to the result, the DEGs in HRNM_T versus HRNM_N were enriched in multiple GO categories, mainly the cell cycle, the IL-17 signaling pathway, and the progesterone-mediated oocyte maturation. The DEGs in LRYM_T versus LRYM_N were enriched mainly in protein digestion and absorption, cytokine-cytokine receptor interaction and ECM-receptor interaction. The DEGs in HRNM_T versus LRYM_T were enriched mainly in the PPAR signaling pathway and the adipocytokine signaling pathway.

Figure 1

GO analysis.

Table 4

KEGG enrichment analyses for both up-regulated and down-regulated genes in HRNM_T vs. HRNM_N.

Category		Genes	p-Value
Up-regulated genes
KEGG_PATHWAY	Cell cycle	E2f1,2,3, CycA, Cdc6, Cdc25A, CDK1, CycB, Bub1, Mad2, BubR1, Mps1	4.17E-06
KEGG_PATHWAY	IL-17 signaling pathway	Eotxin, CXCL10, S100A7, S100A8, S100A9, MMP1, MMP13	7.49E-03
KEGG_PATHWAY	Progesterone-mediated oocyte maturation	PDE3, CycA, Cdc25, Bub1, Mad1/2/3, Cdc2, CycB	2.36E-04
Down-regulated genes
KEGG_PATHWAY	IL-17 signaling pathway	IL-17B	7.49E-03

KEGG – Kyoto encyclopedia of genes and genomes.

Table 5

KEGG enrichment analyses for both up-regulated and down-regulated genes in LRYM_T vs. LRYM_N.

Category		Genes	p-Value
Up-regulated genes
KEGG_PATHWAY	Protein digestion and absorption	FATCD36, FABP, LPL, PGAR, Penlipin, aP2, PEPCK, AQP7	6.11E-08
KEGG_PATHWAY	Cytokine-cytokine receptor interaction	CXCL1, CXCL3, CXCL5, CCL28, LEP, CNTFR, PDGFA, VEGFOP, KIT, VSF11A	5.80E-06
Down-regulated genes
KEGG_PATHWAY	ECM-receptor interaction	Collagen, THBS, Fibionectin	1.38 E-02

KEGG – Kyoto encyclopedia of genes and genomes.

Table 6

KEGG enrichment analyses for both up-regulated and down-regulated genes in HRNM_T vs. LRYM_T.

Category		Genes	p-Value
Up-regulated genes
KEGG_PATHWAY	PPAR signaling pathway	Penlipin	5.82E-05
KEGG_PATHWAY	Adipocytokine signaling pathway	ADIPO	1.13E-03

KEGG – Kyoto encyclopedia of genes and genomes.

Analysis of PPI network and module

The PPI network for the screened DEGs was constructed (Figures 2–4), and contained 129 nodes and 294 edges in HRNM_T versus HRNM_N (Figure 2), 321 nodes and 904 edges in LRYM_T versus LRYM_N (Figure 3), and 16 nodes and 13 edges in HRNM_T versus LRYM_T (Figure 4). Module analysis of PPI networks was performed by MCODE and the most significant module was selected, and was named subnetwork1.

Figure 2

PPI networks for DEGs in HRNM_T versus HRNM_N. The nodes represent proteins and edges represent pair-wise interactions. The red nodes represent the proteins encoded by up-regulated genes and the green represents the proteins encoded by down-regulated genes.

Figure 3

PPI networks for DEGs in LRYM_T versus LRYM_N. The nodes represent proteins and edges represent pair-wise interactions. The red nodes represent the proteins encoded by up-regulated genes and the green represent the proteins encoded by down-regulated genes.

Figure 4

PPI networks for DEGs in HRNM_T versus LRYM_T. The nodes represent proteins and edges represent pair-wise interactions. The red nodes represent the proteins encoded by up-regulated genes and the green represent the proteins encoded by down-regulated genes

We also used DAVID to identify significant enrichment of these genes in multiple GO categories and KEGG pathways of subnetwork1. DEGs in subnetwork1 were mainly enriched in GO terms related to the cell cycle in HRNM_T versus HRNM_N, protein digestion and absorption in LRYM_T versus LRYM_N and the PPAR signaling pathway in LRYM_T versus LRYM_N.

Discussion

Currently, breast cancer is the most frequent cancer in Chinese women and the sixth leading cause of death from cancer in Chinese women. The development of cancer involves multiply levels and stages. Previously, we analyzed clinicopathological characteristics of the study patients. We found Chenfx and Yuay showed high risk for metastasis, but no tumor metastasis occurred. We found that Dongyh and Yumh exhibited low risk of metastasis, but tumor metastasis occurred. Hence, we aimed to study the DEGs in HRNM_T versus HRNM_N, LRYM_T versus LRYM_N, and HRNM_T versus LRYM_T and clarify the differences.

The abnormality of cell proliferation is generally due to the abnormal expression of related genes, which gives rise to the occurrence and development of tumors [11,12]. In this study, we first identified the DEGs in HRNM_T versus HRNM_N, LRYM_T versus LRYM_N, and HRNM_T versus LRYM_T. The results are shown in Table 1. For HRNM_T versus HRNM_N, there were a total of 224 DEGs. There were a total 504 DEGs in LRYM_T versus LRYM_N, and a total of 88 DEGs in LRYM_T versus LRYM_N. DEGs in the three groups were then analyzed by GO (biological process, BP; cellular component, CC; molecular function, MF) and KEGG functional enrichment. For HRNM_T versus HRNM_N, the DEGs were up-regulated mainly in the cell cycle, the IL-17 signaling pathway and the progesterone-mediated oocyte maturation; and DEGs were down-regulated mainly in the IL-17 signaling pathway. Cdc6, Cdc25A, CDK1, and CycB were up-regulated; Cdc6, Cdc25A, CDK1, and CycB are known to be involved in the cell cycle [13,14], and up-regulation of Cdc6, Cdc25A, CDK1, and CycB can result in tumor initiation and progression. Numerous factors contribute to the progression of breast cancer, including inflammation. The IL17 family consist of six protein members, among them IL17B and its receptor; and the IL17RB signaling pathway plays a key role in development and progression of breast cancer [15,16]. Eotaxin, CXCL10, S100A7, S100A8, S100A9, MMP1, and MMP13 were up-regulated and IL-17B was down-regulated. For LRYM_T versus LRYM_N, DEGs were up-regulated mainly in protein digestion and absorption and cytokine-cytokine receptor interaction; and DEGs were down-regulated mainly in ECM-receptor interaction. CXCL1, CXCL3, and CXCL5 were up-regulated; CXCL1, CXCL3, and CXCL5 are known to contribute to the invasion and metastasis [17,18]. Collagen, THBS, and fibronectin are known to contribute to ECM-receptor interaction, which is down-regulated. ECM is an important organizational barrier of tumor metastasis [19-21]. For HRNM_T versus LRYM_T, DEGs were up-regulated mainly in the PPAR signaling pathway and the adipocytokine signaling pathway. Peroxisome proliferator-activated receptor (PPAR) signaling pathways are nuclear hormone receptors that are activated by fatty acids and their derivatives and may be an important predictor of breast cancer response to neoadjuvant chemotherapy [22]. Expression of the nuclear receptor peroxisome proliferator activated receptor delta (PPARdelta) in breast cancer cells is negatively associated with patient survival [23]. ADIPO has been shown to be related to PPAR signaling pathway and AMPK signaling pathway, and up-regulated [24]. Cancer growth and metastasis depends on the availability of energy. Energy-sensing systems are critical in maintaining a balance between the energy supply and utilization of energy for tumor growth. A central regulator in this process is AMPK [24]. In our study, the DEGs were screened to construct PPI networks and significant modules, and named subnetworks were identified. The functional analysis and pathways enrichment showed that DEGs in the subnetworks were enriched in pathways.

Conclusions

This study analyzed the DEGs in HRNM_T versus HRNM_N, LRYM_T versus LRYM_N, and HRNM_T versus LRYM_T. The GO and KEGG functional enrichments of HRNM_T versus HRNM_N, and LRYM_T versus LRYM_N were consistent with earlier studies. For HRNM_T versus LRYM_T, DEGs were up-regulated mainly in the PPAR signaling pathway; DEGs were down-regulated mainly in the adipocytokine signaling pathway.