ORIGINAL RESEARCH published: 26 August 2021 doi: 10.3389/fonc.2021.627811 Circulating microRNAs in Early Breast Cancer Patients and Its Association With Lymph Node Metastases Daniel Escuin 1*†, Laura López-Vilaró 1,2†, Josefina Mora 3, Olga Bell 1, Antonio Moral 4,5, Ignacio Pérez 4, Cristina Arqueros 6, Bárbara Garcı́a-Valdecasas 7, Edited by: Teresa Ramón y Cajal 6, Enrique Lerma 1,2,5 and Agustı́ Barnadas 1,5,6,8* Angela Santoro, Università Cattolica del Sacro Cuore, 1 Clinical Oncology Research Group, Institut d'Investigacions Biomédiques Sant Pau (IIB-Sant Pau), Barcelona, Spain, Italy 2 Department of Pathology, Hospital de la Santa Creu i Sant Pau, Barcelona, Spain, 3 Department of Clinical Biochemistry, Reviewed by: Hospital de la Santa Creu i Sant Pau, Barcelona, Spain, 4 Department of General Surgery, Hospital de la Santa Creu i Sant Mariano Monzo, Pau, Barcelona, Spain, 5 School of Medicine, Universitat Autónoma de Barcelona (UAB), Cerdanyola del Vallès, Spain, University of Barcelona, Spain 6 Department of Medical Oncology, Hospital de la Santa Creu i Sant Pau, Barcelona, Spain, 7 Department of Gynecology, Gian Franco Zannoni, Hospital de la Santa Creu i Sant Pau, Barcelona, Spain, 8 Centro de Investigación Biomédica en Red Cáncer (CIBERONC), Catholic University of the Sacred Madrid, Spain Heart, Italy *Correspondence: MicroRNAs have emerged as important regulators of the metastatic process. In addition, Daniel Escuin descuin@santpau.cat circulating miRNAs appear to be surprisingly stable in peripheral blood making them ideal Agustı´ Barnadas noninvasive biomarkers for disease diagnosis. Here, we performed a proof-of-principle abarnadasm@santpau.cat † study to investigate the expression profile of circulating miRNAs and their association with These authors have contributed equally to this work the metastatic lymph node status in early breast cancer patients. Sentinel lymph node status was detected by one-step nucleic acid (OSNA) analysis. We performed RNA- Specialty section: sequencing in 16 plasma samples and validated the results by qPCR. Gene Ontology term This article was submitted to Women’s Cancer, enrichment and KEGG pathway analyses were carried out using DAVID tools. We found16 a section of the journal differentially expressed miRNAs (q < 0.01) in patients with positive SLNs. Fourteen Frontiers in Oncology miRNAs were down-regulated (miR-339-5p, miR-133a-3p, miR-326, miR-331-3p, Received: 10 November 2020 miR-369-3p, miR-328-3p, miR-26a-3p, miR-139-3p, miR-493-3p, miR-664a-5p, miR- Accepted: 04 August 2021 Published: 26 August 2021 146a-5p, miR-323b-3p, miR-1307-3p and miR-423-3p) and 2 were up-regulated (miR- Citation: 101-3pand miR-144-3p). Hierarchical clustering using differentially expressed miRNAs Escuin D, López-Vilaró L, clearly distinguished patients according to their lymph node status. Gene ontology Mora J, Bell O, Moral A, Pérez I, Arqueros C, Garcı´a-Valdecasas B, analysis showed a significant enrichment of biological processes associated with the Ramón y Cajal T, Lerma E and regulation of the epithelial mesenchymal transition, cell proliferation and transcriptional Barnadas A (2021) Circulating regulation. Our results suggest the potential role of several circulating miRNAs as microRNAs in Early Breast Cancer Patients and Its Association surrogate markers of lymph node metastases in early breast cancer patients. Further With Lymph Node Metastases. validation in a larger cohort of patients will be necessary to confirm our results. Front. Oncol. 11:627811. doi: 10.3389/fonc.2021.627811 Keywords: circulating microRNAs, breast cancer, lymph node, metastasis, surrogate biomarker Frontiers in Oncology | www.frontiersin.org 1 August 2021 | Volume 11 | Article 627811 Escuin et al. miRNAs and BC Locoregional Metastases INTRODUCTION conserved across species (17). MiRNAs regulate the expression of target genes by binding to complementary regions of messenger Breast cancer remains a common disease worldwide and the transcripts to repress their translation or regulate their second cause of cancer death in the US (1). Early diagnosis, degradation. MiRNAs are now recognized as novel post- improvements in treatment and early onset of therapy are transcriptional regulators targeting over 30% of the human important factors determining the prognosis and management genome (18). The overall emerging picture is that of a complex of patients with breast cancer. Various factors such as early age at regulation level of gene expression, in which a single miRNA menarche, late age at first birth and late age at menopause are may control hundreds of targets (19). Many cellular pathways are related to breast cancer risk. However, lymph node (LN) affected by the regulatory function of miRNAs and several affection remains the most important prognosis factor in breast human pathologies, including cancers, have been associated cancer (2). There are a number of factors associated with with misregulation of the miRNAs (16) and their metastases metastases to the LN, including tumor size, presence of (20). Numerous studies have identified widespread alterations in lymphovascular invasion, poor histological grade and age the expression of miRNAs related to human neoplasias. In breast (3, 4). Nevertheless, for a significant number of early-breast cancer, analysis of miRNA expression classified the different cancer patients it is unclear whom will develop metastases. For breast cancer molecular subtypes and correlated these with instance, about 13% of patients with favorable prognostic factors various clinicopathological factors and numerous miRNAs at diagnose will develop metastasis and the percentage increases have been shown to play a pivotal role in various steps of the to 20-30% for LN-negative patients. In contrast, 20-30% of LN- metastatic process. In addition, circulating miRNAs are positive patients will never metastasize (5), therefore it is unclear emerging as prognostic factors in breast cancer (14), but few whether distant metastases arise in a sequential manner from LN studies have correlated their expression with the LN status, the metastases or in parallel through the blood stream and whether occurrence of distant metastases and breast cancer recurrence. other factors such interactions between the tumor and the These studies have analyzed the expression of specific circulating stroma favor locoregional metastases (6). miRNAs by qPCR and have shown promising results (21–24). Most women diagnosed with breast cancer are initially treated Herein, we sought to examine the miRNA content in plasma with surgery to remove the tumor and to determine the presence samples from early breast cancer patients with known SLN and of metastases in the sentinel LNs (SLNs). This is currently the axillary LN metastatic status. We designed a proof-of-principle recommended procedure for axillary staging of early breast study to profile the expression of miRNAs by RNA-sequencing cancer. Our institution use the one-step nucleic amplification using preoperative peripheral blood from patients with early (OSNA) assay (7) to accurately measure total metastatic volume breast cancer who were not previously treated. Our results are in the SLN (8), as an alternative to intraoperative microscopy- preliminary but support the hypothesis of the existence of a based pathological assessment of the SLN. The OSNA assay is a differential miRNA expression profile in the peripheral blood rapid molecular detection of SLN metastasis based on the semi- from breast cancer patients associated with the LN status of their quantification of cytokeratin 19 (CK19) mRNA copy numbers tumors. Our data highlights the potential use of circulating (7). Thus, only patients diagnosed with more than two miRNAs as surrogate markers of locoregional metastases in macrometastatic SLN are further treated with axillary lymph breast cancer. Further studies in a larger number of samples node dissection (ALND) (9), the golden standard procedure for are warranted. invasive breast cancer. However, ALND has been questioned in recent years because of inherent morbidity following the procedure without directly contributing to survival in primary breast cancer patients (10–12) and the recognition that not all MATERIALS AND METHODS patients with nodal disease may require extensive axillary surgery (13). Patients Elucidation of breast cancer’s molecular biological features We studied 16 patients with early breast cancer treated with have had a dramatic effect on how patients are diagnosed and surgery and diagnosed for positive SLNs. All patients had treated. However, effective management of breast cancer is still confirmed diagnosis based on histopathology of tumor biopsy. difficult because of the lack of sensitive and specific biomarkers All tumors were invasive ductal carcinomas (IDCs) with or for early detection and for diseases monitoring. Accumulating without in situ component. In 2 cases, tumors were mixed and evidence in the last years has highlighted the potential use of show presence of invasive lobular carcinoma (ILC) component. peripheral blood circulating nucleic acids in breast cancer Intraoperative SLN were evaluated using the OSNA assay (7). diagnosis, prognosis and for monitoring response to anticancer None of the patients had prior treatment with surgery, therapy. Among these, circulating microRNAs (miRNAs) are chemotherapy or radiation. All patients were hormone increasingly recognized as a promising non-invasive biomarker, receptor (HR) positive, HER2 negative. We collected the given the ease with which miRNAs can be isolated and their following clinical and pathological parameters: age, structural stability under different conditions of sample menopausal status, personal and familiar disease precedents processing and isolation (14–16). and clinical follow-up, tumor stage was be determined MicroRNAs (miRNAs) are a small (19-25 nt) non-coding according to the AJCC/UICC system (25), histological grade RNAs, expressed in a wide variety of organisms and highly was determined using the Elston-Ellis grading system (26), Frontiers in Oncology | www.frontiersin.org 2 August 2021 | Volume 11 | Article 627811 Escuin et al. miRNAs and BC Locoregional Metastases histology (ductal, lobular, special types), presence of associated using the FastQC tool (http://www.bioinformatics.babraham.ac. ductal or lobular carcinoma in situ, presence of vascular and uk/projects/fastqc/). lymphatic invasion, tumor infiltrating lymphocytes, type of invasion (expansive/infiltrating), tumor multifocality, tumor Genome Annotation and Quantification necrosis; proliferation of non-tumoral tissue (ductal of miRNAs hyperplasia, atypical ductal/lobular hyperplasia). Genome annotation was performed using the Exiqon/Xplore RNA pipeline. Following sequencing, Cutadapt (1.9.1) (29) was used to trimmed adaptor sequences. A quality check (QC) was Blood Processing and Isolation of Plasma performed to ensure Q-scores >30 (>99.9% correct) of our data Human plasma samples were collected prospectively from early (30). Reads with correct length were analyzed for the presence of breast cancer patients who have not received any previous UMIs using Cutadapt (1.9.1) and then collapsed by UMIs into treatment. Peripheral blood was withdrawn before surgery. FASTQ files. This approach eliminates library amplification bias Approximately, 10-15ml of peripheral blood was collected for and allows for true identification of the miRNAs. Bowtie2 plasma processing in EDTA tubes. Plasma tubes were processed software (2.2.6) was used for mapping the reads. The mapping within 2 hours of collection and spun at 1200xg for 10 minutes. criterion for aligning reads to spike-ins, abundant sequences and Plasma was aliquot in 1.5 ml tubes and stored at -80C until miRBase_20 was for reads to have perfect match to the reference further processing. All plasma samples used in this study were sequences. To map the genome, one mismatch was allowed in inspected for absence of hemolysis as previously described (27). the first 32 bases of the read. Small insertions and deletions Briefly, the hemolysis score (HS) was determined by ultraviolet- (INDELs) were not allowed. The resulting sequences were visible (UV-Vis) absorbance measurements using a NanoDrop® annotated using the human assembly GRCh37 and the 2000 Spectrophotometer (Thermo Scientific, Barrington, IL, miRBase_20 database. IsomiR analysis was performed USA). Measurements were performed by applying 2 µl of individually for each sample based on the occurrence of count plasma on the micro-volume pedestal after centrifugation at variants for each detected miRNA. Read variants were merged 1000 × g for 5 min at 4°C and using saline (PBS) as a blank. In onto a single count file with a consistent nomenclature across addition, monitoring of hemolysis was conducted by qPCR for samples. Only isomiRs present at a level of 5% of total reads for a all samples by comparing the level of a highly expressed miRNA specific miRNA were retained. Transcripts per million (TPM) in red blood cells (hsa-miR-451a) with a miRNA unaffected by was used as a normalization procedure to correct for differences hemolysis (hsa-miR-23a-3p) as previously described (28). in sequencing depth and to quantified each RNA species. Samples with a DCt > 7 were discarded for further analyses. Differential Expression Analysis RNA Isolation NGS Library Preparation Differential expression analysis was performed using the EdgeR and Next Generation Sequencing statistical software package (Bioconductor, http://www. RNA was isolated from 300ml of plasma samples with the bioconductor.org/). The analysis was performed using the miRNeasy serum/plasma advanced kit (Qiagen, Cat No/ID: trimmed mean of M-values normalization method (TMM) (31), 217204) according to the manufacturer’s instructions. A range based on the log-fold and absolute gene-wise changes in expression of spike-ins were added to the plasma samples prior to RNA between samples. The TMM normalization compensates for isolation. A quality check was performed by qPCR previous to sample specific effects caused by the variation in library sequencing the samples. Sixteen samples were selected to size/sequencing depth between samples and also compensates for p erf orm NGS, i nclu din g 12 pos itive SL Ns (n = 6 under- or over-sampling effects by trimming and scaling factors macrometastasis and n= 6 micrometastasis) and 4 negative that minimize log fold changes between samples across the SLNs. Five ml of total RNA was used to construct the NGS majority of the miRNAs. The isomiR analysis was done using libraries using the QIAseqmiRNA Library Kit (Qiagen, Cat. No: Exiqon in-house scripts (exq_ngs_mircount). Predicted miRNAs 331505). Briefly, after ligation of 3’ and 5’ adapters and Unique analysis was performed based on the read count distribution using Molecular Identifier (UMIs) to miRNAs, complementary DNA the exiqon_ngs_mirpred in house script and the secondary libraries were constructed by reverse transcription followed by 22 structure prediction according to the miRPara classification score cycles of PCR amplification and cDNA cleaned up using QMN (32). Volcano plots were constructed using R programming (33) by beads. A library preparation quality check was performed using plotting the p value (-log10) on the y-axis and the expression fold either Bioanalyzer 2100 (Agilent) or TapeStation 4200 (Agilent). change between the two experimental groups on the x-axis. Based on quality of the inserts and the concentration measurements, libraries were pooled in equimolar ratios and Principal Component Analysis and Heat quantified using the qPCR ExiSEQ LNA™ Quant kit (Exiqon). Map and Unsupervised Clustering The library pools were sequenced with a NextSeq500 platform Principal component analysis (PCA) was performed using R (Illumina) using sequence runs of 75nt single-end reads with an programming and TMM-normalized values as input. The same average number of 10 million reads/sample. Raw data was input was used to generate a heatmap and unsupervised demultiplexed and FASTQ files were generated using the hierarchical clustering by samples and gene expression profile bcl2fastq 2.18.0.12 software (Illumina) and files were checked with R scripts (33). We selected the top 50 miRNAs with the Frontiers in Oncology | www.frontiersin.org 3 August 2021 | Volume 11 | Article 627811 Escuin et al. miRNAs and BC Locoregional Metastases largest coefficient of variation (% CV) across all samples to SLN-positive tumors (76%), of which 6 were OSNA-diagnosed as obtain a cluster of samples. The data was normalized to TMM micrometastasis (38%) and 6 as macrometastasis (38%). and converted to log2 scale. All samples passed the post-sequencing QC, which confirmed that the average read quality and base quality had a Q-score > 30 Gene Ontology Enrichment Analysis (99.9% correct) (30) and the expected read length distribution for Gene Ontology (GO) analyses (34, 35) were done with R TopGO miRNAs (Supplementary Figure 1). All samples were sequenced in package with experimentally verified targets of significantly one excellent run with a median 27.2 million read number. Following differentially expressed miRNAs as input. Two different sequencing and trimming, all reads containing identical insert statistical tests were used and compared. First, a standard sequence and UMI sequence (insert-UMI pair) were collapsed into Fisher’s test was used to investigate enrichment of terms a single read and passed into the analysis pipeline. This allowed for between groups. Second, the Elim method (36) was used to true quantification of the miRNAs by eliminating library incorporate the topology of the GO network and to compensate amplification bias and a better representation of the RNA for local dependencies between GO that could mask significant molecules in the sample. We obtained an average of 1.8 million GO terms. Comparisons from these two methods were used to collapsed reads for each sample and good miRNA mapping reads highlight relevant GO terms. with a very dominant miRNA peak in most of the samples, indicating a good sample/data quality (Supplementary Table 1 and Quantitative Real-Time RT–PCR Analysis Supplementary Figure 2). Overall, we obtained an average Quantitative real-time RT–PCR analysis was done with an ABI genome mapping rate of 46.2% (Figure 1A), which are values well Prism 7500 Sequence Detection System using the miRCURY within the range for plasma samples. After mapping and counting to LNA™ Universal RT cDNA Synthesis Kit (Exiqon). The cDNA relevant entries in mirbase_20, the number of known miRNAs was was diluted 50X and assayed in 10 µl PCR reactions according to calculated using TPM to measure expression. We found comparable the protocol for the miRCURY LNA™ Universal RT microRNA numbers of identified miRNAs using either TPM > 1 (182 miRNAs) PCR System (Exiqon A/S); each microRNA was assayed twice by or TPM > 10 (125 miRNAs) (Figure 1B). We did not identify any qPCR on the plasma Focus microRNA PCR panel. A no- sequences identical to those of known miRNAs in miRBase_20 for template control (NTC) of water was purified with the samples other organisms. However, we were able to predict 80 miRNAs based and profiled like the samples. Analysis of the data was performed on the structural properties of the genome in the indicated locations using the relative miRNA expression according to the resembling those of known miRNAs (Supplementary Table 2). comparative Ct (DDCt) method using negative metastatic Next, we investigated whether the patients were assigned into samples as reference. We used the geNorm (37) or the biological groups based on their miRNA expression. We performed Normfinder algorithm (38) to select the best combination of an unsupervised two-way hierarchical clustering of miRNAs and two reference genes. Data from multiples plates were normalized using UniSp3 spike-in (Exiqon) as interplate calibrators. TABLE 1 | Basic patient and tumor characteristics. Statistics Variable N (%) Differentially expressed miRNAs from RNA-sequencing data were Patients Plasma 16 (100) detected by an exact test based on conditional maximum Age, years Mean +SD 63 ± 13 likelihood (CML) included in the R Bioconductor package Median (range) 63 (46 - 89) edgeR (39). P values from RNA-sequencing were corrected Tumor status T1 10 (62) (q-values) for multiple testing using the Benjamini-Hochberg T2 2 (13) procedure (40). A false discovery rate (FDR) q < 0.05 was T3 4 (24) Node status Negative 4 (25) considered significant. In all group comparisons missing Micrometastasis 6 (38) expression values were treated as zero. Differences in total Macrometastasis 6 (38) numbers of miRNAs between groups were analyzed by two- Axillary Lymph Node Status Negative 3 (50) sided parametric t-tests. For analysis of clinicopathological Positive 3 (50) Tumor grade I 5 (31) parameters, quantitative variables between groups were II 9 (56) compared using the Student’s T-test and qualitative variables III 2 (13) were compared using the X2 or Fisher exact tests. A two-sided Estrogen receptor Negative 0 (0) p-value ≤ 0.05 was considered significant. Positive 16 (100) Progesterone receptor Negative 1 (6) Positive 15 (94) HER2 status Negative 16 (100) RESULTS Positive 0 (0) Ki67 < 20% 15 (94) A total of 25 patients were included in this study. However, only > 20% 1 (6) Surgery Mastectomy 5 (31) samples from 16 patients passed the pre-RNA-sequencing quality Lumpectomy 11 (69) check (QC). The main clinicopathological characteristics of the Lymphovascular invasion Negative 12 (75) patients are described in Table 1. A total of 12 (75%) patients had Positive 4 (25) Frontiers in Oncology | www.frontiersin.org 4 August 2021 | Volume 11 | Article 627811 Escuin et al. miRNAs and BC Locoregional Metastases A B FIGURE 1 | Summary of the mapping results for all samples. (A) Percentage of sequencing reads for each sample. Reads are classified as miRNAs, small RNAs, genome-mapped, outmapped, high abundance (e.g. rRNA, polyA,mtRNA) and unmapped reads. (B) Number of identified known miRNAs with transcripts per million (TPM) normalized numbers of counts >1 (blue bars) or >10 (red bars). samples using the 50 miRNAs with the largest coefficient of multiple testing (q < 0.05) (Table 2). Fourteen miRNAs were down- variation based on TMM counts (Figure 2A). Our results show regulated (miR-339-5p, miR-133a-3p, miR-326, miR-331-3p, miR- that samples did not cluster according to the SLN outcome of the 369-3p, miR-328-3p, miR-26a-3p, miR-139-3p, miR-493-3p, miR- patients, suggesting that other clinicopathological factors are 664a-5p, miR-146a-5p, miR-323b-3pmiR-1307-3p and miR-423-3p) responsible for the variation on the samples. We obtained similar and 2 miRNAs were up-regulated (miR-101-3p and miR-144-3p) results using a principal component analysis. Interestingly, the 2 (Figure 3A). Next, we analyzed the data based on SLN metastasis samples showing the greater variability (p18 and p62) corresponded status subgroups. When we compared patients with macrometastasis to those patients whose tumors had a mixed pathological vs. negative SLNs, we found 42 miRNAs differentially expressed, but component (IDC and ILC) (Figure 2B). only miR-339-5p remained significant after FDR adjustment Despite the unsupervised analysis did not group our samples (p <0.0001, q = 0.0413) (Figure 3B and Table 2). Similar results according to the metastasis status of the patients, we identified were obtained when we compared micrometastasis and negative differentially expressed miRNAs between groups based on the SLN SLNs, which yield 66 miRNAs differentially expressed, but only miR- outcome of the patients. First, we analyzed samples according to the 376c-3p (p = 0.0001, q = 0.046), miR-326 (p = 0.0003, q = 0.049) and positive (n=12) or negative (n=4) SLN metastasis status. We found miR-323b-3p (p = 0.0004, q = 0.049) passed the FDR (Figure 3C and 73 miRNAs with a significant differential expression (p < 0.05). Table 2). Interestingly, we did not find any significantly differentially However, only 16 miRNAs remained significant after correcting for expressed circulating miRNAs between patients with Frontiers in Oncology | www.frontiersin.org 5 August 2021 | Volume 11 | Article 627811 Escuin et al. miRNAs and BC Locoregional Metastases A B FIGURE 2 | Class discovery associated with the SLN metastatic status. The analysis was performed using the 50 miRNAs with the largest coefficient of variation based on trimmed mean of M-values (TMM counts). (A) Heat map and unsupervised hierarchical clustering. Each row represents one miRNA and each column represents one sample. The color represents the relative expression level of a miRNA across all samples. The color scale shows the expression level above (red) or below (green) the mean. (B) Principal Component Analysis (PCA) of samples according to the SLN metastatic status of patients. macrometastasis or micrometastasis SLNs (Figure 3D and Table 2). specific qPCR assays. The down-regulation of 9 out of 14 miRNAs In addition, we used the 16 significant differentially expressed was confirmed in patients with positive SLNs, but we could not miRNAs in patients with positive SLNs to build a heatmap and validate the up-regulation of miR-101-3p and miR-144-3p. hierarchical clustering. Our results show that these miRNAs clearly Furthermore, the degree of down-regulation was higher for those separated patients with negative and positive SLNs (Figure 3E). patients that had additional metastases in their axillar lymph nodes Further validation on the same clinical samples was performed by (Supplementary Figure 3). Frontiers in Oncology | www.frontiersin.org 6 August 2021 | Volume 11 | Article 627811 Escuin et al. miRNAs and BC Locoregional Metastases TABLE 2 | Differentially expressed miRNAs. Names Sequence (5’ – 3’) TMM TMM logFC p value q value Positive Negative hsa-miR-339-5p TCCCTGTCCTCCAGGAGCTCACG 37.4 127.6 -1.8 < 0.0001 0.007 hsa-miR-133a-3p TTTGGTCCCCTTCAACCAGCTG 6.1 28.4 -2.0 < 0.0001 0.008 hsa-miR-326 CCTCTGGGCCCTTCCTCCAG 10.6 50.0 -2.2 < 0.0001 0.008 hsa-miR-331-3p GCCCCTGGGCCTATCCTAGAA 0.9 10.3 -2.8 0.0001 0.011 hsa-miR-369-3p AATAATACATGGTTGATCTTT 7.9 26.0 -1.7 0.0005 0.031 hsa-miR-328-3p CTGGCCCTCTCTGCCCTTCCGT 134.1 350.9 -1.4 0.0005 0.031 hsa-miR-26a-1-3p CCTATTCTTGGTTACTTGCACG 2.1 10.9 -2.5 0.0007 0.034 hsa-miR-139-3p TGGAGACGCGGCCCTGTTGGAGT 30.2 79.8 -1.4 0.0008 0.034 hsa-miR-493-3p TGAAGGTCTACTGTGTGCCAGG 1.9 11.5 -2.1 0.0010 0.034 hsa-miR-664a-5p ACTGGCTAGGGAAAATGATTGGAT 49.7 108.1 -1.1 0.0010 0.034 hsa-miR-101-3p TACAGTACTGTGATAACTGAA 6070.1 3215.4 0.9 0.0011 0.034 hsa-miR-146a-5p TGAGAACTGAATTCCATGGGTT 2960.5 6266.3 -1.1 0.0012 0.034 hsa-miR-144-3p TACAGTATAGATGATGTACT 511.3 295.1 0.8 0.0013 0.035 hsa-miR-323b-3p CCCAATACACGGTCGACCTCTT 10.5 29.9 -1.4 0.0016 0.040 hsa-miR-1307-3p ACTCGGCGTGGCGTCGGTCGTG 150.3 337.4 -1.2 0.0017 0.040 hsa-miR-423-3p AGCTCGGTCTGAGGCCCCTCAGT 353.7 649.1 -0.9 0.0023 0.050 hsa-miR-376c-3p AACATAGAGGAAATTCCACGT 3.8 14.0 -1.8 0.0028 0.056 hsa-miR-1 TGGAATGTAAAGAAGTATGTAT 69.5 168.7 -1.3 0.0037 0.071 hsa-miR-1908-5p CGGCGGGGACGGCGATTGGTC 23.1 59.1 -1.3 0.0042 0.073 hsa-miR-744-5p TGCGGGGCTAGGGCTAACAGCA 138.5 298.1 -1.1 0.0042 0.073 hsa-miR-584-5p TTATGGTTTGCCTGGGACTGAG 419.9 835.4 -1.0 0.0048 0.078 hsa-miR-6721-5p TGGGCAGGGGCTTATTGTAGGAG 2.4 9.3 -1.9 0.0055 0.083 hsa-miR-432-5p TCTTGGAGTAGGTCATTGGGTGG 130.1 432.4 -1.7 0.0055 0.083 hsa-miR-28-3p CACTAGATTGTGAGCTCCTGGA 72.1 150.4 -1.0 0.0058 0.084 hsa-miR-29b-3p TAGCACCATTTGAAATCAGTGTT 292 154 0.93 0.0068 0.094 Macrometastasis Negative hsa-miR-339-5p TCCCTGTCCTCCAGGAGCTCACG 38.9 138.0 -1.8 0.0001 0.041 Micrometastasis Negative hsa-miR-376c-3p AACATAGAGGAAATTCCACGT 1.5 13.7 -3.1 0.0001 0.046 hsa-miR-326 CCTCTGGGCCCTTCCTCCAG 9.8 49.2 -2.3 0.0003 0.049 hsa-miR-323b-3p CCCAATACACGGTCGACCTCTT 7.1 29.5 -2.0 0.0004 0.049 Macrometastasis Micrometastasis hsa-miR-122-5p TGGAGTGTGACAATGGTGTTTG 8948.2 76103.3 -3.1 0.0002 0.062 hsa-miR-125b-2-3p TCACAAGTCAGGCTCTTGGGAC 0.3 8.9 -3.6 0.0006 0.090 Data shows the 25 most significant differentially expressed miRNAs according to the metastatic status of patients. The list includes the average trimmed mean of M-values (TMM) values, logarithmic fold change (logFC), raw p values and Benjamini-Hochberg FDR corrected q values. Next, we sought to understand how our data is related to micrometastasis in their SLNs. As for the cellular component biological functions by performing a gene ontology (GO) (CC), the target miRNAs were significantly located in vesicle and analysis. Selecting Homo sapiens as the background of listed membrane fractions (p < 0.01). Moreover, differentially target genes, we obtained the GO term annotations and KEGG expressed miRNAs were enriched in molecular function (MF) pathway analysis through the functional annotation summaries. terms associated with transcription factors, G protein-related The results are summarized in Figure 4 and Table 3. The top 50 coupled peptide receptor activity, receptor regulator activity and biological process GO terms (p < 0.05) associated with microtubule motor activity (p < 0.01) (Table 3). differentially expressed circulating miRNAs in patients with Our series include 16 patients with early breast cancer and we positive SLNs compared to the reference background (negative reported recurrence in 3 (19%) patients. The median follow-up SLNs samples) are shown in Figure 4A. Our data shows that time was 5.2 years (range 2.2 - 6.4 years). Two patients had differentially expressed miRNAs associated with biological secondary tumors in the colon and 1 patient in the liver. At last processes (BP) markedly focused on epigenetic gene expression follow-up, two patients with recurrences in colon were reported regulation, the epithelial-mesenchymal transition (EMT), alive with disease and we reported 3 deaths in patients due to transcription, cell motility and proliferation processes (p < complications related to the disease. We investigated whether the 0.01) (Figure 4B, Table 3 and Supplementary Data). For differential expressed miRNAs correlated with the patient’s instance, we found that positive regulation of mesenchymal clinico-pathological parameters. The expression of miR-326, cell proliferation term (GO: 0002053) was significantly miR-26a-1-3p, miR-139-3p, miR-101-3p, miR-146a-5p and enriched (p < 0.0028) as well as positive regulation of the miR-144-3p was significantly lower associated in younger histone H3-H4 methylation term (GO:0051571) (p < 0.0017). patients (< 60 years, p < 0.05), the expression of miR-328-3p These two GO terms remained significant even when patients and miR-144-3p was associated with further metastases in the with positive SLNs were sub-classified as having macro- or aLNs (p < 0.05), miR-26a-1-3p,miR-144-3p and miR-323-3p Frontiers in Oncology | www.frontiersin.org 7 August 2021 | Volume 11 | Article 627811 Escuin et al. miRNAs and BC Locoregional Metastases A B C D E FIGURE 3 | Differentially expressed miRNAs according to the SLN status. (A–D) The volcano plots show differentially expressed miRNAs in plasma samples according to the patients’ locoregional metastatic status as indicated. Only significant miRNAs with corrected q values < 0.05 are shown in the plots (red dots). The data show the relationship between non-adjusted p values (y-axis) and the fold change (x-axis) between the experimental groups. (E) Heat map and hierarchical clustering analyzed by samples and miRNAs. The analysis was performed using the 15 miRNAs differentially expressed between patients with positive and negative SLNs. Each row represents one miRNA and each column represents one sample. The color scale shows the expression level above (red) or below (green) the mean. were associated with tumor stage (p < 0.05), miR-664a-5p and DISCUSSION miR-323b-3p showed a non-significant association with tumor status (p = 0.077 and p = 0.069, respectively) and miR-26a-1-3p In recent years, miRNAs have emerged as important regulators showed a non-significant correlation with the recurrence status of the various steps of the metastatic process (41). Currently, (p = 0.067). We did not find any other significant association lymph node affection remains the most important prognosis with other parameters. Due to the low number of events, we were factor in breast cancer (2) and the presence of metastasis in the unable to perform any survival analysis in our cohort of patients. SLNs is still currently the recommended procedure for axillary Frontiers in Oncology | www.frontiersin.org 8 August 2021 | Volume 11 | Article 627811 Escuin et al. miRNAs and BC Locoregional Metastases A B FIGURE 4 | Gene ontology (GO) enrichment analysis for the significant biological processes associated with positive SLNs. (A) Dot plot graph shows the 50 most significant biological process GO terms (y-axis) and the ratio between the number of expressed miRNAs associated to the GO term and the number of significantly differentially expressed genes associated to the GO term (x-axis). The color of the nodes indicates the p-value and the size of the nodes the number of miRNAs associated with a specific GO term. (B) Neural network shows the GO terms for the biological processes associated with patients with positive SLNs. Frontiers in Oncology | www.frontiersin.org 9 August 2021 | Volume 11 | Article 627811 Escuin et al. miRNAs and BC Locoregional Metastases TABLE 3 | Gene Ontology (GO) analysis. GO ID GO Term Counts p value BIOLOGICAL PROCESS GO:0051571 positive regulation of histone H3-K4 methylation 9/15 0.0017 GO:0042462 eye photoreceptor cell development 12/26 0.0017 GO:0070555 response to interleukin-1 17/40 0.0021 GO:0002053 positive regulation of mesenchymal cell proliferation 25/62 0.0028 GO:0071320 cellular response to cAMP 12/23 0.0034 GO:0002407 dendritic cell chemotaxis 4/5 0.0035 GO:0009629 response to gravity 6/8 0.0042 GO:0007097 nuclear migration 6/9 0.006 GO:0051573 negative regulation of histone H3-K9 methylation 6/8 0.0064 GO:0034446 substrate adhesion-dependent cell spreading 8/16 0.0067 CELLULAR COMPONENT GO:0031091 platelet alpha granule 18/44 0.0065 GO:0031983 vesicle lumen 17/46 0.0137 GO:0060205 cytoplasmic membrane-bounded vesicle lumen 17/46 0.0137 GO:0031093 platelet alpha granule lumen 16/41 0.0171 GO:0034774 secretory granule lumen 16/41 0.0171 GO:0044306 neuron projection terminus 9/20 0.0216 GO:1902495 transmembrane transporter complex 14/31 0.0218 GO:1990351 transporter complex 14/31 0.0218 GO:0015030 Cajal body 8/16 0.0239 GO:0034704 calcium channel complex 7/12 0.0245 MOLECULAR FUNCTION GO:0008528 G-protein coupled peptide receptor activity 9/15 0.0013 GO:0001618 virus receptor activity 6/13 0.0018 GO:0030955 potassium ion binding 5/6 0.0064 GO:0005161 platelet-derived growth factor receptor binding 17/39 0.0078 GO:0008798 beta-aspartyl-peptidase activity 4/7 0.0097 GO:0030545 receptor regulator activity 15/36 0.0097 GO:0003777 microtubule motor activity 11/20 0.0115 GO:0046625 sphingolipid binding 5/6 0.0128 GO:0008307 structural constituent of muscle 6/8 0.0131 GO:0017022 myosin binding 6/11 0.0144 Gene set enrichment analysis using GO categories (biological process, cellular component, molecular function) was applied to extract biological meaning from the identified differentially expressed transcripts and predicted mRNA targets. The top 10 GO categories associated with differentially expressed circulating miRNAs in patients with positive SLNs are shown. Counts refers to the ratio between the number of enriched differentially expressed miRNAs and the total number of miRNAs assigned to these terms. P values were calculated with a combination of the Elim method and the Fisher's exact method. GO terms with p values < 0.05 were considered enriched. staging of early breast cancer. The accurate evaluation of patients that showed a 1.9- and 1.7-fold change up-regulation, with involved SLN determines further axillary lymph node respectively. However, we could not confirm the up-regulation dissection (ALND), the golden standard procedure for invasive of these 2 miRNAs and these results agree with the discrepancies breast cancer. However, ALND has been questioned in recent on the direction of the dysregulation for both miRNAs. For years because of inherent morbidity following the procedure instance, dysregulation miR-101-3p has been reported in several without directly contributing to survival. In this study, we sought malignancies, including breast cancer (42, 43). While some reports to gain a better understanding of the role of miRNAs in the indicate up-regulation of miR-101-3p, others indicated the metastatic process and whether specific expression patterns of opposite (42). This due to the fact that mature miR-101-3p miRNAs could predict SLN metastatic status in patients with originates from two different precursors located at different early breast cancer. We performed a proof-of-principle study in chromosomes. One precursor may be processed to 1 or 2 plasma samples from 16 breast cancer patients with known SLN miRNAs and thus, the mature and precursor miRNA levels metastasis status. Importantly, plasma samples were collected might not correlate, and this therefore will influence the clinical prior to any treatment, thus the results using RNA-sequencing interpretation. The same study looked at putative miR-101-3p reflect the basal miRNA expression prior to any therapeutic target genes were analyzed and the most predominant functions intervention in these patients. Our results show a good quality were transcription, metabolism, biosynthesis, proliferation, and sequencing data with mapping rates to miRNAs and comparable transcription factor binding. This result indicated that candidate miRNA discovery across samples. Thus we are confident in the genes have a definitive impact on the pathogenesis of BC (42). accuracy of the reported results. Similar conflicting data has been reported for miR-144-3p. In Our data shows that 16 miRNAs were significantly several human cancers, the expression of miR-144-3p has been differentially expressed in plasma samples from SLN-positive shown decreased (44), but in animal models repression of miR- patients. Overall, we found a general down-regulation of 144 significantly decreased cell proliferation, clonogenicity, miRNAs, with the exception of miR-101-3p and miR-144-3p migration and tumor formation in nude mice (45). Interestingly, Frontiers in Oncology | www.frontiersin.org 10 August 2021 | Volume 11 | Article 627811 Escuin et al. miRNAs and BC Locoregional Metastases one report has shown that up-regulation of miR-144-3p was Among the other circulating miRNAs that were down- associated with families at high-risk for breast cancer (46). regulated in our study, the expression of miR-133a-3p has been These data suggest that the role of miR-144-3p might differ by reported to be down-regulated in paired breast cancer tumor and cancer type and tumor microenvironment. serum samples (64), suggesting the tumor origin of miR-133a-3p. Of those miRNAs down-regulated, miR-339-5p showed a 3.5- In contrast, miR-133a-3p has been found elevated in plasma fold inhibition in patients with positive SLN metastasis. The samples from early-stage BC patients compared to healthy expression of miR-339-5p remained significant when the analysis donors (65, 66) and similar results have been reported for was performed in the subgroup of macrometastatic SLNs and we circulating miR-1307-3p (67). On the other hand, down- observed a non-significant trend towards significance for the regulation of miR-376c-3p (68) and miR-376c-3p have been subgroup of micrometastatic SLNs (q = 0.071). Our results agree linked to breast cancer recurrence (24). MiR-326 has been with previous reports showing that reduced miR-339-5p reported to target B7-H3 in breast cancer, an immunoregulatory expression in breast cancer is associated with increased protein that is overexpressed in several cancers and is often metastasis to lymph nodes (47, 48), high clinical stages and associated with metastasis and poor prognosis (69). worse clinical outcome (47). A similar association with positive Furthermore, its expression has been shown to inhibit LN has been reported in NSCLC patients (49). In addition, miR- tumorigenesis through direct targeting of Nin one binding 339-5p expression is down-regulated in several human cancers protein (NOB1) and the MAPK pathway in glioma cells (52). including NSCLC (49), ovarian carcinoma (50), hepatocellular A main limitation of our study is the small number of samples carcinoma (51), gliomas (52), colorectal cancer (53), analyzed since it was designed as a proof-of-principle study to assess osteosarcoma (54) and breast cancer (48). Mir-339-5p acts as a the feasibility of using circulating miRNAs as potential surrogates of tumor suppressor gene and its expression is required to inhibit the lymph node metastatic status in breast cancer. Therefore, our cell migration and invasion in breast cancer cells (47) in a results are preliminary and must be interpret with caution. mechanism that involves at least the B-cell lymphoma 6 Nonetheless, our data shows several circulating miRNAs that are (BCL6) protein. The authors showed that forced expression of significantly differentially expressed in relation to the SLN metastatic BCL6 results in increased proliferation, anchorage-independent status of the patients. Moreover, we report an overall down- growth, migration, invasion and survival of breast cancer cell regulation of these miRNAs, which in most cases have been lines, whereas knockdown of BCL6 expression reduced these reported to be direct targets of proteins that promote metastasis. oncogenic properties of breast cancer cells (55). Interestingly, Further studies in a larger cohort of patients are warranted to validate miR-339-5p has been shown to inhibit migration and invasion these results and to unveil the molecular mechanisms of the miRNAs by targeting BCL6 in breast cancer (56), ovarian cancer cell lines described here and the various steps of the locoregional metastasis. (50) and in NSCLC (57). In addition, miR-339-5p down- regulation in NSCLC inhibits metastasis of NSCLC by regulating the epithelial-to-mesenchymal (EMT) transition via DATA AVAILABILITY STATEMENT BCL6 (57). A recent report has shown that miR-339-5p regulates EMT through regulation of TGF-b (58) in osteosarcoma (54). The original contributions presented in the study are publicly The EMT and the TGF-b pathways are two of the most available at the Sequence Research Archive under ID important mechanisms underlying the metastatic ability of PRJNA669408 and are available for download here http://www. cancer cells (59, 60). We have previously shown the ncbi.nlm.nih.gov/bioproject/669408. importance of the EMT in breast cancer (61) and here, we show that GO term analysis based on the DE miRNAs showed a significant association with the biological process “positive ETHICS STATEMENT regulation of mesenchymal cell proliferation” (GO:0002053). Another pathway enriched was the positive regulation of H3 This study was conducted according to the Declaration of K4 methylation (GO:0051571), a mark that on a genome-wide Helsinki principles, with approval from the Clinical Research scale is broadly associated with transcriptional regulation, and Ethics Committee at “Institut d’Investigació ’ Biomé dica Sant “negative regulation of H3K9 methylation” (GO:0051573). H3K9 Pau” (IIB Sant Pau). Written informed consent was obtained methylation has been associated with the EMT through from all patients under institutional review board-approved interactions of KDM1A (a H3K9 demethylase) with the protocols. All methods were performed in accordance with the members of the SNAI1 family of zinc finger transcription relevant guidelines and regulations. factors, including SNAI1 (SNAIL) and SNAI2 (SLUG). The expression of SNAI1 and E-cadherin is a hallmark of carcinoma development and metastasis (62). Our data suggest AUTHOR CONTRIBUTIONS that that MiR-101 could be involved in the regulation of these pathways, as it has been shown to directly target the histone DE, LL-V, and AB were involved in the conceptualization of the methyltransferase enhancer of zeste homologue 2 (EZH2), which project. DE, LL-V, OB, JM, IP, AM, CA, BGV, TR and EL were could promote tumor proliferation and invasion (63). involved in resources, investigation and methodology. DE, LL-V, EL, Frontiers in Oncology | www.frontiersin.org 11 August 2021 | Volume 11 | Article 627811 Escuin et al. miRNAs and BC Locoregional Metastases and AB were involved in analysis and interpretation of data. All added to the samples prior to RNA isolation. We observed an excellent correlation of counts corresponding to the spike-ins between the samples. authors contributed to the article and approved the submitted version. Supplementary Figure 2 | Summary of the mapping results for the samples. (A) Total number of reads for each sample sequenced. (B) Read length distribution for each class of RNAs. The peak around 18-23 nucleotides (red) corresponds to FUNDING miRNAs. Instituto de Salud Carlos III and European Regional Supplementary Figure 3 | RNA-sequencing validation. Relative gene expression Development Fund (PI13/00110 and PI19/00362 to D.E.); was performed according to the comparative ddCt (DDCt) method using negative Centro de Investigació n Biomé dica en Red Cá ncer (CB 16/12/ metastatic samples as reference. The geNorm or the Normfinder algorithm were used to select the best combination of two reference genes. Data from multiples 00471 to AB); Red Temá tica de Investigació n Cooperativa de plates were normalized using UniSp3 spike-in as interplate calibrators. Each Cá ncer (RTICC; RD12/0036/0076 to all authors. microRNA was assayed twice by qPCR on the Serum/plasma Focus microRNA PCR panel. (A) Data shows the comparison between patients with negative and positive SLNs. (B) Patients with positive SLNs were divided according to the presence or absence of further axillary lymph nodes (aLNs). Statistical analysis was ACKNOWLEDGMENTS performed using unpaired Student’s t-test. All comparisons shown are statistically significant (p < 0.05). We thank the patients and their caregivers for their participation. Supplementary Table 1 | Mapping and yield results of the sequencing data. Reads are classified into the following classes: miRNA, maps to mirBase_20; small RNA, maps to sRNA database compiled by Exiqon; Predicted miRNA, predicted miRNAs using a prediction software; Genome, reads aligning to the SUPPLEMENTARY MATERIAL reference genome, but not to miRNAs or sRNAs; Outmapped, maps corresponding to poly A and poly homopolymers as well as abundant rRNA and mtRNA; The Supplementary Material for this article can be found online Unmapped, reads not aligning to reference genome. Median values for each type at: https://www.frontiersin.org/articles/10.3389/fonc.2021. of reads are shown. 627811/full#supplementary-material Supplementary Table 2 | Predicted miRNA results of the sequencing data. Supplementary Figure 1 | Data quality checking. (A) Average read quality of the Predicted miRNAs for each sample sorted by total read count. The sequences next generation sequencing (NGS) data. The average read Q-score is plotted on the reported do not match any known miRNA in the miRBase_20. The location is x-axis and the number of reads on the y-axis. A Q-score above 30 is considered formatted as "chromosome: start-stop (strand)". Counts describes the number of high quality data. (B) Base quality of the NGS data. The position in the read is reads which fall onto the location of the predicted miRNA (sum of all samples). Note plotted on the x-axis and the Q-score is plotted on the y-axis. The red line is the that the name and numbering assigned here are project-specific and cannot be median Q-score and the dark blue is the mean value Q-score. A Q-score above 30 compared with the nomenclature found in the miRBase or other databases. (>99.9% correct) is considered high quality data. (C) Read length distribution and adaptor trimming. miRNAs will appear as a peak around 18-23 nucleotides. (D) Supplementary Table 3 | Contributing miRNAs to GO terms for the categories Radar plot showing relative spike-in signal for the samples. A range of spike-ins was biological process (BP), molecular function (MF) and cellular compartment (CC). 9. Peg V, Espinosa-Bravo M, Vieites B, Vilardell F, Antunez JR, de Salas MS, et al. REFERENCES Intraoperative Molecular Analysis of Total Tumor Load in Sentinel Lymph Node: 1. DeSantis CE, Ma J, Gaudet MM, Newman LA, Miller KD, Goding Sauer A, A New Predictor of Axillary Status in Early Breast Cancer Patients. Breast Cancer et al. Breast Cancer Statistics, 2019. CA: Cancer J Clin (2019) 69(6):438–51. Res Treat (2013) 139(1):87–93. doi: 10.1007/s10549-013-2524-z doi: 10.3322/caac.21583 10. Giuliano AE, McCall L, Beitsch P, Whitworth PW, Blumencranz P, Leitch 2. Shek LL, Godolphin W. Model for Breast Cancer Survival: Relative Prognostic AM, et al. Locoregional Recurrence After Sentinel Lymph Node Dissection Roles of Axillary Nodal Status, TNM Stage, Estrogen Receptor Concentration, With or Without Axillary Dissection in Patients With Sentinel Lymph Node and Tumor Necrosis. Cancer Res (1988) 48(19):5565–9. Metastases: The American College of Surgeons Oncology Group Z0011 3. Wildiers H, Van Calster B, van de Poll-Franse LV, Hendrickx W, Roislien J, Randomized Trial. Ann Surg (2010) 252(3):426–32. doi: 10.1097/ Smeets A, et al. Relationship Between Age and Axillary Lymph Node SLA.0b013e3181f08f32 Involvement in Women With Breast Cancer. J Clin Oncol: Off J Am Soc 11. Jagsi R, Chadha M, Moni J, Ballman K, Laurie F, Buchholz TA, et al. Clin Oncol (2009) 27(18):2931–7. doi: 10.1200/JCO.2008.16.7619 Radiation Field Design in the ACOSOG Z0011 (Alliance) Trial. J Clin 4. Rivadeneira DE, Simmons RM, Christos PJ, Hanna K, Daly JM, Osborne MP. Oncol: Off J Am Soc Clin Oncol (2014) 32(32):3600–6. doi: 10.1200/ Predictive Factors Associated With Axillary Lymph Node Metastases in T1a JCO.2014.56.5838 and T1b Breast Carcinomas: Analysis in More Than 900 Patients. J Am Coll 12. Giuliano AE, Ballman KV, McCall L, Beitsch PD, Brennan MB, Kelemen PR, Surge (2000) 191(1):1–6. doi: 10.1016/S1072-7515(00)00310-0 et al. Effect of Axillary Dissection vs No Axillary Dissection on 10-Year 5. Weigelt B, Peterse JL, van ‘t Veer LJ. Breast Cancer Metastasis: Markers and Overall Survival Among Women With Invasive Breast Cancer and Sentinel Models. Nat Rev Cancer (2005) 5(8):591–602. doi: 10.1038/nrc1670 Node Metastasis: The ACOSOG Z0011 (Alliance) Randomized Clinical Trial. 6. Patil N, Allgayer H, Leupold JH. MicroRNAs in the Tumor Microenvironment. JAMA (2017) 318(10):918–26. doi: 10.1001/jama.2017.11470 Adv Exp Med Biol (2020) 1277:1–31. doi: 10.1007/978-3-030-50224-9_1 13. Park KU, Caudle A. Management of the Axilla in the Patient With Breast Cancer. 7. Tsujimoto M, Nakabayashi K, Yoshidome K, Kaneko T, Iwase T, Akiyama F, Surg Clinics North America (2018) 98(4):747–60. doi: 10.1016/j.suc.2018.04.001 et al. One-Step Nucleic Acid Amplification for Intraoperative Detection of 14. Hamam R, Hamam D, Alsaleh KA, Kassem M, Zaher W, Alfayez M, et al. Lymph Node Metastasis in Breast Cancer Patients. Clin Cancer Res (2007) 13 Circulating microRNAs in Breast Cancer: Novel Diagnostic and Prognostic (16):4807–16. doi: 10.1158/1078-0432.CCR-06-2512 Biomarkers. Cell Death Dis (2017) 8(9):e3045. doi: 10.1038/cddis.2017.440 8. Cserni G. Intraoperative Analysis of Sentinel Lymph Nodes in Breast Cancer 15. Heitzer E, Haque IS, Roberts CES, Speicher MR. Current and Future by One-Step Nucleic Acid Amplification. J Clin Pathol (2012) 65(3):193–9. Perspectives of Liquid Biopsies in Genomics-Driven Oncology. Nat Rev doi: 10.1136/jclinpath-2011-200301 Genet (2019) 20:71–88. doi: 10.3390/cells9020276 Frontiers in Oncology | www.frontiersin.org 12 August 2021 | Volume 11 | Article 627811 Escuin et al. miRNAs and BC Locoregional Metastases 16. Condrat CE, Thompson DC, Barbu MG, Bugnar OL, Boboc A, Cretoiu D, Geometric Averaging of Multiple Internal Control Genes. Genome Biol (2002) et al. miRNAs as Biomarkers in Disease: Latest Findings Regarding Their Role 3(7):RESEARCH0034. doi: 10.1186/gb-2002-3-7-research0034 in Diagnosis and Prognosis. Cells (2020) 9(2):276. doi: 10.3390/cells9020276 39. Robinson MD, McCarthy DJ, Smyth GK. Edger: A Bioconductor Package for 17. Bartel DP, Chen CZ. Micromanagers of Gene Expression: The Potentially Differential Expression Analysis of Digital Gene Expression Data. Widespread Influence of Metazoan microRNAs. Nat Rev Genet (2004) 5 Bioinformatics (2010) 26(1):139–40. doi: 10.1093/bioinformatics/btp616 (5):396–400. doi: 10.1038/nrg1328 40. Benjamini Y, Hochberg Y. Controlling the False Discovery Rate: A Practical 18. Croce CM. Causes and Consequences of microRNA Dysregulation in Cancer. and Powerful Approach to Multiple Testing. J R Stat Soc Ser B (Methodol) Nat Rev Genet (2009) 10(10):704–14. doi: 10.1038/nrg2634 (1995) 57(1):289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x 19. Shen J, Stass SA, Jiang F. MicroRNAs as Potential Biomarkers in Human Solid 41. Valastyan S. Roles of microRNAs and Other Non-Coding RNAs in Breast Tumors. Cancer Lett (2013) 329(2):125–36. doi: 10.1016/j.canlet.2012.11.001 Cancer Metastasis. J Mammary Gland Biol Neoplasia (2012) 17(1):23–32. doi: 20. Nicoloso MS, Spizzo R, Shimizu M, Rossi S, Calin GA. MicroRNAs–the Micro 10.1007/s10911-012-9241-9 Steering Wheel of Tumour Metastases. Nat Rev Cancer (2009) 9(4):293–302. 42. Li CY, Xiong DD, Huang CQ, He RQ, Liang HW, Pan DH, et al. Clinical Value of doi: 10.1038/nrc2619 miR-101-3p and Biological Analysis of Its Prospective Targets in Breast Cancer: A 21. Chen W, Cai F, Zhang B, Barekati Z, Zhong XY. The Level of Circulating miRNA-10b Study Based on The Cancer Genome Atlas (TCGA) and Bioinformatics. Med Sci and miRNA-373 in Detecting Lymph Node Metastasis of Breast Cancer: Potential Monit: Int Med J Exp Clin Res (2017) 23:1857–71. doi: 10.12659/MSM.900030 Biomarkers. Tumour Biol (2013) 34(1):455–62. doi: 10.1007/s13277-012-0570-5 43. Liu P, Ye F, Xie X, Li X, Tang H, Li S, et al. Mir-101-3p is a Key Regulator of 22. Si H, Sun X, Chen Y, Cao Y, Chen S, Wang H, et al. Circulating microRNA-92a and Tumor Metabolism in Triple Negative Breast Cancer Targeting AMPK. microRNA-21 as Novel Minimally Invasive Biomarkers for Primary Breast Cancer. Oncotarget (2016) 7(23):35188–98. doi: 10.18632/oncotarget.9072 J Cancer Res Clin Oncol (2013) 139(2):223–9. doi: 10.1007/s00432-012-1315-y 44. Wang W, Peng B, Wang D, Ma X, Jiang D, Zhao J, et al. Human Tumor 23. Shaker O, Maher M, Nassar Y, Morcos G, Gad Z. Role of microRNAs -29b-2, microRNA Signatures Derived From Large-Scale Oligonucleotide Microarray -155, -197 and -205 as Diagnostic Biomarkers in Serum of Breast Cancer Datasets. Int J Cancer (2011) 129(7):1624–34. doi: 10.1002/ijc.25818 Females. Gene (2015) 560(1):77–82. doi: 10.1016/j.gene.2015.01.062 45. Zhang LY, Ho-Fun Lee V, Wong AM, Kwong DL, Zhu YH, Dong SS, et al. 24. Huo D, Clayton WM, Yoshimatsu TF, Chen J, Olopade OI. Identification of a MicroRNA-144 Promotes Cell Proliferation, Migration and Invasion in Circulating microRNA Signature to Distinguish Recurrence in Breast Cancer Nasopharyngeal Carcinoma Through Repression of PTEN. Carcinogenesis Patients. Oncotarget (2016) 7(34):55231–48. doi: 10.18632/oncotarget.10485 (2013) 34(2):454–63. doi: 10.1093/carcin/bgs346 25. Webber C, Gospodarowicz M, Sobin LH, Wittekind C, Greene FL, Mason MD, 46. Chang CW, Wu HC, Terry MB, Santella RM. microRNA Expression in et al. Improving the TNM Classification: Findings From a 10-Year Continuous Prospectively Collected Blood as a Potential Biomarker of Breast Cancer Risk Literature Review. Int J Cancer (2014) 135(2):371–8. doi: 10.1002/ijc.28683 in the BCFR. Anticancer Res (2015) 35(7):3969–77. 26. Elston CW, Ellis IO. Pathological Prognostic Factors in Breast Cancer. I. The 47. Wu ZS, Wu Q, Wang CQ, Wang XN, Wang Y, Zhao JJ, et al. MiR-339-5p Value of Histological Grade in Breast Cancer: Experience From a Large Study Inhibits Breast Cancer Cell Migration and Invasion In Vitro and may be a With Long-Term Follow-Up. Histopathology (1991) 19(5):403–10. doi: Potential Biomarker for Breast Cancer Prognosis. BMC Cancer (2010) 10:542. 10.1111/j.1365-2559.1991.tb00229.x doi: 10.1186/1471-2407-10-542 27. Appierto V, Callari M, Cavadini E, Morelli D, Daidone MG, Tiberio P. A Lipemia- 48. Wang B, Li J, Sun M, Sun L, Zhang X. miRNA Expression in Breast Cancer Independent NanoDrop((R))-Based Score to Identify Hemolysis in Plasma and Varies With Lymph Node Metastasis and Other Clinicopathologic Features. Serum Samples. Bioanalysis (2014) 6(9):1215–26. doi: 10.4155/bio.13.344 IUBMB Life (2014) 66(5):371–7. doi: 10.1002/iub.1273 28. Blondal T, Jensby Nielsen S, Baker A, Andreasen D, Mouritzen P, Wrang 49. Li Y, Zhao W, Bao P, Li C, Ma XQ, Li Y, et al. miR-339-5p Inhibits Cell Teilum M, et al. Assessing Sample and miRNA Profile Quality in Serum and Migration and Invasion In Vitro and may be Associated With the Tumor- Plasma or Other Biofluids. Methods (2013) 59(1):S1–6. doi: 10.1016/ Node-Metastasis Staging and Lymph Node Metastasis of non-Small Cell Lung j.ymeth.2012.09.015 Cancer. Oncol Lett (2014) 8(2):719–25. doi: 10.3892/ol.2014.2165 29. Martin M. Cutadapt Removes Adapter Sequences From High-Throughput 50. Shan W, Li J, Bai Y, Lu X. miR-339-5p Inhibits Migration and Invasion in Sequencing Reads. EMBnetjournal (2011) 17(1):3. doi: 10.14806/ej.17.1.200 Ovarian Cancer Cell Lines by Targeting NACC1 and BCL6. Tumour Biol 30. Cock PJ, Fields CJ, Goto N, Heuer ML, Rice PM. The Sanger FASTQ File (2016) 37(4):5203–11. doi: 10.1007/s13277-015-4390-2 Format for Sequences With Quality Scores, and the Solexa/Illumina FASTQ 51. Wang YL, Chen CM, Wang XM, Wang L. Effects of miR-339-5p on Invasion Variants. Nucleic Acids Res (2010) 38(6):1767–71. doi: 10.1093/nar/gkp1137 and Prognosis of Hepatocellular Carcinoma. Clinics Res Hepatol Gastroenterol 31. Robinson MD, Oshlack A. A Scaling Normalization Method for Differential (2016) 40(1):51–6. doi: 10.1016/j.clinre.2015.05.022 Expression Analysis of RNA-Seq Data. Genome Biol (2010) 11(3):R25. doi: 52. Zhou J, Xu T, Yan Y, Qin R, Wang H, Zhang X, et al. MicroRNA-326 Functions as 10.1186/gb-2010-11-3-r25 a Tumor Suppressor in Glioma by Targeting the Nin One Binding Protein 32. Wu Y, Wei B, Liu H, Li T, Rayner S. MiRPara: A SVM-Based Software Tool (NOB1). PloS One (2013) 8(7):e68469. doi: 10.1371/journal.pone.0068469 for Prediction of Most Probable microRNA Coding Regions in Genome Scale 53. Zhou C, Liu G, Wang L, Lu Y, Yuan L, Zheng L, et al. MiR-339-5p Regulates the Sequences. BMC Bioinf (2011) 12:107. doi: 10.1186/1471-2105-12-107 Growth, Colony Formation and Metastasis of Colorectal Cancer Cells by 33. R_Development_Core_Team. R: A Language and Environment for Statistical Targeting PRL-1. PloS One (2013) 8(5):e63142. doi: 10.1371/journal.pone.0063142 Computing. Vienna, Austria: the R Foundation for Statistical Computing 54. Zhang L, Lu XQ, Zhou XQ, Liu QB, Chen L, Cai F. NEAT1 Induces (2011). Available at: http://www.R-project.org/, ISBN: . Osteosarcoma Development by Modulating the miR-339-5p/TGF-Beta1 34. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene Pathway. J Cell Physiol (2018) 234(4):5097–105. doi: 10.1002/jcp.27313 Ontology: Tool for the Unification of Biology. Gene Ontol Consort Nat Genet 55. Wu Q, Liu X, Yan H, He YH, Ye S, Cheng XW, et al. B-Cell Lymphoma 6 (2000) 25(1):25–9. doi: 10.1038/75556 Protein Stimulates Oncogenicity of Human Breast Cancer Cells. BMC Cancer 35. Gene Ontology C. Gene Ontology Consortium: Going Forward. Nucleic Acids (2014) 14:418. doi: 10.1186/1471-2407-14-418 Res (2015) 43(Database issue):D1049–56. doi: 10.1093/nar/gku1179 56. Yan H, Zhao M, Huang S, Chen P, Wu WY, Huang J, et al. Prolactin Inhibits BCL6 36. Alexa A, Rahnenfuhrer J, Lengauer T. Improved Scoring of Functional Groups Expression in Breast Cancer Cells Through a MicroRNA-339-5p-Dependent From Gene Expression Data by Decorrelating GO Graph Structure. Pathway. J Breast Cancer (2016) 19(1):26–33. doi: 10.4048/jbc.2016.19.1.26 Bioinformatics (2006) 22(13):1600–7. doi: 10.1093/bioinformatics/btl140 57. Li Y, Zhang X, Yang Z, Li Y, Han B, Chen LA. miR-339-5p Inhibits Metastasis 37. Andersen CL, Jensen JL, Orntoft TF. Normalization of Real-Time of Non-Small Cell Lung Cancer by Regulating the Epithelial-to-Mesenchymal Quantitative Reverse Transcription-PCR Data: A Model-Based Variance Transition. Oncol Lett (2018) 15(2):2508–14. doi: 10.3892/ol.2017.7608 Estimation Approach to Identify Genes Suited for Normalization, Applied 58. Meulmeester E, Ten Dijke P. The Dynamic Roles of TGF-Beta in Cancer. to Bladder and Colon Cancer Data Sets. Cancer Res (2004) 64(15):5245–50. J Pathol (2011) 223(2):205–18. doi: 10.1002/path.2785 doi: 10.1158/0008-5472.CAN-04-0496 59. Chanda A, Sarkar A, Bonni S. The SUMO System and TGFbeta Signaling 38. Vandesompele J, De Preter K, Pattyn F, Poppe B, Van Roy N, De Paepe A, Interplay in Regulation of Epithelial-Mesenchymal Transition: Implications for et al. Accurate Normalization of Real-Time Quantitative RT-PCR Data by Cancer Progression. Cancers (2018) 10(8):264. doi: 10.3390/cancers10080264 Frontiers in Oncology | www.frontiersin.org 13 August 2021 | Volume 11 | Article 627811 Escuin et al. miRNAs and BC Locoregional Metastases 60. Brabletz T, Kalluri R, Nieto MA, Weinberg RA. EMT in Cancer. Nat Rev 68. Cuk K, Zucknick M, Madhavan D, Schott S, Golatta M, Heil J, et al. Plasma Cancer (2018) 18(2):128–34. doi: 10.1038/nrc.2017.118 microRNA Panel for Minimally Invasive Detection of Breast Cancer. PloS One 61. Montserrat N, Gallardo A, Escuin D, Catasus L, Prat J, Gutierrez-Avigno FJ, (2013) 8(10):e76729. doi: 10.1371/journal.pone.0076729 et al. Repression of E-Cadherin by SNAIL, ZEB1, and TWIST in Invasive 69. Nygren MK, Tekle C, Ingebrigtsen VA, Makela R, Krohn M, Aure MR, et al. Ductal Carcinomas of the Breast: A Cooperative Effort? Hum Pathol (2011) 42 Identifying microRNAs Regulating B7-H3 in Breast Cancer: The Clinical (1):103–10. doi: 10.1016/j.humpath.2010.05.019 Impact of microRNA-29c. Br J Cancer (2014) 110(8):2072–80. doi: 10.1038/ 62. Ismail T, Lee HK, Kim C, Kwon T, Park TJ, Lee HS. KDM1A bjc.2014.113 Microenvironment, Its Oncogenic Potential, and Therapeutic Significance. Epigenet Chromatin (2018) 11(1):33. doi: 10.1186/s13072-018-0203-3 Conflict of Interest: The authors declare that the research was conducted in the 63. Varambally S, Cao Q, Mani RS, Shankar S, Wang X, Ateeq B, et al. Genomic Loss absence of any commercial or financial relationships that could be construed as a of microRNA-101 Leads to Overexpression of Histone Methyltransferase EZH2 in potential conflict of interest. Cancer. Science (2008) 322(5908):1695–9. doi: 10.1126/science.1165395 64. Chan M, Liaw CS, Ji SM, Tan HH, Wong CY, Thike AA, et al. Identification of Publisher’s Note: All claims expressed in this article are solely those of the authors Circulating microRNA Signatures for Breast Cancer Detection. Clin Cancer and do not necessarily represent those of their affiliated organizations, or those of Res (2013) 19(16):4477–87. doi: 10.1158/1078-0432.CCR-12-3401 the publisher, the editors and the reviewers. Any product that may be evaluated in 65. Kodahl AR, Lyng MB, Binder H, Cold S, Gravgaard K, Knoop AS, et al. Novel this article, or claim that may be made by its manufacturer, is not guaranteed or Circulating microRNA Signature as a Potential non-Invasive Multi-Marker endorsed by the publisher. Test in ER-Positive Early-Stage Breast Cancer: A Case Control Study. Mol Oncol (2014) 8(5):874–83. doi: 10.1016/j.molonc.2014.03.002 Copyright © 2021 Escuin, Loṕ ez-Vilaro,́ Mora, Bell, Moral, Peŕ ez, Arqueros, Garcıa- ́ 66. Shen J, Hu Q, Schrauder M, Yan L, Wang D, Medico L, et al. Circulating miR- Valdecasas, Ramoń y Cajal, Lerma and Barnadas. This is an open-access article 148b and miR-133a as Biomarkers for Breast Cancer Detection. Oncotarget distributed under the terms of the Creative Commons Attribution License (CC BY). The (2014) 5(14):5284–94. doi: 10.18632/oncotarget.2014 use, distribution or reproduction in other forums is permitted, provided the original 67. Shimomura A, Shiino S, Kawauchi J, Takizawa S, Sakamoto H, Matsuzaki J, author(s) and the copyright owner(s) are credited and that the original publication in this et al. Novel Combination of Serum microRNA for Detecting Breast Cancer in journal is cited, in accordance with accepted academic practice. No use, distribution or the Early Stage. Cancer Sci (2016) 107(3):326–34. doi: 10.1111/cas.12880 reproduction is permitted which does not comply with these terms. Frontiers in Oncology | www.frontiersin.org 14 August 2021 | Volume 11 | Article 627811