Next Article in Journal
A Molecular Dynamics Study of a Photodynamic Sensitizer for Cancer Cells: Inclusion Complexes of γ-Cyclodextrins with C70
Previous Article in Journal
Curcumin Mitigates Immune-Induced Epithelial Barrier Dysfunction by Campylobacter jejuni
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Micro RNA Transcriptome Profile in Canine Oral Melanoma

1
Veterinary Teaching Hospital, Joint Faculty of Veterinary Medicine, Kagoshima University, Kagoshima, Kagoshima 890-0065, Japan
2
The United Graduate School of Veterinary Science, Yamaguchi University, Yamaguchi, Yamaguchi 753-8515, Japan
3
Joint Graduate School of Veterinary Medicine, Kagoshima University, Kagoshima, Kagoshima 890-0065, Japan
4
Laboratory of Veterinary Surgery, Graduate School of Agricultural and Life Sciences, The University of Tokyo, Bunkyo City, Tokyo 113-8654, Japan
5
Hygiene and Health Promotion Medicine, Kagoshima University Graduate School of Medicine and Dental Science, Kagoshima, Kagoshima 890-8544, Japan
6
Department of Veterinary Histopathology, Joint Faculty of Veterinary Medicine, Kagoshima University, Kagoshima, Kagoshima 890-0065, Japan
7
Animal Medical Centre, Tokyo University of Agriculture and Technology, Tokyo, Tokyo 183-8538, Japan
*
Author to whom correspondence should be addressed.
Int. J. Mol. Sci. 2019, 20(19), 4832; https://doi.org/10.3390/ijms20194832
Submission received: 5 August 2019 / Revised: 13 September 2019 / Accepted: 27 September 2019 / Published: 28 September 2019
(This article belongs to the Section Molecular Oncology)

Abstract

:
MicroRNAs (miRNAs) dysregulation contribute the cancer pathogenesis. However, the miRNA profile of canine oral melanoma (COM), one of the frequent malignant melanoma in dogs is still unrevealed. The aim of this study is to reveal the miRNA profile in canine oral melanoma. MiRNAs profile of oral tissues from normal healthy dogs and COM patients were compared by next-generation sequencing. Along with tumour suppressor miRNAs, we report 30 oncogenic miRNAs in COM. The expressions of miRNAs were further confirmed by quantitative real-time PCR (qPCR). Pathway analysis showed that deregulated miRNAs impact on cancer and signalling pathways. Three oncogenic miRNAs targets (miR-450b, 301a, and 223) from human study also were down-regulated in COM and had a significant negative correlation with their respective miRNA. Furthermore, we found that miR-450b expression is higher in metastatic cells and regulated MMP9 expression through a PAX9-BMP4-MMP9 axis. In silico analysis indicated that miR-126, miR-20b, and miR-106a regulated the highest numbers of differentially expressed transcription factors with respect to human melanoma. Chromosomal enrichment analysis revealed the X chromosome was enriched with oncogenic miRNAs. We comprehensively analyzed the miRNA’s profile in COM which will be a useful resource for developing therapeutic interventions in both species.

1. Introduction

One person dies every hour from melanomas, new melanoma cases have increased by 53% in the US [1], and the WHO reported that 132,000 new melanoma cases are diagnosed every year in the world. These reports clearly confirm the importance of melanoma studies. Molecular studies have enriched the definition of melanoma sub-types [2,3]. The triple wild-type (TWT) subtype bears the features that underlie non-cutaneous melanoma, including mucosal melanoma [3,4]. Human mucosal melanoma is more aggressive with less favourable prognosis than other subtypes.
Previous studies suggested dog melanoma as a natural model for human melanoma [5,6]. Malignant melanoma is frequent in dog and the majorities are in the oral mucosa. Oral melanoma in dogs is considered a typical model for non-UV or TWT melanoma in human [2,3]. Dog melanoma genes have the same mutations or aberrant expression as human melanoma genes, BRAFV600E, NRAS (Q61) [7], PTEN [5], and KIT [8]. Besides the protein-coding RNAs, non-coding RNAs (ncRNAs) also have important roles in gene regulation. Among them, small non-coding RNAs have widespread regulatory functions in human diseases, and microRNAs (miRNAs) are now in phase I trials to treat human hepatitis, diabetes, lymphoma, scleroderma, and lung cancer [9].
Next-generation sequencing (NGS) has been used widely to study miRNA, including their role in human melanoma. For dogs to be a useful therapeutic preclinical model, knowing the miRNA profile in dog melanoma is important. There are few reports of tumour suppressor miRNAs in canine oral melanoma (COM) [10,11], and studies of miRNA profiles of human TWT or mucosal melanoma are scarce. Moreover, the global deregulated miRNAs expression profile of COM is still unrevealed. So we aimed to use next-generation sequencing to study global aberration of miRNAs in COM bearing following three objectives: (1) Comprehensively profile miRNA expression pattern in COM, (2) Validation of findings in the COM tissue samples by using qPCR, and (3) Investigate the gene-regulatory function and molecular pathways of differentially expressed miRNAs. In our study, we obtained the miRNA profile of eight COM tissue samples by NGS which were further validated by qPCR. We found several miRNAs were differentially expressed. We also explored a new function of miR-450b. The impact of global changes in the miRNA profile was shown by Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. Finally, a common miRNA and transcription factor (TF) network were constructed and analyzed to find the most important miRNAs for the regulation of TFs expression between dog and human melanoma.

2. Results

2.1. Small RNA Profile of Canine Oral Melanoma

To investigate the miRNA profile in COM, RNA from three normal oral tissue samples from healthy dogs (hereafter referred to as ‘‘normal’’) and eight samples from dogs with COM was sequenced (Table S1a). After adapter trimming and quality check, we obtained 51 and 142 million clean reads from normal and melanoma libraries, respectively (Table S1b). Sequences were submitted to the SRA database (PRJNA516252). Length distribution analysis showed 90% and 82% of clean reads in the normal and melanoma libraries, respectively, were 20–24-nt long, indicating an alteration in the small RNA profile (Figure S1a). We annotated the reads using miRBase or Ensembl dog and human ncRNAs (see methods), witch 92.5% and 84.23% of the reads in normal and melanoma, respectively, were annotated (Figure S1b). Among the annotated reads miRNAs were the most abundant small RNAs. Interestingly, the percentage of other ncRNAs reads (Mt-rRNA, Mt-tRNA, linc-RNA, sno-RNA, snRNA, misc-RNA, rRNA, other miR) was two times more in the melanoma group (Figure 1a). SnRNA, snoRNA, and mitochondrial tRNA-derived small RNA fragments were the most altered between the two groups. As miRNAs were the most abundant we analyzed the miRNAs further.
We found significant differences in the chromosome distribution of the annotated miRNAs mapped reads between the normal and melanoma groups (Figure 1b), implying altered global miRNA profiles in the melanoma group. We annotated 542 miRNAs in both groups, among them 64 miRNAs were differentially expressed (Figure 1c). The top 20 highly expressed miRNAs made up >70% of the total reads that were annotated to miRBase (Figure 1d). Among them, 12 miRNAs were common between the groups. We selected the rank of the miRNAs on the basis of their expression. The rank orders of miRNA’s were different in melanoma than normal (Figure 1d). Four of the top 10 highly expressed miRNAs in melanoma were not in the top 10 of the normal group. Importantly, miR-21, which is a well-characterized miRNA oncogene frequently found to be over-expressed in various malignancies, was ranked one in melanoma and 12 in normal. Also, miR-22, miR-148a, and miR-186, all of which have been reported to be oncogenic, but not statistically significant in our study also changed their rank within top 10 in melanoma [12,13,14]. However, the ranks of some known anti-oncogenic miRNAs were much lower in melanoma than in normal. For example, miR-203 and miR-205, which were reported to be anti-oncogenic in melanoma [10] were ranked 6 and 7 in normal, were outside the top 20 in melanoma. So, melanoma reduced expression of anti-oncogenic miRNAs while taking favour of others highly expressed miRNAs by remodelling their expression position according to their target functionality.

2.2. Global miRNAs Expression in Canine Oral Melanoma

Unsupervised hierarchical clustering and principal component analysis were performed for all the differentially expressed miRNAs (FC > 1, without considering the FDR and mean read count) to evaluate similarities between the studied samples at the global level (Figure 2a, and Figure S1c,d). After applying stringent filtering criteria (Fold change >2, FDR <0.05, and miRNA mean read counts in either normal or melanoma >50), we obtained 30 up- and 34 down-regulated miRNAs (Figure 2b, Table S2). The heatmap and clustering tree revealed a distinct miRNA expression pattern between the groups. The principal component analysis and clustering tree showed that the differentially expressed miRNA were enough to distinguish the two groups, and heatmap showed the miRNA expression patterns were similar within a group. Data showed there were significant changes in the miRNA profile in COM.

