Sequencing, alignment, and annotation
As an initial quality check, we mapped all adapter-trimmed reads against the LepWed1.0 genome assembly, with no gaps or mismatches allowed. Approximately 70% of reads obtained from tissues aligned perfectly to the Weddell seal genome and only 43-55% of reads from plasma samples were perfectly aligned (Additional file 1: Fig. S1). These results illuminated a large portion of the read data that was unmapped to the Weddell seal genome (Additional file 1: Fig. S2), ranging from 291,923 unmapped reads in a cortex sample to 2,076,932 in a plasma sample. Plasma, in particular, was over-represented in the unmapped data, with 5 of 6 plasma samples containing >1,000,000 unmapped reads. As several unaligned, putative miRNAs were highly expressed, we further identified reads that did not map to the Weddell seal genome, but were perfectly aligned to a miRNA hairpin sequence annotated in miRBase [26]. For example, the sequence TGAGATGAAGCACTGTAGCT was abundant in our dataset, with 532,097-1,014,236 reads in each sample. This sequence did not map to the Weddell seal genome but was identified as miR-143-3p. Thus, assembly incompleteness appears to reduce the percentage of small RNA reads that could be mapped to the Weddell seal genome, with the greatest impact on plasma samples.
In order to annotate high confidence miRNA loci, we manually curated the set of predictions provided by miRCat2 [27] and miRDeep2 [28]. This led to a final set of 559 loci (union of high confidence predictions for each tool; Additional file 2: Table S1; Additional file 1: Fig. S3). Among these, 329 corresponded to a miRBase annotated hairpin sequence, while the remaining 230 represent novel miRNAs.
We examined variation in the relative abundance of miRNA-5p and miRNA-3p across all samples and did not identify any clear examples of arm switching (i.e. changes of the most abundant miRNA strand). We also used our annotation to identify adenylated/uridylated miRNA isoforms and to investigate their representation in each sample. Plasma displayed the highest rates of 3’ end modification, with up to 12% of all reads corresponding to a modified isoform (Additional file 1: Fig. S1).
Differential miRNA expression among tissues
416 miRNAs were differentially expressed (DE) in a single sample type compared to all others. Among these, 74 were DE in two tissues (Additional file 1: Fig. S4). 50% of tissue-specific DE occurred in the brain, with 31% of significant results in plasma. A smaller set of loci were uniquely elevated or depressed in heart (8%) and muscle (11%; Table 1). These two contractile tissues were also the most similar to each other when expression data was viewed as a heatmap (Fig. 1), although all 4 tissue types were segregated by hierarchical clustering. Principal Component Analysis (PCA) clearly separated samples by tissue along PC1 and PC2 (Fig. 2A). As with the heatmap, PCA points to greater similarity between the contractile tissues, heart and muscle, relative to plasma and brain. Moreover, brain is clearly separated from the other tissues along the PC1 axis, with limited inter-individual variability, especially for PC2. Random forests analyses identified 6 miRNA inputs that were able to separate the four tissue types from each other with zero error, including one novel, unmapped Weddell seal miRNA novel-4-3p, whose expression is significantly upregulated in plasma and downregulated in muscle (Table 2; Fig 3A). Four of these classifiers were upregulated in heart (miR-490-3p, miR-499-5p, miR-30e-5p and miR-30d-5p) and four were downregulated in plasma (miR-95-3p, miR-499-5p, miR-30a-5p and miR-30d-5p; Additional file 3: Table S6).
To shed light on the biological role of miRNAs of interest, we mapped mRNA targets of DE miRNAs to significantly enriched gene pathways. MiRNA downregulation within a tissue could infer higher local expression of target mRNAs. Therefore, the activity of pathways associated with upregulated miRNAs are potentially lowered, while pathways associated with downregulated miRNAs and their targets are potentially enhanced.
MiR-499-5p was the most DE locus in the dataset, a Random Forests classifier, and has the highest expression in heart (Additional file 3: Table S6). mRNA targets of miR-499-5p, predicted to have lower expression in heart than other seal tissues include Jade1, a pro-apoptotic factor. miR-490-3p was also highly expressed in the heart, was a Random Forests classifier, and was among the five most significant DE changes in the dataset (Additional file 3: Table S6). miR-490-3p targets ion channels and transporters, including members of the KCN and SLC gene families (Kcng3, Kcnj15, Kcne2, Kcne4, Kcnmb1, Kcnip3, Kcnk12, Slc1A3, Slc9a1, Slc9a2), associated with pathway enrichments in potassium ion transmembrane transport (GO:0071804, 0071805). High relative expressions of both miRNAs are consistent with prior miRNA biomarker studies in heart [29, 30], and do not differ between pups and adults. Although KEGG disease pathways were excluded from the presented data, it is noteworthy that pathway analysis identified enrichment of three cardiac disease pathways (hsa05410: Hypertrophic cardiomyopathy, hsa05414: Dilated cardiomyopathy, hsa05412: Arrhythmogenic right ventricular cardiomyopathy) associated with mRNA targets of miRNAs that were significantly upregulated in hearts of Weddell seals.
The brain had more DE miRNAs than any other tissue, impacting the largest number of mRNA targets. In general, targets of DE miRNAs in the seal brain were enriched for neuronal and developmental processes (Additional file 3: Table S6). The largest pairwise fold change occurred in miR-488-3p (brain > plasma), which targets steroid hormone receptors including Oxtr, Ghrhr, and Pgrmc2. High relative expression in brain may downregulate the steroid hormone response pathway (GO:0048545).
MiR-296-3p was among 6 downregulated miRNAs in seal swimming muscle. Enhanced pathways associated with targets of this miRNA include those expected for skeletal muscle, related to ryanodine receptor redox state regulation, calcium handling and the sarcoplasmic reticulum (GO:0060314, 0033017, 0016529, 1901019, 0050848), as well as response to extracellular acidic pH (Rab11fip5, Rab11b, Impact, Asic1, Asic2, GO:0010447). We also identified elevated miR-206-3p in muscle, which [31] has been identified as a MyomiR, with expression restricted to this tissue. In contrast to heart, miR-490-3p was lowest in muscle, pointing to an enrichment of target ion channels and potassium transport. Both strands of miR-10 are upregulated in muscle. Although the 5p strand is most abundant (average 5p read counts: 2,233,312 in pups and 1,634,627 in adults; average 3p read counts: 413.6 in pups, 212 in adults), 3p has large tissue-specific fold changes (23-62X) compared to other tissues, downregulating predicted targets related to the regulation of blood circulation (GO:1903522), including Bves, Casq2, and Nos1.
Of the Random Forest classifiers DE in plasma (Table 2), miR-95-3p had the largest sample-specific downregulation (20-53X lower than heart, muscle, and brain; Additional file 3: Tables S3, S4, S6). miR-95-3p targets LDL receptor related protein (Lrp1) as well as the beta-subunit of guanylyl cyclase (Gucy1b1). miR-339-3p, also downregulated in plasma compared to other tissues, targets several key genes in hypoxia sensing (GO:0036293, GO:0070482, GO:0001666) including Epas1, Vhl, Hif1an, and Cygb, but does not change with age (Additional file 1: Fig. S5). Several miRNAs elevated in Weddell seal plasma versus other tissues (miR-18a-5p, miR-221-3p, and miR-34-5p) target pathways regulating blood vessel and vascular remodeling (e.g. GO:0001974), with miR-34-5p specifically targeting Epas1 and angiotensin (Agt).
Differential miRNA expression with development
80 miRNAs were DE between adults and pups across all tissues (Table 1). 188 miRNAs were DE in tissue-specific developmental comparisons, 27 of which were up or downregulated in more than one comparison (Additional file 1: Fig. S7). Heatmap visualization also separates pups from adults within each tissue (Fig. 1). An exception is the miRNA profile of adult skeletal muscle, which clusters more closely with heart samples than with muscle from pups. PC3 and PC4 provide the best developmental clustering (Fig. 2B) but explain only 13.7% of variation. Supervised clustering by Random Forests (Fig.3B) identifies 2 miRNA inputs (miR-29a-3p, miR-542-5p) that separate the samples by developmental stage for all tissues combined (Table 3; Fig. 3B). MiR-29a-3p is significantly upregulated in heart and muscle of adults, whereas miR-542-5p is downregulated in adult brain and muscle (Additional file 3: Table S10). Brain and skeletal muscle miRNAs had large developmental differences among the four sample types investigated (Table 1). Only 13% of significant adult-pup differences occurred in the heart and there were no developmental differences in plasma miRNAs (Table 1).
Predicted targets of miR-29a-3p that would be reduced from a terrestrial pup to a diving adult include components of the mitochondrial electron transport system (Cox4i1, Cox10, Ndufs7, Ndor1, Sdhaf2), genes related to iron regulation (Tfrc, Ireb2) and negative regulators of hypoxia signaling (Hif3a, Vhl). Conversely, higher miR-542-5p in pups is predicted to enhance the abundance of mRNAs that are likely relevant in the adult phenotype. Genes of interest include the mitochondrial citrate transporter Slc25a1, as well as genes known to regulate hematopoiesis (Cd44), anaerobic glycolysis (Ldha), and obesity/lipolysis (Plin1).
High number of lineage specific miRNA orthogroups
Recent studies have highlighted the high rate of novel miRNA gains in mammals [32]. We identified 874 loci clusters in a dataset that included Weddell seal, cow, dog, horse, pig, rabbit, mouse, and human. These were considered miRNA orthogroups and were evaluated to infer the evolutionary patterns of gain and loss across the phylogenetic tree (Fig. 4). High net gain rates in both the dog and the Weddell Seal lineages suggest dynamic miRNA evolution. For the seal, we also observe a high number of lineage specific losses. This may reflect purifying selection on young, selectively neutral miRNAs, however this result might be biased by differences in assembly completeness between the CanFam3.1 dog reference [33] and the LepWed1.0 genome.
Novel Weddell seal miRNAs are more likely to be tissue-specific, as we observed a significantly higher proportion of tissue-specific miRNA families in the novel set compared to the miRBase set (z test, cortex: p=0.0001; heart: p<10-4; muscle: p<10-4; plasma: p<10-4; Fig. 5). This aligns with previous findings, which suggest that young miRNAs are initially expressed in a single or few tissues, then become more broadly expressed later in their evolutionary history [34].
Predicted mRNA targets of novel Weddell seal miRNAs have a functional overlap with targets of annotated miRNAs in this dataset, as evidenced by similarities in pathway enrichments for tissue-specific up- and down-regulated miRNA targets, regardless of whether miRNAs were previously annotated or novel in Weddell seals. For example, VEGF signaling was targeted by both novel and annotated miRNAs in the seal heart, highlighting the importance of vascular development in cardiac tissue. We specifically investigated pathways identified from novel Weddell seal miRNAs, to evaluate the capability for species-specific gene regulation and manifestation of phenotype.
We compared pathways targeted by DE novel versus miRBase-annotated miRNAs, identifying unique pathway targets of novel miRNAs in each tissue. 37 pathways were enriched (p<0.05, see Methods) for the targets of miRNAs that were highest in the brain and 31 pathways were enriched for targets of muscle-elevated miRNAs, which are primarily signaling pathways (Additional file 3: Tables S10, S11). Conversely, only two pathways in each tissue were associated with tissue-specific downregulated miRNAs. Only a limited selection of pathways was enriched by novel miRNAs that were not similarly identified in pathway enrichments derived from annotated miRNAs. Several of these pathways are associated with lipid metabolism (Peroxisome in brain, ABC transporters in heart) and inflammatory signaling (Cytokine-cytokine receptor and Jak-stat signaling in plasma, CAMs in skeletal muscle), both important elements of the Weddell seal phenotype (Fig. 6).
Specific novel miRNAs also indicate potential to control key biochemical and physiological features of the Weddell seal. Novel-15-3p, highly expressed in the Weddell seal heart, is associated with cardiomyopathy through its target dystrophin. Novel-10-3p is upregulated in plasma and is predicted to target 610 mRNAs, including those with vasoactive (Nos1, Nos1ap), iron modulating (Hmox1, Hfe), and lipid metabolic functions (Lep).
Targets of differentially expressed miRNAs control physiological changes associated with elite diving
Within each tissue, age related differences in gene expression, driven by miRNA regulation detected in this dataset likely support established developmental patterns of mammals generally [35, 36]. Beyond this, we examined the 188 miRNAs that are DE between adult and pup in either the brain, heart, or muscle specifically, which may regulate the development of the diving phenotype (Additional file 3: Table S8). No age-specific differences were identified in plasma miRNAs (Table 1).
Several pathway enrichments and specific predicted target mRNAs related to hypoxia signaling were detected in this dataset. The main developmental signal linked to hypoxia responses was the expression of miR-424-3p, which is significantly reduced in the heart overall, and additionally downregulated in the adult seal compared to pup. Egln3 is a prolyl hydroxylase involved in HIF1a degradation under normoxic conditions and is predicted to be targeted by miR-424-3p. Egln3 also hydroxylates pyruvate kinase (PKM) in hypoxia, restricting glycolysis.
Each tissue type revealed DE miRNAs that target nitric oxide (NO) signaling (Fig. 7). Differentially upregulated miRNAs target the beta-1 subunit of guanylyl cyclase (Gucy1b1), which would decrease local expression of the heterodimeric enzyme in the brain, heart and muscle overall, whereas miRNAs targeting Gucy1b1 were downregulated in plasma (Fig. 7). Within the brain, an additional upregulation of novel-25-5p and novel-119-5p was detected in adult seals, which also target Gucy1b1 (Additional file 3: Table S8). Guanylyl cyclase is a major target of NO and has previously been shown to be reduced in seal tissues relative to terrestrial mammals [10]. The alpha-1 subunit (Gucy1a1) is also targeted by several DE miRNAs, including novel-128-5p (Fig. 7). Novel-128-5p is highly expressed in plasma and is also upregulated overall with age (Additional file 3: Tables S1, S7). In addition to miRNAs that may regulate endothelial NO signal transduction via guanylyl cyclase, a number of miRNAs in this dataset target the NO synthase (Nos1) directly (Fig. 7, additional file 3: Table S8), with age-related DE pointing to decreased miRNA expression in adult muscle (miRNAs 377-5p, 214-5p, cfa-543-3p) but increased expression in adult brain (miRNAs 34-5p, 1306-5p).
Contrary to central tissues such as the heart and brain, muscles experience reduced and potentially episodic blood flow during diving [37, 38]. Targets of age-increased miR-218a-3p are enriched for ubiquitin/proteasome pathways, suggesting a less robust response to protein damage with development, as well as reduced apoptosis signaling (miR-29c-3p, miR-129a-3p).
Of 137 novel miRNA loci (133 orthogroups) identified in gain and loss evolution analyses as specific to the Weddell Seal lineage, 9 mature miRNAs are DE between adults and pups in the brain, 4 are DE with age in the heart, and 6 are DE with age in muscle (Additional file 3: Table S8). Four of the novel miRNAs upregulated in the adult brain and 2 that are upregulated in the adult heart are predicted to target Vash1, an angiogenesis inhibitor selective to endothelial cells. Three miRNAs are upregulated in both the heart and skeletal muscle of adult seals and the predicted mRNA targets of these three miRNAs include Ednra, the receptor for endothelin 1 (a potent vasoconstrictor), as well as several members of the lipid transporter ABCA gene family (Abca8, Abca10, Abca12).