Identification of differentially expressed microRNAs in knee anterior cruciate ligament tissues surgically removed from patients with osteoarthritis.

The degradation of cruciate ligaments is frequently observed in degenerative joint diseases, such as osteo-arthritis (OA). The present study aimed to identify the differentially expressed microRNAs (miRNAs or miRs) in knee anterior cruciate ligament (ACL) tissues derived from patients with OA and in health subjects (non-OA). By using Affymetrix miRNA 4.0 microarrays, a total of 22 miRNAs (including let-7f-5p, miR-26b-5p and miR-146a-5p) were found to be upregulated, while 17 (including miR-18a-3p, miR-138-5p and miR-485-3p) were downregulated in the osteoarthritic ACL tissues (fold change ≥2, P-value <0.05). The expression levels of 12 miRNAs were validated by quantitative PCR, and the corresponding results revealed an excellent correlation with the microarray data (R2=0.889). Genes (such as a disintegrin and metalloproteinase domain with thrombospondin type-1 motifs, bone morphogenetic protein-2, runt related transcription factor-2, collagen-1A1 and 2, interleukin-6 and transforming growth factor-β) involved in cartilage development and remodeling, collagen biosynthesis and degradation, inflammatory response and extracellular matrix homeostasis were predicted as potential targets of the dysregulated miRNAs. Moreover, a large set of putative genes were enriched in OA pathogenesis‑associated pathways (such as mitogen-activated protein kinase and vascular endothelial growth factor signaling pathway). Collectively, the data from our study provides novel insight into the ligament injury-related miRNA dysregulation in patients with OA.

Introduction

Osteoarthritis (OA) is a degenerative joint disease characterized by the destruction of articular cartilage, intraarticular inflammation and pathological alterations in peri-articular and subchondral bone (1,2). Various factors are involved in the pathogenesis of OA, including age (3), a history of diabetes, cancer or cardiovascular diseases (4), mechanical influences (5) and genetic factors (6). There is no disease-modifying treatment for the onset or progression of OA and associated structural damage, and the current treatments aim at relieving the symptoms (7). Therefore, the identification of novel molecules involved in the pathogenesis of OA is urgently required, and will provide basis for the development of therapies for OA.

MicroRNAs (miRNAs or miRs) are a category of non-coding RNAs 22–25 nt in length (8). As the key gene regulators, miRNAs directly bind to their target messenger RNAs (mRNAs) in a sequence-specific manner to facilitate degradation of the transcripts and to inhibit the protein translation (8). Differential expression profiles of certain miRNAs in cancers at different stages suggests that miRNAs are novel biomarkers for disease diagnostics (9). The application of microarray technology enables the detection of the expression levels of thousands of miRNAs simultaneously within tens of samples processed in a single experiment (10). The dysregulation of miRNAs has been found in tissue samples derived from patients with OA in a number of previous studies, including let-7 family miRNAs (11), miR-149 (12), miR-21 (13) and miR-24 (14). Most of the earlier studies compared miRNA expression in the injured cartilage and synovium between patients with OA and normal controls (15–17); however, changes in cruciate ligment have been less studied. The cruciate ligament is a collagenous tissue for structural support and provides proprioception to the body by mediating knee kinesthesia (18). Of note, the degradation of the cruciate ligaments frequently occurs in osteoarthritic knees (18,19). The present study was therefore conducted to analyze the miRNA expression profiles in anterior cruciate ligament (ACL) tissues surgically removed from patients with OA and control subjects by using miRNA microarray analysis. In addition, the biological functions and pathways affected by the differentially expressed miRNAs were analyzed.

Materials and methods

Sample recruitment and RNA extraction

Osteoarthritic ACL samples were surgically removed from 3 patients (64.67±3.06 years of age, Kellgren-Lawrence grade III–IV) during knee replacement surgery at Shengjing Hospital of China Medical University, Shenyang, China. Samples derived from 3 patients without OA who encountered ACL rupture were used as controls. The present research protocol was approved by the Institutional Review Board of China Medical University, and written informed consent was obtained from each participant prior to obtaining the samples. Total RNA was extracted from the ACL tissue samples using the total RNA purification kit (Norgen Biotek Corp., Thorold, ON, Canada), quantified on a NanoDrop ND-2100 spectrophotometer (Thermo Fisher Scientific, Inc., Pittsburgh, PA, USA), and assessed on an Agilent 2100 Bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA, USA).

