Whole-Transcriptome Sequence of Degenerative Meniscus Cells Unveiling Diagnostic Markers and Therapeutic Targets for Osteoarthritis

Meniscus plays an important role in joint homeostasis. Tear or degeneration of meniscus might facilitate the process of knee osteoarthritis (OA). Hence, to investigate the transcriptome change during meniscus degeneration, we reveal the alterations of messenger RNA (mRNA), microRNA (miRNA), long noncoding RNA (lncRNA), and circular RNA (circRNA) in meniscus during OA by whole-transcriptome sequence. A total of 375 mRNAs, 15 miRNAs, 56 lncRNAs, and 90 circRNAs were significantly altered in the degenerative meniscus treated with interleukin-1β (IL-1β). More importantly, highly specific co-expression RNA (ceRNA) networks regulated by lncRNA LOC107986251-miR-212-5p-SESN3 and hsa_circ_0018069-miR-147b-3p-TJP2 were screened out during IL-induced meniscus degeneration, unveiling potential therapeutic targets for meniscus degeneration during the OA process. Furthermore, lipocalin-2 (LCN2) and RAB27B were identified as potential biomarkers in meniscus degeneration by overlapping three previously constructed databases of OA menisci. LCN2 and RAB27B were both upregulated in osteoarthritic menisci and IL-1β-treated menisci and were highly associated with the severity of OA. This could introduce potential novel molecules into the database of clinical diagnostic biomarkers and possible therapeutic targets for early-stage OA treatment.

Introduction

Meniscus, a semi-lunar, pale white tissue located between the tibial plateau and femoral condyle, is of great importance to the structure, stability, and biomechanical function of the knee joint (Makris et al., 2011). During attempts to treat common sports injuries such as meniscal tears, previous studies concerning the menisci have been largely focused on applying tissue-engineering techniques to develop qualified meniscus substitutes that possess similar meniscal anatomic structure, biocompatibility, and function (Filardo et al., 2019; Kwon et al., 2019; Shimomura et al., 2019; Zhang et al., 2019). Recently, several studies have discovered that meniscal tears, degeneration, and total meniscectomy are of high relevance with the onset of osteoarthritis (OA), a common chronic disease that mostly occurs in the knee joint (Berthiaume et al., 2005; Hunter et al., 2006; Murphy et al., 2019). The possible “meniscal pathway” to knee OA might be attributed to abnormal biomechanical stress caused by meniscal trauma or dissection (Englund et al., 2012). However, little research has been performed to investigate the specific mechanism underlying meniscal pathology and OA. Therefore, it is necessary to study specific meniscal degenerative mechanisms to cope with the diagnosis and treatment of early OA.

To study the potential regulating mechanism between transcriptional and post-transcriptional modification, bioinformatics has been widely applied. Multiple novel methods have been used to demonstrate the mechanisms responsible for various diseases, for example, with the use of whole-transcriptome sequencing, with which Lei et al. were able to discover the comprehensive circular RNA (circRNA) profile of peripheral blood mononuclear cells in hepatocellular carcinoma (HCC) patients and identified circ_0000798 as a characteristic biomarker for HCC patients (Lei et al., 2019). In orthopedic fields, whole-transcriptome sequencing (Li et al., 2019) and single-cell sequencing (Ji et al., 2019) on human primary OA chondrocytes also assisted with the understanding of the degenerative mechanism of cartilage. Previous studies have also screened out the potential messenger RNA (mRNA), microRNA (miRNA), and long noncoding RNA (lncRNA) in knee cartilage, which might possess cartilage degeneration during OA process (Chen and Chen, 2020; Qi et al., 2020). Recently, with the aid of single-cell sequencing, we were able to identify normal and degenerative meniscus cell types and superficially uncovered that the transition from a meniscus fibrocartilage progenitor (FCP) cell to a meniscus degenerative progenitor cell (DegP) is possibly the key determining factor of the OA meniscus (Sun et al., 2020). However, the molecular mechanism underlying this transition remains unknown.

In this study, we performed whole-transcriptome sequencing on degenerative meniscus with or without interleukin-1β (IL-1β) treatment as inflammatory chemokines, including IL-1β, have been studied as necessary catabolic mediators in OA that are also applicable to the meniscus (McNulty et al., 2013; Cook et al., 2018). Our previous study has also shown that with 48-h IL-1β (5 ng/ml) treatments, the number of FCP cells decreased while DegP increased, possibly causing the transition from FCP to DegP (Sun et al., 2020). Hence, by using IL-1β as an OA inducer in meniscus, we were able to examine the mRNA, miRNA, lncRNA, and circRNA expression profiles in degenerative menisci to identify the characteristics of transcriptional and post-transcriptional differences during the OA degenerative process. Moreover, we overlapped three databases of degenerative menisci to select highly specific biomarkers in the meniscus for diagnosing early-stage OA, including our previous single-cell sequencing on normal menisci and OA degenerative menisci (Sun et al., 2020), whole-transcriptome sequence of IL-1β-abundant menisci, and RNA-seq of control menisci as well as OA degenerative meniscus.

Materials and Methods

Isolation and Culture of Human Meniscus Cells

OA meniscus samples were dissected from 10 OA patients who had the indication of total knee arthroplasty (TKA), and the patients who participated in this program offered written informed consent. Healthy meniscus samples were collected from patients who underwent amputation and did not have OA or rheumatoid arthritis. The enrolled criteria included classic clinical history, pains, signs of dyskinesia, and X-ray imaging. The average age and Kellgren-Lawrence grading scores of the patients are listed in Supplemental Table S2. The exclusion criteria and procedures for sample collection and examination were carried out as described in previous studies (Meng et al., 2018). Afterwards, the menisci were cut into slices and digested with 2 mg/ml of collagenase P for 8–12 h and then implanted into medium containing DMEM/Nutrient Mixture F-12 (Gibco Life Technologies, Grand Island, NY, United States ), 5% fetal bovine serum (FBS; Gibco Life Technologies), and 100 IU/ml of penicillin (PS; Gibco Life Technologies). The meniscus cells were cultured in 6-well plates at 37°C in a humidified atmosphere of 5% CO₂ and 1% oxygen. The cell density was about 1 × 10⁷ per plate.

Inflammatory Stimulation With Interleukin-1β

For whole-transcriptome sequence, three OA meniscus samples dissected from OA patients were collected and plated for cell culture, named OA004, OA006, and OA008. After the meniscus cells fully adhered to the plate and showed 90% cellular confluency in the 6-well plate, we added 5 ng/ml IL-1β in three wells in each sample, named OA004_IL-1B, OA006_IL-1B, and OA008_IL-1B, while simultaneously added refreshed culture medium as control group (OA004_NC, OA006_NC, and OA008_NC). Three samples were treated with 5 ng/ml of IL-1β to simulate OA inflammatory pathology (OA004_IL-1B, OA006_IL-1B, and OA008_IL-1B), while the other three samples were replaced with refreshed medium instead (OA004_NC, OA006_NC, and OA008_NC). All samples were then cultured at 37°C in a humidified atmosphere of 5% CO₂ for 48 h.

Total RNA Extraction

For RNA sequence, four healthy meniscus samples were collected from patients who underwent amputation due to severe femoral fracture who did not have OA or rheumatoid arthritis, and four OA meniscus samples were collected from patients who had the indication of TKA. Total RNA was extracted using TRIzol reagent kit (Invitrogen, Carlsbad, CA, United States) according to the manufacturer’s protocol. We used TRIzol (Invitrogen, Carlsbad, CA, United States) to extract total RNA from each meniscus cell, following the manufacturer’s protocol. The RNA quality was checked by an Agilent 2,200 (Agilent Technologies, Santa Clara, CA, United States) and kept at −80°C, and only samples with an RNA integrity number (RIN) value > 7.0 were used for the cDNA library construction.

cDNA Library Construction

The cDNA libraries were constructed for each pooled RNA sample using the NEBNext^® Ultra™ Directional RNA Library Prep Kit for Illumina (San Diego, CA, United States) according to the manufacturer’s instructions. Generally, the protocol consists of the following steps: depletion of rRNA and fragmented into 150–200 bp using divalent cations at 94°C for 8 min. We further used Dnase I to eradicate contamination after we wiped off ribosome RNA. The cleaved RNA fragments were reverse-transcribed into first-strand cDNA and second-strand cDNA synthesis; fragments were end repaired, A-tailed, and ligated with indexed adapters. Target bands were harvested through AMPure XP Beads (Beckman Coulter, Brea, CA, United States). The products were purified and enriched by PCR to create the final cDNA libraries and quantified by Agilent 2,200. The tagged cDNA libraries were pooled in equal ratio and used for 150-bp paired-end sequencing in a single lane of the Illumina HiSeq X Ten.

