Ex) Article Title, Author, Keywords
pISSN 1598-298X
eISSN 2384-0749
Ex) Article Title, Author, Keywords
J Vet Clin 2023; 40(2): 164-174
https://doi.org/10.17555/jvc.2023.40.2.164
Published online April 30, 2023
Eunpyo Kim1 , Giup Jang2 , Jin-Wook Kim3 , Wan-Hee Kim4 , Geon-A Kim5,*
Correspondence to:*kimgeona2020@gmail.com
Copyright © The Korean Society of Veterinary Clinics.
Canine splenic tumors (STs) are commonly diagnosed during imaging examinations, such as in X-ray and ultrasonography examinations, suggesting their higher prevalence, especially in older dogs. Despite this high prevalence, there are no effective treatment options for STs because of the difficulties in determining therapeutic targets. However, recently, the importance of microRNAs (miRNAs) has evolved owing to their ambivalent characteristics. Biomarkers and novel therapies using miRNAs have been well-studied in human cancer research compared to canine research, except for mammary gland tumors. Therefore, this study aimed to comparatively analyze miRNA expression profiles according to malignancy and biopsy sites to identify novel therapeutic and diagnostic targets. Tissue samples were collected directly from splenic tumor masses and immersed in RNAlater solution for further analysis. To investigate differentially expressed genes (DEGs) between tumor and normal tissues, we used RNA-seq and miRNA microarray analysis. Then, functional analysis based on DEGs was conducted to sort tumor-related DEGs. We found that cfa-miR-150 was upregulated in benign tumors, whereas cfa-miR-134 was upregulated in malignant tumors. Despite limited information on canine miRNAs, we identified two potential biomarkers for the differential diagnosis of STs.
Keywords: spleen tumor, microRNA, dog, benign, malignancy, messenger RNA
Splenic tumors (STs) in dogs are commonly diagnosed during an ultrasound examination. Half of them are found to be malignant (MS), whereas the others are diagnosed as benign (BS) (13). Also, small-breed dogs were more likely to be diagnosed with benign STs than large-breed dogs (13). Hemangiosarcoma is a common malignant tumor of the spleen in older dogs, with a prevalence of 12-18% in all mesenchymal neoplasms (3). Since hemangiosarcoma has a poor prognosis even with surgery and median survival time is known to be about 3 months, early diagnosis or development of biomarker is essential (20). Because the spleen functions in blood storage, total splenectomy is often required to prevent splenic rupture in dogs with STs. The exact causes of STs are unknown, although it has been described that large breeds (e.g. Golden Retrievers, German Shepherds, Labradors) are predisposed (13). Although total or partial splenectomy is suggested to prevent further intra-abdominal hemorrhage and coagulopathy syndrome, perioperative death is common in dogs because of the presence of higher risks, such as thrombocytopenia, anemia, and intra-operative arrhythmia (22). In addition, even if STs are clearly diagnosed by radiological examinations, the definitive diagnosis is always based on histopathologic analysis (21). Thus, a general analysis is needed for more precise prediction of patients’ prognosis and determining therapeutic strategies.
miRNAs, endogenous small non-coding RNAs that function as regulators of gene expression, have been suggested to be dysregulated in cancer (15). Since miRNAs are expressed differently according to the tissue type and malignancy, they have been studied as biomarkers for cancer diagnosis owing to their non-invasive properties (6,11). Furthermore, as altered miRNAs can result in the development of drug resistance, they are being studied as therapeutic targets and newer drug candidates (19). miRNAs can act as both tumor suppressors and oncogenes based on their mRNA targets (25). For example, miR-12/13 targets bcl-2 and oncogenes (e.g., Cyclin D1, MCL1, CDC2, ETS1, and JUN), showing tumor suppressor effects (14). However, one study reviewed the oncogenic properties of miRNAs in malignant transformation, suggesting that they are potential biomarkers and therapeutic targets in cancers (7). Thus, to develop novel biomarkers and drugs, it is important to profile miRNA expression profiles in veterinary oncology.
Compared to the numerous studies on canine mammary gland tumors, miRNA profiling of canine STs is an area that requires further studies. Therefore, this study aimed to comparatively analyze miRNA profiles in canine splenic tumors based on their malignancy and biopsy sites using histological examination to understand the interconnection between altered miRNAs and their mechanism of action and identify potent biomarkers or novel therapeutic target candidates.
The canine participants enrolled in this study was not purchased for experimental usage. This study was approved by the Seoul National University Institutional Animal Care and Use Committee (approval number: SNU-200217-3-2). All samples are acquired under the consent from every owner.
Tissue samples were collected during surgical operation at veterinary clinics in Republic of Korea. The biopsy samples were collected at the center of splenic cancer masses. Right after the resection, the samples were directly immersed in RNAlater solution purchased from Invitrogen (Invitrogen, MA, USA) and formalin solution for histopathological examination to evaluate malignancy. The stored samples are then preserved in 4°C refrigerator until use. All patients are diagnosed without other comorbidities at the time of physical examination. The patients’ information is summarized in Supplementary Table 1.
Total RNA was extracted from each tissue samples using RNeasy mini kit (Qiagen, Doncaster, VIC, Australia) according to the manufacturer’s protocol. RNA purity and integrity were evaluated using ND-2000 Spectrophotometer (NanoDrop, Wilmington, USA) and Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, USA), respectively.
Microarray hybridization was performed using the Affymetrix GeneChip miRNA 4.0 Array (ThermoFisher,MA, USA) as per manufacturer’s protocol. Total RNA samples were labelled using FlashTagTM Biotin RNA Labeling Kit (Genisphere, PA, USA). The labelled RNA samples were then quantified, fractionated, and hybridized to the microarray according to the manufacturer’s protocol. RNA-array hybridization was performed on an Affymetrix® 450 Fluidics Station (ThermoFisher, MA, USA). The arrays were stained using a GeneChip Fluidics Station 450 (Affymetrix, Santa Clara, California, United States), and scanned using an Affymetrix GCS 3000 scanner (Affymetrix, Santa Clara, California, USA). miRNA-mRNA hybridization signals were analyzed using the Affymetrix® GeneChipTM Command Console.
Total RNA concentration was calculated using Quant-IT RiboGreen (Invitrogen, #R11490). To assess the integrity of the total RNA, samples were run on the TapeStation RNA ScreenTape (Agilent, #5067-5576). Only high-quality RNA preparations with RNA integrity number (RIN) >7.0 were used for RNA library construction. A library was independently prepared with 1 µg total RNA of each sample using the Illumina TruSeq Stranded mRNA Sample Prep Kit (RS-122-2101, Illumina Inc., San Diego, CA, USA). The first step in the workflow involved purifying the poly-A containing mRNA molecules using poly-T-attached magnetic beads. Following purification, the mRNA samples were fragmented into small pieces using divalent cations under elevated temperature. The cleaved RNA fragments were then used as templates to generate first strand cDNA using SuperScript II reverse transcriptase (Invitrogen, MA, USA) with random primers. This was followed by second strand cDNA synthesis using DNA Polymerase I, RNase H and dUTP. The generated double stranded cDNA fragments were then subjected to an end repair process involving adenylation followed by adapter ligation. The products were then purified and enriched with PCR to create the final cDNA library. The libraries were quantified using KAPA Library Quantification kits for Illumina Sequencing platforms according to the qPCR Quantification Protocol Guide (KAPA BIOSYSTEMS, Wilmington, USA) and qualified using the TapeStation D1000 ScreenTape (Agilent Technologies, CA, USA). Indexed libraries were then submitted to Illumina NovaSeq (Illumina Inc., CA, USA) to perform paired-end (2 × 100 bp) sequencing.
At the miRNA level, raw data were extracted automatically through the Affymetrix data extraction protocol using the Affymetrix GeneChip® Command Console® Software (AGCC). The CEL files import, miRNA level RMA+DABG-All analyses and results were exported using Affymetrix® Power Tools (APT) Software. Array data were filtered by probes annotated species.
At the mRNA level, we preprocessed the raw reads from the sequencer to remove low quality and adapter sequences before analysis and aligned the processed reads to the Canis lupus familiaris (CanFam3.1) database using HISAT v2.1.0 (10). HISAT utilizes two types of indexes for alignment (a global, whole-genome index and tens of thousands of small local indexes). These two types’ indexes are constructed using the same BWT (Burrows–Wheeler transform)/ a graph FM index (GFM) as Bowtie2. Because of its use in efficient data structures and algorithms, HISAT generates spliced alignments several times faster than Bowtie and BWA. The reference genome sequence of Canis lupus fa-miliaris (CanFam3.1) and annotation data were downloaded from the NCBI database. Transcript assembly of known transcripts was then processed by StringTie v2.1.3b (16). Based on the results, expression abundance of transcripts was calculated as read count or FPKM value (Fragments Per Kilobase of exon per Million fragments mapped) per sample. The expression profiles were used for additional analysis such as that of differentially expressed genes (DEGs). In groups with different conditions, DEGs or transcripts can be filtered through statistical hypothesis testing.
Comparative analysis between test and control samples was carried out using in-dependent t-test and fold change with the null hypothesis indicating no difference exists among groups. False discovery rate (FDR) was controlled by adjusting p-value using Benjamini-Hochberg algorithm. All Statistical tests and visualization of DEGs was conducted using R statistical language (v3.3.2).
We collected miRNA-disease associations in the Human microRNA Disease Database (HMDD) to confirm that the DE-miRNAs are associated with disease (9). HMDD is a database for human-related miRNA-disease associations. However, the analysis was performed considering heterologous miRNAs for interspecies. Because information on spleen tumor does not exist in HMDD, miRNAs related to lymphatic disease, which is a higher concept than spleen tumor, were collected.
We collected miRNA-gene associations from miRDB and performed gene level analysis based on miRNA targets (4). It contains information about predicted target genes and predicted score for miRNAs. We used associations with a predictive score of ≥50 for miRNA-gene relationships, because the amount of known information about spleen tumor is less.
clusterProfiler v3.18.1 in R was used to perform functional analysis of gene set (24). Gene Ontology includes biological process, cellular component, and molecular function, and Kyoto Encyclopedia of Genes and Genomes (KEGG) includes biological pathway. To verify the disease-level for the identified genes, we collected disease-related genes from DISEASES and DisGeNet (17,18). To validate the results of functional analysis, we collected disease-related Gene Ontology and biological pathway from the Comparative Toxicogenomics Database (CTD) (5).
Differentially expressed miRNAs (DE-miRNAs) were identified for three comparisons:1) Case_21 vs. Control_21, 2) Case_22 vs. Control_22, and 3) Case_22 vs. Case_21. Meaning of Case_21 and Control_21 are tumor tissue of benign ST and normal tissue of benign ST, respectively. Meaning of Case_22 and Control_22 are tumor tissue of malignant ST and normal tissue of malignant ST. |Fold change| ≥ 2 and p-value < 0.05 were used to identify DE-miRNAs for each comparison. Fig. 1 shows a volcano plot for identifying DE-miRNAs based on these cut-off values.
Twenty-nine miRNAs (14 upregulated and 15 downregulated) were identified in Case_21/Control_21, 13 miRNAs (5 upregulated and 8 downregulated) in Case_22/Control_22, and 8 miRNAs (4 upregulated and 4 downregulated) in Case_22/Case_21 (Supplementary File 1).We validated the relationship between DE-miRNAs and diseases using the HMDD. We confirmed that 14 of the 29 miRNAs in Case_21/Control_21 and 7 of the 13 miRNAs were associated with lymphatic disease (Table 1).
Table 1 Lymphatic disease-associated miRNAs among DE-miRNAs
Case_21/Control_21 | |||
---|---|---|---|
miRNA | Regulation | Fold change | p-value |
cfa-miR-150 | UP | 2.407 | 0.010 |
cfa-miR-21 | UP | 2.221 | 0.018 |
cfa-miR-221 | UP | 2.022 | 0.039 |
cfa-miR-204 | DOWN | –5.731 | 0.035 |
cfa-miR-222 | UP | 2.536 | 0.002 |
cfa-miR-196b | DOWN | –3.235 | 0.007 |
cfa-miR-224 | DOWN | –14.507 | 0.001 |
cfa-miR-206 | DOWN | –4.508 | 0.003 |
cfa-miR-27b | DOWN | –2.329 | 0.002 |
cfa-miR-194 | DOWN | –2.119 | 0.007 |
cfa-miR-504 | DOWN | –3.728 | 0.035 |
cfa-miR-383* | DOWN | –18.774 | 0.006 |
cfa-miR-133b | DOWN | –14.169 | 0.001 |
cfa-miR-139* | DOWN | –4.927 | 0.023 |
Case_22/Control_22 | |||
cfa-miR-181a | DOWN | –2.568 | 0.011 |
cfa-miR-31 | DOWN | –4.324 | 0.036 |
cfa-miR-326 | DOWN | –4.658 | 0.039 |
cfa-miR-212 | UP | 6.353 | 0.029 |
cfa-miR-197 | DOWN | –2.127 | 0.012 |
cfa-miR-383* | DOWN | –35.603 | 0.005 |
cfa-miR-139* | DOWN | –6.411 | 0.002 |
*Indicates miRNAs that appear in both comparisons.
We confirmed that half of the DE-miRNAs detected in the comparisons were associated with lymphatic disease, and cfa-miR-139 and cfa-miR-383 were common miRNAs. In Case_22/Case_21, miRNAs related to lymphatic disease were not detected; however, three downregulated miRNAs were identified in hematologic neoplasms, a disease at the same level as lymphatic disease (cfa-miR-31, cfa-miR-150, and cfa-miR-708). Thus, three miRNAs were highly expressed in Case_21, and it was confirmed that benign and malignant tumors could be distinguished at the miRNA level. A related study confirmed that hematogenous transmission is the most likely route of transmission, causing splenic metastasis [primary and secondary neoplasms of the spleen]. Differences in miRNA levels between benign and malignant tumors are presumed to be closely related to metastasis and cancer development.
We performed gene-level analysis of DE-miRNAs based on miRDB and miRNA-gene associations collected from the miRNA target prediction database. We collected 5,549 target genes in Case_21/Control_21, 3,038 target genes in Case_22/Control_22, and 2,733 target genes in Case_22/Case_21from their DE-miRNAs. Functional analysis of the gene set was performed using clusterProfiler in the R software. The analysis categories used were GO_BP, GO_CC, GO_MF, and KEGG. GO_BP was a biological process, GO_CC was a cellular component, and GO_MF was a molecular function. To validate whether the results analyzed in each category were related to spleen tumors, disease-related Gene Ontology and KEGG were collected from CTD using the MeSH ID of spleen tumors (D013160). Fig. 2 shows an image of a list associated with spleen tumors among significant Gene Ontologies and KEGG derived from clusterProfiler (p-value < 0.05, Supplementary File 2).
We confirmed that there was no information on spleen tumors in GO and KEGG, but many biological functions related to cancer were identified among the significant categories (Fig. 2). Moreover, we confirmed that the biological functions that prevent the accumulation of damaged biomolecules, such as autophagy and mitophagy, were identified among the significant categories of each comparison.
We collected a list of immune-related pathways from related research to confirm categories related to immunity among the spleen-related function analysis results identified in Fig. 2 [Maternal nematode infection upregulates the expression of Th2 Treg and diapedesis-related genes in the neonatal brain]. Twenty-six KEGG pathway related to immunity are presented in this paper. Table 2 shows a list of genes related to immunity among the spleen-related KEGG genes detected in each comparison.
Table 2 Immune-related KEGG pathway list based on DE-miRNA target genes
Case_21/Control_21 | |
---|---|
KEGG pathway | p-value |
T cell receptor signaling pathway | 0.001219192 |
Endocytosis | 2.68E-08 |
Regulation of actin cytoskeleton | 0.000131578 |
Case_22/Control_22 | |
KEGG pathway | p-value |
Fc epsilon RI signaling pathway | 0.00899355 |
T cell receptor signaling pathway | 0.026774021 |
Endocytosis | 3.26E-10 |
Regulation of actin cytoskeleton | 0.003272488 |
Case_22/Case_21 | |
KEGG pathway | p-value |
Chemokine signaling pathway | 0.049509358 |
Fc epsilon RI signaling pathway | 0.04652657 |
T cell receptor signaling pathway | 0.012534506 |
Endocytosis | 4.23E-06 |
Regulation of actin cytoskeleton | 0.001236387 |
We confirmed that the T cell receptor signaling pathway, endocytosis, and regulation of the actin cytoskeleton appeared in common in the three comparisons. The Fc epsilon RI signaling pathway was confirmed in Case_22/Control_22, and fc epsilon RI and chemokine signaling pathways were detected in Case_22/Case_21.In immune-related research on the bone marrow and spleen, it was confirmed that the T cell receptor signaling pathway, chemokine signaling pathway, and fc epsilon RI signaling pathway were the most enriched pathways [extracellular vesicles mediate radiation-induced systemic bystander signals in the bone marrow and spleen].
We collected disease-related genes from DISEASES, DisGeNet, and identified disease-related genes among the target genes. Ten spleen cancer-related genes and 10 spleen disease-related genes were collected from DISEASES, and 6 splenic neoplasm-related genes were collected from DisGeNet. Because known information about spleen tumors is scarce, several databases and diseases have been used. Table 3 shows a list of disease-related genes among miRNA target genes.
Table 3 Disease-related genes for miRNA target genes
Case_21/Control_21 | ||
---|---|---|
miRNA | Target gene | Database |
cfa-miR-379 | NOTCH2 | DISEASES/Spleen cancer |
cfa-miR-27b | ||
cfa-miR-206 | ||
cfa-miR-299 | CD19 | DISEASES/Spleen cancer |
cfa-miR-150 | ITGAE | DISEASES/Spleen cancer |
cfa-miR-1 | CD5 | DISEASES/Spleen cancer DISEASES/Spleen disease |
cfa-miR-206 | ||
cfa-miR-204 | ||
cfa-miR-543 | LMLN | DISEASES/Spleen disease |
cfa-miR-376a | ||
cfa-miR-299 | THPO | DISEASES/Spleen disease |
cfa-miR-133b | CD8A | DISEASES/Spleen disease |
cfa-miR-1841 | CAST | DisGeNet/Splenic neoplasm |
Case_22/Control_22 | ||
miRNA | Target gene | Database |
cfa-miR-330 | CD5 | DISEASES/Spleen disease |
cfa-miR-376a | LMLN | DISEASES/Spleen disease |
cfa-miR-181a | CD4 | DISEASES/Spleen disease |
cfa-miR-31 | TNF | DISEASES/Spleen disease |
Case_22/Case_21 | ||
miRNA | Target gene | Database |
cfa-miR-342 | CD79B | DISEASES/Spleen cancer |
cfa-miR-665 | CD79A | DISEASES/Spleen cancer |
cfa-miR-150 | ITGAE | DISEASES/Spleen cancer |
cfa-miR-379 | NOTCH2 | DISEASES/Spleen cancer |
cfa-miR-27b | ||
cfa-miR-206 | ||
cfa-miR-31 | TNF | DISEASES/Spleen disease |
CD5 and LMLN were identified as genes commonly identified in benign and malignant tumors, and the remaining genes were identified as genes uniquely expressed in benign or malignant tumors (Fig. 3). In addition, CD79A and CD79B were uniquely derived in benign and malignant tumors compared with the differences between normal and tumor areas. Although the number of genes known to be related to the spleen was small, we confirmed that a large number of disease-related genes were detected among genes affected by differentially expressed miRNAs. Differentially expressed genes (DEGs) at the mRNA level were identified using the same comparisons and conditions (|fold change| ≥ 2 and p-value < 0.05). Fig. 3 shows the volcano plot for identifying DEGs based on these cutoff values.
A total of 1,238 genes (307 upregulated and 931 downregulated) were identified in Case_21/Control_21, 904 genes (369 upregulated and 535 downregulated) in Case_22/Control_22, and 1,156 genes (582 upregulated and 574 downregulated) in Case_22/Case_21 (Supplementary File 3).
Functional analysis of the DEG list was performed using clusterProfiler in R. Similar to the miRNA level, disease-related Gene Ontology and KEGG were collected from CTD using MeSH ID of spleen tumor. Fig. 4 shows an image of a list associated with spleen tumors among significant Gene Ontologies and KEGG derived from clusterProfiler (p-value < 0.05, Supplementary File 4).
Similar to the analysis results at the miRNA level, we confirmed the appearance of multiple cancer-related biological functions in the three comparisons shown in Fig. 4. In particular, the cell cycle, T cell receptor signaling pathway, B cell receptor signaling pathway, and apoptosis were also derived from the difference between benign and malignant tumors as significant functions, along with the previous results.
We identified categories related to immunity among the significant categories related to the spleen (Fig. 4). Table 4 shows a list related to immunity among the spleen-related KEGG genes detected in each comparison at the mRNA level.
Table 4 Immune-related KEGG pathway list based on DEGs
Case_21/Control_21 | |
---|---|
KEGG pathway | p-value |
Chemokine signaling pathway | 0.000612405 |
Regulation of actin cytoskeleton | 1.41E-05 |
Case_22/Control_22 | |
KEGG pathway | p-value |
Fc epsilon RI signaling pathway | 0.003332669 |
Regulation of actin cytoskeleton | 0.000374166 |
Case_22/Case_21 | |
KEGG pathway | p-value |
T cell receptor signaling pathway | 2.87E-09 |
B cell receptor signaling pathway | 0.009193225 |
Regulation of actin cytoskeleton | 0.02036557 |
The pathways identified at the miRNA level also appeared identical at the mRNA level (Table 3). At the mRNA level, the B cell receptor signaling pathway was also derived from the comparison between benign and malignant tumors, and it was confirmed that mature B lymphocytes are closely related to the immune response in the spleen (B cell development in the spleen takes place in discrete steps and is determined by the quality of B cell receptor-derived signals).
Based on the analysis of miRNA and mRNA levels, we examined the correlation between them. First, genes that appeared in common among the target genes of DE-miRNAs and the DEG of mRNA were identified. Among the DEGs identified in the mRNA analysis, only genes corresponding to protein coding were used. Fig. 5 presents a Venn diagram of the genes identified at both levels.
We commonly identified 423 genes in Case_21/Control_21, 173 genes in Case_22/Control_22 and 194 genes in Case_22/Case_21 at the miRNA and mRNA levels (Supplementary File 5). ClusterProfiler was used to confirm the functions of the common genes. Unlike the functional analysis conducted previously at the miRNA and mRNA levels, the analysis results for the top five categories based on p-values in each category were included, not limited to the spleen. Fig. 6 shows the results of functional analysis of the common genes of the three comparisons (Supplementary File 6).
The analysis results based on common genes showed different aspects from the previous analysis results at the miRNA and mRNA levels. Among the significant categories, biological functions related to blood vessels and muscles were identified, and in the comparison between benign and malignant tumors, it was confirmed that significant results were also derived in cellular senescence and dendritic-related functions. Based on the results of the previous DE-miRNA target gene-based functional analysis and mRNA DEG-based functional analysis, a Venn diagram was constructed for categories related to the spleen, as shown in Fig. 7.
Fourteen biological functions in Case_21/Control_21, 26 functions in Case_22/Control_22, and 25 functions in Case_22/Case_21 were identified at the miRNA and mRNA levels (Supplementary File 7). We recognized that more diverse biological functions might appear because the number of miRNA target genes is greater than the number of differentially expressed genes at the mRNA level. The number of functions of the gene clusters identified at the miRNA level was high. Nevertheless, we confirmed that most functions of the genes analyzed at the miRNA and mRNA levels were shared.
We performed various comparisons of canine spleen tumors. We derived various insights into normal and tumor tissues for benign tumors, normal and tumor tissues for malignant tumors, and benign and malignant tumor tissues. By integrating the miRNA and mRNA levels, biological functions that were difficult to grasp at each level were identified, and an organic relationship between miRNA and mRNA was confirmed. There is a lack of information regarding the role of miRNA and mRNA levels in canines. In particular, the amount of information is much smaller, as few studies have been conducted on species other than humans or mice. To confirm the role of the miRNAs and mRNA derived from our study, we conducted related studies based on other species.
We investigated various studies to verify the DE-miRNAs and genes identified in each combination. Beuvink et al. (2) conducted a study and defined a novel expression profile of miRNAs in various mouse organs. The authors confirmed that miR-150 showed higher enrichment than the other miRNAs in the spleen. We confirmed that the expression of cfa-miR-150 was up-regulated in Case_21/Control_21 and downregulated in Case_22/Control_22.This indicates that the expression profile of miR-150 appears specifically in benign tumors. Yan et al. (23) reported that high expression of miR-21 in the mouse spleen promotes T cell homeostasis and inhibits the NF-κB signaling pathway to promote the homeostatic proliferation of lymphocytes. We found that the expression of cfa-miR-21 increased in Case_21/Control_21, suggesting that this immune function was activated in the tumor.
Molnar et al. (12) conducted a study on the chemo-preventive effects of various extracts against carcinogens. The authors confirmed that miR-134 expression was significantly reduced when the extract with a preventive effect was administered to the mouse spleen. In Case_22/Control_22, cfa-miR-134 was identified as an upregulated miRNA, and it was inferred that it could be used as an important marker for identifying malignant tumors. miR-197 plays an important role in various types of cancer and is known to inhibit apoptosis. Andrade et al. (1) reported that a decrease in miR-197 expression was associated with the onset and exacerbation of chronic lymphocytic leukemia. Case_22/Control_22 also showed a pattern of decreased cfa-miR-197 expression, which was confirmed to be a major factor related to immunity.
Grimes et al. (8) conducted RNA-seq analysis to identify miRNAs showing significant differences in hemangiosarcoma and nodule proliferation compared to normal tissues. The authors found 22 miRNAs that were differentially expressed in hemangiosarcoma tissues and derived that these candidate miRNAs have potential effects on the cause of hemangiosarcoma. Among these significant miRNAs, we confirmed that miR-139 appeared in the candidate DE-miRNAs of Case_21 vs Control_21 and Case_22 vs Control_22.
The DE-miRNAs identified in Case_22/Case_21 can be used as biomarkers to distinguish benign tumors from malignant tumors or to treat various diseases. Among the proposed studies to date, few approaches have been proposed to detect changes in tumor stage. miRNAs are important substances that can be used to identify various cellular activities from a biological point of view. If an in-depth analysis is performed on the difference between benign and malignant tumors identified in this study, it can be used as a biomarker that can hinder the metastasis or development of tumors in the canine spleen.
In a comparison between target genes affected by miRNAs and differentially expressed genes at the mRNA level, many genes appeared at both levels. The amount of information currently available about the interactions between miRNAs and genes is very limited. Therefore, additional studies on the relationship between novel miRNAs and genes, especially on the individual analysis and interaction at the miRNA and mRNA levels in canines, are necessary. Based on this information, we inferred that there was an organic interaction between miRNAs and mRNA. In addition, because one miRNA is related to several genes, it is very important to identify miRNAs that can affect a gene group of interest. Further studies are required to determine the relationship between miRNAs and genes that have not yet been identified.
This study was supported by the Cooperative Research Program for Agriculture Science and Technology Development (#PJ014990022021), Research Institute for Veterinary Science, and BK Plus 21 Program.
Canine participants enrolled in this study were not purchased for experimental use. This study was approved by the Institutional Animal Care and Use Committee of Seoul National University (approval number: SNU-200217-3-2).
We specially thank the participated companion dogs and their owners for providing samples.
The authors have no conflicting interests.
J Vet Clin 2023; 40(2): 164-174
Published online April 30, 2023 https://doi.org/10.17555/jvc.2023.40.2.164
Copyright © The Korean Society of Veterinary Clinics.
Eunpyo Kim1 , Giup Jang2 , Jin-Wook Kim3 , Wan-Hee Kim4 , Geon-A Kim5,*
1Department of Theriogenology and Biotechnology, Research Institute for Veterinary Science, College of Veterinary Medicine, Seoul National University, Seoul 08826, Korea
2Research and Development Center, Einocle Inc., Seoul 02841, Korea
3Research and Development Center, ROKIT Healthcare Inc., Seoul 08512, Korea
4Department of General Surgery, College of Veterinary Medicine, Seoul National University, Seoul 08826, Korea
5Department of Biomedical Laboratory Science, School of Healthcare Science, Eulji University, Uijeongbu 11759, Korea
Correspondence to:*kimgeona2020@gmail.com
This is an open access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.
Canine splenic tumors (STs) are commonly diagnosed during imaging examinations, such as in X-ray and ultrasonography examinations, suggesting their higher prevalence, especially in older dogs. Despite this high prevalence, there are no effective treatment options for STs because of the difficulties in determining therapeutic targets. However, recently, the importance of microRNAs (miRNAs) has evolved owing to their ambivalent characteristics. Biomarkers and novel therapies using miRNAs have been well-studied in human cancer research compared to canine research, except for mammary gland tumors. Therefore, this study aimed to comparatively analyze miRNA expression profiles according to malignancy and biopsy sites to identify novel therapeutic and diagnostic targets. Tissue samples were collected directly from splenic tumor masses and immersed in RNAlater solution for further analysis. To investigate differentially expressed genes (DEGs) between tumor and normal tissues, we used RNA-seq and miRNA microarray analysis. Then, functional analysis based on DEGs was conducted to sort tumor-related DEGs. We found that cfa-miR-150 was upregulated in benign tumors, whereas cfa-miR-134 was upregulated in malignant tumors. Despite limited information on canine miRNAs, we identified two potential biomarkers for the differential diagnosis of STs.
Keywords: spleen tumor, microRNA, dog, benign, malignancy, messenger RNA
Splenic tumors (STs) in dogs are commonly diagnosed during an ultrasound examination. Half of them are found to be malignant (MS), whereas the others are diagnosed as benign (BS) (13). Also, small-breed dogs were more likely to be diagnosed with benign STs than large-breed dogs (13). Hemangiosarcoma is a common malignant tumor of the spleen in older dogs, with a prevalence of 12-18% in all mesenchymal neoplasms (3). Since hemangiosarcoma has a poor prognosis even with surgery and median survival time is known to be about 3 months, early diagnosis or development of biomarker is essential (20). Because the spleen functions in blood storage, total splenectomy is often required to prevent splenic rupture in dogs with STs. The exact causes of STs are unknown, although it has been described that large breeds (e.g. Golden Retrievers, German Shepherds, Labradors) are predisposed (13). Although total or partial splenectomy is suggested to prevent further intra-abdominal hemorrhage and coagulopathy syndrome, perioperative death is common in dogs because of the presence of higher risks, such as thrombocytopenia, anemia, and intra-operative arrhythmia (22). In addition, even if STs are clearly diagnosed by radiological examinations, the definitive diagnosis is always based on histopathologic analysis (21). Thus, a general analysis is needed for more precise prediction of patients’ prognosis and determining therapeutic strategies.
miRNAs, endogenous small non-coding RNAs that function as regulators of gene expression, have been suggested to be dysregulated in cancer (15). Since miRNAs are expressed differently according to the tissue type and malignancy, they have been studied as biomarkers for cancer diagnosis owing to their non-invasive properties (6,11). Furthermore, as altered miRNAs can result in the development of drug resistance, they are being studied as therapeutic targets and newer drug candidates (19). miRNAs can act as both tumor suppressors and oncogenes based on their mRNA targets (25). For example, miR-12/13 targets bcl-2 and oncogenes (e.g., Cyclin D1, MCL1, CDC2, ETS1, and JUN), showing tumor suppressor effects (14). However, one study reviewed the oncogenic properties of miRNAs in malignant transformation, suggesting that they are potential biomarkers and therapeutic targets in cancers (7). Thus, to develop novel biomarkers and drugs, it is important to profile miRNA expression profiles in veterinary oncology.
Compared to the numerous studies on canine mammary gland tumors, miRNA profiling of canine STs is an area that requires further studies. Therefore, this study aimed to comparatively analyze miRNA profiles in canine splenic tumors based on their malignancy and biopsy sites using histological examination to understand the interconnection between altered miRNAs and their mechanism of action and identify potent biomarkers or novel therapeutic target candidates.
The canine participants enrolled in this study was not purchased for experimental usage. This study was approved by the Seoul National University Institutional Animal Care and Use Committee (approval number: SNU-200217-3-2). All samples are acquired under the consent from every owner.
Tissue samples were collected during surgical operation at veterinary clinics in Republic of Korea. The biopsy samples were collected at the center of splenic cancer masses. Right after the resection, the samples were directly immersed in RNAlater solution purchased from Invitrogen (Invitrogen, MA, USA) and formalin solution for histopathological examination to evaluate malignancy. The stored samples are then preserved in 4°C refrigerator until use. All patients are diagnosed without other comorbidities at the time of physical examination. The patients’ information is summarized in Supplementary Table 1.
Total RNA was extracted from each tissue samples using RNeasy mini kit (Qiagen, Doncaster, VIC, Australia) according to the manufacturer’s protocol. RNA purity and integrity were evaluated using ND-2000 Spectrophotometer (NanoDrop, Wilmington, USA) and Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, USA), respectively.
Microarray hybridization was performed using the Affymetrix GeneChip miRNA 4.0 Array (ThermoFisher,MA, USA) as per manufacturer’s protocol. Total RNA samples were labelled using FlashTagTM Biotin RNA Labeling Kit (Genisphere, PA, USA). The labelled RNA samples were then quantified, fractionated, and hybridized to the microarray according to the manufacturer’s protocol. RNA-array hybridization was performed on an Affymetrix® 450 Fluidics Station (ThermoFisher, MA, USA). The arrays were stained using a GeneChip Fluidics Station 450 (Affymetrix, Santa Clara, California, United States), and scanned using an Affymetrix GCS 3000 scanner (Affymetrix, Santa Clara, California, USA). miRNA-mRNA hybridization signals were analyzed using the Affymetrix® GeneChipTM Command Console.
Total RNA concentration was calculated using Quant-IT RiboGreen (Invitrogen, #R11490). To assess the integrity of the total RNA, samples were run on the TapeStation RNA ScreenTape (Agilent, #5067-5576). Only high-quality RNA preparations with RNA integrity number (RIN) >7.0 were used for RNA library construction. A library was independently prepared with 1 µg total RNA of each sample using the Illumina TruSeq Stranded mRNA Sample Prep Kit (RS-122-2101, Illumina Inc., San Diego, CA, USA). The first step in the workflow involved purifying the poly-A containing mRNA molecules using poly-T-attached magnetic beads. Following purification, the mRNA samples were fragmented into small pieces using divalent cations under elevated temperature. The cleaved RNA fragments were then used as templates to generate first strand cDNA using SuperScript II reverse transcriptase (Invitrogen, MA, USA) with random primers. This was followed by second strand cDNA synthesis using DNA Polymerase I, RNase H and dUTP. The generated double stranded cDNA fragments were then subjected to an end repair process involving adenylation followed by adapter ligation. The products were then purified and enriched with PCR to create the final cDNA library. The libraries were quantified using KAPA Library Quantification kits for Illumina Sequencing platforms according to the qPCR Quantification Protocol Guide (KAPA BIOSYSTEMS, Wilmington, USA) and qualified using the TapeStation D1000 ScreenTape (Agilent Technologies, CA, USA). Indexed libraries were then submitted to Illumina NovaSeq (Illumina Inc., CA, USA) to perform paired-end (2 × 100 bp) sequencing.
At the miRNA level, raw data were extracted automatically through the Affymetrix data extraction protocol using the Affymetrix GeneChip® Command Console® Software (AGCC). The CEL files import, miRNA level RMA+DABG-All analyses and results were exported using Affymetrix® Power Tools (APT) Software. Array data were filtered by probes annotated species.
At the mRNA level, we preprocessed the raw reads from the sequencer to remove low quality and adapter sequences before analysis and aligned the processed reads to the Canis lupus familiaris (CanFam3.1) database using HISAT v2.1.0 (10). HISAT utilizes two types of indexes for alignment (a global, whole-genome index and tens of thousands of small local indexes). These two types’ indexes are constructed using the same BWT (Burrows–Wheeler transform)/ a graph FM index (GFM) as Bowtie2. Because of its use in efficient data structures and algorithms, HISAT generates spliced alignments several times faster than Bowtie and BWA. The reference genome sequence of Canis lupus fa-miliaris (CanFam3.1) and annotation data were downloaded from the NCBI database. Transcript assembly of known transcripts was then processed by StringTie v2.1.3b (16). Based on the results, expression abundance of transcripts was calculated as read count or FPKM value (Fragments Per Kilobase of exon per Million fragments mapped) per sample. The expression profiles were used for additional analysis such as that of differentially expressed genes (DEGs). In groups with different conditions, DEGs or transcripts can be filtered through statistical hypothesis testing.
Comparative analysis between test and control samples was carried out using in-dependent t-test and fold change with the null hypothesis indicating no difference exists among groups. False discovery rate (FDR) was controlled by adjusting p-value using Benjamini-Hochberg algorithm. All Statistical tests and visualization of DEGs was conducted using R statistical language (v3.3.2).
We collected miRNA-disease associations in the Human microRNA Disease Database (HMDD) to confirm that the DE-miRNAs are associated with disease (9). HMDD is a database for human-related miRNA-disease associations. However, the analysis was performed considering heterologous miRNAs for interspecies. Because information on spleen tumor does not exist in HMDD, miRNAs related to lymphatic disease, which is a higher concept than spleen tumor, were collected.
We collected miRNA-gene associations from miRDB and performed gene level analysis based on miRNA targets (4). It contains information about predicted target genes and predicted score for miRNAs. We used associations with a predictive score of ≥50 for miRNA-gene relationships, because the amount of known information about spleen tumor is less.
clusterProfiler v3.18.1 in R was used to perform functional analysis of gene set (24). Gene Ontology includes biological process, cellular component, and molecular function, and Kyoto Encyclopedia of Genes and Genomes (KEGG) includes biological pathway. To verify the disease-level for the identified genes, we collected disease-related genes from DISEASES and DisGeNet (17,18). To validate the results of functional analysis, we collected disease-related Gene Ontology and biological pathway from the Comparative Toxicogenomics Database (CTD) (5).
Differentially expressed miRNAs (DE-miRNAs) were identified for three comparisons:1) Case_21 vs. Control_21, 2) Case_22 vs. Control_22, and 3) Case_22 vs. Case_21. Meaning of Case_21 and Control_21 are tumor tissue of benign ST and normal tissue of benign ST, respectively. Meaning of Case_22 and Control_22 are tumor tissue of malignant ST and normal tissue of malignant ST. |Fold change| ≥ 2 and p-value < 0.05 were used to identify DE-miRNAs for each comparison. Fig. 1 shows a volcano plot for identifying DE-miRNAs based on these cut-off values.
Twenty-nine miRNAs (14 upregulated and 15 downregulated) were identified in Case_21/Control_21, 13 miRNAs (5 upregulated and 8 downregulated) in Case_22/Control_22, and 8 miRNAs (4 upregulated and 4 downregulated) in Case_22/Case_21 (Supplementary File 1).We validated the relationship between DE-miRNAs and diseases using the HMDD. We confirmed that 14 of the 29 miRNAs in Case_21/Control_21 and 7 of the 13 miRNAs were associated with lymphatic disease (Table 1).
Table 1 . Lymphatic disease-associated miRNAs among DE-miRNAs.
Case_21/Control_21 | |||
---|---|---|---|
miRNA | Regulation | Fold change | p-value |
cfa-miR-150 | UP | 2.407 | 0.010 |
cfa-miR-21 | UP | 2.221 | 0.018 |
cfa-miR-221 | UP | 2.022 | 0.039 |
cfa-miR-204 | DOWN | –5.731 | 0.035 |
cfa-miR-222 | UP | 2.536 | 0.002 |
cfa-miR-196b | DOWN | –3.235 | 0.007 |
cfa-miR-224 | DOWN | –14.507 | 0.001 |
cfa-miR-206 | DOWN | –4.508 | 0.003 |
cfa-miR-27b | DOWN | –2.329 | 0.002 |
cfa-miR-194 | DOWN | –2.119 | 0.007 |
cfa-miR-504 | DOWN | –3.728 | 0.035 |
cfa-miR-383* | DOWN | –18.774 | 0.006 |
cfa-miR-133b | DOWN | –14.169 | 0.001 |
cfa-miR-139* | DOWN | –4.927 | 0.023 |
Case_22/Control_22 | |||
cfa-miR-181a | DOWN | –2.568 | 0.011 |
cfa-miR-31 | DOWN | –4.324 | 0.036 |
cfa-miR-326 | DOWN | –4.658 | 0.039 |
cfa-miR-212 | UP | 6.353 | 0.029 |
cfa-miR-197 | DOWN | –2.127 | 0.012 |
cfa-miR-383* | DOWN | –35.603 | 0.005 |
cfa-miR-139* | DOWN | –6.411 | 0.002 |
*Indicates miRNAs that appear in both comparisons..
We confirmed that half of the DE-miRNAs detected in the comparisons were associated with lymphatic disease, and cfa-miR-139 and cfa-miR-383 were common miRNAs. In Case_22/Case_21, miRNAs related to lymphatic disease were not detected; however, three downregulated miRNAs were identified in hematologic neoplasms, a disease at the same level as lymphatic disease (cfa-miR-31, cfa-miR-150, and cfa-miR-708). Thus, three miRNAs were highly expressed in Case_21, and it was confirmed that benign and malignant tumors could be distinguished at the miRNA level. A related study confirmed that hematogenous transmission is the most likely route of transmission, causing splenic metastasis [primary and secondary neoplasms of the spleen]. Differences in miRNA levels between benign and malignant tumors are presumed to be closely related to metastasis and cancer development.
We performed gene-level analysis of DE-miRNAs based on miRDB and miRNA-gene associations collected from the miRNA target prediction database. We collected 5,549 target genes in Case_21/Control_21, 3,038 target genes in Case_22/Control_22, and 2,733 target genes in Case_22/Case_21from their DE-miRNAs. Functional analysis of the gene set was performed using clusterProfiler in the R software. The analysis categories used were GO_BP, GO_CC, GO_MF, and KEGG. GO_BP was a biological process, GO_CC was a cellular component, and GO_MF was a molecular function. To validate whether the results analyzed in each category were related to spleen tumors, disease-related Gene Ontology and KEGG were collected from CTD using the MeSH ID of spleen tumors (D013160). Fig. 2 shows an image of a list associated with spleen tumors among significant Gene Ontologies and KEGG derived from clusterProfiler (p-value < 0.05, Supplementary File 2).
We confirmed that there was no information on spleen tumors in GO and KEGG, but many biological functions related to cancer were identified among the significant categories (Fig. 2). Moreover, we confirmed that the biological functions that prevent the accumulation of damaged biomolecules, such as autophagy and mitophagy, were identified among the significant categories of each comparison.
We collected a list of immune-related pathways from related research to confirm categories related to immunity among the spleen-related function analysis results identified in Fig. 2 [Maternal nematode infection upregulates the expression of Th2 Treg and diapedesis-related genes in the neonatal brain]. Twenty-six KEGG pathway related to immunity are presented in this paper. Table 2 shows a list of genes related to immunity among the spleen-related KEGG genes detected in each comparison.
Table 2 . Immune-related KEGG pathway list based on DE-miRNA target genes.
Case_21/Control_21 | |
---|---|
KEGG pathway | p-value |
T cell receptor signaling pathway | 0.001219192 |
Endocytosis | 2.68E-08 |
Regulation of actin cytoskeleton | 0.000131578 |
Case_22/Control_22 | |
KEGG pathway | p-value |
Fc epsilon RI signaling pathway | 0.00899355 |
T cell receptor signaling pathway | 0.026774021 |
Endocytosis | 3.26E-10 |
Regulation of actin cytoskeleton | 0.003272488 |
Case_22/Case_21 | |
KEGG pathway | p-value |
Chemokine signaling pathway | 0.049509358 |
Fc epsilon RI signaling pathway | 0.04652657 |
T cell receptor signaling pathway | 0.012534506 |
Endocytosis | 4.23E-06 |
Regulation of actin cytoskeleton | 0.001236387 |
We confirmed that the T cell receptor signaling pathway, endocytosis, and regulation of the actin cytoskeleton appeared in common in the three comparisons. The Fc epsilon RI signaling pathway was confirmed in Case_22/Control_22, and fc epsilon RI and chemokine signaling pathways were detected in Case_22/Case_21.In immune-related research on the bone marrow and spleen, it was confirmed that the T cell receptor signaling pathway, chemokine signaling pathway, and fc epsilon RI signaling pathway were the most enriched pathways [extracellular vesicles mediate radiation-induced systemic bystander signals in the bone marrow and spleen].
We collected disease-related genes from DISEASES, DisGeNet, and identified disease-related genes among the target genes. Ten spleen cancer-related genes and 10 spleen disease-related genes were collected from DISEASES, and 6 splenic neoplasm-related genes were collected from DisGeNet. Because known information about spleen tumors is scarce, several databases and diseases have been used. Table 3 shows a list of disease-related genes among miRNA target genes.
Table 3 . Disease-related genes for miRNA target genes.
Case_21/Control_21 | ||
---|---|---|
miRNA | Target gene | Database |
cfa-miR-379 | NOTCH2 | DISEASES/Spleen cancer |
cfa-miR-27b | ||
cfa-miR-206 | ||
cfa-miR-299 | CD19 | DISEASES/Spleen cancer |
cfa-miR-150 | ITGAE | DISEASES/Spleen cancer |
cfa-miR-1 | CD5 | DISEASES/Spleen cancer DISEASES/Spleen disease |
cfa-miR-206 | ||
cfa-miR-204 | ||
cfa-miR-543 | LMLN | DISEASES/Spleen disease |
cfa-miR-376a | ||
cfa-miR-299 | THPO | DISEASES/Spleen disease |
cfa-miR-133b | CD8A | DISEASES/Spleen disease |
cfa-miR-1841 | CAST | DisGeNet/Splenic neoplasm |
Case_22/Control_22 | ||
miRNA | Target gene | Database |
cfa-miR-330 | CD5 | DISEASES/Spleen disease |
cfa-miR-376a | LMLN | DISEASES/Spleen disease |
cfa-miR-181a | CD4 | DISEASES/Spleen disease |
cfa-miR-31 | TNF | DISEASES/Spleen disease |
Case_22/Case_21 | ||
miRNA | Target gene | Database |
cfa-miR-342 | CD79B | DISEASES/Spleen cancer |
cfa-miR-665 | CD79A | DISEASES/Spleen cancer |
cfa-miR-150 | ITGAE | DISEASES/Spleen cancer |
cfa-miR-379 | NOTCH2 | DISEASES/Spleen cancer |
cfa-miR-27b | ||
cfa-miR-206 | ||
cfa-miR-31 | TNF | DISEASES/Spleen disease |
CD5 and LMLN were identified as genes commonly identified in benign and malignant tumors, and the remaining genes were identified as genes uniquely expressed in benign or malignant tumors (Fig. 3). In addition, CD79A and CD79B were uniquely derived in benign and malignant tumors compared with the differences between normal and tumor areas. Although the number of genes known to be related to the spleen was small, we confirmed that a large number of disease-related genes were detected among genes affected by differentially expressed miRNAs. Differentially expressed genes (DEGs) at the mRNA level were identified using the same comparisons and conditions (|fold change| ≥ 2 and p-value < 0.05). Fig. 3 shows the volcano plot for identifying DEGs based on these cutoff values.
A total of 1,238 genes (307 upregulated and 931 downregulated) were identified in Case_21/Control_21, 904 genes (369 upregulated and 535 downregulated) in Case_22/Control_22, and 1,156 genes (582 upregulated and 574 downregulated) in Case_22/Case_21 (Supplementary File 3).
Functional analysis of the DEG list was performed using clusterProfiler in R. Similar to the miRNA level, disease-related Gene Ontology and KEGG were collected from CTD using MeSH ID of spleen tumor. Fig. 4 shows an image of a list associated with spleen tumors among significant Gene Ontologies and KEGG derived from clusterProfiler (p-value < 0.05, Supplementary File 4).
Similar to the analysis results at the miRNA level, we confirmed the appearance of multiple cancer-related biological functions in the three comparisons shown in Fig. 4. In particular, the cell cycle, T cell receptor signaling pathway, B cell receptor signaling pathway, and apoptosis were also derived from the difference between benign and malignant tumors as significant functions, along with the previous results.
We identified categories related to immunity among the significant categories related to the spleen (Fig. 4). Table 4 shows a list related to immunity among the spleen-related KEGG genes detected in each comparison at the mRNA level.
Table 4 . Immune-related KEGG pathway list based on DEGs.
Case_21/Control_21 | |
---|---|
KEGG pathway | p-value |
Chemokine signaling pathway | 0.000612405 |
Regulation of actin cytoskeleton | 1.41E-05 |
Case_22/Control_22 | |
KEGG pathway | p-value |
Fc epsilon RI signaling pathway | 0.003332669 |
Regulation of actin cytoskeleton | 0.000374166 |
Case_22/Case_21 | |
KEGG pathway | p-value |
T cell receptor signaling pathway | 2.87E-09 |
B cell receptor signaling pathway | 0.009193225 |
Regulation of actin cytoskeleton | 0.02036557 |
The pathways identified at the miRNA level also appeared identical at the mRNA level (Table 3). At the mRNA level, the B cell receptor signaling pathway was also derived from the comparison between benign and malignant tumors, and it was confirmed that mature B lymphocytes are closely related to the immune response in the spleen (B cell development in the spleen takes place in discrete steps and is determined by the quality of B cell receptor-derived signals).
Based on the analysis of miRNA and mRNA levels, we examined the correlation between them. First, genes that appeared in common among the target genes of DE-miRNAs and the DEG of mRNA were identified. Among the DEGs identified in the mRNA analysis, only genes corresponding to protein coding were used. Fig. 5 presents a Venn diagram of the genes identified at both levels.
We commonly identified 423 genes in Case_21/Control_21, 173 genes in Case_22/Control_22 and 194 genes in Case_22/Case_21 at the miRNA and mRNA levels (Supplementary File 5). ClusterProfiler was used to confirm the functions of the common genes. Unlike the functional analysis conducted previously at the miRNA and mRNA levels, the analysis results for the top five categories based on p-values in each category were included, not limited to the spleen. Fig. 6 shows the results of functional analysis of the common genes of the three comparisons (Supplementary File 6).
The analysis results based on common genes showed different aspects from the previous analysis results at the miRNA and mRNA levels. Among the significant categories, biological functions related to blood vessels and muscles were identified, and in the comparison between benign and malignant tumors, it was confirmed that significant results were also derived in cellular senescence and dendritic-related functions. Based on the results of the previous DE-miRNA target gene-based functional analysis and mRNA DEG-based functional analysis, a Venn diagram was constructed for categories related to the spleen, as shown in Fig. 7.
Fourteen biological functions in Case_21/Control_21, 26 functions in Case_22/Control_22, and 25 functions in Case_22/Case_21 were identified at the miRNA and mRNA levels (Supplementary File 7). We recognized that more diverse biological functions might appear because the number of miRNA target genes is greater than the number of differentially expressed genes at the mRNA level. The number of functions of the gene clusters identified at the miRNA level was high. Nevertheless, we confirmed that most functions of the genes analyzed at the miRNA and mRNA levels were shared.
We performed various comparisons of canine spleen tumors. We derived various insights into normal and tumor tissues for benign tumors, normal and tumor tissues for malignant tumors, and benign and malignant tumor tissues. By integrating the miRNA and mRNA levels, biological functions that were difficult to grasp at each level were identified, and an organic relationship between miRNA and mRNA was confirmed. There is a lack of information regarding the role of miRNA and mRNA levels in canines. In particular, the amount of information is much smaller, as few studies have been conducted on species other than humans or mice. To confirm the role of the miRNAs and mRNA derived from our study, we conducted related studies based on other species.
We investigated various studies to verify the DE-miRNAs and genes identified in each combination. Beuvink et al. (2) conducted a study and defined a novel expression profile of miRNAs in various mouse organs. The authors confirmed that miR-150 showed higher enrichment than the other miRNAs in the spleen. We confirmed that the expression of cfa-miR-150 was up-regulated in Case_21/Control_21 and downregulated in Case_22/Control_22.This indicates that the expression profile of miR-150 appears specifically in benign tumors. Yan et al. (23) reported that high expression of miR-21 in the mouse spleen promotes T cell homeostasis and inhibits the NF-κB signaling pathway to promote the homeostatic proliferation of lymphocytes. We found that the expression of cfa-miR-21 increased in Case_21/Control_21, suggesting that this immune function was activated in the tumor.
Molnar et al. (12) conducted a study on the chemo-preventive effects of various extracts against carcinogens. The authors confirmed that miR-134 expression was significantly reduced when the extract with a preventive effect was administered to the mouse spleen. In Case_22/Control_22, cfa-miR-134 was identified as an upregulated miRNA, and it was inferred that it could be used as an important marker for identifying malignant tumors. miR-197 plays an important role in various types of cancer and is known to inhibit apoptosis. Andrade et al. (1) reported that a decrease in miR-197 expression was associated with the onset and exacerbation of chronic lymphocytic leukemia. Case_22/Control_22 also showed a pattern of decreased cfa-miR-197 expression, which was confirmed to be a major factor related to immunity.
Grimes et al. (8) conducted RNA-seq analysis to identify miRNAs showing significant differences in hemangiosarcoma and nodule proliferation compared to normal tissues. The authors found 22 miRNAs that were differentially expressed in hemangiosarcoma tissues and derived that these candidate miRNAs have potential effects on the cause of hemangiosarcoma. Among these significant miRNAs, we confirmed that miR-139 appeared in the candidate DE-miRNAs of Case_21 vs Control_21 and Case_22 vs Control_22.
The DE-miRNAs identified in Case_22/Case_21 can be used as biomarkers to distinguish benign tumors from malignant tumors or to treat various diseases. Among the proposed studies to date, few approaches have been proposed to detect changes in tumor stage. miRNAs are important substances that can be used to identify various cellular activities from a biological point of view. If an in-depth analysis is performed on the difference between benign and malignant tumors identified in this study, it can be used as a biomarker that can hinder the metastasis or development of tumors in the canine spleen.
In a comparison between target genes affected by miRNAs and differentially expressed genes at the mRNA level, many genes appeared at both levels. The amount of information currently available about the interactions between miRNAs and genes is very limited. Therefore, additional studies on the relationship between novel miRNAs and genes, especially on the individual analysis and interaction at the miRNA and mRNA levels in canines, are necessary. Based on this information, we inferred that there was an organic interaction between miRNAs and mRNA. In addition, because one miRNA is related to several genes, it is very important to identify miRNAs that can affect a gene group of interest. Further studies are required to determine the relationship between miRNAs and genes that have not yet been identified.
This study was supported by the Cooperative Research Program for Agriculture Science and Technology Development (#PJ014990022021), Research Institute for Veterinary Science, and BK Plus 21 Program.
Canine participants enrolled in this study were not purchased for experimental use. This study was approved by the Institutional Animal Care and Use Committee of Seoul National University (approval number: SNU-200217-3-2).
We specially thank the participated companion dogs and their owners for providing samples.
The authors have no conflicting interests.
Table 1 Lymphatic disease-associated miRNAs among DE-miRNAs
Case_21/Control_21 | |||
---|---|---|---|
miRNA | Regulation | Fold change | p-value |
cfa-miR-150 | UP | 2.407 | 0.010 |
cfa-miR-21 | UP | 2.221 | 0.018 |
cfa-miR-221 | UP | 2.022 | 0.039 |
cfa-miR-204 | DOWN | –5.731 | 0.035 |
cfa-miR-222 | UP | 2.536 | 0.002 |
cfa-miR-196b | DOWN | –3.235 | 0.007 |
cfa-miR-224 | DOWN | –14.507 | 0.001 |
cfa-miR-206 | DOWN | –4.508 | 0.003 |
cfa-miR-27b | DOWN | –2.329 | 0.002 |
cfa-miR-194 | DOWN | –2.119 | 0.007 |
cfa-miR-504 | DOWN | –3.728 | 0.035 |
cfa-miR-383* | DOWN | –18.774 | 0.006 |
cfa-miR-133b | DOWN | –14.169 | 0.001 |
cfa-miR-139* | DOWN | –4.927 | 0.023 |
Case_22/Control_22 | |||
cfa-miR-181a | DOWN | –2.568 | 0.011 |
cfa-miR-31 | DOWN | –4.324 | 0.036 |
cfa-miR-326 | DOWN | –4.658 | 0.039 |
cfa-miR-212 | UP | 6.353 | 0.029 |
cfa-miR-197 | DOWN | –2.127 | 0.012 |
cfa-miR-383* | DOWN | –35.603 | 0.005 |
cfa-miR-139* | DOWN | –6.411 | 0.002 |
*Indicates miRNAs that appear in both comparisons.
Table 2 Immune-related KEGG pathway list based on DE-miRNA target genes
Case_21/Control_21 | |
---|---|
KEGG pathway | p-value |
T cell receptor signaling pathway | 0.001219192 |
Endocytosis | 2.68E-08 |
Regulation of actin cytoskeleton | 0.000131578 |
Case_22/Control_22 | |
KEGG pathway | p-value |
Fc epsilon RI signaling pathway | 0.00899355 |
T cell receptor signaling pathway | 0.026774021 |
Endocytosis | 3.26E-10 |
Regulation of actin cytoskeleton | 0.003272488 |
Case_22/Case_21 | |
KEGG pathway | p-value |
Chemokine signaling pathway | 0.049509358 |
Fc epsilon RI signaling pathway | 0.04652657 |
T cell receptor signaling pathway | 0.012534506 |
Endocytosis | 4.23E-06 |
Regulation of actin cytoskeleton | 0.001236387 |
Table 3 Disease-related genes for miRNA target genes
Case_21/Control_21 | ||
---|---|---|
miRNA | Target gene | Database |
cfa-miR-379 | NOTCH2 | DISEASES/Spleen cancer |
cfa-miR-27b | ||
cfa-miR-206 | ||
cfa-miR-299 | CD19 | DISEASES/Spleen cancer |
cfa-miR-150 | ITGAE | DISEASES/Spleen cancer |
cfa-miR-1 | CD5 | DISEASES/Spleen cancer DISEASES/Spleen disease |
cfa-miR-206 | ||
cfa-miR-204 | ||
cfa-miR-543 | LMLN | DISEASES/Spleen disease |
cfa-miR-376a | ||
cfa-miR-299 | THPO | DISEASES/Spleen disease |
cfa-miR-133b | CD8A | DISEASES/Spleen disease |
cfa-miR-1841 | CAST | DisGeNet/Splenic neoplasm |
Case_22/Control_22 | ||
miRNA | Target gene | Database |
cfa-miR-330 | CD5 | DISEASES/Spleen disease |
cfa-miR-376a | LMLN | DISEASES/Spleen disease |
cfa-miR-181a | CD4 | DISEASES/Spleen disease |
cfa-miR-31 | TNF | DISEASES/Spleen disease |
Case_22/Case_21 | ||
miRNA | Target gene | Database |
cfa-miR-342 | CD79B | DISEASES/Spleen cancer |
cfa-miR-665 | CD79A | DISEASES/Spleen cancer |
cfa-miR-150 | ITGAE | DISEASES/Spleen cancer |
cfa-miR-379 | NOTCH2 | DISEASES/Spleen cancer |
cfa-miR-27b | ||
cfa-miR-206 | ||
cfa-miR-31 | TNF | DISEASES/Spleen disease |
Table 4 Immune-related KEGG pathway list based on DEGs
Case_21/Control_21 | |
---|---|
KEGG pathway | p-value |
Chemokine signaling pathway | 0.000612405 |
Regulation of actin cytoskeleton | 1.41E-05 |
Case_22/Control_22 | |
KEGG pathway | p-value |
Fc epsilon RI signaling pathway | 0.003332669 |
Regulation of actin cytoskeleton | 0.000374166 |
Case_22/Case_21 | |
KEGG pathway | p-value |
T cell receptor signaling pathway | 2.87E-09 |
B cell receptor signaling pathway | 0.009193225 |
Regulation of actin cytoskeleton | 0.02036557 |