Differences in miRNA differential expression in whole blood between horses with sarcoid regression and progression.

BACKGROUND: Currently no methods are available to predict the clinical outcome of individual horses with equine sarcoid (ES) disease.
OBJECTIVE: To investigate if whole blood microRNA (miRNA) profiles can predict the long-term development of ES tumors. ANIMALS: Five horses with regression and 5 with progression of ES lesions monitored over 5-7 years and 5 control horses free of ES for at least 5 years.
METHODS: For this cohort study, RNA extracted from whole blood samples from the regression, progression, and control groups was used for high throughput sequencing. Known and novel miRNAs were identified using miRDeep2 and differential expression analysis was carried out by the DESeq2 algorithm. Target gene and pathway prediction as well as enrichment and network analyses were conducted using TarBase, mirPath, and metaCore from GeneGo.
RESULTS: Fourteen miRNAs were differentially expressed between regression and progression groups after accounting for the control condition: 4 miRNAs (28.6%) were upregulated and 10 miRNAs (71.4%) were downregulated with >2-fold change. Seven of the 10 downregulated miRNAs are encoded in an miRNA cluster on equine chromosome 24, homologous to the well-known 14q32 cluster in humans. Their target genes show enrichment for pathways involved in viral carcinogenesis. CONCLUSIONS AND CLINICAL IMPORTANCE: Whole blood miRNA expression profiles are associated with long-term ES growth in horses and warrant further validation as prognostic biomarkers in a larger study cohort. Deregulation of miRNAs on equine chromosome 24 might represent a trigger for ES development.

bovine papilloma virus

equine sarcoid

human papilloma virus

likelihood ratio test

microRNA

principal component analysis

RNA integrity number

INTRODUCTION

Equine sarcoid (ES) disease is the most common neoplastic skin condition of horses.1, 2 Its pathogenesis is complex with host genetic susceptibility and bovine papilloma virus (BPV) 1 and 2 infection representing the most important inciting factors. The risk for being affected by ES may be influenced by specific alleles or haplotypes at the major histocompatibility complex,3, 4, 5 but recent population genetics and genome‐wide association studies suggest that susceptibility to ES has a polygenic basis.6, 7, 8 Once affected by ES, the clinical course of the disease can differ dramatically among individuals. Approximately 48% of single occult and verrucous sarcoids in young horses may show complete spontaneous regression over 5‐7 years.9 Such milder forms of ES, however, can transform into aggressive types such as fibroblastic, mixed fibroblastic, or malevolent sarcoids. They may spread extensively or be located in anatomically difficult areas, and thus compromise the health, welfare, and value of the affected animal.10

The diverse clinical forms of ES are still poorly understood. Although it is well accepted that exogenous factors such as trauma, including diagnostic biopsies or surgical intervention, can trigger or exacerbate the disease,11, 12 far less is known about the endogenous factors that may influence disease outcome. Unlike other neoplastic diseases, histological characterization is not useful to discriminate between milder and more aggressive sarcoid phenotypes; increased density of dermal fibroblasts is common to all forms and other microscopic alterations are too variable to differentiate clinical types.13 The DNA of BPV is consistently present in ES, but BPV viral load and abnormal p53 expression do not appear to be associated with ES progression.14 In contrast, an in vitro study suggested that inhibition of BPV gene expression might reverse the transformed phenotype of sarcoid fibroblasts and therefore limit proliferation and invasiveness.15 Recently, abnormal DNA methylation patterns have been detected in sarcoid samples and these epigenetic changes might help to assess growth behavior of the tumors.16 In summary, it remains unknown why some horses show severe exacerbations and others are clinically stable or even experience spontaneous regression.

MicroRNAs (miRNAs) are small, 20‐23 nucleotide long, noncoding RNA molecules that mainly downregulate the expression of their target genes. MicroRNA expression is dysregulated in various malignancies,17 and tumor‐type specific miRNA fingerprints have been proposed as promising diagnostic and prognostic cancer biomarkers.18 Changes in miRNA expression have been observed in ES tissue and BPV‐transformed equine fibroblast cell lines.19, 20 Furthermore, miRNA fingerprints have prognostic potential in human papilloma virus (HPV)‐associated cancers.21, 22, 23

Thus, our aim was to compare miRNA expression profiles in whole blood of horses with regression versus progression of ES disease and ES‐free control horses, assess their potential as noninvasive predictors of the course of ES disease, and predict their possible biological function.

MATERIAL AND METHODS

Study cohort

