Comprehensive characteristics of microRNA expression profile of equine sarcoids

Biochimie 137 (2017) 20e28

Contents lists available at ScienceDirect

Biochimie journal homepage: www.elsevier.com/locate/biochi

Research paper

Comprehensive characteristics of microRNA expression profile of equine sarcoids Klaudia Pawlina a, *, Artur Gurgul a, Tomasz Szmatoła a, Christoph Koch b, €hlmann b, c, Maciej Witkowski d, Monika Bugno-Poniewierska a Kathrin Ma a Laboratory of Genomics, Department of Animal Genomics and Molecular Biology, National Research Institute of Animal Production, Krakowska 1, 32-083, Balice, Poland b €nggassstrasse 124c, Postfach 8466, CH-3001, Swiss Institute of Equine Medicine ISME, Faculty of Veterinary Medicine, University of Bern and Agroscope, La Bern, Switzerland c Equine Clinic: Surgery and Radiology, Department of Veterinary Medicine, Free University of Berlin, Oertzenweg 19b, 14163, Berlin, Germany d Department of Animal Reproduction, Institute of Veterinary Medicine, University of Agriculture, Al. Mickiewicza 21, 31-120, Krakow, Poland

a r t i c l e i n f o

a b s t r a c t

Article history: Received 13 January 2017 Accepted 28 February 2017 Available online 1 March 2017

Equine sarcoids are the most common neoplasms occurring in horses. Despite frequent occurrence, they are still not well described at the molecular level. Thus, in the present study, we performed a comprehensive comparative analysis of sarcoid miRNAome profile to identify aberrantly expressed microRNAs, along with their structural variants, potentially useful as biomarkers and, in a wider perspective, broaden the knowledge about this tumor and underlying mechanisms. To this end, we conducted next generation sequencing and as a result we identified both known and potentially novel miRNAs. Differential expression analysis revealed the existence of almost one hundred miRNAs being over- or underexpressed in sarcoids in comparison to healthy tissue (p-adj<0.05), of which many are known for their involvement in processes crucial for neoplastic transformation. Among upregulated miRNAs there were those associated with decreased cell adhesion abilities as well as engaged in global protein production, while downregulation of some miRs i.a. increased cell expansion abilities. Moreover, we identified altered expression levels of miRNA variants (isomiRs) between the investigated tissues. Further analysis revealed that 50 isomiRs comprise different seed sequences leading to target gene switching followed by activation of different biological pathways. Our results are the first which revealed the complexity of microRNA profiles in equine sarcoids and skin tissue, along with the dynamism of their growing in importance concomitants, namely isomiRs. They also showed miRNA molecules and biological pathways important from the sarcoid oncogenesis point of view. © 2017 Published by Elsevier B.V.

Keywords: microRNA isomiR Sarcoid Horse NGS

1. Introduction microRNAs, belonging to the class of small non-coding RNAs [1], have an ability to orchestrate gene expression through binding to gene transcripts. It further leads to translational repression, degradation [2] or in some cases up regulation of target gene expression [3].

* Corresponding author. E-mail addresses: [email protected] (K. Pawlina), artur.gurgul@ izoo.krakow.pl (A. Gurgul), [email protected] (T. Szmatoła), [email protected] (C. Koch), [email protected] (K. M€ ahlmann), [email protected] (M. Witkowski), monika.bugno@izoo. krakow.pl (M. Bugno-Poniewierska). http://dx.doi.org/10.1016/j.biochi.2017.02.017 0300-9084/© 2017 Published by Elsevier B.V.

Numerous miRNA-focused studies [4e7] showed their importance and engagement in the regulation of a wide range of vital biological processes such as development, metabolism or apoptosis [8,9]. Considering that these processes are essential for proper cell functioning, it is unsurprising that alterations of microRNA expression profiles were identified in the course of different diseases, such as cancer [10,11]. In the carcinogenesis process, miRNAs can function as oncogenes (e.g. miR-155) or tumor suppressors (e.g. miR-15a) [12]. A unique combination of oncogenic and tumor suppressor microRNA expression profiles constitutes tumor signature [8,12]. These signature microRNAs may serve as potential biomarkers for early diagnosis, staging, prognosis as well as response to treatment [12,13]. For example, serum levels of miR-21 were correlated with diagnosis and prognosis of colorectal cancer

K. Pawlina et al. / Biochimie 137 (2017) 20e28

