(C) PLOS One
This story was originally published by PLOS One and is unaltered.
 . . . . . . . . . .



Charting the equine miRNA landscape: An integrated pipeline and browser for annotating, quantifying, and visualizing expression [1]

['Jonah N. Cullen', 'Department Of Veterinary Population Medicine', 'College Of Veterinary Medicine', 'University Of Minnesota', 'St. Paul', 'Minnesota', 'United States Of America', 'Jakub Cieslak', 'Department Of Genetics', 'Animal Breeding']

Date: 2025-09 

Abstract MicroRNAs (miRNAs) are essential regulators of gene expression, yet few comprehensive databases exist for miRNA expression in non-model species, limiting our ability to characterize their roles in gene regulation, development, and disease. Similarly, isomiRs - length and sequence isoforms of canonical miRNAs with potentially altered regulatory targets and functions - have received even less attention in non-model species, including the horse, leaving a critical gap in our understanding of their biological significance. To address these challenges, we developed an open-source, containerized pipeline for identifying and quantifying miRNAs and isomiRs (FARmiR: Framework for Analysis and Refinement of miRNAs), and an associated interactive browser (AIMEE: Animal IsomiR and MiRNA Expression Explorer). AIMEE was developed to make miRNA expression data more accessible and user-friendly, a feature often lacking from other expression atlases. These tools were developed using equine data but can be readily extended to other species. Using these tools, we aggregated 461 small RNA-seq datasets, spanning 61 distinct tissues, integrating data from public repositories, an American Quarter Horse cohort, and the Functional Annotation of ANimal Genome (FAANG) consortium Thoroughbred samples, predicting 5,781 miRNAs and isomiRs. This work represents the largest systematically curated atlas of equine miRNA expression to date, providing a valuable resource that will enhance our understanding of miRNA and isomiR functions in tissue-specific regulation and ultimately improve biomarker discovery, functional genomics, and precision veterinary medicine.

Author summary MicroRNAs (miRNAs) play a critical role in gene regulation, influencing many biological processes including development, metabolism, and disease. Despite their importance, miRNA research in non-model species like the horse has remained limited. Moreover, what research has been done has been more narrowly focused on canonical miRNA expression in specific tissues under various disease states. To address this gap, we developed a pipeline and interactive browser for reproducibly quantifying and visualizing miRNA expression. Using these tools, we aggregated and analyzed 461 small RNA sequencing datasets from 61 equine tissue types, generating the most comprehensive equine miRNA expression atlas to date. Our findings highlight the complexity of miRNA and isomiR expression across tissues and delivers a user-friendly resource for researchers studying gene regulation. Our work enhances the ability to investigate miRNA-regulated biological processes in the horse and provides a framework for similar studies in other species.

Citation: Cullen JN, Cieslak J, Petersen JL, Bellone RR, Finno CJ, Kalbfleisch TS, et al. (2025) Charting the equine miRNA landscape: An integrated pipeline and browser for annotating, quantifying, and visualizing expression. PLoS Genet 21(9): e1011835. https://doi.org/10.1371/journal.pgen.1011835 Editor: Androniki Psifidi, Royal Veterinary College, UNITED KINGDOM OF GREAT BRITAIN AND NORTHERN IRELAND Received: February 17, 2025; Accepted: August 6, 2025; Published: September 5, 2025 Copyright: © 2025 Cullen et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability: All presented data are either already publicly available through NCBI’s Sequence Read Archive (SRA), the FAANG data portal, or EMBL-EBI’s European Nucleotide Archive. No unpublished results from other individuals are included. Primary: PRJNA1224819 (https://dataview.ncbi.nlm.nih.gov/object/PRJNA1224819?reviewer=o1n5qkg09j03k37uc8783480h3) FAANG: https://data.faang.org/dataset?sortTerm=accession&sortDirection=desc PRJEB14485: https://www.ebi.ac.uk/ena/browser/view/PRJEB14485 PRJEB20494: https://www.ebi.ac.uk/ena/browser/view/PRJEB20494 PRJEB27174: https://www.ebi.ac.uk/ena/browser/view/PRJEB27174 PRJNA573500: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA573500 PRJNA593513: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA593513 PRJNA637872: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA637872 PRJNA681404: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA681404 PRJNA705083: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA705083 PRJNA726388: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA726388 PRJNA306540: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA306540 PRJNA352827: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA352827 PRJNA169102: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA169102 PRJNA393269: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA393269 PRJNA553103: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA553103. Funding: This project was supported by the United States Department of Agriculture (USDA) National Institute of Food and Agriculture (NIFA) Grant 016-10133 (M.E.M, J.R.M.), University of Minnesota College of Veterinary Medicine NIFA Multistate competitive funds 2017-2020, 2020-2022 (M.E.M., J.R.M.), USDA NIFA Agriculture and Food Research Initiative (AFRI) Predoctoral Fellowship 2020-67034-31739 (J.N.C.), Grayson-Jockey Club Research Foundation (C.J.F., J.L.P, R.R.B.), USDA NIFA Animal Breeding and Functional Annotation of Genomes (A1201) Grant 2019-67015-29340 (J.L.P, C.J.F., R.B., T.K., E.G.) as well as NRSP-8 Species Coordinator Funds from the USDA NIFA (C.J.F., J.L.P, R.R.B.), and the the UC Davis Center for Equine Health (C.J.F., J.L.P, R.R.B.). SDA salary was supported by National Institutes of Health’s National Center for Advancing Translational Sciences, grants K12TR002492 and UL1TR002494. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health’s National Center for Advancing Translational Sciences. J.N.C was funded by USDA NIFA Fellowships, 2020-67034-31739 and 2024-67012-41761. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: The authors have declared that no competing interests exist.

