Molecular Medicine Reports Molecular Medicine Reports 1791-2997 1791-3004 D.A. Spandidos 10.3892/mmr.2018.9685 mmr-19-01-0093 Articles Microarray analysis of the molecular mechanisms associated with age and body mass index in human meniscal injury HuangPeiyan GuJun WuJunguo GengLei HongYang WangSiqun WangMinghai Department of Orthopedic Surgery, Shanghai Fifth People's Hospital Affiliated to Fudan University, Shanghai 200240, P.R. China Correspondence to: Dr Minghai Wang, Department of Orthopedic Surgery, Shanghai Fifth People's Hospital Affiliated to Fudan University, 128 Ruili Road, Minhang, Shanghai 200240, P.R. China, E-mail: wwmmhai2@163.com 012019 21112018 19 1 93 102 04022018 04102018 Copyright: © Huang et al. 2019 This is an open access article distributed under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made.

The aim of the present study was to identify genes and functional pathways associated with meniscal injuries affected by age or body mass index (BMI) using microarray analysis. The GSE45233 gene expression dataset with 12 injured meniscus samples associated with age and BMI and GSE66635 dataset with 12 injured and 12 normal meniscus samples were downloaded from the Gene Expression Omnibus database. Differentially expressed genes (DEGs) were identified based on age or BMI in GSE45233. DEGs between injured and normal meniscus samples in GSE66635 were also identified. Common DEGs between GSE45233 and GSE66635 were identified as feature genes associated with age or BMI, followed by protein-protein interaction (PPI) network and functional pathway enrichment analyses for the feature genes. Finally, the GSE51588 genome-wide expression profile was then downloaded from the GEO database to validate the results. A total of 1,328 DEGs were identified. Of these, 28 age-associated and 20 BMI-associated meniscal injury genes were obtained. B-cell lymphoma-2 (Bcl-2) and matrix metalloproteinase-14 were identified as hub genes in the PPI networks. Functional pathway enrichment analysis revealed that vascular endothelial growth factor A (VEGFA), transferrin (TF) and Bcl-2 were involved in the hypoxia-inducible factor 1 signaling pathway. TF was involved in the mineral absorption function pathway associated with BMI. Additionally, TF and VEGFA were identified to be overlapping candidate genes of GSE45233 and GSE66635, and DEGs in GSE51588. Therefore, VEGFA, TF, and Bcl-2 may be important genes for human meniscal injuries. Additional evaluations of these results are required.

 meniscal injury body mass index age differentially expressed genes protein-protein interaction network functional enrichment analysis Introduction

Meniscal injury is a common disorder that may cause pain, swelling or mechanical symptoms in patients (1). As articular cartilages of the knee joint, menisci are crucial for normal function of the knee joint and protection of the articular surfaces (2). Articular cartilages, as specialized connective tissues of diarthrodial joints, are affected by overall fitness levels (3,4).

Previously, the effect of factors including aging and obesity on articular cartilages has been widely studied. With aging, changes occur in articular cartilages, including degenerative changes in the morphology, loss of cartilage matrix proteins, cleavage of type II collagen, worsening of mechanical properties and dysregulated expression of associated genes (4–7). An inverse correlation between cartilage thickness within the knee joint and age has been demonstrated (8). The mechanisms of obesity affecting articular cartilages are complex, and include biomechanical effects via increased loading and meta-inflammation (9–14). Abnormalities in the articular cartilage of the knee are also associated with obesity, as revealed by the investigation of obese adolescent patients suffering from knee pain (15).

Details of the gene expression profile in human meniscal injury have been described previously (16–18). Expression of vascular endothelial growth factor A (VEGFA), tumor necrosis factor-α and matrix metalloproteinases (MMPs) were demonstrated to be varied by age in patients with meniscal tears (16,19). Neuronal apoptosis-inhibitory protein, apoptosis inhibitor of macrophage, carbohydrate sulfotransferase 15 and MMP28 were suggested to be associated with obesity (17). However, gene expression changes associated with injuries and aging or obesity in the human meniscus have not been fully investigated.

Microarray analysis of gene expression changes and functional pathways may be useful to improve understanding of the associations between aging or obesity and meniscal injuries. In previous studies, GSE45233 microarray data was used to study variation in gene expression signatures in human injured meniscus with age or body mass index (BMI) and the degree of chondrosis in the knee (16,17). GSE66635 microarray data was used to compare transcriptome signatures between the meniscus and articular cartilage from knees undergoing arthroscopic partial meniscectomy (18). In the present study, the GSE45233 and GSE66635 datasets were downloaded and analyzed to identify genes and functional pathways associated with age or BMI in human injured meniscus. Differentially expressed genes (DEGs) were identified and feature genes were obtained by comparing the DEGs from the 2 datasets. Protein-protein interaction (PPI) networks of these feature genes were constructed to obtain hub genes. Functional pathway enrichment analysis for feature genes was also performed to identify key pathways in the development of age or BMI-associated meniscal injuries. Furthermore, the GSE51588 genome-wide expression profile was downloaded from the Gene Expression Omnibus (GEO) database to validate the results. The results concerning gene expression changes in human injured meniscus may provide an improved understanding of the associations between meniscal injuries and aging or obesity.

 Materials and methods <sec> <title>Affymetrix microarray data