This study was fully approved by the Ethical Committee of the Canton of Bern (VD 2227.0, VD 2227.1 und VD 2227.2, BE110/15). Animals were client‐owned (regression and progression groups) or clinic‐owned horses (control group) that were cared for according to local regulations. Whole blood samples collected into EDTA tubes for a previous longitudinal study9 were available for discovery of prognostic miRNA fingerprints. They had been taken from 61 Franches‐Montagnes horses that had been diagnosed with ES at the age of 3 years (1st examination and blood collection). After 5‐7 years, a standardized clinical follow‐up with a focus on ES lesion behavior was performed on all horses. Horses were defined as suffering from ES based on the characteristic clinical appearance of skin lesions and excluded when the diagnosis of ES was unclear because of atypical appearance or location of the lesions. Animals were then grouped based on progression versus regression of lesions. For this study, 5 horses with marked regression (complete disappearance of lesions at 2nd examination) and 5 horses with marked progression (exacerbation at 2nd examination with more aggressive phenotype[s]) were selected. Whole blood samples also were collected from a control group consisting of 5 clinic‐owned horses ≥12 years of age that were free of sarcoid lesions for at least 5 years and free of other signs of systemic disease at the time of examination. In control horses, whole blood collection into EDTA tubes was performed by the end of the 5‐year observation period. Because of the high costs associated with next generation sequencing, a low number of study subjects (5 horses per group) were selected for this exploratory cohort study. Because no general consensus currently exists on what the sample size should be for miRNA biomarker studies, a previously recommended approach was used.24

RNA sample preparation and sequencing of multiplexed small RNA samples

Whole blood samples used in this study were stored for long term (9‐11 years) for horses in the ES regression and progression groups and for short term (6 weeks) for control horses at −80°C before the analysis for this study. Samples were handled according to a previously established protocol for long‐term stored equine whole blood.25 The blood was gently thawed on ice and transferred to a PAXgene blood RNA tube 16 hours before the start of RNA extraction with a PAXgene Blood miRNA Kit (Qiagen, Hombrechtikon, Switzerland). The RNA concentrations were measured with Qubit RNA BR assay kit (ThermoFisher, Basel, Switzerland). The RNA integrity number (RIN) was assessed using Fragment Analyzer with the Standard Sensitivity RNA Analysis Kit (Advanced Analytical, Heidelberg, Germany). A total of 14.6‐54.6 ng of small RNA was used for stranded single‐end library preparation (NEBNext Multiplex Small RNA Library Prep Set for Illumina) and sequenced (HiSeq3000, Illumina).

Raw data preprocessing

An average of approximately 17 million reads was collected for each of the 15 libraries. Raw sequencing data quality was assessed using FastQC software.26 Truncated adapter sequences were reported by FastQC to be over‐represented. Cutadapt (v 1.8)27 was used to trim the full and truncated adapter sequences and to filter out reads <15 bp with the following options: ‐‐trim‐n a FASTA file> ‐m 15 aq 20, in which “‐m” denotes the minimum fragment length in bp and “‐q” denotes the minimum read quality score encoded in asci (phred quality + 33) format. The raw sequence dataset is available at the European Bioinformatics Institute http://www.ebi.ac.uk/ena/data/view/PRJEB27174.

MicroRNA annotation

Sequencing reads that passed quality control were mapped to the horse reference genome (EquCab2) available from NCBI using the miRDeep2 (v0.0.7) package.28 The processed reads from all 15 FASTQ files were mapped by the mapper.pl script of the miRDeep2 tool using a configuration file containing a list of all FASTQ file names and a 3‐letter code for each sample (used to label each read according to the phenotype [regression, progression, and control]). The module mapped the reads to the genome with Bowtie129 keeping only the alignments with 0 mismatches (option ‐n) in the seed region and reads that did not map to >5 different loci in the genome (option ‐m). The module was run with the following parameters: ‐d ‐v ‐u ‐e ‐h ‐o 16 ‐m ‐n ‐q ‐p. This procedure generated nonredundant reads with the counts, source, and genomic location of the mapped reads. The miRDeep2 then was used to annotate putative novel miRNAs as well as validate predicted miRNAs from miRBase30 using (1) mapped reads from the mapper.pl output, (2) equine mirBase predicted miRNAs, and (3) human miRBase predicted miRNAs.28 The putative novel miRNAs predicted by miRbase were assigned a score based on the 5p and 3p support and secondary structure that conformed to biogenesis miRNAs. Putative novel miRNAs predicted by miRDeep2 with a score of 4 were retained for further analysis because the signal over noise ratio for this cutoff was >10 at 14.3. Known equine and human miRNAs were obtained from the miRBase database, release 21.30 Furthermore, overlapping novel miRNAs (by at most 1 nucleotide) were merged by BEDtools.31 The code used to implement this approach is outlined in Supporting Information File S1.