[14]. Despite a significant improvement in the field of identification of new miRNA molecules in recent years, the coverage of microRNA profiles of livestock species is still lacking. The horse can be a good example with its 690 records in miRBase Release 21.0 [15,16]. When it comes to characterization of horse miRNAome profile in pathological states, the data are even more scarce. Research conducted so far concerned osteochondrosis [17] and bovine papillomavirus type 1 transformed equine cells [18]. The latter studies were focused on a cell culture model of sarcoid tumor, which is the most common skin cancer occurring in horses. It is estimated that it may constitute as much as 40% of all diagnosed cancers [19e21]. Sarcoids are classified as locally malignant, benign fibroblastic tumors of fibrous tissue [22], differing in aggressiveness [23e25]. Despite the fact that they do not metastasize, they may recur if treatment is not carried out properly [26]. In sporadic cases, complete and constant remission is observed [27] but, on the other hand, protracted, unsuccessful treatment may lead to animal euthanasia. BPV (Bovine Papillomavirus) is believed to be the causative agent, however, the mechanism of infection and its oncogenic activity have not been fully understood yet [22,28,29]. Since horses are of great importance in livestock production, it is vital to acquire knowledge about mechanisms which orchestrate gene expression and may be disrupted in the course of carcinogenesis. Therefore, taking into account the fact that different types of cancers exhibit perturbed expression profiles of microRNAs which are involved in the oncogenesis process, we hypothesized that sarcoid cells are also characterized by altered expression of oncomiRs and their structural variants. We applied next generation sequencing along with a subsequent bioinformatical analysis. The obtained results will shed light on the horse miRNAome profile and its disruptions caused by neoplastic transformation to sarcoid. This will further facilitate understanding of microRNA functioning and mechanisms engaged in the oncogenesis process occurring in sarcoids. 2. Materials and methods 2.1. Small RNA library preparation Tissue samples from equine sarcoids as well as healthy skin tissue from five horses were collected in Equine Clinic in Bern, _ Switzerland as well as in Horse Clinic Słuzewiec in Warsaw, Poland. The samples were provided from clinical patients treated for equine sarcoid disease in those clinics. Sample collection for this study was approved by the Animal Experimentation Ethics Committee of the Canton of Berne, Switzerland (BE 30/11; 11 April 2011). For the _ samples obtained in Horse Clinic Słuzewiec in Warsaw, Poland, Part 1, Division 1.2, Paragraph 1 of “Legislation for the protection of animals used for scientific or educational purposes” of Poland, stating that no Ethic Committee approval is needed when providing veterinary services, was applied. Moreover, all efforts were made to minimize suffering of animals. Immediately after excision, the samples were immersed in RNAlater Stabilization Solution (Thermo Fisher Scientific) and stored in 20  C. Total RNA was isolated using Direct-zol RNA kit (Zymo Research) according to the manufacturer protocol. The RNA quality control was carried out with the use of 2200 TapeStation instrument (Agilent). Then, one mg of total RNA was used to construct libraries from each sample following the TruSeq Small RNA Sample Preparation Kit protocol (Illumina, protocol 15004197 Rev.F, February 2014). In brief, the samples were ligated with 30 and 50 adapters, respectively. Then, RNA-adapter ligation products were reverse transcribed and PCR amplified using one of 12 different indexed primers for each sample. Next, the libraries were gel-purified on Novex 6% TBE PAGE gel (Invitrogen).

21

The quality of the obtained libraries was checked using 2200 TapeStation instrument (Agilent) whereas concentration was assessed using Qubit 2.0 Fluorometer (Life Technologies). Furthermore, the libraries were mixed with the PhiX control library (Illumina). Then, they were clustered on Illumina Flowcell_v3 in cBot cluster station and sequenced with the use of HiScanSQ (Illumina) according to the manufacturer protocol, applying 36 cycles. 2.2. Data analysis The resulting reads were analyzed with the use of UEA sRNA Workbench V3.2. [30]. First, 30 adapter sequences were trimmed off and then the obtained sequences were filtered using the following parameters: 17e35 nt in length, minimum abundance of at least 3 supporting reads, tRNA and rRNA sequences excluded from the data. Then, miRNA sequences were identified applying miRCat tool with default parameters for animals, except for: minimum abundance (3 or greater), minimum length (17 nt) and maximum length (25 nt). The obtained miRNAs were then checked for the concordance with the miRBase 21.0. Next, to profile miRNAs' abundance, miRProf tool was used applying the following parameters: minimum length of 17 nt, maximum length of 25 nt. The potentially new sequences were then checked for the presence of other sRNA species, such as snoRNAs, using RNAcentral database (Release 5) [31]. 2.3. Differential expression analysis To characterize differential expression (DE) of isomiRs, first, we prepared a library of reference sequences on the basis of miRBase 21.0. Then, we assigned a unique number to each isomiR of a miRNA and carried out the DE analysis using two algorithms applying different methodological approaches, namely: DESeq2 (a parametric method based on the negative binomial distribution) [32] and SAMSeq (a nonparametric method based on the Wilcoxon test) [33]. The paired tests were used for both approaches. For further analyses, only miRNAs showing statistically significant differences between the tested samples (p  0.05 for both algorithms, corrected for multiple testing) were chosen to augment the likelihood of the estimation. 2.4. Pathway analysis In order to identify target genes and pathways potentially modulated by the differentially expressed microRNAs, the mirPath DIANA Tools web application [34] was used with DIANA e TarBase v7.0 as a reference database of target genes. When it comes to the analysis of pathways in which potentially new miRNA are engaged, first, we employed MR-microT software [35,36] to predict target genes for four chosen novel miRNAs. Then, the target genes were subjected to pathway analysis using KEGG Pathway Database [37]. 2.5. IsomiRs' analysis First, the analysis of alternative seed sequences was carried out. 50 isomiRs of all miRNAs detected in the study were extracted from the obtained data, and then prediction of their target genes using TargetRank software was performed [38]. Pathway analysis was carried out with the use of KEGG Pathway Database [37]. Next, the analysis of modification types of isomiR sequences was conducted. Estimations of occurrences of particular nucleotides at both ends of isomiRs were carried out and the obtained results were then subjected to Chi-square test.

22

K. Pawlina et al. / Biochimie 137 (2017) 20e28