The gene expression profile datasets GSE66635, deposited by Rai et al (18), and GSE45233, deposited by Rai et al (16,17), were downloaded from the GEO (www.ncbi.nlm.nih.gov/geo/) database. GSE66635 and GSE45233 were based on the GPL16686 Human Gene 2.0 ST Array (Affymetrix Inc.; Thermo Fisher Scientific, Inc., Waltham, MA, USA) and GPL10558 HumanHT-12 V4.0 expression Beadchip (Illumina, Inc., San Diego, CA, USA) platforms, respectively. A total of 24 samples were available in the GSE66635 dataset, including 12 injured meniscus and 12 normal articular cartilages from patients undergoing partial meniscectomy. The second microarray profile, GSE45233, included 12 isolated injured meniscus tissues from patients undergoing arthroscopic partial meniscectomy.

 Data preprocessing and feature genes screening

The probe-level data in the CEL files were converted to expression measures, and idle probes were deleted. When multiple probes corresponded to the same gene, the expression values of those probes were averaged. The 2 datasets were aggregated according to gene names, followed by background correction and quartile data normalization by the robust multiarray average algorithm in R language (20).

Based on the GSE66635 gene expression profile data, DEGs in injured menisci from patients undergoing partial meniscectomy were identified using the Linear Models for Microarray Data (Limma) package (21) in R language compared with normal articular cartilage. The cut-off values were P<0.05 and fold change (FC) >1.5.

Samples of GSE45233 were grouped by age and BMI according to World Health Organization standards (lean, BMI <18.5; normal, BMI 18.5–25; obese, BMI >25 kg/m2) (22). The age group included 6 samples of patients >40 years of age (old group) and 6 samples of patients <40 years of age (young group). DEGs associated with age were then identified in the old group using the Limma package and compared with the young group.

The BMI group contained 8 samples of lean patients (lean group) and 4 samples of obese patients (obese group); no patients were in the normal category according to this dataset. DEGs associated with age were identified in the obese group using the Limma package and compared with the lean group.

Common DEGs between GSE66635 and the age group in GSE45233 were defined as ‘feature’ genes associated with age-associated meniscal injuries. Similarly, common DEGs identified from GSE66635 and the BMI group in GSE45233 were BMI-associated feature genes.

 Integrative analysis of feature gene expression levels

A total of 9 samples of patients >40 years of age and 3 samples of patients <40 years of age were selected from the injured meniscus samples in GSE66635. Combining these samples with those in the age group from GSE45233, a new age group with 15 samples of patients >40 years of age and 9 samples of patients <40 years of age was obtained. Similarly, 8 obese samples from GSE66635 were selected and combined with those in the BMI group from GSE45233, obtaining a new BMI group with 8 lean samples and 12 obese samples. The significances of average expression level differences for feature genes extracted from the two new groups were analyzed by t-test using the Limma package (23). Heat maps were generated using Pheatmap (24) in R language, to exhibit the results of bidirectional hierarchical clustering based on the expression value of the feature genes.

 Classification of the sample classification model

Support vector machines (SVM) (25) are a useful technique for two-group classification problems. To verify the identifiability and classifiability of the feature genes extracted from the two new groups, the SVM classifier model was performed based on the expression value of these feature genes. SVMs were initially optimized using the DEGs in the training set and the feature genes were examined in the test set to separate the samples in the new age and BMI groups.

 PPI network construction

The Human Protein Reference Database (HPRD; www.hprd.org/) (26) is a protein database that contains information on human protein functions, including PPIs, post-translational modifications, and enzyme-substrate and disease associations. PPI networks for age-associated feature genes and BMI-associated feature genes were respectively constructed by mapping the two groups of feature genes to the PPI pairs downloaded from HPRD, which were then visualized using Cytoscape (27). Finally, hub genes with the highest node degree in the PPI networks were identified using CytoNCA (28) plugin (network without weight) in Cytoscape.

 Pathway enrichment analysis

The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database (www.genome.jp/kegg/pathway.html) (29) is used for the systematic analysis of gene functions and associated genomic information with higher order functional information. The feature genes were examined using the KEGG pathway database to identify functional pathways. The enrichment process was performed by the Fisher algorithm using KOBAS2.0 (30) as follows:

P=1-∑i=0x-1(Mi)(N-MK-i)(NK)

where N is the number of total genes in the genome, M is the number of genes in the pathway and K is the number of DEGs. The Fisher score indicates the probability of X functional pathway genes among K DEGs.

 Data validation

