Identification of differentially expressed genes and biological pathways in bladder cancer.

The purpose of the present study was to identify key genes and investigate the related molecular mechanisms of bladder cancer (BC) progression. From the Gene Expression Omnibus database, the gene expression dataset GSE7476 was downloaded, which contained 43 BC samples and 12 normal bladder tissues. GSE7476 was analyzed to screen the differentially expressed genes (DEGs). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed for the DEGs using the DAVID database, and a protein‑protein interaction (PPI) network was then constructed using Cytoscape software. The results of the GO analysis showed that the upregulated DEGs were significantly enriched in cell division, nucleoplasm and protein binding, while the downregulated DEGs were significantly enriched in 'extracellular matrix organization', 'proteinaceous extracellular matrix' and 'heparin binding'. The results of the KEGG pathway analysis showed that the upregulated DEGs were significantly enriched in the 'cell cycle', whereas the downregulated DEGs were significantly enriched in 'complement and coagulation cascades'. JUN, cyclin‑dependent kinase 1, FOS, PCNA, TOP2A, CCND1 and CDH1 were found to be hub genes in the PPI network. Sub‑networks revealed that these gene were enriched in significant pathways, including the 'cell cycle' signaling pathway and 'PI3K‑Akt signaling pathway'. In summary, the present study identified DEGs and key target genes in the progression of BC, providing potential molecular targets and diagnostic biomarkers for the treatment of BC.

Introduction

Bladder cancer (BC) is the fourth most common cancer in men and the seventh most common solid tumor in women worldwide, with an estimated 430,000 new cases diagnosed in 2012 (1,2). While the incidence rate is stable or declining in men, it exhibits an increasing trend in women (3). BC has a complex biological behavior, with frequent relapse and metastasis (4). Previous data shows that about one-third of initial BC cases will exhibit local progression and distant metastasis, and the 5-year survival rate is <62% (5). However, the mechanism underlying BC is not clear, and the mechanisms of occurrence, recurrence and metastasis are still unknown. Therefore, it is of great value to explore the molecular mechanisms involved in the apoptosis, proliferation, metastasis and invasion of BC for the improvement of prevention, diagnosis and therapy.

The histopathology and molecular pathways in BC pathogenesis have been described. Somatic copy number alterations in multiple regions have been identified in previous studies, including amplification of PPARG and E2F3, with loss of CDKN2A and RB1 (6,7). The Cancer Genome Atlas (TCGA) project reported that potential therapeutic targets had been identified in 69% of the bladder tumors investigated; 42% of the tumors were reported to have targets in the phosphatidylinositol-3-OH kinase/AKT/mTOR pathway, and 45% were reported to have targets in the RTK/MAPK pathway (8). So far, knowledge of the molecular biology of BC has lagged behind that of other cancers. No molecular or gene-targeting agents have been approved for the treatment of the disease. Therefore, understanding the molecular mechanism of BC is vital for the development of more precise diagnostic and effective therapeutic strategies.

With the continuous development of bioinformatics and molecular biology, it is possible to explore the mechanism of carcinogenesis and development at the molecular level. In previous decades, a large number of important signaling pathways in tumorigenesis were identified through analysis of the expression profiles of gene microarrays. This technology has also been used for genomic analysis, which may aid in the discovery of key genes that are interrelated with tumorigenesis (9).

In the present study, a gene expression profile (GSE7476) was downloaded from the Gene Expression Omnibus database (GEO). The differentially expressed genes (DEGs) between the controls and BC samples were analyzed. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses in the DAVID database were applied to analyze the functional enrichment and significant pathways associated with the DEGs. In addition, we constructed a PPI network to identify the critical DEGs and significant modules. This study aimed to investigate the involvement of genes critical to BC, and to promote the development of novel targeted agents for BC therapeutic intervention.

Materials and methods

Microarray data

