Integrative approach uncover microRNA interactome dysregulation in osteoarthritis cartilage

      Purpose: MicroRNAs (miRNAs) are small non-coding RNAs (∼22 nucleotides) that play an important role in gene expression regulation by binding to multiple mRNAs regulating their inhibition or degradation. Moreover, miRNAs are known to be altered in disease conditions including osteoarthritis (OA). As such, miRNAs are recognized as important targets for treatment options. Here, we uncover the dysfunctional miRNA interactome in OA articular cartilage using an integrated analysis of miRNAs and mRNA next generation sequencing datasets.
      Methods: This study includes a total of N=134 primary OA patients who underwent a joint replacement surgery, from which macroscopically unaffected (preserved) and lesioned OA articular cartilage were collected (RAAK-study). From these cartilage samples, n=64 (32 paired preserved and lesioned) were used for small non-coding RNA-sequencing, and n=70 samples (35 paired preserved and lesioned) for mRNA-sequencing (Hiseq 2500 and 4000, Illumina). In n=21 of these samples, both mRNA and miRNA sequencing data were obtained. Sequencing alignments, mapping and annotation were done with a custom in-house pipeline. The DESeq2 R-package was used for normalization and statistical analysis. Differential expression analysis of miRNA and mRNA was performed between pairs of preserved and lesioned OA cartilage and molecules with false discovery rate (FDR) < 0.05 and absolute fold change (FC) > 1.5 were considered significantly differentially expressed (DE). To construct a miRNA-mRNA co-expression network, we calculated the Pearson correlation between DE miRNAs and DE genes. Only correlations with |r2| > 0.5 and P-value < 0.05 were considered for further analysis. To prioritize correlated miRNA-mRNA target genes interactions, we integrated predicted interactions from TargetScan and microT-CDS with experimentally validated miRNA-mRNA interactions from miRTarbase and Tarbase v7 using the miRComb R package. Finally, pathway enrichment analysis was performed on correlated-target genes using DAVID. Pathways with FDR < 0.05 were considered significantly enriched.
      Results: We identified 132 miRNAs and 385 mRNAs with significant DE (FDR < 0.05, FC > 1.5.) between preserved and lesioned OA articular cartilage (Figure 1A). Of the miRNAs, 31 were previously reported to be involved in OA, including, miR-502-3p (P = 0.04, FC = 0.7) and miR-23a-3p (P = 0.03, FC = 1.51). Additionally, we identified 101 DE miRNAs that were not previously associated to OA articular cartilage. For example, miR-107, which was previously only linked to obesity and to insulin sensitivity, is significantly upregulated in lesioned OA cartilage (P = 4.4x10-5, FC = 1.68). To uncover the miRNA-interactome we identified 98 mRNAs to be significantly correlated (FDR<0.05, |r2| >0.5) with 55 DE miRNAs. The miRNA-interactome included miRNAs that target multiple mRNAs, such as miR-107 which is negatively correlated and predicted to target the same 6 DE genes (Figure 1B). Vice versa, BMPR1B gene (P = 1.27 x 10-7, FC = 0.3) is regulated by miR-502-3p and miR-23a-3p and miR-339-5p (Figure 1B). Finally, upon performing pathway analyses on the genes targeted by miRNAs in cartilage, we revealed that these targeted genes were significantly (FDR < 0.05) enriched in 12 pathways, including extracellular matrix (FDR = 0.007) and chondrogenesis (FDR = 0.04).