Discovery and characterization of the feline miRNAome

The domestic cat is an important human companion animal that can also serve as a relevant model for ~250 genetic diseases, many metabolic and degenerative conditions, and forms of cancer that are analogous to human disorders. MicroRNAs (miRNAs) play a crucial role in many biological processes and their dysregulation has a significant impact on important cellular pathways and is linked to a variety of diseases. While many species already have a well-defined and characterized miRNAome, miRNAs have not been carefully studied in cats. As a result, there are no feline miRNAs present in the reference miRNA databases, diminishing the usefulness of medical research on spontaneous disease in cats for applicability to both feline and human disease. This study was undertaken to define and characterize the cat miRNAome in normal feline tissues. High-throughput sequencing was performed on 12 different normal cat tissues. 271 candidate feline miRNA precursors, encoding a total of 475 mature sequences, were identified, including several novel cat-specific miRNAs. Several analyses were performed to characterize the discovered miRNAs, including tissue distribution of the precursors and mature sequences, genomic distribution of miRNA genes and identification of clusters, and isomiR characterization. Many of the miRNAs were regulated in a tissue/organ-specific manner.

This study was undertaken to define the cat miRNAome and the expression and potential regulation of microRNAs (miRNAs) in normal feline tissues. Many other species have already been analyzed for miRNA expression (such as the dog), but miRNAs have not been carefully analyzed in cats. Previous work includes computational prediction of miRNAs within the cat genome 9, 10 and miRNA high-throughput sequencing performed on a feline kidney cell line 13 . However, no high-throughput sequencing and validation has been performed on normal cat tissues, thus far. As a result, there are no feline miRNAs present in miRBase, the reference microRNA database (www.mirbase.org) 14 . In this study, high-throughput sequencing was performed on 12 different normal cat tissues. 271 candidate feline miRNA precursors, encoding a total of 475 mature sequences were identified, including several novel cat-specific miRNAs. Many of these miRNAs are regulated in a tissue-specific manner.

Results
Generation of RNA libraries from cat tissues and processing of sequencing data. To define and characterize the cat miRNAome, we first generated twenty-seven small-RNA libraries from different tissue samples of 9 different cats (Supplementary Table S1.1). High-throughput sequencing was performed and the raw reads were pre-processed by the SOLiD software Lifescope. Processed high-quality sequences were analyzed by miRDeep2, a computational tool to map, analyze and score deep sequencing data for the identification of known and novel miRNAs 15 (See Materials and Methods section for extended details).
MiRDeep2 analysis returned 1182 putative miRNA precursors, but only 378 of them had a score above the chosen minimum threshold of 5 and these were further processed by our in-house pipeline (see Materials and Methods; scripts are available upon request). We first removed candidates that were poorly expressed or that matched potential repetitive sequences and other types of short RNAs. We thus obtained 273 candidate miRNA precursors, encoding a total of 475 mature sequences. We analyzed conservation, tissue distribution, genomic location, arm preference and sequence variation of these novel miRNAs. The results of this analysis are summarized in Supplementary Table S1.1-S1.5 and further described in the following sections. Figure 1 illustrates the computational pipeline employed for data analysis.