Gene expression profiles for BC (GSE7476) were downloaded from the GEO on the NCBI website (http://www.ncbi.nlm.nih.gov/geo). The probe-level data were converted into the corresponding gene symbols to detect the expression of gene transcript levels, according to the annotation information downloaded from the platform GPL570 (Affymetrix Human Genome U133 Plus 2.0 Array), which contains 54,675 probes. The gene expression profiles consisted of 12 urothelial samples from patients with prostatic hyperplasia or renal failure with no evidence of bladder malignancy, and 43 tumor samples from different BC risk groups. The mean age of the BC risk groups, which consisted of 39 males and 4 females, was 77 years (10). In total, 15 low-grade superficial tumor samples, 13 high-grade superficial tumor samples and 15 high-grade muscle-invasive tumors samples were assigned to the BC risk group (10). The healthy control (HC) group, which comprised 12 males, had a median age of 59 years (10). The datasets from the 12 HC and 43 BC samples were analyzed.

Data processing and screening of DEGs

The CEL file data of GSE7476, downloaded from the GEO database, were read using the affy package in the R programming language (R). The original probe-level data were converted into gene symbols. Then, the expression values of multiple probes for the same gene were transformed into a single value by taking the mean expression value. The RMA method (robust multi-array average) was applied to carry out data pre-processing, including background correction, normalization and expression calculation of the original array data. The Limma package in R (11) was used to identify the DEGs between BC and HC samples. The Benjamini-Hochberg (BH) method (12) was introduced to adjust the raw P-values into a false discovery rate (FDR) to avoid the multi-test problem, which might produce too many false positive results. P<0.05 and |log2 fold change (FC)|≥1 were set as the thresholds for identifying DEGs.

Functional and pathway enrichment analysis of DEGs

GO and KEGG analyses were applied for the functional annotation and pathway analysis, using the Database for Annotation Visualization and Integrated Discovery (DAVID; https://david.ncifcrf.gov/) (13). The human genome was selected as the background parameter. P<0.05 and a count ≥2 were set as the thresholds to indicate a statistically significant difference.

PPI network construction and analysis of modules

PPI analyses may be helpful in identifying the generic organizational principles of functional networks, and to provide novel insights into protein function (14). In order to reveal the functional associations between proteins on a genome-wide scale, the STRING database (http://string-db.org/) online software (15,16) was used to construct a PPI network.

PPI networks were created after all DEGs were imported into the Cytoscape plugin. Confidence score ≥0.4 was set as the cut-off criterion. Molecular Complex Detection (MCODE) (17) was then applied to conduct module analysis in the resulting PPI network with the following parameters: Node score cutoff, ≥2; degree cutoff, ≥2; max depth, 100; and K-core, ≥2.

Results

Data preprocessing and DEG screening

The RNA was isolated from the tissue from BC and HC samples, respectively, for use in the microarray studies. A total of 20,487 gene symbols were discerned and the gene expression matrix of the samples was obtained. Based on the R analysis, a total of 1,173 DEGs were identified in BC compared with HC samples, including 859 upregulated genes and 314 downregulated genes. P<0.05 and |FC|≥2.0 were set as the threshold criteria. The top 10 upregulated DEGs and top 10 downregulated DEGs are listed in Table I. A volcano plot of the DEGs is presented in Fig. 1.

Table I.

Top ten upregulated and downregulated differentially expressed genes between bladder cancer and normal tissues.

A, The top 10 upregulated differentially expressed genes
Gene symbol	logFC	AveExpr	t	P-value	adj.P.Val	B
DPP3	1.367807192	6.823383737	−15.62834787	1.34E-09	6.41E-07	12.55923385
PAFAH1B3	2.465494694	7.611594612	−14.08412651	4.65E-09	1.26E-06	11.35210194
TFPT	1.064195524	6.593293201	−12.69386019	1.58E-08	2.55E-06	10.14074068
RANGAP1	1.069859311	6.672884639	−12.47693921	1.93E-08	2.81E-06	9.939786046
IGFBP3	2.382533743	10.63365369	−12.1139983	2.73E-08	3.49E-06	9.59578721
PVRL4	2.113260416	7.562648025	−11.09448808	7.52E-08	7.52E-06	8.574037157
SEC61A1	1.001492679	8.052241965	−10.79227308	1.03E-07	9.34E-06	8.254355505
MTFP1	1.022940775	6.246895347	−10.71273652	1.12E-07	9.91E-06	8.168871656
ESRP1	1.760663741	8.297155088	−10.59987707	1.27E-07	1.10E-05	8.046588009
ABRACL	1.327112508	7.101266084	−10.57221373	1.30E-07	1.12E-05	8.016436726
B, The top 10 downregulated differentially expressed genes
Gene symbol	logFC	AveExpr	t	P-value	adj.PVal	B
SCARA5	−3.36100537	4.768657761	30.97352259	3.13E-13	5.84E-09	19.86271169
LINC01082	−3.753683567	5.294466313	29.16744989	6.59E-13	5.84E-09	19.29373062
OLFML1	−2.95576203	4.434015195	28.55851363	8.55E-13	5.84E-09	19.08956481
TMEM100	−2.533343457	3.70073062	26.68801945	1.97E-12	1.01E-08	18.41959434
MIR100HG	−5.32118364	4.747251266	25.79518302	3.01E-12	1.23E-08	18.07515932
CFD	−4.881864287	8.124954765	24.85704639	4.74E-12	1.62E-08	17.69452213
SLIT2	−4.050964305	5.040374213	23.90645097	7.67E-12	2.09E-08	17.28787436
PRDM6	−3.036050363	4.838278315	23.58237834	9.07E-12	2.09E-08	17.14414695
MRGPRF	−4.045214224	5.668562294	23.55835407	9.18E-12	2.09E-08	17.13338551
LRFN5	−1.962030571	3.960439049	22.38800885	1.72E-11	3.51E-08	16.59060679

Figure 1.

Volcano plot of the DEGs. The abscissa represents logFC and the ordinate represents-log10 (P-value). The red and green dots indicate DEGs, while the black dots represent genes that are not differentially expressed between bladder cancer and healthy control tissues. Red, upregulation; green, downregulation. DEG, differentially expressed gene; FC, fold-change.

To further explore the systematic characterization and biological functions of the identified DEGs, functional annotation and pathway analysis, including GO and KEGG, were performed using DAVID.

In this study, the three GO categories [cellular component (CC), biological process (BP) and molecular function MF)] were detected, respectively, using DAVID. The top 15 GO terms of the upregulated and downregulated DEGs are shown in Table II and Fig. 2, respectively. The upregulated DEGs were involved in the different GO terms, such as ‘cell division’ (ontology: BP), ‘nucleoplasm’ (ontology: CC) and ‘protein binding’ (ontology: MF) (Table IIA and Fig. 2A). The most significantly downregulated DEGs were related to the GO terms ‘extracellular matrix organization’ (ontology: BP), ‘proteinaceous extracellular matrix’ (ontology: CC) and ‘heparin binding’ (ontology: MF) (Table IIB and Fig. 2B).

Table II.

The top 15 enrichedgene ontology terms of up-regulated DEGs and downregulated DEGs.

A, The top 15 enriched gene ontology terms of the upregulated DEGs
Category	Term	Count	P-value
BP	Cell division	47	2.13E-28
BP	Mitotic nuclear division	28	7.12E-15
BP	sSster chromatid cohesion	19	6.92E-14
BP	G1/S transition of mitotic cell cycle	17	9.83E-12
BP	Cell proliferation	28	7.54E-11
CC	Nucleoplasm	97	6.57E-14
CC	Condensed chromosome kinetochore	17	6.24E-13
CC	Midbody	19	3.04E-12
CC	Spindle pole	16	2.50E-10
CC	Kinetochore	14	5.65E-10
MF	Protein binding	198	6.15E-12
MF	Protein kinase binding	27	4.95E-10
MF	ATP binding	51	5.69E-07
MF	Microtubule binding	14	3.52E-05
MF	Microtubule motor activity	8	3.11E-04
B, The top 15 enriched gene ontology terms of the down-regulated DEGs
Category	Term	Count	P-value
BP	Extracellular matrix organization	45	1.92E-19
BP	Platelet degranulation	27	2.96E-13
BP	Cell adhesion	58	1.86E-12
BP	Muscle organ development	20	9.72E-09
BP	Muscle contraction	21	4.34E-08
CC	Proteinaceous extracellular matrix	60	9.31E-26
CC	Extracellular matrix	62	5.68E-25
CC	Extracellular space	132	2.99E-19
CC	Extracellular region	145	6.81E-18
CC	Sarcolemma	26	1.01E-14
MF	Heparin binding	37	1.33E-16
MF	Integrin binding	26	1.77E-12
MF	Extracellular matrix binding	11	5.86E-08
MF	Actin binding	34	1.34E-07
MF	Calcium ion binding	62	2.85E-07

DEG, differentially expressed gene; BP, biological process; CC, cellular component; MF, molecular function.

Figure 2.

The top 15 enriched Gene Ontology terms of upregulated (A) and downregulated (B) differentially expressed genes. DEG, differentially expressed gene.

Subsequently, KEGG pathway analysis demonstrated that the upregulated DEGs were enriched in five key pathways (Table IIIA and Fig. 3A), including ‘cell cycle’, ‘DNA replication’ and ‘p53 signaling pathway’, whereas the downregulated DEGs were enriched in five different pathways (Table IIIB and Fig. 3B), including ‘complement and coagulation cascades’, ‘focal adhesion’ and ‘hypertrophic cardiomyopathy (HCM)’.

Table III.

The top five enriched pathways of upregulated differentially expressed genes and downregulated differentially expressed genes.

A, The top 5 enriched pathways of the upregulated DEGs
Pathway ID	Name	Count	P-value	Genes
hsa04110	Cell cycle	21	9.58E-14	CDC7, CDC6, CDK1, SKP2, TTK, ESPL1, CDC20, PTTG1, MCM2, MCM4, MCM5, CCNB1, CCNE1, CCND1, MAD2L1, MCM7, CCNB2, PLK1, PCNA, BUB1B, CCNA2
hsa03030	DNA replication	10	1.47E-08	RFC5, PRIM1, RFC4, MCM7, PCNA, MCM2, RNASEH2A, MCM4, MCM5, FEN1
hsa04115	p53 signaling pathway	8	2.53E-04	BID, CCNB1, CCNE1, CDK1, CCND1, CCNB2, RRM2, IGFBP3
hsa05219	Bladder cancer	6	9.83E-04	CCND1, FGFR3, VEGFA, CDH1, MMP1, DAPK1
hsa04114	Oocyte meiosis	9	1.03E-03	CCNE1, CDK1, PPP1CA, MAD2L1, PLK1, AURKA, CDC20, ESPL1, PTTG1
B, The top 5 enriched pathways of the downregulated DEGs
Pathway ID	Name	Count	P-value	Genes
hsa04610	Complement and coagulation cascades	19	1.96E-09	C7, A2M, C5AR1, C3, F13A1, F8, C1R, SERPING1, BDKRB2, C1S, PLAUR, CD55, F3, SERPINE1, CFH, TFPI, CFI, CFD, PROS1
hsa04510	Focal adhesion	31	3.28E-08	CAV2, CAV1, PPP1R12B, TNC, MYL9, VCL, LAMB2, COL6A3, ILK, PPP1R12A, COL6A2, PDGFC, ZYX, PDGFD, THBS1, ITGA1, IGF1, ACTN1, FLNC, COL5A1, COL4A6, FLNA, LAMA2, VEGFC, ITGA5, FYN, JUN, ITGA8, ITGA7, MYLK, PARVA
hsa05410	Hypertrophic cardiomyopathy	15	1.55E-05	ACTC1, IL6, CACNA2D1, ITGA1, IGF1, TPM2, TPM1, DES, SGCG, ITGA5, DMD, ITGA8, ITGA7, SGCA, SGCB
hsa05205	Proteoglycans in cancer	25	2.53E-05	FGFR1, CAV2, CAV1, LUM, PPP1R12B, DCN, TIMP3, SDC2, ANK2, PPP1R12A, RRAS, FAS, THBS1, FGF2, IGF1, FLNC, FZD7, FLNA, ITPR1, WNT2B, PLAUR, FZD10, ITGA5, HSPB2, HBEGF
hsa05414	Dilated cardiomyopathy	15	3.70E-05	CACNA2D1, ACTC1, ITGA1, IGF1, TPM2, TPM1, DES, SGCG, ITGA5, DMD, PLN, ITGA8, ITGA7, SGCA, SGCB

DEG, differentially expressed gene.

Figure 3.

The top five enriched pathways of upregulated (A) and downregulated (B) differentially expressed genes. DEG, differentially expressed gene.

PPI network construction and module selection

STRING was applied to construct the PPI network of the DEGs. This PPI network consisted of 959 nodes interacting via 6,400 edges. Seven hub genes appeared in the top 10 genes list in terms of degree, betweenness and closeness, simultaneously. Among these genes, Jun proto-oncogene (JUN) showed the highest node degree, which was 144. The others included cyclin-dependent kinase 1 (CDK1, degree=125), Fos proto-oncogene (fos, degree=122), proliferating cell nuclear antigen (PCNA, degree=101), topoisomerase (DNA) II alpha (TOP2A, degree=100), cyclin D1 (CCND1, degree=98) and cadherin 1 (CDH1, degree=98).

Moreover, 27 functional clusters were selected from the PPI network using MCODE. The top 3 significant modules were selected (Fig. 4), and the pathway enrichment annotation of the genes involved in the modules was analyzed using KEGG pathway analysis, which revealed that the genes in modules 1–3 were mainly associated with the ‘cell cycle’ signaling pathway (both appearing in module 1 and module 2), and ‘PI3K-Akt signaling pathway’.

Figure 4.

Top 3 modules from the PPI network. (A) The top 3 significant modules were selected. (B) Module 1 of DEGs from the PPI network; (C) module 2 of DEGs from the PPI network; and (D) module 3 of DEGs from the PPI network. PPI, protein-protein interaction; DEG, differentially expressed gene.

Discussion

BC is one of the most common types of malignant cancer in China and has a high mortality rate (18). BC is the most common form of urinary tract malignant tumor. Approximately 95% of bladder tumors are urothelial, and their treatment mainly centers around surgery; however, relapse and metastasis after surgery are common (4). The key genes and pathways associated with BC were identified in the present study using bioinformatics methods.

In the present study, R was used to extract the genetic information from GSE7476, and a total of 1,173 genes were identified to be differentially expressed between BC and HC samples, among which 314 were upregulated and 859 were downregulated in BC. The upregulated DEGs were mainly enriched in ‘cell division’, ‘nucleoplasm’ and ‘protein binding’, while the downregulated DEGs were mainly involved in ‘extracellular matrix organization’, ‘proteinaceous extracellular matrix’ and ‘heparin binding’. Moreover, the KEGG pathway enrichment analysis results showed that the upregulated DEGs were associated with the ‘cell cycle’, ‘DNA replication’ and the ‘p53 signaling pathway’, whereas the downregulated DEGs were mainly enriched in the ‘complement and coagulation cascades’, ‘focal adhesion’ and ‘HCM’.

Previous studies have demonstrated that tumor development is associated with the activation of the coagulation cascade. The exact mechanism through which coagulation proteins promote tumorigenesis remains unclear; however, it is possibly associated with hemostatic factor changes and peritumoral deposition of fibrin (19–21). The cell cycle is the series of events that occur between cell duplication and division, and is closely associated with cell growth, anabolism and proliferation (22). Uncontrolled cell proliferation and DNA replication comprise one of the hallmarks of cancer (23). p53 is known to be mutated in >50% of all human cancers, including bladder carcinoma (24). Alterations in p53 expression levels are correlated with tumor recurrence, lower survival rates (25) and poor prognosis in BC patients (26). Therefore, investigating these signaling pathways may aid in elucidating the carcinogenic mechanism of BC.

In addition, a PPI network was constructed to identify the key DEGs. We used proteins that corresponded to genes to construct the PPI network, and found that seven hub genes (JUN, CDK1, FOS, PCNA, TOP2A, CCND1 and CDH1) appeared in each of the top 10 gene lists in terms of degree, betweenness and closeness. JUN and FOS both exhibited downregulated expression, and were identified as main hub genes, with degree values of 144 and 122, respectively. FOS and JUN are proto-oncogenes belonging to the family of activator protein 1 (AP1) transcription factors (27,28). Ye et al (29) reported that AP-1 plays a vital role in cellular migration, metastasis, proliferation, transformation, apoptosis and inflammation. C-FOS, a major member of the FOS family, has been demonstrated to be involved in the regulation of cell growth, differentiation, proliferation, transformation and apoptosis (30). Previous studies have shown that the level of C-FOS in BC tissues is significantly higher than that in adjacent non-cancer and normal tissues (31,32). Most of the research on c-Jun (a major member of the JUN family), indicates that it may contribute to tumor initiation and invasiveness (33,34). Huhe et al (35) revealed that high c-Jun expression served a vital role in tumor progression, and may be a diagnostic and therapeutic biomarker in urothelial carcinoma of the bladder.

The results of our study also showed that CCND1 and CDH1 were enriched in several pathways in BC. CCND1, a cell cycle regulatory factor, promotes the progression of the cell cycle through the G1/S phase limit points. Overexpression of the CCND1 gene can result in uncontrolled cell proliferation and tumor occurrence by shortening the G1 phase (36,37). Many patients with cancer have been found to overexpress CCND1, and thus to have a poor prognosis (38,39). Xu et al (40) reported that the expression of CCND1 is associated with the progression of BC; therefore, CCND1 may be considered as an auxiliary diagnostic factor and potential prognostic marker for BC patients. CDH1 encodes a classical cadherin of the cadherin superfamily. A previous report (41) showed that CDH1 plays an important role in suppressing the invasive phenotype of urothelial BC cells. Many studies have shown that the classical cadherins and related molecular pathways may be attractive therapeutic targets to restrain tumor progression in patients with BC (42–46).

Cell cycle progression is controlled by cyclin-dependent kinases (CDKs) and cyclins. Cell cycle dysregulation may lead to uncontrolled cell proliferation and the subsequent development of cancer (47). CDK1 regulates the G1-S transition in the cell cycle, a process that is important for the development of centrosome mutation (48). CDK1 is a vital regulator in cell proliferation, and overexpression may lead to high tumor aggressiveness and poor prognosis (49–51). Some antibodies, including anti-CDK1, have been used to investigate cell proliferation (52). One study revealed that determination of the specific activity of CDK1 may be useful in the prediction of outcomes in breast cancer patients (53). Therefore, CDK1 may also play an important role in BC tumorigenesis, and further study is required to identify whether it may serve as a potential molecular marker associated with BC.

PCNA, which encodes a nuclear protein that functions as a cofactor of DNA polymerase delta, serves as an important proliferative marker in carcinogenesis (54). The synthesis rate of PCNA has a direct impact on the proliferative rate of cells (55). An early study reported that significant clinical information obtained from immunohistochemical staining for PCNA may be helpful in the initial selection of therapies and the evaluation of chemotherapeutic effects (56). Malkas et al (57) reported that polyclonal antibodies against cancer-associated PCNA (caPCNA), which can serve as a diagnostic marker of breast cancer, have been developed. Therefore, further investigation is necessary to clarify the underlying biological links between PCNA and BC.

TOP2A is an essential nuclear enzyme involved in DNA replication, and its expression is decreased at the end of mitosis and increased during the S to G2/M phases in bladder urothelial carcinoma (58). Many studies have reported that the expression of TOP2A is increased in skin, breast, brain, ovarain and small cell lung cancers, and such increased expression is associated with shortened survival (59–63). Overexpression of TOP2A has been demonstrated to be related to recurrence and increased risk of death (64), and with late-stage BC (65). Lindén et al (66) reported that TOP2A could serve as a vital urinary biomarker candidate for BC. However, further investigation is required to elucidate the exact mechanism of action of TOP2A in the development and progression of BC.

The module analysis in the PPI network demonstrated that the development of BC was associated with the cell cycle signaling pathway and PI3K-Akt signaling pathway. It is well known that the cell cycle signaling pathway plays a key role in controlling the normal progression of the cell cycle. In the entire cell cycle regulatory network, abnormalities in various types of molecules can affect cell proliferation and apoptosis, potentially leading to uncontrolled cell growth and ultimately causing tumors. Akt plays a central role in the signaling pathways involved in cell growth, proliferation, angiogenesis, metabolism, apoptosis and migration (67,68), and has already been found to be associated with cancer (69). Many studies have identified PI3K/Akt overexpression and activation in a variety of tumor tissues, such as ovarian cancer, colorectal cancer, lymphoma, pancreatic cancer, non-small cell lung cancer, lymphoma and gastric cancer (70–75). Therefore, blocking the cell cycle and inhibiting the PI3K-AKT signaling pathway are promising approaches for therapeutic intervention in BC patients.

In conclusion, the current study aimed to identify DEGs involved in the progression of BC via comprehensive bioinformatics analysis. This study provides several key genes and pathways for future investigation into the mechanisms and biomarkers of BC. However, a lack of experimental verification is a limitation of this study. Further experimental research is necessary to investigate the pathogenic mechanism of BC.