To validate the results, the GSE51588 genome-wide expression profile was downloaded from the GEO database. This dataset was obtained from 40 subchondral bone tissue samples from patients with osteoarthritis and 10 samples from normal subjects based on the GPL13497 Agilent-026652 Whole Human Genome Microarray 4×44 K v2 platform. The original data in TXT format were preprocessed using the Limma package (version 3.30.13; www.bioconductor.org/packages/2.9/bioc/html/limma.html) (21). The preprocessing methods included background correction and data normalization. DEGs were identified using the Bayes method in Limma with a threshold P<0.05.

The overlapping DEGs between GSE66635 and GSE45233 were selected as candidate genes, which were then compared with the DEGs obtained from GSE51588.

 Results <sec> <title>Screening of candidate feature genes

Following data normalization, genes exceeding the difference threshold (P<0.05 and FC >1.5) were screened. A total of 364 DEGs were identified in injured meniscus samples compared with normal meniscus samples in GSE66635, among which 135 DEGs were down-regulated and 229 DEGs were upregulated. A total of 424 DEGs were screened in the old group compared with the young group in GSE45233, including 222 downregulated genes and 202 upregulated genes. In the obese group, compared with the lean group, 540 genes were differentially expressed, including 340 downregulated and 200 upregulated genes. Comparison of the DEGs in the GSE66635 dataset with the DEGs of the age group in GSE45233 identified 28 age-associated meniscal injury feature genes (Table I). Of these, 15 feature genes were downregulated, which included S100 calcium-binding protein A1 (S100A1) and BARX homeobox 2, and 13 feature genes were upregulated, which included cluster of differentiation 36 and endomucin (EMCN). A total of 20 common genes between the GSE45233 dataset and the BMI group in GSE45233 were identified as BMI-associated meniscal injury feature genes (Table I). Of these, 11 feature genes were downregulated, which included EMCN and furry, and 9 feature genes were upregulated, which included IGF binding protein 1 and S100 calcium-binding protein A8 (S100A8).

 Integrative analysis of expression levels of feature genes

Following integration of the GSE66635 and GSE45233 datasets, a new age-associated meniscal injury group including 24 samples and a new BMI-associated meniscal injury group including 20 samples were obtained. Intra-group differences of the 28 age-associated meniscal injury feature genes and 20 BMI associated meniscal injury feature genes were examined by t-test. For the age group, the average expression levels of the young group were significantly increased compared with those of the old group (P=0.021); for the BMI group, the average expression levels of the lean group were increased compared with those of the obese group, but the difference was not significant (P=0.3663; Fig. 1A). Heat maps (Fig. 1B and C) of bidirectional hierarchical clustering indicated that the samples within the age group were completely separated into the young and old groups, while 1 obese sample (obese-12) in the BMI group was assigned to the lean group.

 Classification of SVM classifier model

To observe the classifiability of the feature genes from the age and BMI classifications, the SVM classifier model was utilized. Samples in the age group were divided into two groups, but 1 old sample was included with the young samples (Fig. 2A). The samples in the BMI group were completely divided into lean and obese groups (Fig. 2B).

 PPI network analysis

Analysis of the HPRD database identified 215 unique PPI pairs for the age group and 90 unique PPI pairs for the BMI group. The PPI network for age-associated feature genes contained 222 nodes, including 18 feature genes (Fig. 3A). The PPI network for BMI-associated feature genes contained 102 nodes, including 15 feature genes (Fig. 3B). Network topology analysis indicated that the constructed PPI networks obeyed scale-free network attributes (Fig. 4) and the node degree of the network followed the distribution, obtaining y = 64.6 * x−2.57 for the age group and y = 47.82 * x−2.46 for the BMI group, where × represents node degree. The top 5 nodes with the higher node degree of the age group were B-cell lymphoma-2 (Bcl-2), CyclinB1, S100A1, VEGFA and Transferrin (TF) (Table II). The top 5 nodes of the BMI group were MMP14, S100A8, TF, chemokine (C-X-C motif) ligand 12 and endothelin receptor type A (Table II).

 Pathway enrichment analysis

KEGG pathway enrichment analyses for the top 10 feature genes were performed. Feature genes associated with age-associated meniscal injury were significantly enriched in the hypoxia-inducible factor 1 (HIF-1) signaling pathway, including VEGFA, TF, and Bcl-2 (Table III). Feature genes associated with BMI-associated meniscal injury were significantly associated with the mineral absorption function pathway, in which TF was enriched.

 Data validation

From GSE51588, a total of 3,403 DEGs (1,853 upregulated and 1,550 downregulated) were identified (data not shown). Following comparison of GSE66635 and GSE45233, a total of 28 and 20 common DEGs (candidate genes) associated with age and BMI, respectively, were identified. Finally, 11 overlapping genes were obtained between candidate genes and DEGs in GSE51588 (Table IV). These included TF and VEGFA.

 Discussion