2.6. Validation of the results The expression levels of 10 selected microRNAs were validated with the use of RT-qPCR method. The validation was conducted with the use of eight commercially available as well as two custom synthesized LNA primer sets (Locked Nucleic Acids) (Exiqon), namely NA1 and NA2 e complementary to potentially new miRNAs identified in this study (S1 Table). The reverse transcription reactions were carried out using Universal cDNA synthesis kit II (Exiqon), according to the manufacturer protocol. Subsequent quantitative PCR reactions were conducted with the use of ExiLENT SYBR Green master mix (Exiqon), according to the protocol. Each reaction was run in triplicate, and for each pair of primers NTC (Non-Template Control) was also included. Moreover, to exclude “plate effect”, on every plate an additional qPCR reaction with „spike in” was included, as specified in the manufacturer protocol. The amplification was carried out on Eco Real-Time PCR System (Illumina), including melting curve analysis. The efficiency of amplification reactions was estimated for each pair of primers using standard curve method. Obtained data were analyzed with the use of NormFinder software [39], to choose the most stable miRNAs which could serve as endogenous control. In order to evaluate relative expression levels DDCt method including reaction efficiency (E) [40] was employed. 3. Results 3.1. MiRNAome profile of sarcoid tissue To identify new miRNAs and characterize miRNA expression profile of equine sarcoids as well as healthy skin tissue, we applied next generation sequencing using Illumina HiScanSQ platform. As a result, we obtained ~124.7 million sequence reads from all constructed libraries, encompassing ~4.8 million unique sequences (GEO accession number GSE87901). Then, the sequences were mapped to the reference genome (Ensembl e EquCab2.0) and compared to the miRBase (21.0), using the miRCat algorithm. The analysis revealed the expression of 249 distinct known miRNAs, 49 microRNA* complementary to microRNAs deposited so far in the miRBase 21.0 and 523 unique potentially novel microRNAs, of which 81 (~16%) were classified as other sRNA species according to the RNAcentral database. The great majority of the identified microRNAs were 22 bp in length, followed by 21 and 23 bp. Around 100 of potentially new miRNA molecules were detected solely in the sarcoid tissue samples. For some microRNAs, one or a few miRNA* sequences were identified as well as genomic duplications of precursor sequences and alternative precursors (S3 Table). 3.2. MicroRNAs are deregulated in equine sarcoid tissue Firstly, we conducted a cluster analysis on the basis of the expression profiles of the investigated samples using DESeq2 software. This analysis revealed grouping of four sarcoid samples with the healthy tissue samples, which is why they were excluded from further analysis. Then, we repeated clustering and as a result we distinguished two distinct clusters, namely the cluster of the sarcoid samples and the cluster of the control samples (Fig. 1). Next, we carried out Differential Expression (DE) analysis applying DESeq2 and SAMseq algorithms. As a result, we identified 40 upregulated miRNA sequences (consisting of 108 isomiRs) and 34 downregulated microRNAs (including 121 isomiRs) (S4 Table). For some miRNAs, DE was observed for only one isomiR, whereas in the case of other microRNAs, the level of expression varied for up to 12 isomiRs (miR-199a-3p) (S4 Table). Moreover, all of differentially

Fig. 1. Heatmap of the sample-to-sample distances, on the basis of the microRNA next generation sequencing expression data. S e sarcoid samples, C e control samples.

expressed isomiRs stemming from one microRNA showed the same direction of expression changes between the examined tissues. Ten microRNAs showing the greatest differences in the expression levels between the normal tissue and the sarcoids for each used algorithm are depicted in Table 1. Log2FC e the binary logarithm of the Fold Change parameter. IsomiRs are denoted as “is.X”; where “X” means subsequent isomiRs. Letter „m” means that an isomiR was the major one, among other isomiRs of the same microRNA, in a sample. isomiRs with additional nucleotides in reference to the canonical miRNA sequence or having an internal substitution (maximally 1 nucleotide) were marked with “(1)”. “NA” stands for potentially new microRNAs, detected in this study. Bolded letters mean that given microRNA was identified by both algorithms. 3.3. Deregulated miRNAs are engaged in cancer-related pathways In order to shed light on roles of the identified DE miRNAs in cellular functioning, we carried out the pathway analysis on the basis of microRNA target genes. As a result, we identified 86 pathways (S5 Table). The largest number of microRNAs was involved in “microRNAs in cancer” pathway, whereas the largest number of target genes was engaged in “pathways in cancer”. Moreover, other pathways associated with carcinogenesis or viral infection were identified including PI3K-Akt and MAPK signaling pathways, viral carcinogenesis, transcriptional misregulation in cancer and hepatitis B. Among the identified potentially new miRNAs were those expressed only in sarcoid tissue. To elucidate their roles in tumor cells we investigated the engagement of a few of them in biological pathways. The analysis showed that they are engaged mainly in cancer related pathways, such as MAPK, PI3-Akt, Ras signaling pathways, miRNAs in cancer, transcriptional misregulation in cancer, viral infection as well as different pathways related to signaling and cell interactions (S6 Table). 3.4. MicroRNA isoforms are abundantly expressed We identified over 1000 isomiRs which are miRNA variants

K. Pawlina et al. / Biochimie 137 (2017) 20e28

23

Table 1 Ten statistically significantly differentially expressed (p < 0.05) microRNAs, along with their isomiRs, characterized by the greatest differences in expression levels between the sarcoid and control samples, identified by the two applied algorithms. DESeq2

SAMSeq

upregulated

Log2FC

downregulated

Log2FC

upregulated

Log2FC

downregulated

Log2FC

miR-450a-is.3m miR-450a-is.5 miR-450a-is.4 miR-136-3p*(NA) miR-409-3p-is.(1)5 miR-381-is.3 miR-127-is.2 miR-409-3p-is.(1)4 miR-409-3p-is.7 miR-134-is.2 miR-299-is.(1)2 miR-450a-is.(1)1 miR-450a-is.(1)2 miR-432-is.2 miR-409-3p-is.6 miR-381-is.(1)1 miR-134-is.3m miR-370-is.5m miR-432-is.5 miR-214-is.3

8.51 8.36 7.92 7.88 7.36 6.89 6.78 6.45 6.34 5.98 5.95 5.95 5.77 5.62 5.29 5.26 5.26 5.15 5.11 5.10

miR-200c-is.4 miR-200a-is.5m NA/hsa-miR-375 miR-141-is.5 miR-141-is.(1)2 miR-205-is.8 miR-200b-is.(1)1 miR-200b-is.5 miR-141-is.(1)6 miR-200a-is.4 miR-141-is.2 miR-200a-is.(1)2m miR-200b-is.6 miR-200a-is.2 miR-141-is.(1)4 miR-205-is.(1)5 miR-205-is.11m miR-205-is.10 miR-141-is.4m miR-141-is.1 miR-96-is.3m miR-200b-is.3 miR-205-is.9 miR-708-is.5 miR-200b-is.(1)3m miR-429-is.3m miR-182-is.3

10.57 10.36 10.36 9.63 9.44 9.37 9.25 9.09 9.03 9.03 8.78 8.77 8.70 8.68 8.68 8.66 8.42 8.37 8.34 8.12 8.07 7.9 7.89 7.82 7.73 7.57 7.54

miR-136-3p*(NA) miR-450a-is.3m miR-450a-is.5 miR-450a-is.4 miR-450a-is.(1)1 miR-432-is.5 miR-299-is.4m miR-450a-is.(1)2 miR-127-is.2 miR-494-is.(1)m miR-381-is.3 miR-409-3p-is.7 miR-758-is.2m miR-214-is.3

