We selected 35 MEC samples from patients treated at the A.C.Camargo Cancer Center and Barretos Cancer Hospital, Barretos, São Paulo, Brazil. Two experienced pathologists (FPS and VCA) in salivary gland tumors reviewed the diagnosis of all tumor cases and graded according to Auclair et al., 1992 (4). Demographic, clinical, pathological, therapeutic, and follow-up data were obtained from the patients’ medical records (Table 1). A reference RNA (Human Universal Reference Total RNA, Clontech, Mountain View, California, USA) was used and hybridized with both tumor RNA and normal salivary gland RNA. Six surrounding normal salivary glands (NSG/control) tissues were removed during surgical procedures of six MEC patients, and they were hybridized with reference RNA to further compare their mRNA and miRNA expressions with MEC’s (tumor) mRNA and miRNA expressions. All samples were collected from treatment-naive patients. Written informed consent was obtained from all patients before the sample collection. The National Human Research Ethics Committee approved the study (Protocol #1.380.762/2015).

miRNA expression analyses were performed in 25 out of the 35 fresh-frozen MEC samples and six NSG; no tissue or total RNA was available for analyses in the remaining 10 samples (Table 1). Hybridizations were performed using a one-color SurePrint 8X60K Human miRNA platform (G4870A, Agilent Technologies, Santa Clara, CA, USA), as recommended by the supplier. Background correction, quantile normalization, log2 transformation, and statistical tests were conducted using BRB ArrayTools software v. 4.4.0 (Biometric Research Branch, National Cancer Institute, Bethesda, MD, USA - https://brb.nci.nih.gov/BRB-ArrayTools/index.html). Sequences with more than 10% of MEC and NSG samples presenting undetectable expression (below background signal) were removed. The mean of the probes representing the same miRNA was used in the subsequent steps. miRNAs DE between MEC and NSG groups were identified with a p-value <0.05 (random variance t-test), false discovery rate (FDR) <0.05, and fold change (FC) ≥ 2 and ≤ -2. Supervised hierarchical clustering analysis was performed using 1-minus correlation distance and complete linkage (BRB array tools). Robustness of hierarchical clustering analyses was confirmed using pvclust package (R program) (Supplementary Figure S1). Data were deposited in the Gene Expression Omnibus (GEO) database with the accession number GSE199692.

Array-based gene expression analysis was performed in 34 out of the 35 fresh-frozen MEC samples and five NSG; one MEC and one NGS sample were excluded based on inferior RNA quality (Table 1). Hybridizations were performed using Two-color SurePrint G3 Human Gene Expression Microarray 8x60K (G4851B, Agilent) platform, as previously described (27). Data processing and analyses were carried out using similar parameters described for miRNA profiling (BRB array tools). Identification of DE mRNAs (p-value < 0.001, FDR < 0.05, FC ≥ 2 and ≤ -2) and supervised hierarchical clustering analysis were performed as described above. Pvclust package (R program) was used to confirm the robustness of hierarchical clustering analyses (Supplementary Figure S1). The data were deposited in the GEO database (accession number GSE169754).

Target transcripts from the disrupted miRNAs were predicted using the miRWalk 2.0 tool (http://www.umm.uni-heidelberg.de/apps/zmf/mirwalk/), considering only the interactions predicted by at least three of four different bioinformatic algorithms (miRWalk, miRanda, RNAhybrid, and Targetscan). miRNA and mRNA expression data from 24 MEC samples tested by both procedures were integrated based on a significant negative correlation (Pearson correlation, p-value < 0.05) between predicted miRNA-mRNA interactions. Experimentally validated interactions were additionally obtained from the miRTarBase database (28).

Pathway enrichment analysis was performed with KOBAS 3.0 (http://kobas.cbi.pku.edu.cn) and pathDIP (http://ophid.utoronto.ca/pathDIP) tools, comprising PANTHER, Reactome, and KEGG databases. Default parameters were adopted in KOBAS 3.0, and only experimentally detected protein-protein interactions were considered in PathDIP. The threshold used in both in silico tools was defined as p-value < 0.001 (hypergeometric test) and adjuscted p-value < 0.05 (Benjamini and Hochberg method).

Human Mucoepidermoid Carcinoma (UM-HMC-2) cells were isolated from the intermediate grade (stage IVb) parotid gland MEC of a 59-year-old Caucasian female and cultured according to Warner et al. (29).

miRNA-22 and miR-205 expression in UM-HMC-2 was confirmed using qRT-PCR (data not shown). Then, UM-HMC-2 cells were transfected with pSpCas9(BB)-2A-GF (PX458 expression vector, Addgene plasmid # 48138) expressing CRISPR-Cas9 and sgRNA targeting either miR-22 or miR-205 using Fugene HD transfection reagent (Promega, Madison, WI, USA). This resulted in transient expression of Cas9-sgRNA. Cells transfected with an empty plasmid were used as a control. After 72 hours, cells were sorted for GFP (Green fluorescent protein) positive population using a Sony SH800 cell sorter (Sony Biotechnology, San Jose, CA, USA), and were cloned as single cells per well in a flat bottom 96-well plate. Successfully expanded clones were then screened by capillary sequencing to detect nonhomologous end-joining CRISPR-Cas9 induced gene editing. Clones with predicted out of frame insertions and deletions (indels) were selected and expanded. The predicted effect of the CRISPR editing on miRNAs was assessed using the TIDE tool (30). Details of all sgRNAs and primers used in the experiments, as well as the CRISPR knockout efficiency, are summarized in the supplementary information (Supplementary Figure S2A and Supplementary Table S1).

In addition to sequencing, CRISPR knockout of miR-22 and miR-205 was confirmed using qRT-PCR. The miRNA was extracted with miRNeasy Tissue/Cell Advanced Mini Kit (Qiagen, Hilden, Germany) and transcripted to cDNA using miScript II RT Kit (Qiagen) following manufacturer’s instructions. The miScript universal primer and miRNA-specific primers for Hsa-miR-22-3p (MS00003220) and hsa-miR-205-5p (MS00003780) were purchased from Qiagen. The relative quantitative expressions were normalized to the endogenous control human RNU-6 (MS00033740) purchased also from Qiagen. Quantitative real-time PCR was performed on Applied Biosystems QuantStudio 5 Real-Time PCR System. qRT-PCR results are summarized in the supplementary information (Supplementary Figure S2B).

In order to study the effect of miR-knockouts, 11 genes were selected and evaluated by qRT-PCR: PTEN, LAMC1, CADM1, HER3, MYCBP, SNAI1, YAP1, CD147, SMAD4, ESR1 (ESR1) and ZEB2. One thousand ng of the total RNA was used for cDNA synthesis. Synthesis was done using iScript cDNA Synthesis Kit (Bio-Rad Laboratories, Hercules, CA, USA) according to the manufacturer’s instructions. Two nanograms of cDNA was used for performing qRT-PCR with the Fast SYBR Green Master Mix (Thermo Fisher Scientific) as per the manufacturer’s instructions. The relative quantitative expression was normalized to the endogenous control GAPDH. The primers were purchased from Metabion (Planegg, Germany) and the sequences are summarized in the supplementary information (Supplementary Table S1). Quantitative real-time PCR was performed on Applied Biosystems QuantStudio 5 Real-Time PCR System.

A CellTiter-Glo (CTG) 2.0 Luminescent Cell Viability Assay (Promega, Madison, WI, USA) was used to determine the effect of miR-22 and miR-205 on the cells’ viability. Briefly, 100 μL of cell suspension was dispensed in the Perkin Elmer ViewPlate-96 microplate with a clear flat bottom and black well walls (Perkin Elmer Inc., Waltham, MA, USA) for a final concentration of 1000 cells per well. After 72 hours, 100 μL of the CellTiter-Glo reagent was dispensed into the wells, and the luminescence reads were measured using a PHERAstar plate reader (BMG Labtech, Ortenberg, Germany).

IncuCyte 96-well ImageLock Microplate wells (Sartorius, Göttingen, Germany) were coated with 300 μg/mL Myogel for migration and invasion assays (31). The cells were seeded at a density of 25,000 cells per well in 100 μL of complete medium for both assays. After 24 hours at 37°C, a 96-pin IncuCyte WoundMaker Tool (Sartorius) was used to make uniform wounds on the confluent monolayer of the cell. The wells were washed two times with media, and 100 μL of complete medium was added. For the invasion plate, 50 μL of Myogel-collagen gel (2.4 mg/mL Myogel, 0.8 mg/mL type I rat tail collagen) (Corning Incorporated, Corning, NY, USA) was added on top of the cells. After the gel was solidified, 50 μL of media was added, and the plates were transferred to an incubator. The wound closing was monitored automatically every 2 hours for two days using IncuCyte S3 Live-Cell Imaging System (Sartorius). Analysis of wound closing (width of the wound) was performed using Matlab. Mathematical function decorrelation was used to make the cells’ intensity substantially higher than the background.

The spheroid invasion assay was done according to Naakka et al. (32). The UM-HMC-2 cells were seeded at a concentration of 1000 cells per well in 50 µL of the complete medium using a U-shaped ultra-low attachment 96-well plate (Corning, New York, USA) and incubated for four days. Next, the spheroids were embedded in 50 μL Myogel-fibrin gel containing 0.5 mg/mL Myogel, 0.3 U/mL thrombin (Sigma-Aldrich), 33.3 mg/mL aprotinin (Sigma-Aldrich), and 0.5 mg/mL fibrinogen (Merck). After the Myogel-fibrin matrix (30 min) solidification, 100 μL of complete medium was added to the wells.

Images of the spheroids were captured daily using Nikon Eclipse TS100 Inverted Microscope (Nikon, Minato, Tokyo, Japan) at 4x magnification. Analysis of the spheroid invasion area and length of the longest branch was performed using ilastik (freeware) and Fiji ImageJ 1.51 software (33).

All in vitro assays were repeated at least three times, each performed at least in triplicate. Statistical analyses were carried out with SPSS v.25.0 (IBM Corporation, Chicago, IL, USA) and the GraphPad Prism software (v. 6.0; GraphPad Software Inc., La Jolla, CA, USA). Student’s T-test or One-way ANOVA followed by Bonferroni correction was used in posthoc analysis. A p-value < 0.05 was considered as statistically significant. Figures were created with Origin 2018b graphing software (OriginLab Corporation, Northampton, MA, USA).

Overall survival analysis was performed using the Kaplan-Meier estimator with the log-rank test in IBM SPSS Statistics for Windows, Version 25.0 (IBM Corp. Armonk, NY, USA). The miRNA expression values were dichotomized below and above the median (p-value < 0.05). A random variance t-test using BRB ArrayTools software (v. 4.4.0) was applied to investigate differences in the miRNA expression in relation to the histological grade, lymph node, and distant metastasis (p-value < 0.05, FDR < 0.05).

After excluding uniformly low expressed miRNAs in MEC and NSG samples, 530 miRNAs and 19,911 mRNAs were considered for further analysis. We found 46 DE miRNAs (18 overexpressed and 28 underexpressed) in MEC (Supplementary Table S2). The most significant (P adjusted ≤ 0.005) overexpressed miRNAs included miR-21-5p (FC=10.2), miR-22-3p (FC=2.0), miR-181a-5p (FC=2.9), miR-205-3p (FC=14.7), and miR-224-3p (FC=5.3). The miR-363-3p (FC =-16.3), miR-625-5p (FC =-18.5), miR-885-5p (FC =-10.7), miR-892b (FC =-2.7), and miR-1288-3p (FC =-2.7) were significantly underexpressed (Figure 2). A similar approach used for mRNAs unveiled 3,162 mRNAs differentially expressed in MEC (1,488 overexpressed and 1,674 underexpressed (Supplementary Table S3).

Supervised hierarchical clustering analysis based on the DE transcripts revealed two clusters in both miRNA (Figure 3A) and mRNA (Figure 3B). Although these two main clusters completely separated MEC from NSG samples, no association was observed when comparing the clinical-pathological parameters (histological grade, lymph node involvement, and distant metastasis) with the clusters generated by both miRNA and mRNA analysis.

The miRNA-target prediction analysis resulted in 20,816 miRNA-mRNA putative interactions. The integrative analysis revealed a miRNA-mRNA network comprising 696 negatively correlated interactions and inverted FCs (44 miRNAs and 444 mRNAs) (Supplementary Table S4). The main biological pathways uncovered by miRNA targets corroborated by the integrative analysis were cell signaling, cell cycle, and cancer-related pathways (Table 2).

Lower expression levels of miR-582-5p, miR-3125, and miR-4324 were found in high-grade MEC compared to low and intermediate grades (Supplementary Figure S3). Increased expression levels of miR-205-5p and miR-224-5p (both overexpressed in MEC) and decreased expression levels of miR-139-3p, miR-145-3p, miR-148a-3p, miR-186-5p, miR-338-3p, miR-363-3p and miR-4324 were significantly related to worse overall survival in MEC patients (Figure 4).

Among the list of differentially expressed miRNAs, we selected miR-22 and miR-205 for functional assays for the following reasons: they were significantly overexpressed (adjuscted p-value <0.005) (Figure 2), presented high interactivity in the integrative analysis (>10 underexpressed mRNA predicted targets negatively correlated with the miRNA expression) (Supplementary Table S4), and showed clinical association with worse prognosis (increased miR-205 expression was associated with shorter overall survival) (Figure 4).

We explored the use of the CRISPR-Cas9-based method to knockout miR-22 and miR-205 in MEC. Cell viability, migration, and invasion assays were performed in the MEC cell line UM-HMC-2. The cell viability was the lowest in the miR-205 knockout, followed by miR-22-knockout cells, but with no statistical significance (Figure 5A).

The scratch wound migration assay showed that miR-22 and miR-205 knockouts reduced cell migration. The effect was the same in both knockout cell lines compared to the empty vector, but miR-22-knockout cells migrated significantly slower than the control cells (Figures 5B–D and Supplementary Figure S4A). The cell lines showed different invasion speeds, and the knockouts invaded slower in both scratch wound invasion (Figures 5E–F and Supplementary Figure S4A) and spheroid invasion assays (Figures 5G–K and Supplementary Figure S4B). Both miR-22 and miR-205 knockout cell lines invaded slower than the cell line with empty gRNA vectors. However, the effect was statistically significant only when miR-22-knockout cells were compared to the empty vector in the spheroid invasion assay.

In order to understand the mechanism behind the effect of miR-22 and miR-205 knockout on MEC cell behaviour, we studied the expression of specific molecules: PTEN, LAMC1, CADM1, HER3, MYCBP, SNAI1, YAP1, CD147, SMAD4, ESR1 and ZEB2 which, based on the literature, are known to be targets either for miR-22 or miR-205. We reported significant differences in two of the targets: estrogen receptor alpha (ESR1) for miR-22 and zinc finger E-box-binding homeobox 2 (ZEB2) and miR-205 (Figure 6). These molecules influence cell proliferation, migration and invasion (34–36). As expected, miR-22 knockout cells have significantly higher expression of ESR1, and miR-205 knockout cells have significantly higher expression of ZEB2 compared with the empty vector.