Meniscal injury is a common disease caused by athletic events and activities in daily life. Patient age and BMI are associated meniscal injuries (17,18). An understanding of the molecular mechanism of the associations between meniscal injuries and age or BMI is important. In the present study, 2 datasets (GSE45233 and GSE66635) associated with meniscal injury were downloaded from the GEO database and analyzed. Feature genes associated with age-associated meniscal injury, including VEGFA, TF and Bcl-2 were involved in the HIF-1 signaling pathway. In addition, the feature gene TF was also associated with BMI-associated meniscal injury and was significantly enriched in the mineral absorption function pathway.

Bcl-2 is a regulator of apoptosis from the B-cell lymphoma-2 family, and its function is to inhibit cell death, rather than promote proliferation (18,31–33). Overexpression of Bcl-2 in mouse bones suppresses apoptosis of bone cells (34) and suppression of Bcl-2 increases apoptosis in leukemic cells (35). Meniscal injuries are directly associated with the development of osteoarthritis (36). Emerging evidence indicates that apoptosis serves an important role in osteoarthritis pathology (37). Iwata et al (38) demonstrated that the levels of apoptosis in chondrocytes in mice consuming a high-fat diet was increased. The present study identified that Bcl-2 was downregulated in the injured meniscus samples compared with the normal meniscus, meaning that apoptosis may be increased in injured menisci tissues, which was in accordance with data from Iwata et al (38). The expression of Bcl-2 was downregulated in the old group compared with the young group, in accordance with a previous study (16). Bcl-2 was the hub gene with a node degree of 74 in the PPI network for the feature genes associated with age in the injured meniscus. All these data imply that aging may accelerate apoptosis of the injured menisci of older patients.

VEGFA is an important angiogenetic protein with a selective mitogenic effect on vascular endothelial cells (39). As a vascular endothelial growth factor, VEGFA was demonstrated to be necessary for the survival of chondrocytes during skeletal development (40) and as a regulator of osteoblast differentiation during bone development (41). VEGFA may increase the osteogenic healing capacity by promoting osteogenic and endothelial differentiation (42). A previous study revealed that VEGF levels were increased following the creation of meniscal lesions in rabbits (43). Chen et al (39) recently identified that VEGFA was downregulated in osteoarthritis chondrocytes compared with normal chondrocytes. In the present study, the expression of VEGFA was downregulated in the injured meniscus compared with the normal meniscus tissues, which is consistent with previous data that the healing capacity of menisci tissue may be damaged in injured meniscus, as the angiogenesis function of VEGFA is be inhibited in the downregulated condition (44,45). Furthermore, the expression of VEGFA was downregulated in the old group compared with the young group, consistent with a previous study (16). Therefore, the healing capacity for injured meniscus of older patients may be not as efficient compared with younger patients.

In the present study, TF was identified as a meniscal injury feature gene associated with age and BMI. TF is an iron transport protein, which delivers iron from absorption centers and storage sites to all tissues as it circulates in the plasma (46). TF serves a major role in angiogenesis during endochondral bone formation and is produced by hypertrophic chondrocytes (47). Hypertrophy and neovascularization are specific signs of healing of injured menisci (48), through which TF may be produced and transported. Based on these data, the upregulation of TF in the injured meniscus indicates that the healing of the injured meniscus may have occurred, and that TF may have a role in the healing progress. However, this hypothesis requires experimental verification. TF was identified to be downregulated in the old group compared with the young group in the present study, which is consistent with previous data (16), indicating that recovery of the injured meniscus in old patients may be slower. Patient age has a detrimental effect on the healing potential of the injured meniscus (48). The reason why the expression of TF was upregulated in the obese group compared with the lean group is not clear, as a tendency for lower TF saturation in obesity has been suggested (49). Additional investigation is required to verify the association between TF and healing of injured meniscus.

One of the key results of the present study was that the feature genes in the age group were enriched in the HIF-1 signaling pathway and the feature genes in the BMI group were enriched in the mineral absorption function pathway. HIF-1 is a transcriptional activator that targets genes encoding proteins that increase O2 delivery and mediate adaptive responses to O2 deprivation (50,51). HIF-1α may protect articular cartilage by promoting the chondrocyte phenotype, maintaining chondrocyte viability and supporting metabolic adaptation to a hypoxic environment (52). Articular cartilage is a hypoxic tissue, in which HIF-1 is essential for survival, growth, energy generation and matrix synthesis of chondrocytes (33,44). HIF-1 is the major regulator of VEGFA in the context of angiogenesis (53,54). VEGFA may promote angiogenesis (44,45), which is an indication of the healing process in injured menisci (48). Expression of oxygen-regulated TF is positively mediated by HIF-1 (55), which occurs in the hypoxic meniscal tissues. In addition, Bcl-2 may promote HIF-1-mediated VEGF expression under hypoxia by increasing the expression level of HIF-1 (56). In the present study, VEGFA and Bcl-2 were downregulated in the injured meniscus compared with the normal meniscus tissue. This indicates that the HIF-1 signaling pathway was inhibited. These results indicate that apoptosis may be increased, and angiogenesis may be decreased in meniscal injury tissues. The expression of VEGFA, TF and Bcl-2 was downregulated in the old group compared with the young group, which may explain the longer time generally required for older patients to recover from meniscal injuries compared with younger patients (57). However, TF was upregulated in the injured meniscus, which was in contrast to the data that TF is positively mediated by HIF-1 under hypoxic condition (55). This may be a result of co-regulation by the HIF-1 signaling pathway and the mineral absorption function pathway. An additional contradiction is that TF was upregulated in the obese group and, contrarily, TF saturation in obesity is usually decreased (49). Experimental verification of these contradictions is important for additional understanding of the mechanisms of meniscal injuries.