Differential expression analysis of miRNA

The expression level quantification of novel (including merged ones) and known miRNA, predicted by miRDeep2, was performed by the quantifier.pl script applying the parameters ‐k ‐j. This script produced a raw counts file that then was used as input into the DESeq R/Bioconductor package for further statistical analysis.32 Differential expression tests in DESeq were conducted only for miRNA genes with read counts in ≥1 of the samples. For DESeq2, the “DESeq” function was used to estimate size factors, dispersions, and finally fit a model. The DeSeq2 likelihood ratio test (LRT) was used to contrast the conditions (regression ‐ control) ‐ (progression ‐ control) and obtain differential expressed miRNA among the 3 levels. The bioconductor package DegReport was used to visualize the results.33 MiRNAs were considered as differentially expressed if the absolute value of log2‐fold change was >1 and the P value adjusted for multiple testing was <.05.

MicroRNA target genes, enrichment, pathway, and network analysis

Enrichment analysis with GeneCodis34, 35, 36 was performed using target genes of differentially expressed miRNAs determined by TarBase v7.0 (Diana Tools),37, 38 which catalogs only published, experimentally validated miRNA:gene interactions. Pathway prediction of the union of genes regulated by the differentially expressed miRNAs was performed by mirPath v.3 (Diana Tools).38 MetaCore v.6.32 (Thomson Reuters) software was used to build networks of the differentially expressed miRNAs. Expanded networks (auto expand algorithm default: 50 nodes) were built using known and predicted targets of the differentially expressed miRNAs.

RESULTS

The regression group consisted of 3 mares and 2 stallions, and the progression group of 3 mares and 2 geldings. All horses were Franches‐Montagnes horses and 3 years of age at the time of sampling. Sarcoid types consisted only of occult and verrucous lesions in both groups at the 1st presentation and verrucous, nodular, and fibroblastic lesions in the progression group at the 2nd examination. Anatomical sites and types of the lesions as well as treatment attempts are presented in Supporting Information Table S1. The control group consisted of 2 mares, 2 stallions, and 1 gelding aged 12‐23 years (mean ± SD, 16 ± 4.36 years). The control group comprised 3 Franches‐Montagnes and 2 Swiss Warmblood horses.

RNA quality and quantity control

The concentration of RNA used for library preparation ranged from 14.0 to 59.6 ng/μL (mean ± SD, 34.1 ± 14.5 ng/μL). The RNA concentration in whole blood was 33.7 ± 16.5 ng/μL in ES‐affected horses and 35.1 ± 11.2 ng/μL in control horses (Supporting Information Table S1). The RIN ranged from 1.0 to 8.3 (3.2 ± 2.4). Mean RIN was 2.2 ± 2.4 in whole blood of ES‐affected horses and 5.2 ± 0.6 in control horses (Supporting Information Table S1). Because the correlation between the coefficient of variation of expressed miRNAs and RIN quality of the input RNA was 0.2, control for RIN in differential gene expression analysis was not performed.

Read sequence data quality control

Approximately 17 million reads from the 15 libraries were collected and after trimming for adapters and low‐quality bases, approximately 11 million reads remained. The number of unique reads mapped by Bowtie1 to the horse reference genome was 981 012. A summary of miRNA read sequencing is shown in Supporting Information Table S2. Figure 1 displays the read length distribution for adapter and quality‐trimmed reads for the 3 groups of samples.

Figure 1

Read length distribution of all miRNA raw read libraries of control horses (CTL), horses with progressive sarcoids (PRG), and horses with regressive sarcoids (RGR). len, length

Read mapping and putative novel miRNA identification

The miRDeep2 algorithm reported 786 putative novel miRNAs with miRDeep2 score between 0 and 10 (Supporting Information Table S3) and detected 271 novel miRNAs and 313 known miRNAs that passed the relatively stringent score cutoff of 4. These then were used for studying the differential expression with the DeSeq2 tool. The novel miRNAs also were blasted against the Rfam database for potential hits to rRNA or tRNA, which did not result in any significant matches. These resulting high confidence novel miRNAs are listed in Supporting Information Table S4.

Differential expression analysis with DeSeq2

The 522 miRNAs with raw counts values quantified by miRDeep (from 15 samples) were used for differential analysis. After filtering rows with zero expression from the 522 miRNAs, 514 miRNA genes remained. A principal component analysis (PCA) plot was generated to determine the overall sample clustering based on the 200 most variable miRNAs. Three outliers were removed based on PCA clustering (Supporting Information Figures S1 and S2). The 3 outliers were ES_1595_E belonging to the control group, ES_1660_E belonging to the progression group, and ES_1677_E belonging to the regression group. The miRNA read counts from the remaining 12 samples were used for the differential expression.