The sequencing library of miRNA was prepared from total RNA by using NEBNext Small RNA Library Prep Set for Illumina (NEB) according to the manufacturer’s instructions. Briefly, RNA was ligated with 5′-RNA and 3′-RNA adapters, reversely transcribed into cDNAs, and PCR amplified. The PCR products were size selected and sequenced on HiSeq X Ten platform.

RNA Sequencing Mapping

Mapping of paired-end reads: Before read mapping, clean reads were obtained from the raw reads by removing the adaptor sequences, reads with >5% ambiguous bases (noted as N), and low-quality reads containing more than 20% of bases with qualities of <20. The clean reads were then aligned to human genome [version: GRCh38 National Center for Biotechnology Information (NCBI)] using the hisat2 (Kim et al., 2015). HTseq (Anders et al., 2015) was used to get gene counts, and RPKM method was used to determine the gene expression. The clean reads of miRNA library were mapped to Human miRNA database (miRBase v22.0) to achieve the miRNA expression.

Circular RNA Identification and Quantification

The pipeline “acfs,” which was publicly available at https://code.google.com/p/acfs/, was used to identify circRNA in each sample including the following steps (You et al., 2015): Unmapped Reads Collection: BOWTIE2 version 2.2.5 (Langmead and Salzberg, 2012) was used as the mapping method to the respective reference genome [GRCH37.p13 NCBI] utilizing the parameter bowtie2 --end-to-end --sensitive --mm --phred33 --fr --rg-id S13171 --rg SM:S13171 --rg LB:S13171 --rg PL:Illumina -p 8 -X 500 -k 4 -x.)

Circular RNA Identification

Unmapped reads were collected to identify the circRNA utilizing BWA mem (bwa mem -t 1 -k 16 -T 20): partial alignments of segments within a single read that mapped to 1) regions on the same chromosome and no more than 1 Mb away from each other 2) on the same strand 3) but in reverse order were retained as candidates supporting head-to-tail junction. The strength of potential splicing sites supported by these candidate head-to-tail junction reads was then estimated using MaxEntScan33. The exact junction site was determined by selecting the donor and acceptor sites with the highest splicing strength score. Candidate circRNAs were reported if the head-to-tail junction was supported by at least two reads and the splicing score was greater than or equal to 10.

Expression Analysis

To estimate the expression of circRNA, we re-aligned all the unmapped reads to the circRNA candidates by using the BWA-mem under the following parameter: bwa mem -t 1 -k 16 -T 20. As for most of the circRNAs, there is no direct evidence for their exact sequence: we filled in the sequence using existing exon annotation. Sequence at the 5′ end was concatenated to the 3′ end to form circular junctions. Reads that mapped to the junction (with an overhang of at least 6 nt) were counted for each candidate.

Dif-Gene-Finder

We applied EBSeq (Leng et al., 2013) algorithm to filter the differentially expressed genes, after the significant analysis, p-value, and false discovery rate (FDR) analysis under the following criteria (Benjamini et al., 2001):

MiRNA under the following criteria: 1) fold change ≥2 or ≤0.5; 2) FDR ≤ 0.05.

mRNA under the following criteria: 1) fold change ≥2 or ≤0.5; 2) FDR ≤ 0.05.

NcRNA under the following criteria: 1) fold change ≥2 or ≤0.5; 2) FDR ≤ 0.05.

CircRNA under the following criteria: 1) fold change ≥2 or ≤0.5; 2) FDR ≤ 0.05.

Gene Ontology Analysis