miRNA microarray procedures

The RNA samples were tailed with Poly(A) and labeled with biotin using FlashTag™ biotin HSR ligation mix (Affymetrix, Inc., Santa Clara, CA, USA) according to the manufacturer's instructions. The labeled RNA samples were hybridized onto the Affymetrix miRNA 4.0 arrays on a hybridization oven 645, washed and stained on fluidics station 450, and then scanned with a Scanner 3000 (all from Affymetrix, Inc.).

Data analysis

Array images were analyzed with GeneChip Command Console software (version 4.0; Affymetrix, Inc.) to generate raw data. The obtained raw data were first normalized with robust multi-array average (RMA) using Expression Console software (version 1.3.1; Affymetrix, Inc.) and then analyzed with GeneSpring software (version 12.5; Agilent Technologies, Inc.). Principal component analysis (PCA) is a mathematical algorithm that is performed to reduce the data dimensionality while retaining most of the variation in the data set (20). Differentially expressed miRNAs were identified by evaluating the fold change (FC). miRNAs with an FC ≥2 and a P-value <0.05 (t-test) were considered as differentially expressed. Hierarchical clustering was performed to analyze the distinguishable miRNA expression patterns among the samples. Genes targeted by the identified differentially expressed miRNAs were shown as the intersection of Targetscan, PITA and microRNA.org databases (GeneSpring software, version 12.5). These putative target genes were subjected to Gene Ontology (GO) biological process annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis using FunNet algorithm. The P-value was calculated with a unilateral Fisher's exact test and corrected by false discovery rate (FDR). GOs and pathways with P-value <0.05 and FDR <0.05 were considered as significant. To display the combinatorial interactions between miRNA pairs and their shared targets, genes with P-value <0.05 (calculated by hypergeometric distribution) were presented using Cytoscape software.

Quantitative PCR

Quantitative PCR was performed on cDNA synthesized from the same RNA samples used in the prior microarray analysis. Primers used in this study were listed in Table I. The expression levels of 6 upregulated miRNAs (hsa-let-7f-5p, hsa-let-7g-5p, hsa-miR-146a-5p, hsa-miR-146b-3p, hsa-miR-26b-5p and hsa-miR-335-5p) and 6 downregulated miRNAs (hsa-miR-18a-3p, hsa-miR-485-3p, hsa-miR-665, hsa-miR-675-5p, hsa-miR-1207-5p and hsa-miR-138-5p) were determined on the Exicycler™ 96 (Bioneer, Daejeon, Korea) using SYBR-Green (Solarbio, Beijing, China). Triplicate reactions were performed. The data were analyzed by the comparative threshold cycle (Ct) method. U6 was used as the endogenous control.

Table I

Primers used in this work.

MiRBase accession number	Name	Sequence information (5′→3′)
MIMAT0000067	hsa-let-7f-5p
	RT primer	GTTGGCTCTGGTGCAGGGTCCGAGGTATTCGCACCAGAGCCAACAACTAT
	Forward real-time PCR primer	CGCGGCTGAGGTAGTAGATTGT
MIMAT0000414	hsa-let-7g-5p
	RT primer	GTTGGCTCTGGTGCAGGGTCCGAGGTATTCGCACCAGAGCCAACAACTGT
	Forward real-time PCR primer	CGGTCGTGAGGTAGTAGTTTGT
MIMAT0000449	hsa-miR-146a-5p
	RT primer	GTTGGCTCTGGTGCAGGGTCCGAGGTATTCGCACCAGAGCCAACAACCCA
	Forward real-time PCR primer	GCGAGGTGAGAACTGAATTCCA
MIMAT0004766	hsa-miR-146b-3p
	RT primer	GTTGGCTCTGGTGCAGGGTCCGAGGTATTCGCACCAGAGCCAACCCAGAA
	Forward real-time PCR primer	GACTGCCCTGTGGACTCAGTTC