Introduction MicroRNAs (miRNAs) are a class of endogenous small (~22 nucleotides) noncoding RNAs that primarily regulate gene expression through binding messenger RNA (mRNA), triggering mRNA degradation or the inhibition of translation [1]. Over 60% of human protein-coding genes contain at least one miRNA binding site, highlighting the importance of miRNA-mRNA pairing and miRNA-based regulation [2]. MicroRNAs play critical roles in numerous biological processes, including development [3], stem-cell maintenance [4], metabolism [5], apoptosis [6], and cell-cell communication [7]. Given their large regulatory footprint, miRNA dysregulation has been implicated in many diseases [8–12], and miRNA profiling has demonstrated both ubiquitous and tissue-specific expression across species [13–16]. High-throughput sequencing (HTS) for profiling of human small RNA expression has identified miRNA isoforms (isomiRs) that differ from the reference miRNA (Ref-miR) [17]. Originally believed to be sequencing or alignment errors, isomiRs are now thought to be the result of imprecise cleavage or post-transcriptional modification [18–21], resulting in differences in sequence composition and/or length from the reference form [22–24]. IsomiRs can be categorized into five classes based on compositional differences from their Ref-miR counterparts: reference or canonical (identical to the miRNA database sequence [25,26]), 5’ isomiR (length changes to the 5’ end), 3’ isomiR (length changes to the 3’ end), polymorphic (sequence modifications with identical length), and mixed (length and sequence changes) [27,28]. The region of the miRNA critical for mRNA target recognition is referred to as the “seed sequence”. IsomiRs with edited or shifted seed regions (e.g., 5’-isomiRs or polymorphisms within residues 2–8) are of particular interest [29,30], as they may have altered target sets that impact similar or different biological pathways as compared to their Ref-miR counterpart [31–36]. Similar to Ref-miRs, isomiRs can be highly expressed [28] with tissue-dependent expression patterns [37,38] and have been associated with various disease states, highlighting their potential as biomarkers [39–41]. Despite the availability of many isomiR-capable processing tools [42], “gold standards” do not exist for isomiR profiling. As a result, isomiRs are often excluded from miRNA studies and databases, introducing profiling biases and missing biological significance [43]. While several miRNA databases exist, only a few include isomiR expression [38,44], and, with a few exceptions, these are limited to humans and mice [45]. Characterization of the tissue-specific miRNA transcriptome in companion animals and livestock has lagged behind that of humans and model organisms for several reasons. First, miRNA profiling is often based on the available Ref-miR sequences from databases such as miRBase (v22) [25], and the more sparse annotation of miRNAs in these species may have a significant impact on expression profiling. For example, the current version of miRBase contains 2,656 mature miRNA sequences in humans and 1,978 in mice, compared to only 690 in horses. Of the miRBase horse miRNAs, over half (359) are based on a single in silico study [46], and the remainder (331) are from small RNA-sequencing (RNA-seq) of testes from one horse [47]. Moreover, while earlier versions of MirGeneDB (2.1) [26] did not include horse, the latest release (3.0) [48] provides a database of manually curated, high-confidence miRNAs for EquCab3. While this represents a major improvement in annotation quality over previous resources, it is more conservative in scope, containing 438 equine mature miRNA sequences from 417 precursors (69 of which are not in miRBase). As the goal of this study was to catalog all potentially expressed miRNAs, we included the broader set of miRBase sequences. Second, few large-scale miRNA expression profiles of normal equine tissues exist. Most of the equine miRNA efforts have focused on profiling select tissues from horses with or without disease [49–53], during exercise [54–56], or at different developmental stages [50]. MiRNA expression profiling across healthy tissues has been previously conducted with only limited tissue sets [57]. The largest study targeted nine tissues, including gluteus medius muscle tissues (the largest muscle in horses, involved in hip extension) from normal horses and horses with polysaccharide storage myopathy type I (PSSM1) [58]. While these studies identified tissue-specific reference and putative novel miRNAs, isomiRs were not thoroughly characterized. To the best of our knowledge, only two equine miRNA analyses have included isomiRs [59,60]. Third, miRNA data can be difficult to access for animal researchers intending to explore and utilize previous findings, and analyses of miRNA data requires bioinformatics proficiency. Finally, the majority of available computational tools and pipelines for characterizing isomiRs were developed for processing human data, many of which require additional informatic manipulations to configure for non-model organisms [42]. To address these limitations, we developed the Framework for Analysis and Refinement of miRNAs (FARmiR) pipeline and Animal IsomiR and MiRNA Expression Explorer (AIMEE) browser. We used these tools to catalog the most comprehensive collection of miRNA expression at the isoform resolution in normal equine tissues, resulting in the identification of 5,463 isomiRs. These resources allow for the reproducible processing and analysis of high-throughput small RNA-seq datasets, and the export of raw and normalized data for additional downstream analyses. The aggregated expression data, together with the pipeline and browser, will be a significant resource for researchers investigating the composition and functional roles of miRNAs in species that lack expression databases like those available for humans and mice. Finally, our pipeline-browser framework will enable researchers to process their own RNA-seq samples in a consistent and reproducible manner, facilitating comparisons across different disease states, phenotypes of interest, and developmental processes with our comprehensive catalog.