2.3. Validation of miRNA Expression

We selected 20 differentially expressed miRNAs for validation by qPCR by 12 normal oral and 17 melanoma tissue samples (Figure 2c,d and Figure S2). We selected these 20 miRNAs on the basis of three different criteria: (1) miRNAs those were not reported or validated previously (miR-450b, 301a, 140, 542, 223, 190, 429, 96, 375,183, 200b, 141), (2) miRNAs those were reported or validated previously (miR-21, miR-122, 383 and 143) [15], (3) miRNAs those were beyond (numerically close) to the stringent filtering criteria (miR-122, miR-34a, miR-338, miR-182, and miR-1). Among the up-regulated miRNAs, miR-450b, miR-223, miR-140, miR-542, and miR-383 showed >10-fold change and miR-301a, miR-21, and miR-190 showed >5 fold-change (Figure 2c and Figure S2a–c). Among the down-regulated miRNAs, miR-429, miR-200b, miR-141, and miR-375 showed <−10-fold change, miR-96 showed <−3 fold-change, and miR-183 and miR-143 showed <−2-fold change (Figure 2d and Figure S2d). There was significant positive correlation of fold differences between the NGS and qPCR results (Figure 2e), and the expression levels of miR-122, miR-34a, miR-338, miR-182, and miR-1 (Figure 2c,d and Figure S2) which were beyond our stringent filtering criteria, showed similar trends between NGS and qPCR indicates the strength of our filtering criteria. These results confirm the expression of several oncogenic and tumour suppressor miRNAs in COM revealed by NGS and validated by qPCR.

2.4. Gene Regulatory Function of Oncogenic miRNAs

MiRNAs do their function by negatively regulating the gene expression in the mRNA and protein level. To know the differentially expressed miRNAs function, we need to focus on the expression of their respective target genes. However, it is unfeasible to explore all the differentially expressed miRNAs targets in a single study. So, we selected three miRNAs (miR-450b, miR-301a, and miR-223) to know their gene regulatory function because to our knowledge miR-450b has no report in melanoma and other two miRNAs are less studied in human melanoma and no report in COM. We selected PAX9, NDRG2, and ACVR2A as targets of miR-450b, miR-301a, and miR-223, respectively, from previous studies where they validated miRNA-mRNA binding experimentally [16,17,18]. The binding sites in the 3′ un-translated regions of these genes predicted by TargetScan 7.2 were conserved between human and dog (Figure S3). So, we hypothesized that similar phenomenon (miRNA-mRNA binding) may exist in the canine oral melanoma. As in our experiment miR-450b, 301a and 223 was up-regulated so if our hypothesis is true the target genes must be down-regulated and they should have a strong negative correlation with the respective pairs. Our qPCR results showed significant down-regulation of PAX9, NDRG2, and ACVR2a in melanoma compared with normal, and the relative expression of the respective miRNA–mRNA pairs showed significant negative correlation (Figure 3a,b). This inverse relationship indicates the miRNAs may bind the respective mRNA targets like previous human studies and suppress their expression in melanoma. In addition, study reported that PAX9 down-regulation decrease BMP4 expression which can increase MMP9 expression [19,20]. So we investigated the PAX9-BMP4-MMP9 axis in our study. Our qPCR results showed that BMP4 was down-regulated and MMP9 was up-regulated in melanoma tissue samples (Figure 3c). MMP9 is required for the degradation of the extracellular matrix, which is a prerequisite for tumour invasion and positively correlates with tumour metastasis. So we expected high MMP9 expression in metastatic cells along with miR-450b. Therefore, we further investigated the relative expressions of miR-450b, PAX9, BMP4, and MMP9 in two COM cell lines: KMEC established from primary oral melanoma and LMEC from metastatic mandibular lymph node of oral melanoma [21]. qPCR analysis showed that miR-450b was up-regulated and PAX9 and BMP4 were significantly down-regulated in LMEC compared with KMEC. Conversely, as expected, MMP9 was significantly up-regulated in LMEC compared with KMEC, as shown in Figure 3d. So, we concluded that up-regulation of miR-450b and down-regulation of PAX9 and BMP4 were interconnected with the high MMP9 expression in metastatic melanoma cells (Figure 3e).

2.5. Gene Ontology and KEGG Pathway Analysis of the Differentially Expressed miRNAs

To determine the global function of the differentially expressed miRNAs, we predicted their target genes by overlaying the results obtained using TargetScan and miRDB. We detected 2555 and 2464 target genes of the down- and up-regulated miRNAs, respectively (Table S3). We functionally annotated the target genes by assigning gene ontology (GO) terms and KEGG pathways. Target genes of the down-regulated miRNAs were analysed against cancer and other databases (Table S4). We found oncogenic genes related terms were enriched which indicates these genes have an oncogenic function. From GO analysis we found protein modification (e.g., phosphorylation, transcription, or repression from DNA), extracellular matrix, and receptor signalling GO terms were assigned for the target genes of down-regulated miRNAs. It indicates down-regulated miRNAs inhibit their target genes to maintain target genes terms related function in normal condition which was disrupted in melanoma due to the miRNA down-regulation. Target genes of the up-regulated miRNAs were involved mainly in the protein ubiquitination system because ubiquitin-dependent protein catabolic process and ubiquitin-protein ligase activity terms were enriched in GO (Table 1). Ubiquitylation is long known for driving cell cycle transition and ubiquitin ligase has relation to the cell cycle control. This indicates that up-regulated miRNAs may be involved in cell cycle control by ubiquitination system.
The KEGG pathway analysis of the target genes revealed that the down-regulated miRNAs were involved in tuning of several signalling pathways that are known to be disrupted in diseases, and the up-regulated miRNAs were related to remodelling of extracellular matrix organization, packaging, circadian entrainment, recycling and modification of receptors, proteins, chemokines and enzymes in favour of disease progression (Table 2).

2.6. miRNA–Transcription Factor Interaction Network between Dog and Human Melanoma

Previous human melanoma studies indicate that several differentially expressed miRNAs have similar trend of expression with COM [22,23,24,25,26]. Transcription factors like MITF can play a critical role in melanoma [27,28]. So, we were interested to find a common miRNA–TFs co-regulatory network in human and dog for melanoma. Moreover, miRNAs those can target maximum number of TFs will be considered to be important for melanoma-related TFs regulation. We considered the same seed sequences miRNAs between human and dog from our study and the miRNAs target (genes) orthologues TFs that were differentially expressed in human melanoma (GSE31909) were selected for network construction. A total of 34 up-regulated and 33 down-regulated TFs were obtained from GSE31909 (Table S5). We constructed two networks, one using down-regulated miRNAs and up-regulated TFs, and the other using up-regulated miRNAs and down-regulated TFs (Figure 4a,b). See methods for details. We measured the degree and betweenness centrality of the networks to detect the key miRNAs that can regulate maximum TFs in both groups. Nodes that had higher centrality values than average were considered to influence the network biologically.
In the down-regulated miRNA–TF regulatory network (Figure 4a), the miR-126, miR-183, and miR-200 families, let-7 family members had higher degree centralities than average. Among them, miR-126 had the highest centrality (Table S6a).
In the up-regulated miRNA–TF regulatory network (Figure 4b), miR-130 family, and miR-9, miR-20b, miR-371, miR-106a, miR-450b, miR-21, and miR-424, had higher degree centrality than average. Among the miRNAs, miR-20b and miR-106a had the highest centralities (Table S6b). So it indicates that among the differentially expressed miRNAs, miR-126, miR-20b and miR-106a are the most potent to regulate the maximum number of TFs in melanoma.

2.7. Differential miRNA Chromosomal Enrichments

