Objectives: To present an unbiased approach to identify positional transcript single nucleotide polymorphisms (SNPs) of osteoarthritis (OA) risk loci by allelic expression imbalance (AEI) analyses... Show moreObjectives: To present an unbiased approach to identify positional transcript single nucleotide polymorphisms (SNPs) of osteoarthritis (OA) risk loci by allelic expression imbalance (AEI) analyses using RNA sequencing of articular cartilage and subchondral bone from OA patients. Methods: RNA sequencing from 65 articular cartilage and 24 subchondral bone from OA patients was used for AEI analysis. AEI was determined for all genes present in the 100 regions reported by the genome-wide association studies (GWAS) catalog that were also expressed in cartilage or bone. The count fraction of the alternative allele (phi) was calculated for each heterozygous individual with the risk SNP or with the SNP in linkage disequilibrium (LD) with it (r(2) > 0.6). Furthermore, a meta-analysis was performed to generate a meta-phi (null hypothesis median phi = 0.49) and P-value for each SNP. Results: We identified 30 transcript SNPs (28 in cartilage and two in subchondral bone) subject to AEI in 29 genes. Notably, 10 transcript SNPs were located in genes not previously reported in the GWAS catalog, including two long intergenic non-coding RNAs (lincRNAs), MALAT1 (meta-phi = 0.54, FDR = 1.7x10(-4)) and ILF3-DT (meta-phi = 0.6, FDR = 1.75x10(-5)). Moreover, 12 drugs were interacting with seven genes displaying AEI, of which seven drugs have been already approved. Conclusions: By prioritizing proxy transcript SNPs that mark AEI in cartilage and/or subchondral bone at loci harbouring GWAS signals, we present an unbiased approach to identify the most likely functional OA risk-SNP and gene. We identified 10 new potential OA risk genes ready for further translation towards underlying biological mechanisms. Show less
Houtman, E.; Tuerlings, M.; Suchiman, H.E.D.; Lakenberg, N.; Cornelis, F.M.F.; Mei, H.L.; ... ; Meulenbelt, I. 2022
Objectives To investigate whether the deiodinase inhibitor iopanoic acid (IOP) has chondroprotective properties, a mechanical stress induced model of human aged explants was used to test both... Show moreObjectives To investigate whether the deiodinase inhibitor iopanoic acid (IOP) has chondroprotective properties, a mechanical stress induced model of human aged explants was used to test both repeated dosing and slow release of IOP. Methods Human osteochondral explants subjected to injurious mechanical stress (65%MS) were treated with IOP or IOP encapsulated in poly lactic-co-glycolic acid-polyethylene glycol nanoparticles (NP-IOP). Changes to cartilage integrity and signalling were determined by Mankin scoring of histology, sulphated glycosaminoglycan (sGAG) release and expression levels of catabolic, anabolic and hypertrophic markers. Subsequently, on a subgroup of samples, RNA sequencing was performed on 65%MS (n = 14) and 65%MS+IOP (n = 7) treated cartilage to identify IOP's mode of action. Results Damage from injurious mechanical stress was confirmed by increased cartilage surface damage in the Mankin score, increased sGAG release, and consistent upregulation of catabolic markers and downregulation of anabolic markers. IOP and, though less effective, NP-IOP treatment, reduced MMP13 and increased COL2A1 expression. In line with this, IOP and NP-IOP reduced cartilage surface damage induced by 65%MS, while only IOP reduced sGAG release from explants subjected to 65%MS. Lastly, differential expression analysis identified 12 genes in IOP's mode of action to be mainly involved in reducing metabolic processes (INSIG1, DHCR7, FADS1 and ACAT2) and proliferation and differentiation (CTGF, BMP5 and FOXM1). Conclusion Treatment with the deiodinase inhibitor IOP reduced detrimental changes of injurious mechanical stress. In addition, we identified that its mode of action was likely on metabolic processes, cell proliferation and differentiation. Show less
Ramos, Y.F.M.; Almeida, R.C. de; Lakenberg, N.; Suchiman, E.; Mei, H.L.; Kloppenburg, M.; ... ; Meulenbelt, I. 2021
Objective: To identify and validate circulating micro RNAs (miRNAs) that mark gene expression changes in articular cartilage early in osteoarthritis (OA) pathophysiology process. Methods: Within... Show moreObjective: To identify and validate circulating micro RNAs (miRNAs) that mark gene expression changes in articular cartilage early in osteoarthritis (OA) pathophysiology process. Methods: Within the ongoing RAAK study, human preserved OA cartilage and plasma (N = 22 paired samples) was collected for RNA sequencing (respectively mRNA and miRNA). Spearman correlation was determined for 114 cartilage genes consistently and significantly differentially expressed early in osteoarthritis and 384 plasma miRNAs. Subsequently, the minimal number of circulating miRNAs serving to discriminate between progressors and non-progressors was assessed by regression analysis and area under receiver operating curves (AUC) was calculated with progression data and plasma miRNA sequencing from the GARP study (N = 71). Results: We identified strong correlations (rho >= |0.7|) among expression levels of 34 unique plasma miRNAs and 21 genes, including 4 genes that correlated with multiple miRNAs. The strongest correlation was between let-7d-5p and EGFLAM (rho = -0.75, P = 6.9 x 10(-5)). Regression analysis of the 34 miRNAs resulted in a set of 7 miRNAs that, when applied to the GARP study, demonstrated clinically relevant predictive value with AUC > 0.8 for OA progression over 2 years and near-clinical value for progression over 5 years- (AUC = 0.8). Conclusions: We show that plasma miRNAs levels reflect gene expression levels in cartilage and can be exploited to represent ongoing pathophysiological processes in articular cartilage. We advocate that identified signature of 7 plasma miRNAs can contribute to direct further studies toward early biomarkers predictive for progression of osteoarthritis over 2 and 5 years. Show less
Tuerlings, M.; Hoolwerff, M. van; Houtman, E.; Suchiman, E.H.E.D.; Lakenberg, N.; Mei, H.L.; ... ; Meulenbelt, I. 2021
Objective To identify key determinants of the interactive pathophysiologic processes in subchondral bone and cartilage in osteoarthritis (OA).Methods We performed RNA sequencing on macroscopically... Show moreObjective To identify key determinants of the interactive pathophysiologic processes in subchondral bone and cartilage in osteoarthritis (OA).Methods We performed RNA sequencing on macroscopically preserved and lesional OA subchondral bone from patients in the Research Arthritis and Articular Cartilage study who underwent joint replacement surgery due to OA (n = 24 sample pairs: 6 hips and 18 knees). Unsupervised hierarchical clustering and differential expression analyses were conducted. Results were combined with data on previously identified differentially expressed genes in cartilage (partly overlapping samples) as well as data on recently identified OA risk genes.Results We identified 1,569 genes that were significantly differentially expressed between lesional and preserved subchondral bone, including CNTNAP2 (fold change [FC] 2.4, false discovery rate [FDR] 3.36 x 10(-5)) and STMN2 (FC 9.6, FDR 2.36 x 10(-3)). Among these 1,569 genes, 305 were also differentially expressed, and with the same direction of effect, in cartilage, including the recently recognized OA susceptibility genes IL11 and CHADL. Upon differential expression analysis with stratification for joint site, we identified 509 genes that were exclusively differentially expressed in subchondral bone of the knee, including KLF11 and WNT4. These genes that were differentially expressed exclusively in the knee were enriched for involvement in epigenetic processes, characterized by, e.g., HIST1H3J and HIST1H3H.Conclusion IL11 and CHADL were among the most consistently differentially expressed genes OA pathophysiology-related genes in both bone and cartilage. As these genes were recently also identified as robust OA risk genes, they classify as attractive therapeutic targets acting on 2 OA-relevant tissues. Show less
Almeida, R.C. de; Mahfouz, A.; Mei, H.L.; Houtman, E.; Hollander, W. den; Soul, J.; ... ; Meulenbelt, I. 2021
Objective. To identify OA subtypes based on cartilage transcriptomic data in cartilage tissue and characterize their underlying pathophysiological processes and/or clinically relevant... Show moreObjective. To identify OA subtypes based on cartilage transcriptomic data in cartilage tissue and characterize their underlying pathophysiological processes and/or clinically relevant characteristics.Methods. This study includes n = 66 primary OA patients (41 knees and 25 hips), who underwent a joint replacement surgery, from which macroscopically unaffected (preserved, n = 56) and lesioned (n = 45) OA articular cartilage were collected [Research Arthritis and Articular Cartilage (RAAK) study]. Unsupervised hierarchical clustering analysis on preserved cartilage transcriptome followed by clinical data integration was performed. Protein-protein interaction (PPI) followed by pathway enrichment analysis were done for genes significant differentially expressed between subgroups with interactions in the PPI network.Results. Analysis of preserved samples (n = 56) resulted in two OA subtypes with n = 41 (cluster A) and n = 15 (cluster B) patients. The transcriptomic profile of cluster B cartilage, relative to cluster A (DE-AB genes) showed among others a pronounced upregulation of multiple genes involved in chemokine pathways. Nevertheless, upon investigating the OA pathophysiology in cluster B patients as reflected by differentially expressed genes between preserved and lesioned OA cartilage (DE-OA-B genes), the chemokine genes were significantly downregulated with OA pathophysiology. Upon integrating radiographic OA data, we showed that the OA phenotype among cluster B patients, relative to cluster A, may be characterized by higher joint space narrowing (JSN) scores and low osteophyte (OP) scores.Conclusion. Based on whole-transcriptome profiling, we identified two robust OA subtypes characterized by unique OA, pathophysiological processes in cartilage as well as a clinical phenotype. We advocate that further characterization, confirmation and clinical data integration is a prerequisite to allow for development of treatments towards personalized care with concurrently more effective treatment response. Show less
Houtman, E.; Hoolwerff, M. van; Lakenberg, N.; Suchiman, E.H.D.; Zwaag, E.V.V. van der van der; Nelissen, R.G.H.H.; ... ; Meulenbelt, I. 2021
Introduction: Likely due to ignored heterogeneity in disease pathophysiology, osteoarthritis (OA) has become the most common disabling joint disease, without effective disease-modifying treatment... Show moreIntroduction: Likely due to ignored heterogeneity in disease pathophysiology, osteoarthritis (OA) has become the most common disabling joint disease, without effective disease-modifying treatment causing a large social and economic burden. In this study we set out to explore responses of aged human osteochondral explants upon different OA-related perturbing triggers (inflammation, hypertrophy and mechanical stress) for future tailored biomimetic human models.Methods: Human osteochondral explants were treated with IL-1 beta (10 ng/ml) or triiodothyronine (T3; 10 nM) or received 65% strains of mechanical stress (65% MS). Changes in chondrocyte signalling were determined by expression levels of nine genes involved in catabolism, anabolism and hypertrophy. Breakdown of cartilage was measured by sulphated glycosaminoglycans (sGAGs) release, scoring histological changes (Mankin score) and mechanical properties of cartilage.Results: All three perturbations (IL-1 beta, T3 and 65% MS) resulted in upregulation of the catabolic genes MMP13 and EPAS1. IL-1 beta abolished COL2A1 and ACAN gene expression and increased cartilage degeneration, reflected by increased Mankin scores and sGAGs released. Treatment with T3 resulted in a high and significant upregulation of the hypertrophic markers COL1A1, COL10A1 and ALPL. However, 65% MS increased sGAG release and detrimentally altered mechanical properties of cartilage.Conclusion: We present consistent and specific output on three different triggers of OA. Perturbation with the pro-inflammatory IL-1 beta mainly induced catabolic chondrocyte signalling and cartilage breakdown, while T3 initiated expression of hypertrophic and mineralization markers. Mechanical stress at a strain of 65% induced catabolic chondrocyte signalling and changed cartilage matrix integrity. The major strength of our ex vivo models was that they considered aged, preserved, human cartilage of a heterogeneous OA patient population. As a result, the explants may reflect a reliable biomimetic model prone to OA onset allowing for development of different treatment modalities. Show less
Hoolwerff, M. van; Metselaar, P.I.; Tuerlings, M.; Suchiman, H.E.D.; Lakenberg, N.; Ramos, Y.F.M.; ... ; Meulenbelt, I. 2020
Objective To identify robustly differentially expressed long noncoding RNAs (lncRNAs) with osteoarthritis (OA) pathophysiology in cartilage and to explore potential target messenger RNA (mRNA) by... Show moreObjective To identify robustly differentially expressed long noncoding RNAs (lncRNAs) with osteoarthritis (OA) pathophysiology in cartilage and to explore potential target messenger RNA (mRNA) by establishing coexpression networks, followed by functional validation. Methods RNA sequencing was performed on macroscopically lesioned and preserved OA cartilage from patients who underwent joint replacement surgery due to OA (n = 98). Differential expression analysis was performed on lncRNAs that were annotated in GENCODE and Ensembl databases. To identify potential interactions, correlations were calculated between the identified differentially expressed lncRNAs and the previously reported differentially expressed protein-coding genes in the same samples. Modulation of chondrocyte lncRNA expression was achieved using locked nucleic acid GapmeRs. Results By applying our in-house pipeline, we identified 5,053 lncRNAs that were robustly expressed, of which 191 were significantly differentially expressed (according to false discovery rate) between lesioned and preserved OA cartilage. Upon integrating mRNA sequencing data, we showed that intergenic and antisense differentially expressed lncRNAs demonstrate high, positive correlations with their respective flanking sense genes. To functionally validate this observation, we selectedP3H2-AS1, which was down-regulated in primary chondrocytes, resulting in the down-regulation ofP3H2gene expression levels. As such, we can confirm thatP3H2-AS1regulates its sense geneP3H2. Conclusion By applying an improved detection strategy, robustly differentially expressed lncRNAs in OA cartilage were detected. Integration of these lncRNAs with differential mRNA expression levels in the same samples provided insight into their regulatory networks. Our data indicates that intergenic and antisense lncRNAs play an important role in regulating the pathophysiology of OA. Show less
Raz, Y.; Akker, E.B. van den; Roest, T.; Riaz, M.; Rest, O. van de; Suchiman, H.E.D.; ... ; Slagboom, P.E. 2020
Skeletal muscles control posture, mobility and strength, and influence whole-body metabolism. Muscles are built of different types of myofibers, each having specific metabolic, molecular, and... Show moreSkeletal muscles control posture, mobility and strength, and influence whole-body metabolism. Muscles are built of different types of myofibers, each having specific metabolic, molecular, and contractile properties. Fiber classification is, therefore, regarded the key for understanding muscle biology, (patho-) physiology. The expression of three myosin heavy chain (MyHC) isoforms, MyHC-1, MyHC-2A, and MyHC-2X, marks myofibers in humans. Typically, myofiber classification is performed by an eye-based histological analysis. This classical approach is insufficient to capture complex fiber classes, expressing more than one MyHC-isoform. We, therefore, developed a methodological procedure for high-throughput characterization of myofibers on the basis of multiple isoforms. The mean fluorescence intensity of the three most abundant MyHC isoforms was measured per myofiber in muscle biopsies of 56 healthy elderly adults, and myofiber classes were identified using computational biology tools. Unsupervised clustering revealed the existence of six distinct myofiber clusters. A comparison with the visual assessment of myofibers using the same images showed that some of these myofiber clusters could not be detected or were frequently misclassified. The presence of these six clusters was reinforced by RNA expressions levels of sarcomeric genes. In addition, one of the clusters, expressing all three MyHC isoforms, correlated with histological measures of muscle health. To conclude, this methodological procedure enables deep characterization of the complex muscle heterogeneity. This study opens opportunities to further investigate myofiber composition in comparative studies. Show less
Insights into individual differences in gene expression and its heritability (h(2)) can help in understanding pathways from DNA to phenotype. We estimated the heritability of gene expression of 52... Show moreInsights into individual differences in gene expression and its heritability (h(2)) can help in understanding pathways from DNA to phenotype. We estimated the heritability of gene expression of 52,844 genes measured in whole blood in the largest twin RNA-Seq sample to date (1497 individuals including 459 monozygotic twin pairs and 150 dizygotic twin pairs) from classical twin modeling and identity-by-state-based approaches. We estimated for each gene h(total)(2), composed of cis-heritability (h(cis)(2), the variance explained by single nucleotide polymorphisms in the cis-window of the gene), and trans-heritability (h(res)(2), the residual variance explained by all other genome-wide variants). Mean h(total)(2) was 0.26, which was significantly higher than heritability estimates earlier found in a microarray-based study using largely overlapping (>60%) RNA samples (mean h(2) = 0.14, p = 6.15 x 10(-258)). Mean h(cis)(2) was 0.06 and strongly correlated with beta of the top cis expression quantitative loci (eQTL, rho = 0.76, p < 10(-308)) and with estimates from earlier RNA-Seq-based studies. Mean h(res)(2) was 0.20 and correlated with the beta of the corresponding trans-eQTL (rho = 0.04, p < 1.89 x 10(-3)) and was significantly higher for genes involved in cytokine-cytokine interactions (p = 4.22 x 10(-15)), many other immune system pathways, and genes identified in genome-wide association studies for various traits including behavioral disorders and cancer. This study provides a thorough characterization of cis- and trans-h(2) estimates of gene expression, which is of value for interpretation of GWAS and gene expression studies. Show less
X-inactivation is a well-established dosage compensation mechanism ensuring that X-chromosomal genes are expressed at comparable levels in males and females. Skewed X-inactivation is often... Show moreX-inactivation is a well-established dosage compensation mechanism ensuring that X-chromosomal genes are expressed at comparable levels in males and females. Skewed X-inactivation is often explained by negative selection of one of the alleles. We demonstrate that imbalanced expression of the paternal and maternal X-chromosomes is common in the general population and that the random nature of the X-inactivation mechanism can be sufficient to explain the imbalance. To this end, we analyzed blood-derived RNA and whole-genome sequencing data from 79 female children and their parents from the Genome of the Netherlands project. We calculated the median ratio of the paternal over total counts at all X-chromosomal heterozygous single-nucleotide variants with coverage ≥10. We identified two individuals where the same X-chromosome was inactivated in all cells. Imbalanced expression of the two X-chromosomes (ratios ≤0.35 or ≥0.65) was observed in nearly 50% of the population. The empirically observed skewing is explained by a theoretical model where X-inactivation takes place in an embryonic stage in which eight cells give rise to the hematopoietic compartment. Genes escaping X-inactivation are expressed from both alleles and therefore demonstrate less skewing than inactivated genes. Using this characteristic, we identified three novel escapee genes (SSR4, REPS2, and SEPT6), but did not find support for many previously reported escapee genes in blood. Our collective data suggest that skewed X-inactivation is common in the general population. This may contribute to manifestation of symptoms in carriers of recessive X-linked disorders. We recommend that X-inactivation results should not be used lightly in the interpretation of X-linked variants. Show less
Bonder, M.J.; Luijk, R.; Zhernakova, D.V.; Moed, M.; Deelen, P.; Vermaat, M.; ... ; BIOS Consortium 2017