MIMAT0000083	hsa-miR-26b-5p
	RT primer	GTTGGCTCTGGTGCAGGGTCCGAGGTATTCGCACCAGAGCCAACACCTAT
	Forward real-time PCR primer	CGCGGCTTCAAGTAATTCAGG
MIMAT0000765	hsa-miR-335-5p
	RT primer	GTTGGCTCTGGTGCAGGGTCCGAGGTATTCGCACCAGAGCCAACACATTT
	Forward real-time PCR primer	CGCAGCTCAAGAGCAATAACGA
MIMAT0002891	hsa-miR-18a-3p
	RT primer	GTTGGCTCTGGTGCAGGGTCCGAGGTATTCGCACCAGAGCCAACCCAGAA
	Forward real-time PCR primer	CGACTACTGCCCTAAGTGCTC
MIMAT0002176	hsa-mir-485-3p
	RT primer	GTTGGCTCTGGTGCAGGGTCCGAGGTATTCGCACCAGAGCCAACAGAGAG
	Forward real-time PCR primer	CTGCTGTCATACACGGCTCTC
MIMAT0004952	hsa-miR-665
	RT primer	GTTGGCTCTGGTGCAGGGTCCGAGGTATTCGCACCAGAGCCAACAGGGGC
	Forward real-time PCR primer	CAGTTAACCAGGAGGCTGAGG
MIMAT0004284	hsa-miR-675-5p
	RT primer	GTTGGCTCTGGTGCAGGGTCCGAGGTATTCGCACCAGAGCCAACCACTGT
	Forward real-time PCR primer	CTATAATGGTGCGGAGAGGGCC
MIMAT0005871	hsa-miR-1207-5p
	RT primer	GTTGGCTCTGGTGCAGGGTCCGAGGTATTCGCACCAGAGCCAACCCCCTC
	Forward real-time PCR primer	CTTATTGGCAGGGAGGCTG
MIMAT0000430	hsa-miR-138-5p
	RT primer	GTTGGCTCTGGTGCAGGGTCCGAGGTATTCGCACCAGAGCCAACCGGCCT
	Forward real-time PCR primer	CGGTGCAGCTGGTGTTGTGAAT
	Universal reverse primer	GTGCAGGGTCCGAGGTATTC

Results

PCA distinguishes patients with OA from control subjects

PCA was performed in the 6 knee ACL tissues based on the microarray data, and the corresponding results revealed that the ACL samples from the patients with OA could be distinguished from those of the control subjects (Fig. 1). The above results suggested that the specimens used in this study were properly prepared and could be classified into two distinct groups.

Figure 1

PCA of the miRNA profiles in knee ACL samples from patients with OA and control subjects. Each sample was assayed using an Affymetrix miRNA 4.0 microarray. The red dots represent samples from patients with OA, whereas the blue dots represent samples from the control subjects. PCA, principal component analysis; miRNA, microRNA; ACL, anterior cruciate ligament; OA, osteoarthritis.

Identification of differentially expressed miRNAs in knee ACL tissues

Our data indicated that 22 miRNAs were upregulated and 17 miRNAs were downregulated in the osteoarthritic ACL tissues (FC ≥2 and P-value <0.05). Twelve miRNAs were selected for further validation regarding their expression levels (such as hsa-miR-138-5p, hsa-miR-26b-5p, hsa-miR-665), their previous correlations with OA (such as hsa-miR-146a-5p, hsa-let-7g-5p, hsa-let-7f-5p), or data from the following analysis (such as hsa-miR-146b-3p, hsa-miR-1207-5p). Log2 results of the miRNA expression levels from the microarray analysis and the quantitative PCR analysis revealed an excellent correlation (R2=0.889; Fig. 2A). All analyzed human miRNAs were presented in a volcano plot (Fig. 2B), and the dysregulated ones were assessed via hierarchical clustering analysis (Fig. 2C). Results from hierarchical clustering analysis revealed that the ACL tissue samples were divided into two distinct clusters based on their pathological statuses. Collectively, these results implied that the microarray data reflected the reliable miRNA expression patterns in ACL tissues and that the samples from same situation clustered together.