Conservation analysis.
We performed a preliminary BLAST analysis against the cat genome and the Rfam and Refseq-RNA databases, in order to filter out highly repetitive sequences and sequences matching other types of short RNA. This analysis removed two candidates, as they were found to be part of longer rRNA sequences. In addition, one of these two candidates had more than five significant hits against different chromosomal areas of the cat genome (E < 2e-08).
We then performed a BLAST search for the predicted miRNA precursors across miRBase (v. 21), to assess their conservation and identify potential orthologs. For 246 out of the 271 sequences, the analysis reported a significant match with sequences from at least one other species (EValue < 1e-5), with an average of 32 orthologous miRNAs per candidate. Human, cow and dog were the species with most matches. The predicted mature sequences were highly conserved, in most cases identical to their closest homolog, and were assigned to miRBase families according to their seed sequence (see Supplementary Table S1.2 columns E, F, K and AK and Supplementary Table S1.3).
The 25 unmatched precursors were considered putative feline specific miRNAs. They encoded a total of 33 mature sequences, 28 of which had a conserved seed and were assigned to the corresponding miRBase families. It should be noted that most of these cat specific miRNAs had rather low counts and, with a few exceptions, appeared to be tissue specific. We assigned them a temporary ID based on their chromosome.
Finally, we performed an in-silico analysis to investigate whether there was evidence in the cat genome of the conserved miRNAs that were missed by our analysis. For this task, we considered only the 268 miRNA families that were conserved in at least five species including human. These families contained 515 different miRNA genes, 214 of which were recovered by our deep sequencing analysis (41%) (see Supplementary Table S1.2, column F). Moreover, for 16 miRNAs missed by our analysis, we found at least one other miRNA from the same family (See Supplementary Table S12, column F). We then performed a BLAST search for the remaining 301 miRNA genes against the cat genome using human mature and precursor sequences as probes. We obtained a significant match for 164 distinct mature miRNAs encoded by 124 miRNA genes. However, only 33 out of the 124 miRNA genes exhibited conservation of their precursor sequences. Overall, our deep sequencing and in-silico analyses recovered 216 out of 268 conserved miRNA families (Supplementary Table S1.2, column BQ and S1.12, column C). A more thorough computational analysis would be necessary for a more accurate estimate. However, this goes beyond the purpose of our study.
Tissue distribution. Although most miRNAs were ubiquitously expressed, we observed significant differences in their distribution across the examined tissues. Figure S1 reports the number of different miRNA precursors and mature sequences across the analyzed samples. Unsupervised hierarchical clustering based on normalized expression of the most highly expressed mature miRNAs from each of the 271 precursors accurately distinguished between the different tissues (Fig. 2). We performed differential expression analysis of the mature miRNAs in order to identify signatures of tissue-specific miRNAs. We applied DESeq 2, a method for differential analysis of RNAseq count data, to perform pairwise comparisons of each tissue against all the others 16 . Uneven group sizes may result in lower power with the groups containing fewer samples. Nevertheless, our analysis identified groups of organ-enriched miRNAs in brain, liver, pancreas, lymph nodes, kidney, testis, spleen and lung (BH-adjusted p-value < 0.05). Table 1 summarizes all the significant tissue/organ-enriched and tissue/organ-specific miRNAs identified by our analysis (Details in Supplementary Table S1.2 and Supplementary  Tables S3). Arm selection preference. A miRNA precursor may encode one or two functional mature miRNAs, one from each arm of the hairpin (5p and 3p) and the most expressed product is commonly referred to as the Genomic Distribution and miRNA clusters. The miRNA precursors identified by our analysis were widely distributed throughout the genome with the exception of chromosome Y where no miRNAs were located. In humans, only two miRNAs have been found on chromosome Y. We discovered an average of 14.4 miRNAs per chromosome. Chromosomes B3 and B2 had the highest and the lowest number of miRNAs (if chromosome Y is excluded), 46 and 3, respectively. One miRNA was also assigned to the "unknown chromosome" chrUn_ JH409706, and four were assigned to the fragment chrC1_JH408690_random. (See Supplementary Table S1.4).
Our analysis also identified 31 different miRNA clusters. A cluster is a group of precursors with an inter-miRNA gene distance of less than 10 kb on the same genomic strand 19 . Clustered miRNAs are generally Computational pipeline of data analysis. The figure illustrates the four steps of the computational pipeline employed to analyze the RNAseq data. Pre-Processing: raw data were processed by the SOLiD software Lifescope in order to obtain good quality mappable reads. miRNA detection: this step was carried out by applying miRDeep2 to the mappable reads. Post-processing: the output of miRDeep2 was further analyzed by BLAST against different databases in order to assess conservation of the predicted miRNAs and remove sequences matching other kinds of small RNAs. Data Analysis: this step consisted of the application of a series of ad-hoc scripts for the extraction of descriptive statistics. The tool IPA was used to perform the functional enrichment analysis of potential targets of the identified miRNAs, which were predicted by the software miRiam.
transcribed as polycistronic primary transcripts (pri-miRNAs) and then processed into shorter pre-miRNAs to generate distinct mature miRNAs. It has been demonstrated that clustered miRNAs may regulate gene expression either individually or in combination in a coordinated manner and thus they may be functionally related. Figure 3 depicts the genomic distribution of the miRNA genes and highlights their organization in clusters. The largest cluster was identified on chromosome B3 and consisted of 29 distinct miRNA genes. We observed a large area of about 10 Mb on chromosome X encompassing 5 close clusters containing a total of 16 different miR-NAs. The size of the other clusters ranged from 2 to 5 miRNAs. Table 2 and Supplementary Table S1.5 provide a detailed report of the clusters, including conservation in human and dog.
Of the 31 clusters, 28 and 27 were completely or partially conserved in human and dog, respectively. Our analysis also revealed a cluster conserved in horse on chromosome X, consisting of miR-514 and miR-8908n. Among the non-conserved clusters, we identified a cluster with two cat-specific miRNAs on chromosome X (mir-chrX_38640 and mir-chrX_38642), a cluster containing miR-153 and a cat-specific miRNA on chromosome A2, and a cluster on chromosome E2 containing mir-295, mir-302d-1 and mir-371. miR-295 has not been characterized either in human or in dog, but orthologs were reported for mouse and rat, while miR-371 was conserved in many species, including human and dog, though clustered with other miRNAs. On the other hand, miR-302d-1 was a non-conserved precursor expressing a mature sequence similar to that of human miR-302d-3p.
IsomiRs. The term isomiR refers to variations in size and sequence from the canonical reference miRNA sequence annotated in miRBase. Initially considered to be sequencing artifacts, isomiRs were recently reported as functional variants with a specific biological role, probably generated by variation in processing by Drosha and/or Dicer 20 . The majority of miRNA genes encode mature isomiRs and it appears that in some cases human miRNA genes express isomiRs as the dominant transcript in specific cell types. IsomiRs differ from canonical miRNAs by the addition or the removal of one or more bases at the 5′ and/or 3′ end of the sequence 21 . The additional bases usually match the reference precursor sequence, thus these isomiRs are called templated. Non-templated and polymorphic isomiRs are rarer isoforms which harbor nucleotide changes from the precursors at either the 5′/3′ ends or internally. These variants are probably the result of post-transcriptional modifications, i.e. editing. A-to-I editing is the most common form of RNA editing 22 .
Our analysis classified the detected mature miRNA sequences into 5′/3′ templated and non-templated isomiRs, polymorphic isomiRs and canonical forms. Figure 4 shows the distribution of the different types of variants overall and in the different tissues. Overall, the most common products observed were the canonical form (46.3%), and the 3′ templated isoforms (33.4%), accounting for roughly 80% of the total number of reads. This was expected and in line with previous findings in other animal models, as 3′ isomiRs are much more likely to preserve the function of the canonical form. Target recognition, indeed, is mostly mediated by the seed sequence at the 5′ end of the miRNA, thus length modifications at the 3′ end are likely to not have any functional impact. The third most frequently observed type of variant was the polymorphic isomiR (11.2%), followed by 3′ non-templated and 5′ templated isomiRs (6.0% and 2.4%, respectively). The latter is particularly relevant, as 5′ modifications of the miRNA sequence may significantly affect target recognition and, consequently, be functionally important, as reported by recent work 23 . Non-templated 5′ isomiRs and simultaneous 5′/3′ templated and non-templated variants all together accounted for less than 1% of the total reads. Details on isomiRs are given in Supplementary Table S1.7.