Studies reported that more than 50% of miRNA genes are encoded in the cancer-associated regions or fragile sites or chromosomal breakpoints which are frequently absent, amplified or rearranged in tumour specimens[29,30]. Thereby, the dysregulation of miRNA expression occurred in tumour. For example, miR-15a/16-1 is located in genomic locus containing tumour suppressor that is frequently deleted in leukaemia [31]. So to investigate the melanoma-associated regions or fragile sites it’s important to understand in which chromosome the differential miRNAs reside. We analyzed the chromosomal locations of all 542 annotated miRNAs from miRBase. Stem-loop or mature sequences were mapped against the dog genome to obtain locations for the hsa-miRs (miRNAs that were annotated by the human sequence). Among the 542 miRNAs, 70 (12.84%) were in the X chromosome, and 14 (2.57%) and 24 (4.40%) were in chromosomes 4 and 5, respectively. Three chromosomes were focused because most of differentially expressed miRNAs were encoded in these chromosomes (Figure 5a,b). Out of 30 up-regulated miRNAs, 11 (34.4%) were in X, an about 2.8-fold significant enrichment (p = 5.6 × 10−4). Among the down-regulated miRNAs, four in chromosome 4 (p = 0.008, enrichment = 4.55-fold) and six in chromosome 5 (p = 0.001, enrichment = 4.55-fold) were enriched.
Among the 11 up-regulated miRNAs in Chromosome X, nine miRNAs encode from two clusters: mir-106a/18b/20b/19b-2/92a-2/363 and mir-424/503/542/450a-2/450a-1/450b. Other members in the clusters (miR-19b-2, miR-92a-2, and miR-503) are not listed among the up-regulated miRNA list because they did not meet our stringent criteria, but the changes in their expression were similar to other members of the clusters, except miR-92a-2. Two miRNAs, miR-223 and miR-421, are encoded separately as single genes.
Among the four miRNAs in Chromosome 4 two from miR-143/145 cluster, and rest miR-1271 and miR-1260a encode as a single gene. Rest six down-regulated miRNAs from chromosome 5 encode from two clusters, the miR-200 family and miR-99a-2/let-7a-2, and miR-101 as a single gene. All the family miRNAs had similar expression patterns.
These results are consistent with that study, found clustered miRNAs stay and act together [32]. It also indicates that in chromosome X, 4 and 5 most of the differentially expressed miRNAs were from clusters and their other cluster members also have a similar trend of expression. This may recognize COM specific chromosomal fragile sites.

3. Discussion

Despite dogs being considered as a preclinical model for human melanoma [6], until now, the global miRNA profile was not fully revealed. In this study, we comprehensively analyzed the miRNA profile in COM. The expression levels of miRNAs studied previously [10,11] and our recently reported oncogenic miRNA [15] expression were consistent with those of the present study. Moreover, we detected several differentially expressed miRNAs that have not been reported previously (Table S2).
Some of the differentially expressed miRNAs (up-regulated miR-301a, 130, 383, 21, 454, 335, 132, 423,424, 146b, 9, 20b and down-regulated let-7a, 7b miR-126, 125a, 183, 26b, 29c, 152, 31, 145, 141, 205, 203, 200, 101) were reported in human melanoma [22,23,24,25,26]. The expression trends of these miRNAs correlated well between human melanoma and COM. This indicates an overlap of miRNomes between the species and can be used as a model for human miRNA therapeutics development. It also affirms that dogs share much of their ancestral DNA with humans [33]. To further understand the functions of the miRNAs compared to humans, the targets of miR-450b (PAX9), miR-301a (NDRG2), and miR-223 (ACVR2A) which were reported previously [16,17,18] in human were analyzed. We found, PAX9 and NDRG2, which were down-regulated in human and mouse melanoma, also were down-regulated in COM [34,35]. The expression of the miR-450b–PAX9 and miR-301a–NDRG2 pairs was significantly negatively correlated, which supports two studies that reported miR-450b and miR-301a can bind and suppress PAX9 and NDRG2 activity, respectively [16,17]. ACVR2A is reported here for the first time in melanoma with significant negative correlation with miR-223. These results suggest that the tumour suppressive function of PAX9, NDRG2, and ACVR2a were disrupted by miR-450b, miR-301a, and miR-223, respectively, to maintain oncogenic characteristics in COM as like the human studies. Moreover, as previous study reported that PAX9 is down-regulated and miR-301a, 223 is up-regulated in human melanoma [22,35,36]. It also indicates that there are similarities between human and canine melanoma in respect of PAX9, miR-301a and 223 expressions.
The predicted binding site of miR-450b-PAX9 is conserved between human and dog (Figure S3a). BMP4, a downstream of PAX9, was suppressed when miR-450b degraded the function of PAX9, resulting in an increase in MMP9. Previous studies also showed that suppression of PAX9 decreased BMP4 expression and subsequently increased MMP9 [19,37]. Our study revealed that, in COM, this axis is also maintained and interconnected with high expression of miR-450b. Additionally, high MMP9 expression in metastatic melanoma cells may be maintained by this axis.
Axon guidance, endocytosis, and regulation of actin cytoskeleton and pathways in cancer are common pathways between human and dog melanoma regulated by miRNAs [26]. With canine cutaneous melanoma, only PI3K-Akt signalling, focal adhesion, and ECM-receptor interaction pathways are common [38]. This is not surprising because there are molecular differences between cutaneous and mucosal melanomas, so different pathways are likely to be affected [2,3].
TFs can regulate single or multiple gene expressions, so investigation of miRNAs that influence TFs could be more meaningful. MiR-126 has maximum influence over eight TFs that were up-regulated in melanoma. Although, low miR-126 expression was found to have poor prognostic value in several cancers [39], up-regulated miR-20b and miR-106a influenced 11 and 10 TFs, respectively. The miR-20b seed sequence is similar to that of human miR-17-5p, and miR-106a belongs to the miR-17-92 family. Over-expression of hsa-miR-17 and miR-106a is a good predictor of poor overall survival in several human cancers [40], indicating these miRNAs may be a prognostic marker and also a good therapeutic option in both species.
Chromosomal enrichment showed that the X chromosome harboured up-regulated miRNAs. In human melanoma, X-linked miRNAs are also enriched. Women have consistent advantageous prognosis in melanoma compared with men [41]. However, in mucosal melanoma, the incidence is higher in female [42]. Also, breast cancer has X chromosome-linked differential miRNA enriched in women [43]. To our knowledge, until now, the correlation between X-linked miRNA and poor survival has not been explained. However, in humans, a high number of miRNAs related to cancer located in X chromosome compared to Y [44]. Also study shows 52.5% of human miRNA genes are encoded in the cancer-associated regions or fragile sites or chromosomal breakpoints [30]. As a result miRNAs are frequently absent, amplified or rearranged in tumour specimens [30,31]. So, the speculation that miRNA clusters or family members are co-regulated with melanoma-related genes to achieve a regulatory net outcome on a cell or environment is a reasonable explanation of the enrichment of X chromosome-linked differentially expressed miRNAs. On the other hand, down-regulated miRNAs enriched in chromosome 4 and 5. Previous studies showed that miR-15a/16-1 is down-regulated in leukaemia, due to the deletion of genomic locus containing a putative tumour suppressor-containing region, [31]. Also let-7g/mir-135-1 deletion are reported in varies human malignancies [30]. So, further study regarding the down-regulation of two cluster miR-143/145 and 200 families may be interesting to find putative region in the respective chromosomes.
One drawback in our experiment might be the use of less normal samples in NGS screening. However, we overcome the issue by the qPCR validation of 20 differentially expressed miRNAs within 12 normal and 17 melanoma samples. Moreover, all the normal samples were from healthy laboratory beagle dogs which minimize the limitation of diversity regarding normal samples.
Our study comprehensively established a miRNA profile of COM that has not been previously implicated. We have also shown the significance of miR-450b over-expression in melanoma metastatic cells and future studies are necessary to evaluate the others. Furthermore, we are able to report some melanoma-related miRNAs that are also important in human. Besides, as dog oral melanoma is considered as a typical model for non-UV or TWT melanoma in human our findings will give an insight into the miRNA expression of TWT and mucosal melanoma.

4. Materials and Methods

4.1. Clinical Samples and Canine Melanoma Cell Lines

COM tissue specimens were acquired from tumours excised from dogs that had undergone surgery at the Veterinary Teaching Hospital of Kagoshima University. Informed consents were obtained from dog owners. COM patient’s information is presented in Table S1a. Normal oral tissues were collected from healthy laboratory beagle dogs (age range 8–10 years) at Kagoshima University. Experimental conditions and design were approved by Kagoshima University and Veterinary Teaching Hospital ethics committee (KV004; 18.04.2011). All experimental methods were carried out in accordance with the approved guidelines and regulation.
Tissue samples were collected immediately after excision from dogs that had undergone surgery. The diagnosis was confirmed histopathologically by the hospital. The tissue specimens were placed in RNAlater (AM7021, Invitrogen, Carlsbad, CA, USA) immediately after isolation and stored at −80 °C after overnight incubation at 4 °C.
Dog melanoma cell lines KMEC and LMEC were stored in freezing medium (039-23511, CultureSure, Fujifilm Wako Pure Chemical Corporation, Osaka, Japan). Cell lines were cultured according to the procedure described previously [21]. Cells were grown until confluence and then RNA was extracted for evaluation.