Discussion From 462 samples across 61 tissue types, AIMEE includes the largest collection of uniformly processed equine miRNA expression data to date. This expression atlas and explorer contains nearly 6,000 putative and canonical miRNA sequences, ~ 95% of which are isomiRs. The inclusion of isomiRs is of particular value as the importance of these molecules has rapidly grown in the last decade [28]. Consistent with previous studies, the equine atlas contains far fewer 5’-isomiRs compared to 3’-isomiRs [32]. Due to the well-known batch effects inherent in RNA-seq data [71,72], robust rank aggregation was used to circumvent the challenges associated with batch correction methods and integrate samples from heterogeneous sources. Also, in the context of minimizing unwanted technical variation, the tissue specificity index (TSI) was calculated for each unique miRNA that was sufficiently expressed (i.e., mean RPM > 20 in at least one organ) independently for the included primary and FAANG samples. From the FAANG samples, for example, a 3’-isomiR of eca-miR-7 and a 5’-isomiR of eca-miR-215 were exclusively expressed in the pituitary and cecum, respectively. The miR-7 miRNA family has been associated with regulating gene expression in the pituitary of mice [73], rats [74], and cattle [75]. Similarly, the miR-215 family has been repeatedly identified as a critical tumor suppressor in gastric and colorectal cancer [76,77], as well as associated with ulcerative colitis [78]. The scope of this miRNA expression atlas has implications for both basic research and future clinical applications. The tissue-specific miRNA and isomiR profiles generated as part of this work advances our understanding of gene regulatory mechanisms across equine tissues. In particular, the isomiR cataloging of various muscle types, cartilage, and blood, creates opportunities for biomarker discovery, improving precision in the diagnosis and monitoring of various equine conditions. Data from previous studies investigating the role of miRNAs in equine osteochondritis and osteoarthritis [49,79] can be analyzed along with the uniformly processed tissue data from this study, possibly contributing to the identification of additional biomarkers for early osteoarthritis detection. Furthermore, understanding the distribution of specific miRNAs across tissues is essential for developing targeted therapeutic interventions, particularly for metabolic and inflammatory conditions where miRNA-mediated regulation plays a significant role [80]. The characterization of isomiRs in individual tissues highlights new possibilities for personalized medicine approaches. For example, the isomiR catalog we developed for miRNAs like eca-miR-7 (174 seed variants) and eca-miR-140-3p (125 seed variants) underscores the potential for individual variation in isomiR expression. This type of fine-tuned control of the regulatory machinery could have significant functional consequences that may warrant careful consideration in therapeutic strategies. This equine miRNA expression resource does have some limitations. First, publicly available samples that met the inclusion criteria were based entirely on the descriptions reported by the authors in the original published main text, supplemental materials, or data repositories. Consequently, we cannot be as confident in the health assessments of these public samples compared to the primary or FAANG sets which included extensive ante and postmortem phenotyping. Second, we only identified 252 of the 690 miRBase-defined canonical miRNAs [25]. This is in contrast to previous studies that identified greater proportions, for example, in the testes [50,81]. This could be for a couple of reasons. Identifying miRNAs in earlier studies, especially in studies that were not using tools designed for miRNA identification but instead repurposing general aligners like BWA [82], can over-inflate miRNA counts. Even when using tools designed for miRNA profiling, there are significant differences in sensitivity and accuracy and poorly overlapping miRNA sets between tools [83]. Additionally, misannotated or spurious miRNAs in the reference set may hinder accurate identification and quantification. Many profiling tools utilize miRBase [25] or miRCarta [84] to define the miRNA space of a given species. For the horse, this may be suboptimal, as the 690 miRBase reference is based on two studies, with over half (359) being based on in silico prediction [46], and the remainder (331) from small RNA of testes of one horse of unknown breed and health status [47]. Both of these studies relied on the previous horse assembly [85]. In contrast, MirGeneDB v3.0 provides a curated set of 438 mature miRNAs (387 unique sequences) based on the current EquCab3 assembly [86]. We observed 259 of the 387 unique sequences (66.9%) exactly matched the mature sequences identified in this atlas, with an additional 2,150 exhibiting full-length containment. These findings underscore the potential biological complexity introduced by isomiRs and the challenges of maintaining comprehensive miRNA annotations. Importantly, FARmiR was not designed for this task but instead to catalog expression using the known miRNA space and predict additional miRNA loci beyond that space. As a result, it is possible that the parental miRNA precursors (and locations) have not been filtered against the current repeat element annotations, leading to hairpin-shaped precursors that are not true miRNAs [63]. Third, we did not include any in silico miRNA target prediction nor enrichment analyses based on predicted targets. This decision was made for two reasons: 1) it has been previously reported that computational target prediction (without experimental validation) in non-model organisms is not reliable [87]. 2) A follow-up to this study will couple prediction and validation to predict the RNA targets more accurately. Finally, this atlas profiled miRNA expression at the tissue-level, using bulk small RNA-seq. Characterization at the cellular level has demonstrated cell-specific expression of both human Ref-miRs and isomiRs [38,88,89]. Although still in its infancy and technically challenging, recent advances in HTS technologies have demonstrated the feasibility and potential for single-cell miRNA-seq protocols [90,91]. To date, there have been no rigorous characterizations of the equine miRNA landscape using single-cell RNA-seq methodology, thus representing a natural extension of this current work. Using FARmiR and AIMEE represents a unified approach to miRNA data processing, analysis, and visualization. As an open-source, containerized pipeline for miRNA identification and quantification, FARmiR aligns with the principles of sustainable data analysis, promoting transparency, adaptability to other species, and in silico reproducibility [92,93]. AIMEE provides researchers a user-friendly and accessible option for exploring miRNA tissue profiles without the substantial effort and cost associated with preparing small RNA-seq data. AIMEE also includes a filtering option to consider only the mature sequences from this atlas that overlap with the curated set provided by MirGeneDB v3.0 [48]. The primary, FAANG, and public data, uniformly processed with FARmiR and accessible via AIMEE, contribute valuable information about the expression of equine miRNAs and, ultimately, the tissue-specific regulatory functions contributing to equine health and performance.