There are certain limitations in the present study. The key genes identified were not validated in clinical or in vivo studies, although data validation was performed. Additionally, the sample size was small. Therefore, future clinical and in vivo studies will be performed to validate the results of the present study, and to investigate the underlying molecular mechanisms of the effects of aging and obesity on meniscal injuries.

Aging may affect the development of meniscal injuries through the HIF-1 signaling pathway, in which VEGFA, TF, and Bcl-2 are involved. Obesity may affect the mineral absorption function pathway in injured menisci by regulating the expression of TF. However, there remains a lack of experimental evidence to confirm the hypotheses in the present study, and the molecular mechanisms of the effects of aging and obesity on these pathways remain unclear. In addition, why and how TF is involved in the HIF-1 signaling pathway and the mineral absorption function pathway requires additional study.

 Acknowledgements

Not applicable.

 Funding

The present study was funded by the Minhang District Natural Science Project (grant no. 2017MH77).

 Availability of data and materials

The analyzed data sets generated during the study are available from the corresponding author on reasonable request.

 Authors' contributions

PH contributed to data collection, analysis and interpretation, obtained funding and wrote the manuscript. SW and YH contributed to the study design. LG performed the experiments and contributed to the data collection. JG performed data analysis and interpretation. JW conducted statistical analysis. MW contributed to data collection and interpretation. All authors read and approved the final manuscript.

 Ethics approval and consent to participate

Not applicable.

 Patient consent for publication

Not applicable.

 Competing interests