As recommended by DESeq2 for this type of experimental design, a LRT was performed to test for differences in gene expression among the sample groups (in this case, a condition factor with 3 levels: regression, progression, and control). Fourteen differentially expressed miRNAs were detected between the regression and progression group which were accounted for by control condition (Table 1). The differential expression was identified with an adjusted P value <.05 and a fold change of approximately >2. Four miRNAs (28.6%) were upregulated and 10 miRNAs (71.4%) were downregulated in the regression group (Table 1). The degPlot function from DEGreport was used to view the expression of the top 10 differentially expressed miRNAs (Figure 2). The degPatterns function from the R package DEGreport was used to extract and plot genes that had a similar trend across the 3 conditions (control [CTL], progression [PRG], and regression [RGR] on the x‐axis; Supporting Information Figure S3).

Table 1

Differentially expressed known and putative novel miRNAs in horses with regressive (regression‐control) and progressive sarcoids (progression‐control)

	baseMean	log2‐FoldChange	P value	P adj
Eca‐miR‐381	63.6763	−3.4913	2.23E‐09	9.97E‐07
Eca‐miR‐134	39.5413	−3.5363	3.13E‐07	6.99E‐05
Eca‐miR‐127	138.0183	−3.6306	6.20E‐07	9.24E‐05
Eca‐miR‐382	15.4181	−4.2427	1.56E‐05	.00174842
Eca‐miR‐379	10.6011	−5.0780	3.49E‐05	.003
20_15133	37.1299	1.7902	.0001	.009
Eca‐miR‐107b	218.9891	−1.3709	.0002	.01
Eca‐miR‐146a	1252.0312	−1.7017	.0002	.01
Eca‐miR‐125a‐5p	25.2910	1.1047	.0007	.03
Eca‐miR‐432	2.9624	−5.3878	.0007	.03
Eca‐miR‐24	1371.2551	0.0419	.0009	.03
Eca‐miR‐323‐5p	14.6754	−3.2733	.001	.04
Eca‐miR‐1249	2192.4562	0.7270	.001	.04
Eca‐miR‐378	1234.2905	−1.5619	.001	.04

Abbreviation: P adj = adjusted P value.

Figure 2

The top 10 differentially expressed miRNA patterns in control horses (CTL) and horses with progressive (PRG), and regressive sarcoids (RGR). The y‐axis shows normalized counts of the respective miRNA

MicroRNA target genes, enrichment, and pathway analysis

Using TarBase v7.0 (Diana Tools), 265 target genes with prediction score ≥0.9 (suggested threshold for stringent analysis) were predicted for the 14 differentially expressed miRNAs in the regression versus progression groups (Supporting Information Table S5). These target genes then were used for enrichment and pathway analysis. The results of the singular enrichment analysis of target genes using GeneCodis for gene ontology biological process and molecular function showed enrichment for cancer and cancer‐related pathways (Supporting Information Table S6). MirPath v.3 (Diana Tools) predicted 67 pathways for all differentially expressed miRNAs in the regression versus progression group (P < .05) with the pathway “viral carcinogenesis” being between the top 2 pathways (Table 2).

Table 2

Pathway prediction for the differentially expressed miRNAs in horses with regressive sarcoids versus progressive sarcoids according to mirPath v.3 (Diana Tools) (P ≤ .05). The number of miRNAs associated with each pathway and number of genes involved in each pathway are indicated. The top 10 pathways are listed

KEGG pathway	P value	Number of genes	Number of miRNAs
Proteoglycans in cancer	4.38E‐18	106	10
Viral carcinogenesis	3.44E‐09	103	10
Fatty acid metabolism	3.52E‐09	22	4
Prion diseases	3.52E‐09	15	7
Hippo signaling pathway	3.55E‐07	70	9
Chronic myeloid leukemia	7.46E‐07	45	8
Cell cycle	1.44E‐06	72	9
Fatty acid biosynthesis	2.20E‐06	6	4
Lysine degradation	2.20E‐06	28	6
Hepatitis B	2.20E‐06	70	9

Biomarker‐disease associations and network analysis for downregulated miRNAs in the regression versus progression group

To determine the comprehensive functional character of the differentially expressed miRNAs, ontology enrichment analysis was performed with target genes of up‐ and downregulated miRNAs by the MetaCore software algorithm of GeneGo. The direct target genes of downregulated miRNAs in the regression group were top sarcoma and connective/soft tissue neoplasm biomarkers according to GeneGo disease ontology (Figure 3).