Figure 2

Differentially expressed miRNAs in knee ACL tissues of patients with OA. (A) Log2 FC for microarray data (y-axis) is plotted against log2 FC for quantitative PCR data (x-axis) for each gene. (B) Volcano plots of miRNA microarray data (x-axis, log2 FC; y-axis, negative log10 P-value). The red plots depict differentially expressed miRNAs with an FC ≥2 and P<0.05 between the two groups. (C) Hierarchical cluster analysis of the miRNA microarray data. Each column depicts a single tissue sample and each row represented a miRNA. Red and green indicate high and low expression levels, respectively, whereas black indicates the mean expression levels. Distance metric, pearson centered; linkage rule, average. miRNA, microRNA; ACL, anterior cruciate ligament; FC, fold change; OA, osteoarthritis.

Microarray-based GO and KEGG pathway annotations

Three data bases predicted a total of 5,356 genes as putative targets for the differentially expressed miRNAs (Fig. 3A). Genes involved in cartilage remodeling (21), collagen biosynthesis (22), extracellular matrix (ECM) homeostasis (23) and inflammation (24) are summarized in Table II. Moreover, the GO annotation (at the biological process level) and KEGG pathway analysis of all putative genes revealed that these genes were enriched in 41 GO items and 23 KEGG pathways (data not shown). As indicated in Table III, several essential biological processes, including DNA-dependent regulation of transcription, signal transduction, multicellular organismal development, were affected by the differentially expressed miRNAs. Additionally, a large set of genes implicated in mitogen-activated protein kinase (MAPK), vascular endothelial growth factor (VEGF), protein kinase C (PKC)-mitogen-activated protein kinase (MEK), phosphatidylinositol 3-kinase (PI3K)-protein kinase B (AKT), as well as the WNT signaling pathway were mediated by the dysregulated miRNAs (Fig. 3B and C; pathway ID, 05200 for Fig. 3C). The above results provided useful information for us to understand how cruciate ligament injuries develop in OA patients.

Figure 3

KEGG pathway annotations of the putative target genes. (A) Genes targeted by the differentially expressed miRNAs are predicted with Targetscan, PITA and microRNA.org databases, and the corresponding numbers are shown in the venn diagram. (B) Putative targeted genes are annotated into several KEGG pathways. Pathways with P-value <0.05 and FDR <0.05 are selected and presented here (the top 10 pathways). (C) Pathways in cancer (ID, 05200). Genes mediated by the differentially expressed miRNAs are marked in red. miRNAs, microRNAs; FDR, false discovery rate; OA, osteoarthritis; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Table II

Prediction of target genes potentially related to osteoarthritis.

miRNAs	Gene symbols	Function
hsa-let-7f/7g-5p	BMP2, COL1A1, COL1A2, COL3A1, COL4A1, COL4A6, COL5A2, COL14A1, COL15A1, COL24A1, TGFBR1, ADAMTS8, IL6, IL10, IL13, HMGA1, HMGA2	Cartilage development and remodeling
hsa-miR-10a-5p	COL4A4, COL24A1, ADAMTS4
hsa-miR-26b-5p	CILP, COL1A2, COL9A1, COL10A1, COL11A1, ADAMTS19, HMGA1, HMGA2, IL6, IL1RAP	Collagen biosynthesis and degradation
hsa-miR-138-5p	MMP16, ADAMTSL3, IL6R, IL1RAP
hsa-miR-146a-5p	CCL5, CXCR7, ADAMTS3, ADAMTS18
hsa-miR-146b-3p	COL11A1, MMP24, VEGFA
hsa-miR-335-5p	COL5A1, COL6A3, COL19A1, ADAMTS19, CCL5	ECM homeostasis
hsa-miR-542-5p	ADAMTS8
hsa-miR-485-3p	COL12A1, TGFB3, MMP20, ADAMTS3	Inflammatory response
hsa-miR-572	COL7A1
hsa-miR-665	COL8A2, TGFBR1, TGFBR2, ADAMTS8, CXCL11, CXCL12
hsa-miR-1207-5p	COL9A2, TGFBR1, ADAMTS10, ADAMTS19
hsa-miR-1254	RUNX2, TGFBR3, ADAMTSL5