Materials and methods . Ethics statement. Tissues were collected from horses donated for use in this study at the time of euthanasia. All procedures were approved by the Institutional Animal Care and Use Committee at the University of Minnesota (1712-35369A). All methods were performed in accordance with the IACUC guidelines and regulations. Three small RNA-seq data sources were utilized for this atlas, one primary (i.e., the 12 horses enrolled as part of this study) and two public sources. The latter two sources consisted of four Thoroughbreds from the publicly available Functional Annotation of Animal Genomes (FAANG) consortium (referred to as FAANG), and systematic identification of compatible datasets from public databases. Primary data source Horses. A total of 12 healthy horses (11 American Quarter Horses and 1 American Paint Horse) donated to the University of Minnesota provided the primary source of data for this study. Physical examinations were performed by a boarded large animal internal medicine specialist to ensure horses were clinically normal. Horses were housed on a dry lot as a group and were fed free choice hay and water. Isolation and twice daily temperature checks were performed for two weeks to safeguard against communicable diseases. The primary source consisted of two cohorts with the first in the fall of 2013 (3 mares, one gelding, and one stallion 5–15 years of age) and the second in the fall of 2020 (7 mares 12–14 years of age). Tissue collection and RNA isolation. Horses were sedated with 0.05 mg/kg intravenous (IV) injection of xylazine, an alpha-2 adrenergic agonist, and humanely euthanized with ≥ 100 mg/kg injection of pentobarbital sodium (Fatal-Plus, Vortech Pharmaceutical) following the American Association of Equine Practitioners Euthanasia Guidelines (https://aaep.org/resource/euthanasia-guidelines/). Immediately following euthanasia, all tissues were collected in triplicate and flash frozen in liquid nitrogen or placed in RNAlater Stabilization Solution (ThermoFisher Scientific). Both flash-frozen and RNAlater aliquots were stored at -80°C until further processing. A total of 31 tissues were collected from the first cohort and up to 70 tissues from the second cohort (S1 Table). All tissues were collected within 45 min of euthanasia. The following tissues (n = 31) were selected for further processing: articular cartilage, lamina, subchondral bone, synovium, multiple muscles (gluteal, longissimus, and sacrocaudalis dorsalis), duodenum, jejunum, cecum, liver, spleen, thyroid, kidney (cortex and medulla), pituitary, cerebellum, hypothalamus, lung, heart (left and right atrium and ventricle), multiple fat depots (nuchal, omental, retroperitoneal, shoulder, tailhead, and visceral), ovary, and uterus. Only tissues with at least two biological replicates were included for library preparation and sequencing. Tissues were selected to capture a range of potential expression diversity across tissue types, as well as to represent those actively being investigated by our laboratory and others. The remaining tissues are biobanked in the University of Minnesota Equine Genetics and Genomics Laboratory tissue repository. RNA was isolated using tissue-type specific protocols to maximize RNA quality and yield. Due to wide diversity of tissue types and composition, no single extraction protocol was appropriate for all sampled tissues. Tissue type specific protocols were developed and refined based on composition and assessment of RNA quality via NanoDrop 8000 spectrophotometer (Thermo Scientific, cat. no ND-8000-GL, Waltham, MA, USA) and the Agilent TapeStation system (Agilent Technologies, Santa Clara, CA, USA). See Supplemental Materials for complete descriptions of isolation protocols. Tissues isolated and sequenced as part of this study are referred to as primary throughout the remainder of this manuscript. Small RNA library preparation and sequencing. Sample QC, library preparations, and sequencing reactions were conducted at GENEWIZ, LLC/Azenta US, Inc (South Plainfield, NJ, USA). The RNA samples were quantified using a Qubit 2.0 Fluorometer (ThermoFisher Scientific, Waltham, MA, USA), and RNA integrity (RIN) was checked using the Agilent TapeStation (Agilent Technologies, Santa Clara, CA, USA). The Small RNA sequencing libraries were prepared using the NEBNext Small RNA Library Prep Kit for Illumina using the manufacturer’s instructions (New England Biolabs, Ipswich, MA, USA). Briefly, Illumina 3’ and 5’ adapters were added to RNA molecules with a 5’-phosphate and a 3’-hydroxyl group sequentially. A reverse transcription reaction was used to create single-stranded cDNA. The cDNA was then PCR amplified using a common primer and a primer containing an index sequence. The amplified cDNA construct was purified by polyacrylamide gel electrophoresis, and the correct band (~145–160 bp) was excised from the gel and eluted with water. The eluted cDNA was concentrated by EtOH precipitation generating the sequencing libraries. The sequencing libraries were multiplexed and clustered onto a flowcell. After clustering, the flowcell was loaded onto the Illumina HiSeq 4000 or NovaSeq 6000 instrument according to the manufacturer’s instructions. The samples were sequenced using a 2 × 150 bp paired-end (PE) configuration. Image analysis and base calling were conducted by the HiSeq Control Software (HCS). Raw sequence data (.bcl files) were converted into FASTQ files and de-multiplexed using Illumina bcl2fastq 2.20 software. One mismatch was allowed for index sequence identification. Equine Functional Annotation of ANimal Genomes (FAANG) data Ethics statement. All procedures were approved by the Institutional Animal Care and Use Committee at the University of California - Davis (protocols #19037 and 21033). All methods were performed in accordance with the IACUC guidelines and regulations. Horses. As part of the equine FAANG project, small RNA-seq libraries were constructed and sequenced using tissue samples collected from two healthy Thoroughbred mares (aged 4 and 5 years) [61] and two healthy Thoroughbred stallions (aged 3 and 4 years) [62]. Tissue collection and RNA isolation. Horses were sedated and humanely euthanized prior to sample collection and preservation as described [61,62]. The following tissues (n = 39) were selected for RNA isolation and miRNA sequencing: adrenal cortex, amygdala, articular cartilage (fetlock and stifle), bone marrow, brain, cecum, cerebellum vermis, deep digital flexor tendon, dorsal root ganglia, epididymis (caput, cauda, and corpus), fat depots (abdominal and loin), heart (left and right ventricle, right atrium), hippocampus, kidney (cortex and medulla), lamina, liver, lung, lymph node, muscle types (gluteal and longissimus dorsi), ovary, parietal cortex, peripheral blood mononuclear cells (PBMC), pituitary, sesamoid bone, two skin sections (neck and unpigmented), superficial digital flexor tendon, T8 thoracic vertebrae, temporal cortex, testis, and third metacarpal bone (S2 Table). Small RNA library preparation and sequencing. Sample QC, library preparations, and sequencing reactions were conducted at the University of Minnesota Genomics Center (UMGC) (n = 40) or Admera Health (n = 46). At UMGC, total RNA isolates were quantified using a fluorimetric RiboGreen assay on a BioTek Synergy 2 (BioTek Instruments, Winooski, VT, USA). Total RNA integrity was assessed using capillary electrophoresis (Agilent BioAnalyzer 2100, Agilent Technologies, Santa Clara, CA, USA) to generate RNA Integrity Numbers (RIN). While UMGC recommends total RNA yields of at least 500ng and RINs > 8, TruSeq Small RNA libraries were prepared for samples (n = 36) with a range of RINs (5.4 – 9.5). Libraries for four bone-derived samples (two sesamoid and two third metacarpals) with the lowest RINs (3.7 – 4.0) were prepared using the Clontech smRNA-Seq kit (Takara Bio USA, San Jose, CA, USA) for Illumina following the manufacturer’s instructions, which is more accommodating for low-input RNA. The sequencing libraries were multiplexed and clustered onto a flowcell. After clustering, the flowcell was loaded onto the Illumina HiSeq 2500 and sequenced using a 2 × 50 bp PE configuration. Raw sequence data (.bcl files) were de-multiplexed and converted into FASTQ files using CASAVA v1.8.2 software. At Admera Health, isolated RNA sample quality was assessed by BioAnalyzer RNA Eukaryotic Nano Assay (Agilent Technologies, Santa Clara, CA, USA) and quantified by Invitrogen Broad Range RNA Qubit Assay (Invitrogen, Thermo Fisher Scientific, Waltham, MA, USA). Library construction was performed with TruSeq Small RNA kit (Illumina, San Diego, CA, USA). Final libraries quantity was assessed by Qubit 2.0 Fluorometer (ThermoFisher Scientific, Waltham, MA, USA) and quality was assessed by TapeStation HSD1000 ScreenTape (Agilent Technologies, Santa Clara, CA, USA). Mean final library size was about 200 bp with an insert size of about 50 bp. Illumina 6-nt single-indices were used. Equimolar pooling of libraries was performed based on QC values and sequenced on an Illumina NovaSeq 6000 with a read length configuration of 150 PE for 20M PE reads per sample (10M in each direction) and data was trimmed to 1x50 sequencing configuration. Tissue samples sequenced across different runs or lanes were combined prior to processing. In the primary dataset, both fragment ends were sequenced for all libraries. However, only the first read of the FAANG data was available and the included public samples were sequenced from single-end (SE) libraries. Therefore, only the first read from the primary dataset was included to mitigate library configuration as a source of unwanted technical variation or batch effect, an unavoidable and inevitable challenge with aggregating multi-center RNA-seq data [71,94]. Public data sources Database search and validation. The National Center for Biotechnology Information (NCBI) BioProject and PubMed databases were queried using “(equine OR horse) AND (miRNA OR microRNA)” on April 4, 2023, and May 3, 2023, respectively, to identify potentially relevant studies and datasets. The search was limited to studies and datasets describing Illumina-based small RNA-seq of equine tissue samples with publicly available sequencing data. Studies that profiled non-blood derived fluids/cell types (e.g., semen, milk, synovial fluid), extracellular vesicles, or other culture-expanded cell types (e.g., skeletal muscle satellite cells, chondrocytes) were excluded to focus on native tissue-derived expression. These sample types were omitted due to the potential for increased technical variability associated with the additional, type-specific handling and preparation steps, the evaluation of which was beyond the scope of this study. In contrast, blood derived sample types were included due to their widespread usage and established relevance in miRNA expression studies despite variability in processing. From each study, samples were considered eligible if described as healthy without experimental treatment. For example, data from the pre-challenge samples in a challenge study [95] were included, while data from horses described as “non-performers” in an exercise or endurance trial were excluded [56]. Similarly, if horses were not described as “non-performers,” both pre- and post-race samples were considered eligible [55]. Samples were excluded if health status was not described or could not be inferred from available metadata. Since it has been suggested that sex and breed-specific miRNA profiles may exist in the horse [58] and other species [96,97], samples lacking sex or breed information were also omitted. Mixed breeds were only included if at least one parental breed was described (i.e., horses labeled “cross” were excluded, whereas horses labeled “Thoroughbred-cross” were included). Samples described as male were classified as male-unknown unless expressly indicated as stallion or gelding. In the case of females, if age was reported, a distinction was made between mare (>3 years) and filly. Similarly, males under 3 years were labelled colt. All sample data passing inclusion criteria were downloaded from the NCBI Sequence Read Archive (SRA) with parallel-fastq-dump v0.6.6 (https://github.com/rvalieris/parallel-fastq-dump) (a fastq-dump wrapper from the SRA toolkit v2.11.0 [98]) or directly from the European Nucleotide Archive (ENA). Study details and sample data for included public samples are reported in S3 Table. Small RNA processing and analysis We developed the Framework for Analysis and Refinement of miRNAs (FARmiR; https://github.com/jonahcullen/FARmiR) pipeline to catalog miRNA expression in a reproducible, transparent, and adaptable manner, adhering to the principles of sustainable data analysis [92]. FARmiR is an open-source, Snakemake-based [92] containerized pipeline [93] developed for quality assessment, pre-processing, and miRNA profiling (Fig 11). Although designed for use with any animal species, we illustrate its application and inputs using the current equine reference genome EquCab3 [86]. PPT PowerPoint slide

PNG larger image

TIFF original image Download: Fig 11. Overview of the FARmiR pipeline. Created with BioRender.com. https://doi.org/10.1371/journal.pgen.1011835.g011 Read quality assessment and pre-processing. Raw sequence reads from all three sources were first assessed for quality without any filtering using fastp v1.3.1 [99]. Adapters and low-quality reads (> 40% of bases with phred quality scores < 15 or at most 5 missing base calls) were removed. For SE samples, the initial approach employs adapt_find to identify and remove adapters [100]. If adapters were not detected, DNApi is used for adapter prediction followed by cutadapt [101,102]. Although not included here, adapter detection and removal with fastp is used for samples with PE layouts [99]. The adapter prediction functionality of fastp was not included for SE samples due to observed inconsistencies in adapter prediction across datasets. We instead implemented this sequential approach to improve robustness without requiring prior knowledge of the adapter sequence, which is often not readily available in the publication. If included, PE reads will be joined (no mismatches allowed and a minimum overlap of six bases required) [103]. By convention, miRNA profiling typically only considers SE-sequenced libraries (or the first read of PE layouts). We designed FARmiR to handle PE layouts as there is recent evidence that PE sequencing of small RNA-seq libraries may improve isomiR identification and accuracy [104]. Sequence error correction, a critically important consideration for isomiR discovery and tissue-specific quantification, was conducted with miREC [105], followed by length filtering (minimum length of 17) [99]. To remove likely non-miRNA molecules, cleaned reads were first aligned against the EquCab3 [86] non-coding reference (excluding miRNA sequences) from Ensembl (release 103). Unmapped reads were then aligned against non-coding RNA (excluding miRNA sequences) from the Rfam v14.8 database [106]. Both alignments use bowtie v1.3.1 [107] with zero mismatches permitted (“-n 0”), only one reportable alignment allowed (“-m 1”), and without alignment against the reverse-complement (“--norc”). Reportable alignments were limited to one as there are existing miRNA sequences that have significant overlap with other non-miRNA regions, such as piRNAs [108]. Reads unmapped against Rfam with lengths greater than 25 bases were removed (referred to as post-filtered reads) [99]. Pre- and post-filtered reads were aligned against EquCab3 using bowtie v.1.3.1 [107]. Default parameters were used except allowing for no mismatches (“-n 0”) within the first 8 bases (“-l 8”), up to 6 multiple alignments (“-m 6”) and reporting only the best stratum alignment (“-a --best --strata”). A maximum of 6 multiple alignments were permitted as there are at most 6 multimaps between the miRBase mature and hairpin sequences. We note however there exist up to 11 multimaps between human mature and hairpin sequences. Mapped reads and the Ensembl (release 103) EquCab3 annotation were input to featureCounts v2.0.1 [109], allowing for multi-mapped reads (“-M”) with per-exon counting (“-t exon”). MultiQC was used to aggregate pre-processing steps [110]. Sample filtering. Samples with log 2 post-filtering counts less than three times the unbiased median absolute deviation were excluded [111]. Following this, samples with reads consisting of less than 10% miRNA biotypes (as estimated by featureCounts) were also discarded. Finally, only tissues with at least two independent samples (i.e., regardless of age, sex, or breed) were retained for quantification. This approach was taken to improve overall data quality and reduce the influence of batch effects. Updated miRNA annotation. The current version of the equine miRNA annotation from miRBase is based on EquCab2 [85], and to the best of our knowledge, an EquCab3 annotation is not publicly available. To address this, we used a combination of Liftoff [112] and LiftOver [113] to map hairpins and mature sequence loci from EquCab2 to EquCab3. Of the 715 non-redundant miRNA hairpin genes (4 hairpins are encoded at more than one locus: eca-mir-703 at 5 loci, eca-mir-9128–2 at 3, and both eca-mir-1842 and eca-mir-9150 at 2 loci each), 7 hairpins were unmapped during conversion to EquCab3. All 7 unmapped hairpin loci (eca-mir-8921–1, eca-mir-8921–2, eca-mir-8923–4, eca-mir-8923–5, eca-mir-8934–1, eca-mir-8935–6, and eca-mir-8946) were predicted to be on unplaced EquCab2 contigs and are not included in either the Ensembl 103 or most recent 109 annotations. However, of the 690 mature miRNA sequences in EquCab2 (59 miRNAs are encoded by more than one hairpin), only one (eca-miR-8946) was not identified in EquCab3. Discovery and quantification of miRNAs Novel miRNA loci. Samples that passed filtering were used for novel miRNA prediction with mirPRo v.1.1.4 [114]. MirPRo uses Novoalign v3.02.12 [115] and a prediction algorithm similar to miRDeep2 [116] with modifications to reduce false positives. Novel prediction was performed using default parameters with mature and hairpin miRNA sequences of human, mouse, dog, cow, and pig from miRBase v22 as related species [25]. miRBase v21 family information was used as miRBase v22 does not contain miRNA family information. There is no apparent difference between miRBase v21 and v22 families for equine miRNAs. BLAST+ v2.13.0 [117] was used to query predicted hairpins against equine RNA sequences from RNAcentral [118] to remove other possible hairpin-forming RNAs (e.g., tRNA and rRNA). Genomic loci of predicted hairpins with a maximum overlap of 17 (based on a minimum mature sequence length of 17 as above) were merged with BEDtools [119]. Predicted hairpin sequences were then queried against the EquCab3 reference genome, and sequences with at most 15 alignments were retained. Predicted hairpin counts were then normalized using the post-filter counts from above to reads per million (RPM) and filtered to include only hairpins with RPMs ≥ 20 in at least 2 samples. In the cases of overlapping and merged loci, the hairpin with the highest total counts across all samples was retained. If total counts were equal between two or more hairpins, one was randomly selected to define the genomic locus. Expressed candidate loci that overlap (maximum of 17) with the updated miRNA annotation hairpin loci were removed. Retained expressed candidate loci were then used to select the associated mature representative. Novel miRNAs were named following a modified strategy previously established for sheep miRNA profiling [120]. Briefly, retained candidate loci were queried (BLAST+ v2.13.0 with default settings and -task blastn) against the miRBase precursor sequences for cattle, sheep, pig, mouse, and human. Alignments with a q-value (minimum false discovery rate (FDR)) below 0.01 and sufficient query coverage (> 80%) were retained and the subject(s) with the highest percent identity was selected. For queries with valid alignments to more than one species, the evolutionarily closest species was chosen for naming. Note the original mirPRo naming convention (i.e., eca-novel-mir-NNNN) was preserved for downstream analyses and within AIMEE to avoid potential confusion. Orthologs of predicted miRNAs are noted where applicable and available in the supplemental data (S4 Table). The candidate hairpins and mature miRNAs were converted to BED format and merged with the updated miRNA annotation to be utilized for quantification. Ref-miRs and isomiRs. Quantification was conducted with isoMiRmap [63]. IsoMiRmap requires species-specific “mapping bundles.” Currently, the developers only provide a human miRNA mapping bundle. To address this, we generated a species mapping bundle prior to quantification, starting with the putative miRNA annotation (see Loher et al. (2021) [63] for details and https://github.com/jonahcullen/FARmiR for implementation). The post-filtered and sized reads from each sample, along with the mapping bundle, were then quantified by isoMiRmap. This procedure returns (among other outputs) a Ref-miR/isomiR miRGFF3 annotation [121], a proposed standard from the miRNA TRanscriptomic Open Project (miRTOP), which utilizes a previously described naming convention [122] and nucleotide-based unique identifier per isomiR [123]. RPM (normalized by the final post-filtered and sized read counts) and count data were extracted from these annotations and assembled into an expression matrix. Ref-miR/isomiRs were retained only where at least 2 samples have RPMs greater than 20. Following this, samples consisting of less than 10% non-zero RPMs were excluded, and ultimately, tissues were dropped if less than 2 independent samples remained. While this filtering improves interpretability, it may impact detection of lowly expressed miRNAs, particularity in sample types with high compositional variability across RNA biotypes. Various metrics were then calculated for each Ref-miR (e.g., numbers of hairpins, isomiRs, and seed-shifted isomiRs). Shiny app and browser. To promote accessibility and maximize insights derived from the collected data, we developed a Shiny web application built with the golem framework [124]. This app, Animal IsomiR and MiRNA Expression Explorer (AIMEE; https://github.com/jonahcullen/AIMEE), was designed to provide user-friendly accessibility and point-and-click functionality for data exploration and visualization, simple analyses, and exporting quantification data for subsequent investigations. AIMEE utilizes expression matrices generated by FARmiR and associated sample metadata as input for multiple visualization and analysis modules. These include simple quality control, read loss through FARmiR processing, isolating Ref-miRs/isomiRs of interest, and upset analysis with ComplexUpset [125]. Assessing inter- or intra-tissue Ref-miR/isomiR rank aggregation and tissue-specificity via the organ enrichment index (OEI) is also possible with AIMEE. Robust rank aggregation. Unbiased integration of small RNA-seq data is critical when attempting to make inferences amidst unknown or unquantified sources of technical variation. One such approach to reduce these biases is to represent each sample as an independent prioritized miRNA gene list (i.e., miRNAs ranked by expression at a user-defined threshold). These gene lists can then be aggregated using order statistics, with p-values assigned to each miRNA denoting how its rank in each sample list surpasses what would be expected by chance. AIMEE implements the R package RobustRankAggreg (RRA) [67,68], a probabilistic aggregation model that is robust to the variations associated with high-throughput data from disparate origins, to assign significance probabilities to each miRNA from user-selected tissue or sample sets. As RRA operates on within-sample rankings as opposed to raw expression values, it helps reduce – but does not eliminate – the impact of technical variation across samples. A more traditional batch effect correction (e.g., ComBat from the sva R package [126]) was not feasible given the diverse and unbalanced nature of the included samples. Organ-enrichment index. The OEI is a measurement of miRNA specificity [14,127], analogous to the mRNA tissue specificity index (TSI) ‘tau’ [128]. OEI values range from 0 to 1, with scores close to 0 indicating a miRNA that is expressed at similar levels across many or most organs (i.e., ‘housekeepers’), and scores close to 1 corresponding to miRNAs expressed in only one specific organ (i.e., organ-enriched). To calculate OEI, tissue types originating from the same organ were grouped together. Specifically, heart chambers were collapsed into “heart”, kidney regions into “kidney”, and skin samples into “skin”. Brain regions (parietal and temporal cortex, hippocampus, hypothalamus, amygdala, cerebellum and cerebellum vermis, and dorsal root ganglia) were combined into a single “brain” grouping. All other tissues were retained as separate groupings to preserve potential site-specific expression and enrichment. The mean RPM was calculated for each miRNA per organ or combined organ grouping. Only miRNAs with a mean RPM > 20 in at least one organ or group were retained. Then, the OEI was calculated for each miRNA j as where x j,i represents miRNA j in organ i normalized by the maximum RPM across all organs, and N is the total number of organs considered. Importantly, due to the likely presence of batch effects, AIMEE contains OEI values calculated separately using the primary or FAANG collection.

Acknowledgments We would like to thank the following individuals for providing assistance preparing tissue samples and RNA extractions: Em Adam, Samantha Beeson, Kendall Blanchard, Lauren Hughes, Elaine Norton, Kelsey Springer, and Matthew Sturos. We would also like to thank Ernest F Bailey, Vinzenz Gerber, Vidhya Jagannathan, Tosso Leeb, Leslie A Lyons, and Cecilia M Penedo.

[END]
---
[1] Url: https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1011835

Published and (C) by PLOS One
Content appears here under this condition or license: Creative Commons - Attribution BY 4.0.

via Magical.Fish Gopher News Feeds:
gopher://magical.fish/1/feeds/news/plosone/