36.69 34.22 33.87 33.53 32.99 32.38 32.28 32.18 32.14 31.95 31.84 31.82 31.78 31.69

miR-205-is.10 miR-200c-is.5m miR-200a-is.5m miR-200a-is.(1)2m miR-429-is.3m miR-205-is.9 miR-200b-is.6 miR-200c-is.4 miR-141-is.(1)6 miR-200a-is.4 miR-429-is.(1)4 miR-141-is.5 miR-141-is.(1)2 miR-205-is.8 miR-200b-is.(1)1 miR-200b-is.5 miR-708-is.4m miR-141-is.2 miR-141-is.(1)4 miR-200a-is.2 miR-96-is.3m miR-141-is.1 miR-200b-is.3 miR-95-is.1m miR-205-is.7 miR-429-is.2 miR-200c-is.1 miR-708-is.5 miR-200c-is.3 miR-182-is.3

39.13 36.63 36.49 36.47 36.38 36.36 36.21 35.69 35.66 35.66 35.13 34.98 34.56 34.35 34.19 34.04 33.80 33.76 33.59 33.59 33.46 32.89 32.87 32.68 32.66 32.58 32.57 32.49 32.46 32.23

differing in sequence length (30 or/and 50 end) and/or having an internal nucleotide polymorphism in respect to the miRBase 21.0 deposited reference particle. For 142 microRNAs, it was an isomiR particle which had the highest expression level. For single samples, the number of isomiRs with the highest expression levels ranged from 41 to 88, which constituted from 25% to 49% of the detected known microRNAs. For 12.5% of the identified known microRNAs, no reference sequence, deposited in the miRBase 21.0 miRNA, was detected, solely its isomiRs. The most common isomiR type in our data was an alteration of the 30 end (664 isomiRs), then the 50 end modifications (522 isomiRs), and the least frequent was an internal nucleotide polymorphism (27 isomiRs).

3.5. IsomiRs comprise different seed sequences leading to target gene switching Alterations in microRNA sequence may lead to changes in the seed region which is responsible for target gene recognition. In our research, one seed sequence was expressed for almost half of the identified miRNAs while the rest of them had two or even more alternative seed sequences. However, the increase of the number of alternative seeds was reflected in the decrease of the number of miRNAs (Fig. 2). IsomiRs of only three miRNAs comprised the highest number of seed regions, namely eight (Table 2). Next, to shed light on the roles and significance of different seed sequences, we carried out a computational analysis to identify target genes specific for alternative seeds. The analysis showed that the number of target genes for each alternative seed differed considerably (Table 2). It also revealed that alternative seeds introduced their unique target genes, not targeted if other isoform was expressed. Those unique genes constituted on average even up to 51% of genes targeted by all isomiRs of one miRNA.

To further confirm the potential of alternative seed sequences to exhibit different roles, we conducted pathway analysis on the basis of their unique target genes. It revealed that, even though most of the pathways were common, they were organized differently in the hierarchy as well as there also were pathways specific for one target gene set (S7 Table). 3.6. MicroRNA isoforms are not equally modified To elucidate the nature of isomiRs, we carried out further in silico analysis aimed at assessing the frequency of occurrences of particular nucleotides at the 30 and 50 termini and as internal substitution. The analysis showed that there are no statistically significant differences with regard to nucleotide substitution (pvalue>0.05). However, when it comes to the length isoforms, the ends are not modified equally frequently (Fig. 3) and nucleotides are not distributed evenly (Fig. 4). The differences were statistically significant at the 30 end in general. However, having taken into consideration the type of a modification, that is elimination or addition of a nucleotide, those differences were statistically significant at both ends, except for 50 addition (Table 3). Moreover, statistically significant differences were also observed for particular nucleotides (Table 3). 3.7. Validation of the results using RT-qPCR The results of differential expression analysis obtained using next generation sequencing were subjected to validation applying PCR method. Of the 10 selected microRNAs, two with the most stable expression were used as endogenous controls. For the rest, the level of expression in the sarcoid tissue in comparison to the healthy tissue was estimated. In seven of eight cases, the direction

24

K. Pawlina et al. / Biochimie 137 (2017) 20e28

2

Fig. 2. Distribution of seed regions per microRNA.

Table 2 Data on the number of genes potentially targeted by different seed sequences of microRNAs with the highest number of alternative seeds. miRNA

Number of seeds

Mean number of target genes

let-7a let-7c let-7f miR-10b miR-21 miR-22 miR-148a miR-181a miR-205

8 7 7 7 8 7 7 8 7

316 306 349 445 544 518 459 751 981

± ± ± ± ± ± ± ± ±

94 85 82 142 162 79 182 630 393

of expression changes determined with the use of SAMseq software was consistent with the results of RT-qPCR method, which means that the concordance of the results was 87.5%. The results of the RNA-Seq analyzed with the use of DESeq2 algorithm coincided with the results of RT-qPCR technique for six miRNA sequences, meaning 75% concordance. 4. Discussion The aim of our study was to investigate whole miRNAome expression profile, including both known and potentially novel sequences, in equine sarcoids and tumor-distant skin tissue to characterize miRNA activity in this tumor. Data obtained from next generation sequencing allowed to identify 249 known and 484 potentially novel miRNAs. The majority of the identified microRNAs in this study was 22 nt in length, which is the most commonly observed length of miRNAs in various animal species [41,42]. To give insight into processes taking place in sarcoid neoplastic transformation we identified microRNAs differentially expressed between the investigated samples. The analysis revealed many known oncomiRs, both up- and downregulated. Among the 10

Mean number of seed-unique target genes 100 161 173 201 195 241 222 292 283

± ± ± ± ± ± ± ± ±

61 (38%) 80 (51%) 66 (49%) 73 (46%) 67 (37%) 80 (46%) 128 (46%) 274 (39%) 216 (27%)

