Role of Bclaf1 in Promoting Adrenocortical Carcinoma Proliferation: A Study Combining the Use of Bioinformatics and Molecular Events

Introduction

The available treatments for adrenocortical carcinoma (ACC) are limited, and among them, surgical resection is basically the only curative option.^1,2 Nevertheless, up to 85% of patients who undergo radical resection exhibit poor outcomes, such as local disease recurrence and concurrent metastases.^3–6 Not surprisingly, the prognosis of ACC is poor; the 5-year overall survival (OS) rate is only 38%.⁷

Therefore, the accurate identification of adrenal carcinoma at an early stage of the disease^8,9 and understanding the mechanism of ACC are particularly urgent for exploring new drug targets.

The Oncomine database indicates that Bcl-2-associated transcription factor 1 (Bclaf1) mRNA is more highly expressed in ACC than in normal tissues. The Cancer Genome Atlas (TCGA) provides a gene map of human cancer and contains an abundant amount of complete clinical follow-up information.^10,11 Our study aimed to identify core genes associated with BCLAF1 by analysing the TCGA database. G2/mitotic-specific cyclin B1 (Cyclin B1) and cyclin-dependent kinase 1 (CDK1) play important roles in controlling the transition of cells from the G2 phase to the M phase, which is a key event in the initiation of mitosis,^12,13 and their disorder can lead to uncontrolled cell proliferation and contribute to the development of tumours.¹⁴ Previous studies have explored CDK1 and Cyclin B1 as drug targets.¹⁵

Bclaf1 is a powerful molecule that influences mRNA transcription and splicing to participate in tumour progression^16,17 and might be a new target of CDK1 inhibitors. In our previous research, we found that Bclaf1 promotes the proliferation and clonal formation of hepatocellular carcinoma (HCC) cell lines.¹⁸ Thus, the mechanism through which Bclaf1 affects ACC remains unclear. Based on an analysis of RNA sequencing (RNA-seq) data of TCGA-ACC, we assessed the impact of Bclaf1 on ACC and explored its molecular mechanism.

Materials and Methods

Data Sources

