miRNA profiles of canine cutaneous mast cell tumours with early nodal metastasis and evaluation as potential biomarkers

Cutaneous mast cell tumours (MCTs) are common skin neoplasms in dogs. MicroRNAs (miRNAs) are post-transcriptional regulators involved in several cellular processes, and they can function as tumour promoters or suppressors. However, the role of miRNAs in canine MCTs has not yet been elucidated. Thus, the current study aimed to characterize miRNA profiles and to assess their value as biomarkers for MCTs. miRNA expression profiles were assessed in formalin-fixed, paraffin-embedded samples by next-generation sequencing. Ten samples were MCT tissues, and 7 were healthy adjacent tissues. Nine dysregulated miRNAs (DE-miRNAs) were then validated using RT-qPCR in a larger group of MCT samples, allowing the calculation of ROC curves and performance of multiple factor analysis (MFA). Pathway enrichment analysis was performed to investigate miRNA biological functions. The results showed that the expression of 63 miRNAs (18 up- and 45 downregulated) was significantly affected in MCTs. Five DE-miRNAs, namely, miR-21-5p, miR-92a-3p, miR-338, miR-379 and miR-885, were validated by RT-qPCR. The diagnostic accuracy of a panel of 3 DE-miRNAs—miR-21, miR-379 and miR-885—exhibited increased efficiency in discriminating animals with MCTs (AUC = 0.9854) and animals with lymph node metastasis (AUC = 0.8923). Multiple factor analysis revealed clusters based on nodal metastasis. Gene Ontology and KEGG analyses confirmed that the DE-miRNAs were involved in cell proliferation, survival and metastasis pathways. In conclusion, the present study demonstrated that the miRNA expression profile is changed in the MCT microenvironment, suggesting the involvement of the altered miRNAs in the epigenetic regulation of MCTs and identifying miR-21, miR-379 and miR-885 as promising biomarkers.

a non-denaturing acrylamide gel (≈ 141 bp bands). After elution from the gel, the size-selected libraries were quantified and sequenced on a NextSeq 500 sequencer (Illumina). The resulting number of reads per sample varied from 53,000 to 29,000,000. Eight MCT samples with insufficiently high mapping rates were excluded from further analysis. For two dogs (numbers 8 and 15 in Table 1), the results for the tumours and matched healthy adjacent margins were reported.
A count table was used to identify differentially expressed miRNAs via DESeq2 analysis 15 . The DESeq threshold was set by discarding low-expression miRNAs having an average count of 2 or less. This analysis revealed the expression of 246 and 116 Canis familiaris (cfa) miRNAs in healthy margins and MCTs, respectively.
The expression profiles of sequenced samples were used to carry out cluster analysis. Samples were grouped into two clusters: healthy tissues and MCTs (Fig. 1a).
To rank the most differentially expressed miRNAs (DE-miRNAs) between healthy and MCT samples, the results of differential expression analysis performed with DESeq2 were further filtered a more stringent cut-off criteria of an adjusted P value of 0.01 and an absolute log2FC of 2.4. This filtering allowed the identification of sixty-three miRNAs whose abundance differed significantly between MCT and healthy samples, demonstrating that 45 miRNAs were downregulated with a log2FC of between − 2.4 and − 13.4 and 18 miRNAs were upregulated with a log2FC of between 2.4 and 6.9 (Fig. 1b).

