The mRNA sequencing data from 19, 131 samples were downloaded from TCGA, TARGET, and GTEx projects hosted at the University of California at Santa Cruz (UCSC) Cancer Genome Browser (https://xenabrowser.net/datapages/?cohort=TCGA%20TARGET%20GTEx). Additional data were obtained from the literature (Malta et al., 2018). The data were cleaned, leading to the retention of 38 cancer types. Data on clinical information, copy number variation (CNV), tumor mutation burden (TMB), and single-nucleotide polymorphisms (SNPs) were downloaded from the UCSC database at the same time.

The transcriptome data were converted into log2 (TPM+1) to assess differences between normal and tumor tissues. For protein expression, we downloaded mass spectrometry data from the Clinical Proteomic Tumor Analysis Consortium (CPTAC) database for further analysis (Chandrashekar et al., 2017; Chandrashekar et al., 2022). Z-values represented standard deviations from the median across the given cancer types.

The codes used in this study are available in the GitHub repository (https://github.com/Honglei-Zhou/BOLA2B.git). The genes involved in our study were listed as Appendix A1 at the end of article.

Cox proportional hazard ratio regression analysis was performed with the ‘survival’ package in R to evaluate the dependency of the survival data such as overall survival (OS), disease-specific survival (DSS), and progression-free interval (PFI) on gene expression. The data were divided into groups according to the median level of BOLA2B expression; this was applied to all cancer types. Log-rank tests were used to determine the statistical significance of the hazard ratios. Forest plots and Kaplan–Meier curves were drawn with the R packages “ggplot” and “survminer,” respectively.

Data on mRNA modifications were collected and analyzed. Writer, eraser, and reader genes that were highly associated with the m1A, m5C, and m6A modifications were identified. Pearson’s correlation coefficients were used to determine the associations between the identified genes and BOLA2B expression in various cancers.

The BOLA2B gene expression profiles were extracted separately for each tumor type and were then mapped to the gene symbol, and the infiltration scores of B cells, CD4 T cells, CD8 T cells, neutrophils, macrophages, and dendritic cells (DCs) were determined using the “IOBR” package in R.

The expression profiles of immune-related genes were divided into inhibitory and stimulatory categories and were analyzed using the same process. The corr. test function in the ‘psych’ R package was used to calculate Pearson’s correlation coefficients. The ESTIMATE score for each tumor was calculated using the R package “ESTIMATE.”

The TMB for each tumor was calculated using the “TMB” function in “MFT” tools. The TMB and gene expression data of the samples were integrated into a single matrix, and a log2 (x+0.001) transformation was applied to the expression data.

The BOLA2B CNV data were downloaded from the Genomic Data Commons (GDC) Data Portal. The data were processed using GISTIC software. Sample barcodes were used to combine the clinical survival data with the CNV data. The ‘survival’ package in R was used to fit data on survival time and survival status with the CNV data, and log-rank tests were used to examine differences in survival rates.

Exon sequencing data for nearly 30 tumors are available in TCGA database. Missense_Mutation, Frame_Shift_Del, Non-sense_Mutation, In_Frame_Ins, Splice_Site, In_Frame_Del, and Frame_Shift_Ins were considered as deleterious mutations in this analysis.

Pearson’s correlation coefficients were used to assess interactions between microsatellite instability (MSI) and BOLA2B gene expression. RNA-based stemness scores (RNAss) indicated that the RNA expression-based score drives the main stemness (Malta et al., 2018).

The Search Tool for the Retrieval of Interacting Gene (STRING) database (https://string-db.org/) is a web application that analyzes potential protein–protein interaction networks using multiple algorithms. We used the single protein name ‘BolA2B’ and the organism ‘Homo sapiens’ to obtain a map of predicted interacting proteins. To obtain more reliable results, we set specific parameters with a minimum required interaction score set as ‘medium confidence (0.4)’ and the maximum number of first-shell interactors to ‘not more than 50 interactors’. Functional enrichment of the genes was analyzed using Gene Ontology (GO), including information on biological process, molecular function, and cellular components, and Kyoto Encyclopedia of Genes and Genomes (KEGG).

BOLA2B shRNA1 and shRNA2 were amplified by PCR and subcloned into the pLKO.1 vector. The sequence of BOLA2B shRNA1 was GGC​ACG​TGA​GCG​ACA​GAA​ATG​TTC​AAG​AGA​CAT​TTC​TGT​CGC​TCA​CGT​GCC, and the sequence of BOLA2B shRNA2 was GCG​AGA​AGC​TGC​AGC​GGG​ACC​TTC​AAG​AGA​GGT​CCC​GCT​GCA​GCT​TCT​CGC. The Ctrl group used the same plasmid vector as the experimental group, except for the shRNA sequence. Lentivirus vectors were used to transfect MDA-MB-231 and SKBR3 cells for 24 h, followed by selection with 2 μg/ml puromycin to obtain cells with stable expression.

The human breast cancer cell lines MDA-MB-231 and SKBR3 were cultured in Dulbecco’s modified Eagle medium (DMEM) supplemented with 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin at 37°C with 5% CO₂. For the measurement of cell viability, cells of different groups were seeded in 96-well plates at the same density of 5 × 10³/mL, and the Cell Counting Kit-8 (Beyotime, #C0038) solution was added at the specific time points. The cells were then incubated for 2 h at 37°C, and the absorbance was read at 450 nm. Cell proliferation was measured using the EdU Cell Proliferation Assay Kit (RiboBio, #C10310-1), according to the manufacturer’s instructions. Nuclei were counterstained with 1 μg/ml of DAPI (Beyotime, #C1002). The proportion of cells stained with EdU was examined by fluorescence microscopy (Carl Zeiss AG, Axio Vert. A1).

Cell cycle changes were monitored by flow cytometry. Cells were harvested and fixed with 70% cold ethanol at -20°C overnight. The DNA was stained with propidium iodide (PI), according to the manufacturer’s instruction (Multi Sciences. #CCS01).

Specimens of different cancers were collected from the First Affiliated Hospital of Nanjing Medical University. The project was approved by the Ethics Committee of Jiangsu Province People’s Hospital, with the code number 2020-SR-477. Immunohistochemistry (IHC) was performed using kits from MXB Biotechnologies, #KIT-9710 and #DAB-2031, in accordance with the manufacturer’s instructions. IHC scores for BolA2B staining were determined by light microscopy and calculated as the staining intensity score multiplied by the staining percentage score. The staining intensities were classified as follows: no staining, 0; weak staining, 1; moderate staining, 2; strong staining, 3; and very strong staining, 4. The staining percentages were classified as 0%–10%, 0; 10%–25%, 1; 25%–50%, 2; 50%–75%, 3; and 75%–100%, 4.

Targeted cell proteins were extracted with the RIPA buffer. Proteins of different molecular weights were separated by SDS-PAGE. After electransferring onto the PVDF membrane, 5% skim milk was used to block extra places for 2 h. Then, the membrane was incubated with targeted primary antibodies at 4°C overnight. The antibodies used in the study were against BolA2B (Proteintech, # 26080-1-AP), tubulin (Proteintech, #11224-1-AP), CDK6 (Cell Signaling Technology [CST], #13331), cdc6 (CST, #3387), cyclin B1 (CST, #12231), cyclin D1 (CST, #55506), cyclin E1 (CST, #20808), p21 (CST, #2947), p-AKT (CST, #4060), AKT (CST, #9272), p-mTOR (CST, #5536), mTOR (CST, #2983), p-ERK (CST, #4370), and ERK (CST, #4695). The next day, the membrane was incubated with secondary antibodies (CST, #7074P2). The bands could be observed by Immobilon™ Western Chemiluminescent HRP Substrate (Millipore, United States)

Continuous data were analyzed by t-tests, with paired t-tests used for comparing differences between groups. Categorical data were analyzed using χ2 tests. Pearson’s correlation analysis was used for comparing associations between groups of continuous data. Two-sided p-values less than 0.05 were considered statistically significant, and asterisks were used in some figures to indicate the level of significance: *= (p < 0.05), **= (p < 0.01), and *** = (p < 0.001).

We analyzed BOLA2B mRNA expression between tumor and normal tissues in 33 cancer types. As shown in Figure 1A, BOLA2B was upregulated in BLCA, BRCA, CESC, CHOL, COAD, ESCA, HNSC, KIRC, KIRP, LIHC, LUAD, LUSC, PCPG, PRAD, READ, STAD, and UCEC tumor samples compared to normal tissues. Paired comparisons also verified that BOLA2B was a cancer-related factor (Figure 1B). BLCA, BRCA, CHOL, COAD, HNSC, KIRC, KIRP, LIHC, LUAD, LUSC, PRAD, READ, STAD, THCA, and UCEC tumor samples all showed higher expression compared to their paired adjacent normal tissues. We then analyzed the BolA2B protein levels in different cancers. Analysis of data from the CPTAC database showed that COAD, OV, KIRC, UCEC, LUAD, PAAD, HNSC, GBM, and LIHC tumor tissues all expressed increased levels of BolA2B; the exception was BRCA, as shown in Figure 1C. These results indicated the association of increased BOLA2B expression with multiple tumor types.

We then evaluated the prognostic significance of BOLA2B in pan-cancer analysis. In the overall survival (OS) assessment, the forest plot showed that the expression was significantly correlated with NB (HR = 2.55, 95%CI: 1.80–3.61), GBMLGG (HR = 1.90, 95%CI: 1.47–2.44), UVM (HR = 5.83, 95%CI: 2.67–12.72), KIRC (HR = 1.62, 95%CI: 1.33–1.96), LIHC (HR = 1.27, 95%CI: 1.08–1.49), KIPAN (HR = 1.27, 95%CI: 1.07–1.51), LGG (HR = 1.68, 95%CI: 1.15–2.45), ACC (HR = 2.13, 95%CI: 1.03–4.44), and OV (HR = 0.86, 95%CI: 0.79–0.93) (Figure 2A). Kaplan–Meier curves were used to evaluate differences in OS between the high- and low-BOLA2B subgroups in different types of cancer. As shown in Figure 2B, high expression of BOLA2B predicted better OS in OV patients, while in other cancers, such as NB, GBM, UVM, KIRC, LIHC, KIPAN, LGG, and ACC, patients with low BOLA2B expression showed better OS rates.

Similar results were obtained in the disease-specific survival (DSS) analysis. The hazard ratios of BOLA2B were significant for GBMLGG (HR = 1.87, 95%CI: 1.43–2.44), LGG (HR = 1.63, 95%CI: 1.10–2.42), KIRP (HR = 2.36, 95%CI: 1.44–3.89), KIRAN (HR = 1.69, 95%CI: 1.38–2.08), KIRC (HR = 2.06, 95%CI: 1.64–2.58), THYM (HR = 3.88, 95%CI: 1.09–13.75), LIHC (HR = 1.26, 95%CI: 1.03–1.55), UVM (HR = 5.39, 95%CI: 2.42–11.99), and OV (HR = 0.85, 95%CI: 0.78–0.93), as shown in Supplementary Figure 1A. Kaplan–Meier plots showed the specific survival probability for every significant cancer influenced by BOLA2B expression (Supplementary Figure 1B).

Finally, analysis of the progression-free interval (PFI) indicated that high levels of BOLA2B mRNAs were also strongly associated with worse prognosis in GBMLGG (HR = 1.81, 95%CI: 1.42–2.32), LGG (HR = 1.41, 95%CI: 1.05–1.90), KIRP (HR = 2.06, 95%CI: 1.21–3.53), KIPAN (HR = 2.17, 95%CI: 1.49–3.16), KIRC (HR = 3.36, 95%CI: 2.12–5.32), LIHC (HR = 1.67, 95%CI: 1.22–2.30), SKCM-M (HR = 1.32, 95%CI: 1.02–1.71), UVM (HR = 4.55, 95%CI: 1.99–10.42), KICH (HR = 3.99, 95%CI: 1.17–13.64), UCS (HR = 3.04, 95%CI: 1.45–6.36), and BLCA (HR = 0.56, 95%CI: 0.38–0.82) (Supplementary Figure 2A). Kaplan–Meier curves confirmed these results (Supplementary Figure 2B). These findings demonstrated the risks associated with abnormally increased BOLA2B expression and suggest that BOLA2B could be used as a potential biomarker to predict patient survival in multiple cancers.

mRNA modifications such as m1A, m5C, and m6A, are known to regulate protein translation. We analyzed the correlations between BOLA2B mRNA expression and the expression of related genes. As shown in Figure 3, correlations between genes involved in m1A, m5C, and m6A modifications and BOLA2B mRNA expression were present in numerous cancers.

BOLA2B m1A modifications (Figure 3) in DLBC, THYM, WT, KIPAN, PRAD, UCEC, THCA, CESC, KIRP, PCPG, BRCA, HNSC, READ, CHOL, MESO, and SARC were found to be negatively correlated with a reader gene, while positive correlations were seen with writer genes in cancers such as LUAD, LUSC, PAAD, COAD, COADREAD, TGCT, ALL, LAML, LIHC, GBMLGG, LGG, SKCM, NB, UCS, ESCA, STAD, STES, ACC, and BLCA. In terms of m5C modifications, the results showed that all the investigated cancers except OV were negatively associated with the eraser gene TET2. All cancers were highly positively associated with the writer gene NSUN5 and reader gene ALYREF synchronization. Regarding m6A, the writer gene KIAA1429 was positively correlated with BOLA2B across all cancers, while other m6A writer genes showed negative correlations with BOLA2B in most cancers. Of these cancers, OV was unique as it was positively associated with m1A, m6A, and m5C, regardless of their modifier regulators. These findings suggest that there are three different methods of post-transcriptional regulation of BOLA2B in different tumors.

Next, we analyzed interactions between BOLA2B expression and immune cell infiltration. As shown in Figure 4A, BOLA2B expression was negatively correlated with immune cell infiltration in most tumors. However, in LIHC, BOLA2B expression was positively associated with CD4⁺ T cells, neutrophils, macrophages, and DCs, all of which are involved in humoral immunity. Similar results were observed in KICH and SARC tumors. In UVM, BOLA2B expression showed a positive correlation with CD8⁺ T cells but was negatively correlated with other immune cells such as CD4⁺ T cells, neutrophils, macrophages, and DCs, suggesting that BOLA2B may affect cellular immunity in UVM.

Analysis of the association between BOLA2B and immune-related proteins and cytokines indicated that BOLA2B was negatively correlated with several immune-associated proteins such as EDNRB, ENTPD1, and TLR4 in most tumors, as shown in Figure 4B. We also found that OV, UVM, and SARC had specific differences in terms of immunity; also, in OV, BOLA2B was highly positively correlated with the regulation of immune cytokines.

In addition, we evaluated the effects of BOLA2B on the immune response in 33 cancers using the ESTIMATE algorithm. The ESTIMATE score consists of the purity score within the tumor and the immune score of the tumor. As shown in Figure 4C, the ESTIMATE score of UVM tumor was strongly associated with BOLA2B expression, in contrast to other tumors, such as DLBC, READ, and LUSC, where negative correlations were observed. Further analysis revealed that of all cancers, UVM showed the highest correlation with BOLA2B expression, in both the immune and stromal scores (Supplementary Figures 3A, B). DLBC, on the other hand, showed the most negative correlation with BOLA2B in the ESTIMATE score, which was largely due to the immune score (Supplementary Figures 3A, B). For both READ and LUSC tumors, BOLA2B expression influenced both immune and stromal aspects (Supplementary Figures 3A, B).

It is well known that tumor tissues consist of both tumor and non-tumor cells, including immune cells, stromal cells, and mesenchymal cells. Tumor purity is significantly associated with the clinical characteristics, genomic expression, and biological properties of cancers. We analyzed the effects of BOLA2B expression on tumor purity. As shown in the lollipop plot in Figure 4D, the tumor purity of LUAD, GBM, and SKCM cancers was strongly associated with BOLA2B. In addition, the interactions between the TMB and BOLA2B expression were investigated. As seen in Figure 4E, the lollipop plot shows that BOLA2B expression in DLBC, READ, and KIPAN was positively correlated with TMB but was negatively correlated in SARC, THYM, and UVM.

Mutations in genes often affect tumor occurrence and progression. Here, we explored several types of mutations associated with BOLA2B and their influence on various types of tumors. As shown in Figure 5A, the CNVs of BOLA2B revealed the presence of copy number amplifications of BOLA2B in several tumors, including CESC, LUAD, COAD, COADREAD, BRCA, ESCA, STES, SARC, STAD, HNSC, LUSC, LIHC, OV, and BLCA. These led to the increased expression of BOLA2B relative to WT and Del groups. The percentages of each type of BOLA2B CNVs in the 33 cancers were then plotted and illustrated in a pie chart. This showed that heterozygous copy number amplification was most prevalent in ACC, KIRP, and BRCA, while there were fewer CNVs in THCA, THYM, LGG, LAML, and PCPG (Figure 5B). We also found that survival prognosis, including OS, PFS, and DSS, was affected by CNVs in ACC, COAD, PCPG, and UCEC (Figure 5C). The relationship between single-nucleotide variants (SNVs) and BOLA2B expression was also examined in the different cancers. This showed that mutations in TP53, TTN, and MUC16 played key roles in influencing the expression of BOLA2B, and simultaneous mutations in TTN and MUC16 were presented in three different cancers (Figure 5D). This suggests that further exploration of the relationship between BOLA2B and mutations in these genes is warranted.

In addition, we explored the relationship between microsatellite instability (MSI) and BOLA2B. The occurrence of MSI in tumor tissues is due to a functional breakdown in DNA mismatch repair. The highest correlation between MSI and BOLA2B was observed in DLBC (Figure 5E). Finally, we depicted the relationship between tumor stemness and BOLA2B. As shown in Figure 5F, the scatter plot indicated that the stemness of TGCT was significantly positively correlated with BOLA2B, while that of LAML showed the most negative association with BOLA2B. The stemness of most other tumors was positively correlated with BOLA2B.

In this section, we analyzed the function of BOLA2B from a molecular biology perspective. The GO cellular component analysis showed that BOLA2B was associated with the TOR complex and histone acetylation (Figure 6A). In the biological process category, BOLA2B was associated with the aggregation of iron–sulfur clusters and cell-cycle arrest (Figure 6B), while in the molecular function category, BOLA2B was mainly involved in the formation of iron–sulfur clusters of both 2-valent and 4-valent iron and could function as a transcription factor (Figure 6C). KEGG pathway analysis showed that BOLA2B was mainly involved in the mTOR, autophagy, and PI3K-AKT pathways (Figure 6D).

To study the function of BOLA2B in cancer, the gene was knocked down in the MDA-MB-231 and SKBR3 breast cancer cell lines using shRNA. As shown in Figure 7A, shRNA1 and shRNA2 both downregulated the BolA2B protein level in MDA-MB-231 and SKBR3 cells. The CCK-8 assay showed that BOLA2B knockdown significantly reduced the proliferative capacity of breast cancer cells, as shown in Figure 7B. EdU assays were used to confirm the effects of BOLA2B on cell proliferation. As shown in Figure 7C, the percentage of proliferating cells in the BOLA2B-knockdown group was lower than that in the control group. These results demonstrated the involvement of BOLA2B in the proliferation of breast cancer cells.

Since knockdown of BOLA2B reduced cancer cell proliferation, we assumed that it would affect the cell cycle. Flow cytometry showed that BOLA2B knockdown in MDA-MB-231 and SKBR3 cells resulted in G2/M arrest (Figure 8A). Western blotting showed a significant upregulation of cyclin D1 and p21, compared with the control group. The levels of Cdc6 and cyclin E1 were slightly increased, while there was little influence on those of CDK6 and cyclin B1 (Figure 8B).

Previous studies have suggested that BOLA2B may be associated with the mTOR and PI3K-AKT pathways. We used western blotting to verify bioinformatics analysis. The results showed that in breast cancer, the signaling activation pattern may be dependent on the cell type. Both AKT and mTOR were upregulated in the BOLA2B-knockdown MDA-MB-231 cells, although there was a little change in their phosphorylation status. In SKBR3 cells, on the other hand, the mTOR pathway was significantly activated, while the MAPK pathway was inactivated. No change was seen in the AKT pathway (Figure 8C). We speculate that this phenomenon could be explained by different BRCA subtypes represented by MDA-MB-231 and SKBR3 cells.

Last, we performed the IHC staining of BolA2B in pan-cancer tissue samples. Representative images of THCA, LIHC, READ, and BRCA samples are shown in Figure 8D. Essentially, it was observed that BOLA2B was highly expressed in tumor tissues, except for gastric cancer and LIHC. This suggests that BOLA2B may have different functions in different cancer types.