Figure 3

The top GeneGo disease ontology terms associated with direct targets of downregulated miRNAs in horses with regressive versus progressive sarcoids. The y‐axis represents ranks according to the highest log of the P value in the x‐axis

The 10 miRNAs found to be downregulated in whole blood of the regression versus progression group are likely to be associated with a more favorable prognosis. Downregulation of some of their human analogs in plasma also has been associated with a more favorable prognosis in human cancers (Supporting Information Table S7). Literature‐based network analysis using GeneGo for downregulated miRNAS in whole blood of horses with regression of ES disease predicted miR‐146a, miR‐381, and miR‐107 as central regulatory elements (Figure 4). MiR‐146a, miR‐381, and miR‐107 deregulation also has been associated with HPV‐associated cancers (Supporting Information Table S7). Furthermore, the expression of c‐Myc, a major oncogene, is predicted to be regulated by 2 of the downregulated miRNAs (miR‐134 and miR‐382), which in human beings are described to be part of a large miRNA cluster on chromosome 14q32. This miRNA cluster is dysregulated in several types of cancers in humans (Supporting Information Table S8). Further miRNAs downregulated in whole blood of horses with ES regression and belonging to this cluster are miR‐381, miR‐127, miR‐379, miR‐432, and miR‐323‐5p. In the horse, these 7 miRNAs all are encoded by 3 polycistronic miRNA clusters on equine chromosome 24.39

Figure 4

Network for downregulated miRNA genes in horses with regressive versus progressive sarcoids. This figure illustrates the complex interactions of the downregulated miRNA genes in the regression versus progression group. MiR‐146a, miR‐381, and miR‐107 seem to represent central hubs that may influence the expression of various, cancer‐relevant proteins, among them c‐Myc, a major oncogene, and p53, one of the most important tumor suppressor genes

Enrichment and network analysis for upregulated miRNAs in the regression versus progression group

The direct targets of the upregulated genes show enrichment in network processes such as cell cycle, DNA damage, and DNA mismatch repair (Figure 5).

Figure 5

Top GeneGo disease ontology terms associated with direct targets of upregulated miRNAs in horses with regressive versus progressive sarcoids. The y‐axis represents ranks according to the highest log of the P value in the x‐axis. BER, base excision repair; JAK/STAT, janus kinases/signal transducer and activator of transcription proteins; MAPK, mitogen‐activated protein kinase; MMR, mismatch repair; NER, nucleotide excision repair

In our study, increased concentrations of novel miRNA 20_15133, eca‐miR‐125a, eca‐miR‐24, and eca‐miR‐1249 are associated with more favorable outcome of ES disease (regression). Plasma concentrations of miR‐24 and miR‐125a were associated with outcome in various cancers of humans and turned out to be valuable prognostic biomarkers (Supporting Information Table S7). According to the network analysis, miR‐24 and miR‐125a seem to be major effector miRNAs. They are predicted to regulate the expression of E2F3, which interacts directly with the retinoblastoma protein (pRB) to control the expression of genes involved in cell cycle regulation according to GeneCards human gene database, and of STAT3, which is aberrantly activated in many cancers in humans40 (Figure 6).

Figure 6

Network for upregulated miRNA genes in horses with regressive versus progressive sarcoids. According to network analysis using GeneGo miR‐24 and miR‐125a‐5p are predicted to be major effector miRNAs in the group of upregulated miRNAs in the regression group by controlling expression of E2F3 and STAT3

Literature research of deregulated miRNAs

A literature research of the 14 dysregulated miRNAs is depicted in Supporting Information Tables S7 and S8. Supporting Information Table S7 summarizes the current scientific knowledge of the 14 dysregulated miRNAs in cancers of humans and horses, their possible oncogenic and tumor‐suppressive properties depending on tissue and tumor type, and their suitability as noninvasive diagnostic and prognostic biomarkers based on selected publications. Supporting Information Table S8 summarizes which miRNAs are part of the miRNA cluster on equine chromosome 24 homologous to the human 14q32 cluster and depicts their deregulation in different cancer types in humans.

DISCUSSION

Whole blood miRNA expression profiles of horses with ES regression resembled those of control horses and were clearly distinct from those of the progression group (Figure 2). Fourteen miRNAs were differentially expressed in horses with ES regression versus progression and are potential prognostic biomarkers that may allow noninvasive prediction of clinical outcome in ES‐affected horses using only a blood sample. These miRNAs may possess oncogenic or tumor‐suppressing properties in a tissue‐dependent manner (Supporting Information Table S7), and their biological functions as circulating miRNAs need to be further elucidated.

