Identication And Validation of Pivotal Genes Related To Meniscus Senescence Based On Gene Expression Proling Analysis And In Vivo And In Vitro Models Detection

Background: The compositional change in the meniscus with aging would increase the tissue vulnerability of the meniscus, which would induce meniscus tearing. Here, we investigated the molecular mechanism of age-related meniscus degeneration with gene expression proling analysis, and validate pivotal genes in vivo and in vitro models. Methods: The GSE45233 dataset, including 6 elderly meniscus samples and 6 younger meniscus samples, was downloaded from the Gene Expression Omnibus (GEO) database. To screen the differential expression of mRNAs, identify the miRNAs targeting hub genes, and forecast the potentially toxic drugs, we completed a series of bioinformatics analyses, including functional and pathway enrichment, protein-protein interaction network, hub genes screening, construction of a lncRNA–miRNA–mRNA network, and molecular docking of potential drugs. Furthermore, crucial genes were examined in human senescent menisci, mouse senescent meniscus tissues and mouse meniscus cells stimulated by IL-1β. Results: In total, the most signicant 4 hub genes (RRM2, AURKB, CDK1, and TIMP1), 5 miRNAs (hsa-miR-6810-5p, hsa-miR-4676-5p, hsa-miR-6877-5p, hsa-miR-8085, and hsa-miR-6133) that regulated such 4 hub genes, and potential toxic drugs (Cladribine, Danusertib, Barasertib, Riviciclib, and Dinaciclib) that had a targeting effect on these genes, were nally identied. Moreover, these hub genes were decreased in meniscus cells in vitro and meniscus tissues in vivo, which indicated that hub genes were related to meniscus senescence and could serve as potential biomarkers for age-related meniscus tearing. Conclusions: In short, the integrated analysis of gene expression prole, co-expression network, and models detection identied pivotal genes, which elucidated the possible molecular basis underlying the senescence meniscus and also provided prognosis clues for early-onset age-related meniscus tearing.


Introduction
The meniscus is an important structure of knee joint, playing a vital role in load transmission, shock absorption, joint lubrication and nutrition, proprioception, and stability [1]. The position of the meniscus and its function in load transfer and shock absorption make the meniscus prone to traumatic and degenerative injuries. Therefore, the meniscus tearing is one of the most common intra-articular injuries of knee joint, and also one of the major causes of gonyalgia and lameness [2,3]. Other studies have also suggested that meniscus tearing injury is not only one of the high-risk factors for inducing osteoarthritis (OA), but also further aggravates the process of osteoarthritis, and ultimately affects people's basic daily activities [4,5].
The population-based cohort study found that the age is a major risk factor for meniscus tearing, especially in people over 40 years old [6,7]. Previous studies have shown that aging could cause compositional changes in the meniscus, which includes degeneration of collagens, cells, and proteoglycans [8,9]. Meanwhile, these age-related tissue changes, such as aggregation of advanced glycation end products, thickening of collagen bers, and mucoid degeneration, would increase the tissue vulnerability of the meniscus, which would eventually induce meniscus tearing [10]. Tsujii et al. found that age-related tissue homeostasis changes in the meniscus might be caused by senescence of meniscus cells [8]. However, there are few studies on meniscus aging, and also the speci c molecular mechanism and key genes of meniscus cell senescence are still unclear [11,12].
Bioinformatics is an emerging interdisciplinary eld which combines molecular biology and information technology. It is of great importance to reveal that the differentially expressed genes (DEGs) are involved in molecular mechanisms, speci c pathways, protein-protein interactions (PPI), and associations between disease and genes [13]. In order to study the molecular mechanism of age-related meniscus degeneration and further explore potential biological therapeutic targets, we carried out extensive bioinformatics methods to identify key DEGs and functional pathways in the aged meniscus tissues.
Then, we constructed a lncRNA-miRNA-mRNA network and performed pharmacogenomics analysis to search potential miRNAs and drugs that may act on the screened key DEGs. Moreover, we further con rmed the expressions of these key DEGs in human senescent menisci, mouse aging meniscus tissues, and IL-1β-stimulated mouse meniscus cells.