Data retrieved from multiple research centres were used for the integrated analysis performed in this study. The Giordano Adrenal and Giordano Adrenal 2 datasets were selected from the Oncomine database¹⁹ (https://www.oncomine.org). The RNA-seq data of patients with ACC from the TCGA database, including data from 76 tumour tissues and the corresponding clinical information, were derived from the Genomic Data Commons (GDC) API (https://portal.gdc.cancer.gov), and samples with duplicated TCGA barcodes and incomplete survival data were excluded. GSE10927 and GSE90713 were then downloaded from the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database. We used Gene Expression Profiling Interactive Analysis (GEPIA)²⁰ to conduct a series of online analyses, including box plots and survival analysis. ChIP-seq data were acquired from the Cistrome Data Browser (http://cistrome.org).

Gene Expression Analysis

For the RNA-seq data, the expression levels of mRNA were computed as fragments per kilobase per million (FPKM). After transcripts per million (TPM) transformation, the median value of BCLAF1 expression was used to divide the patients into BCLAF1^high and BCLAF1^low groups. Using R language (limma package, R version 3.6.1), we identified key differentially expressed genes (DEGs) between the BCLAF1^high and BCLAF1^low groups. Genes with abs(log2(fold-change)) > 1 (absolute value) and a false discovery rate (FDR) < 0.01 were selected as candidate signatures for further analysis. Volcano plots were generated based on the DEGs and drawn using the ggplot2 and pheatmap R packages.

Weighted Gene Co-Expression Network Analysis (WGCNA) and Protein-Protein Interaction (PPI) Network

Gene co-expression network analysis was specifically performed using the R package WGCNA.²¹ The screened genes were uploaded to the Search Tool of the Retrieval of Interacting Genes/Proteins (STRING) database (http://string-db.org) to construct a PPI network with an upper confidence threshold of 0.400. The results from the STRING analysis were visualized in a degree-attributed circle layout using Cytoscape v.3.7.2 (http://cytoscape.org/).

Gene Set Enrichment Analysis and Receiver Operating Characteristic (ROC) Curve Analysis

WebGestalt²² (http://www.webgestalt.org/) was utilized for the GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of genes. ROC curves were analysed by employing the ggplot2 R package, and a Cox proportional hazards model was used to predict the 1–5-year survival states.

Gene Expression Analysis

For the RNA-seq data, the expression levels of mRNA were computed as fragments per kilobase per million (FPKM). After transcripts per million (TPM) transformation, the median value of BCLAF1 expression was used to divide the pa tients into BCLAF1^high and BCLAF1^low groups. We applied R language (limma package, R version 3.6.1) to identify key differentially expressed genes (DEGs) between the BCLAF1^high and BCLAF1^low groups. Genes with abs(log2(fold-change)) > 1 (absolute value) and an FDR < 0.01 were selected as candidate signatures for further analysis. Volcano plots were generated based on the DEGs and drawn using the ggplot2 and pheatmap R packages.

Cell Culture

SW-13 (RRID: CVCL_0542) and NCI-H295R (RRID: CVCL_0458) cell lines were purchased from Center for Excellence in Molecular Cell Science, CAS (Shanghai, China). SW-13 were maintained in L-15 medium supplemented with 10% foetal bovine serum (Gibco, USA), 1 mM glutamine, and 100 U/mL penicillin-streptomycin (Solarbio, Beijing, China). NCI-H295R were maintained in special medium for NCI-H295R cells. (Procell, Wuhan, China). The cell lines were authenticated by STR profiling.

Oligonucleotide and Plasmid Transfection and Stable Cell Line Generation

Effective siRNA oligonucleotides that target Bclaf1 were purchased from GenePharma (Suzhou, China). Full-length Bclaf1-overexpressing plasmids were kindly provided by Dr. Tang (College of Veterinary Medicine, China Agricultural University). SW-13 cells were transfected with oligonucleotides or plasmids using Lipofectamine 3000 (Invitrogen, CA, USA) and were harvested for assays 48 h after transfection. The synthesized short hairpin RNAs targeting Bclaf1 (shBclaf1) or nonspecific control RNAs (shNC) were designed and cloned into LV3 lentiviral vectors, which were constructed by GenePharma (Suzhou, China). The siRNA and shRNA sequences are listed in Table S1. SW-13 cell lines stably expressing shBclaf1 or the vector were generated after infection with the lentiviral particles and were selected with 0.8 μg/mL puromycin (MP Biomedicals, USA).

Western Blotting (WB)

The cells were lysed with lysate (KeyGEN, China; KGP2100), and proteins were extracted to measure the protein concentrations with Bradford reagent (Beyotime, China; P0006C). At least 50 μg of sample was then separated in 10% SDS-PAGE gels and transferred onto PVDF membranes (Merck Millipore, Billerica, MA, USA). The antibodies used in the Western blotting assay are detailed in Table S2.

RNA Extraction, Reverse Transcription and Quantitative PCR

RNAiso (Takara, Tokyo, Japan) was used to extract total RNA from SW-13 cells. Total RNA was used for reverse transcription with the Prime Script RT Reagent Kit (Takara). Then, real-time quantitative PCR was performed using SYBR Premix Ex Taq (Takara) on the CFX Maestro sequence detection system (Bio-Rad, Hercules, CA, USA). The relative expression levels were calculated using the 2-ΔΔCT method based on the normalization control β-actin for the SW-13 cell line. The primers for amplification of human genes are presented in Table S1.

Cell Proliferation

To measure cell proliferation, we conducted a series of cell proliferation experiments, including Cell Counting Kit-8 (Dojindo, Tabaru, Mashikimachi, Japan) assays, crystal violet staining, and EdU assays. The experimental methods are described below. After transfection, cells were incubated for 2 h at 37 °C with Cell Counting Kit-8 solution, and the absorbance at 450 nm was measured. For crystal violet staining, 5 × 10³ cells in 2 mL of medium were plated in six-well plates and cultured for 5 days. The cells were then fixed with methanol, stained with 0.1% crystal violet for 15 min and counted. For the EdU assay, immunofluorescence (IF) staining was performed with an EdU kit (RiboBio, Guangzhou, China) according to the manufacturer’s protocol. The results were then quantified with Image-Pro Plus 6.0 software.

Flow Cytometry

To analyse the cell cycle, the cells were trypsinized and fixed in 70% cold ethanol at −20 °C. Prior to flow cytometry, the cells were washed three times in phosphate-buffered saline (PBS) and treated with a final concentration of 0.2 mg/mL RNaseA and 50 μg/mL propidium iodide (PI) for 30 min at 37 °C while protected from light. Flow cytometry was performed with an Accuri C6 Plus flow cytometer (Becton, Dickinson and Company, USA).

Immunohistochemistry (IHC)

Paraffin sections of tissues were subjected to IHC according to the standard protocol (Huayin Health, Guangzhou, China) using primary antibodies against CDK1 (1:100) and Cyclin B1 (1:200). The normalized expression of CDK1 and Cyclin B1 (integrated optical density [IOD]) in tissues was statistically analysed using Image-Pro Plus 6.0.

Statistical Analysis and Other Software

All of the data are presented as the means ± standard deviations (SDs) from at least three independent experiments, and representative results are presented. The differences between groups were analysed by one-way ANOVA, Dunnett’s test, or Student’s t-test, and a simple correlation analysis was performed; all the analyses were performed with SPSS software version 22.0 (IBM, Chicago, IL, USA) and GraphPad Prism 7 software according to the data features. The significance level was set to p < 0.05. To identify the potential function of Bclaf1, GSEA was performed using GSEA v4.0.3 software.

Results

WGCNA Revealed That the Turquoise Gene Modules Were Related to Tumour Prognosis and the Expression of Bclaf1

Current situation should be considered with caution, as several of the patients were diagnosed with late-stage adrenal adenocarcinoma with distant metastasis in our hospital. Unfortunately, they had lost their best chance of surgical treatment by the time of diagnosis. Therefore, we aimed to explore this rare and highly malignant tumour. In the Giordano Adrenal 2 dataset from the Oncomine database, BCLAF1 was more highly expressed in ACC tissues (29 cases) than in normal tissues (Figure S1A). OS or disease-free survival (DFS) analysis performed using the GEPIA database revealed that high BCLAF1 expression in patients was significantly correlated with poor prognosis (Figure S1B and C). To further explore the role of Bclaf1, we used the median value of BCLAF1 expression obtained from an analysis of the TCGA dataset and divided the patients into BCLAF1^high and BCLAF1^low groups. Atotal of 2323 differentially expressed genes (DEGs) were identified from the comparison of the BCLAF1^high group with the BCLAF1^low group (Figure S1D). We then used the DEGs to construct gene co-expression networks using WGCNA based on complete clinical information from the TCGA database. After performing clustering, we excluded outliers (Figure S1E) and obtained 65 ACC samples. Several clinical traits of each sample are shown in the heat map, including age, sex, TNM stage classification, OS, and progress free survival (PFS) (Figure1A). To explore the relationship between BCLAF1 and each module, we used the mRNA expression level of BCLAF1 as characterization information to construct aWGCNA. Six was subsequently selected as the soft thresholding power (Figure S1F), and the correlated modules (r < 0.25) were then merged (Figure S1G and H). According to their co-expression relationships, we identified 14 modules and revealed their relevant clinical traits (Figure1B). The correlation coefficients (r) and pvalues are illustrated in aheatmap. Anegative correlation between turquoise gene modules and survival time (r=−0.46) indicated that the turquoise gene modules were strongly correlated with the prognosis of ACC. We found that BCLAF1 is closely related to the turquoise gene modules (r=0.53), which indicates that BCLAF1 plays an important role in the genetic mapping changes that lead to the prognosis of ACC.

Figure 1 Identification of gene module related to Bclaf1. (A) Clustering of samples based on the gene expression (upper panel) and baseline characteristics of each sample was described by a heat map (lower panels). Deeper one refers to higher levels of clinical traits. The TNM staging system classifies tumours by using the tumour (T), lymph node (N), and metastasis (M) attributes, where “T” refers to the extent of the primary tumour (T) and “N” and “M” describe the absence or presence of lymph node(s) and distant metastasis, respectively. (B) Heatmaps of the correlations between different gene modules and clinical traits. Each cell contains the corresponding correlation and p-value. The colour of the right axis depicts the correlation coefficients of each cell.

Bclaf1 Plays an Important Role in the Proliferation of ACC

To further underscore the biological role of Bclaf1 in ACC, two cell lines were used as an ACC model for the in vitro experiments. We knocked down (siBclaf1) and overexpressed Bclaf1 with a plasmid expressing the Bclaf1-encoding gene (FL) in cells using cell transfection techniques. The transfection efficiency was determined by Western blotting (Figure 2A and B). The downregulation of Bclaf1 by siRNA transfection reduced the proliferative capacity of NCI-295R cells in comparison with that of siRNA control cells (Figure S1I). For SW-13 cells, cell lines with stable knockdown of Bclaf1 were constructed (Figure S2A). In addition, the effect of si Bclaf1 and sh-Bclaf1 on SW-13 cell proliferation was the same as that on NCI-295R (Figure S2B). The inhibition of Bclaf1 expression reduced the growth rate and DNA replication activity of NCI-295R cells, whereas transient transfection with a plasmid expressing the Bclaf1-encoding gene (FL) yielded the opposite results (Figure 2C and D). Consistent data were obtained for the SW-13 cell line (Figure 2E and F). These results indicate that Bclaf1 might be an important molecule that promote tumour proliferation of ACC, as shown in Figure 2.

Figure 2 Effect of Bclaf1 on the proliferation of ACC cells. (A and B) Bclaf1 was downregulated and overexpressed by cell transfection. (C and E) The proliferation of NCI-H295R and SW-13 cells transfected as described above was determined by CCK-8 assay. (D and F) Cell proliferation following transfection was assessed by EdU IF staining. *p < 0.05, **p < 0.01, and ***p < 0.001 vs the control.

Bclaf1 Contributes to Cell Cycle Progression in ACC

To determine the functions related to the DEGs, we performed an enrichment analysis using WebGestalt. A KEGG analysis revealed that BCLAF1^high-specific genes were mainly involved in signalling pathways such as the cell cycle (Figure 3A). Notably, a GO analysis demonstrated that the mitotic cell cycle phase transition BP (Figure 3B) was significantly associated with the genes identified from the BCLAF1^high samples. Furthermore, we performed GSEA to detect pathways connected with BCLAF1. Various gene sets were found to be upregulated in samples from patients with high BCLAF1 expression, and these included genes related to the mitotic spindle and G2/M checkpoint (Figure 3C–F). The results from the analysis of the GSE10927 dataset also showed that the expression level of BCLAF1 was correlated with the mitotic rate of tumours (Figure S2C). To further confirm the results from the bioinformatics analysis, we performed a cell cycle assay based on flow cytometry. More NCI-295R cells were arrested at the G2/M phase after cotreatment with siRNA transfection than after control treatment, whereas overexpressed Bclaf1 had the opposite effect (Figure 3G). Consistent data were obtained for the SW-13 cell line (Figure 3H). Based on the above-mentioned results, we surprisingly found that BCLAF1 in ACC is positively correlated with the cell cycle.

Figure 3 Bclaf1 is involved in cell cycle progression. (A and B) KEGG (A) pathway enrichment analysis and GO biological process analysis (B) of BCLAF1high- and BCLAF1low-specific mRNA expression. Screening was conducted at FDR < 0.05. (C–F) GSEA plot; several key pathways, including the mitotic spindle and G2/M checkpoint, which are related to the cell cycle, were visualized. Parameters: p-value < 0.05, FDR q-value < 0.25. (G and H) Images of the cell cycle results (left panels) and stacked bar charts (right panels) of the cell cycle phase distribution of SW-13 cells. *p < 0.05 and **p < 0.01 vs the control.

CCNB1 and CDK1 are the Key Genes in the Turquoise Module That Correlate with Bclaf1

We then paid close attention to the meaningful clinical information (r > 0.4) in Figure 1B, including pathological tumour, OS (in months), and disease-free survival (PFS, in months). After overlapping these factors (Figure S2D), we identified 53 genes with the following criteria: module membership > 0.8, gene significance of pathological tumour > 0.2, and p value < 0.01. We further explored the prognostic value of these 53 genes by combining them to predict their impact on OS using R. We used Cox proportional hazards regression to analyse the risk of the patients based on the expression of the above-mentioned 53 genes and then divided the patients into high- and low-risk groups. The relationship between gene expression and patient risk score is plotted in Figure S2E. As seen from the figure, the high-risk group was on the right side of the figure and correspondingly had poor survival status and high gene expression. An ROC curve analysis revealed that these 53 genes were valuable with an AUC value of 1 to predict the survival status of patients after years of follow-up (Figure S2F). The survival curve indicated that high-risk patients had a poor prognosis (Figure 2G). These data demonstrated that BCLAF1^high-specific genes are associated with the prognosis of ACC.

To identify highly connected hub genes among these 53 genes, we also constructed a PPI network of these 53 genes and sorted them by degree to evaluate the number of edges to which nodes were connected. CCNB1 was identified as the hub gene with the highest degree (Figure S3A), and CDK1 exhibited a high degree. This result was consistent with the scatter diagram (Figure 4A–C), in which we focused on the genes in the top right corner of the plot because these were correlated with both clinical traits and module membership. CCNB1 has been reported to affect the cell cycle together with CDK1.^23,24 Importantly, OS or DFS analysis performed using the GEPIA database revealed that high expression of these two genes in patients was significantly correlated with poor prognosis (Figure S3B and C). To compare gene expression in different pathological types, we performed differential analysis by GEPIA, as shown in Figure S3D and E. CDK1 and CCNB1 were more highly expressed in ACC tissues than in normal tissues. We then used the Oncomine database for further verification (Figure S3F and G). Microarray-based IHC data, which included 10 ACC specimens, eight ACA specimens and four normal specimens, also confirmed higher levels of these two proteins in ACC tissues than in ACA and normal tissues, as illustrated in Figure 4D. These data suggest that CDK1 and CCNB1 are the hub genes among the BCLAF1^high-specific genes and is highly expressed in tumour tissue.

Figure 4 CCNB1 and CDK1 are the key genes in the turquoise module. (A–C) Significance of the genes in the turquoise module for three clinical characteristics. The significant genes were identified based on module membership > 0.8 and gene significance > 0.2. (D) Representative images of IHC staining of CDK1 and Cyclin B1 in different types of tissue specimens; the scale represents 100 μm (left panel). Quantification of protein levels according to IHC integrated optical density in ACC, ACA and normal tissues; *p < 0.05 and **p < 0.01 (right panel).

Bclaf1 Promotes the Proliferation of ACC by Influencing CDK1 and Cyclin B1

Based on the above-described studies, we focused on two hub genes, CDK1 and CCNB1, which can play important roles in ACC progression and prognosis. By analysing the data from the TCGA and GEO databases, we found that BCLAF1 was linearly correlated with CDK1 and CCNB1 (Figure 5A–C). Consistent with the aforementioned results, we wondered whether Bclaf1 is involved in the cell cycle biological process of ACC by regulating CDK1 and Cyclin B1. We then verified our conjecture through molecular biology experiments with NCI-H295R and SW-13 cells. Consistent with our hypothesis, the abundance of CDK1 and CCNB1 mRNA in SW-13 cells significantly decreased after Bclaf1 downregulation (Figure S3H and I), which suggested that Bclaf1 acts on these at the mRNA level. Interestingly, the CDK1 and Cyclin B1 protein levels showed consistent changes: the downregulation of Bclaf1 protein levels by siRNA transfection reduced the protein levels of CDK1 and Cyclin B1, and the overexpression of Bclaf1 after transfection yielded the reverse effect (Figure 5D). For SW-13 cells, the downregulation of Bclaf1 levels by transient siRNA transfection or stable shRNA transfection reduced the levels of key proteins of the cell cycle (Figure 5E). Taken together, these findings support the hypothesis that Bclaf1 promote tumour proliferation and cell cycle transition by participating in the regulation of CDK1 and Cyclin B1.

Figure 5 Bclaf1 promotes the proliferation of ACC by modulating CDK1 and Cyclin B1. (A–C) Correlation plots showing positive associations between Bclaf1 expression and the expression of CDK1 or CCNB1. Pearson correlation coefficient r values and p-values are shown. (D and E) Regulation of CDK1 and CCNB1 via Bclaf1 knockdown or overexpressed after cell transfection. The protein levels were detected by Western blotting (right panel). *p < 0.05 and **p < 0.01 vs the control.

Discussion

Although ACC is a rare disease, it tends to progress rapidly and relapse easily and is thus associated with a poor prognosis. This study explores strategies for enabling ACC patients to receive accurate effective treatment as early as possible.

It is difficult to conduct prospective cohort studies or cross-sectional studies and even harder to conduct clinical trials because there are a limited number of ACC patients, and they are not easy to recruit. Fortunately, combining complete clinical follow-up information, we built a bridge using bioinformatics methods between macroscopic disease occurrence and microscopic molecular events to begin our study, which is the beginning consistent with the objective clinical reality.

Through a bioinformatics analysis, we found an association between BCLAF1 and CDK1. We identified the critical CDK1 and CCNB1 genes related to the development of ACC, which is in line with other strategic bioinformatics research results.²⁵

The current treatment for ACC is mainly early surgical treatment with mitotane adjuvant therapy to prevent recurrence.^1,2 Radiotherapy combined with chemotherapy is often selected after metastasis for patients with advanced disease,²⁶ and the overall treatment plan is limited, particularly because the effect of treatment on patients with advanced disease is not good. CDK1 is a key protein regulating the cell cycle, whereas cyclin B1 binds to CDK1 to drive mitosis and promote cell passage through the M phase.^23,24 Therefore, CDK1 has long been considered a good target for drug development against tumours and other proliferative disorders. Present studies show that CDK1 is part of the extensive cyclin protein family, and the selectivity and specificity of ATP-competitive CDK inhibitors are not ideal.¹⁵ Although some compounds have significant preclinical effects, their development and application are limited. It is thought that improving the selectivity of CDK1 inhibitors through further research might be key to their successful development. Many new approaches for CDK1 inhibitors have emerged recently.²⁷

Interestingly, the experimental data investigated in this study demonstrated that Bclaf1 influences the cell cycle in ACC and participates in the regulation of CDK1 and cyclin B1 specifically through gene regulation. Bclaf1 is a promising prospective drug target for ACC. Furthermore, Bclaf1 is a powerful molecule that can both affect the proliferation of cancer cells and regulate the transcription of Hif1α,¹⁷ as indicated in our previous research. The cell adhesion molecule binding molecular function (Figure S4A) and the chromosomal region cellular component (Figure S4B) were significantly associated with the genes identified in BCLAF1^high samples. Using the Cistrome Data Browser, we found that the binding of Bclaf1 to the promoter regions of CDK1 and CCNB1, particularly that of CDK1, shows strong specificity (Figure S4C and D). Bclaf1 might promptly adjust the transcription of CDK1 and CCNB1. Although this molecule might thus be a new potential therapeutic target of ACC, further work is needed to thoroughly interpret the molecular mechanism. Further confirmation is needed in in vivo and in vitro trials related to treatment.

Conclusions

In conclusion, the disease burden of ACC cannot be ignored. Our study indicated that tumour proliferation and cell cycle transition are influenced by Bclaf1, which has been shown to be involved in the regulation of CDK1 and Cyclin B1.