Among the set of differentially expressed and potentially prognostic miRNAs, 7 miRNAs (eca‐miR‐127, eca‐miR‐134, eca‐miR‐323‐5p, eca‐miR‐379, eca‐miR‐381, eca‐miR‐382, and eca‐miR‐432) belong to members of the human 14q32 miRNA network and all were downregulated in whole blood of horses with ES regression. In human beings, the 14q32 locus harbors a large bipartite miRNA aggregate which is recognized as the largest identified miRNA cluster to date comprising more than 50 members.41 It frequently is silenced in cancer41, 42, 43, 44, 45, 46, 47, 48 but can be upregulated in certain tumor types.49, 50, 51 The miRNA network on the 14q32 locus has been implicated in a variety of cancer types in humans and deregulations have been linked to clinical outcome in osteosarcoma,47 ovarian cancer,48 lung adenocarcinoma,49 and melanoma.41 In the horse, the genomic organization of the miRNA genes corresponding to the human 14q32 locus is not perfectly conserved. However, 3 corresponding polycistronic miRNA clusters on equine chromosome 24 encode all of the 7 previously mentioned equine miRNAs.39 In total, equine chromosome 24 harbors 40 miRNA genes in 4 polycistronic clusters.39 Of the 56 miRNA genes in the human 14q32 cluster,41 37 homologous equine miRNAs are found in the equine chromosome 24 region. Thus, we consider the miRNA cluster in the equine chr24: 42937011‐42745108 region homologous to the human 14q32 miRNA cluster. Interestingly, human chromosome 14q32 is 1 of the most common HPV integration sites in cervical cancer. The corresponding miRNA cluster on equine chromosome 24 might represent a predilection site for BPV integration in ES disease.52 The differential expression of miRNAs belonging to the miRNA cluster on chromosome 24 in whole blood of horses with sarcoids also might reflect a dysregulation at the tumor tissue level: Eca‐miR‐381 and eca‐miR‐134 both are upregulated in sarcoid tissue of horses requiring treatment.20 Ours is the first study to report an association of dysregulation of an miRNA cluster on chromosome 24 with cancer in horses.

Some of the differentially expressed miRNAs in our study have been linked to HPV‐induced carcinogenesis in cancers of humans, as also is predicted by pathway analysis: miR‐381 belongs to HPV core miRNAs and is downregulated in cervical squamous cell carcinoma and HPV‐associated head and neck squamous cell carcinoma.53, 54 In human oropharyngeal carcinoma, HPV positivity is associated with miR‐107 downregulation, possibly because of HPV‐enhanced immune response, whereas overexpression of miR‐107 was associated with decreased overall and disease‐free survival.23 The main HPV oncoproteins, E6 and E7, can strongly influence miRNA expression in HPV‐induced cancers: miR‐24 is upregulated by E6 and E7 expression and may promote cell proliferation by targeting the cell cycle inhibitor p27.55 In contrast, miR‐1249 is downregulated in HPV‐positive human kerationcytes in response to HPV16 E6/E7 oncoprotein expression.56 It needs to be further evaluated whether these mechanisms of miRNA dysregulation also apply to BPV oncoproteins in ES disease. Recently, miRNA deregulation has been attributed to the presence of BPV genomes in BPV‐transformed equine fibroblasts.19

According to network analysis, eca‐miR‐24 and eca‐miR‐125a‐5p turned out to be central regulatory elements, presumably by regulating the expression of E2F3 and STAT3, and their upregulation has been associated with advantageous outcome of ES disease. Both miRNAs already have been shown to be dysregulated at tissue or cell culture level: eca‐miR‐24 has been found to be downregulated in BPV‐transformed equine fibroblasts.19 Eca‐miR‐125a‐5p is highly expressed in sarcoid tissue and downregulated in serum of horses with aggressive, fibroblastic ES lesions (unpublished data—manuscript in preparation).

In our study, miRNA fingerprints in horses with ES regression versus progression were evaluated in whole blood. Typically, serum is the preferred body fluid for noninvasive miRNA profiling.57 The value, advantages, and limitations of whole blood samples as input material for miRNA studies are controversial, because changes in blood cell counts can strongly influence miRNA expression (such as to an immune response against neoplastic disease).58 Other studies, however, highlight the potential of whole blood compared to serum: whole blood contains higher miRNA concentrations59 and important sources of pre‐analytical errors related to serum processing, such as hemolysis, can be avoided.60 Furthermore, tumor‐specific exosomal miRNA might only be detectable in serum if a well‐perfused neoplastic mass of substantial size is present, whereas cellular blood cell‐derived miRNA fingerprints may already be detectable in sufficient amounts at earlier stages of tumor development,61 and with masses that are poorly perfused, such as ES lesions. Moreover, sarcoids in horses are not accompanied by changes in blood cell counts unless marked superficial ulceration is present, which was not the case in any of the initially mild ES lesions in our study. Thus, the use of whole blood as input material for miRNA analysis seemed superior to serum, particularly in light of the presence of only small and mild ES lesions in the regression and progression groups at the time of sample collection.

