Human breast cancer cell lines were purchased from the Cell Resource Center of Beijing Xiehe (Beijing, China) and cultivated at 37°C in an atmosphere containing 5% CO₂. MCF-7 and MDA-MB-231 cells were maintained in high-glucose Dulbecco's modified Eagle's medium (Gibco; Thermo Fisher Scientific, Inc., Waltham, MA, USA) supplemented with 10% fetal bovine serum (Hyclone; GE Healthcare, Logan, UT, USA) and 100 U/ml penicillin (Thermo Fisher Scientific, Inc.). Cells were seeded in 6-well plates and were grown overnight, in order to reach 60–70% confluence prior to transfection. Transfection of MCF-7 and MDA-MB-231 cells with Myc-USP18 plasmid (250 ng/μl, 8 μl; RC202374; OriGene Technologies, Inc., Beijing, China) or vector plasmid (500 ng/μl, 4 μl) was conducted at room temperature for 30 min using Lipofectamine^® 2000 transfection reagent (Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer's protocol. Prior to detecting USP18 expression, cells were treated with 2 μM 5′Aza (Sigma-Aldrich; Merck KGaA, Darmstadt, Germany) for 3 days at 37°C; medium and 5′Aza were replenished every day.

A total of 24 h post-transfection, breast cancer cells were serum starved overnight, and were then treated with 100 ng/ml EGF (Sigma-Aldrich; Merck KGaA) for 5 min at 37°C, in order to detect the effects of USP18 overexpression on the EGF-mediated AKT signaling pathway. Furthermore, LY294002 stock solution (10 mM; Merck KGaA) was diluted with dimethyl sulfoxide, and 24 h post-transfection, cells were treated with 10 μM LY294002 for 24 h at 37°C.

The CpG island of the USP18 gene was predicted using The University of California Santa Cruz (UCSC) Genome Browser database (http://genome.ucsc.edu) and the bisulfite-conversion-based methylation PCR primers were designed using the MethPrimer database (http://www.urogene.org/methprimer/). Genomic DNA was extracted for methylation analysis from MCF-7 and MDA-MB-231 cells in culture using the Genomic DNA Miniprep kit (Sigma-Aldrich; Merck KGaA). Genomic DNA (1 μg) was then modified with sodium bisulfite using the DNA Bisulfite Conversion kit (Tiangen Biotech Co., Ltd., Beijing, China) according to the manufacturer's protocol. MSP was run in a total volume of 20 μl; MSP reactions were subjected to initial incubation at 95°C for 5 min, followed by 35 cycles at 94°C for 20 sec, annealing at 60°C for 30 sec and 72°C for 20 sec. Final extension was performed at 72°C for 5 min. MSP products were separated by 2% agarose gel electrophoresis and were visualized following ethidium bromide (Thermo Fisher Scientific, Inc.) staining. The following primers were used: Unmethylated forward, 5′-TGTTTTTTATTTAGTGGAAAATGA-3′ and reverse 5′-ACCAAAAACCAAACTAAAACCAA-3′; and methylated forward, 5′-CGTTTTTTATTTAGTGGAAAACGA-3′ and reverse, 5′-GAAAACCGAACTAAAACCGAA-3′. Normal blood lymphocyte DNA acted as a negative control and was obtained from peripheral blood of three healthy volunteers (two female volunteers aged 45 and 46, and one male volunteer aged 37). The present study was approved by the Research Ethics Committee of Qilu Hospital of Shandong University (Qingdao, China) and was conducted according to the Declaration of Helsinki. All volunteers provided written informed consent.

Microarray data were obtained from seven GEO datasets that were accessed from the National Centers for Biotechnology Information GEO database (http://www.ncbi.nlm.nih.gov/geo/), which serves as a public repository for gene expression datasets. The accession numbers were as follows: GSE66695, GSE72653, GSE61500, GSE23938, GSE2528, GSE61499 and GSE61500. Differentially expressed genes were determined using GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/), which is an interactive web tool that compares two groups of samples under the same experimental conditions and can analyze almost any GEO series.

Total RNA was isolated using TRIzol^® (Invitrogen; Thermo Fisher Scientific, Inc.) and cDNA was synthesized using a reverse transcription kit (Bio-Rad Laboratories, Inc., Hercules, CA, USA), according to the manufacturer's protocol. qPCR was performed using SYBR Green master mix on a HT7500 system (Applied Biosystems; Thermo Fisher Scientific, Inc.). The following primers were used: USP18 forward, 5′-AACGTGCCCTTGTTTGTCCAA-3′ and reverse 5′-GAGTCCTTCACCCGGATCGTA-3′; Skp2 forward, 5′-ATGCCCCAATCTTGTCCATCT-3′ and reverse 5′-CACCGACTGAGTGATAGGTGT-3′; and GAPDH forward, 5′-CTGGGCTACACTGAGCACC-3′ and reverse 5′-AAGT GGTCGTTGAGGGCAATG-3′. qPCR was conducted as follows: 95°C for 10 min, followed by 40 amplification cycles at 95°C for 15 sec, 60°C for 15 sec and 72°C for 30 sec. The FC (fold change) of gene expression was calculated using the quantification cycle (Cq) method with the following formula: FC = 2^−ΔΔCq, where ΔΔCq = (Cq_target − Cq_reference) sample − (Cq_target − Cq_reference) control (12).

Cells were solubilized in lysis buffer [50 mM Tris (pH 7.4), 150 mM NaCl, 10% glycerol, 1 mM EDTA, 1 mM EGTA, 1% SDS, 1 mM Na₃VO₄, 1 mM DTT and 10 mM NaF]. Protein samples were obtained from whole cell lysates, mixed with 2X loading buffer [1:1; 20 mM Tris-HCl (pH 6.8), 100 mM DTT, 2% SDS, 20% glycerol and 0.016% bromophenol blue] and incubated at 99°C for 15 min. Protein content determination within cell lysates was conducted using the bicinchoninic acid (BCA) method (Pierce BCA protein assay kit; Thermo Fisher Scientific, Inc.). After quantification, protein samples (~10 μl/300 μl whole cell lysates, ~20 μg) were separated by 10% SDS-PAGE and transferred to nitrocellulose membranes. Subsequently, the blots were blocked in blocking buffer [5% non-fat dry milk in Tris-buffered saline containing 0.1% Tween (TBST)] for 1 h at room temperature, and were then incubated with primary antibodies (1:1,000) in blocking buffer overnight at 4°C. The blots were washed three times with TBST and were incubated for 1 h at room temperature with horseradish peroxidase-conjugated anti-rabbit (A4914) immunoglobulin G secondary antibodies (1:3,000; Sigma-Aldrich; Merck KGaA) in blocking buffer. Finally, the blots were washed three times with TBST and visualized by enhanced chemiluminescence (ECL) using an ECL kit (Applygen Technologies, Inc., Beijing, China). The following primary antibodies were used in the present study: Anti-Myc (#2278, 1:3,000), anti-USP18 (#4813, 1:1,000), anti-AKT (#4685, 1:1,000), anti-phosphorylated (p) AKT_308 (#4056, 1:500), anti-pAKT_473 (#3787, 1:1,000), anti-EGFR (#4267, 1:1,000), anti-Skp2 (#2652, 1:1,000) and anti-GAPDH (#5174, 1:5,000; all from Cell Signaling Technology, Inc., Danvers, MA, USA).

The mRNA (RNA Seq v2) and protein (reverse phase protein array) data were downloaded from https://www.synapse.org, and the clinical information for patients in TCGA_breast cancer (BRCA) dataset was downloaded from cBioPortal database (www.cbioportal.org); these data were used to analyze differential mRNA expression, correlation, prognosis and gene set enrichment. Methylation data and the results of correlation analysis between methylation percentage and mRNA level were downloaded from MethHC database (http://methhc.mbc.nctu.edu.tw/php/index.php). The protein expression data were obtained from The Cancer Proteome Atlas database (http://www.tcpa-portal.org/tcpa/).

Cell proliferation was assessed using the Cell Counting Kit-8 (CCK8) assay (Promega Corporation, Madison, WI, USA), according to the manufacturer's protocol. The transfected cells were plated in 96-well plates (2,000 cells/well) and cell proliferation was detected every 24 h. Briefly, 10 μl CCK8 solution was added to each well and incubated for 1 h at 37°C. The solution was then measured spectrophotometrically at 450 nm.

A total of 2 days post-transfection with Myc-USP18 or empty vectors, 500 cells were subsequently seeded into each well of 6-well plates with full culture medium. After 7–10 days, cells were fixed with 70% ethanol for 10 min and stained with 0.5% crystal violet solution. Nikon Eclipse 600 photomicroscope (Nikon Corporation, Tokyo, Japan) was used to observe colony formation and to capture images. The experiment was conducted in three independent triplicates.

The association between phenotypes, pathways and USP18/Skp2 expression was analyzed using GSEA v2.2 (http://www.broad.mit.edu/gsea/). GSEA calculates a gene set Enrichment Score that estimates whether genes from a pre-defined gene set (obtained from the Molecular Signatures Database, MSigDB, http://software.broadinstitute.org/gsea/msigdb/index.jsp) are enriched among the highest- or lowest-ranked genes or distributed randomly. Default settings were used. Thresholds for significance were determined by permutation analysis (1,000 permutations). false discovery rate (FDR) was calculated. A gene set was considered significantly enriched when the FDR score was <0.05.

Cells transfected with USP18 or empty vectors were fixed in 70% ethanol for 12 h at 4°C, treated with 0.1 mg/ml RNase A for 30 min at 37°C, and stained with 50 mg/ml propidium iodide for 30 min at 37°C in the dark (PI; BD Biosciences, Franklin Lakes, NJ, USA). The cells were then analyzed by flow cytometry using FACSCalibur (BD Biosciences) (13), and cell cycle distributions were calculated using CellQuest software version 5.1 (BD Biosciences). These experiments were performed three times in triplicate.

To explore the biological function of USP18 in vivo, the Database for Annotation Visualization and Integrated Discovery database (https://david.ncifcrf.gov/home.jsp) was used to map the KEGG pathways and GO terms associated with the differentially expressed genes in USP18 knockout (KO) mice (GEO dataset, GSE61500).

Search Tool for the Retrieval of Interacting Genes (STRING; https://string-db.org/) is a database of known and predicted protein interactions, which may aid in the comprehensive description of cellular mechanisms and functions. A PPI network of the important differentially expressed genes with GO annotations in USP18 KO mice was constructed using the STRING database.

JASPAR database (http://jaspar.binf.ku.dk/) was used to predict transcription factor binding sites in DNA sequences.

All experiments were conducted in three independent triplicates and data are presented as the means ± standard error of the mean. Statistical analyses were performed using SPSS 19.0 (IBM Corp., Armonk, NY, USA) and GraphPad Prism 5 (GraphPad Software, Inc., La Jolla, CA, USA). Group distributions were compared using Student's t-test, or one-way analysis of variance followed by post hoc Bonferroni tests for the analysis of multiple comparisons. Correlations were analyzed using the Pearson correlation analysis. In addition, log-rank test for the generated Kaplan-Meier curves was conducted to evaluate the association between gene expression and survival rate. P<0.05 was considered to indicate a statistically significant difference.

Aberrant promoter methylation always affects transcriptional gene expression. Sequence analysis using the UCSC database (http://genome.ucsc.edu/index.html) revealed a CpG island located within the USP18 promoter region (Fig. 1A). As shown in Fig. 1B, the promoter region of USP18 exhibited poor methylation in MCF-7 and MDA-MB-231 cells, and the USP18 methylation level was significantly downregulated in breast cancer compared with in normal tissues (Fig. 1C), and primer efficiency was verified by a negative control (normal blood lymphocyte DNA; data not shown). In addition, demethylation was induced by treatment with 5′Aza, which successfully increased the mRNA and protein expression levels of USP18 in MCF-7 and MDA-MB-231 cells (Fig. 1D and E). Subsequently, the expression levels of USP18 were detected in tumor tissues compared with in adjacent tissues using public datasets; the results indicated that USP18 was clearly enhanced in breast cancer tissues (Fig. 1F and G). To further evaluate whether a correlation existed between promoter DNA methylation and USP18 mRNA expression in breast cancer, the present study further detected USP18 promoter DNA methylation and mRNA expression data from the MethHC database (http://methhc.mbc.nctu.edu.tw/php/index.php). Linear regression analysis demonstrated a significant negative correlation between mRNA expression and promoter methylation in USP18 (Fig. 1H). These findings suggested that reduced methylation of USP18 in the promoter region may induce USP18 upregulation in breast cancer.

To elucidate the biological function of USP18 in MCF-7 and MDA-MB-231 cells, the cells were transfected with a human USP18 vector. Myc protein expression was confirmed in Myc-USP18-transfected MCF-7 and MDA-MB-231 cells by western blotting (Fig. 2A). Subsequently, the effects of USP18 on cell proliferation were examined by CCK8 and colony formation assays. As shown in Fig. 2B, cell proliferation was significantly promoted by USP18 in MCF-7 and MDA-MB-231 cells (P<0.0001, Fig. 2B). In addition, USP18 expression promoted colony-forming ability in MCF-7 and MDA-MB-231 cells (P<0.0001, Fig. 2C). In order to evaluate the role of USP18 in clinical patients with breast cancer, GSEA results revealed that USP18 mRNA expression was positively correlated with the epithelial cell proliferation gene set (Fig. 2D). These results indicated the enhancing effects of USP18 on breast cancer growth.

The present study further examined the effects of USP18 on cell cycle progression by flow cytometry. Following PI staining, the proportion of MCF-7 cells in S phase was increased from 7.15 to 9.22% in response to USP18 overexpression (Fig. 3A). Similar findings were observed in MDA-MB-231 cells in response to USP18 overexpression; the number of cells in S phase was elevated from 14.53 to 18.35% (Fig. 3B). GSEA results revealed that USP18 expression was negatively correlated with the cell cycle arrest gene set and the apoptosis gene set in patients with breast cancer (Fig. 3C). These results demonstrated that USP18 may participate in breast cancer development by promoting cell cycle progression and inhibiting apoptosis; these findings are similar to those presented in a previous study (14).

The present study demonstrated that USP18 may regulate cell cycle progression and apoptosis in breast cancer cells; therefore, further investigation was conducted to determine the physiological and pathological roles of USP18. KEGG pathway enrichment analysis and GO enrichment analysis were performed in USP18 KO mice and mouse mammary tumor models. As shown in Fig. 4A, the AKT/FoxO signaling pathway was one of the most dysregulated pathways, as determined by KEGG pathway analysis, in USP18 KO mice. In addition, as shown in Fig. 4B, the majority of the enriched GO terms belonged to the biological process category. Cell death, cell cycle and cell proliferation were the most relevant biological processes. Furthermore, to obtain PPI data, differentially expressed genes, which were involved in cell death, cell cycle and cell proliferation biological processes, were uploaded to the STRING database. These samples with a PPI score >0.7 were used to construct the PPI network (Fig. 4C). Skp2, transforming growth factor β receptor 1, interleukin 6, polo-like kinase 2, Cbl proto-oncogene, MDM2 proto-oncogene, kinase insert domain receptor and aurora kinase B were the main nodes in this network. In addition, these genes were enriched in the KEGG pathway and GO terms. To further elucidate the function of USP18 in tumor progression, MMTV-PymT (GSE23938) and MMTV-neu (GSE2528) mouse mammary tumor models were used. As shown in Fig. 4D and E, USP18 was significantly upregulated in mouse models of mammary tumors. These results demonstrated that USP18 may serve an important role in breast cancer progression.

To investigate the main pathways in breast cancer, KEGG pathway analysis was conducted on differentially expressed genes in breast cancer tissues compared with matched normal tissues. As shown in Fig. 5A, the AKT signaling pathway was the most significantly enriched pathway, and the cell cycle process was dysregulated in clinical patients. KEGG pathway analysis validated that the AKT signaling pathway serves an important role in the occurrence of breast cancer. To evaluate the association between USP18 and the AKT pathway in breast cancer, GSEA was conducted; the results revealed that USP18 expression was positively correlated with the AKT signaling pathway (Fig. 5B). In previous studies (14,15), USP18 was reported to upregulate EGFR, which is the main upstream kinase of AKT that promotes activation of AKT. Consistent with the results of the present study, and those of previous studies (14,15), it was demonstrated that the protein expression levels of EGFR were upregulated in the USP18 high expression group from TCPA_BRCA (reverse phase protein array) dataset (Fig. 5C). Furthermore, pAKT expression was positively correlated with expression of the downstream target gene Skp2 in TCPA_BRCA dataset. It is well known that Skp2 is a key cell cycle regulatory factor (Fig. 5D). In addition, GSEA results demonstrated that activation of the AKT signaling pathway was positively correlated with Skp2 expression (Fig. 5E). In addition, Skp2 was negatively associated with p27 protein expression, which is the specific substrate of Skp2 and a cell cycle inhibitor in TCPA_BRCA dataset (Fig. 5F). These results demonstrated that USP18 may have an important role in clinical breast cancer progression, which is associated with the AKT/Skp2 pathway.

To confirm the effects of USP18 on EGFR expression and the AKT/Skp2 pathway in breast cancer cells, MCF-7 and MDA-MB-231 cells were transiently transfected with USP18 and empty vector plasmids. As shown in Fig. 6A, USP18 overexpression upregulated the protein expression levels of EGFR and Skp2, and increased EGF-mediated AKT phosphorylation in MCF-7 and MDA-MB-231 cells. Following 5 min of EGF stimulation, the activation levels of AKT in USP18-transfected cells were higher than those in empty vector plasmid-transfected cells. In addition, Skp2 mRNA expression was significantly downregulated in USP18 KO mice (Fig. 6B). In order to further evaluate the association between USP18 and Skp2 in breast cancer, USP18 and Skp2 mRNA expression, and p27 protein expression data were retrieved from TCGA dataset. Linear regression analysis demonstrated a significant positive correlation between USP18 and Skp2 expression (Fig. 6C), and a significant negative correlation between USP18 and P27 in TCPA_BRCA dataset (Fig. 6D). Notably, the effects of USP18 on Skp2 expression were attenuated following treatment with the phosphoinositide 3-kinase inhibitor LY294002 (Fig. 6E). In a previous study, FoxO3, which is a negative regulator of Skp2, was reported to be negatively regulated by AKT (16). As presented in Fig. 6F, the results of a bioinformatics analysis demonstrated that the promoter region of Skp2 contains a FoxO3 binding site. In addition, among patients with high FoxO3 expression, USP18 was positively associated with Skp2 expression; however, among patients with low FoxO3 expression, USP18 was not correlated with Skp2 expression, thus indicating that USP18 may upregulate Skp2 in a FoxO3 (AKT)-dependent manner in breast cancer. These findings suggested that USP18 may increase EGFR expression and promote AKT/Skp2 pathway activation in breast cancer.

In order to gain further support for the role of USP18 in breast cancer progression, TCGA dataset was used. Notably, the expression levels of USP18 were positively associated with increasing TNM stages of breast tumors (Fig. 7A). To further analyze the prognostic value of USP18 in breast cancer tissues, a Kaplan-Meier analysis was conducted on TCGA dataset. As shown in Fig. 7B, high USP18 expression was significantly associated with a worse disease-free survival (log-rank P<0.05) of patients. In addition, the expression levels of USP18 in breast cancer tissues from various breast cancer subtypes were analyzed from TCGA dataset. As shown in Fig. 7C, compared with in matched normal tissues, USP18 was significantly upregulated in HER2-positive patients, which is consistent with the role of Skp2, as reported in a previous study (17). In addition, the expression of USP18 in HER2-positive breast cancer was increased compared with in other types of breast cancer. Further stratification of the patient groups, based on the expression of USP18/Skp2, demonstrated that patients with high expression of both USP18 and Skp2 had the worst survival rate (Fig. 7D). These results suggested that USP18/Skp2 may act as a potential biomarker in breast cancer.