Gene Ontology (GO) analysis was performed to facilitate elucidating the biological implications of unique genes in the significant or representative profiles of the target gene of the differentially expressed miRNA in the experiment. We downloaded the GO annotations from NCBI (http://www.ncbi.nlm.nih.gov/), UniProt (http://www.uniprot.org/), and the GO (http://www.geneontology.org/). Fisher’s exact test was applied to identify the significant GO categories, and FDR was used to correct the p-values.

Pathway Analysis

Pathway analysis was used to find out the significant pathway of the differential genes according to Kyoto Encyclopedia of Genes and Genomes (KEGG) database. We turn to Fisher’s exact test to select the significant pathway, and the threshold of significance was defined by p-value and FDR.

GO-Tree

The GO is structured as a directed acyclic graph, and each term has defined relationships to one or more other terms. GO-Tree is built based on the GO directed acyclic graph to provide user-friendly data navigation and visualization. We selected the significant GO-Term (p-value < 0.01) in GO analysis based on the up and down differentially expressed genes to construct the GO-Tree to summarize the function affected in the experiment (Zhang et al., 2004).

Path-Act-Network

KEGG (Ogata et al., 1999) included metabolism, membrane transport, signal transduction, and cell cycle pathways. We picked the genes in enriched biological pathway and using Cytoscape (Shannon et al., 2003) for graphical representations of pathways.

Target Analysis

We utilized the miRanda (Enright et al., 2003) and the tools for predicting differentially expressed miRNA target on circRNA, lncRNA, and mRNA.

qRT-PCR and Immunohistochemistry

The RNA extracted from the meniscus cells was reverse-transcribed into cDNA using Super-Script™ III Reverse Transcriptase (Invitrogen). Each primer was designed based on the sequence displayed in the NCBI database. We performed quantitative reverse transcriptase PCR (qRT-PCR) using 2× SYBR Green Master Mix (Arraystar, Rockville, MD, United States ) on an Applied Biosystems (Foster City, CA, United States ) ViiA 7 Real-time PCR System. The final reaction system consisted of 1 µl of cDNA, 3.2 µl of double-distilled water, 0.4 µl of forward and backward primers, and 5 µl of 2× SYBR Green PCR Master Mix. Gene expression levels were measured using the 2^−ΔΔCt method. The primer sequences are listed in Supplemental Table S1. In addition, for miRNA validation, total RNA was extracted by miRNeasy Mini Kit (Qiagen, Venlo, Netherlands), and cDNA was synthesized using PrimeScript™ RT Master Mix (Takara, Shiga, Japan). qRT-PCR was performed on a CFX96 system (Bio-Rad, Hercules, CA, United States). GAPDH was used as a housekeeping gene for mRNA, lncRNA, and circRNA, while U6 was applied for miRNA as internal reference genes.

Immunohistochemical analysis was also performed according to previous methods (Sun et al., 2020). For antigen retrieve, sections in 0.1% EDTA were incubated with moderate heat in microwave for 10 min. For staining, sections were treated with 3% normal goat serum for 1 h and incubated with antibodies specific to LCN2 (#26991-1-AP; ProteinTech, Chicago, IL, United States) and RAB27B (#13412-1-AP; ProteinTech).

Statistical Analysis

Statistical analyses were performed using the Statistical Package for the Social Sciences (SPSS), version 25.0 software (SPSS Inc., Chicago, IL, United States). Data are presented as the mean ± SD of the results of at least three independent experiments. Student’s t-test and the Mann–Whitney U test were applied to identify significant differences between groups, where appropriate. Spearman’s rank correlation analysis was used to examine the correlation between two variables (Figure 6D). A p-value < 0.05 was considered statistically significant for all tests. In addition, in order to correct the batch effect, RUVseq package from R language was applied for batch correction. Furthermore, the heatmaps and volcano plots were exported by R language Heatmap package 2, and the scatter plots were exported by ggplot2 package.

Results

Interleukin-1β Might Facilitate Meniscus Degeneration During Osteoarthritis

To test if IL-1β possesses the effect of meniscus degeneration, we treated menisci with IL-1β (5 ng/ml) for 48 h. As a result, meniscus markers like COL1A1, COL2A1, COL3A1, COL6A1, and ACAN were significantly downregulated after inflammatory stimulation, while inflammatory markers like MMP1, MMP3, and ADAMTS5 were upregulated (Figure 1A). Therefore, we suggest that IL-1β might acquire degenerative effect on meniscus, which is similar with chondrocyte during OA.

FIGURE 1

FIGURE 1. Differential expression profile of messenger RNA (mRNA) between degenerative menisci with and without IL-1β stimulation. (A) The expression pattern of meniscus marker genes and inflammatory marker genes in meniscus cells treated with IL-1β (5 ng/ml) as determined by qRT-PCR analysis. GAPDH was used as the internal reference gene for qRT-PCR relative expression. Error bars reveal the standard deviation or the standard error of the data. Student’s t test and Mann–Whitney U test were used to identify significant differences between groups, where appropriate. *p < 0.05, **p < 0.01, ***p < 0.001. (B) Hierarchical clustering illustrates distinguished expression difference of mRNA between the two groups and homogeneity between groups. (C) Volcano plot of differentially expressed mRNAs. (D) Scatter plot of differentially expressed mRNAs. (E) The 20 most enriched Gene Ontology (GO) terms for differentially expressed mRNA in degenerative menisci treated with IL-1β. (F) The 20 most enriched pathway terms for differentially expressed mRNA in degenerative menisci treated with IL-1β. (G) Relative expression levels of selected mRNAs in negative control versus IL-1β-treated osteoarthritis (OA) meniscus. GAPDH was used as the internal reference gene for qRT-PCR relative expression. Error bars reveal the standard deviation or the standard error of the data. Student’s t-test and Mann–Whitney U test were used to identify significant differences between groups, where appropriate. *p < 0.05, **p < 0.01, ***p < 0.001.

Differential Messenger RNA Expression Profile

A total of 14,800 mRNAs were identified in OA meniscus samples. The hierarchical clustering heatmap, volcano plots, and scatter plots revealed the distinguishable gene expression mapping of each sample (Figures 1B–D). After IL-1β stimulation, 145 mRNAs were significantly downregulated (${\log }_2\mathbf{FC}$ < 1, FDR < 0.05), and 230 mRNAs were significantly upregulated (${\log }_2\mathbf{FC}$ > 1, FDR < 0.05) compared with those in degenerative meniscus without IL-1β treatment. Among these, aggrecan (ACAN) (${\log }_2\mathbf{FC}$ = −2.348, FDR = 0) was markedly downregulated, and a disintegrin metallopeptidase with thrombospondin type 1 motif, 5 (ADAMTS5) (${\log }_2\mathbf{FC}$ = 1.093, FDR = 0.011), cholesterol 25-hydroxylase (CH25H) (${\log }_2\mathbf{FC}$ = 27.594, FDR = 0), cytochrome P450, family 7, subfamily B, polypeptide 1 (CYP7B1) (${\log }_2\text{FC}$ = 12.014, FDR = 0), and matrix metalloproteinase 3 (MMP3) were significantly upregulated (${\log }_2\mathbf{FC}$ = 4.917, FDR = 0.030). As both of them were largely studied in OA cartilage, we further validated the sequencing results using qRT-PCR, and the expression trend was concurrent with the sequencing results (Figure 1G). GO and KEGG pathway analyses were performed to uncover the related functions and signaling pathways of the differentially expressed genes (DEGs). The top 20 enriched GO terms and pathways are listed in Figures 1E,F. DEGs were significantly enriched for inflammatory response (FDR = 5.937E−21) and chemotaxis (FDR = 7.175E−14). Inflammatory signaling pathways such as cytokine–cytokine receptor interactions (FDR = 2.129E−14), TNF (FDR = 2.354E−15), and NOD-like receptor signaling pathways (FDR = 3.248E−15) were remarkably enriched with DEGs upon IL-1β treatment. Interestingly, rheumatic arthritis pathway enrichment was also observed.

Differentially Expressed microRNA Profile and Its Target Gene Prediction

MiRNA expression was evaluated in OA menisci with or without IL-1β treatment. In total, 1,145 miRNAs were examined, and only 15 differentially expressed microRNA (DEMs) were identified in the hierarchical clustering heatmap (${\log }_2\mathbf{FC}$ <1 or >1, FDR < 0.05) (Figure 2A). The most upregulated gene was hsa-miR-147b-5p (${\log }_2\mathbf{FC}$ = 3.929, FDR = 3.114E−09). Intriguingly, only one miRNA, hsa-miR-3065-5p, was specifically downregulated by IL-1β stimulation (${\log }_2\mathbf{FC}$ = −1.038, FDR = 0.006) (Table 1). qRT-PCR confirmed several upregulated miRNAs expressed in degenerative menisci with or without IL-1β treatment (Figure 2B).

FIGURE 2

FIGURE 2. Differential expression profile of microRNA (miRNA) between degenerative menisci with and without IL-1β stimulation. (A) Hierarchical clustering illustrates distinguished expression difference of miRNA between the two groups and homogeneity between groups. (B) Relative expression level of selected mRNAs in negative control versus IL-1β-treated osteoarthritis (OA) menisci. U6 was used as the internal reference gene for qRT-PCR relative expression. Error bars reveal the standard deviation or the standard error of the data. The statistical methods are described above. *p < 0.05, **p < 0.01, ***p < 0.001.

TABLE 1

TABLE 1. Fourteen differentially expressed miRNAs (OA menisci with vs. without IL-1β treatment;${\log }_2\mathbf{FC}$ < 1, FDR < 0.05).

Expression Profile of Long Noncoding RNAs and Long Noncoding RNA–MicroRNA–Messenger RNA Network Prediction

In total, 5,997 lncRNAs were identified in this study, with eight significantly downregulated lncRNAs (${\log }_2\text{FC}$ < 1, FDR < 0.05) and 48 upregulated lncRNAs (${\log }_2\text{FC}$ > 1, FDR < 0.05) after IL-1β stimulation. LncRNA LOC105379771 (${\log }_2\mathbf{FC}$ = 5.482, FDR = 8.689E−05) was the most upregulated lncRNA with nearly a 45-fold change, whereas lncRNA DNM1P9 was the most downregulated (${\log }_2\mathbf{FC}$ = −5.002, FDR = 0.0133). Notably, the upregulated lncRNAs were evidently more than the downregulated ones. Hierarchical clustering analysis, volcano plots, and scatter plots revealed distinguishable lncRNA expression profiles between control groups and IL-1β treatment groups (Figures 3A–C), and qRT-PCR validation of four predicted upregulated and three predicted downregulated lncRNAs further confirmed the authenticity (Figure 3D).

FIGURE 3

FIGURE 3. Differential expression profile of long noncoding RNA (lncRNA) and lncRNA LOC107986251 ceRNA network prediction. (A) Hierarchical clustering illustrates distinguished expression difference of lncRNA between the two groups and homogeneity between groups. (B) Volcano plots of differentially expressed lncRNAs. (C) Scatter plots of differentially expressed lncRNAs. (D) Relative expression level of selected lncRNAs in negative control versus IL-1β-treated osteoarthritis (OA) menisci. GAPDH was used as the internal reference gene for qRT-PCR relative expression. Error bars reveal the standard deviation or the standard error of the data. The statistical methods are described above. *p < 0.05, **p < 0.01, ***p < 0.001. (E) LncRNA LOC107986251 network consists of one lncRNA, eight microRNAs (miRNAs), and 97 mRNAs (RNAhybrid_Energy < −25). The red diamond represents downregulated lncRNA LOC107986251. The orange arrows represent upregulated miRNAs. The purple circles represent suppressed mRNAs. (F) Venn diagram of the predicted lncRNA LOC107986251 ceRNA networks by miRanda and RNAhybrid algorithms. (G) qRT-PCR validation of LOC107986251, hsa-miR-212-5p, and SESN3 ceRNA regulation pattern upon IL-1β stimulation in degenerative menisci. GAPDH was used as the internal reference gene for qRT-PCR relative expression. Error bars reveal the standard deviation or the standard error of the data. The statistical methods are described above. *p < 0.05, **p < 0.01, ***p < 0.001.

Additionally, to determine the potential lncRNA regulation mechanism in the menisci during OA, we performed lncRNA–miRNA–mRNA network analysis using the RNAhybrid algorithm (RNAhybrid_Energy < −25). We identified 1,077 different lncRNA–miRNA–mRNA networks regulated by six different lncRNAs, including LOC100130417, LOC101232810, LOC107986251, LOC200772, and DNM1P9. LncRNA LOC107986251 (${\log }_2\mathbf{FC}$ = −1.767, FDR = 1.303E−12) was predicted to have 252 RNA interaction networks in IL-1β-treated degenerative menisci, which possessed the highest amount of RNA interaction networks in the six lncRNAs described earlier (Figure 3E). To discover the most specific and reliable co-expressed RNA (ceRNA) regulatory pathways of lncRNAs, we overlapped the miRanda and RNAhybrid algorithm results (miRanda_Score ≥150, miRanda_Energy < −20, and RNAhybrid_Energy < −25). Consequently, we screened out six ceRNA networks concerning lncRNA LOC107986251 (Figure 3F). A total of 36 different biological processes were identified by GO analysis, and the most enriched were related to regulation of response to oxygen species (FDR = 0.0217) and amyloid precursor protein catabolic processes (FDR = 0.0217) (Supplemental Figure S2A). Only three pathways were confirmed to be enriched in the predicted network during KEGG pathway analysis, in which steroid synthesis (FDR = 0.0177) was the most enriched (Supplemental Figure S2B). Among these, the lncRNA LOC107986251-miR-212-5p-SESN3 network was selected for further qRT-PCR validation in OA menisci with IL-1β treatment, as the downregulation of Sestrin3 (SESN3) in OA cartilage has been described as one of the causes of deficiency in cellular homeostasis, thereby leading to OA. Consequently, validation results were consistent with overlapping prediction (Figure 3G).

Differential Circular RNA Expression Profile and Circular RNA–MicroRNA–Messenger RNA Network Prediction

A total of 13,715 circRNAs were analyzed concerning differentially expressed circRNA (DECs). The heatmap, volcano plots, and scatter plot results illustrated the distinct circRNA variation between degenerative menisci with and without IL-1β cultivation (Figures 4A–C). A total of 55 circRNAs were significantly upregulated, and 34 circRNAs were significantly downregulated in the IL-1β group compared with those in the no IL-1β group. Further, 73 circRNAs had already been identified in the CircBase database, including 46 upregulated circRNAs and 27 downregulated circRNAs; and qRT-PCR confirmed several circRNA expression patterns (Figure 4F). Among these, hsa_circ_0094044 was the most upregulated (${\log }_2\mathbf{FC}$ = 5.926, FDR = 5.288E−07), whereas hsa_circ_0000277 expression was the most evidently suppressed (${\log }_2\mathbf{FC}$ = −4.716, FDR = 9.706E−05). Additionally, 17 circRNAs were not marked in the CircBase database, suggesting several novel circRNAs for further investigation. GO analysis indicated the top 20 highly enriched GO terms and proposed that the parental gene of DECs were largely enriched for cAMP catabolic process (FDR = 0.0900), regulation of nucleic acid-templated transcription (FDR = 0.0900), negative regulation of phosphatase activity (FDR = 0.1159), and Rab GTPase activity (FDR = 0.1158) (Figure 4D). Pathway analysis also revealed the 20 most enriched pathways (Figure 4E). The downregulated transcripts were notably enriched in the case of morphine addiction. For the upregulated transcripts, lysine degradation (FDR = 0.2918) was remarkably enriched. In addition, purine metabolism (FDR = 0.2918), which is highly associated with the pathophysiology of urarthritis, was also prominently enriched with DECs upon IL-1β stimulation.

FIGURE 4

FIGURE 4. Differential expression profile of circular RNA (circRNA) and potential ceRNA prediction. (A) Hierarchical clustering illustrates distinguished expression difference of circRNA between the two groups and homogeneity between groups. (B) Volcano plots of differentially expressed circRNAs. (C) Scatter plots of differentially expressed circRNAs. (D) The 20 most enriched Gene Ontology (GO) terms for the parental genes of differentially expressed circRNA in degenerative menisci treated with IL-1β. (E) The 20 most enriched pathway terms for the parental genes of differentially expressed mRNA in degenerative menisci treated with IL-1β. (F) Relative expression levels of selected circRNAs in negative control versus IL-1β-treated osteoarthritis (OA) menisci. GAPDH was used as the internal reference gene for qRT-PCR relative expression. Error bars reveal the standard deviation or the standard error of the data. The statistical methods are described above. *p < 0.05, **p < 0.01, ***p < 0.001. (G) Hsa_circ_0018069 ceRNA network consists of one circRNA, seven miRNAs, and 97 mRNAs (RNAhybrid_Energy < −25). The red diamond represents downregulated lncRNA LOC107986251. The purple arrows represent upregulated miRNAs. The orange circles represent suppressed mRNAs. (H) Venn diagram of the predicted hsa_circ_0018069 ceRNA networks by miRanda and RNAhybrid algorithms. (I) qRT-PCR validation of hsa_circ_0018069, hsa-miR-147-3p, and TJP2 ceRNA regulation patterns upon IL-1β stimulation in degenerative menisci. GAPDH was used as the internal reference gene for qRT-PCR relative expression. Error bars reveal the standard deviation or the standard error of the data. The statistical methods are described above. *p < 0.05, **p < 0.01, ***p < 0.001.

By merely applying the RNAhybrid algorithm (RNAhybrid_Energy < −25), 1,024 circRNA–miRNA–mRNA networks were predicted to be involved in IL-1β-stimulated degenerative menisci concerning 13 DECs; and hsa_circ_0018069 (${\log }_2\mathbf{FC}$ = −3.030, FDR = 0.0135) was established to be involved in the regulation of 246 ceRNA networks (Figure 4G), which possesses the highest amount of ceRNA relation networks. Only one ceRNA regulatory pathway, hsa_circ_0018069-miR-147b-3p-TJP2 network, was screened out by overlapping the miRanda and RNAhybrid algorithm results (miRanda_Score ≥150, miRanda_Energy < −20, and RNAhybrid_Energy < −25) (Figure 4H). qRT-PCR also confirmed the expression pattern in degenerative menisci with IL-1β stimulation (Figure 4I). Furthermore, GO and KEGG pathway analyses showed that the ceRNA network was highly enriched in biological pathway like the regulation of membrane permeability (FDR = 0.0019) (Supplemental Figures S3A,B).

qRT-PCR Validations in Normal and Degenerative Menisci and Screening for Potential Diagnostic Messenger RNA Biomarkers in the Menisci During Early-Stage Osteoarthritis

To further confirm the results of whole-transcriptome sequencing, we selected previous qRT-PCR-verified DEMs, differentially expressed lncRNAs (DELs), and DECs to validate their expression patterns between normal and degenerative menisci by qRT-PCR (Figure 5A). Consequently, three out of five miRNAs, five out of six lncRNAs, and four out of six circRNAs were shown to have concurrent expression trends with the sequencing results and qRT-PCR validation between the OA menisci with and without IL-1β stimulation (Figures 5B–D). In addition, the lncRNA LOC107986251-miR-212-5p-SESN3 network and hsa_circ_0018069-miR-147b-3p-TJP2 network were also confirmed (Figures 5E,F). Interestingly, miR-147-3p and miR-212-5p expression were not significantly altered between normal and degenerative menisci, unlike in IL-1β-treated groups. This phenomenon may be largely attributed to stronger and immediate exogenous chemokine inflammatory-stimulated effects on fragile degenerative meniscus tissue. Inflammatory stimulation of endogenous chronic chemokines might not be strong enough for the significant alterations of miR-147-3p, miR-212-5p, and other unchanged ncRNAs.

FIGURE 5

FIGURE 5. qRT-PCR certification on control and degenerative meniscus and the selection of potential osteoarthritis (OA) biomarkers in meniscus. (A) Hematoxylin–eosin staining and arthroscopy of the meniscus from normal and degenerative OA meniscus. (B) Relative expression levels of selected differentially expressed microRNAs (DEMs) in normal versus degenerative meniscus. U6 was used as the internal reference gene for qRT-PCR relative expression. Error bars reveal the standard deviation or the standard error of the data. The statistical methods are described above. *p < 0.05, **p < 0.01, ***p < 0.001. (C) Relative expression levels of selected differentially expressed circRNA (DECs) in normal versus degenerative meniscus. GAPDH was used as the internal reference gene for qRT-PCR relative expression. Error bars reveal the standard deviation or the standard error of the data. The statistical methods are described above. *p < 0.05, **p < 0.01, ***p < 0.001. (D) Relative expression levels of selected differentially expressed lncRNAs (DELs) in normal versus degenerative meniscus. GAPDH was used as the internal reference gene for qRT-PCR relative expression. Error bars reveal the standard deviation or the standard error of the data. The statistical methods are described above. *p < 0.05, **p < 0.01, ***p < 0.001. (E) Relative expression patterns of lncRNA LOC107986251hsa_circ_0018069 in normal versus degenerative meniscus. GAPDH was used as the internal reference gene for qRT-PCR relative expression. Error bars reveal the standard deviation or the standard error of the data. The statistical methods are described above. *p < 0.05, **p < 0.01, ***p < 0.001. (F) Relative expression patterns of hsa_circ_0018069 in normal versus degenerative meniscus. GAPDH was used as the internal reference gene for qRT-PCR relative expression. Error bars reveal the standard deviation or the standard error of the data. The statistical methods are described above. *p < 0.05, **p < 0.01, ***p < 0.001. (G) The scatter plots of differentially expressed mRNAs between normal and degenerative menisci. (H) The Venn diagram of single-cell sequence data of normal meniscus versus degenerative meniscus, whole-transcriptome sequence data of control versus IL-1βstimulated OA degenerative meniscus, and RNA sequence data of normal meniscus versus OA degenerative meniscus. (I) qRT-PCR confirmation of the screening mRNA (LCN2, RAB27B, PRDM1, and SERPINB2). GAPDH was used as the internal reference gene for qRT-PCR relative expression. Error bars reveal the standard deviation or the standard error of the data. The statistical methods are described above. *p < 0.05, **p < 0.01, ***p < 0.001.

Simultaneously, we performed RNA-seq on four healthy control menisci and four degenerative menisci in order to select potential diagnostic biomarkers for early-stage OA (Figure 5G). LCN2, RAB27B, SERPINB2, and PRDM1 were screened out by overlapping previously constructed single-cell sequencing data on normal and OA menisci (Sun et al., 2020), and whole-transcriptome sequence data on IL-stimulated meniscus (Figure 5H).

LCN2 and RAB27B Might Act as Biomarkers in Meniscus for OA Severity Predictors and Early OA Diagnosis

We further examined whether meniscus-specific LCN2 and RAB27B possess the potential of predicting OA severity. qRT-PCR confirmed LCN2 and RAB27B expression patterns, both of which show significant upregulation in OA degenerative menisci, while also time-dependently upregulated in inflammatory chemokine-stimulated menisci (Figures 5I, 6A,B). Interestingly, LCN2 and RAB27B expression showed robust correlation with patients’ OA severity based on OARSI Osteoarthritis Cartilage Histopathology Assessment System (Waldstein et al., 2016; Figure 6C). LCN2 and RAB27B were also examined in spontaneous aging C57BL/6J mouse model to validate if meniscus-specific LCN2 and RAB27B could act as biomarkers for early-stage OA. Both of them were discovered to be significantly upregulated at the age of 26 weeks, which is approximately 40 years old in human lifespan (Figure 6D). This suggests that LCN2 and RAB27B might be potential diagnostic biomarkers in meniscus for OA severity prediction and early-stage OA diagnosis.

FIGURE 6

FIGURE 6. LCN2 and RAB27B might act as biomarkers in the meniscus for early osteoarthritis (OA) diagnosis. (A) Relative expression levels of LCN2, RAB27B, PRDM1, and SERPINB2 in IL-1β-treated menisci. GAPDH was used as the internal reference gene for qRT-PCR relative expression. Error bars reveal the standard deviation or the standard error of the data. The statistical methods are described above. *p < 0.05, **p < 0.01, ***p < 0.001. (B) Relative expression levels of LCN2 and RAB27B in menisci treated with time-dependent IL-1β stimulation (0, 12, 24, 48, 72, and 96 h). GAPDH was used as the internal reference gene for qRT-PCR relative expression. Error bars reveal the standard deviation or the standard error of the data. (C) Safranin O/Fast Green staining of patients’ knee cartilage and immunohistochemical (IHC) staining of patients’ menisci with antibody against LCN2 and RAB27B. (D) Correlation coefficient between LCN2 and RAB27B expression quantified by IHC positive cell percentage and OARSI score in patients’ samples (n = 5). (E) Safranin O and IHC staining of LCN2 and RAB27B in mice knee at 1, 2, 4, 8, and 26 weeks and quantification of positive cells (n = 3). Error bars reveal the standard deviation or the standard error of the data. Scale bar, 50 μm.

Discussion

Whole-transcriptome sequencing is a novel bioinformatics analysis method to test the differential expression levels of mRNA, miRNA, lncRNA, and circRNA between normal and pathological tissues. This technique has already been widely applied in the field of oncology (Zheng et al., 2016). A recent study illustrated the comprehensive transcriptome map of normal and OA cartilage and identified four DELs and six DEGs targeted by lncRNAs during OA (Li et al., 2019). Potential OA-associated genes, pathways, competing endogenous RNA networks, and co-expression networks in knee cartilage were further identified in lately studies, thereby offering a better understanding of OA mechanism (Chen and Chen, 2020; Qi et al., 2020). However, a comprehensive analysis of the expression patterns of mRNA, miRNA, lncRNA, and circRNA in OA meniscus, another important knee joint anatomic structure, remains unknown. A previous study had already described that IL-1β stimulation on chondrocytes could act as an in vitro model for OA (Kapoor et al., 2011). Simultaneously, IL-1β performed similar effects on menisci in our study. Therefore, we systematically analyzed the expression profile in degenerative menisci obtained from patients with last-stage OA with or without IL-1β treatment. As a result, we identified 14,800 genes, 1,145 miRNAs, 5,997 lncRNAs, and 13,715 circRNAs. Among these, 375 mRNAs, 15 miRNAs, 56 lncRNAs, and 56 circRNAs were significantly modified subsequent to IL-1β treatment. Following principal component analysis (PCA), we have discovered that sample OA006_NC exhibited high heterogeneity as compared with OA004_NC and OA008_NC (Supplemental Figure S1). This phenomenon might contribute to slight influence on the following sequence results, and we will discuss it in our limitations.

A total of 375 DEGs were examined, and upregulated genes were remarkably more pronounced than downregulated genes. With this, our study confirmed several DEGs that were previously discussed in previous research on OA cartilage, including MMP3 (Shi et al., 2016), superoxide dismutase 2 (SOD2) (Fu et al., 2016), ADAMTS5 (Mokuda et al., 2019), CH25H, cytochrome P450, family 7, subfamily B, polypeptide 1 (CYP7B1) (Choi et al., 2019), and bone morphogenetic protein 2 (BMP2) (Blaney Davidson et al., 2015). Nonetheless, several genes that were found to be differentially expressed in degenerative menisci, such as COL1A1 and COL10A1 (Brophy et al., 2017), were not significantly altered in our study. The lack of sample abundance might contribute to this phenomenon. In terms of GO and KEGG pathway analyses, most enriched genes were highly associated with biological processes implicated in inflammation, such as inflammatory response, chemokine-mediated signaling pathways, chemotaxis, and response to lipopolysaccharide, potentially contributing to meniscus inflammation during the degenerative process. Based on these data, it is possible that IL-1β might contribute to the initiation of general chronic knee joint inflammation within menisci.

The attempt to test the DEMs permitted the discovery of the possible co-expression RNA (ceRNA) regulation networks of lncRNAs and circRNAs. However, we only identified 15 DEMs through sequencing, possibly because of batch effect. In order to screen more DEMs, we performed batch-correction methods to eliminate the effect as much as possible. Consequently, we only screened significantly upregulated miRNAs. As Brophy et al. (Brophy et al., 2018) also predicted relatively low DEMs in the menisci dissected from TKA patients compared with those in arthroscopic partial meniscectomy (APM)-derived menisci, it is possible that only a few DEMs can be detected in degenerative menisci. Interestingly, miR-146-5p was specifically upregulated in OA006_IL-1β (46-fold changes). The differences between the sequences might contribute to meniscus sample heterogeneity between patients as we discussed before, and the inflammatory cytokine treatment might act diversely between different primary meniscus cells. However, after qRT-PCR validation, miR-146-5p was upregulated in all other three samples, suggesting that miR-146-5p is actually upregulated upon IL-1β stimulation. Therefore, we believe that a meniscus database for OA patients needs to be constructed in the future in order to cut down errors brought by sample heterogeneity.

LncRNAs over 200 nucleotides in length are also known to be derived from mammalian genomes and have been studied as a decoy for miRNA to combine with and inhibit expression (Ponting et al., 2009; Wang and Chang, 2011). For instance, Wang et al. (2019) demonstrated that lncRNA FOXD2-AS1 increased the expression levels of TLR4 by sponging with miR-27a-3p, thereby inducing chondrocyte proliferation. On the other hand, knockdown of lncRNA-like lncRNA MF12-AS1 leads to miR-130a-3p upregulation and therefore interferes with the expression of TCF4, which results in increased chondrocyte viability and inhibition of apoptosis, inflammatory response, and extracellular matrix (ECM) degeneration in OA (Luo et al., 2020). All these studies suggest that the sponging function of lncRNA is an important mechanism within OA cartilage. In our present work, we screened out 56 DELs in IL-1β-treated degenerative menisci versus non-IL-1β-treated degenerative menisci. A previous study identified 10 DEL results using TKA to obtain degenerative menisci versus APM to garner a traumatic meniscus (Brophy et al., 2018). LncRNA expression differences might possibly be based on the divergence of OA patients or the conspicuous inflammatory effect of IL-1β. Based on our DEL results, we performed lncRNA–miRNA–mRNA network prediction by applying the RNAhybrid algorithm, and lncRNA LOC107986251 possessed the greatest amount of ceRNA networks in degenerative menisci with IL-1β treatment. Furthermore, we overlapped miRanda and RNAhybrid results to screen out the most specific lncRNA regulatory network. Six lncRNA–miRNA–mRNA ceRNA networks are potentially regulated in the pathogenesis of meniscus OA. Among these, SESN3, which was previously investigated for supporting chondrocyte homeostasis and is suppressed in OA cartilage (Shen et al., 2017), was also downregulated by the modulation of the LOC107986251-hsa-miR-212-5p-SESN3 network in OA-induced degenerative menisci. The qRT-PCR validation supported this result. Therefore, the downregulation of lncRNA LOC107986251 might induce miR-212-5p expression and inhibit SESN3 expression, leading to the meniscus and cartilage degenerative process, suggesting a potential crosslink between menisci and cartilage during OA. Nonetheless, deeper mechanistic validation is needed to confirm this hypothesis.

CircRNAs are novel regulatory RNAs that have been recently investigated in OA chondrocytes. Recently, CircSERPINE2-miR-1271-5P-E26 transformation-specific-related gene axis was found to have a strong protective effect against OA by inhibiting chondrocyte apoptosis and ECM degeneration (Shen et al., 2019). This finding suggests that the miRNA-sponging function of circRNA is a potential source of OA development. However, few studies on circRNA have been employed for OA-induced degenerative menisci. Hence, as whole-transcriptome sequencing can simultaneously detect circRNA expression as well, we identified novel DECs in degenerative menisci with or without IL-1β stimulation, including 56 significantly upregulated and 34 significantly downregulated circRNAs. Similarly, we predicted the circRNA–miRNA–mRNA network using the same method as the prediction of the lncRNA ceRNA network, discovering that hsa_circ_0018069 possessed the highest number of networks. Specifically, only one ceRNA regulation network was identified by overlapping the miRanda and RNAhybrid outputs. The downregulation of hsa_circ_0018069 retarded its sponging effect with hsa-miR-147b-3p so that the expression of TJP2 was dysregulated during the OA process in the menisci. Conversely, the hsa_circ_0018069-miR-147b-3p-TJP2 axis might also serve as protection against meniscus degeneration in OA, like circSERPINE2 on the cartilage. It has been reported that hsa_circ_0018069 expression is inhibited in bladder cancer tissues and may serve as a clinical biomarker for early bladder cancer. TJP2 has been studied previously. However, none of these network components have been further evaluated in menisci or cartilage during OA, which might suggest a possible novel regulatory pathway in meniscus degeneration. Our qRT-PCR validation confirmed the predicted expression pattern of hsa_circ_0018069-miR-147B-3p-TJP2 in the menisci with IL-1β stimulation, yet this axis still requires further verification in vitro and in vivo.

Menisci have been largely reported to have an important role in OA progress, and destabilization of the medial meniscus (DMM) model is a common OA model for mice (Berthiaume et al., 2005; Hunter et al., 2006; Murphy et al., 2019). However, whether meniscus degeneration or meniscus-specific biomarkers forecast the onset of OA or the severity of OA remains unknown. Hence, aside from ncRNAs, we screened mRNAs by overlapping the three constructed meniscus databases. LCN2 and RAB27B exhibited significant upregulation during OA in menisci. LCN2, also known as neutrophil gelatinase-associated lipocalin, has recently been identified as a pro-inflammatory adipokine in OA chondrocytes. However, LCN2 overexpression via adenovirus injection into the murine joint did not trigger OA pathology, and global LCN2 knockout mice showed no restoration of cartilage in DMM-induced mice (Choi and Chun, 2017). This might imply that early-stage OA was triggered via LCN2 activation in the menisci but not in the chondrocytes. On the other hand, no studies have been conducted concerning RAB27B in OA; however, a recent investigation revealed that RAB27B acts as a downstream mediator of HIF-2α to regulate the formation of the vascular network (Bhurke et al., 2020). Intriguingly, RAB27B had also been predicted to be highly expressed in DegP in a previous study (Sun et al., 2020). Since meniscus degeneration after trauma or tear is of high relevance to the avascular traits in the white zone of the menisci during OA, it is reasonable to believe that RAB27B might contribute to increasing avascular traits during meniscus degeneration in OA. Following immunohistochemical studies on menisci derived from human suggested that these two meniscus-specific biomarkers correlated with OA severity. In vivo study showed meniscus-specific LCN2 and RAB27B remarkably upregulated at the age of 26 weeks (6 months) in mice and specifically distributed in the internal zone of menisci. Strangely, mice at the age of 52 weeks (1 year) did not show highly positive LCN2 and RAB27B in menisci. We hypothesize that LCN2 and RAB27B might act as a warning and protective signal for OA in murine knee. As the aging and OA develops, their expression begins to fade away. In OA patients, on the other hand, meniscus-specific LCN2 and RAB27B remain consistently expressed and play the part of OA severity prediction. Since menisci or cartilage from early-stage OA patients were usually not able to be obtained, the relation between LCN2 and RAB27B and the period of OA prediction in human remain blurry and require further analysis. Anyway, both of these results are promising for the study of the mechanism underlying meniscus degeneration during OA.

The main strength of this study is to use the advanced high throughout sequence methods—whole-transcriptome sequence to predict the potential mRNA and noncoding RNA, which is more comprehensive than mere RNA sequence. Furthermore, based on the whole-transcriptome sequence data, we overlapped miRanda and RNAhybrid predicting algorithm, and we were able to predict two specific RNA regulatory axis—lncRNA LOC107986251-miR-212-5p-SESN3 and hsa_circ_0018069-miR-147b-3p-TJP2—which could be a novel target for the early treatment of degenerative menisci. More importantly, by combining different databases, we were also able to discover two highly specific markers, LCN2 and RAB27B, which are also highly specific since these two biomarkers were both significantly altered in three different databases of degenerative meniscus.

Although several novel findings were proposed in the OA-induced degenerative meniscus, this study has several limitations. First of all, IL-1β diluent was not used as an exact positive control, while we applied refreshed medium instead. Moreover, following PCA, we have discovered that sample OA006_NC exhibited heterogeneity compared with OA004_NC and OA008_NC (Supplemental Figure S1). This phenomenon might contribute to slight influence on the following sequence results, and we will discuss it in our limitations. Hence, a larger database of degenerative menisci from OA patients and even normal menisci should be built in order to offer a global understanding of distinct genes and ncRNA expression during meniscus degeneration, so that further investigation of meaningful clinical biomarkers for OA patients can be efficiently performed. It could also cut down some examination errors brought by sample heterogeneity as we mentioned above. Another limitation is the highly rigorous selection for lncRNA and circRNA target prediction by overlapping miRanda, RNAhybrid algorithm, and miRNA sequencing, which might contribute to relatively less ceRNA network results. Nevertheless, it also helps us to identify highly specific ceRNA regulatory pathways during meniscus degeneration during OA. In addition, we performed simple validation on the differential expression of each ncRNA and mRNA using qRT-PCR. To further confirm their specific mechanism and function in the degenerative process of OA menisci, more in-depth research into significantly upregulated and downregulated ncRNAs should be performed.

In summary, this study illustrated a transcriptome profile of OA menisci by a whole-transcriptome sequencing method and specifically identified two highly specific ceRNA networks regulated by lncRNA LOC107986251 and hsa_circ_0018069, which possibly play an important role during the meniscal degeneration process, and two potential mRNA biomarkers, LCN2 and RAB27B, in the meniscus for future OA diagnosis. All these bioinformatics results could be of value to researchers seeking to understand the underlying mechanism of meniscus degeneration in OA, hence exploiting new diagnostic biomarkers for early-stage OA clinical lab tests and potential therapeutic targets for treating OA. Furthermore, the relationship between degenerative knee menisci and cartilage during the OA process could also be explored based on the present study. However, additional efforts are merited for DEGs, DEMs, DELs, DECs, and ceRNA networks to achieve the aforementioned goals.

Data Availability Statement

The data presented in the study are deposited in the GEO repository, accession number GSE185064, accession number GSE171652.

Ethics Statement

The animal study was reviewed and approved by the Ethical Committee of The First Affiliated Hospital of Sun Yat-sen University. Written informed consent was obtained from the individual(s), and minor(s)’ legal guardian/next of kin, for the publication of any potentially identifiable images or data included in this article.

Author Contributions

ZJ and XD designed and performed the experiments, and they also wrote the article. XW and HS conceived the research and collected the samples. HL, SH, and QH collected and analyzed the data. WL and ZZ conceived and supervised the study and also contributed to writing the article.

Funding

This study was funded by the National Natural Science Foundation of China (81972049, 82172467), the First Affiliated Hospital of Sun Yat-sen University Ke Ling Funding program for Novel and Distinguished talents (R07005), the Guangdong Natural Science Funds for Distinguished Young Scholar of China (2021B1515020008), the Science and Technology Project of Guangzhou City, China, Number: 201710010164, the Postdoctoral Science Foundation of China (2020M683085), the Guangdong Provincial Natural Science Foundation of China (2020A1515110943, 2021A1515010454), Sun Yat-sen University’s College Basic Research Service Fee Project-Young Teacher Cultivation Project (19ykpy63), and the Natural Science Foundation of Guangdong Province (grant number: 2016A030310156).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors thank Professor Xuerong Li (Department of the Zhongshan School of Medicine, Sun Yat-sen University) for technical assistance.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.754421/full#supplementary-material

References

Anders, S., Pyl, P. T., and Huber, W. (2015). HTSeq--a Python Framework to Work with High-Throughput Sequencing Data. Bioinformatics 31, 166–169. doi:10.1093/bioinformatics/btu638

PubMed Abstract | CrossRef Full Text | Google Scholar

Benjamini, Y., Drai, D., Elmer, G., Kafkafi, N., and Golani, I. (2001). Controlling the False Discovery Rate in Behavior Genetics Research. Behav. Brain Res. 125, 279–284. doi:10.1016/s0166-4328(01)00297-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Berthiaume, M.-J., Raynauld, J. P., Martel-Pelletier, J., Labonte, F., Beaudoin, G., Bloch, D. A., et al. (2005). Meniscal Tear and Extrusion Are Strongly Associated with Progression of Symptomatic Knee Osteoarthritis as Assessed by Quantitative Magnetic Resonance Imaging. Ann. Rheum. Dis. 64, 556–563. doi:10.1136/ard.2004.023796

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhurke, A., Kannan, A., Neff, A., Ma, Q., Laws, M. J., Taylor, R. N., et al. (2020). A Hypoxia-Induced Rab Pathway Regulates Embryo Implantation by Controlled Trafficking of Secretory Granules. Proc. Natl. Acad. Sci. USA 117, 14532–14542. doi:10.1073/pnas.2000810117

PubMed Abstract | CrossRef Full Text | Google Scholar

Blaney Davidson, E. N., Vitters, E. L., Bennink, M. B., van Lent, P. L. E. M., van Caam, A. P. M., Blom, A. B., et al. (2015). Inducible Chondrocyte-specific Overexpression of BMP2 in Young Mice Results in Severe Aggravation of Osteophyte Formation in Experimental OA without Altering Cartilage Damage. Ann. Rheum. Dis. 74, 1257–1264. doi:10.1136/annrheumdis-2013-204528

PubMed Abstract | CrossRef Full Text | Google Scholar

Brophy, R. H., Sandell, L. J., and Rai, M. F. (2017). Traumatic and Degenerative Meniscus Tears Have Different Gene Expression Signatures. Am. J. Sports Med. 45, 114–120. doi:10.1177/0363546516664889

PubMed Abstract | CrossRef Full Text | Google Scholar

Brophy, R. H., Zhang, B., Cai, L., Wright, R. W., Sandell, L. J., and Rai, M. F. (2018). Transcriptome Comparison of Meniscus from Patients with and without Osteoarthritis. Osteoarthritis Cartilage 26, 422–432. doi:10.1016/j.joca.2017.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, H., and Chen, L. (2020). An Integrated Analysis of the Competing Endogenous RNA Network and Co-expression Network Revealed Seven Hub Long Non-coding RNAs in Osteoarthritis. Bone Jt. Res. 9, 90–98. doi:10.1302/2046-3758.93.bjr-2019-0140.r2

PubMed Abstract | CrossRef Full Text | Google Scholar

Choi, W.-S., and Chun, J.-S. (2017). Upregulation of Lipocalin-2 (LCN2) in Osteoarthritic Cartilage Is Not Necessary for Cartilage Destruction in Mice. Osteoarthritis Cartilage 25, 401–405. doi:10.1016/j.joca.2016.07.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Choi, W.-S., Lee, G., Song, W.-H., Koh, J.-T., Yang, J., Kwak, J.-S., et al. (2019). The CH25H-CYP7B1-RORα axis of Cholesterol Metabolism Regulates Osteoarthritis. Nature 566, 254–258. doi:10.1038/s41586-019-0920-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Cook, A., Cook, J., and Stoker, A. (2018). Metabolic Responses of Meniscus to IL-1β. J. Knee Surg. 31, 834–840. doi:10.1055/s-0037-1615821

PubMed Abstract | CrossRef Full Text | Google Scholar

Englund, M., Roemer, F. W., Hayashi, D., Crema, M. D., and Guermazi, A. (2012). Meniscus Pathology, Osteoarthritis and the Treatment Controversy. Nat. Rev. Rheumatol. 8, 412–419. doi:10.1038/nrrheum.2012.69

PubMed Abstract | CrossRef Full Text | Google Scholar

Enright, A. J., John, B., Gaul, U., Tuschl, T., Sander, C., and Marks, D. S. (2003). MicroRNA Targets in Drosophila. Genome Biol. 5, R1. doi:10.1186/gb-2003-5-1-r1

PubMed Abstract | CrossRef Full Text | Google Scholar

Filardo, G., Petretta, M., Cavallo, C., Roseti, L., Durante, S., Albisinni, U., et al. (2019). Patient-specific Meniscus Prototype Based on 3D Bioprinting of Human Cell-Laden Scaffold. Bone Jt. Res. 8, 101–106. doi:10.1302/2046-3758.82.bjr-2018-0134.r1

PubMed Abstract | CrossRef Full Text | Google Scholar

Fu, Y., Kinter, M., Hudson, J., Humphries, K. M., Lane, R. S., White, J. R., et al. (2016). Aging Promotes Sirtuin 3-Dependent Cartilage Superoxide Dismutase 2 Acetylation and Osteoarthritis. Arthritis Rheumatol. 68, 1887–1898. doi:10.1002/art.39618

PubMed Abstract | CrossRef Full Text | Google Scholar

Hunter, D. J., Zhang, Y. Q., Niu, J. B., Tu, X., Amin, S., Clancy, M., et al. (2006). The Association of Meniscal Pathologic Changes with Cartilage Loss in Symptomatic Knee Osteoarthritis. Arthritis Rheum. 54, 795–801. doi:10.1002/art.21724

PubMed Abstract | CrossRef Full Text | Google Scholar

Ji, Q., Zheng, Y., Zhang, G., Hu, Y., Fan, X., Hou, Y., et al. (2019). Single-cell RNA-Seq Analysis Reveals the Progression of Human Osteoarthritis. Ann. Rheum. Dis. 78, 100–110. doi:10.1136/annrheumdis-2017-212863

PubMed Abstract | CrossRef Full Text | Google Scholar

Kapoor, M., Martel-Pelletier, J., Lajeunesse, D., Pelletier, J.-P., and Fahmi, H. (2011). Role of Proinflammatory Cytokines in the Pathophysiology of Osteoarthritis. Nat. Rev. Rheumatol. 7, 33–42. doi:10.1038/nrrheum.2010.196

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, D., Langmead, B., and Salzberg, S. L. (2015). HISAT: a Fast Spliced Aligner with Low Memory Requirements. Nat. Methods 12, 357–360. doi:10.1038/nmeth.3317

PubMed Abstract | CrossRef Full Text | Google Scholar

Kwon, H., Brown, W. E., Lee, C. A., Wang, D., Paschos, N., Hu, J. C., et al. (2019). Surgical and Tissue Engineering Strategies for Articular Cartilage and Meniscus Repair. Nat. Rev. Rheumatol. 15, 550–570. doi:10.1038/s41584-019-0255-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Langmead, B., and Salzberg, S. L. (2012). Fast Gapped-Read Alignment with Bowtie 2. Nat. Methods 9, 357–359. doi:10.1038/nmeth.1923

PubMed Abstract | CrossRef Full Text | Google Scholar

Lei, B., Zhou, J., Xuan, X., Tian, Z., Zhang, M., Gao, W., et al. (2019). Circular RNA Expression Profiles of Peripheral Blood Mononuclear Cells in Hepatocellular Carcinoma Patients by Sequence Analysis. Cancer Med. 8, 1423–1433. doi:10.1002/cam4.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Leng, N., Dawson, J. A., Thomson, J. A., Ruotti, V., Rissman, A. I., Smits, B. M. G., et al. (2013). EBSeq: an Empirical Bayes Hierarchical Model for Inference in RNA-Seq Experiments. Bioinformatics 29, 1035–1043. doi:10.1093/bioinformatics/btt087

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, H., Yang, H. H., Sun, Z. G., Tang, H. B., and Min, J. K. (2019). Whole-transcriptome Sequencing of Knee Joint Cartilage from Osteoarthritis Patients. Bone Jt. Res. 8, 290–303. doi:10.1302/2046-3758.87.bjr-2018-0297.r1

CrossRef Full Text | Google Scholar

Luo, X., Wang, J., Wei, X., Wang, S., and Wang, A. (2020). Knockdown of lncRNA MFI2-AS1 Inhibits Lipopolysaccharide-Induced Osteoarthritis Progression by miR-130a-3p/TCF4. Life Sci. 240, 117019. doi:10.1016/j.lfs.2019.117019

PubMed Abstract | CrossRef Full Text | Google Scholar

Makris, E. A., Hadidi, P., and Athanasiou, K. A. (2011). The Knee Meniscus: Structure-Function, Pathophysiology, Current Repair Techniques, and Prospects for Regeneration. Biomaterials 32, 7411–7431. doi:10.1016/j.biomaterials.2011.06.037

PubMed Abstract | CrossRef Full Text | Google Scholar

McNulty, A. L., Rothfusz, N. E., Leddy, H. A., and Guilak, F. (2013). Synovial Fluid Concentrations and Relative Potency of Interleukin-1 Alpha and Beta in Cartilage and Meniscus Degradation. J. Orthop. Res. 31, 1039–1045. doi:10.1002/jor.22334

CrossRef Full Text | Google Scholar

Meng, F., Li, Z., Zhang, Z., Yang, Z., Kang, Y., Zhao, X., et al. (2018). MicroRNA-193b-3p Regulates Chondrogenesis and Chondrocyte Metabolism by Targeting HDAC3. Theranostics 8, 2862–2883. doi:10.7150/thno.23547

PubMed Abstract | CrossRef Full Text | Google Scholar

Mokuda, S., Nakamichi, R., Matsuzaki, T., Ito, Y., Sato, T., Miyata, K., et al. (2019). Wwp2 Maintains Cartilage Homeostasis through Regulation of Adamts5. Nat. Commun. 10, 2429. doi:10.1038/s41467-019-10177-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Murphy, C. A., Garg, A. K., Silva-Correia, J., Reis, R. L., Oliveira, J. M., and Collins, M. N. (2019). The Meniscus in Normal and Osteoarthritic Tissues: Facing the Structure Property Challenges and Current Treatment Trends. Annu. Rev. Biomed. Eng. 21, 495–521. doi:10.1146/annurev-bioeng-060418-052547

PubMed Abstract | CrossRef Full Text | Google Scholar

Ogata, H., Goto, S., Sato, K., Fujibuchi, W., Bono, H., and Kanehisa, M. (1999). KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res. 27, 29–34. doi:10.1093/nar/27.1.29

PubMed Abstract | CrossRef Full Text | Google Scholar

Ponting, C. P., Oliver, P. L., and Reik, W. (2009). Evolution and Functions of Long Noncoding RNAs. Cell 136, 629–641. doi:10.1016/j.cell.2009.02.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Qi, X., Yu, F., Wen, Y., Li, P., Cheng, B., Ma, M., et al. (2020). Integration of Transcriptome-wide Association Study and Messenger RNA Expression Profile to Identify Genes Associated with Osteoarthritis. Bone Jt. Res. 9, 130–138. doi:10.1302/2046-3758.93.bjr-2019-0137.r1

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a Software Environment for Integrated Models of Biomolecular Interaction Networks. Genome Res. 13, 2498–2504. doi:10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, T., Alvarez-Garcia, O., Li, Y., Olmer, M., and Lotz, M. K. (2017). Suppression of Sestrins in Aging and Osteoarthritic Cartilage: Dysfunction of an Important Stress Defense Mechanism. Osteoarthritis Cartilage 25, 287–296. doi:10.1016/j.joca.2016.09.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, S., Wu, Y., Chen, J., Xie, Z., Huang, K., Wang, G., et al. (2019). CircSERPINE2 Protects against Osteoarthritis by Targeting miR-1271 and ETS-Related Gene. Ann. Rheum. Dis. 78, 826–836. doi:10.1136/annrheumdis-2018-214786

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, J., Zhang, C., Yi, Z., and Lan, C. (2016). Explore the Variation of MMP3, JNK, P38 MAPKs, and Autophagy at the Early Stage of Osteoarthritis. IUBMB Life 68, 293–302. doi:10.1002/iub.1482

PubMed Abstract | CrossRef Full Text | Google Scholar

Shimomura, K., Rothrauff, B. B., Hart, D. A., Hamamoto, S., Kobayashi, M., Yoshikawa, H., et al. (2019). Enhanced Repair of Meniscal Hoop Structure Injuries Using an Aligned Electrospun Nanofibrous Scaffold Combined with a Mesenchymal Stem Cell-Derived Tissue Engineered Construct. Biomaterials 192, 346–354. doi:10.1016/j.biomaterials.2018.11.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, H., Wen, X., Li, H., Wu, P., Gu, M., Zhao, X., et al. (2020). Single-cell RNA-Seq Analysis Identifies Meniscus Progenitors and Reveals the Progression of Meniscus Degeneration. Ann. Rheum. Dis. 79, 408–417. doi:10.1136/annrheumdis-2019-215926

PubMed Abstract | CrossRef Full Text | Google Scholar

Waldstein, W., Perino, G., Gilbert, S. L., Maher, S. A., Windhager, R., and Boettner, F. (2016). OARSI Osteoarthritis Cartilage Histopathology Assessment System: A Biomechanical Evaluation in the Human Knee. J. Orthop. Res. 34, 135–140. doi:10.1002/jor.23010

CrossRef Full Text | Google Scholar

Wang, K. C., and Chang, H. Y. (2011). Molecular Mechanisms of Long Noncoding RNAs. Mol. Cel. 43, 904–914. doi:10.1016/j.molcel.2011.08.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Cao, L., Wang, Q., Huang, J., and Xu, S. (2019). LncRNA FOXD2-AS1 Induces Chondrocyte Proliferation through Sponging miR-27a-3p in Osteoarthritis. Artif. Cell Nanomed. Biotechnol. 47, 1241–1247. doi:10.1080/21691401.2019.1596940

PubMed Abstract | CrossRef Full Text | Google Scholar

You, X., Vlatkovic, I., Babic, A., Will, T., Epstein, I., Tushev, G., et al. (2015). Neural Circular RNAs Are Derived from Synaptic Genes and Regulated by Development and Plasticity. Nat. Neurosci. 18, 603–610. doi:10.1038/nn.3975

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, B., Schmoyer, D., Kirov, S., and Snoddy, J. (2004). GOTree Machine (GOTM) a Web-Based Platform for Interpreting Sets of Interesting Genes Using Gene Ontology Hierarchies. BMC Bioinf. 5, 16. doi:10.1186/1471-2105-5-16

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Z. Z., Chen, Y. R., Wang, S. J., Zhao, F., Wang, X. G., Yang, F., et al. (2019). Orchestrated Biomechanical, Structural, and Biochemical Stimuli for Engineering Anisotropic Meniscus. Sci. Transl Med. 11 (487), eaao0750. doi:10.1126/scitranslmed.aao0750

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, Q., Bao, C., Guo, W., Li, S., Chen, J., Chen, B., et al. (2016). Circular RNA Profiling Reveals an Abundant circHIPK3 that Regulates Cell Growth by Sponging Multiple miRNAs. Nat. Commun. 7, 11215. doi:10.1038/ncomms11215

PubMed Abstract | CrossRef Full Text | Google Scholar