Introduction

Breast cancer (Br Ca) is the most common cancer in women and the fifth leading cause of cancer death among females with cancer around the world1,2. Breast cancer accounts for one in eight diagnosed cancer cases out of an estimated 2.3 million new cancer cases, according to the latest WHO data3,4. In recent decades, the quality of life and survival of breast cancer patients have improved greatly due to the development of early screening, early diagnosis, early treatment, and the introduction of new drugs and technologies, but the incidence remains high5,6. Therefore, it is of great significance to explore the molecular mechanism of breast cancer occurrence and to search for related genes for targeted therapy.

Transcription factors (TFs) specifically recognize the DNA sequence (transcription factor binding sites) through their own DNA-binding domains to regulate gene expression7. Zinc finger protein 32 (ZNF32) is a transcription factor that belongs to the Krüppel family and has been linked to breast cancer progression8,9. At present, there are few reports on the role of this gene in tumors. ZNF32 protects cancer cells against oxidative stress-induced apoptosis via regulating C1QBP transcription, according to a previous study10. ZNF32 has also been shown to be able to regulate autophagy and protect breast cancer cells from the death triggered by stimuli11. Furthermore, Zfp637 as the mouse homologue of the ZNF32 could markedly increase mTERT expression and telomerase activity and maintain telomere length12. As recently reported, ZNF32 is involved in multidrug resistance in lung adenocarcinoma13. However, there are still many unanswered questions regarding the biological function of this gene, especially its exact mechanism for the proliferation, adhesion and survival of breast cancer cells.

Molecular targeted therapy, as a new means of cancer treatment, interferes with specific molecules to prevent cancer growth, progression and metastasis. At the same time, it is becoming a hot spot in cancer treatment research because of its low toxicity and high efficiency14,15. RNA-sequencing (RNA-Seq) technology which is a high-throughput sequencing technology has developed rapidly in recent years. Currently, RNA-seq has been widely used in basic medical research, which provides a powerful tool for finding target genes for cancer therapy. In this study, RNA-sequencing was used for the first time to compare the difference between breast cancer cells with overexpression of ZNF32 and the control cells. This study may provide an experimental basis for finding new therapeutic targets for breast cancer.

Materials and methods

Cell culture

The human breast carcinoma cell lines ZR-75-30 (ER+; PR-; HER2+) were purchased from American Tissue Culture Collection (ATCC, Manassas, VA, USA). Maintained in RPMI 1640 medium containing 10% FBS (Gibco, USA) in humidified atmosphere at 37℃, 5% CO2, and regularly tested for mycoplasma to verify their negative status.

Lentiviral infection of BC cells

ZNF32 overexpression lentivirus and its negative control (LV5 EF-1a/GFP&Puro) were purchased from Genepharma (Shanghai, China). Procedures was done followed the manufacturer’s instruction16,17. Stable overexpression ZR-75-30 cell lines were selected with puromycin and fluorescent intensity was detected under an Olympus IX-73 fluorescence microscope (Olympus, Tokyo, Japan). When the fluorescence in breast cancer cells with overexpression of ZNF32 and control cells (NC) reached 100%, we considered that the cells were infected and screened successfully.

RNA extraction, quality control and sequencing

When the cells are cultured to 90% in a 10 cm culture dish, discard the consumed culture medium, wash twice with PBS, add TRlzol to lyse the cells in the culture dish, and then transfer to an EP tube. The total RNA was extracted according to the instruction manual of the TRlzol Reagent (Life technologies, California, USA). Then the concentration and purity of RNA were measured by NanoDrop 2000 (Thermo Fisher Scientific, Wilmington, DE). In addition, the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA) is used to evaluate RNA integrity. The sequencing library was successfully constructed by using Hieff NGS Ultima Dual-mode mRNA Library Prep Kit for Illumina (Yeasen Biotechnology (Shanghai) Co., Ltd.). After the quality of the library was evaluated by AMPure XP system, the library was sequenced on the Illumina NovaSeq platform to produce a double-ended read length of 150 bp. The experimental strategy is shown in Fig. 1. The sequencing data generated in this study have been deposited in the NCBI Sequence Read Archive (SRA) database under Bioproject number PRJNA1092341. And can be reviewed via at the following URL: https://www.ncbi.nlm.nih.gov/sra/?term=PRJNA1092341.

Fig. 1
figure 1

The workflow of data collection and technical analysis for this study. The solid arrow shows the flow of the experiment. Diagram detailing the experiment including the sample model establishment, data acquisition and data analysis.

Full size image

Data filtering and mapping to reference genomes

Before the data analysis, the high quality Reads is obtained by filtering and removing the Reads containing joints and low quality Reads. All subsequent studies used high-quality data. These clean reads were then mapped to the reference genome sequence. Soft Hisat2 techniques were utilized for mapping with the reference genome. Then use StringTie to assemble the aligned reads for subsequent analysis.

Differential expression quantification and analysis

The expression level of each gene was normalized by the reads per kilobase per million (RPKM)18. DESeq2 provide statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution19. The resulting P values were adjusted using the Benjamini and Hochberg’s approach for controlling the false discovery rate. To identify DEGs, the DEseq2 package was used to filter the DEGs for the NC and OE groups. After the statistical analysis, we screened the DEGs with a fold change ≥ 2 and the false discovery rate (FDR) threshold was FDR < 0.01.

Gene ontology and pathway analysis

Gene Ontology (GO) enrichment analysis of the differentially expressed genes (DEGs) was implemented by the cluster profiler packages based Wallenius non-central hyper-geometric distribution20, which can intuitively understand the functional items related to differential genes. Moreover, the differentially expressed genes were annotated in the KEGG database21,22, and the annotation results were classified according to the pathway types in KEGG to draw a KEGG classification diagram. In addition, a hypergeometric test was applied to find pathways that were significantly enriched in differentially expressed genes compared with the entire genome background.

RT-qPCR validation

To validate the DEGs discovered by transcriptome sequencing, reverse transcription quantitative real-time PCR (RT-qPCR) was carried out. Primer Premier 5 software was used to design sequence-specific primers for randomly chosen genes (Table 1). Afterwards, real-time RT-qPCR was performed with qPCR SYBR Green SuperMix according to its manufacturer’s instructions (bimake, USA). 18 S rRNA gene as the endogenous reference gene to normalize the relative mRNA expression.

Table 1 Primers for RT-qPCR in this study.
Full size table

Western blot

Proteins from cells were extracted using RIPA Lysis Buffer (Biosharp, #BL504A, Beijing, China) supplemented with 1% protease inhibitor (Biosharp, #BL612A, Beijing, China). 25 μg proteins were loaded and separated on 12% SDS-PAGE and then transferred electrophoretically to a 0.45 μm polyvinyl difluoride membranes. β-actin antibody was ordered from Sangon (Shanghai, #D110001 China). Anti-ZNF32 antibody was ordered from Sigma (St. Louis, # SAB1406611 USA). All above antibodies were used in a dilution ratio of 1:1000. Horseradish peroxidase-conjugated secondary antibody to rabbit IgG (1:5000, Santa), and horseradish peroxidase-conjugated secondary antibody to mouse IgG (1:8000, Santa). The membrane was developed using Immobilon™ Western Chemiluminescent HRP Substrate (Millipore).

Statistical analysis

Statistical analysis was performed with GraphPad Prism 8.8 (GrpahPad Software, San Diego, CA, USA). To evaluate the significant differences between two groups, the means were compared using Student’s t-test. P values < 0.05 were considered statistically significant. *P < 0.05, **P < 0.01, ***P < 0.001 and ****P < 0.0001.

Results

Construction of breast cancer cell line stably overexpressing ZNF32

In order to study the role of ZNF32 in breast cancer development, we first generated ZR-75-30 breast cancer cells stably overexpressing ZNF32 (named as OE) by lentivirus. Because of the presence of fluorescent tags in the vector, when the fluorescence of cells in both the overexpression and the control reached 100%, we considered that the cells overexpressing ZNF32 were successfully constructed (Fig. 2A). WB and RT-qPCR were used to detect the overexpression efficiency, the results showed that the expression of ZNF32 in overexpression cells was significantly higher than that in the control cells (Fig. 2B C). In summary, the stably overexpressing ZNF32 breast cancer cell line ZR-75-30 was successfully generated.

Fig. 2
figure 2

ZNF32 overexpression breast cancer cells construction. (A) Representative images of ZR-75-30 after infected by lentivirus. (B) ZNF32 expression levels were detected by qPCR in ZNF32 over-expressing ZR-75-30 cells. (C) ZNF32 expression levels were detected by WB in ZNF32 over-expressing ZR-75-30 cells. The blots were cut based on the standard band positions, then incubated with the appropriate antibodies and in Additional File 1 the full-length blots are presented.

Full size image

RNA-seq of the ZNF32 overexpression and control cells transcriptome

To better understand the molecular changes in response to ZNF32 overexpression in BC cells, three cDNA libraries from control cells (NC-1, NC-2, and NC-3) and overexpression cells (OE-1, OE-2, and OE-3) were constructed. After sequencing quality control, a total of 44.21Gb Clean Data was obtained, and the percentage of Q30 bases in each sample was not less than 94.80%. The output statistics of sequencing sample data are shown in Table 2. In general, the sequencing results were highly reliable.

In addition, the correlation of gene expression level between samples is an important index to test the repeatability of biological experimental operations, evaluate the reliability of differentially expressed genes and assist the screening of abnormal samples. We used Pearson correlation coefficient r as the evaluation index of biological repetitive correlation. After calculation, we found that the correlation was high in the control group and overexpression group, indicating that the three biological repetitive samples in each group had strong correlation and good repeatability. The heat map of the correlation between samples in this study is shown in Fig. 3A. Principal component analysis (PCA) can reduce multiple variables to a few independent variables (principal components). In transcriptome analysis, the expression of genes will be reduced to find the law of sample distribution from complex data, so as to evaluate the degree of discretization of samples. The first principal component (PC1) can explain 40.09% of the characteristics of the original data set, and the second principal component (PC2) can explain 15.48% of the characteristics of the original data set (Fig. 3B). In this study, different treatments were significantly separated in the results of principal component analysis, which indicated that there were significant changes in the transcriptional level between the overexpression group and the control group.

Fig. 3
figure 3

Sample correlation analysis and PCA. (A) Sample correlation analysis, the value on each color block on the heat map in the figure represents the correlation between the two samples corresponding to the horizontal and vertical axes of the color block. The larger the value, the higher the correlation. (B) PCA, In PCA, different coordinates represent different principal components, percentage represents the contribution value of corresponding principal components to the difference of samples, each point represents a sample, the negative control group showed a red rectangle, and the overexpression group showed a blue circle.

Full size image
Table 2 Data output statistics of sequencing samples.
Full size table

Mining of DEGs

Differences in transcription levels between the overexpressed and control groups could be analyzed by clustering heat maps. The results of cluster heat map analysis showed that there were significant differences between the control group and the overexpressed group, and 3 biological copies in each group were clustered together, indicating that the biological copies had good consistency (Fig. 4A). The DEGs (Genes with significantly different expression levels) of the experimental group and the control cells were screened with FDR < 0.01 and Fold Change ≥ 2 as the threshold. A total of 604 DEGs were screened, of which 303 genes were up-regulated and 301 genes were down-regulated (Fig. 4B). All DEGs were listed in Supplementary Table 1. Table 3 lists the top 20 up-regulated and down-regulated genes with the largest fold change, and we use ZNF32 as a positive control.

Fig. 4
figure 4

RNA-seq results of differential expression genes. (A) Gene cluster expression heat map. (B) Volcano map to show the distribution of DEGs. The X axis represents Log2 fold change, and the Y axis represents -log10 FDR. Black is the genes of nonsignificant differences, red and blue are the genes of significant differences. The red color represents up-regulated genes, the blue color down-regulated genes and the black color represents unchanged genes.

Full size image
Table 3 The top 20 of up-regulated and down-regulated gene.
Full size table

ZNF32 is closely related to the development of breast cancer

In order to intuitively understand the main biological functions of differentially expressed genes, we classify differentially expressed genes at the secondary classification level of GO database. The results of GO classification (Fig. 5A) annotated 20 biological processes, 3 cellular components and 11 molecular functions. The main biological processes involved in differential genes are cellular processes, biological regulation, stress response, metabolic processes, multicellular biological processes, cellular signaling, development, localization, immune system processes, biological adhesion, development and reproduction, etc. Among the cellular components, the number of differential genes expressing cellular anatomical entities, intracellular and protein complexes is the largest. Among the molecular functions, the number of genes involved in binding, catalytic activity and transcriptional regulation is the largest. The results of gene ontology classification of genes related to the corresponding GO terms have been attached in supplementary Table 2. In addition, in order to understand the significant enrichment of GO entries compared with the whole genome background, we found that cell communication, positive regulation of cell population proliferation, positive regulation of receptor signal pathway through STAT, phosphorylation, signal transduction, regulation of cell migration and BMP signal pathway were significantly enriched in the biological process (Fig. 5B). Among the cellular components, the most abundant items are extracellular regions, extracellular spaces and presynaptic active areas (Fig. 5C). In molecular function, the top three are insulin-like growth factor binding, growth factor activity and calcium binding (Fig. 5D). These GO enrichment categories provide clues about the role of ZNF32 in breast cancer cells, which will help to screen target genes that contribute to these categories.

Fig. 5
figure 5

GO classification of DEGs. (A) The result of GO classification according to the number of genes enriched in each item. The abscissa is the GO classification, the left side of the ordinate is the percentage of the number of genes, and the right side is the number of genes. Biological Process (B), Cellular Component (C) and Molecular Function (D) were annotated according to the degree of differences.

Full size image

The results of KEGG pathway classification (Fig. 6A) showed that the signaling pathway with the most enriched differentially expressed genes in the cellular process classification was the focal adhesion, with a total of 6 up-regulated DEGs and 14 down-regulated DEGs enriched. The top 3 in the environmental information processing classification are PI3K-Akt signaling pathway, cytokine-cytokine receptor interaction and MAPK signaling pathway. Pathway in cancer rank first in the classification of human diseases. The corresponding KEGG enrichment results are attached in supplemental Table 3. The statistical results of pathway enrichment showed that the most significantly enriched pathway was focal adhesion, including ECM-receptor interaction, PI3K-AKT, HIPPO and TNF signaling pathways, which are closely related to cancer development (Fig. 6B). They participated in a variety of physiological process such as cell growth, development, proliferation, inflammatory response, occurrence and differentiation of tumor cells23,24,25. In the annotation results of focal adhesion pathway (Fig. 6C), the relationship between labeled proteins and corresponding DEGs is shown in Table 4. ECM protein is affected by the up-regulated gene FN1 and down-regulated gene LAMA1, COL1A1, COL4A5, LAMA5, NewGene1547 and SPP1, and then the changes of ITGA protein and ITGB protein are affected by up-regulated gene ITGA9 and up-regulated gene TTGB3, respectively. GF is affected by the down-regulated gene EGF, VEGFC, PDGFC, PDGFB, which then affects the PTK protein by down-regulating the gene KDR. Other proteins affected by DEGs include pl30Cas, Parvin, MLC, MLCK, Actin and PAK. These changes ultimately affect cell motility, proliferation and survival. In cancer pathway annotation results (Supplementary Fig. 1), DEGs involved Wnt signaling pathway, cAMP signaling pathway, ECM-receptor interaction, Focal adhesion, P13K-Akt signaling pathway, MAPK signaling pathway, PPAR signaling pathway, P53 signaling pathway, TGF-β signaling pathway, VEGF signaling pathway, HIF-1 signaling pathway and other signaling pathways. Understanding the differential genes in the above pathways is crucial for studying the biological function of ZNF32 in breast cancer. Taken together, ZNF32 may regulate the occurrence and development of tumor through the above signal pathways.

Fig. 6
figure 6

KEGG pathway classification and functional enrichment for DEGs. (A) KEGG Pathway annotation. (B) KEGG Enrichment analysis. Abscissa is the proportion of interested genes annotated in this entry to the number of differentially expressed genes, ordinate each KEGG annotation entry. The size of the dot represents the number of differentially expressed genes annotated in the pathway, and the color of the dot represents the Q value of the hypergeometric test. (C) Focal adhesion (ko04510) pathway annotation result map. Compared to the control cells, the proteins labeled in the red box were associated with up-regulated genes, and the proteins labeled in the green box were associated with down-regulated genes. The proteins labeled in blue are associated with both up-regulated and down-regulated genes. The KEGG pathway map is copyrighted by Kanehisa laboratories, license number 241,008.

Full size image
Table 4 Annotation proteins and their corresponding DEGs in focal adhesion.
Full size table

qRT-PCR validation of the RNA-seq results

The consistency and repeatability of DEGs discovered by transcriptome sequencing was validated using qRT-PCR. We selected the DEGs in the positive regulation of cell population proliferation and cell migration regulation in GO term (Fig. 5B) and the most significantly enriched DEGs in focal adhesion in KEGG (Fig. 6B) to draw heat maps respectively (Fig. 7). We first selected the up-regulated and down-regulated genes with the largest fold change among all differentially expressed genes as CA9 and LYPD6B, respectively. Secondly, we selected the up-regulated genes with the largest fold change in the focal adhesion pathway: PAK6, FN1, ITGA9 and the down-regulated gene KDR (Fig. 7A). Among the DEGs related to cell population proliferation, the two most up-regulated genes CRLF1, TFF2 and the most down-regulated gene IL15 were selected (Fig. 7B). Among the DEGs related to cell migration, the three most up-regulated genes were selected: TFF2, SNAI2, ENPP2 (Fig. 7C). All these genes were significantly differentially expressed and were consistently upregulated or downregulated with the gene expression changes based on RNA-Seq (Fig. 8A B). Table 5 showed that the fold change of the selected DEGs. All above suggesting that the DEGs obtained from transcriptome sequencing were reliable.

Fig. 7
figure 7

Clustering heat map of different functionally DEGs. (A) Focal adhesion (ko04510). (B) Positive regulation of cell population proliferation (GO :0008248). (C) Regulation cell migration (GO :0030334).

Full size image
Fig. 8
figure 8

Verification of the reliability of sequencing results. (A) qRT-PCR analysis showed the relative expression of CA9, LYPD6B, PAK6, FN1, ITGA9, KDR, CRLF1, TFF2, IL15, SNAI2 and ENPP2 in the overexpression cells and the control cells. (B) Comparison between differentially expressed genes and qRT-PCR results. X-axis shows 11 selected genes, and Y-axis represents relative multiple change. The data of this study is expressed as the mean standard deviation of three parallel measurements.

Full size image
Table 5 The Fold change of the selected DEGs.
Full size table

ZNF32 target genes prediction in breast cancer

ZNF32 belongs to the zinc finger protein super family and performs its biological functions as a transcription factor. ZNF32 binding sequence G (A/C/T) ATTT have been reported before26. CA9, LYPD6B, PAK6, FN1, ITGA9, KDR, CRLF1, TFF2, IL15, SNAI2 and ENPP2 were closely related to the proliferation, metastasis and invasion of tumor cells, meanwhile they were the differential genes after overexpression of ZNF32 in breast cancer cells. Then, we obtained the promoter sequences of above genes from the NCBI database and searched the binding sequence of ZNF32 in the promoter region these genes (-2000 to + 100). The target gene FASTA file containing the promoter sequence in the Additional file 2. The potential binding sequence of ZNF32 were found in CA9, PAK6, ENPP2, FN1, ITGA6, SNAI2, IL15, KDR and LYPD6B promoter area (Fig. 9). These results suggest that ZNF32 may regulate the biological processes such as cell proliferation, metastasis and invasion in breast cancer cells through transcriptional regulation of CA9, PAK6, ENPP2, FN1, ITGA9, SNAI2, IL15, KDR, LYPD6B (Fig. 10).

Fig. 9
figure 9

Prediction of binding sequences of ZNF32 in the target genes. (A) ZNF32 binding sequences predicted in CA9 promoter. (B) ZNF32 binding sequences predicted in PAK6 promoter. (C) ZNF32 binding sequences predicted in FN1 promoter. (D) ZNF32 binding sequences predicted in ENPP2 promoter. (E) ZNF32 binding sequences predicted in SNAI2 promoter. (F) ZNF32 binding sequences predicted in LYPD6B promoter. (G) ZNF32 binding sequences predicted in ITGA9 promoter. (H) ZNF32 binding sequences predicted in KDR promoter. (I) ZNF32 binding sequences predicted in IL15 promoter.

Full size image
Fig. 10
figure 10

Effects of ZNF32 and target genes on breast cancer tumor progression.

Full size image

Discussion

The reason why cancer is terrible is that tumor cells have biological characteristics such as dysregulation proliferation and differentiation, invasion and metastasis27,28,29. Breast cancer is initially a localized disease, but it can metastasize to lymph nodes and distant organs30. According to statistics, most breast cancer deaths are not due to the primary tumor itself, but from metastasis31. Not only that, drug resistance mediated by cancer cell metastasis limits the success of cancer treatment and is a huge obstacle that needs to be overcome clinically32. Due to the limitations of traditional methods such as radiotherapy and chemotherapy in breast cancer treatment, molecular targeted therapy has gradually become an indispensable part of modern breast cancer treatment. Therefore, the search for new molecular therapeutic targets for breast cancer is now a hot research topic. In recent years, scholars have found that ZNF32 plays an important regulatory role in the occurrence and development of breast cancer33. In this study, RNA-seq technique was used for the first time to compare the breast cancer cells with overexpression of ZNF32 and the control group, so as to determine the global expression changes related to the overexpression of zinc finger protein 32 in breast cancer cells.

In this study, we identified a total of 604 DEGs, including 303 up-regulated genes and 301 down-regulated genes. We observed that in breast cancer cells overexpressing ZNF32, some genes involved in key processes such as tumor proliferation, metastasis and adhesion were significantly up- or down-regulated. The up-regulated gene with the largest fold change was CA9. Carbonic anhydrase (CAs) are a superfamily of metalloenzymes found in all living kingdoms, balancing the reaction between CO2, bicarbonate and H. CA9, a membrane-associated α-CA, has been a drug target for various cancers34. Recent studies have found that the upregulation of CA9 expression can promote the progression of colorectal cancer35. It can change the acid-base homeostasis of tumor cells, resulting in pH gradient and driving cancer cell migration36,37,38,39. In addition, CA9 can also promote tumor environmental acidification, resulting in the acquisition of metastatic phenotype and chemical resistance of anticancer drugs. After overexpression of ZNF32, LYPD6B was the gene with the largest fold change among the down-regulated DEGs. Studies have shown that the invasiveness of melanoma cells is negatively correlated with the expression of hypermethylation gene LYPD6B40, so we speculate that the significant down-regulation of LYPD6B may increase the invasiveness of cancer cells.

GO enrichment analysis in the category of biological processes, we found that differentially expressed genes were significantly enriched in cell communication, positive regulation of cell population proliferation, signal transduction, and regulation of cell migration. In cell population proliferation, the two up-regulated genes with the largest fold change were CRLF1 and TFF2. In addition, the down-regulated gene with the largest fold change was IL15. Cytokine receptor-like factor 1 (CRLF1) is a secretory protein that promotes the malignant phenotype of thyroid carcinoma by activating MAPK/ERK and PI3K/AKT pathways41. Previous studies have shown that miR-3065-3p promotes dryness and metastasis of colorectal cancer by targeting CRLF142. Trefoil factor (TFF) family is a small secretory peptide characterized by trefoil structure composed of TFF1, TFF2 and TFF3. It has been reported that the expression of TFFs is up-regulated in a variety of malignant tumors, including gastric cancer43, colon cancer, hepatocellular carcinoma (HCC), breast cancer, prostate cancer, etc. In addition, TFF levels in these tumor tissues or patient serum are closely related to malignancy and prognosis. TFF promotes tumorigenesis and metastasis by enhancing cell proliferation, invasion, metastasis and angiogenesis and inhibiting apoptosis. Interleukin-15 (IL-15) is a pleiotropic cytokine with multiple effects that can improve the immune response to tumor cells to achieve the effect of cancer treatment. IL-15 has the potential to reject and destroy cancer cells in tumor microenvironment by expanding and activating immune-related cells44. It has been found that cytotoxic intrinsic lymphoid cells can sense IL15 expressed by cancer cells to inhibit malignant tumors in humans and mice45. Many existing studies have confirmed that the up-regulated expression of IL15 can inhibit the proliferation of cancer cells and help improve the prognosis of breast cancer46,47,48. In cell migration, the top three differentially expressed genes TFF2, SNAI2 and ENPP2 had the largest up-regulation fold change. siRNA-mediated silencing of SNAI2 inhibits the activation of Snail/Slug, thereby inhibiting gastric cancer cell proliferation, invasion and migration, EMT, tumor growth and lymph node metastasis49. The protein product of ectonucleotide pyrophosphates phosphodiesterase 2 (ENPP2) is Autotaxin (ATX), and ATX is related to the pathogenesis of various cancers, especially breast cancer50. Aberrant expression of ATX affects tumor progression, metastatic potential, and invasiveness51. In fact, ATX has been recognized as a potential diagnostic biomarker and drug target for cancer and chronic inflammatory diseases, ENPP2 is considered one of the top 40 genes that promote the metastatic process52,53. In addition, studies have shown that the tumor suppressor gene microRNA-101-3p can bind to the 3′UTR of ENPP2 mRNA and inhibit its translation, thereby reducing migration, invasion, and proliferation54.

Cancer cells undergo a series of processes such as isolation, migration, invasion, extravasation and proliferation in a certain part of the body, which will eventually lead to tumor metastasis. Transcriptome analysis showed that high expression of ZNF32 is accompanied by changes in downstream focal adhesion, ECM-receptor interaction. The extracellular matrix (ECM) supports multiple functions and plays an important role in regulating the tumor microenvironment. When normal cells become cancerous, if the ECM does not support these tumor cells, further development, invasion, and metastasis cannot occur. In the focal adhesion pathway, ECM proteins are affected by upregulation of FN1. Fibronectin 1 (FN1) drives the initiation and progression of multiple cancers through metabolic reprogramming. Recent studies have found that silencing FN1 inhibits breast cancer cell proliferation, invasion and migration, and promotes apoptosis55. Since triple-negative breast cancer (TNBC) is characterized by metastasis and invasion and has a poor prognosis, some researchers screened for differentially expressed genes between TNBC and normal cancer-free tissues, including fibronectin 1 (FN1), and through the protein interaction predicts it to be a hub gene in cell adhesion56. Integrins are a type of transmembrane glycoprotein composed of α subunit (ITGA) and β subunit (ITGB) that serve as cell surface adhesion molecules. The up-regulated expression of ITGA9 and ITGB3 affects the levels of integrin proteins. When extracellular matrix proteins serve as ligands and bind to receptors such as integrins, they subsequently recruit focal adhesion proteins and cytoskeletal elements to form stable adhesion. Research shows that ITGA9 is related to cell migration, cell invasion and epithelial-mesenchymal transition (EMT) of tumor cells57. ITGA9 levels in TNBC were significantly higher than those in other breast cancer subtypes and higher ITGA9 levels were associated with lower survival rates in patients with TNBC. Also significantly up-regulated is PAK6, which is a member of the p21-activated kinase (PAK) family of serine/threonine kinases and is closely related to the poor prognosis of various cancers such as liver cancer, prostate cancer, cervical cancer, and colorectal cancer58,59,60,61. Some studies have found that overexpression of PAK6 makes cells more resistant to chemotherapy drugs, while cells downregulated by PAK6 show sensitivity to drugs62. Compared with normal tissues, the mRNA expression of PAK6 in breast cancer is significantly up-regulated, and it is important in the regulation of cancer cell adhesion63,64. Cancer stem-like cells (CSCs) are a subset of tumor cells, which have remarkable abilities of self-renewal, proliferation, adhesion and metastasis65,66,67. Their existence is considered to be the key factor leading to tumor progression, metastasis and cancer drug resistance68,69. Previous studies have found that ZNF32 imparts stem-like properties to breast cancer by participating in signal transduction33. According to the results of this study, ZNF32 are closely related to focal adhesion, PI3K-Akt signaling pathway, Hippo signaling pathway and so on. Coincidentally, these pathways are also associated with tumor stem-like properties and cancer progression, and these findings provide strong data support for previous studies.

To sum up, we speculate that ZNF32, as a transcription factor, can transcriptional activation of CA9, PAK6, CRLF1, ENPP2, FN1, ITGA9, SNAI2 and TFF2, while inhibiting the transcriptional activities of IL15, KDR and LYPD6B. Through sequence search, we found potential binding sequences of ZNF32 in the promoter regions of CA9, PAK6, ENPP2, FN1, ITGA9, SNAI2, IL15, KDR and LYPD6B. Therefore, we speculate that transcription factor ZNF32 can bind to the promoter regions of these target genes, thereby promoting the proliferation, adhesion and metastasis of breast cancer cells (Fig. 10). As for the specific way to regulate it, further experimental exploration is needed. In addition to regulating these target genes which are closely related to the proliferation, migration and adhesion of cancer cells at the transcriptional level, ZNF32 can alter cancer-related signaling pathways in cells, such as PI3K-AKT signaling pathway, Hippo signaling pathway, TNF signaling pathway, etc. The occurrence and development of breast cancer is driven by dysregulation of cell signaling pathways, and ZNF32 may regulate the occurrence and development of tumors through these signaling pathways.

Conclusions

In this study, we constructed the stable overexpression of ZNF32 in breast cancer cell lines, and detected the transcriptomic changes associated with ZNF32 overexpression in breast cancer cells by RNA-Seq for the first time. Transcriptome analysis found a total of 604 differentially expressed genes, including 301 down-regulated genes and 303 up-regulated genes. Through KEGG enrichment analysis, we found that they are closely related to focal adhesion, ECM-receptor interaction, PI3K-Akt signaling pathway, Hippo signaling pathway and TNF signaling pathway. We screened out 11 differentially expressed genes with biological functions such as cell proliferation, migration, and invasion. The qRT-PCR results showed that the expression trends of the selected genes were consistent with the sequencing results. Sequence analysis revealed that ZNF32 transcriptional binding sequences exist in the promoter regions of CA9, PAK6, ENPP2, FN1, ITGA9, SNAI2, IL15, KDR and LYPD6B. Therefore, we speculate that transcription factor ZNF32 can transcriptionally regulate downstream target genes, thereby affecting the proliferation, metastasis, and invasion of breast cancer cells. This study deepens the understanding of ZNF32 function in breast cancer cells and provides a basis for finding new therapeutic targets for breast cancer.