most overexpressed microRNAs in the sarcoid samples we identified i.a. miR-450a, which is also upregulated in a type of skin cancer, namely oral squamous cell carcinoma and is believed to regulate cell adhesion, a process essential for cancer development [11]. Taking into account the same tissue of origin of squamous cell carcinoma and sarcoid, upregulation of miR-450a could also increase cell detachment abilities in sarcoid tissue. The influence of miRNAs on sarcoid carcinogenesis may be exerted by other upregulated oncomiRs, such as miR-10a and 10b. Those miRNAs, whose aberrant expression was also identified in liver cancer, colorectal cancer and leukemia, target genes which play roles in translation and its regulation [43]. Therefore, miR-10a and 10b have influence on the global cell protein production, which makes them important factors of the oncogenesis process enabling the domination and expansion of cancer phenotype [43]. The obtained results also suggest their involvement in the regulation of sarcoid cellular processes. Disruptions of miR expression in the sarcoid samples were also manifested by downregulation of miR-200a, miR-200b and miR141. These miRNAs are also expressed in other virus-related cancers, namely human gastric cancer caused by Epstein-Barr virus

K. Pawlina et al. / Biochimie 137 (2017) 20e28

25

Fig. 3. Occurrences of particular nucleotides at the 50 and 30 ends of isomiRs.

Fig. 4. Occurrences of different types of modifications of isomiRs' ends. “-” stands for elimination of a nucleotide; “þ” stands for addition of a nucleotide.

[44,45] and HPV-induced cervical cancer [46]. Among genes targeted by them are ZEB1 and ZEB2 encoding transcription repressors of crucial importance in maintaining the epithelial cell phenotype [47]. The ZEB-induced repression leads to epithelial-mesenchymal

transition (EMT), manifested by cell morphology disruptions, which in turn increase cell migration and invasion abilities [47]. Therefore, apart from its physiological roles EMT is believed to play a vital part in neoplastic transformation [48,49]. Downregulation of

26

K. Pawlina et al. / Biochimie 137 (2017) 20e28

Table 3 Chi-square test results with respect to the ends, modification types and nucleotides. End p-value

5′ 0.052

Type of modification p-value

5′ elimination 0.04

5′ addition 0.051

3′ elimination 0.00337

3′ addition 0.00000883

Nucleotide p-value

A 0.035

G 0.0096

U 0.00001

C 0.019

expression of miR-200 family members may be reflected in ZEBinduced EMT occurrence and, as a result, increased cell invasiveness. It can be indirectly evidenced by research on invasive breast cancer, where miR-200 family was underexpressed [47,48]. Moreover, underexpression of these miRs reported in this study also supports this idea. One might postulate that in sarcoids EMT could be one of mechanisms responsible for invasion of tumor cells into deeper layers of skin. One of factors responsible for altered expression of miR-200 family may be p53 protein, which increases expression levels of this family, thus leading to EMT inhibition [50]. It is noteworthy that in the course of BPV infection, a causative agent of sarcoids, p53 is deactivated [51,52]. It may be followed by a decrease in expression levels of miR-200 family and as a result increase the likelihood of EMT. What is more, viral proteins may themselves influence miRNA expression profiles [10], which could also take place in sarcoid tumors. To further confirm the engagement of the identified DE microRNAs in oncogenesis, we carried out the analysis of biological pathways. It revealed that the largest number of miRNA target genes belongs in the category “pathways in cancer”, while many other pathways, such as viral carcinogenesis, HTLV-I infection or p53 signaling pathway, were involved in viral infections and carcinogenesis. Moreover, pathways involving genes responsible for the regulation of vital cellular processes were identified, such as regulation of actin cytoskeleton, cell cycle, ribosome or apoptosis. Actin cytoskeleton abnormalities are commonly observed in cancer cells, due to its importance in maintaining cell integrity. They may disrupt cell morphology and lead to the acquisition of excessive migratory characteristics, which promotes cancer progression [53]. The corrupted machinery of cancer cells also influences cell cycle, ribosomal processes and apoptosis, owing to which cells gain uncontrolled proliferation and survival abilities. Altogether, the results stay in agreement with the nature of sarcoid that is a neoplasm caused by a viral infection and confirm the engagement and importance of microRNAs in sarcoid carcinogenesis. Among the identified microRNAs, there also were previously unobserved molecules, of which over 100 were detected solely in the tumor tissue. The analysis of biological pathways revealed the occurrence of cancer related pathways (e.g. MAPK, PI3-Akt) as well as those connected with miRNAs or viral infections. It may imply that those miRs are characteristic for this type of cancer and may be synthesized during its progression to promote and support the invasion and expansion of sarcoid cancer cells. However, in order to verify it and obtain empirical evidence, further research is necessary. MicroRNAs may be expressed as length or sequence variants, so called isomiRs, which are a common phenomenon [54e59]. They may be the most prevalent sequences in up to almost 50% of the identified miRNA molecules [57,58]. Our results are no exception, isomiRs were identified in both tumor and normal skin tissue and constituted from 25% to 49% of the most highly expressed microRNAs in single samples, in reference to the miRBase deposited canonical sequence. So far, no such estimations have been made for

3′ 0.00259