Validation of candidate miRNAs by Stem-Loop qRT-PCR.
We confirmed the expression of 19 different microRNAs from our high-throughput sequencing data, including nine novel, cat-specific miRNAs. We chose representative miRNAs that were enriched in different tissues/organs to confirm that our high throughput sequencing data was quantitative (see Fig. 5). We also confirmed by qRT-PCR that fca-miR-chrD4_30107-3p and fca-miR-chrE3_34323-5p are enriched in the brain (data not shown). Each of these miRNAs was significantly enriched in specific tissues/organs, except for fca-miR-26a, fca-miR-151-3p, fca-miR-361-5p and fca-miR-chrC2_23051-3p, which were constitutively expressed in all tissues/organs examined. Heat maps of both the high-throughput sequencing and the stem-loop qRT-PCR data were virtually superimposable ( Fig. 5 and Supplementary Table S1.10). Furthermore, a Pearson's correlation coefficient analysis between the deep-sequencing and TaqMan data revealed that most miRNAs had a very high correlation between the two assays (see Supplementary Table S1.11). Only miR-26a-5p and miR-23051-3p were not significant, probably due to the fact that they were constitutively expressed. However, the correlation was still high (0.46 and 0.55).

Comparison with previous studies. A recent study by Sun et al. performed high throughput sequencing
on a feline kidney cell line (F81) before and after infection with mink enteritis virus (MEV) 13 . A careful comparison between our study and the Sun study revealed 186 common miRNAs (Supplementary Table S1.8). In addition, there were 264 unique miRNAs in our study and, surprisingly, 78 unique miRNAs in the Sun study.
In another study by Tamazian et al., a total of 370 unique feline miRNA genes were detected and mapped based upon homology to miRNA sequences from 36 species described in the miRBase database 10 . We have done a careful comparison between the lists of miRNAs generated in Tamazian's and in our study (see Supplementary  Table S1.9). We identified 213 miRNAs that had the same name and location between the two studies, confirming our assignment of these 213 miRNAs. In addition, there were 58 and 157 miRNAs that were unique to our study and Tamazian et al., respectively.
Functional analysis of cat-specific miRNAs. Our analysis identified 25 novel, cat-specific miRNA precursors encoding 33 mature miRNAs. With the exception of seven miRNAs that were expressed in at least two different tissues, most novel miRNAs appeared to be tissue specific. In order to investigate the potential roles for the novel cat-specific miRNA candidates, we performed functional enrichment analysis of their predicted targets. In particular, we focused on 21 tissue-specific miRNAs expressed in brain, kidney, testis, spleen and ovary (see
Similarly, the results for the three kidney-specific miRNAs, including fca-chrE3_33626-5p, which was validated by RT-PCR, showed significant enrichment for pathways that are relevant for kidney function, such as AMPK signaling, cell cycle checkpoint control, Estrogen receptor signaling, Renin-Angiotensin signaling and Renal cell carcinoma signaling (BH-adjusted P < 0.05).
The results for the spleen-specific miRNAs showed significant enrichment for pathways relevant for spleen function, such as EGF Signaling, FLT3 Signaling in Hematopoietic Progenitor Cells, IGF-1 Signaling and CD40 Signaling (BH-adjusted P < 0.05).

Discussion
The domestic cat is not only a wonderful companion animal, but it also serves as a translational model for many human genetic, metabolic, and degenerative diseases and forms of cancer. In order to improve feline medicine and translate findings to humans, it is necessary that we have knowledge of the expression and regulation of miR-NAs in normal feline tissues, since miRNAs play a crucial role in many biological processes and pathways and are implicated in most if not all diseases. In this study, we performed a genome-wide high-throughput sequencing analysis to define and characterize the feline miRNAome, which consists of 271 candidate feline miRNA precursors, encoding a total of 475 mature sequences, some of which appear to be tissue/organ-enriched and tissue/ organ-specific.
Several of these miRNAs have also been shown to be tissue-enriched in other species, thus supporting the results of our analysis. For example, miR-1, miR-206 and miR-133a are highly enriched in cardiac and skeletal muscle and are collectively known as myomiRs 24 . They play a fundamental role in the regulation of muscle cell differentiation, development and maintenance 25 . Our analysis confirmed that they are also highly enriched in cat muscle tissue (lip and tongue samples). Other oral tissue-enriched miRNAs that we identified included miR-24, which has been reported as up-regulated in oral squamous cell carcinoma 26 , and miR-184, which was shown to be associated with anti-apoptotic and proliferative processes in tongue carcinoma 27 .
Our analysis identified 88 miRNAs specifically enriched in brain tissue. Among them, miR-219, miR-124, miR-153, miR-128, miR-132 and miR-139 are known to be brain enriched or brain-specific in other species and to be involved in several brain-specific functions 28 . For example, miR-219 has been reported to inhibit the proliferation, anchorage-independent growth and migration of glioma cells and to promote oligodendrocyte differentiation and myelination, while miR-124 has been shown to promote neurogenesis, inhibit proliferation of glioblastoma multiform cells and induce differentiation of brain tumor stem cells. Another important brain-enriched miRNA that was also found highly expressed in brain tissue by our analysis was miR-128. This miRNA has been reported to be involved in synaptogenesis, reduce neuroblastoma cell motility and invasiveness, and regulate apoptosis and inhibit proliferation and self-renewal in glioma.
Other tissue-specific miRNAs identified by our analysis that were also reported to be tissue-specific in other species include: miR-122, a liver-specific miRNA that functions as a tumor-suppressor gene in hepatocellular carcinoma 29 ; miR-216, a pancreas-enriched miRNA, which has been reported as a marker for acute phase pancreatic injury and whose down-regulation is thought to be crucial in the development of pancreatic cancer 30 ; miR-150, which our analysis reported as enriched in lymph nodes, is known to be expressed in mature B and T cells and, in particular, to regulate differentiation and the cytolytic effector function in CD8+ T cells 31 . Finally, our analysis identified a cluster of testis-enriched miRNAs located on chromosome X, including miR-506, miR-507, miR-508 and miR-514, that was previously reported as preferentially expressed in testis 32 .
Several studies have shown that arm selection preference may be tissue-specific 33 . We found that, based on the total number of reads per tissue, the choice of the dominant arm was consistent across all tissues for 252 out of 271 precursors (93%), while 19 of them exhibited arm switches. In some cases, arm switching appeared to be tissue specific, such as miR-378, whose dominant form was 5p in the ovary and 3p in all other tissues, and miR-493 and miR-140, whose dominant product was 5p in all tissues except skin, where the 3p form was expressed at a higher level. This is consistent with the fact that miR-140-3p was reported to be the dominant form in melanocytes and melanoma cells by previous studies 34 . In other cases, we observed less specific arm switching events across tissues. For example, the dominant product for let-7i was 5p in testis, ovary, lymph node, spleen, lung and oral tissue and 3p in brain, pancreas, skin, liver and kidney.
Arm switching can also be species specific 35 . For example, the dominant miR-10 sequences in fly (Drosophila melanogaster) and beetle (Tribolium castaneum) are processed from opposite arms 18 . Arm usage appears to be regulated through sequences in the primary miRNA sequence, outside the mature miRNA duplex, and the targets of miRNAs encoded from opposite arms may differ significantly. This strongly suggests that changes in arm preference throughout nature may have relevant functional consequences 36 .
Cluster analysis showed that, of the 31 clusters found in cat, 28 and 27 were completely or partially conserved in human and dog, respectively. For example, the miR-221/222 cluster on chromosome X was conserved in both human and dog also on chromosome X. These miRNAs, which also constitute a family, are known to be involved in cancer in humans 37 . They act either as oncogenes or tumor suppressors, depending on the type of tumor, and have been reported to have a role in drug resistance 38 . Recent studies, indeed, have revealed that targeting miR-221/222 may enhance sensitivity of cancer cell lines to drugs 39 . The miR-221/222 cluster is also considered a key player in vascular biology, as it contributes to vascular remodeling and plays a prominent role in atherosclerosis and in metabolic diseases, being involved in the regulation of insulin resistance 40 . Another important conserved cluster was the miR-17/92 cluster on chromosome A1. It is one of the most studied miRNA clusters and it is known to be involved in normal development and homeostasis, as well as in the pathogenesis of cancer, immune, cardiovascular and neurodegenerative diseases 41,42 . Our analysis also revealed a cluster conserved in horse on chromosome X, consisting of miR-514 and miR-8908n. The human ortholog of miR-514 (miR-514a) was reported to be involved in cancer by initiating melanocyte transformation and promoting melanoma growth 43 , while miR-8908n was characterized in horse only and its function remains unknown.
A careful comparison between our study and a study on a feline kidney cell line 13 revealed 186 common miRNAs. In addition, there were 264 unique miRNAs in our study and, surprisingly, 78 unique miRNAs in the Sun study. Our study had a broader scope and included 12 different tissues, thus the detection of a much larger number of miRNAs in our samples compared to just one cell line was reasonable to expect. However, it was surprising that Sun et al. identified 78 miRNAs that were not detected in our samples, especially since our study included three highly correlated kidney samples. Thus, these 78 miRNAs are probably induced during kidney cell immortalization/transformation and/or tissue culturing. Finally, there were 16 miRNAs that exhibited arm switches between the two studies. Five of these switches were very dramatic, going from > 90% of one arm in one study to > 90% of the opposite arm in the other study (fca-miR-1307, fca-miR-144, fca-miR-197, fca-miR-199 and fca-miR-33a). The most likely explanation for the differences observed between the two studies is that several miRNA genes are either turned off or on when cells are taken out of their primary environment and grown in tissue culture.
In another study, a total of 370 unique feline miRNA genes were detected and mapped based upon homology to miRNA sequences from 36 species described in the miRBase database 10 . We identified 213 miRNAs that were identical between the two studies, as well as 58 and 157 miRNAs that were unique to our study and Tamazian et al., respectively. The 157 miRNAs that we didn't observe in our study probably include miRNAs that are not expressed in cat and/or are expressed in cat tissues/organs that we did not examine.
Our analysis also identified 25 novel, cat-specific miRNA precursors encoding 33 mature miRNAs, some of which were expressed in a tissue-specific manner and may be involved in the regulation of important biological processes and pathways. Functional analysis of their predicted targets confirmed significant enrichment for several pathways that are specific or relevant in the enriched tissues (Tables S2.1-S2.5), supporting our findings. However, further investigations will be necessary to confirm such roles and functions. Despite the remarkable advances made in recent years, computational prediction of miRNA targets still represents a challenge, as target recognition by miRNAs is a very complex and dynamic mechanism, still only partially understood 44 . Moreover, the UTR-ome of the cat has not yet been fully characterized, thus our analysis was limited to only ~8,300 UTR sequences, which is considerably less than the 20,285 putative genes present in cat 9 . Hence, the missing UTRs could affect the analysis. Nevertheless, this analysis provided useful information about the potential role of the candidate cat-specific miRNAs and support for their tissue specificity.
It has been argued that species-specific and lineage-specific microRNAs constitute a large proportion of miR-NAs, and many of them have been recently identified through deep sequencing technologies in a variety of animals 45 . Although conserved miRNAs have been shown to be significantly more associated with disease than non-conserved miRNAs 46 , recent evidence suggests functional roles for species-specific miRNAs in a variety of biological processes, including apoptosis 47 , inflammation 48 and nervous system development 49 . Silkworm-specific miRNAs were reported to have a functional preference for regulating genes involved in life-cycle-associated traits 50 . A more recent study has shown that a set of cattle-specific miRNAs may regulate the expression of target genes involved in insulin signaling through both conserved and cattle-specific binding sites, thus delineating gene expression divergence between ruminant and non-ruminant species 51 .
In conclusion, this study filled the gap in our knowledge regarding feline miRNAs by characterizing the feline miRNAome through deep sequencing of multiple cat tissues and organs. Our analysis identified 271 candidate feline miRNA precursors encoding a total of 475 mature sequences, including several novel cat-specific miRNAs. Further analyses were performed to characterize the discovered miRNAs in terms of genomic distribution, tissue expression, arm preference, isomiRs and function. Our comprehensive miRNA dataset provides the scientific community with a novel and useful resource for basic biology and veterinary clinical research, which may also inform and translate research findings for the benefit of humans.

Materials and Methods
Tissues were obtained from control research cats and client-owned cats (with consent from the owners), collected at the time of euthanasia. The use of these research cats, as well as the client consent form, were approved by the Institutional Animal Care and Use Committee (IACUC) at The Ohio State University. All methods were performed in accordance with the IACUC guidelines and regulations.

Quantitative Stem-Loop
Bioinformatics and Statistical Analysis. Raw reads obtained from the Applied Biosystems SOLiD 5500 system were pre-processed by the SOLiD software Lifescope, which filtered out low quality reads and performed mapping to the cat reference genome (felCat5, downloaded from UCSC: https://genome.ucsc.edu/cgi-bin/ hgGateway?db=felCat5). High quality mapped reads were then processed by miRDeep2, a computational tool that identifies candidate miRNA precursors and their mature sequences from small RNA sequencing data. miRD-eep2 makes use of other tools, such as the short read aligner Bowtie and the RNA secondary structure prediction tool RNAfold from the ViennaRNA package 2.0 53, 54 . Short reads were aligned to the reference genome, then all candidates whose structure and read signature were inconsistent with Drosha/Dicer processing were filtered out. Potential hairpin precursors were assigned a score according to a Bayesian probabilistic model of miRNA biogenesis.
No cat miRNAs were included in miRBase 21, thus all the detected miRNAs were considered novel. Candidate miRNAs detected by miRDeep2 were filtered by their scores. Then we applied filters based on score, expression and conservation by using our in-house developed ad-hoc scripts (Python and R) and a local version of BLAST. Similarly to what has been described in previous literature, we set a score cut-off threshold corresponding to the lowest score that yielded a signal-to-noise ratio higher than 10:1 (21.5:1) 55, 56 . Such value was 5, thus all the candidates with score below 5 were discarded. We kept all the candidate miRNAs expressed at detectable levels in at least one tissue, i.e. read count ≥10, this being a reasonable and commonly used count threshold 57 .
BLAST was used to filter out highly repetitive sequences and sequences matching other types of short RNA, and to evaluate conservation at the precursor and the mature sequence level.
Tissue distribution was evaluated by pairwise average-linkage unsupervised hierarchical clustering of the most highly expressed mature sequences from each precursor using Spearman's rank correlation as a distance measure. Pairwise differential expression analysis was performed between different tissues by DESeq 2 in order to identify significant signatures of tissue-enriched miRNAs. Ad-hoc scripts were used to evaluate arm preference and detect arm switching across the samples, to evaluate the genomic distribution of the identified precursors and identify miRNA clusters, and to perform isomiRs analysis. Chromosome ideograms were created using the NCBI Genome Decoration Tool (http://www.ncbi.nlm.nih.gov/genome/tools/gdp).
We performed an in-silico analysis to investigate whether there was evidence in the cat genome of the conserved miRNAs that were missed by our analysis. We then performed a BLAST search for the remaining 301 miRNA genes against the cat genome using human mature and precursor sequences as probes. We discarded precursors exhibiting a match of less than 70% of their nucleotide bases and mature sequences without a complete seed match and EValue > 0.1.
The scripts used for the analysis are available upon request.