4.2. RNA Extraction and Sequencing

A mirVana RNA Isolation kit (AM1560, Thermo Fisher Scientific, Waltham, MA, USA) was used for RNA isolation according to the Manufacture’s standard protocol. The total RNA concentration was measured using a NanoDrop 200c spectrophotometer (ND2000C, Thermo Fisher Scientific). The RNA quality and integrity were assessed with an Agilent 2100 Bioanalyzer (G2939BA, Agilent Technologies, Santa Clara, CA, USA). The RNA Integrity Number (RIN) mean value was 8.8 (range 7–10) for tissue samples and 9.9 (range 9.6–10) for the KMEC and LMEC cell lines.
Following RNA isolation and quality measurement, samples were sequenced by the Hokkaido System Sciences Company (Hokkaido, Japan). Briefly, small RNA (sRNA) libraries were constructed using 1 µg of total RNA with the TruSeq Small RNA Library Preparation kit (Illumina, San Diego, CA) following the manufacturer’s protocol. After obtaining the sRNAs (18–30 nt) from the total RNA, 5′ and 3′ adaptors were ligated to the sRNAs. Then, reverse transcription followed by amplification was performed to create cDNA constructs. A gel purification step was applied to purify the amplified cDNA constructs for cluster generation and Illumina/Hiseq2500 sequencing analysis by the Hokkaido System Science Co., Ltd. (Hokkaido, Japan). The high-quality cleaned reads that we obtained from the company are shown in Table S1b (Phred score > 34). The raw sequences have been submitted to NCBI sequence read archive (SRA) database under accession number PRJNA516252

4.3. Bioinformatics Analysis of Small RNA Reads