animal cancers, but our results corroborate the human cancer data [57,58] and reveal another aspect of miRNA conservation in mammal species. Most of the identified isomiRs were altered in length at the 30 end in comparison to the canonical sequence, which is in agreement with other research [58e60]. Alterations of the 50 end length are less frequent because they may affect affinity for target genes, caused by the displacement of the active region and changes in the thermodynamic stability of miRNA molecules [60]. Consequences may be substantial, namely alterations in the set of targeted genes. Thus, we paid particular attention to 50 isomiRs and decided to test their potential for seed region changes. The analysis revealed that almost half of the analyzed sequences had one seed region, but the rest of them had up to eight alternative seeds, which corroborates results obtained by Wojcicka et colleagues [57]. Having in mind possible target gene switching triggered by seed conversion, we predicted genes potentially targeted by alternative seeds of a few miRs. Obtained results showed high variability in the number of target genes for single seed regions with the standard deviation reaching as many as 630 target genes for miR-181a. What is more, estimated in the next step percentage of seed-unique target genes ranged from 27% to 51%, which means that even half of target genes identified for one isomiR may be regulated solely by that particular isoform. This illustrates how significant and powerful may be even single nucleotide changes in miRNA sequences. To further confirm the importance and different roles of alternative seed regions, we carried out the pathway analysis. On the one hand it showed many common pathways, but, on the other hand, there also were those isomiR specific. Recently conducted research on human cell cultures [60] confirmed that miRNA isoforms play distinct roles in a wide variety of cellular processes, and our results suggest that other species and tissues are no different. These observations are additionally supported by the results of the DE analysis, which showed that not all isomiRs deriving from one miRNA are equally expressed. Disrupted expression of isomiRs in the sarcoid tissue suggests that this is a dynamic process which is subjected to changes in response to various stimuli. These results are supported by research on fruitfly [54], Atlantic halibut [56] as well as in humans [60], in dependence to gender [58,61], race and population [58]. Such a strategy would allow for selective, temporary gene silencing, depending on organism needs. Thus, one may hypothesize that altered expression of miRNA isoforms in sarcoids contribute to their survival and expansion. In order to shed some light on microRNA isoforms' properties, we looked closer at isomiRs' modifications and nucleotide composition. The analysis showed that the 30 end is significantly more frequently modified, which may be explained by the fact that there is no active region there so nucleotide changes are not eliminated. Interestingly, having taken into account modification type we observed that also a 50 modification, namely elimination, became statistically significant. We hypothesize that elimination events are not avoided as much as addition events because the latter may result in seed region shifts leading to changes of targeted genes. We went further in our analysis and estimated the frequency

K. Pawlina et al. / Biochimie 137 (2017) 20e28

of occurrences of single nucleotides, which revealed substantial differences. High frequency of A and U nucleotides at the 50 end and the rare occurrence of G nucleotides was also observed by other researchers [62,63] and may be explained by different preferences of Argonaute protein [64]. Moreover, recently conducted research confirmed another possibility, namely Dicer cleavage preferences towards releasing miRNAs with U nucleotide at the 50 end and avoiding those with G nucleotide [65]. Our estimations support those results. 5. Conclusions The research reported here is the first one to give insight into the miRNAome profile of equine sarcoid and skin tissues. It elucidates the roles of microRNAs in the pathogenesis of sarcoids and sheds light on the emerging importance of single isomiRs in regulation of gene expression. It also shows how dynamic and rich is the repertoire of miRNA isoforms, which should be taken into account when analyzing causes and effects of neoplastic transformation. Moreover, one should not forget potential involvement of viral proteins in miRNA expression disruptions. Altogether, similarities with data on miRNA functions in human neoplasms suggest the existence of common, interspecies mechanisms, leading to changes in miRNA expression profiles during oncogenesis. Acknowledgments The research was funded from statutory activity of National Research Institute of Animal Production in Poland, number 04e008.1. Appendix A. Supplementary data Supplementary data related to this article can be found at http:// dx.doi.org/10.1016/j.biochi.2017.02.017. References [1] G. Meister, T. Tuschl, Mechanisms of gene silencing by double-stranded RNA, Nature 431 (2004) 343e349. [2] L. He, G.J. Hannon, MicroRNAs: small RNAs with a big role in gene regulation, Nat. Rev. Genet. 5 (2004) 522e531. [3] S. Vasudevan, Y. Tong, J.A. Steitz, Switching from repression to activation: microRNAs can up-regulate translation, Science 318 (2007) 1931e1934. [4] S.M. Johnson, H. Grosshans, J. Shingara, M. Byrom, R. Jarvis, A. Cheng, et al., RAS is regulated by the let-7 microRNA family, Cell 20 (2005) 635e647. [5] L. Ma, J. Teruya-Feldstein, R.A. Weinberg, Tumour invasion and metastasis initiated by microRNA-10b in breast cancer, Nature 449 (2007) 682e688. [6] M.F. Segura, D. Hanniford, S. Menendez, L. Reavie, X. Zou, S. Alvarez-Diaz, et al., Aberrant miR-182 expression promotes melanoma metastasis by repressing FOXO3 and microphthalmia-associated transcription factor, Proc. Natl. Acad. Sci. U. S. A. 106 (2009) 1814e1819. n, T. Oskarsson, D. Padua, Q. Wang, P.D. Bos, et al., [7] S.F. Tavazoie, C. Alarco Endogenous human microRNAs that suppress breast cancer metastasis, Nature 451 (2008) 147e152. [8] S. Volinia, G.A. Calin, C.-G. Liu, S. Ambs, A. Cimmino, F. Petrocca, et al., A microRNA expression signature of human solid tumors defines cancer gene targets, Proc. Natl. Acad. Sci. U. S. A. 103 (2006) 2257e2261. [9] S. Deng, G.A. Calin, C.M. Croce, G. Coukos, L. Zhang, Mechanisms of microRNA deregulation in human cancer, Cell Cycle 7 (17) (2008) 2643e2646. pez-Urrutia, V. García-Castillo, N. Jacobo-Herrera, [10] A. Pedroza-Torres, E. Lo L.A. Herrera, O. Peralta-Zaragoza, et al., MicroRNAs in cervical cancer: evidences for a miRNA profile deregulated by HPV and its impact on radioresistance, Molecules 19 (2014) 6263e6281. [11] E.-W. Hsing, S.-G. Shiah, J.-R. Hsiao, P.-C. Lyu, J.-Y. Chang, The role of miR-450a in oral squamous cell carcinoma carcinogenesis, FASEB J. 29 (1) (2015). Supplement 711.22. [12] G.A. Calin, C.M. Croce, MicroRNA signatures in human cancers, Nat. Rev. Cancer 6 (2006) 857e866. [13] S. Sethi, S. Ali, D. Kong, P.A. Philip, F.H. Sarkar, Clinical implication of MicroRNAs in molecular pathology, Clin. Lab. Med. 33 (2013) 773e786. [14] Y. Toiyama, M. Takahashi, K. Hur, T. Nagasaka, K. Tanaka, Y. Inoue, et al., Serum miR-21 as a diagnostic and prognostic biomarker in colorectal cancer,