Quantification of DE-miRNAs in healthy versus MCT samples by RT-qPCR.
To validate the NGS results and measure the abundances of DE-miRNAs in MCTs, RT-qPCR was performed on both the sequenced samples and a separate group of 11 MCTs and related healthy adjacent (normal) tissue samples. To validate the sequencing results, 9 DE-miRNAs-miR-370, miR-379, miR-92a, miR-21, miR-26a, miR-342, miR-885, miR-375 and miR-338-were selected based on the fold change and read count values. MiR-122, miR-128 and miR- Table 1. Summary of clinical and histopathological data. a MCT samples sequenced using NGS. b Healthy margins sequenced using NGS. c $amples in which miRNAs selected for the RT-qPCR validation step were not detected. d Samples for which healthy margins were not collected. e Classification system proposed by Weishaar et al. 7 . HN histological node, NGS next-generation sequencing. Diagnostic value of DE-miRNAs in dogs with MCT. Receiver operating characteristic (ROC) analysis was used to assess the diagnostic value of DE-miRNAs as biomarkers to further discriminate between MCT and healthy adjacent tissue. To confirm the diagnostic efficacy of each miRNA, the associated area under the curve (AUC) was calculated. Table 2 shows a summary of the diagnostic performance of each DE-miRNA and shows combinations of three DE-miRNAs. The AUC was fair for miR-92a (AUC = 0.7427) and miR-338 (AUC = 0.7339) and excellent for miR-21 (AUC = 0.9825), miR-379 (AUC = 0.9211) and miR-885 (AUC = 0.9181) (Fig. 3). Discriminant analysis was used to further investigate the potential for improving diagnostic performance by analysing multiple DE-miRNAs. Statistical analysis was performed to examine the weighted average relative quantification (RQ) values of the miRNAs with an AUC of > 0.9 (miR-21, miR-379 and miR-885) (Fig. 4). The median expression levels were 0.0301 (range 0.0069-0.9334) and 0.99998 (range 0.3485-1) in healthy margin and MCT samples, respectively (Fig. 4a). The predicted probability of being able to discriminate a sample as positive based on the logit model [logit = 1/(1 + exp (− (− 4.92611-1.31822 × expression level of miR-885 + 0.40746 × expression level of miR-379 + 0.86787 expression level of miR-21)))] was used to construct the ROC curve (Fig. 4b). The AUC for the panel of these three DE-miRNAs was 0.9854 (95% CI 0.9854-0.9854) with a cut-off value of 0.1654, and a sensitivity of 100% and a specificity of 94.4%.
Potential of the miRNA panel for the detection of nodal metastases. Excised lymph nodes were categorized in accordance with the Weishaar classification system for nodal metastases (2014) 7 , and the potential of the three-miRNA panel to discriminate patients with and without lymph node involvement was evaluated. Two groups, namely, HN0/1, including non-metastatic and pre-metastatic samples, and HN2, including early metastatic samples, were included for further analysis ( Table 1). The weighted average relative quantification (RQ) values of the miRNA panel (miR-21, miR-379 and miR-885) were calculated (Fig. 5). The median expression levels were 0.3179 (range 0.0071-0.8858) and 0.9424 (range 0.3741-1) in the HN0/1 and HN2 groups, respectively (Fig. 5a). The predicted probability of being able to discriminate a sample as metastatic based on the logit model [logit = 1/(1 + exp (−(− 1.58980-7.91569 × expression level of miR-885 + 0.13130 × expression level of miR-379 + 0.05084 expression level of miR-21)))] was used to construct the ROC curve (Fig. 5b). The AUC for the panel of these three DE-miRNAs was 0.8923 (95% CI 0.759-1.000) with a cut-off value of 0.5528, a sensitivity of 92.3% and a specificity of 80%. MFA identifies individuals with similar profiles who are close to each other on the factor map. Collectively, the components F1 and F2 explained 68.53% of the total variance in the data (Fig. 5c). The first component (F1) explained 42.76% of the variance and separated the HN0/1 group from the HN2 group according to lymph node involvement. The second component (F2) explained 25.77%, discriminating non-metastatic HN0 samples (samples 4 and 10) in the upper right panel from high-grade early metastatic HN2 samples (samples 16 and 9) in the lower-left panel (Supplementary Table S1).
Gene Ontology and pathway enrichment analysis of miRNAs. The MiRWalk 3.0 and DAVID databases were searched to retrieve the candidate target genes of DE-miRNAs and to perform mRNA enrichment analysis, respectively. Of the predicted mRNA targets of downregulated miRNAs, 196 were in the 3′UTR, 45 were in the 5′UTR and 171 were in the CDS. Of the predicted mRNA targets of upregulated miRNAs, 16 were in the 3′UTR, 3 were in the 5′UTR and 11 were in the CDS. The list of candidate target genes is provided in Table 3.
Gene Ontology (GO) analysis was performed using DAVID for three categories: biological process (BP), cellular component (CC) and molecular function (MF). An overview of the top 10 terms significantly enriched with target genes for each of the above GO categories is presented in Fig. 6. The enriched GO BP terms mainly included regulation of transcription from RNA polymerase II promoter and protein ubiquitination involved in ubiquitin-dependent protein catabolic process; the CC terms were related to the cytoplasm, nucleus and nucleoplasm, while the MF terms focused on transcription factor activity and sequence-specific DNA binding. KEGG pathway analysis was performed on candidate targets of DE-miRNAs. Figure 7 shows the top 10 significantly enriched KEGG pathways, with PI3K-Akt signalling pathway, small cell lung cancer, viral carcinogenesis and microRNAs in cancer at the top of the list.
Although the prognostic role of both the Patnaik 5 and Kiupel 4 grading systems in canine MCTs is widely accepted, histological grading alone cannot accurately predict the risk of local and distant metastases [3][4][5]9,49 . Nodal metastases have been reported in 20-49% of cutaneous MCTs at first presentation, and identification of lymph node involvement is crucial for accurate tumour staging and prognosis 3,49,50 . Recently, a novel classification system for the evaluation of nodal metastasis in canine MCTs has been proposed and correlated with the clinical outcome, providing evidence that dogs diagnosed with early metastatic/overt metastatic (HN2-HN3) nodes have a shorter life expectancy 7 . In our study, a three-miRNA signature (miR-379-miR-21-miR-885) accurately discriminated between healthy adjacent tissue and MCT tissue (AUC = 0.9854) and identified patients with early nodal metastases (AUC = 0.8923). Since the number of enrolled patients did not allow us to perform discriminant analysis of parameters such as survival time and progression-free interval, the present results provide a background to investigate new biomarkers of MCT outcome in different matrices, including blood, to support clinical decision making.
In conclusion, the present study demonstrated that the expression levels of miR-21, miR-379, miR-92a, miR-885 and miR-338 in the tumour microenvironment are changed compared to those in healthy adjacent tissues and differ in dogs with early nodal metastases compared to those without nodal involvement, suggesting that these miRNAs may epigenetically modulate genes involved in MCT progression and metastasis. Our study provides insights into the emerging roles of miRNAs in veterinary oncology, although more efforts are required to establish the role and molecular targets of the investigated DE-miRNAs. Since the sample size influences the clinical sensitivity and specificity of the test, further studies are necessary to confirm the diagnostic value of miRNAs by increasing the number of patients. . Patients were recruited after the owner provided written informed consent. All experiments were performed following the relevant guidelines and regulations. Samples were trimmed and processed according to currently recommended guidelines 51 (Table 1). Cutaneous MCTs at first presentation without distant metastasis 52 and sentinel/ regional lymph nodes were surgically removed and histologically classified 53 and graded 4,5 . In addition, neoplastic involvement of sentinel lymph nodes was categorized as previously described 7,52 . For all samples, after bright field microscopy observation of the haematoxylin-eosin-stained slide, the corresponding paraffin block was penetrated using a biopsy punch with a plunger (Miltex) to collect a portion of the tumour (21 MCTs) or a portion of the healthy dermal connective tissue (16 margins); the latter samples were used as controls. For MCT samples, areas of necrosis, haemorrhage or inflammation were avoided, if present. miRNA extraction and next-generation sequencing (NGS). MiRNAs were extracted using an miRNeasy FFPE kit (Qiagen, Cat. No. 217504) following the manufacturer's instructions. The RNA quality and quantity were verified according to MIQE guidelines 54 . The RNA concentration was determined in a Qubit 2.0 fluorometer with a Qubit microRNA Assay Kit (Invitrogen, Cat. No. Q32880). A pilot NGS study was performed on 10 MCTs and 7 healthy adjacent tissue samples (Table 1). Small RNA transcripts were converted into barcoded cDNA libraries. Library preparation was performed as previously reported 55 using an NEBNext Multiplex Small RNA Library Prep Set (Cat. No. NEB#E7560) for Illumina, and sequencing was performed in a NextSeq 500 sequencer (Illumina Inc., USA).
Computational analysis. The output of the NextSeq500 Illumina sequencer was demultiplexed using bcl2fastq Illumina software embedded in the docker4seq package 56 . miRNA expression quantification was performed using the workflow and implementation previously described 57 . In brief, after adapter trimming with cutadapt 58 , sequences were mapped using SHRIMP 59 to Canis familiaris precursor miRNAs available in miR-Base 22.0-March 2018 (https ://www.mirba se.org/). Using GenomicsRanges 60 , an R script, was used to identify the number of reads on precursor miRNAs mapping to the expected location on mature miRNAs. The detected counts were organized in a table including all analysed samples. For visualization purposes, only CPM (counts per million reads) values were used. Differential expression analysis was conducted using the DESeq2 Biocon-  To obtain cDNA, reverse transcription was performed using a TaqMan Advanced miRNA cDNA Synthesis Kit (Cat. No. A28007, Applied Biosystems) following the manufacturer's instructions.
Quantitation was performed in a scaled down reaction volume (15 µl) in a CFX Connect Real-Time PCR Detection System (Bio-Rad) using 7.5 μl of 2X TaqMan Fast Advanced Master Mix (catalogue number 4444557), 0.75 μl of miRNA-specific TaqMan Advance assay reagent (20X), 1 μl of cDNA and water to make up the remaining volume. The thermal cycling profile was as follows: 50 °C for 2 min, 95 °C for 3 min and 40 cycles at 95 °C for 15 s and 60 °C for 40 s. No-RT controls and no-template controls were included. The geometric mean of the reference miRNA abundance was used for normalization. Relative quantification of target miRNAs was carried out after sample normalization using the geometric mean of the reference miRNA abundance in Bio-Rad CFX Maestro Software using the 2 -ΔΔCq method. miRNA target prioritization. The target genes of DE-miRNAs were retrieved using MiRWalk 3.0 61 , which includes 3 miRNA target prediction programs (miRDB 62 , miRTarBase 63 and TargetScan 64 ). Analysis was performed on the entire gene sequence (including the 5′UTR, CDS, and 3′UTR). The list of target genes predicted by at least two of the three tools was included in further analysis, mRNA functional enrichment analysis was performed using the DAVID bioinformatic resource 65,66 , and biological pathways in the KEGG database 67 were examined for enrichment. Statistical analysis. Statistical analysis was performed using XLStat software for Windows (Addinsoft, New York, USA). Data were tested for normality using the Shapiro-Wilk test; as the data were not normally distributed, the nonparametric Wilcoxon signed-rank test for paired samples was applied. Quantitative (miRNA quantification and tumour size) and qualitative (lymph node HN classification 7 ) variables were used for ordination analysis using the 'Multiple Factor Analysis' (MFA) function. Receiver operating characteristic (ROC) analysis was performed as previously reported to determine the diagnostic accuracyy 68 . Statistical significance was accepted at a P value of ≤ 0.05.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.