The RNA sequencing data were imported into the CLC Genomics Workbench (CLC Bio, Qiagen, Germany) as recommended in the manufacturer’s manual (http://resources.qiagenbioinformatics.com). Normalization of reads, quality, ambiguity, and adapter trimming as well as quality control was performed using the CLC Genomics Workbench (versions 9 and 10). Briefly, the sequencing generated about 103 and 266 million reads from the normal and melanoma libraries, respectively, with single-end reads (Table S1b). We performed a two-step trimming process to remove adapters and other contaminants. In step one, we aimed to remove low quality, ambiguous nucleotides, 3′ adapters, and short (>15 nt) and long reads (>29 nt). In step two, we removed contaminated sequences, 5′ adapters, and the Illumina stop oligo sequence (5′-GAATTCCACCACGTTCCCGTGG-3′). Finally, we obtained about 51 and 142 million reads from the normal and melanoma libraries, respectively, for further analysis of the small RNA reads (Table S1b). Clean reads were analyzed according to the small RNA analysis guideline of the CLC Genomics Workbench. Briefly, the CLC Genomics Workbench was used to extract and count the small RNA from the clean reads and then compare them to databases of miRNAs and other small RNA databases for annotation. Sequence/fragment counts were used as the expression values for the miRNAs/small RNAs in the libraries.
To annotate the small RNA other than miRNA, CLC bio uses two other reference databases (Canis_familiris.canfam3.1.ncrna and Homo_sapiens.GRCh37.ncrna) from ensemble to annotate sequences that had no matches in miRBase [45]. Differential expression between the two groups was followed the EDGE (empirical analysis of differential gene expression) analysis within the CLC bio.

4.4. Edge Analysis

EDGE follows the exact test developed by Robinson and Smyth for two-group comparisons [46]. The exact test counts data that follow a negative binomial distribution and compares the counts in one set of count samples against the counts in another set of count samples. The variability of each group also is taken into account. The original count data are used because the algorithm assumes that the counts on which it operates are negative binomially distributed. We used the default parameters throughout the analysis. Fold change was calculated from the estimated average count per million (cpm) from each group. The estimated average cpm is derived internally in the exact test of the algorithm. Fold change indicates the difference in average cpm values between the groups. The FDR is based on the p-value of the exact test.

4.5. Expression Analysis by qPCR

To measure the expressions of miRNAs and mRNAs by qPCR we used TaqMan microRNA and gene expression assays (Thermo Fisher Scientific). Total RNA (2 ng) was reverse transcribed to cDNA using a TaqMan MicroRNA Reverse Transcription kit (4366597, Thermo Fisher Scientific) according to the manufacturer’s protocol. qPCR was performed using a TaqMan First Advanced Master Mix kit and a one-step plus real-time PCR system (Thermo Fisher Scientific). Thermal cycling was performed according to the manufacturer’s instructions. All experiments were performed in duplicate. The Cq values of RNU6B in the normal and melanoma samples were consistent between the groups, so RNU6B was used as an internal control to calculate miRNA expression. ΔCq was calculated by subtracting the Cq values of RNU6B from the Cq value of the target miRNA. ΔΔCq was calculated by subtracting the mean target miRNA ΔCq value from the ΔCq value of the normal and melanoma samples. The expression level was evaluated using the 2−∆∆Cq method [47]. qPCR reactions of undetermined Cq were assigned Cq = 36 cycle. TaqMan MicroRNA assays were used in this study can quantitate mature miRNAs. MiRNA primer IDs were as follows: RNU6B (ID: 001093), miR-450b (ID: 006407), miR-301a (ID: 000528), miR-140 (ID: 007661), miR-383 (ID: 000573), miR-542 (ID: 001284), miR-223 (ID: 000526), miR-190 (ID: 000489), miR-21 (ID: 000397), miR-122 (ID: 002245), miR-34a (ID: 000426), miR-338 (ID: 000548), miR-429 (ID: 001077), miR-96 (ID: 000186), miR-375 (ID: 000564), miR-183 (ID: 002269), miR-182 (ID: 002334), miR-1 (ID: 000385), miR-143 (ID: 002249), miR-200b (ID: 002251), and miR-141 (ID: 000463).
For the target gene mRNAs, 250 ng RNA was reverse transcribed to cDNA using ReverTra Ace qPCR RT master mix with gDNA Remover (FSQ-301, Toyobo, Osaka, Osaka Prefecture, Japan). The qPCR procedure was the same as that used for the miRNA experiments. The 2−∆∆Cq method also was used to calculate the expression. GAPDH was used as an internal control. The TaqMan gene expression assay was used in the experiments. The gene IDs were as follows: GAPDH (ID: Cf04419463_gH), PAX9 (ID: Cf02705737_m1), MMP9 (ID: Cf02621845_m1), BMP4 (ID: Cf01041266), NDRG2 (ID: Cf02631635_m1), and ACVR2A (ID: Cf02664427_m1).

4.6. Gene Ontology, Pathway Analysis and Network Construction

Gene Ontology and pathway analysis of miRNA target genes were done using the Database for Annotation, Visualization and Integrated Discovery (DAVID) [48]. A common miRNA–TF interaction network was constructed between human and dog by analyzing the differentially expressed TFs from GSE31909. Briefly, we used TargetScan 7.2 [49] and miRDB [50] to predict miRNA targets. The common target genes between the two predictions were considered as targets for the respective miRNAs. A low FDR was considered to indicate a strong relationship between the annotation and the gene.
To construct a common miRNA–mRNA interaction network between humans and dogs we analyzed the BioProject GSE31909 datasets using the GEO2R tool (https://www.ncbi.nlm.nih.gov/geo/info/geo2r.html#background) to get the differentially expressed genes in human melanoma. We picked the target genes of the differentially expressed miRNAs from the differentially expressed genes in GSE31909. From the differentially expressed target genes, we only considered the TFs for network construction. We also considered the miRNAs that had the same seed sequences as the orthologous human miRNA sequences. The MSigDb gene families (http://software.broadinstitute.org/gsea/msigdb/index.jsp) were used to select the transcription factors (TFs) from the miRNA target genes. Since the expression of miRNAs and their targets are inversed, we built miRNA–TF co-regulatory networks with the inversely expressed TFs using Cytoscape v3.5 (http://www.cytoscape.org/) [51]. That means we built two networks, one with down-regulated miRNAs and up-regulated TFs and the other with up-regulated miRNAs and down-regulated TFs. Since TFs can regulate each other, we also performed a STRING v.10.5 (confidence score 0.700) (http://string-db.org/) [52] network analysis within each group of TFs. The final miRNA–TF regulatory network was obtained after merging the STRING TFs network with the respective miRNA–TF regulatory network in Cytoscape. The degree and betweenness of the network were measured using CentiScaPe 2.2 [53].

4.7. Statistical Analysis

We used GraphPad Prism 7 (www.graphpad.com) for the statistical analysis. Hierarchical clustering analysis was performed using the log10 value that was converted from the expression value of every miRNA from each sample. Unsupervised hierarchical clustering was done with Euclidean distance metric and complete linkage. Comparison of the qPCR data was done using Mann-Whitney test followed by Tukey’s test where appropriate. p values < 0.05 were considered significant. Correlation analysis was performed using Spearman’s correlation coefficient. For the miRNA chromosomal enrichment analysis, we used hyper-geometric test.

5. Conclusions

To the best of our knowledge, this study is the first report of comprehensively studied global miRNA profile about COM. Important miRNAs with respect to human melanoma was also explored in this study. As dogs are considered as models for human melanoma, further study will better explain the pathogenesis of melanoma in both species. Also key therapeutic option may reveal by the in-depth future study.

Supplementary Materials

Supplementary materials can be found at https://www.mdpi.com/1422-0067/20/19/4832/s1.

Author Contributions

Conceptualization, N.M. (Naoki Miura) and M.M.R.; methodology, N.M. (Naoki Miura), M.M.R., and Y.-C.L.; validation, N.M. (Naoki Miura) and M.M.R.; formal analysis, N.M. (Naoki Miura), M.M.R., Y.-C.L., A.A.H., and H.-w.C.; investigation, N.M. (Naoki Miura), M.M.R., Y.-C.L, A.A.H., H.-w.C., H.K., N.M. (Naoki Miura), Y.T., T.N., and R.F.; resources, N.M. (Naoki Miura), Y.T., H.K., N.M. (Noriaki Miyoshi), T.N., and R.F., data curation, N.M. (Naoki Miura) and M.M.R.; writing—original draft preparation, M.M.R.; writing—review and editing, N.M. (Naoki Miura), M.M.R, Y.T., H.K., N.M. (Noriaki Miyoshi), T.N., and R.F.; supervision, N.M. (Naoki Miura); project administration, N.M. (Naoki Miura); funding acquisition, N.M. (Naoki Miura).

Funding

This study was supported by the Japan Society for the Promotion of Science KAKENHI (grant nos. 17H03926, 15H14872, 25292180, and 22780283).

Acknowledgements

We thank Ayako Masuda for their help with the experiments. We also thank Margaret Biswas, from Edanz Group (www.edanzediting.com/ac) for editing a draft of this manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

MiRNAsMicroRNAs
COMCanine oral melanoma
qPCRQuantitative real-time PCR
NGSNext generation Sequencing

References

  1. Siegel, R.L.; Miller, K.D.; Jemal, A. Cancer Statistics, 2018. CA Cancer J. Clin. 2018, 67, 7–30. [Google Scholar] [CrossRef] [PubMed]
  2. Akbani, R.; Akdemir, K.C.; Aksoy, B.A.; Albert, M.; Ally, A.; Amin, S.B.; Arachchi, H.; Arora, A.; Auman, J.T.; Ayala, B.; et al. Genomic classification of cutaneous melanoma. Cell 2015, 161, 1681–1696. [Google Scholar] [CrossRef] [PubMed]
  3. Hayward, N.K.; Wilmott, J.S.; Waddell, N.; Johansson, P.A.; Field, M.A.; Nones, K.; Patch, A.M.; Kakavand, H.; Alexandrov, L.B.; Burke, H.; et al. Whole-genome landscapes of major melanoma subtypes. Nature 2017, 545, 175–180. [Google Scholar] [CrossRef] [PubMed]
  4. Hernandez, B.; Adissu, H.A.; Wei, B.R.; Michael, H.T.; Merlino, G.; Mark Simpson, R. Naturally occurring canine melanoma as a predictive comparative oncology model for human mucosal and other triple wild-type melanomas. Int. J. Mol. Sci. 2018, 19, 394. [Google Scholar] [CrossRef] [PubMed]
  5. Gillard, M.; Cadieu, E.; De Brito, C.; Abadie, J.; Vergier, B.; Devauchelle, P.; Degorce, F.; Dréano, S.; Primot, A.; Dorso, L.; et al. Naturally occurring melanomas in dogs as models for non-UV pathways of human melanomas. Pigment Cell Melanoma Res. 2014, 27, 90–102. [Google Scholar] [CrossRef] [PubMed]
  6. Simpson, R.M.; Bastian, B.C.; Michael, H.T.; Webster, J.D.; Prasad, M.L.; Conway, C.M.; Prieto, V.M.; Gary, J.M.; Goldschmidt, M.H.; Esplin, D.G.; et al. Sporadic naturally occurring melanoma in dogs as a preclinical model for human melanoma. Pigment Cell Melanoma Res. 2014, 27, 37–47. [Google Scholar] [CrossRef] [PubMed]
  7. Mochizuki, H.; Kennedy, K.; Shapiro, S.G.; Breen, M.B. BRAF mutations in canine cancers. PLoS ONE 2015, 10, e0129534. [Google Scholar] [CrossRef] [PubMed]
  8. Chu, P.Y.; Pan, S.L.; Liu, C.H.; Lee, J.; Yeh, L.S.; Liao, A.T. KIT gene exon 11 mutations in canine malignant melanoma. Vet. J. 2013, 196, 226–230. [Google Scholar] [CrossRef]
  9. Rupaimoole, R.; Slack, F.J. MicroRNA therapeutics: Towards a new era for the management of cancer and other diseases. Nat. Rev. Drug Discov. 2017, 16, 203–222. [Google Scholar] [CrossRef]
  10. Noguchi, S.; Mori, T.; Hoshino, Y.; Yamada, N.; Maruo, K.; Akao, Y. MicroRNAs as tumour suppressors in canine and human melanoma cells and as a prognostic factor in canine melanomas. Vet. Comp. Oncol. 2013, 11, 113–123. [Google Scholar] [CrossRef]
  11. Noguchi, S.; MORI, T.; Hoshino, Y.; Nami, Y.; Nakagawa, T.; Sasaki, N.; Akao, Y.; Maruo, K. Comparative study of anti-oncogenic microRNA-145 in canine and human malignant melanoma. J. Vet. Med. Sci. 2012, 74, 1108090601. [Google Scholar] [CrossRef] [PubMed]
  12. Islam, F.; Gopalan, V.; Vider, J.; Wahab, R.; Ebrahimi, F.; Lu, C.T.; Kasem, K.; Lam, A.K.Y. MicroRNA-186-5p overexpression modulates colon cancer growth by repressing the expression of the FAM134B tumour inhibitor. Exp. Cell Res. 2017, 357, 260–270. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  13. Friedrich, M.; Pracht, K.; Mashreghi, M.F.; Jäck, H.M.; Radbruch, A.; Seliger, B. The role of the miR-148/-152 family in physiology and disease. Eur. J. Immunol. 2017, 47, 2026–2038. [Google Scholar] [CrossRef] [Green Version]
  14. Wang, J.; Li, Y.; Ding, M.; Zhang, H.; Xu, X.; Tang, J. Molecular mechanisms and clinical applications of MIR-22 in regulating malignant progression in human cancer (Review). Int. J. Oncol. 2017, 50, 345–355. [Google Scholar] [CrossRef] [PubMed]
  15. Ushio, N.; Rahman, M.M.; Maemura, T.; Lai, Y.C.; Iwanaga, T.; Kawaguchi, H.; Miyoshi, N.; Momoi, Y.; Miura, N. Identification of dysregulated microRNAs in canine malignant melanoma. Oncol. Lett. 2019, 17, 1080–1088. [Google Scholar] [CrossRef] [PubMed]
  16. Sun, M.M.; Li, J.M.F.M.; Guo, L.L.; Xiao, H.T.; Dong, L.; Wang, F.; Huang, F.B.; Cao, D.; Qin, T.; Yin, X.H.; et al. TGF-β1 suppression of microRNA-450b-5p expression: A novel mechanism for blocking myogenic differentiation of rhabdomyosarcoma. Oncogene 2013, 33, 2075. [Google Scholar] [CrossRef]
  17. Guo, Y.J.; Liu, J.X.; Guan, Y.W. Hypoxia induced upregulation of miR-301a/b contributes to increased cell autophagy and viability of prostate cancer cells by targeting NDRG2. Eur. Rev. Med. Pharmacol. Sci. 2016, 20, 101–108. [Google Scholar]
  18. Yang, L.; Li, Y.; Wang, X.; Mu, X.; Qin, D.; Huang, W.; Alshahrani, S.; Nieman, M.; Peng, J.; Essandoh, K.; et al. Overexpression of miR-223 tips the balance of pro- and anti-hypertrophic signaling cascades toward physiologic cardiac hypertrophy. J. Biol. Chem. 2016, 291, 15700–15713. [Google Scholar] [CrossRef]
  19. Zhou, J.; Gao, Y.; Lan, Y.; Jia, S.; Jiang, R. Pax9 regulates a molecular network involving Bmp4, Fgf10, Shh signaling and the Osr2 transcription factor to control palate morphogenesis. Development 2013, 140, 4709–4718. [Google Scholar] [CrossRef] [Green Version]
  20. Laulan, N.B.; St-Pierre, Y. Bone morphogenetic protein 4 (BMP-4) and epidermal growth factor (EGF) inhibit metalloproteinase-9 (MMP-9) expression in cancer cells. Oncoscience 2015, 2, 309–316. [Google Scholar] [CrossRef]
  21. Inoue, K.; Ohashi, E.; Kadosawa, T.; Hong, S.-H.; Matsunaga, S.; Mochizuki, M.; Nishimura, R.; Sasaki, N. Establishment and characterization of four canine melanoma cell Lines. J. Vet. Med. Sci. 2004, 66, 1437–1440. [Google Scholar] [CrossRef] [PubMed]
  22. Cui, L.; Li, Y.; Lv, X.; Li, J.; Wang, X.; Lei, Z.; Li, X. Expression of MicroRNA-301a and its functional roles in malignant melanoma. Cell Physiol. Biochem. 2016, 40, 230–244. [Google Scholar] [CrossRef] [PubMed]
  23. Larsen, A.C.; Mikkelsen, L.H.; Borup, R.; Kiss, K.; Toft, P.B.; Von Buchwald, C.; Coupland, S.E.; Prause, J.U.; Heegaard, S. MicroRNA expression profile in conjunctival melanoma. Investig. Ophthalmol. Vis. Sci. 2016, 57, 4205–4212. [Google Scholar] [CrossRef] [PubMed]
  24. Mirzaei, H.H.R.H.; Gholamin, S.; Shahidsales, S.; Sahebkar, A.; Jafaari, M.R.; Mirzaei, H.H.R.H.; Hassanian, S.M.; Avan, A. MicroRNAs as potential diagnostic and prognostic biomarkers in melanoma. Eur. J. Cancer 2016, 53, 25–32. [Google Scholar] [CrossRef] [PubMed]
  25. Xu, Y.; Brenn, T.; Brown, E.R.S.; Doherty, V.; Melton, D.W. Differential expression of microRNAs during melanoma progression: MiR-200c, miR-205 and miR-211 are downregulated in melanoma and act as tumour suppressors. Br. J. Cancer 2012, 106, 553–561. [Google Scholar] [CrossRef] [PubMed]
  26. Fattore, L.; Costantini, S.; Malpicci, D.; Ruggiero, C.F.; Ascierto, P.A.; Croce, C.M.; Mancini, R.; Ciliberto, G. MicroRNAs in melanoma development and resistance to target therapy. Oncotarget 2015, 8, 22262–22278. [Google Scholar] [CrossRef] [PubMed]
  27. Bai, X.; Fisher, D.E.; Flaherty, K.T. Cell-state dynamics and therapeutic resistance in melanoma from the perspective of MITF and IFNγ pathways. Nat. Rev. Clin. Oncol. 2019, 16, 549–562. [Google Scholar] [CrossRef] [PubMed]
  28. Garraway, L.A.; Widlund, H.R.; Rubin, M.A.; Getz, G.; Berger, A.J.; Ramaswamy, S.; Beroukhim, R.; Milner, D.A.; Granter, S.R.; Du, J.; et al. Integrative genomic analyses identify MITF as a lineage survival oncogene amplified in malignant melanoma. Nature 2005, 436, 117–122. [Google Scholar] [CrossRef] [PubMed]
  29. Laganà, A.; Russo, F.; Sismeiro, C.; Giugno, R.; Pulvirenti, A.; Ferro, A. Variability in the incidence of miRNAs and genes in fragile sites and the role of repeats and CpG islands in the distribution of genetic material. PLoS ONE 2010, 5, e11166. [Google Scholar] [CrossRef] [PubMed]
  30. Calin, G.A.; Sevignani, C.; Dumitru, C.D.; Hyslop, T.; Noch, E.; Yendamuri, S.; Shimizu, M.; Rattan, S.; Bullrich, F.; Negrini, M.; et al. Human microRNA genes are frequently located at fragile sites and genomic regions involved in cancers. Proc. Natl. Acad. Sci. USA 2004, 101, 2999–3004. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  31. Calin, G.A.; Dumitru, C.D.; Shimizu, M.; Bichi, R.; Zupo, S.; Noch, E.; Aldler, H.; Rattan, S.; Keating, M.; Rai, K.; et al. Frequent deletions and down-regulation of micro-RNA genes miR15 and miR16 at 13q14 in chronic lymphocytic leukemia. Proc. Natl. Acad. Sci. USA 2002, 99, 15524–15529. [Google Scholar] [CrossRef] [PubMed]
  32. Dambal, S.; Shah, M.; Mihelich, B.; Nonn, L. The microRNA-183 cluster: The family that plays together stays together. Nucleic Acids Res. 2015, 43, 7173–7188. [Google Scholar] [CrossRef] [PubMed]
  33. Lindblad-Toh, K.; Wade, C.M.; Mikkelsen, T.S.; Karlsson, E.K.; Jaffe, D.B.; Kamal, M.; Clamp, M.; Chang, J.L.; Kulbokas, E.J.; Zody, M.C.; et al. Genome sequence, comparative analysis and haplotype structure of the domestic dog. Nature 2005, 438, 803–819. [Google Scholar] [CrossRef] [PubMed]
  34. Kim, A.; Yang, Y.; Lee, M.S.; Yoo, Y.D.; Lee, H.G.; Lim, J.S. NDRG2 gene expression in B16F10 melanoma cells restrains melanogenesis via inhibition of Mitf expression. Pigment Cell Melanoma Res. 2008, 21, 653–664. [Google Scholar] [CrossRef] [PubMed]
  35. Hata, S.; Hamada, J.-I.; Maeda, K.; Murai, T.; Tada, M.; Furukawa, H.; Tsutsumida, A.; Saito, A.; Yamamoto, Y.; Moriuchi, T. PAX4 has the potential to function as a tumor suppressor in human melanoma. Int. J. Oncol. 2008, 33, 1065–1071. [Google Scholar]
  36. Sand, M.; Skrygan, M.; Sand, D.; Georgas, D.; Gambichler, T.; Hahn, S.A.; Altmeyer, P.; Bechara, F.G. Comparative microarray analysis of microRNA expression profiles in primary cutaneous malignant melanoma, cutaneous malignant melanoma metastases, and benign melanocytic nevi. Cell Tissue Res. 2013, 351, 85–98. [Google Scholar] [CrossRef] [PubMed]
  37. Shon, S.K.; Kim, A.; Kim, J.Y.; Kim, K., II; Yang, Y.; Lim, J.S. Bone morphogenetic protein-4 induced by NDRG2 expression inhibits MMP-9 activity in breast cancer cells. Biochem. Biophys. Res. Commun. 2009, 385, 198–203. [Google Scholar] [CrossRef] [PubMed]
  38. Brachelente, C.; Cappelli, K.; Capomaccio, S.; Porcellato, I.; Silvestri, S.; Bongiovanni, L.; De Maria, R.; Verini Supplizi, A.; Mechelli, L.; Sforna, M. Transcriptome Analysis of Canine Cutaneous Melanoma and Melanocytoma Reveals a Modulation of Genes Regulating Extracellular Matrix Metabolism and Cell Cycle. Sci. Rep. 2017, 7, 6386. [Google Scholar] [CrossRef]
  39. Dong, Y.; Fu, C.; Guan, H.; Zhang, Z.; Zhou, T.; Li, B.S. Prognostic significance of miR-126 in various cancers: A meta-analysis. Onco Targets Ther. 2016, 9, 2547–2555. [Google Scholar] [CrossRef]
  40. Liu, F.; Zhang, F.; Li, X.; Liu, Q.; Liu, W.; Song, P.; Qiu, Z.; Dong, Y.; Xiang, H. Prognostic role of miR-17-92 family in human cancers: Evaluation of multiple prognostic outcomes. Oncotarget 2017, 8, 69125–69138. [Google Scholar] [CrossRef]
  41. Scoggins, C.R.; Ross, M.I.; Reintgen, D.S.; Noyes, R.D.; Goydos, J.S.; Beitsch, P.D.; Urist, M.M.; Ariyan, S.; Sussman, J.J.; Edwards, M.J.; et al. Gender-related differences in outcome for melanoma patients. Ann. Surg. 2006, 243, 693–700. [Google Scholar] [CrossRef] [PubMed]
  42. Mihajlovic, M.; Vlajkovic, S.; Jovanovic, P.; Stefanovic, V. Primary mucosal melanomas: A comprehensive review. Int. J. Clin. Exp. Pathol. 2012, 5, 739–753. [Google Scholar] [PubMed]
  43. Muñoz-Rodríguez, J.L.; Vrba, L.; Futscher, B.W.; Hu, C.; Komenaka, I.K.; Meza-Montenegro, M.M.; Gutierrez-Millan, L.E.; Daneri-Navarro, A.; Thompson, P.A.; Martinez, M.E. Differentially expressed MicroRNAs in postpartum breast cancer in Hispanic women. PLoS ONE 2015, 10, e0124340. [Google Scholar] [CrossRef] [PubMed]
  44. Ghorai, A.; Ghosh, U. miRNA gene counts in chromosomes vary widely in a species and biogenesis of miRNA largely depends on transcription or post-transcriptional processing of coding genes. Front. Genet. 2014, 5, 100. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  45. Kozomara, A.; Griffiths-Jones, S. MiRBase: Annotating high confidence microRNAs using deep sequencing data. Nucleic Acids Res. 2014, 42, D68–D73. [Google Scholar] [CrossRef] [PubMed]
  46. Robinson, M.D.; Smyth, G.K. Small-sample estimation of negative binomial dispersion, with applications to SAGE data. Biostatistics 2008, 9, 321–332. [Google Scholar] [CrossRef] [PubMed]
  47. Schmittgen, T.D.; Livak, K.J. Analyzing real-time PCR data by the comparative CT method. Nat. Protoc. 2008, 3, 1103–1108. [Google Scholar] [CrossRef]
  48. Huang, D.W.; Sherman, B.T.; Tan, Q.; Kir, J.; Liu, D.; Bryant, D.; Guo, Y.; Stephens, R.; Baseler, M.W.; Lane, H.C.; et al. DAVID Bioinformatics Resources: Expanded annotation database and novel algorithms to better extract biology from large gene lists. Nucleic Acids Res. 2007, 35, 169–175. [Google Scholar] [CrossRef] [PubMed]
  49. Lewis, B.P.; Burge, C.B.; Bartel, D.P. Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are microRNA targets. Cell 2005, 120, 15–20. [Google Scholar] [CrossRef]
  50. Wong, N.; Wang, X. miRDB: An online resource for microRNA target prediction and functional annotations. Nucleic Acids Res. 2015, 43, 146–152. [Google Scholar] [CrossRef]
  51. Lopes, C.T.; Franz, M.; Kazi, F.; Donaldson, S.L.; Morris, Q.; Bader, G.D.; Dopazo, J. Cytoscape Web: An interactive web-based network browser. Bioinformatics 2011, 27, 2347–2348. [Google Scholar] [CrossRef] [PubMed]
  52. Szklarczyk, D.; Franceschini, A.; Wyder, S.; Forslund, K.; Heller, D.; Huerta-Cepas, J.; Simonovic, M.; Roth, A.; Santos, A.; Tsafou, K.P.; et al. STRING v10: Protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 2015, 43, D447–D452. [Google Scholar] [CrossRef] [PubMed]
  53. Scardoni, G.; Petterlini, M.; Laudanna, C. Analyzing biological network parameters with CentiScaPe. Bioinformatics 2009, 25, 2857–2859. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Figure 1. Profile of small RNA reads in canine oral melanoma by next-generation sequencing: (a) Percentages of the clean reads annotated under the different small RNA categories. Normal (n = 3), Melanoma (n = 8), (b) The miRNA reads were mapped to the canine genome (Canfam3.1) and the mapped reads ratio distributed in each chromosome was analyzed. Multiple t-test (many t-tests at once-one per row), * p < 0.05, ** p < 0.01, (c) Venn diagram showing the total numbers of identified and differentially expressed miRNAs in melanoma. First Venn (blue) indicates the total number of miRNAs identified in both groups. Other three Venn indicates the number of differentially expressed miRNAs (up (+) and down-regulated (−)). fc, fold change, (d) Top 20 highly expressed miRNAs in the normal and melanoma. Twelve miRNAs were common between the groups and rests eight were left blank in melanoma. In the bar graph, the number in each cell represents the rank of the miRNA with respect to the normal group. Arrows represent the changed position of miRNAs in melanoma group, Red arrows indicate the positions of significant deferentially expressed miRs. miRNA/s (miR/s), * statistically significant differentially expressed miRs between normal and melanoma.
Figure 1. Profile of small RNA reads in canine oral melanoma by next-generation sequencing: (a) Percentages of the clean reads annotated under the different small RNA categories. Normal (n = 3), Melanoma (n = 8), (b) The miRNA reads were mapped to the canine genome (Canfam3.1) and the mapped reads ratio distributed in each chromosome was analyzed. Multiple t-test (many t-tests at once-one per row), * p < 0.05, ** p < 0.01, (c) Venn diagram showing the total numbers of identified and differentially expressed miRNAs in melanoma. First Venn (blue) indicates the total number of miRNAs identified in both groups. Other three Venn indicates the number of differentially expressed miRNAs (up (+) and down-regulated (−)). fc, fold change, (d) Top 20 highly expressed miRNAs in the normal and melanoma. Twelve miRNAs were common between the groups and rests eight were left blank in melanoma. In the bar graph, the number in each cell represents the rank of the miRNA with respect to the normal group. Arrows represent the changed position of miRNAs in melanoma group, Red arrows indicate the positions of significant deferentially expressed miRs. miRNA/s (miR/s), * statistically significant differentially expressed miRs between normal and melanoma.
Ijms 20 04832 g001
Figure 2. Differential miRNA expression in COM: (a) Unsupervised euclidean hierarchical clustering by the miRNA normalized expression values in the normal and melanoma libraries, (b) Heatmap visualizes the expression of statistically significant differentially expressed miRNAs in the normal and melanoma libraries. The colour scale (upper right) indicates the expression values. Up- and down-regulated miRNAs are shown in red to green, respectively. Colour saturation indicates the deviation from the median, (c) Relative expression of oncogenic miRNAs and (d) tumour suppressor miRNAs selected from the next-generation sequencing confirmed by qPCR. The Y-axes indicates the relative miRNA expression levels normalized against RNU6B (normal n = 12, melanoma n = 17, Mann-Whitney test followed by Tukey’s test, * p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001), (e) Correlation of fold change between next-generation sequencing (NGS) and qPCR of the up- and down-regulated miRNAs. Each black circle in the graph represents a miRNA.
Figure 2. Differential miRNA expression in COM: (a) Unsupervised euclidean hierarchical clustering by the miRNA normalized expression values in the normal and melanoma libraries, (b) Heatmap visualizes the expression of statistically significant differentially expressed miRNAs in the normal and melanoma libraries. The colour scale (upper right) indicates the expression values. Up- and down-regulated miRNAs are shown in red to green, respectively. Colour saturation indicates the deviation from the median, (c) Relative expression of oncogenic miRNAs and (d) tumour suppressor miRNAs selected from the next-generation sequencing confirmed by qPCR. The Y-axes indicates the relative miRNA expression levels normalized against RNU6B (normal n = 12, melanoma n = 17, Mann-Whitney test followed by Tukey’s test, * p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001), (e) Correlation of fold change between next-generation sequencing (NGS) and qPCR of the up- and down-regulated miRNAs. Each black circle in the graph represents a miRNA.
Ijms 20 04832 g002aIjms 20 04832 g002b
Figure 3. Gene regulatory function of miR-450b, miR-301a, and miR-223: (a) Relative expression of the target genes PAX9, NDRG2, and ACVR2A. Y-axes indicates the relative mRNA expression normalized against GAPDH, (b) Spearman correlation of the expression of the miR-450b-PAX9, miR-301a-NDRG2, and miR-223-ACVR2A pairs, (c) Relative expression of BMP4 and MMP9 in melanoma tissue samples, (d) Relative expression of miR-450b, PAX9, BMP4, and MMP9 in the canine melanoma cell lines KMEC and LMEC. Mann-Whitney test followed by Tukey’s test, * p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001. (e) Schematic representation of the miR-450b regulatory function. MiR-450b inhibits PAX9 and, as a result, BMP4-MMP9 regulation is disrupted (↑–up-regulation, ⊺- Inhbition).
Figure 3. Gene regulatory function of miR-450b, miR-301a, and miR-223: (a) Relative expression of the target genes PAX9, NDRG2, and ACVR2A. Y-axes indicates the relative mRNA expression normalized against GAPDH, (b) Spearman correlation of the expression of the miR-450b-PAX9, miR-301a-NDRG2, and miR-223-ACVR2A pairs, (c) Relative expression of BMP4 and MMP9 in melanoma tissue samples, (d) Relative expression of miR-450b, PAX9, BMP4, and MMP9 in the canine melanoma cell lines KMEC and LMEC. Mann-Whitney test followed by Tukey’s test, * p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001. (e) Schematic representation of the miR-450b regulatory function. MiR-450b inhibits PAX9 and, as a result, BMP4-MMP9 regulation is disrupted (↑–up-regulation, ⊺- Inhbition).
Ijms 20 04832 g003aIjms 20 04832 g003b
Figure 4. Common miRNA–transcription factor (TF) regulatory network between human and dog: (a) Regulatory network for down-regulated miRNAs and its up-regulated target TFs. MiR-126 influences eight TFs and has the highest centrality, (b) Regulatory network for up-regulated miRNAs and its down-regulated target TFs. MiR-20b and miR-106a have highest centralities. A miRNA (V-shaped) or TF (oval-shaped) is considered a node and line between nodes considered edge. Green and red indicate degree scores less and above than average and saturation shows deviation. Colour gradient of the nodes (miRNAs/mRNAs) indicates the strength of network regulation. Edge width represents edge betweenness. Node’s Edges with highest degree scores are in red.
Figure 4. Common miRNA–transcription factor (TF) regulatory network between human and dog: (a) Regulatory network for down-regulated miRNAs and its up-regulated target TFs. MiR-126 influences eight TFs and has the highest centrality, (b) Regulatory network for up-regulated miRNAs and its down-regulated target TFs. MiR-20b and miR-106a have highest centralities. A miRNA (V-shaped) or TF (oval-shaped) is considered a node and line between nodes considered edge. Green and red indicate degree scores less and above than average and saturation shows deviation. Colour gradient of the nodes (miRNAs/mRNAs) indicates the strength of network regulation. Edge width represents edge betweenness. Node’s Edges with highest degree scores are in red.
Ijms 20 04832 g004
Figure 5. Chromosomal enrichment of differentially expressed miRNAs (miRs): (a) Chromosome enrichment of the down-regulated miRNAs, (b) Chromosome enrichment of the up-regulated miRNAs. Hypergeometric test, * p < 0.05, ** p < 0.01, *** p < 0.001.
Figure 5. Chromosomal enrichment of differentially expressed miRNAs (miRs): (a) Chromosome enrichment of the down-regulated miRNAs, (b) Chromosome enrichment of the up-regulated miRNAs. Hypergeometric test, * p < 0.05, ** p < 0.01, *** p < 0.001.
Ijms 20 04832 g005
Table 1. Gene ontology (GO) functional analysis of the target genes of differentially expressed miRNAs.
Table 1. Gene ontology (GO) functional analysis of the target genes of differentially expressed miRNAs.
Down-Regulated miRNA’s Target Genes Ontology
Biological Process
TermsFE 1FDR 2
GO:0018105~peptidyl-serine phosphorylation 2.3230.001
GO:0045944~positive regulation of transcription from RNA polymerase II promoter 1.4610.025
Cellular Component
GO:0005634~nucleus 1.2771.65 × 10−6
GO:0005654~nucleoplasm 1.3742.49 × 10−5
GO:0005737~cytoplasm 1.2321.09 × 10−4
GO:0005911~cell-cell junction 2.110.068146
Molecular Function
GO:0004702~receptor signaling protein serine/threonine kinase activity 2.9329.98 × 10−4
GO:0005201~extracellular matrix structural constituent 3.3730.002697
GO:0003682~chromatin binding 1.7090.002892
GO:0003714~transcription corepressor activity 2.1040.057312
Up-Regulated miRNA’s Target Genes Ontology
Biological Process
TermsFE 1 FDR 2
GO:0042787~protein ubiquitination involved in ubiquitin-dependent protein catabolic process 2.2080.004
Cellular Component
GO:0005654~nucleoplasm 1.4411.35 × 10−7
GO:0005737~cytoplasm 1.2282.87 × 10−4
GO:0005794~Golgi apparatus 1.5908.03 × 10−4
GO:0005634~nucleus 1.2280.001041
Molecular Function
GO:0061630~ubiquitin protein ligase activity 2.0270.012
GO:0044212~transcription regulatory region DNA binding 2.1270.019
1 Fold enrichment, 2 False discovery rate.
Table 2. KEGG pathway analysis of the target genes of the differentially expressed miRNAs.
Table 2. KEGG pathway analysis of the target genes of the differentially expressed miRNAs.
Down-Regulated miRNAs Target Genes Pathway
TermsFE 1FDR 2
cfa05206:MicroRNAs in cancer2.5226.72 × 10−7
cfa04010:MAPK signaling pathway1.9171.20 × 10−4
cfa04151:PI3K-Akt signaling pathway1.7621.65 × 10−4
cfa04360:Axon guidance2.3074.96 × 10−4
cfa05205:Proteoglycans in cancer1.9818.64 × 10−4
cfa04910:Insulin signaling pathway2.2098.88 × 10−4
cfa04152:AMPK signaling pathway2.2300.003
cfa04510:Focal adhesion1.9010.003
cfa04722:Neurotrophin signaling pathway2.1640.010
cfa04012:ErbB signaling pathway2.3840.013
cfa04512:ECM-receptor interaction2.3840.013
Up-regulated miRNAs target genes pathway
TermsFE 1FDR 2
cfa04144:Endocytosis2.4241.31 × 10−11
cfa04810:Regulation of actin cytoskeleton2.2732.78 × 10−7
cfa05200:Pathways in cancer1.6040.010
cfa04710:Circadian rhythm3.5360.060
cfa05410:Hypertrophic cardiomyopathy (HCM)2.3980.063
cfa05414:Dilated cardiomyopathy2.3460.065
cfa04713:Circadian entrainment2.2100.098
1 Fold enrichment. 2 False discovery rate.

Share and Cite

MDPI and ACS Style

Rahman, M.M.; Lai, Y.-C.; Husna, A.A.; Chen, H.-w.; Tanaka, Y.; Kawaguchi, H.; Miyoshi, N.; Nakagawa, T.; Fukushima, R.; Miura, N. Micro RNA Transcriptome Profile in Canine Oral Melanoma. Int. J. Mol. Sci. 2019, 20, 4832. https://doi.org/10.3390/ijms20194832

AMA Style

Rahman MM, Lai Y-C, Husna AA, Chen H-w, Tanaka Y, Kawaguchi H, Miyoshi N, Nakagawa T, Fukushima R, Miura N. Micro RNA Transcriptome Profile in Canine Oral Melanoma. International Journal of Molecular Sciences. 2019; 20(19):4832. https://doi.org/10.3390/ijms20194832

Chicago/Turabian Style

Rahman, Md. Mahfuzur, Yu-Chang Lai, Al Asmaul Husna, Hui-wen Chen, Yuiko Tanaka, Hiroaki Kawaguchi, Noriaki Miyoshi, Takayuki Nakagawa, Ryuji Fukushima, and Naoki Miura. 2019. "Micro RNA Transcriptome Profile in Canine Oral Melanoma" International Journal of Molecular Sciences 20, no. 19: 4832. https://doi.org/10.3390/ijms20194832

APA Style

Rahman, M. M., Lai, Y. -C., Husna, A. A., Chen, H. -w., Tanaka, Y., Kawaguchi, H., Miyoshi, N., Nakagawa, T., Fukushima, R., & Miura, N. (2019). Micro RNA Transcriptome Profile in Canine Oral Melanoma. International Journal of Molecular Sciences, 20(19), 4832. https://doi.org/10.3390/ijms20194832

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop