Identifying liver metastasis-related hub genes in breast cancer and characterizing SPARCL1 as a potential prognostic biomarker

Discovery and validation datasets

In the GEO (http://www.ncbi.nlm.nih.gov/geo/), datasets GSE124648 (Sinn et al., 2019) and GSE58708 (McBryan et al., 2015) were obtained with the following keywords: breast neoplasms, breast cancer, liver metastasis, expression profiling by array, attribute name tissue, and Homo sapiens. Array data of liver metastases (N = 16) and primary tumors (control, N = 130) from the GPL96 (HG-U133A; Affymetrix Human Genome U133A Array) platform were selected for differential expression analysis. The array data for GSE58708 consisted of three patients with BC liver metastasis vs. three controls as a validation dataset (Table 1).

DEG identification

R software (version 4.0.3; R Core Team, 2020) was employed for bioinformatics analysis. The “limma” package was adopted to identify DEGs between BC and liver metastasis (Ritchie et al., 2015), with the following cut-off criteria: —log₂FC—>2.0 and adj.P.Val<0.05. Similarly, the R package “TCGAbiolinks” was used to obtain biological data for BC in the TCGA database (Colaprico et al., 2016), and hub gene expression was verified by analyzing the differentially expressed mRNAs in BC. “ggplot2” and “heatmap” packages were used for visualizing the DEGs.

Enrichment analysis of DEGs

The “clusterProfile” package was employed to conduct gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of the DEGs (Yu et al., 2012), and a gene set enrichment analysis (GSEA) method was used for the KEGG enrichment analysis. p < 0.05 was the criteria set for statistical significance.

Integration of PPI networks and identification of hub genes

The STRING (http://www.string-db.org/) database provides comprehensive information on protein–protein interactions (PPI) (Szklarczyk et al., 2021). A PPI network of DEGs was integrated using the STRING database. Then, the results were visualized and analyzed using Cytoscape (version 3.7.2). The important modules were extracted via the plug-in MCODE, with the cut-off criteria being an MCODE score>5 and nodes>5. The other default parameters were set as Max. Depth = 100, K-Core = 2, Node score cut-off = 0.2, Degree Cut-off = 2. The plug-in cytoHubba was employed to identify the top 30 hub genes in the PPI network, with the following cut-off criteria: degree.layout ≥25 (MCC algorithm) (Shannon et al., 2003). GEPIA (http://gepia.cancer-pku.cn/) was used to obtain the hub gene expression between BC and normal tissues, with the following default parameters: p ≤ 0.01 and —Log₂FC— ≥1 (Tang et al., 2017). The Kaplan–Meier plotter (https://kmplot.com/analysis/) was used to evaluate the correlation between the expression of the hub genes and patient survival (Lanczky & Gyorffy, 2021). The median expression was used to divide the patients into high and low expression groups. Then, the Kaplan–Meier curves of overall survival (OS) and recurrence-free survival (RFS) were drawn. Finally, hub genes (SPARCL1 and SERPINA1) with prognostic values were identified. p < 0.05 using the log-rank test was considered statistically significant.

Validation of hub genes and clinicopathological correlation analysis

A box plot was plotted to verify SPARCL1 expression in liver metastasis based on the GSE58708 dataset. The correlations between SPARCL1 expression and clinicopathological characteristics of patients with BC were analyzed using TCGA-BRCA clinical data. Then, based on patient survival data and SPARCL1 levels, the optimal expression threshold was calculated using the “survminer” package. The correlation between gene expression and clinicopathological parameters was tested using Pearson’s chi-squared test. p < 0.05 was considered statistically significant.

GSEA

GSEA can be used to evaluate whether a predefined gene set shows statistically significant differences between two groups with different phenotypes. Expression (.gct) and phenotype (.cls) information files for SPARCL1 were uploaded to the GSEA software (version 4.1.0) to conduct enrichment analysis. The chip platform and gene set database selected were “Human_ENSEMBL_Gene_ID_MSigDB.v7.0” and “c2.cp.kegg.v7.1. symbols.gmt,” respectively. The normalized enrichment score was calculated. Both false discovery rate (FDR q-val) and nominal p-value (NOM p-val) < 0.05 were considered statistically significant.

Clinical specimens and cell lines

Tissue samples were obtained from patients with BC who underwent radical mastectomy at Shanghai Tenth People’s Hospital (Shanghai, China). A total of 30 pairs of tumor and matched normal adjacent breast tissues were collected and preserved in liquid nitrogen. All clinical samples were obtained with the written informed consent of participants, and the project was approved by the Ethics Committee of Shanghai Tenth People’s Hospital (No. 2020-KN174-01). BC cells (MDA-MB-231, BT549, and MCF-7) and mammary epithelial cells (MCF-10A) were purchased from the Cell Bank of Type Culture Collection of Chinese Academy and grown in mammary epithelial basal medium (Cambrex, East Rutherford, NJ, USA) and Gibco Dulbecco’s modified Eagle medium (DMEM; Thermo Fisher Scientific, Waltham, MA, USA) supplemented with 10% fetal bovine serum (FBS; Gibco) and 1% penicillin-streptomycin (Enpromise, Shanghai, China). All cells were maintained in a 5% CO₂ incubator at 37 °C.

Cell transfection, RNA extraction, and RT-qPCR

Lentiviral vector pLKO.1 and lentiviral packaging kit (Genomeditech, Shanghai, China) was used for construct stable SPARCL1 knockdown cells. sh-NC and shRNAs targeting SPARCL1 (SPARCL1 sh-1 and SPARCL1 sh-2) were purchased from IBSbio (Shanghai, China) and the sequences were as follows: SPARCL1 sh-1: 5′- CCCACAATGATAACCAAGAAA-3′, SPARCL1 sh-2: 5′- GCAGAGAAATAAAGTCAAGAA-3′. The cells transduced by lentivirus were selected with 2 µg/ml puromycin for 3 days. The pcDNA3.1 vector was used to construct SPARCL1 overexpression plasmid (IBSbio, Shanghai, China). Plasmids were transfected into BC cells using Invitrogen Lipofectamine 3000 (Thermo Fisher Scientific).

A lentiviral short hairpin RNA (shRNA) construct allows for stable and long-term knockdown of the targeted gene. Invitrogen TRIzol reagent (Thermo Fisher Scientific) was used to extract total RNA from cell lines and tissues. HiScript III RT SuperMix for qPCR (Vazyme, Nanjing, China) was used to generate the cDNAs, and Hieff qPCR SYBR Green Master Mix (Yeasen, Shanghai, China) was used to conduct real-time quantitative-polymerase chain reaction (RT-qPCR) following the manufacturer’s protocol, with ACTB as internal control for normalizing SPARCL1 expression. We used the following primers to conduct RT-qPCR. SPARCL1 (forward: 5′-CCAACTGAAGGTACATTGGACAT-3′, reverse: 5′-CTGTGAAGGAACTAACACCAGG-3′) and ACTB (forward: 5′-CATGTACGTTGCTATCCAGGC-3′, reverse: 5′-CTCCTTAATGTCACGCACGAT-3′).

MTT, colony formation, wound healing, and Transwell assays

To evaluate the effects of SPARCL1 on BC cells, post-transfection BC cell lines were cultured in 96-well plates at a density of 1,500 cells. Next, 20 µL methylthiazolyldiphenyl-tetrazolium bromide (MTT; Yeasen) was added to each well at 0, 24, 48, 72, and 96 h after inoculation, and the samples were incubated for 4 h at 37 °C in a 5% CO₂ incubator to assess cell viability. Then, 150 µL DMSO was added to each well after removing the supernatant. The absorbance was measured using a microplate spectrophotometer (BioTek Instruments, Winooski, VT, USA) at 490 nm. Cell proliferation curves were plotted according to the absorbance value.

Transfected BC cell lines were prepared as single-cell suspensions, inoculated at a density of 750 cells/well until prominent colonies were formed. The cells were fixed with 95% ethanol and stained with 0.1% crystal violet (Yeasen) to detect cell colony formation ability. Representative pictures were recorded, and clone colonies were counted.

To detect cell mobility, the cell monolayer was scratched with RNA free tips when the transfected BC cell density reached 90% or more. In the following step, DMEM supplemented with 2% FBS was used as the culture medium. The healing of scratches was observed at 0 and 12 h using the same field of view to calculate cell mobility.

The transfected BC cells and 500 µL of DMEM containing 10% FBS were added to Transwell’s upper and lower chambers (Corning, Corning, NY, USA), respectively. After culturing for 18 h, migrated cells were fixed with 4% paraformaldehyde and stained with 0.1% crystal violet to assess the migratory ability. Representative images were captured using an inverted microscope.

Statistical analysis

Wilcoxon matched-pairs signed-rank test was used to compare the expression of SPARCL1 between BC and control samples. Unpaired Student’s t-test was used to compare the expression of SPARCL1 between MCF-10A and BC cell lines. The results of the MTT assay were analyzed using a two-way analysis of variance. All experiments were repeated three times. The experimental data were analyzed and plotted using Prism v8.3.0 (GraphPad Software, San Diego, CA, USA). p < 0.05 was considered statistically significant.

DEG identification

We identified 332 DEGs in the GSE124648 dataset, of which 116 and 216 were upregulated and downregulated genes, respectively (Fig. 1A). A heatmap was used to visualize the top 50 liver metastasis-related genes in the dataset (Fig. 1B). Among the top 50 DEGs, 60% showed significant down expression in liver metastasis compared to BC tissues. Another 40% of highly expressed DEGs showed higher expression compared to BC tissues. These results suggest that downregulated DEGs in BC liver metastasis probably play a more critical role than upregulated DEGs.

GO and KEGG enrichment analysis of DEGs

We further performed functional and pathway enrichment analyses on the 332 DEGs. GO annotation divides gene function into three categories: cellular components (CC), molecular function (MF), and biological process (BP) (Sinn et al., 2019). The top eight enriched GO terms for each category are shown in Fig. 2. CC analysis indicated that these DEGs were particularly related to the extracellular matrix, endoplasmic reticulum lumen, collagen-containing extracellular matrix, collagen trimer, and blood microparticles (Fig. 2A). MF analysis showed that the DEGs were mainly involved in extracellular matrix structural constituents, glycosaminoglycan binding, extracellular matrix structural constituents conferring tensile strength, heparin-binding, and collagen binding (Fig. 2B). The BP category was mainly enriched in extracellular structure organization, extracellular matrix organization, wound healing, humoral immune response, and complement activation (Fig. 2C). The KEGG pathway enrichment analysis results further showed that the upregulated DEGs were significantly enriched in neutrophil extracellular trap formation, alcoholic liver disease, and neuroactive ligand-receptor interaction, whereas the downregulated DEGs were significantly enriched in malignancy-related pathways, including pathways in cancer, BC, PI3K-Akt signaling pathway, MAPK signaling pathway, and focal adhesion (Fig. 2D). The above results suggest that these DEGs were mainly enriched in entries associated with extracellular matrix remodeling, and remarkably, downregulated DEGs were significantly enriched in various signaling pathways involved in the regulation of malignant tumor pathogenesis.

PPI network and SPARCL1 identified as a prognostic-related hub gene

The PPI network contained 332 nodes and 2,972 edges (Fig. S1), and we selected the top three hub modules identified by the MCODE plug-in for display (Figs. 3A–3C), where each node represents one DEG and the edges between nodes represent interactions. The top 30 hub genes with a high degree of connectivity were identified from the network (Fig. 3D). The expression of these genes and their modules are listed in Table 2. Among the 30 hub genes, SERPINA1 and VCAN were found to be upregulated in BC, whereas ALB, IGFBP3, SPARCL1, and FSTL1 were downregulated and no significant difference was discovered in the expression of other hub genes. Survival analyses showed favorable OS and RFS in patients with BC with upregulated SERPINA1 and SPARCL1 expression (Fig. 4). However, it was paradoxical that SERPINA1 showed high expression in both breast cancer and liver metastasis, but predicted a better patient prognosis. Therefore, we selected SPARCL1 for further evaluation, as the above results also suggested that it was likely to be a diagnostic and prognostic biomarker for BC.

Validation of SPARCL1 and clinicopathological correlation analysis

To affirm the accuracy of the above results, we externally validated the independent dataset GSE58708. Following differential expression analysis of the GSE58708 dataset, 683 DEGs were obtained, including 410 upregulated and 273 downregulated genes (Fig. S2). SPARCL1 expression was also significantly downregulated in liver metastasis compared with that in BC tissues (log₂FC = −2.617; Fig. S3A). In the TCGA–BRCA dataset, SPARCL1 expression in BC tissues was also significantly downregulated compared with that in normal mammary gland tissue (Fig. S3B). The correlation between SPARCL1 levels and the clinical characteristics of patients with BC was further analyzed and showed that the lower SPARCL1 expression was significantly related to age, tumor-node-metastasis (TNM) stage, estrogen receptor (ER) status, progesterone receptor (PR) status, histological type, molecular type, and living status, whereas no significant difference in node stage and Her-2 status was observed (Table 3). The above findings suggest that SPARCL1 expression was downregulated both in BC and liver metastasis, and its low expression correlates with younger age, high TNM stages, and HR negative status of patients, implying that SPARCL1 may play a tumor suppressive role in BC patients.

GSEA

From the viewpoint of enrichment of gene sets, finding the effects of subtle changes on biological pathways or functions is easier in theory (Subramanian et al., 2005). Therefore, we performed GSEA based on the different expression levels of SPARCL1 in BC. The results suggested that 50 gene sets were upregulated in the low expression SPARCL1 phenotype, 47 gene sets were significantly enriched at FDR <25%, and 40 gene sets were significantly enriched at nominal p-value <5%. Figure 5 displays six signaling pathways (DNA replication, cell cycle, oxidative phosphorylation, homologous recombination, spliceosome, and proteasome) that were significantly enriched when SPARCL1 was downregulated in BC, revealing the potential molecular mechanisms by which SPARCL1 participates in BC occurrence and progression.

SPARCL1 downregulation in BC tissues and cells

Total RNA was extracted from tissues and cell lines for RT-qPCR to validate SPARCL1 expression in BC. In comparison with the paired adjacent normal breast tissues, SPARCL1 was significantly downregulated in BC tissues (N = 30; Fig. 6A). In comparison with the normal breast epithelial cell line, SPARCL1 expression was also decreased in the BC cell lines (Fig. 6B). These results are consistent with those from our bioinformatics analysis, suggesting that SPARCL1 is indeed expressed significantly less in BC compared to in normal breast tissue.

SPARCL1 knockdown-induces proliferation and migration of BC cells in vitro

To explore the biological function of SPARCL1, we used shRNAs to knockdown SPARCL1 in BC cells. The results suggest that shRNAs could significantly downregulate the level of SPARCL1 in BC cells compared to sh-NC (Fig. 7A). The results of colony formation and MTT assays suggested that SPARCL1 inhibition promoted the colony formation and proliferation of BC cells (Figs. 7B and 7C). Cell migration is an essential step in tumor progression and metastasis. In comparison with the control group (sh-NC), the healing ability of SPARCL1-knockdown MDA-MB-231 cells was significantly enhanced (Fig. 7D) and cell migration was significantly increased (Fig. 7E). Collectively, these results suggested that downregulation of SPARCL1 promoted the colony formation, proliferation, and migration of BC cells.

Overexpression of SPARCL1 inhibited the proliferation and migration of BC in vitro

To further evaluate the effects of overexpression of SPARCL1 on the biological function of BC cell lines, SPARCL1 overexpression plasmid was transfected into BC cells. The results of RT-qPCR showed that the expression of SPARCL1 was significantly upregulated in the SPARCL1 OE group compared to the control (vector) group (Fig. 8A). MTT assay revealed that the proliferation of BC cells was significantly decreased after the overexpression of SPARCL1 (Fig. 8C), and the results of the colony formation assay showed that overexpression of SPARCL1 inhibited the colony formation of BC cells (Fig. 8B). The results of the wound-healing assay showed that there was a significant decrease in cell migration area after SPARCL1 overexpression compared to the control group (Fig. 8D). The results of Transwell assay showed that the number of cells in the lower chamber was significantly reduced after SPARCL1 overexpression (Fig. 8E). These findings demonstrated that the colony formation, proliferation, and migration of BC cells were inhibited after elevating SPARCL1 expression.