The authors declare that they have no competing interests.

 References GreisPEBardanaDDHolmstromMCBurksRTMeniscal injury: I. Basic science and evaluationJ Am Acad Orthop Surg10168176200210.5435/00124635-200205000-0000312041938 NobleJTurnerPGThe function, pathology, and surgery of the meniscusClin Orthop Relat Res626819863757377 AnandacoomarasamyALeibmanSSmithGCatersonIGiuffreBFransenMSambrookPNMarchLWeight loss in obese people has structure-modifying effects on knee articular cartilageAnn Rheum Dis712632201210.1136/ard.2010.14472522135412 MusumeciGSzychlinskaMAMobasheriAAge-related degeneration of articular cartilage in the pathogenesis of osteoarthritis: Molecular markers of senescent chondrocytesHistol Histopathol30112201510.7243/2055-091X-2-125010513 ArmstrongCGMowVCVariations in the intrinsic mechanical properties of human articular cartilage with age, degeneration and water contentJ Bone Joint Surg Am648894198210.2106/00004623-198264010-000137054208 BuckwalterJARoughleyPJRosenbergLCAge-related changes in cartilage proteoglycans: Quantitative electron microscopic studiesMicrosc Res Tech28398408199410.1002/jemt.10702805067919527 BankRABaylissMTLafeberFPMaroudasATekoppeleJMAgeing and zonal variation in post-translational modification of collagen in normal human articular cartilage. The age-related increase in non-enzymatic glycation affects biomechanical properties of cartilageBiochem J330345351199810.1042/bj330034594615291219146 DallaPLCovaMPozzi-MucelliRSMRI appearance of the articular cartilage in the knee according to ageJ Belge Radiol80172019979103709 JonesDGArticular cartilage degeneration: Etiologic association with obesityOchsner J91371392009216034293096275 SowersMRKarvonen-GutierrezCAThe evolving role of obesity in knee osteoarthritisCurr Opin Rheumatol22533537201010.1097/BOR.0b013e32833b4682204851733291123 MezhovVCiccutiniFMHannaFSBrennanSLWangYYUrquhartDMWlukaAEDoes obesity affect knee cartilage?A systematic review of magnetic resonance imaging data. Obes Rev151431572014 BlazekKFavreJAsayJErharthledikJAndriacchiTAge and obesity alter the relationship between femoral articular cartilage thickness and ambulatory loads in individuals without osteoarthritisJ Orthop Res32394402201410.1002/jor.2253024281940 TravascioFEltoukhyMCamiSAsfourSAltered mechano-chemical environment in hip articular cartilage: Effect of obesityBiomech Model Mechanobiol13115201410.1007/s10237-013-0545-523483348 DenningWMWinwardJGPardoMBHopkinsJTSeeleyMKBody weight independently affects articular cartilage catabolismJ Sports Sci Med142902962015259835774424457 WidhalmHKMarlovitsSWelschGHDirisamerANeuholdAvan GriensvenMSeemannRVécseiVWidhalmKObesity-related juvenile form of cartilage lesions: A new affliction in the knees of morbidly obese children and adolescentsEur Radiol22672681201210.1007/s00330-011-2281-521947483 RaiMFPatraDSandellLJBrophyRHTranscriptome analysis of injured human meniscus reveals a distinct phenotype of meniscus degeneration with agingArthritis Rheum6520902101201310.1002/art.3798423658108 RaiMFPatraDSandellLJBrophyRHRelationship of gene expression in the injured human meniscus to body mass index: A biologic connection between obesity and osteoarthritisArthritis Rheumatol6621522164201410.1002/art.38643246921314116431 RaiMFSandellLJZhangBWrightRWBrophyRHRNA microarray analysis of macroscopically normal articular cartilage from knees undergoing partial medial meniscectomy: Potential prediction of the risk for developing osteoarthritisPLoS One11e0155373201610.1371/journal.pone.0155373271710084865200 BrophyRHRaiMFZhangZTorgomyanASandellLJMolecular analysis of age and sex-related gene expression in meniscal tears with and without a concomitant anterior cruciate ligament tearJ Bone Joint Surg Am94385393201210.2106/JBJS.K.00919223624943284860 IrizarryRAHobbsBCollinFBeazer-BarclayYDAntonellisKJScherfUSpeedTPExploration, normalization, and summaries of high density oligonucleotide array probe level dataBiostatistics4249264200310.1093/biostatistics/4.2.24912925520 SmythGKLimma: Linear models for microarray data. In: Bioinformatics and computational biology solutions using R and bioconductorGentlemanRCareyVJHuberWIrizarryRADudoitSSpringerNew York, NY3974202005 AppiahCAWorld Health Organisation. Global Database on Body Mass IndexWorld Health Organisation: Geneva 20062014 RitchieMEPhipsonBWuDHuYLawCWShiWSmythGKlimma powers differential expression analyses for RNA-sequencing and microarray studiesNucleic Acids Res43e47201510.1093/nar/gkv007256057924402510 KoldeRPheatmapPretty Heatmaps2015 CortesCVapnikVSupport-vector networksMach Learn20273297199510.1023/A:1022627411411 Keshava PrasadTSGoelRKandasamyKKeerthikumarSKumarSMathivananSTelikicherlaDRajuRShafreenBVenugopalAHuman protein reference database-2009 updateNucleic Acids Res37D767D772200910.1093/nar/gkn89218988627 ShannonPMarkielAOzierOBaligaNSWangJTRamageDAminNSchwikowskiBIdekerTCytoscape: A software environment for integrated models of biomolecular interaction networksGenome Res1324982504200310.1101/gr.123930314597658403769 TangYLiMWangJPanYWuFXCytoNCA: A cytoscape plugin for centrality analysis and evaluation of protein interaction networksBio Systems1276772201510.1016/j.biosystems.2014.11.00525451770 OgataHGotoSSatoKFujibuchiWBonoHKanehisaMKEGG: Kyoto encyclopedia of genes and genomesNucleic Acids Res272934199910.1093/nar/27.1.299847135148090 XieCMaoXHuangJDingYWuJDongSKongLGaoGLiCYWeiLKOBAS 2.0: A web server for annotation and identification of enriched pathways and diseasesNucleic Acids Res39W316W322201110.1093/nar/gkr483217153863125809 CorySAdamsJMThe Bcl2 family: Regulators of the cellular life-or-death switchNat Rev Cancer2647656200210.1038/nrc88312209154 VauxDLCorySAdamsJMBcl-2 gene promotes haemopoietic cell survival and cooperates with c-myc to immortalize pre-B cellsNature335440442198810.1038/335440a03262202 Fernández-TorresJMartínez-NavaGAGutiérrez-RuízMCGómez-QuirozLEGutiérrezMRole of HIF-1α signaling pathway in osteoarthritis: A systematic reviewRev Bras Reumatol Engl Ed571621732017(In Portuguese)28343622 YamashitaJDattaNSChunYHYangDYCareyAAKreiderJMGoldsteinSAMcCauleyLKRole of Bcl2 in osteoclastogenesis and PTH anabolic actions in boneJ Bone Miner Res23621632200810.1359/jbmr.07121118086008 CimminoACalinGAFabbriMIorioMVFerracinMShimizuMWojcikSEAqeilanRIZupoSDonoMmiR-15 and miR-16 induce apoptosis by targeting BCL2Proc Natl Acad Sci USA1021394413949200510.1073/pnas.050665410216166262 WilliamsLBAdesidaABAngiogenic approaches to meniscal healingInjury49467472201810.1016/j.injury.2018.01.02829395218 ZamliZSharifMChondrocyte apoptosis: A cause or consequence of osteoarthritis?Int J Rheum Dis14159166201110.1111/j.1756-185X.2011.01618.x21518315 IwataMOchiHHaraYTagawaMKogaDOkawaAAsouYInitial responses of articular tissues in a murine high-fat diet-induced osteoarthritis model: Pivotal role of the IPFP as a cytokine fountainPLoS One8e60706201310.1371/journal.pone.0060706235932893625196 ChenHTianYMiR-15a-5p regulates viability and matrix degradation of human osteoarthritis chondrocytes via targeting VEGFABiosci Trends10482488201710.5582/bst.2016.0118727916780 ZelzerEMamlukRFerraraNJohnsonRSSchipaniEOlsenBRVEGFA is necessary for chondrocyte survival during bone developmentDevelopment13121612171200410.1242/dev.0105315073147 DuanXMurataYLiuYNicolaeCOlsenBRBerendsenADVegfa regulates perichondrial vascularity and osteoblast differentiation in bone developmentDevelopment14219841991201510.1242/dev.117952259773694460736 BehrBTangCGermannGLongakerMTQuartoNLocally applied VEGFA increases the osteogenic healing capacity of human adipose-derived stem cells by promoting osteogenic and endothelial differentiationStem Cells29286296201110.1002/stem.581217324863400547 RuizIbán MÁComellas MeleroNMartinez-BotasJOrtizADiaz HerediaJGrowth factor expression after lesion creation in the avascular zone of the meniscus: A quantitative PCR study in rabbitsArthroscopy3011311138201410.1016/j.arthro.2014.04.08924951133 DeckersMMvan BezooijenRLvan der HorstGHoogendamJvan Der BentCPapapoulosSELöwikCWBone morphogenetic proteins stimulate angiogenesis through osteoblast-derived vascular endothelial growth factor AEndocrinology14315451553200210.1210/endo.143.4.871911897714 WeijtsBGBakkerWJCornelissenPWLiangKHSchaftenaarFHWestendorpBde WolfCAPaciejewskaMScheeleCLKentLE2F7 and E2F8 promote angiogenesis through transcriptional activation of VEGFA in cooperation with HIF1EMBO J3138713884201210.1038/emboj.2012.231229030623463843 BonkovskyHLIron and the liverAm J Med Sci3013243199110.1097/00000441-199101000-000061847276 CarlevaroMFAlbiniARibattiDGentiliCBenelliRCermelliSCanceddaRCanceddaFDTransferrin promotes endothelial cell migration and invasion: Implication in cartilage neovascularizationJ Cell Biol13613751384199710.1083/jcb.136.6.137590874502132523 SenanVSucheendranJPrasad and BalagopalKHistological features of meniscal injuryKerala J Orthop2430362011 ChengHLBryantCCookRO'ConnorHRooneyKSteinbeckKThe relationship between obesity and hypoferraemia in adults: A systematic reviewObes Rev13150161201210.1111/j.1467-789X.2011.00938.x21981048 SemenzaGSignal transduction to hypoxia-inducible factor 1Biochem Pharmacol64993998200210.1016/S0006-2952(02)01168-112213597 SemenzaGLHypoxia-inducible factor 1: Oxygen homeostasis and disease pathophysiologyTrends Mol Med7345350200110.1016/S1471-4914(01)02090-111516994 ZhangFJLuoWLeiGHRole of HIF-1α and HIF-2α in osteoarthritisJoint Bone Spine82144147201510.1016/j.jbspin.2014.10.00325553838 PagèsGPouysségurJTranscriptional regulation of the vascular endothelial growth factor gene-a concert of activating factorsCardiovasc Res65564573200510.1016/j.cardiores.2004.09.03215664382 LiaoDJohnsonRSHypoxia: A key regulator of angiogenesis in cancerCancer Metastasis Rev26281290200710.1007/s10555-007-9066-y17603752 RolfsAKvietikovaIGassmannMWengerRHOxygen-regulated transferrin expression is mediated by hypoxia-inducible factor-1J Biol Chem2722005520062199710.1074/jbc.272.32.200559242677 TrisciuoglioDGabelliniCDesideriMRagazzoniYDe LucaTZiparoEDel BufaloDInvolvement of BH4 domain of bcl-2 in the regulation of HIF-1-mediated VEGF expression in hypoxic tumor cellsCell Death Differ1810241035201110.1038/cdd.2010.175212338463131942 RoosHAdalberthTDahlbergLLohmanderLSOsteoarthritis of the knee after injury to the anterior cruciate ligament or meniscus: The influence of time and ageOsteoarthritis Cartilage3261267199510.1016/S1063-4584(05)80017-28689461

