In the BC group, plasma samples were collected from patients with Stage I breast cancer (T1N0M0) (1), who did not have other systemic diseases or cancer at the time of their initial diagnosis, prior to receiving any treatment. BC tissues from these patients were collected at the time of surgery. All the patients involved were women. The age range of patients was 20–63 years.

In the control group, plasma samples were collected from patients with benign breast disease, including breast fibroadenoma and mammary adenosis; these patients also did not have other systemic diseases or cancer at the time of collection. The plasma samples were collected at the time of diagnosis prior to surgery, and normal breast tissue from the same patient was collected during surgery, which served as the benign lesion sample. The diagnoses of the patients were confirmed by clinical and pathological findings. A summary of the clinicopathological characteristics of control and BC groups is presented in Table I. All experiments were performed in accordance with the approved guidelines and regulations of the Shenzhen People's Hospital. All experimental protocols were approved by the Ethics Committees of Shenzhen People's Hospital. Written informed consent was obtained from all participants, in accordance with the Declaration of Helsinki.

Venous blood samples (6 ml/patient) were collected in tubes containing 10.8 mg EDTA as the anticoagulant (BD Biosciences). Plasma was separated from blood cells by centrifugation at 1,000 × g for 10 min at 4°C. The plasma layer was further centrifuged at 16,000 × g for 10 min at 4°C prior to extracellular vesicle extraction. Tissue samples were identified by a clinical pathologist and cut into 1 cm³ sections for storage at −80°C until further processing.

Extracellular vesicles were isolated from the plasma using ExoQuick precipitation (System Biosciences, LLC), according to the manufacturer's protocol. Briefly, the extracellular vesicles were precipitated by incubation; the ExoQuick extracellular vesicle precipitation reagent was added at 4°C for 60 min. The extracellular vesicle pellet was collected by centrifugation at 1,500 × g for 10 min at 4°C and resuspended in 10 mM PBS to four times the original plasma volume.

A copper mesh was placed on a clean wax plate and 100 µl of the extracellular vesicle suspension was added. After 4 min, the copper mesh was removed and placed in 2% phosphotungstic acid for 5 min at room temperature. The mesh was dried on a filter paper. The mesh was examined under a transmission electron microscope (HT7700; Hitachi, Ltd.) and images were captured.

The extracellular vesicle pellet was dissolved in RIPA protein lysis buffer (Sangon Biotech, Co., Ltd.), and the protein concentration was determined using the bicinchoninic acid protein assay kit (Thermo Fisher Scientific, Inc.). The proteins were separated on an SDS-PAGE gel (the concentrations of separating gel and stacking gel were 12 and 5% respectively). A total of 30 µg of protein was loaded per lane, while the protein ladder (Thermo Fisher Scientific, Inc.; 26616) was loaded 3 µl per lane. The proteins were electrotransferred onto a PVDF membrane. Membranes were blocked with 5% skimmed milk in TBST at room temperature for 1 h. The membrane was then incubated with the following primary antibodies at 4°C overnight: Tumor susceptibility gene 101 protein (TSG101; Abcam; ab83, dilution 1:1,000), CD9 (Abcam; ab92762, dilution 1:1,000) and heat shock cognate 70 protein (HSC70; Santa Cruz; sc7298, dilution 1:300). The membrane was then incubated with the HRP-conjugated corresponding secondary antibodies (Proteintech; SA00001-1, dilution 1:1,000 and CST; 7074P2, dilution 1:10,000) at room temperature for 1 h and treated with enhanced chemiluminescence detection reagents (4A Biotech; 4AW012-100). The specific protein bands were visualized and recorded on film using a chemiluminescence imaging system (2000; Tanon).

Total RNA was extracted from the extracellular vesicles and tissues using TRIzol^® reagent (Life Technologies; Thermo Fisher Scientific, Inc.), according to the manufacturer's protocol. The purity of the isolated RNA was determined according to the ratio of radiation absorbance at 260 nm to that at 280 nm, using a Nanodrop ND-1000 system (Thermo Fisher Scientific, Inc.); 200 ng total RNA was used for the construction of the miRNA library using the NEBNext Multiplex Small RNA Library Prep Set for Illumina kit (New England BioLabs, Inc.). Each RNA sample underwent adaptor ligation, cDNA synthesis and polymerase chain reaction amplification. The quality of the RNA library was verified by Agilent 2100 Bioanalyzer (Agilent Technologies, Inc.). The library was sequenced on the Hiseq 2500 (Illumina Inc.,). RNA-seq datasets were generated with 250 bp paired-end reads with a 10 million-read sequencing depth of the multiplexed samples.