27

J. Natl. Cancer Inst. 105 (12) (2013) 849e859. [15] S. Griffiths-Jones, R.J. Grocock, S. van Dongen, A. Bateman, A.J. Enright, miRBase: microRNA sequences, targets and gene nomenclature, Nucleic Acids Res. 34 (Database issue) (2006) D140eD144. [16] S. Griffiths-Jones, H.K. Saini, S. van Dongen, A.J. Enright, miRBase: tools for microRNA genomics, Nucleic Acids Res. 36 (Database issue) (2008) D154eD158. [17] C. Desjardin, A. Vaiman, X. Mata, R. Legendre, J. Laubier, S.P. Kennedy, et al., Next-generation sequencing identifies equine cartilage and subchondral bone miRNAs and suggests their involvement in osteochondrosis physiopathology, BMC Genomics 15 (2014) 798. [18] N. Terron-Canedo, W. Weir, L. Nicolson, C. Britton, L. Nasir, Differential expression of microRNAs in bovine papillomavirus type 1 transformed equine cells, Vet. Comp. Oncol. 4 (APR 2016), http://dx.doi.org/10.1111/vco.12216. Version of Record online. [19] W. Ragland, G.H. Keown, G.R. Spencer, Equine sarcoid, Equine Vet. J. 2 (1) (1970) 2e11. [20] E. Marti, S. Lazary, D.F. Antczak, H. Gerber, Report of the first international workshop on equine sarcoid, Equine Vet. J. 25 (1993) 397e407. [21] J.J. Bertone, Equine Geriatric Medicine and Surgery, Saunders Elsevier, 2006. [22] G. Chambers, V.A. Ellsmore, P.M. O'Brien, S.W.J. Reid, S. Love, M.S. Campo, et al., Association of bovine papillomavirus with the equine sarcoid, J. Gen. Virol. 84 (2003) 1055e1062. [23] D.C. Knottenbelt, S.E.R. Edwards, E.A. Daniel, Diagnosis and treatment of the equine sarcoid, Pract. 17 (1995) 123e129. [24] R.R. Pascoe, D.C. Knottenbelt, Neoplastic conditions, in: R.R. Pascoe, D.C. Knottenbelt (Eds.), Manual of Equine Dermatology, Saunders, London, 1999, pp. 244e252. [25] D.C. Knottenbelt, A suggested clinical classification for the equine sarcoid, Clin. Tech. Equine Pract. 4 (2005) 278e295. [26] L. Bogaert, M. van Poucke, C. De Baere, J. Dewulf, L. Peelman, R. Ducatelle, et al., Bovine papillomavirus load and mRNA expression, cell proliferation and p53 expression in four clinical types of equine sarcoid, J. Gen. Virol. 88 (2007) 2155e2161. [27] D.C. Knottenbelt, in: International Conference „Equine Oncology and Chosen Aspects of the Medical Law”. Book of Abstracts: 61e70, Wrocław, 2014, 7 VI. [28] C. Olson Jr., R.H. Cook, Cutaneous sarcoma-like lesions of the horse caused by the agent of bovine papilloma, Proc. Soc. Exp. Biol. Med. 77 (1951) 281e284. [29] W.L. Ragland, C.A. McLaughlin, G.R. Spencer, Attempts to relate bovine papilloma virus to the cause of equine sarcoid: horses, donkeys and calves inoculated with equine sarcoid extracts, Equine Vet. J. 2 (1970) 168e172. [30] M.B. Stocks, S. Moxon, D. Mapleson, H.C. Woolfenden, I. Mohorianu, L. Folkes, et al., The UEA sRNA workbench: a suite of tools for analysing and visualizing next generation sequencing microRNA and small RNA datasets, Bioinformatics 28 (2012) 2059e2061. [31] The RNAcentral Consortium, RNAcentral: an international database of ncRNA sequences, Nucleic Acids Res. 43 (Database issue) (2015) D123eD129. [32] M.I. Love, W. Huber, S. Anders, Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2, Genome Biol. 15 (2014) 550. [33] J. Li, R. Tibshirani, Finding consistent patterns: a nonparametric approach for identifying differential expression in RNA-Seq data, Stat. Methods Med. Res. 22 (5) (2013) 519e536. [34] I.S. Vlachos, N. Kostoulas, T. Vergoulis, G. Georgakilas, M. Reczko, M. Maragkakis, et al., DIANA miRPath v.2.0: investigating the combinatorial effect of microRNAs in pathways, Nucleic Acids Res. 40 (Web Server issue) (2012) W498eW504. [35] I. Kanellos, T. Vergoulis, D. Sacharidis, T. Dalamagas, A.G. Hatzigeorgiou, S. Sartzetakis, T. Sellis, MR-microt: a MapReduce-based MicroRNA Target Prediction Method, SSDBM, 2014, p. 47. [36] M. Reczko, M. Maragkakis, P. Alexiou, I. Grosse, A.G. Hatzigeorgiou, Functional microRNA targets in protein coding sequences, Bioinformatics 28 (6) (2012) 771e776. [37] M. Kanehisa, S.K.E.G.G. Goto, Kyoto encyclopedia of genes and genomes, Nucleic Acids Res. 28 (2000) 27e30. [38] C.B. Nielsen, N. Shomron, R. Sandberg, E. Hornstein, J. Kitzman, C.B. Burge, Determinants of targeting by endogenous and exogenous microRNAs and siRNAs, RNA 13 (2007) 1e18. [39] C.L. Andersen, J.L. Jensen, T.F. Ørntoft, Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets, Cancer Res. 64 (2004) 5245e5250. [40] M.W. Pfaffl, A new mathematical model for relative quantification in real-time RT-PCR, Nucleic Acids Res. 29 (2001) e45. [41] M. Nielsen, J.H. Hansen, J. Hedegaard, R.O. Nielsen, F. Panitz, C. Bendixen, et al., MicroRNA identity and abundance in porcine skeletal muscles determined by deep sequencing, Anim. Genet. 41 (2010) 159e168. [42] L. Guo, Y. Zhao, S. Yang, H. Zhang, F. Chen, An integrated analysis of miRNA, lncRNA, and mRNA expression profiles, Biomed. Res. Int. 2014 (2014) 345605. [43] A.H. Lund, miR-10 in development and cancer, Cell Death Differ. 17 (2010) 209e214. [44] A. Shinozaki, T. Sakatani, T. Ushiku, R. Hino, M. Isogai, S. Ishikawa, et al., Downregulation of MicroRNA-200 in EBV-associated gastric carcinoma, Cancer Res. 70 (2010) 4719e4727. [45] Y. Du, Y. Xu, L. Ding, H. Yao, H. Yu, T. Zhou, et al., Down-regulation of miR-141 in gastric cancer and its involvement in cell growth, J. Gastroenterol. 44