Data collection
The GSE45233 dataset in Gene Expression Omnibus was downloaded for subsequent analysis.

Identi cation of DEGs in elderly and younger samples
We used the robust multiarray average algorithm to perform background correction and quartile data normalization on the GSE45233 dataset, and then identi ed DEGs associated with age by means of the Limma package in R language (version3.6.3). P value < 0.05 and log 2 fold change (log 2 FC) < − 1.5 or log 2 FC > 1.5 were de ned as the cutoff standard. The volcano plot and heatmap were completed in virtue of Ggplot2 and Pheatmap packages.
KEGG and GO analysis annotate the functions of DEGs DAVID (version 6.8; https://david.ncifcrf.gov/) is an online bioinformatics database that integrates biological data and analysis tools to provide systematic and comprehensive biological annotation information for extensive lists of genes and proteins. The concept of Gene Ontology (GO) is widely used in the eld of bioinformatics, which contains three notions of functional information: biological processes (BP), cellular components (CC), and molecular functions (MF). Currently, Kyoto Encyclopedia of Genes and Genomes (KEGG; https://www.kegg.jp/) is a database that integrates genomic, chemical, and system functional information for signaling pathway analysis. We employed the DAVID online tool to perform GO and KEGG analysis of DEGs.
Construction of the PPI network STRING (Version11.0; https://string-db.org/), a search tool for interacting genes/proteins, was applied for predicting the PPI network and discovering the possible relationships. Cytoscape (version 3.8.0; https://cytoscape.org/), a powerful network building software, was used to establish the PPI network.

Identi cation and analysis of signi cant hub genes
The hub genes were screened through the PPI network utilizing the Cytoscape plugin cytoHubba. Since the singlular algorithm is sometimes biased, we employed a four-fold algorithm to collectively recognize hub genes, and these algorithms were adopted respectively. The algorithms used to recognize hub genes include Edge Percolated Component (EPC), Degree, Maximal Clique Centrality (MCC), and Maximum Neighborhood Component (MNC). Then the common hub genes were obtained through the intersection of Venn diagram websites (http://bioinformatics.psb.ugent.be/webtools/Venn). Metascape, as a forceful gene function analysis and annotation tools, was applied for batch enrichment analysis and network construction of genes and proteins to understand their functions. The common hub genes were further analysed through Metascape analysis.
Construction of lncRNA-miRNA-mRNA network miRWalk (http://mirwalk.umm.uni-heidelberg.de/), as a bioinformatics suit to forecast miRNA-mRNA interactions, could predict miRNAs of the common Hub gene and then obtain the miRNA-mRNA interaction pairs with the score > 0.9.
[16] The DIANA-miRPath (https://www.microrna.gr/miRPathv3) is an online tool dedicated to assessing the regulatory role of miRNAs and forecasting relevant regulation pathways. The DIANA-LncBase (http://carolina.imis.athena-innovation.gr/diana_tools/) was conducted to search the lncRNAs targeting at miRNAs. In the end, a lncRNA-miRNA-mRNA regulatory network was accomplished by means of Cytoscape software.

Pharmacogenomics analysis of hub genes
The Pharmacogenomics and Drug-Gene Interaction database (DGIdb, http://dgidb.org/) work on collecting and categorizing information about the relationship between genetic variation and drug response. In this research, the network database was devoted to forecast drugs that might target at hub genes. Then, we downloaded 2D structures of small molecule drugs from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/) before being prepared for docking. The software PyMol (version 2.4; https://pymol.org/2/) was applied to determine the lowest energy conformations for docking. With the help of AutoDock software (Version 4.2.6; http://autodock.scripps.edu/), we completed the exible molecular docking between the macromolecular proteins and the drugs to determine the relative position and orientation in the process of a ligand binding to a receptor. Finally, the images were captured from software PyMoL.
Human/Mouse meniscus tissues from normal and aging knees Human meniscus tissues were obtained from 12 patients undergoing arthroscopic partial meniscectomy, including 6 patients older than 40 years and 6 patients younger than 40 years. This human research protocol was authorized by the Medical Ethics Committees (approval number: 2019K-K011). All patients have signed informed consent. Animal protocols used in this study were approved by the Animal Welfare week, where standard laboratory chow and water were freely consumed. Then, after the mice were anesthetized by intraperitoneal injection of 1% pentobarbital (60mg/kg), the mice were euthanatized by spinal dislocation and the femurs were collected for subsequent detection. We provided an ARRIVE checklist in the Supplementary Material to reveal that the ARRIVE guidelines were adhered to in this research.

Histological analysis in the meniscus tissues
After xation with 4% paraformaldehyde for 2 days, the entire mouse knee joints and human meniscus were decalci ed, dehydrated and embedded in para n. The sections of mouse knee joints were stained with Safranin-O Fast Green and Toloniumchloride Blue to observe histopathological changes.
For immunohistochemistry analysis, knee joint sections were depara nized, hydrated, antigen retrieval, and blocked in 3% BSA (Servicebio, Wuhan, China) at room temperature for 1 h. Then, the sections incubated with primary antibodies for AURKB (1:200, dilution), CDK1 (1:100, dilution), RRM2 (1:100, dilution), and TIMP1(1:250, dilution, Abcam) overnight at 4°C, followed by incubating with biotinylated secondary antibody (1:100 dilution). Finally, after the sections incubating with an avidin-biotinylated HRP complex solution, these were immersed in diaminobenzidine (DAB) to reveal peroxidase activity. Observations and photographs were taken under the Photo Imaging System (Nikon H550S, Japan). For immunohistochemical-negative controls, the primary antibodies were replaced by non-immunized IgG in immunostaining. Then, we measured the mean optical density (MOD) of 5 random elds in each section using Image Pro-plus (version 6.0) to determine the dyeing intensity. The percentage of positive cells in each eld was the ratio of the number of positive cells to the total number of meniscus cells in the corresponding region.

SA-β-gal staining
The mouse meniscus cells were derived from meniscus tissues of SPF 3-month-old C57BL/6J mice. In brief, the full thickness meniscus tissues were cut into pieces, and digested with serum-free DMEM/antibiotics containing collagenase D (Roche, Basel, Switzerland) for 6h at 37°C [17]. The meniscal cells were cultured on DMEM containing 10% FBS and inoculated in 12-well plates. Then, treatment of the cells was conducted with IL-1β (10ng/ml) for 2d to achieve the induction of cell senescence injury by in ammatory factors [18]. The meniscus cells were xed at room temperature for 15 min, after washing three times with PBS, and stained with β-galactosidase staining solution for 12 h at 37°C without CO 2 environment. After that, these were rinsed again with PBS and were observed via an inverted microscope (Nikon, Japan).

Apoptosis assay
For apoptosis analysis, the meniscus cells treated with IL-1β (10ng/ml) for 2 days. After cells resuspended with 100 µL1x binding buffer, 5 µL Annexin V-FITC and 5 µL 7-ADD were added and incubated at room temperature for 15 min without light. The apoptosis rate was detected by using a FACS Aria III Flow cytometry (BD Biosciences, USA).

Cell immuno uorescence analysis
For immuno uorescence staining, the meniscus cells were treated with IL-1β (10ng/ml) for 2 days, and then the cells were xed in 4% paraformaldehyde for 15 minutes at room temperature. After washing with PBS, the cells were permeabilized with 0.3% Triton X-100 for 15 minutes and blocked with 3%BSA at room temperature for 40 minutes. Then, the cells were washed with PBS and incubated with primary antibodies of CDK1(1:200, dilution), AURKB (1:200, dilution), RRM2 (1:100, dilution), and TIMP1(1:150, dilution, Abclonal) at 4℃ overnight. The meniscus cells were rinsed three times with PBS and incubated with secondary antibodies (1:100; ThermoFisher, USA) for 1 h at room temperature the next day. In the end, we stained the nuclei with 4',6-diamidino-2-phenylindole (DAPI, Thermo Fisher, USA) for 5 minutes. All of the images were captured under the uorescence microscope (Nikon, Japan).
qRT-PCR con rmation of the hub genes Total RNA from meniscus tissue was extracted with TRIzol reagent following the manufacturer's protocol. After the purity and concentration of isolated RNA were determined by the nano-drop-2000 nucleic acid analyzer (ThermoFisher, USA), the RNA reverse transcription kits were devoted to reverse transcribe total RNA into cDNA. Quantitative real-time PCR (qRT-PCR) was performed with the ABI StepOnePlus cycler (Applied Biosystems, USA) in a 10-µl reaction system. To accurately quantify the gene expression levels, mRNA levels of the internal reference gene glyceraldehyde phosphate dehydrogenase (GAPDH) was normalized by the 2 −ΔΔCT method. The primer sequences of genes applied in this study were listed in the Table 1.

Statistical analysis
The expression levels of hub genes that were differentially expressed between young and aging tissues or cells were determined by implementing Student's t-test with Prism 6.0 (GraphPad Software, USA), and the data of the experimental results were presented as mean ± standard error of mean (SEM). In addition, for GO functional annotation, DAVID calculated the P value of the false discovery and enrichment rate via Benjamini correction. The P < 0.05 was considered to be statistically signi cant difference.

Identi cation of DEGs
A total of 631 DEGs, including 341 up-regulated genes and 290 down-regulated genes, were identi ed from the GSE45233 dataset by comparing the gerontic and young meniscus samples. The heatmap and volcano plot of total DEGs were shown in Figure. 1A and B, respectively.

GO function and KEGG pathway enrichment analysis of the DEGs
The results of the GO analysis showed that DEGs related to CC were notably concentrated in extracellular region, extracellular space, integral component of plasma membrane, and proteinaceous extracellular matrix ( Figure. 2A). Variations in DEGs connected with BP were mainly concentrated in signal transduction, regulation of transcription from RNA polymerase II promoter, and regulation of cell proliferation ( Figure. 2B). Regarding MF, DEGs were prominently gathered in calcium ion binding, chromatin binding, sequence-speci c DNA binding, and heparin-binding ( Figure. 2C). KEGG pathway analysis manifested that the top typical pathways correlated with DEGs were PI3K-Akt signaling pathway, transcriptional misregulation in cancers, and Jak-STAT signaling pathway ( Figure. 2D).

Construction of PPI network and hub genes analysis
The PPI network of the DEGs was constructed by employing STRING (Figure. 3A). The Degree, MCC, MNC, and EPC algorithms were applied to ltrate hub genes shown in Figure. 3B, C, D, and E, respectively. The mutual hub genes (RRM2, AURKB, CDK1, and TIMP1) were discerned by the Venn diagram ( Figure. 3F). According to Metascape analysis, GO enrichment showed that the signi cant hub genes related to BP were mainly concentrated in aging and regulation of the cell cycle process ( Figure. 3G).
Construction of a lncRNA-miRNA-mRNA network and molecular docking of potential drugs with hub genes As for the screening of miRNAs, since there were no miRNAs that could act on all Hub genes, miRNAs working on three Hub genes were selected as ideal targets. Through the miRWalk website predictive analysis, we screened 5 miRNAs (hsa-miR-6810-5p, hsa-miR-4676-5p, hsa-miR-6877-5p, hsa-miR-8085, and hsa-miR-6133) ( Figure. 4A, Table 2). As illustrated in Figure. 4B, these miRNAs are referred to some pathways that include glycosaminoglycans biosynthesis. To further understand the function of the hub genes in meniscal degeneration, we employed Cytoscape to combine lncRNA/miRNA interactions with miRNA/mRNA interactions, thus constructing a lncRNA-miRNA-mRNA network ( Figure. 4C).
Patients undergoing meniscal resection may have di culty regaining previous motion abilities. [19,20] Therefore, preventing meniscus aging is one of the treatment strategies to reduce the incidence of meniscus tearing. Through probing the DGIdb website, we selected the above hub genes for molecular docking analysis in order to ascertain potential drugs that could act on these genes ( Figure. 5D). The predicted results showed the high-a nity hydrogen binding events between drugs, including Cladribine, Danusertib, Barasertib, Riviciclib, and Dinaciclib, and residues of hub genes, so it was necessary to consider the possible senescence and degeneration of meniscus caused by these drugs in clinical application. The expression of AURKB, RRM2, CDK1, and TIMP1 were decreased in senescent meniscus tissues SafraninO-Fast green and Toloniumchloride Blue staining results of the meniscus of 18-month-old aging mice revealed more severe degenerative histopathological changes than the meniscus in 3-month-old young mice due to decreased collagen ber content ( Figure.5A). Immunohistochemistry showed that the expression levels of AURKB, RRM2, CDK1, and TIMP1 protein were higher in the meniscus of young mice than that of the senescent mice ( Figure.5B and C). In order to further explore the expression levels of hub genes in senescent meniscus, we applied human meniscus for immunohistochemical detection. The results illustrated that the Hub genes were also poorly expressed in senescent human meniscus ( Figure.5D and E). Moreover, the mRNA levels of AURKB, RRM2, CDK1, and TIMP1 were also decreased both in human and mouse menisci, respectively ( Figure.5F). These founds suggested that the decreased gene expressions of AURKB, RRM2, CDK1, and TIMP1 might be a trigger or indicator that participating in meniscus aging and degeneration.

Hub genes expression in cultured mouse meniscus cells
As an additional method to further investigate the expressions of hub genes in senescent meniscus tissues, we treated normal mouse meniscus cells with IL-1β (10ng/ml) for 48h to induce degeneration in mouse meniscus cells in vitro. Take into consideration that the overexpression and accumulation of endogenous lysosomal β-galactosidase are one of the characteristics of senile cells, SA-β-gal is often taken as a biomarkers of senile cells. [21,22] More positive stained cells, the blue ones, were observed in the IL-1β group than the control by the SA-β-gal staining ( Figure. 6A). Meanwhile, more apoptotic cells were detected in the IL-1β group than the control by ow cytometry ( Figure. 6B and D). In addition, the suppressed mRNAs expressions of AURKB, RRM2, CDK1, and TIMP1 were detected in meniscus cells treated with IL-1β (10ng/ml) when compared with the control group ( Figure. 6C). Then, the immuno uorescence assay further supported that these Hub genes were poor expressions in senescent meniscus cells (Figure. 6E). These consequences further con rmed the decreased AURKB, RRM2, CDK1, and TIMP1 expressions in the senescent meniscus.

Discussion
The incidence of meniscal injuries is on the rise, which might be partly attributed to the general aging of the population [23]. Previous studies have found that age-related meniscus denaturation induces an increase in fragility of meniscus tissue, thereby increasing the incidence of meniscus tearing in the elderly population [24]. However, the molecular pathologic mechanism of meniscus age-related degeneration is still poorly studied [25]. The high-throughput microarray technology combined with bioinformatics analysis has been widely used to predict the pathogenesis and potential molecular therapeutic targets of many diseases.
In the current study, a microarray dataset GSE45233 containing young and gerontic meniscus samples was obtained and a bioinformatics analysis was completed. In view of the results of GO terms enrichment analysis, we found a relationship between the DEGs and cell senescence. According to the CC analysis results of DEGs, we found that the major changes in the components of the cells were located in the extracellular region. The extracellular domain of meniscus cells mainly refers to the extracellular matrix (ECM), which plays a crucial role in sustaining structural integrity and mechanical properties of meniscus [26]. Histopathological analysis showed that meniscus aging was associated with ECM structure loss [27]. Variations in DEGs related to MF were primarily concentrated in calcium ion binding, and chromatin binding. In fact, imbalance of calcium ion homeostasis and chromatin rearrangement are also often associated with cell senescence and ECM degeneration [28,29]. Based on the KEGG database, we discovered that the DEGs were principally gathered in the PI3K-Akt signaling pathway, which could act on downstream targets including forkhead box O transcription factors (FoxO) and mammalian target of rapamycin (mTOR), and also activate gene P53 to regulate cell growth [30][31][32]. Herranz et al. found that mTOR controls senescence secretory by indirectly regulating the expressions of senescence signaling molecules [33]. Furthermore, FoxO may play a protective role against meniscus aging by regulating autophagy [34]. Such reports further con rmed our nding.
The PPI network was carried out to predict the connections of protein functions of DEGs, and the common hub genes-AURKB, CDK1, RRM2, and TIMP1-were screened. CDK1 is a catalytic subunit of mitogenic promotors and plays a key role in cell cycle control, such as mitosis and G2-M transformation [35]. Saito et al. demonstrated that CDK1 was essential for chondrocyte proliferation and differentiation [36]. AURKB is a type of Aurora kinase involved in the regulation of chromosome alignment and separation in mitosis and meiosis [37]. Thus, both AURKB and CDK1 could re ect the meniscus cell cycle to a certain degree. RRM2 encodes one of two non-identical subunits for ribonucleotide reductase, and catalyzes the formation of deoxyribonucleotides from ribonucleotides [38]. Previous studies have found that the decrease in RRM2 levels could lead to DNA damage accumulation by inhibiting deoxyribonucleotide triphosphates (dNTPs) expression, thereby mediating cell senescence [39]. TIMP1, natural inhibitors of the matrix metalloproteinases (MMPs), is a group of peptidases involved in degradation of the extracellular matrix [40]. The balance between matrix degradation and synthesis is impaired during meniscus aging. Since elastolysis or collagenolysis is put down to the joint effect of several members of MMPs, TIMP1 may also play an anti-meniscus aging role to some extent [41]. In the current study, the decrease in the expressions of AURKB, CDK1, RRM2, and TIMP1 was observed in both the human and mouse senescent menisci, as well as meniscus cells treated with exogenous IL-1β. Such nding indicated that both cell cycle arrest and ECM degeneration contribute to the meniscus senescence.
For the ve selected key miRNAs targeting the hub genes, there were no relevant studies reporting that they were linked to meniscus before. Further pathways enrichment analysis of these miRNA revealed that these miRNAs were involved in the synthesis of glycosaminoglycans, an important component of the extracellular matrix of the meniscus, of which, the reduction was often accompanied by the occurrence of OA [42]. With regard to the drugs (Cladribine, Danusertib, Barasertib, Riviciclib, and Dinaciclib) targeting hub genes based on pharmacogenomics tools prediction, their toxic and side effects on meniscus aging and damage should also be taken into account in the treatment of related diseases [43]. Cellular senescence is a complex phenotype characterized by two aspects: persistent cell cycle arrest and the production of pro-in ammatory molecules known as senescence associated secretion phenotype (SASP) [44]. Considering that the senescence of articular chondrocytes acting as hyporeplicative cells was mainly driven by SASP factors, we utilized IL-1β (10ng/ml) to construct an in vitro cell senescence model [45]. Our current study has shown that cell senescence was increased in IL-1β-treated mouse meniscus cells. The decreased expressions of screened hub genes in vitro meniscus cells supported in vivo experiments and further identi ed these genes were related to cell senescence.
However, our conclusions are only originated from theoretical prediction and models detection, and the role of pivotal genes in meniscus senescence still needs further clinical veri cation. Still, our current study would provide a basis for nding markers of the aging meniscus to a certain extent.

Declarations
Ethics approval and consent to participate This human research protocol was authorized by the Medical Ethics Committees of Renmin Hospital of Wuhan University (approval number: 2019K-K011). All patients have signed informed consent. Animal protocols used in this study were approved by the Animal Welfare Committee in Wuhan University (License number: 14016). These approvals included consents to the collection, storage, and testing of samples for research conducted in this manuscript.

Not Applicable
Availability of data and materials The protein sequences and structures from this study were deposited in RCSB PDB database and were accessible through PDB ticket number 6C2R, 2B9R, and 6LE1, at: https://www.rcsb.org/structure/6C2R, https://www.rcsb.org/structure/2B9R, and https://www.rcsb.org/structure/6LE1. The DNA and RNA sequencing data were downloaded from Genebank database through gene names RRM2 (Gene ID:

Competing interests
No competed or con icts of interest reported. No bene ts in any form have been received or will be received from a commercial party related directly or indirectly to the subject of this article.