Fastx_toolkit (http://hannonlab.cshl.edu/fastx_toolkit/) was used to remove low quality reads from the sequencing data. A quality statistic was performed using FastQC software (version 0.13.0) (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) as a quality control tool for the high-throughput sequencing data. All reads were compared with the reference sequence, human mature miRNA from miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/php/index.php), using the FastQC high precision sequence alignment algorithm (15). Gene expression was quantified in reads per kilobase of transcript per million mapped reads (RPKM) (16,17).

Differential miRNA expression analysis was performed. EdgeR was used to determine the differentially expressed miRNAs (18). P<0.05 and the absolute value of log2 (FoldChange) <1 (|log2 (FoldChange)|>1) were considered to indicate a statistically significant difference. Clustering analysis was performed using a hierarchical clustering method based on the expression level of differentially expressed miRNAs [log10(RPKM+1)]. miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/php/index.php) was used to identify potential human miRNA target genes. Gene Ontology (GO) analysis was performed using topGO (version 2.18.0) (https://bioconductor.org/packages/release/bioc/html/topGO.html) to investigate the biological processes, cellular components and specific molecular functions of the identified differentially expressed coding genes. Pathway analysis was performed using KOBAS (kobas2.0-20150126) (18) to determine the involvement of co-expressed genes in different biological pathways according to the Kyoto Encyclopedia of Genes and Genomes (KEGG). The hypergeometric test was used for statistical analysis. miRNA-target gene-network analysis was performed using mirWalk (http://mirwalk.umm.uni-heidelberg.de/) (19). A receiver operating characteristic (ROC) curve was generated to calculate the relationship between sensitivity and specificity for the disease group compared with the healthy controls. Data analysis and ROC curve analysis was performed using SPSS (version 19.0; IBM Corp.) and GraphPad Prism (version 7; GraphPad Software, Inc.).

The Kaplan-Meier method was used to compute the survival analyses using OncoLnc (http://oncolnc.org/) and Oncomir (http://www.oncomir.org). OncoLnc is a tool for interactively exploring survival correlations, and for downloading clinical data coupled to expression data for mRNAs, miRNAs or long noncoding RNAs (lncRNAs). OncoLnc contains survival data for 8,647 patients from 21 cancer studies performed by The Cancer Genome Atlas (TCGA; http://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga), along with miRNAs from TCGA. OncomiR is an online resource for exploring miRNA dysregulation in cancer. Using combined miRNA-seq, RNA-seq and clinical data from TCGA, it is able systematically perform statistical analyses to identify dysregulated miRNAs that are associated with tumor development and progression in the majority of cancer types (20).

The analysis of differentially expressed miRNAs was evaluated using R studio software (version 3.4.2) with the edgeR package (21,22). Differential expression was assessed for each gene using an exact test analogous to Fisher's exact test, but adapted for over-dispersed data (22). P<0.05 was considered to indicate a statistically significant difference. Data are expressed as the median with range, and were analyzed using GraphPad Prism software (version 7.0). Pearson correlation analysis was performed to assess the correlation between age and miRNA expression profile using SPSS (version 19). The scatter plot was drawn using SPSS (version 19).

The present study initially assessed the extracellular vesicles that were isolated from the plasma samples. Transmission electron microscopy revealed spherical vesicles of 30–100 nm in diameter in all samples (Fig. 1A), suggesting that the present study had successfully purified extracellular vesicles from the plasma. According to the standard curve, the mean protein concentration of extracellular vesicles in the breast cancer group was 24.16 µg/µl, and was 23.84 µg/µl in the control group. Western blotting was performed to detect three conventional extracellular vesicle protein markers, namely TSG101, CD9 and HSC70 (15). With an equal amount of protein loaded in each lane, these three extracellular vesicle protein markers were all highly enriched in the isolated extracellular vesicles relative to the plasma (Fig. 1B). These results confirmed the successful purification of intact extracellular vesicles from all plasma samples.

The optical density (OD)260/OD280 value range of total RNA extracted from extracellular vesicles in each samples ranged between 1.70 and 2.00, which was indicative of high RNA purity. Differential miRNA expression analysis was performed to identify the candidate extracellular vesicle-packaged miRNAs. A series of bioinformatics analyses were performed to reveal the relationship between extracellular vesicle-packaged miRNAs and the pathogenesis and progression of early-stage breast cancer. The edgeR method, based on negative binomial distribution, was adopted. When compared with the control group, there were 41 differentially expressed miRNAs out of 1,111 miRNAs in the blood samples of the BC group, of which 21 were significantly upregulated, while 20 were significantly downregulated miRNAs (P<0.05). In addition, when compared with the control group, there were 181 differentially expressed miRNAs in the tissue samples of the BC group, of which 102 were significantly upregulated, while 79 were significantly downregulated (P<0.05) (Table II). A volcano plot illustrating the distribution of differential miRNA between breast cancer and control group is presented in Fig. 2A.

Hierarchical clustering was performed using R software based on the RPKM values of differentially expressed miRNAs. Hierarchical clustering analysis revealed differences in the extracellular vesicle-packaged miRNA expression in the blood samples between the breast cancer and control groups (Fig. 2B).

miRTarBase was used to identify potential human miRNA target genes of differentially expressed miRNAs. GO was used to classify the functions of the target genes from three structured networks, including biological processes, cellular components and molecular function.

In the present study, GO enrichment was performed in blood and tissue samples comparing the early-stage breast cancer group with the control group (Fig. 3). In the biological process analysis, the 20 most prominent terms in the blood and tissue samples were compared, and 9 out of 20 terms were identical, including ‘cellular macromolecule metabolic process’, ‘regulation of gene expression’ and ‘RNA metabolic process’. In the cellular components analysis, the present study contrasted the 20 most prominent terms in the blood and tissue samples, and 19 out of 20 terms were identical, including ‘nucleoplasm’, ‘cytosol’ and ‘nuclear lumen’. In the molecular function analysis, the present study contrasted the 20 most prominent terms in the blood and tissue samples, and 8 out of 20 terms were observed to be identical, including ‘enzyme binding’, ‘DNA binding’ and ‘nucleic acid binding’.

KEGG pathway analysis was performed to identify significantly affected pathways. The results indicated that in the blood samples, the target genes of differentially expressed miRNAs were enriched in ‘ErbB signaling pathway’, ‘cell cycle’, ‘signaling pathways regulating pluripotency of stem cells’ and ‘pathways in cancer’. In the tissue samples, the target genes of differentially expressed miRNAs were enriched in ‘PI3K-Akt signaling pathway’, ‘MAPK signaling pathway’, ‘mTOR signaling pathway’ and ‘pathways in cancer’. Comparison of the 20 most prominent terms of the KEGG pathway analysis in the blood and tissue samples revealed that 11 out of 20 terms were the same, including ‘signaling pathways regulating pluripotency of stem cells’, ‘pathways in cancer’, ‘cell cycle’, ‘proteoglycans in cancer’ and the ‘mTOR signaling pathway’ (Fig. 4).

Comparison of the early BC group with the control group revealed that six miRNAs were significantly altered in the blood and tissue analyses (Table III and Fig. 5). miRNA-205-5p was significantly upregulated in the blood samples (P=0.011650), yet was significantly downregulated in the tissue samples (P=0.000304). Excluding miRNA-205-5p, a total of five miRNAs exhibited identical changing trends in the blood and tissue analyses. These miRNAs included miRNA-375, miRNA-24-2-5p, miRNA-548b-5p, miRNA-655-3P and miRNA-376b-5p. Network analysis of miRNAs and target genes was performed based on the five miRNAs, and it was revealed that these five miRNAs were associated with the same target genes and regulated each other, thereby affecting biological processes (Fig. 6). At the same time, the statistical analysis of patient's age showed that the mean age between the breast cancer group and control group was significantly different (P=0.019). However, bivariate correlation analysis revealed that there was no significant correlation between these 5 miRNAs (miR-375, miR-24-2-5p, miR-548b-5p, miR-655-3p and miR-376b-5p) and age in the blood groups (Table SI and Fig. S1).

To evaluate the utility of extracellular vesicle-packaged miRNA levels for discriminating cases of early-stage BC from controls, ROC curve analysis was performed. A total of five miRNAs identified in the primary analyses were highly differentially expressed between patients with early-stage BC and controls (Fig. 7). The highest area under the curve (AUC) for a single miRNA was achieved with miRNA-548b-5p [AUC=0.785; 95% confidence interval (CI)=0.585–0.984; P=0.022].

In addition, the different combinations of these miRNAs were investigated (these results are presented in Table SI), and the highest AUC was achieved via the combination of miRNA-375, miRNA-24-2-5p, miRNA-548b-5p and miRNA-655-3p (AUC=0.808; 95% CI=0.629–0.986; P=0.013). An increase in the AUC was observed when miRNA-375 was combined with other miRNAs (Fig. 7), such as miRNA-375 + miRNA-655-3p (AUC=0.762; 95% CI=0.557–0.966; P=0.035), miRNA-375 + miRNA-655-3p + miRNA-548b-5p (AUC=0.792; 95% CI=0.605–0.979; P=0.018).

Because the highest AUC was achieved via the combination of miRNA-375, miRNA-24-2-5p, miRNA-548b-5p and miRNA-655-3p, KEGG pathway analysis of these four miRNAs was subsequently performed. As illustrated in Fig. 8, the Venn diagram had little overlap, indicating that the four miRNAs were involved in independent pathways. Taken together, these data demonstrated that extracellular vesicle-packaged miRNAs exhibit potential as diagnostic markers for early-stage breast cancer.

Based on online data from TCGA, the Kaplan-Meier curves were constructed, and log-rank test analysis was performed. Low expression of miRNA-24-2-5p was associated with a high survival rate over the period between 3,000 and 6,000 days. A similar trend was also observed in miRNA-375; low expression of miRNA-375 was associated with high survival over the period between 2,000 and 6,000 days. For miRNA-655-3p, a decrease in this miRNA in breast cancer was identified, although its lower expression was also associated with a higher survival rate. miRNA-548b-5p was decreased in the breast cancer group, and the high expression of miRNA-548b-5p was associated with a high survival rate (Fig. 9). However, it was not possible to include miRNA-376b-5p in the Kaplan-Meier analysis, as there were no expression data for miRNA-376b-5p in breast cancer within TCGA database, but the Kaplan-Meier curve and log-rank test analysis of miRNA-376b-5p was available based on cancer of the bladder, colon, pancreas, rectum and stomach. Apart from bladder cancer, the results demonstrated that the increased expression of miRNA-376b-5p was associated with a higher survival rate in the other cancer types examined (Fig. S2).