Comparison of the average expression levels between groups of feature genes. (A) A t-test of the average expression level of feature genes was performed. (B) Heat map of bidirectional hierarchical clustering for feature genes in the age group. Comparison of the average expression levels between groups of feature genes. (C) Heat map of bidirectional hierarchical clustering for feature genes in the body mass index group.

Scatter plots of support vector machines classification of the samples. (A) Age and (B) body mass index groups are presented. The orange color represents the young and obese samples. The green color represents old and lean samples. The x-axis and y-axis represent the sample position vector coordinates according to the classifier.

PPI networks of feature genes and network topology analysis. (A) PPI network constructed based on the feature genes in the age group. (B) PPI network constructed based on the feature genes in the BMI group. The yellow color represents upregulated genes and the blue color represents downregulated genes. The young group serves as the control in the age group, and the lean group serves as the control in the BMI group. The pink color represents the remaining interactional genes, with the exception of the feature genes. PPI, protein-protein interaction; BMI, body mass index.

Node degree distribution of the feature genes. (A) Node degree distribution of the age group network. (B) Node degree distribution of the body mass index group network.

Feature genes in age and body mass index groups.

A, Age group
Gene P-value log2FC
S100A1 0.0004 −5.64
BARX2 0.0096 −2.95
IGFBP1 0.0341 −2.46
HAPLN1 0.0031 −2.21
DIXDC1 0.0123 −2.01
FGFBP2 0.0381 −1.92
FBLN7 0.0076 −1.81
VEGFA 0.0173 −1.66
GREM1 0.0171 −1.56
TF 0.0281 −1.40
BCL2 0.0345 −1.30
DCHS1 0.0020 −1.29
SMOC1 0.0042 −1.23
CAPN6 0.0167 −1.20
PTGES 0.0352 −1.08
VIT 0.0459 1.00
LYZ 0.0258 1.23
FGL2 0.0033 1.24
CCNB1 0.0030 1.31
PBK 0.0043 1.39
BST2 0.0102 1.46
CTSS 0.0222 1.48
HTRA4 0.0068 1.51
DDHD1 0.0093 1.55
TSPAN7 0.0007 1.57
CALCRL 0.0198 2.61
EMCN 0.0171 6.40
CD36 0.0010 9.33
B, Body mass index group
Gene P-value log2FC
EMCN 0.0169 −6.40
FRY 0.0353 −4.20
EDNRA 0.0084 −3.18
MFAP2 0.0304 −2.26
CXCL12 0.0113 −2.19
ABHD2 0.0010 −1.86
EYA4 0.0339 −1.79
ARHGAP11B 0.0145 −1.57
MKL2 0.0260 −1.41
MMP14 0.0398 −1.26
ISLR 0.0121 −1.18
DIXDC1 0.0150 1.03
BTC 0.0212 1.14
RCAN1 0.0191 1.34
TF 0.0325 1.45
CST6 0.0397 1.48
RAD51AP2 0.0007 1.58
KAL1 0.0098 2.40
S100A8 0.0042 2.74
IGFBP1 0.0157 2.95