As for any retrospective study, our study has some limitations. Diagnosis of ES in our study was clinical, even though a definitive diagnosis of ES requires histopathology. However, clinical diagnosis has good reliability and because trauma of any nature carries the risk of exacerbating ES,10, 12 biopsies were not performed. Furthermore, surgical excision of the relatively mild lesions in our study (Supporting Information Table S1) was not indicated in any of the horses at the time of first presentation.

Because of the retrospective and longitudinal nature of our study, blood samples of different quality collected at different time points throughout the study period had to be utilized. Long‐term stored whole blood samples collected at the beginning of the study period were used in the ES‐affected horses, in contrast to the short‐term stored samples from control horses collected by the end of the study period. This difference in sample quality explains the higher RIN values in the control group. Still, we considered having a reliable control group consisting of daily monitored clinic‐owned horses, which did not develop sarcoids at any point throughout the 5‐year study period, more important than including a control group from the previously described Franches‐Montagnes study population,9 from which a blood sample collected at the beginning of the study period would have been available, which was however only examined on 2 occasions throughout the 5‐7 year study period. In the latter group, smaller transiently occurring sarcoid lesions might have been easily missed and could have substantially biased the results of the study. Comparison of miRNA differential expression in blood of horses with ES collected at the beginning, at several time points during, and at the end of the study period might have shed further light on changes of miRNA fingerprints during the course of ES disease and might have allowed us to draw additional conclusions regarding their possible biological functions. Further longitudinal studies therefore should investigate miRNA differential expression repeatedly and closely during the clinical course of ES.

Some reports state that miRNA quantification is influenced by overall RNA integrity,62 and that RNA degradation compromises the reliability of miRNA analysis.63 More recent studies, however, show that miRNA profiles themselves are not affected by the extent of RNA degradation.64, 65 In our study, we could confirm that RIN did not have any impact on differential gene expression analysis. Thus, our study shows that despite the compromised quality of a subset of blood samples because of long‐term storage, a robust miRNA study still could be conducted. This observation encourages the use of valuable blood samples from bioarchives for future longitudinal studies.

The mean age of the control horses was significantly higher than the mean age of horses in the ES groups. Sarcoids often develop in young horses,4 thus only ES‐free middle‐aged and old horses (≥12 years) were included in the control group. We did not want to include younger age‐matched controls because of the risk of developing ES lesions and becoming false‐negative controls. Age may affect miRNA expression profiles in human beings.66, 67 Its effect on the results of the present study was considered to be less important than the need for a reliable control group.

Our study provides a set of 14 blood‐derived miRNAs significantly associated with prognosis of ES disease in a limited number of samples with markedly favorable (regression) versus aggressive (progression) clinical courses. These miRNA fingerprints warrant further validation as prognostic biomarkers in a larger study cohort. Admittedly, to date, miRNA differential expression analysis remains an elaborate diagnostic technique and for now may not serve as a practical test and prognostic biomarker in equine practice. However, in the future, qRT‐PCR may replace next generation sequencing to assess miRNA expression levels, which would render the procedure less cost prohibitive and time consuming. Additional studies focusing on correlation of whole blood miRNA expression and miRNA and target gene expression at the tissue level as well as monitoring of changes of miRNA expression over time may shed more light on the speculated biological functions of these blood cell‐derived miRNAs, particularly, if dysregulation of the miRNA cluster on equine chromosome 24 by integration of the BPV genome is a possible mechanism that drives disease progression. Regardless of the biological function of blood‐derived miRNAs, these potential biomarkers might offer novel perspectives in the clinical setting, such as for prepurchase examinations or therapeutic decision‐making.

CONFLICT OF INTEREST DECLARATION

Authors declare no conflict of interest.

OFF‐LABEL ANTIMICROBIAL DECLARATION

Authors declare no off‐label use of antimicrobials.

INSTITUTIONAL ANIMAL CARE AND USE COMMITTEE (IACUC) OR OTHER APPROVAL DECLARATION

This study was fully approved by the Ethical Committee of the Canton of Bern (VD 2227.0, VD 2227.1 und VD 2227.2, BE110/15).