OA, osteoarthritis; miRNAs, microRNAs; BMP, bone morphogenetic protein; RUNX, runt related transcription factor; COL, collagen; CILP, cartilage intermediate layer protein; HMGA, high mobility group AT-hook; IL, interleukin; ADAMTS, a disintegrin and metalloproteinase domain with thrombospondin type-1 motifs; CCL, chemokine (C-C motif) ligand; CXCL, chemokine (C-X-C motif) ligand; CXCR, chemokine (C-X-C motif) receptor; TGFB, transforming growth factor β; TGFBR, TGFB receptor; MMP, matrix metalloproteinase.

Table III

Identified biological process GO terms for the differentially expressed miRNAs (top 10).

GO ID	GO term	List hits	P-value
GO:0006355	Regulation of transcription, DNA-dependent	491	5.38E-09
GO:0007165	Signal transduction	289	8.11E-05
GO:0007275	Multicellular organismal development	237	5.38E-05
GO:0006351	Transcription, DNA-dependent	166	1.24E-08
GO:0006468	Protein phosphorylation	153	1.25E-08
GO:0007155	Cell adhesion	151	3.54E-05
GO:0045944	Positive regulation of transcription from RNA polymerase II promoter	138	4.43E-06
GO:0007399	Nervous system development	131	8.66E-10
GO:0007049	Cell cycle	121	5.31E-06
GO:0007411	Axon guidance	114	4.42E-12

GO, Gene Ontology; miRNAs, microRNAs; list hits, the number of genes annotated by the GO biological process category or annotation cluster within the analyzed list of target genes; P-value, the significance P-value of the gene enrichment in the GO biological process category or annotation cluster, calculated with the unilateral Fisher's exact test and corrected with the false discovery rate (FDR).

Establishment of miRNA-gene regulatory network

A large set of genes were predicated as targets for the differentially expressed miRNAs through the Targetscan, PITA and microRNA.org data bases. In order to visualize and integrate the interactions between the dysregulated miRNAs and their targets, a miRNA-gene regulatory network was established by using Cytoscape software (Fig. 4). Our results revealed that the differentially expressed miRNAs may function in combination to exert effects on their target genes.

Figure 4

The microRNA (miRNA)-gene network. The association between miRNA pairs with their targeted genes is obtained on the basis of the cumulative hypergeometric statistic, and presented by the Cytoscape software.

Discussion

miRNAs play crucial roles in mediating chondrogenesis, and are considered to link to the pathogenesis of cartilage-related diseases, including OA (25). In this study, miRNA microarray was performed to compare the miRNA expression levels in knee ACL tissues from patients with OA to those of the controls. Appropriate grouping of the 6 ACL samples was confirmed by PCA and heatmap data. We found that 22 miRNAs were upregulated and 17 were downregulated in the osteoarthritic ACL tissues. Additional bioinformatics was performed to analyze the biological processes and pathways that were affected by the identified differentially expressed miRNAs. The obtained data enhanced our understanding of the roles of the dysregulated miRNAs in OA pathogenesis.

Reportedly, let-7 miRNAs can regulate skeletal development by orchestrating the proliferation and differentiation of chondrocytes (11). The enforced overexpression of Lin28a, a let-7 inhibitor, has been shown to accelerate cartilage regrowth in a model of tissue injury (26). It is likely that the abnormal upregulation of let-7 miRNAs contributes to the degeneration of articular cartilage. In this study, to the best of our knowledge, we demonstrate for the first time that the expression of let-7f-5p and let-7g-5p was increased by 2.04- and 1.68-fold (log2FC) in the osteoarthritic ACL tissues, respectively. Though all let-7 family members share the identical seed region (GAGGUAG) (27), only these two let-7 members were identified to be dysregulated in OA-affected ligaments. Bone morphogenetic protein (BMP)2 has been reported to promote osteogenesis (28). Apart from its role in bone formation, the pre-injection of recombinant human BMP2 in the semitendinosus tendon enables successful ACL reconstruction following injury (29), suggesting a beneficial role of BMP2 in ligament injury. miR-140-5p is a potent regulator of BMP2 (30), and its expression is markedly reduced in osteoarthritic articular cartilage tissues (31), but not in ACL tissues, as evidenced by our microarray data. Of note, we found that BMP2 is a possible target for let-7f/7g-5p (Table II), although the interaction between them has not been entirely clarified. Moreover, apart from BMP2, other factors related to collagen biosynthesis and degradation and inflammatory response, such as transforming growth factor β receptor 1 (TGFβR1), various types of collagens (COL1A1 and COL1A2) and interleukins (IL)-6 were also putative targets for let-7f/7g-5p. To address the roles of let-7f/7g in osteoarthritic ligament lesion, their targets should also be taken into consideration.

