High expression of CCNB1 driven by ncRNAs is associated with a poor prognosis and tumor immune infiltration in breast cancer

Data collection

GSE4922 was acquired from the GEO database, which includes gene expression and clinical data on 347 BC patients. The DEGs of BC samples are compared to normal samples using R’s limma software package. The following criteria were used to define DEGs: adjusted P < 0.05 and | logfc | > 0.1. For Pan-cancer analysis, 33 kinds of tumor and normal samples data from TCGA (https://cancergenome.nih.gov/) and GTEx database (https://commonfund.nih.gov/GTEx) were collected using R's ggpubr package (version 3.6.3). The PAM50 procedure was used to classify BC and normal samples into subtypes HER2, Basal like, luminal A, luminal B, and normal [9, 10].

Pan-cancer expression of the target gene

The expression of the target gene was extracted from the transcriptome data of 33 distinct types of cancer and normal tissues in TCGA database to find out which cancers and their corresponding normal tissues besides cancer had a significant differential expression of the target gene.

Analysis of TMB, DNA methylation, immune cell infiltration and immunohistochemistry

CCNB1 was compared with 33 distinct types of cancer data in the downloaded TCGA and GTEx data to explore the relationship between the CCNB1 and TMB, and the connection between the CCNB1 and immune cells in BC tissue. The RNAseq data in the medium level 3 HTSeq-FPKM format and the illumina human DNA methylation 450 DNA methylation data and BC DNA methylation data were acquired from TCGA (Spearman correlation coefficient r = −0.140, P < 0.001 as the cut-off criterion). The expression levels of CCNB1 in normal breast tissues and BC tissues were validated utilizing Human Protein Atlas (THAP) database (https://www.proteinatlas.org/).

Survival analysis filter, clinical correlation analysis and independent prognosis analysis

The Kaplan-Meier (KM) and Cox regression method were used to filter, brush, and identify genes with significant differences in this survival evaluation. The set filtering criteria must satisfy KM < 0.05 and Cox value < 0.05. HR > 1 suggested that the gene was upregulated in BC, while HR < 1 suggested that the gene was downregulated in BC. After that, the clinical correlations of these genes with ER, lymph node, p53, age, and tumor size were analyzed, and the genes with the highest connection being chosen for further investigation. The survival time and survival state were extracted, and the gene was visualized by independent prognosis analysis survival analysis.

Survival analysis

Survival analysis and expression validation of CCNB1 using KM analysis combined with Cox regression were carried out by using the survival and survminer package by R (version 3.6.3).

Identification of DEGs, function and pathway enrichment analysis

Based on CCNB1 expression, BC cases were separated into low and high groups, with the median expression serving as the cut-off value. Specifically, BC cases were placed in the high group if their CCNB1 expression was higher than the median, and in the low group if their CCNB1 expression was lower than the median. The DEGs were identified by using P < 0.05 and | logfc | > 0.1 as the cut-off criterion. High expression samples and low expression samples were used to divide the target genes into two categories.

Cluster profiler

The Gene Ontology (GO) pathway analysis of DEGs were identified using the R’clusterprofile package. The KEGG pathway analysis of DEGs was predicted using the web program Metascape (https://metascape.org/gp/index.html#/main/step1). The KEGG signaling pathway analysis was carried out by using internet program KEGG (https://www.genome.jp/kegg/) [11].

Gene mapping and PPI network construction, module analysis of target gene

The first 20 most significant genes were analyzed according to their interaction. The PPI information was then evaluated using string (search tool) that was used to retrieve the interaction genes. We constructed a PPI network of DEGs by string and applied a confidence score ≥ 0.4 as the cut-off criteria to examine the relevant link between these DEGs. The PPI subnet centered on CCNB1 is then graphically generated using the Cytoscape program.

Candidate miRNA and lncRNA prediction

To find CCNB1 upstream binding miRNAs, we applied miRWalk and miRMap. For the following studies, only the forecasted miRNAs that occurred often in the two algorithms indicated above were used. These projected miRNAs were considered CCNB1 candidate miRNAs.

StarBase (http://starbase.sysu.edu.cn/) is a database dedicated to miRNA research [11]. In BC, starBase was set up to do miRNA-CCNB1, lncRNA-miR-139-5p, or lncRNA-CCNB1 expression correlation analysis. StarBase was used to explore at the expression of miR-139-5p in BC and normal controls. In addition, starBase was utilized to identify candidate lncRNAs that could bind to miR-139-5p. The association between lncRNA expression and CCNB1 was assessed by Spearman correlation.

Cell culture

MCF-7 (Luminal A), BT474 (Luminal B), SKBR3 (Her2+, i.e., Her2 enriched), MDA-MB-231 (Basal like) and MCF 10A cell lines (Normal) [12] were from the Chinese Academy of Science (Shanghai, China). The three cell lines were cultured at 37°C and 5% CO² in a sterile cell culture incubator. The cell lines except MCF-7 were cultured using MEM (Procell, Wuhan, China) containing fetal bovine serum (FBS, Procell, Wuhan, China) and penicillin (100 U/mL) and streptomycin (100 μg/mL) (P/S, Procell, Wuhan, China) to make up a concentration of 10% FBS and 1% P/S. MCF-10A was cultured with the following medium: DMEM+5% HS+20 ng/ml EGF+0.5 μg/ml hydrocortisone+10 μg/ml insulin+1% NEAA+1% P/S.

qRT-PCR

Total RNA was isolated using RNA-easy Isolation Reagent (Vazyme), and cDNA was generated with HiScript^® III SuperMix (Vazyme) according to the kit’s instructions. qRT-PCR was performed in a StepOnePlus (Applied Biosystems) device (Vazyme) using ChamQTM SYBR^® qPCR Master Mix. The specific primers used for qRT-PCR are listed in Table 1. Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as an internal reference gene. CCNB1 was compared and analyzed by 2−ΔΔct method.

Availability of data and materials

All data and materials used in this work were publicly available and also available based on request.

CCNB1 expression in pan-cancer

The transcriptome data of 33 types of cancer and normal tissues in the TCGA database was used to extract CCNB1 expression. CCNB1 has been discovered to be strongly expressed in 19 cancer types, including BC, while low expressed in normal tissues (Figure 1A).

Figure 1. Expression of CCNB1 in normal breast tissue and tumor tissue. (A) Expression of CCNB1 in pan cancers. CCNB1 expression was shown to be higher in cancer tissues than in normal tissues in the TCGA study. (B) TMB of CCNB1 in pan cancer. A positive association with cancer is shown by a red dot bigger than 0; a negative correlation with cancer is indicated by a red dot less than 0. (C) The THPA website was used to examine the expression of the CCNB1 protein in BC specimens and noncancerous breast tissue. Three images of malignant and non-cancerous breast tissues were shown. Abbreviations: N: normal; T: tumor; TCGA: The Cancer Genome Atlas; THPA: The Human Protein Atlas; TMB: tumor mutation burden. ^*P < 0.05, ^**P < 0.01, ^***P < 0.001.

Analysis of TMB, immunohistochemistry and immune cell infiltration of CCNB1

In BC, CCNB1 expression was positively correlated with TMB (Figure 1B). In the analysis of immune cell infiltration, the expression of CCNB1 was positively correlated with the expression of Macrophages M0, Macrophages M1, Dendritic cells activated, T cells follicular helper, T cells CD4 memory activated, it was negatively correlated with the expression of B cells naive, Mast cells resting, T cells CD4 memory resting, Monocytes (Figure 2). Eventually, the upregulation of CCNB1 in BC tumor tissues was confirmed by immunohistochemistry in the Human Protein Atlas (THPA) database (Figure 1C).

Figure 2. The correlation between CCNB1 and immune cell infiltration. Each point (in black) in the figure represents a sample. The abscissa is the content of an immune cell, and the expression of a gene in the ordinate. R > 0 suggests that the gene is positively correlated with the content of an immune cell, R < 0 suggests that the gene is negatively correlated with the content of an immune cell.

Survival analysis

The CCNB1 high expression group had a substantially lower survival time (Including: OS, overall survival; RFS, recurrence free survival; DMFS, distant metastasis free survival) than the CCNB1 low expression group, demonstrating that CCNB1 expression in BC was adversely linked with survival time (Figure 3A). CCNB1 was statistically significant (P < 0.001) in both univariate and multivariate analyses of variables linked to survival, indicating that it may be considered as an independent prognostic factor for BC (Figure 3B).

Figure 3. Survival and clinical association of CCNB1 as an independent prognostic factor in BC. (A) The survival analysis was performed for OS, RFS and DMFS. The CCNB1 high expression group had a considerably shorter survival duration than the CCNB1 low expression group. (B) Univariate and multivariate analyses of factors related to survival. (C) Clinical characteristics associated with CCNB1. It can be seen that CCNB1 is highly expressed in ER (−), lymphnode (+), and p53 (+) groups. Abbreviations: OS: overall survival; RFS: recurrence free survival; DMFS: distant metastasis free survival.

Independent prognostic analysis, clinical correlation analysis and single gene selection

BRIP1, CCNB1, CDT1, CENPW, CSRNP3, DIAPH3, ENKD1, LINC00472, LOC101928767, SPC24, SYCP3 are the genes with the significant difference, according to survival analysis filtering (Table 2). Analysis of the clinical associations of these genes with ER, lymph node, p53, age, and tumor size revealed a high correlation between the CCNB1 with ER, lymph node, p53, and tumor size (Table 3). It can be seen that CCNB1 is highly expressed in ER (−), lymphnode (+), and p53 (+) groups (Figure 3C).

CCNB1 upstream miRNA prediction and analysis

The role of ncRNAs in gene regulation has increasingly been known. To determine whether CCNB1 was influenced by certain ncRNAs, we first anticipated upstream miRNAs that may potentially bind to CCNB1 and eventually discovered 24 miRNAs. There should be a negative association between miRNA and CCNB1 based on miRNA’s action mechanism in the control of target gene expression. In BC, CCNB1 had the strongest negative correlation with miR-139-5p. The expression of miR-139-5p in BC was assessed, as well as its prognostic value. MiR-139-5p was significantly downregulated in BC, as shown in Figure 4B, 4C, and its overexpression was associated to patients’ prognosis. All of these data point to miR-139-5p as the most effective CCNB1 regulating miRNA in BC.

Figure 4. Identification of miR-139-5p as a potential upstream miRNA of CCNB1, and AC026401.3 as a potential upstream lncRNA of miR-139-5p. (A) The expression correlation between predicted miRNAs and CCNB1 in BC. (B) The expression of miR-139-5p in BC and control normal samples was determined. (C) The prognostic value of miR-139-5p in BC assessed plotter. (D) The expression of AC026401.3 in BC and control normal samples was determined. (E) Between CCNB1 and AC026401.3, there was a substantial positive association.

Upstream lncRNAs of AC026401.3 prediction and analysis

By using starBase database, the upstream lncRNAs of miR-139-5p were projected. A total of 75 lncRNAs have been predicted. GEPIA was used to assess the expression levels of these lncRNAs in BC. Only RP11-553L6.5 was highly elevated and showed prognostic significance in BC when compared to normal controls, as well as being strongly related with CCNB1. According to the competing endogenous RNA (ceRNA) theory, lncRNA can improve mRNA expression by binding to shared miRNAs in a competitive manner. As a result, lncRNA and miRNA should have a negative connection, while lncRNA and mRNA should have a positive association. As suggested in Figure 4, compared with normal breast tissue, RP11-553L6.5 was highly expressed in BC (Figure 4D) and had a significant positive correlation with CCNB1 (Figure 4E).

RP11-553L6.5 may be the most potential upstream lncRNAs of the miR-139-5p/CCNB1 axis in BC, based on expression and correlation analyses.

TMB and DNA methylation of CCNB1

We then calculated the BC TMB associated with CCNB1 for each patient in the TCGA and GTEx databases, and we found that the expression of CCNB1 in luminal A BC cell lines was positively correlated with TMB, P < 0.05; while in BC cell lines the expression in luminal B, Her2+, and Basal like BC cell line was not statistically significant (Supplementary Figure 1A). Correlation analysis results of DNA methylation indicating that CCNB1 and DNA methylation have a negative correlation (Supplementary Figure 1B).

Expression of CCNB1 in different molecular subtyping of BC

We use TCGA and GTEx database to classify BC and normal cases into luminal A, luminal B, Her2+, basal like and normal groups. We found that the expression of CCNB1 in the four molecular subtyping (luminal A, luminal B, Her2+, basal like) of BC was higher than that in the normal group (Figure 5A).

Figure 5. Expression of CCNB1 in different molecular subtyping of BC. (A) The relative expression of CCNB1 in luminal A, luminal B, Her2+, basal like and MCF 10A cell lines. CCNB1 expression was higher in the four molecular subtyping of BC (luminal A, luminal B, Her2+, basal like) than in the normal group. (B) Compared with the MCF 10A cell line, CCNB1 was up-regulated in luminal A, luminal B and basal-like cell lines, but down-regulated in Her2+ cell line, as verified by qRT-PCR. ^*p < 0.05, ^**p < 0.01, ^***p < 0.001. Abbreviation: BC: breast cancer.

qRT-PCR

The expression of CCNB1 in luminal A (MCF-7), Luminal B (BT474), Her2+ (SKBR3), Basal like (MDA-MB-231) BC cell lines and human normal breast epithelial cell lines (MCF-10A) was verified by qRT-PCR.

The results showed that the expression of CCNB1 in the luminal A BC cell line was higher than that in the human normal breast epithelial cell line (Figure 5B).

Identification of DEGs, function and pathway enrichment analysis

In this study, 347 BC samples were collected from GSE4922 of GEO database. There were 109 DEGs identified, with 58 up-regulated genes and 51 down-regulated genes. DEGs are visualized on the heat map (Figure 6A). DEGs impacted extracellular matrix organization, leukocyte migration, and extracellular structure organization, according to Go analyses (Figure 6B-a). The cell cycle, PPAR signaling route, and TCF signaling beta pathway were the most enriched in KEGG pathways (Figure 6B-b). Among them, CCNB1 was most enriched in the cell cycle, which was involved in the most significantly enriched signaling pathway in BC tumorigenesis and pathogenesis (Figure 7B).

Figure 6. Related genes and signaling pathways affected by CCNB1. (A) DEGS: heat map. According to the different expression levels of CCNB1, they are divided into two groups: high expression (Green in the right half of the image) and low (Pink in the left half of the image) expression. The abscissa is the sample, and the ordinate is the gene name. (B) Go and KEGG analysis. (a) Go analysis demonstrated that CCNB1 and genes significantly associated with CCNB1 alterations mainly enriched extracellular matrix organization, leukocyte migration and extracellular structure organization. (b) The major enriched pathways were those connected with cell cycle, AMPK, PPAR signaling pathway and TGF signaling beta pathway according to the KEGG analysis. (C) The genes that are strongly linked to CCNB1.

Figure 7. PPI network and KEGG pathway analysis. (A) PPI network of the significant genes in BC. (a) The PPI network was constructed using 40 genes selected from the string web database. (b) Construction of CCNB1 subnetwork of PPI. (B) KEGG pathway analysis; CCNB1 are involved in the cell cycle. Abbreviations: DEGs: differentially expressed genes; GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; PPI: protein-protein interaction.

Gene mapping and PPI network construction, module analysis of target gene

The interaction between the top 20 genes and the bottom 20 genes was used to perform correlation analysis (Figure 6C). Positive correlation genes with CCNB1: ZNF367, UHRF1, FBXO5, CDCA5, KIF14, ANLN, SPC24, DEPDC1, UBE2T, NUF2, ASPM, MAD2L1, HAUS8, MND1. Negative correlation genes with CCNB1: SCNN1G, IQCA1, PRO1082, ANKRD30A, LINC00993, DNALI1, CAPN8, NEK10, MGP, TEKT3, RP11−250B2.3, FGF9, PCDHB5, DIO3OS, OSR1, SCN4B, ADCY4, ADAM23, SLC28A3, CYP4Z1, AK5. To explore the highly likely correlation between these DEGs, we built a PPI network of DEGs via string and set the confidence score ≥ 0.4 (Figure 7A-a). Finally, we depict the PPI sub network centered on CCNB1 by using Cytoscape (Figure 7A-b).