HUMAN ETHICS APPROVAL DECLARATION

Authors declare human ethics approval was not needed for this study.

Supporting Information Table S1 Overview of the study cohort and RNA quality and quantity after extraction from the individual whole blood samples. Minimal RNA integrity number (RIN) is 1, maximal RIN is 10. Abbreviations: ID = identification number, f = female, m = male, mc = male castrated, RGR = regression, PRG = progression, CTL = control, RIN = RNA integrity number

Click here for additional data file.

Supplemental Table 2 MicroRNA sequencing read summary

Click here for additional data file.

Supplemental Table 3 MiRDeep2 output table. The table includes the miRDeep2 score, the number of novel miRNAs reported by miRDeep2, the estimated false positives, estimated true positives, the number of known miRNAs in horses, the number of known miRNAs in the dataset, the number and percentage of known miRNAs detected by miRDeep2, the estimated signal‐to‐noise ratio and the excision gearing

Click here for additional data file.

Supplemental Table 4 Novel miRNAs predicted by miRDeep2

Click here for additional data file.

Supplemental Table 5 Target gene prediction for the 14 differentially expressed miRNAs in the regression vs progression group including the total number of gene targets per miRNA and a list of all target genes according to TarBase v7.0 (Diana Tools) with a prediction score of ≥0.9

Click here for additional data file.

Supplemental Table 6 Singular enrichment analysis for target genes of the differentially expressed miRNAs in the regression vs progression group [of gene ontology (GO) biological process, GO molecular function and KEGG pathways according to GeneCodis]

Click here for additional data file.

Supplemental Table 7 Literature research (selected papers) and current state of knowledge regarding the biological function of the 14 differentially expressed equine miRNAs in the regression vs progression group and their human homologues as well as their validity as diagnostic and prognostic biomarkers in different types of cancers, including human HPV‐ and equine BPV‐induced cancers or cancer models. Abbreviations: HN SCC = head and neck squamous cell carcinoma, O(P) SCC = oro(pharyngeal) squamous cell carcinoma, LSCC = laryngeal squamous cell carcinoma, CC = cervical carcinoma, SNP = single nucleotide polymorphisms, HPV+ = HPV positive, BC = breast cancer, ES = equine sarcoid, CRC = colorectal cancer, RCC = renal cell carcinoma, L(A)C = lung (adeno)carcinoma, HNSCC = head and neck squamous cell carcinoma, EPS = epitheloid sarcoma, OS = osteosarcoma, OC = ovarian cancer, NSCLC = non‐small cell lung cancer, HCC = hepatocellular carcinoma, HBV HCC = HBV related hepatocellular carcinoma, GC = gastric cancer, ESCC = esophageal squamous cell carcinoma, MTC = medullary thyroid carcinoma, BL = Burkitt lymphoma, PC = prostate cancer, MPM = malignant pleural mesothelioma, PaC = pancreatic cancer, PTC = papillary thyroid carcinoma, GL = glioma, LC = lung cancer, MM = multiple myeloma, MEN1 = multiple endocrine neoplasia type 1, NPC = nasopharyngeal carcinoma

Click here for additional data file.

Supplemental Table 8 Overview of differentially expressed miRNAs in whole blood of horses with sarcoid regression vs progression, whose human homologues belong to the 14q32 miRNA network and whose dysregulation is associated with cancer. The human 14q32 miRNA cluster corresponds to three polycistronic miRNA clusters on chromosome 24 in the equine species. Abbreviations: ME = melanoma, OS = osteosarcoma, GIST = gastrointestinal stromal tumor, BLC = bladder cancer, LAC = lung adenocarcinoma, OC = ovarian cancer

Click here for additional data file.

Supplemental Figure 1 Principal component analysis (PCA) clustering of samples based on expression level of top 200 genes. ES_Numbers represent the different individuals. Abbreviations: CTL = control, PRG = progression, RGR = regression, PC = principal component

Click here for additional data file.

Supplemental Figure 2 Principal component analysis (PCA) clustering of samples after removal of outliers based on expression level of top 200 genes. ES_Numbers represent the different individuals. Abbreviations: CTL = control, PRG = progression, RGR = regression, PC = principal component

Click here for additional data file.

Supplemental Figure 3 Clustered and scaled expression patterns for the top 100 miRNAs ordered by adjusted P‐values. Boxplots are shown for the expression patterns of each miRNA within the group to give a better idea of how well the groupings fit the expression data. Abbreviations: CTL = control, PRG = progression, RGR = regression

Click here for additional data file.

Supplemental File 1 Code used for miRNA discovery and expression

Click here for additional data file.