Formation and degradation of collagens and ECM proteins are mediated by miR-26 family members (32). miR-26b is suggested to contribute to rheumatoid arthritis regarding to its elevation in IL-17 producing T cells (33). Such findings indicate that miR-26b may participate in inflammatory diseases in the joints. A significant upregulation of miR-26b-5p (previously miR-26b) was found in osteoarthritic ACL tissues in the present study. Although several putative targets of miRNA-26b, such as high mobility group AT-hook 1 (HMGA1 and HMGA2), cartilage intermediate layer protein (CILP), as well as a variety of collagens (COLs) are implicated in the development and progression of OA (34–36), the direct correlation of miRNA-26b dysregulation with OA has not been fully elucidated, and requires for further exploration.

miR-146a controls knee joint homeostasis by balancing inflammatory responses in cartilage (37). Its expression is increased in articular cartilage and/or synovium derived from patients with OA (38,39). Our results were consistent with these earlier findings by showing a significant upregulation of miR-146a-5p (previously miR-146a) in osteoarthritic knee ACL tissues. Studies on the correlation between OA pathogenesis and miR-146b-3p are limited. Our data indicated that miR-146b-3p was also overexpressed in osteoarthritic ACL tissues. Several genes associated to ECM homeostasis and inflammation such as matrix metalloproteinase (MMP)24 and VEGFA in OA (40,41) were predicted as targets for miR-146b-3p.

A disintegrin and metalloproteinase domain with thrombospondin type-1 motifs (ADAMTS) are a new family of metalloproteases that play important roles in physiological and pathological conditions (42,43). Previous studies have demonstrated that ADAMTS7 overexpression leading to the increased expression of tumor necrosis factor (TNF)-α and MMPs contributes OA development (44), while the knockdown or knockout of ADAMTS4 and/or ADAMTS5 prevents OA progression (45,46). These studies suggest that ADAMTS may be the potential molecular targets for the prevention and treatment of OA. In this study, we found that ADAMTS3, 4, 8, 10, 18, 19, ADAMTS-like-3, -5 were the putative targets for several differentially expressed miRNAs, including let-7f/7g-5p, miR-146a-5p, miR-1207-5p (Table II). Investigations of the interaction between these dysregulated miRNAs and their target ADAMTS will help to understand the mechanisms through which OA develops and progresses.

Several essential biological processes, such as the DNA-depen dent regulation of transcription, signal transduction and multicellular organismal development, are affected by miRNAs with differential expression levels in OA as indicated in GO annotation. To provide an overall understanding of the association between the dysregulated miRNAs and OA pathogenesis, KEGG pathway analysis was further performed. We found that several pathways enriched by the putative target genes were essential for OA pathogenesis. For instance, a study from Prasadam et al demonstrated that p38 MAPK phosphorylation was decreased in OA-affected chondrocytes as compared to normal chondrocytes, and that the inactivation of p38 signaling leads to OA-like changes in rats (47). In addition, activation of VEGF signaling has been suggested to contribute to synovial inflammation during the progression of OA (48).

In conclusion, our study revealed that 39 miRNAs were differentially expressed in knee ACL tissues from patients with OA. The functional bioinformatic analyses suggest that the dysregulated miRNAs may regulate cartilage development and remodeling, collagen biosynthesis and degradation, ECM homeostasis and pathology by interacting with their targets. Collectively, our study provides novel insight into the ligament injury-related miRNA dysregulation in patients with OA.