FC, fold change.

List of top 10 highest degree nodes of the age and body mass index groups.

A, Age group
Gene P-value log2FC Degree
BCL2 0.0345 −1.30 74
CCNB1 0.0030 1.31 30
S100A1 0.0004 −5.64 18
VEGFA 0.0173 −1.66 18
TF 0.0281 −1.40 14
CD36 0.0010 9.33 13
PBK 0.0043 1.39 12
CTSS 0.0222 1.48 8
HAPLN1 0.0031 −2.21 7
CALCRL 0.0198 2.61 6
B, Body mass index group
Gene P-value log2FC Degree
MMP14 0.0398 −1.26 20
S100A8 0.0042 2.74 16
TF 0.0325 1.45 14
CXCL12 0.0113 −2.19 7
EDNRA 0.0084 −3.18 6
IGFBP1 0.0157 2.95 5
MFAP2 0.0304 −2.26 5
RCAN1 0.0191 1.34 5
BTC 0.0212 1.14 4
MKL2 0.0260 −1.41 3

Pathways associated with feature genes in the age and BMI groups.

Group Pathway ID P-value Genes
Age HIF-1 signaling pathway hsa04066 0.002013 VEGFA, TF, BCL2
BMI Mineral absorption hsa04978 0.043813 TF

HIF-1, hypoxia-inducible factor 1; VEGFA, vascular endothelial growth factor A; TF, transferrin; BCL2, B-cell lymphoma 2; BMI, body mass index.

Overlapping genes between candidate genes and differentially expressed genes in GSE51588.

Differentially expressed genes Candidate genes in GSE51588
Gene name P-value logFC Type P-value logFC
TF 2.81×10-2 −1.40 Age 3.20×10-3 −0.56
VEGFA 1.73×10-2 −1.66 Age 7.95×10-3 −0.62
BST2 1.02×10-2 1.46 Age 6.25×10-3 0.64
EMCN 1.71×10-2 6.40 Age 1.23×10-5 0.65
VIT 4.59×10-2 1.00 Age 4.17×10-5 0.70
ABHD2 9.91×10-4 −1.86 BMI 2.09×10-2 −0.56
ARHGAP11B 1.45×10-2 −1.57 BMI 1.19×10-2 −0.72
BTC 2.12×10-2 1.14 BMI 4.12×10-3 1.13
DIXDC1 1.50×10-2 1.03 BMI 2.27×10-7 0.63
KAL1 9.84×10-3 2.40 BMI 1.39×10-3 0.95
RAD51AP2 6.60×10-4 1.58 BMI 4.94×10-2 0.37

BMI, body mass index.