28

K. Pawlina et al. / Biochimie 137 (2017) 20e28

(2009) 556e561. [46] X. Hu, J.K. Schwarz, J.S. Lewis Jr., P.C. Huettner, J.S. Rader, J.O. Deasy, et al., A microRNA expression signature for cervical cancer prognosis, Cancer Res. 70 (2010) 1441e1448. [47] G. Muralidhar, M. Barbolina, The miR-200 family: versatile players in epithelial ovarian cancer, Int. J. Mol. Sci. 16 (2015) 16833e16847. [48] P.A. Gregory, A.G. Bert, E.L. Paterson, S.C. Barry, A. Tsykin, G. Farshid, et al., The miR-200 family and miR-205 regulate epithelial to mesenchymal transition by targeting ZEB1 and SIP1, Nat. Cell Biol. 10 (5) (2008) 593e601. [49] M. Korpal, E.S. Lee, G. Hu, Y. Kang, The miR-200 family inhibits epithelialmesenchymal transition (EMT) and promotes mesenchymal-epithelial transition (MET) by direct targeting of E-cadherin transcriptional repressors ZEB1 and ZEB2, J. Biol. Chem. 283 (22) (2008) 14910e14914. [50] T. Kim, A. Veronese, F. Pichiorri, T.J. Lee, Y.-J. Jeon, S. Volinia, et al., p53 regulates epithelialemesenchymal transition through microRNAs targeting ZEB1 and ZEB2, J. Exp. Med. 208 (2011) 875e883. [51] L. Rapp, Y. Liu, Y. Hong, E.J. Androphy, J.J. Chen, The bovine papillomavirus type 1 E6 oncoprotein sensitizes cells to tumor necrosis factor alpha-induced apoptosis, Oncogene 18 (1999) 607e615. € tzler, A. Bielak-Zmijewska, [52] L. Potocki, A. Lewinska, J. Klukowska-Ro W. Grabowska, I. Rzeszutek, et al., Sarcoid-derived fibroblasts: links between genomic instability, energy metabolism and senescence, Biochimie 97 (2014) 163e172. [53] A. Hall, The cytoskeleton and cancer, Cancer Metastasis Rev. 28 (2009) 5. [54] S.L. Fernandez-Valverde, R.J. Taft, J.S. Mattick, Dynamic isomiR regulation in Drosophila development, RNA 16 (2010) 1881e1888. [55] C.T. Neilsen, G.J. Goodall, C.P. Bracken, IsomiRsethe overlooked repertoire in the dynamic microRNAome, Trends Genet. 28 (2012) 544e549. [56] T.T. Bizuayehu, C.F. Lanes, T. Furmanek, B.O. Karlsen, J.M. Fernandes,

[57]

[58]

[59]

[60]

[61] [62]

[63]

[64]

[65]

S.D. Johansen, et al., Differential expression patterns of conserved miRNAs and isomiRs during Atlantic halibut development, BMC Genomics 13 (2012) 11. A. Wojcicka, M. Swierniak, O. Kornasiewicz, W. Gierlikowski, M. Maciag, M. Kolanowska, et al., Next generation sequencing reveals microRNA isoforms in liver cirrhosis and hepatocellular carcinoma, Int. J. Biochem. Cell Biol. 53 (2014) 208e217. P. Loher, E.R. Londin, I. Rigoutsos, IsomiR expression profiles in human lymphoblastoid cell lines exhibit population and gender dependencies, Oncotarget 5 (18) (2014) 8790e8802. K. Pawlina, A. Gurgul, M. Oczkowicz, M. Bugno-Poniewierska, The characteristics of the porcine (Sus scrofa) liver miRNAome with the use of next generation sequencing, J. Appl. Genet. 56 (2) (2015) 239e252. A.G. Telonis, P. Loher, Y. Jing, E. Londin, I. Rigoutsos, Beyond the one-locusone-miRNA paradigm: microRNA isoforms enable deeper insights into breast cancer heterogeneity, Nucleic Acids Res. 43 (19) (2015) 9158e9175. L. Guo, T. Liang, J. Yu, Q. Zou, A comprehensive analysis of miRNA/isomiR expression with gender difference, PLoS One 11 (5) (2016) e0154955. H. Seitz, J.S. Tushir, P.D. Zamore, A 5’-uridine amplifies miRNA/miRNA* asymmetry in Drosophila by promoting RNA-induced silencing complex formation, Silence 2 (2011) 4. H.Y. Hu, Z. Yan, Y. Xu, H. Hu, C. Menzel, Y.H. Zhou, et al., Sequence features associated with microRNA strand selection in humans and flies, BMC Genomics 10 (2009) 413. F. Frank, N. Sonenberg, B. Nagar, Structural basis for 5’-nucleotide basespecific recognition of guide RNA by human AGO2, Nature 465 (2010) 818e822. J. Starega-Roslan, P. Galka-Marciniak, W.J. Krzyzosiak, Nucleotide sequence of miRNA precursor contributes to cleavage site selection by Dicer, Nucleic Acids Res. 43 (22) (2015) 10939e10951.