Assessment of DNA Methylation Patterns in the Bone and Cartilage of a Nonhuman Primate Model of Osteoarthritis

Abstract

Objective

Osteoarthritis (OA) affects humans and several other animals. Thus, the mechanisms underlying this disorder, such as specific skeletal tissue DNA methylation patterns, may be evolutionary conserved. However, associations between methylation and OA have not been readily studied in nonhuman animals. Baboons serve as important models of disease and develop OA at rates similar to those in humans. Therefore, this study investigated the associations between methylation and OA in baboons to advance the evolutionary understanding of OA.

Design

Trabecular bone and cartilage was collected from the medial condyles of adult female baboon femora, 5 with and 5 without knee OA. The Infinium HumanMethylation450 BeadChip (450K array) was used to identify DNA methylation patterns in these tissues.

Results

Approximately 44% of the 450K array probes reliably align to the baboon genome, contain a CpG site of interest, and maintain a wide distribution throughout the genome. Of the 2 filtering methods tested, both identified significantly differentially methylated positions (DMPs) between healthy and OA individuals in cartilage tissues, and some of these patterns overlap with those previously identified in humans. Conversely, no DMPs were found between tissue types or between disease states in bone tissues.

Conclusions

Overall, the 450K array can be used to measure genome-wide DNA methylation in baboon tissues and identify significant associations with complex traits. The results of this study indicate that some DNA methylation patterns associated with OA are evolutionarily conserved, while others are not. This warrants further investigation in a larger and more phylogenetically diverse sample set.

Keywords

osteoarthritis animal model articular cartilage bone DNA methylation

Introduction

Osteoarthritis (OA) is a complex degenerative joint disease, and OA of the knee is one of the leading causes of disability across the globe.¹ Thus, research endeavors to describe the molecular mechanisms that contribute to this disorder are underway. Both genetic and environmental factors have some effect,^2,3 as well as epigenetic factors that bridge the gap between genetics and the environment. In particular, DNA methylation, which regulates gene expression, is thought to play an influential role in the development of degenerative skeletal disorders like OA.^4-10 Animal models, such as mice, rats, rabbits, guinea pigs, dogs, sheep, goats, and horses, have been essential in discerning some of the processes inherent to OA development.^11,12 Nevertheless, all these animals are limited in their abilities to fully inform our understanding of human OA, so the search to find a gold standard animal model for OA is still ongoing.¹³ Finally, although variation in skeletal tissue DNA methylation patterns are thought to be involved in the development and progression of OA,^4-10 this epigenetic mechanism has not been readily studied in animal models. In part, this may be due to the compatibility constraints of DNA methylation assays that are designed specifically for humans, as well as the limited efforts to optimize these assays for nonhuman animals.^14,15

Nonhuman primates can serve as important models of disease for humans because they are phylogenetically close to humans. Baboons (Papio sp.) are a particularly good model of disease, especially OA,¹⁶ as they naturally develop OA at rates similar to those observed in humans.^3,16 Additionally, because of their evolutionary proximity to humans, further investigation of the molecular processes innate to OA development and progression in baboons as compared to these mechanisms in humans will advance the evolutionary understanding of this disease. Finally, the relative genetic conservation between baboons and humans makes the optimization and use of standardized DNA methylation assays possible. Specifically, the Infinium HumanMethylation450 BeadChip (450K array), which is a cost-efficient application for assessing genome-wide DNA methylation patterns in humans, has been successfully used for some nonhuman primate species. These and other nonhuman primate DNA methylation studies have primarily used DNA extracted from blood or other soft tissues.^14,15,17-20 However, this technique has not yet been used to study DNA methylation variation in baboon skeletal tissues or how it relates to the development of OA in a nonhuman primate species.

For this study, we used the 450K array to identify DNA methylation patterns in femur bone and cartilage of age-matched female baboons, 5 with and 5 without knee OA. This was done to validate that the 450K array can be used for nonhuman primate skeletal tissue DNA extracts. Additionally, this study was performed to determine if DNA methylation variation is associated with OA in baboons and in a manner similar to that observed in humans.

Methods

Ethics Statement

Nonhuman primate tissue samples included were opportunistically collected at routine necropsy of these animals. No animals were sacrificed for this study, and no living animals were used in this study.

Samples

Baboon (Papio sp.) samples come from captive colonies at the Southwest National Primate Research Center in Texas. These samples are ideal because many environmental factors that influence skeletal development and maintenance (e.g., diet and exposure to sunlight, which influences vitamin D production) are controlled and consistent across individuals. Femora were opportunistically collected at routine necropsy of these animals and stored in −20°C freezers at the Texas Biomedical Research Institute after dissection. These preparation and storage conditions ensured the preservation of skeletal DNA methylation patterns. Samples include skeletally healthy adult baboons (n = 5) and adult baboons with severe OA (n = 5). Age ranges are comparable between each group (19.30 ± 1.70 and 19.24 ± 1.73 years, respectively), and only females were included in this study.

Assessment of Osteoarthritis

Classification of adult baboons as having healthy or OA knees was determined through visual examination of the distal femora and macroscopic inspection of the distal articular surface cartilage. Each specimen was assigned an OA severity score. Briefly, grade 1 is unaffected, grade 2 is mild OA as indicated by cartilage fibrillation, grade 3 is moderate OA as indicated by cartilage lesions, and grade 4 is advanced OA as indicated by eburnation.³ From this, binary classifications were made such that all healthy adult baboons have 100% grade 1 on one or both distal femora, and all OA adult baboons have a variable percentage of grades 3 or 4 on one or both distal femora ( Fig. 1 ).

Figure 1.

Representative examples of baboon knees (distal femora) that are healthy or have OA.

DNA Extraction

From the distal femoral condyles, cartilage scrapings were collected using scalpels and processed with a homogenizer, and trabecular bone samples from the interior of the condyle were obtained using a Dremel and pulverized into bone dust using a BioPulverizer. Tissues were collected from this region because OA legions in both humans and baboons generally appear first on the cartilage of the medial condyle. Both tissues are included in this project because human skeletal epigenetic studies are based on trabecular bone and cartilage, and it is important to standardize tissue type for comparative purposes. DNA was extracted from these processed tissues using a phenol-chloroform protocol optimized for skeletal tissues²¹ and quantified using both Nanodrop and Qubit machines (Supplementary Table S1, available in the online version of the article).

Genome-Wide DNA Methylation Profiling

Genome-wide DNA methylation was assessed using Infinium HumanMethylation450 BeadChips (450K array). These arrays analyze the methylation status of more than 485,000 sites throughout the genome, covering 99% of RefSeq genes and 96% of the UCSC-defined CpG islands and their flanking regions. For each sample, approximately 500 ng of genomic DNA (Supplementary Table S1, available in the online version of the article) was bisulfite converted using the EZ DNA Methylation Gold Kit, according to the manufacturer’s instructions (Zymo Research), with modifications described in the Infinium Methylation Assay Protocol. Following manufacturer guidelines (Illumina), this processed DNA was then whole-genome amplified, enzymatically fragmented, hybridized to the arrays, and imaged using the Illumina iScan system. The array data discussed in this publication have been deposited in NCBI’s Gene Expression Omnibus and are accessible through GEO Series accession number GSE101733.

Processing of Methylation Data

Raw fluorescent data were normalized to account for the noise inherent within and between the arrays themselves. Specifically, we performed a normal-exponential out-of-band (Noob) background correction method with dye-bias normalization²² to adjust for background fluorescence and dye-based biases and followed this with a between-array normalization method (functional normalization),²³ which removes unwanted variation by regressing out variability explained by the control probes present on the array as implemented in the minfi v1.20.2 package in R v3.3.1,²⁴ which is part of the Bioconductor project.²⁵ This method has been found to outperform other existing approaches for studies that compare conditions with known large-scale differences,²³ such as those assessed in this study.

After normalization, methylation values (β values) for each site were calculated as the ratio of methylated probe signal intensity to the sum of both methylated and unmethylated probe signal intensities. These β values range from 0 to 1 and represent the average methylation levels at each site across the entire population of cells from which DNA was extracted (0 = completely unmethylated sites, 1 = fully methylated sites). Every β value in the Infinium platform is accompanied by a detection P value, and those with failed detection levels (P > 0.05) in greater than 10% of samples were removed from downstream analyses. Because β values have high heteroscedasticity, they are not statistically valid for use in differential methylation analyses.²⁶ Thus, M values were calculated as the log transformed ratio of methylated signal to unmethylated signal and used in the statistical analyses described below.

The probes on the arrays were designed to hybridize specifically with human DNA, so our use of nonhuman primate DNA required that probes nonspecific in the baboon genome, which could produce biased methylation measurements, be computationally filtered out and excluded from downstream analyses. This was accomplished using 2 different methods modified from Hernando-Herraez et al.¹⁴ and Ong et al.¹⁵ For both methods, we used blastn²⁷ to map the 485,512 50bp probes onto the Papio anubis genome (Assembly: Panu_2.0, Accession: GCF_000264685.2) using an e-value threshold of e⁻¹⁰. We retained probes that successfully mapped to the baboon genome, had only one unique BLAST hit, and targeted CpG sites (Supplemental Probe Annotation File, available in the online version of the article). These 213,858 probes that aligned to the baboon genome were then further filtered using 1 of 2 methods.

The first method used criteria based on sequence alignment,¹⁴ and the second method used criteria based on gene symbol similarities.¹⁵ For the first method (alignment filter), we only retained probes that had 0 mismatches in the 5 bp closest to and including the CpG site and that had 0 to 2 mismatches in the 45 bp not including the CpG site.¹⁴ This more stringent sequence similarity filtering retained 133,264 probes. For the second method (gene symbol filter), we identified the closest baboon gene to each probe alignment site and checked for corresponding gene name matches between humans and baboons.¹⁵ For baboons, this information was obtained from GFF and Ensembl BioMart data. Only those probes with partial or complete gene matches were retained. This more lenient filtering retained 130,307 probes. Additionally, cross-reactive probes,²⁸ probes containing single nucleotide polymorphisms (SNPs) at the CpG site, probes detecting SNP information, probes detecting methylation at non-CpG sites, and probes targeting sites within the sex chromosomes were removed using the minfi v1.20.2 package in R v3.3.1.²⁴ This filtering retained a final set of 120,305 probes for the alignment filter criteria and 112,760 probes for the gene symbol criteria ( Fig. 2 ).

Figure 2.

(A and B) Density plots of β values after normalization and probe filtering using the alignment filter criteria (A) or the gene symbol filter criteria (B). (C and D) Multidimensional scaling plots showing the first 2 principal components that describe genome-wide methylation variation after normalization and filtering using the alignment filter criteria (C) or the gene symbol filter criteria (D). Each point represents one sample that is either from healthy bone, healthy cartilage, OA bone, or OA cartilage. In the multidimensional scaling plots, these categories do not form distinct clusters.

Statistical Analysis of Differential Methylation

To identify sites that were significantly differentially methylated across comparative groups, we designed and tested generalized linear mixed models (GLMMs) which related the variables of interest to the DNA methylation patterns for each site, while accounting for latent variables.²⁹ Sites found to have significant associations were classified as significant differentially methylated positions (DMPs). Specifically, a GLMM was used to estimate differences in methylation levels for each of the following contrasts: (1) between bone and cartilage in OA baboons, (2) between bone and cartilage in healthy baboons, (3) between OA and healthy baboon bone, (4) between OA and healthy baboon cartilage, (5) among all 4 combinations of tissue type and disease state (healthy bone vs. healthy cartilage vs. OA bone vs. OA cartilage).

Additional variables included in this GLMM were unknown latent variables calculated using the iteratively reweighted least squares approach in the sva v3.22.0 package in R v3.3.1.^30,31 The four latent variables estimated were included to help mitigate any unknown batch and cell heterogeneity effects on methylation variation at each site. No predefined batch effects for the arrays were included because these did not appear to have large effects on the data (Supplementary Figure S1, available in the online version of the article). Alternative methods to account for cell heterogeneity exist, but they are specific to whole blood,^30,32 require reference epigenetic data, or are reference free methods³³ that are comparable to the sva method.³⁴ Out of the known cell types in skeletal tissues, only chondrocytes and osteoblasts have reference epigenomes available on the International Human Epigenomics Consortium, and these are only for humans, not nonhuman primates. Thus, because no standard method is available to correct for the heterogenous cell structure in nonhuman primate skeletal tissue, we chose the described sva method.

The GLMM design matrix was fit to the M value array data by generalized least squares using the limma v3.30.13 package in R v3.3.1,^25,35 and the estimated coefficients and standard errors for the defined tissue type and disease status contrasts were computed. Finally, for each coefficient, an empirical Bayes approach³⁶ was used to compute moderated t statistics, log-odds ratios of differential methylation, and associated P values adjusted for multiple testing.³⁷ Significant DMPs for the effect of tissue type and disease status contrasts were defined as those having log fold changes in M values corresponding to an adjusted P value of less than 0.05.

Results

The aim of this study was to use the 450K array to identify DNA methylation patterns in femur bone and cartilage of age-matched female baboons, 5 with and 5 without knee OA. To do this, we first assessed the effectiveness of the 450K array in identifying DNA methylation patterns in baboon DNA and of different probe filtering methods.

Alignment of 450K Array Probes with the Baboon Genome

Probes from the 450K array were aligned to the baboon genome (Supplemental Probe Annotation File, available in the online version of the article).^14,15 Out of the 485,512 50bp probes on the array, 213,858 probes (44%) map to the baboon genome with e-values less than e⁻¹⁰, have only unique BLAST hits, and target a CpG site ( Fig. 3A ). Out of these reliably mapped probes, 133,264 probes (62%) were retained after the alignment filter criteria ( Fig. 3B ), while 130,307 probes (61%) were retained after the gene symbol filter criteria ( Fig. 3C ). A total of 83,142 probes overlapped between both filtering methods (62% for the alignment filter criteria and 64% for the gene symbol filter criteria, Supplementary Figure S2, available in the online version of the article). Probes that reliably mapped to the baboon genome, that met the alignment filter criteria, or that met the gene symbol criteria covered approximately 18,800 genes with an average coverage of nine, six, or eight probes per gene, respectively (Supplementary Table S2, available in the online version of the article). Additionally, the retained probes covered a range of locations with respect to genes and CpG islands (Supplementary Table S2, available in the online version of the article), indicating that these filtered probes maintain a wide distribution throughout the genome.

Figure 3.

Filtering effects on the 450K array probes for baboons. (A) Pie chart showing the percent of 450K array probes that map to the baboon (Papio anubis) genome with e-values less than e⁻¹⁰, have only unique BLAST hits, and target a CpG site. Out of 485,512 probes total, 213,858 probes (44%) meet these criteria. (B) Pie chart showing the percent of probes, out of those that successfully mapped to the baboon genome, that contain 0 mismatches in 5 bp of the probe by and including the targeted CpG site and 0 to 2 mismatches in 45 bp of the probe not including the CpG site. Out of the 213,858 mapped probes, 133,264 probes (62%) meet these criteria. (C) Pie chart showing the percent of probes, out of those that successfully mapped to the baboon genome, with gene symbol matches to humans. Out of the 213,858 mapped probes, 130,307 probes (61%) meet these criteria.

Effectiveness of 450K Array Probes Using Baboon DNA

To determine how effectively the 450K array probes measured DNA methylation in baboon DNA, we performed Spearman correlation tests between the hybridization efficiency of each probe and parameters defining the alignment quality of each probe to the baboon genome. Specifically, both probe alignment bitscores and percent identity were significantly negatively correlated with probe hybridization efficiency, and probe alignment e-values were significantly positively correlated with probe hybridization efficiency, regardless of filtering criteria (Supplementary Table S3, available in the online version of the article). However, filtering probes using the alignment filter criteria retained proportionally more successfully hybridized probes than filtering probes using the gene symbol filter criteria (Supplementary Figure S3, available in the online version of the article). Thus, filtering probes using the alignment filter criteria likely produces more reliable results.

Differential Methylation and Osteoarthritis

Significant DMPs were only identified between healthy and OA individuals in cartilage tissues (Supplementary Table S4, available in the online version of the article). All these DMPs displayed decreased methylation in OA cartilage samples as compared to healthy cartilage samples, and some of these patterns overlapped with those previously identified in humans. Conversely, no DMPs were found between tissue types or between disease states in bone tissues. When filtering probes using the alignment filer criteria, 6 significant DMPs were identified between OA cartilage samples and healthy cartilage samples, while only two DMPs were identified when filtering probes using the gene symbol criteria ( Table 1 ). One locus, RFXAP, matched between these filtering methods. RUNX1 has previously been found to be differentially methylated in OA and healthy cartilage in humans, with OA cartilage having lower methylation as compared to healthy cartilage.⁵ The other genes associated with these probes (KLHL26, RFXAP, MIR497, MIR195, ELF1, ACSL1, and CMIP) have not previously been associated with OA-related differential methylation in humans.^4-10

Table 1.

Significant DMPs Between OA Cartilage and Healthy Cartilage.^a

	Probe ID	ΔM	P	Human Gene Symbol	Human CpG Position	Human Gene Group	Human CpG Island	Relation to CpG Island	Other Human OA Genes	Baboon Gene Symbol	Baboon CpG Position
Alignment Filter Probes	cg05295841	−3.6136	0.0247	KLHL26	chr19:18771126	Body	chr19:18771078-18771340	Island	CRTC1; CRLF1; COMP	CRTC1	chr19:17216922
	cg17890983	−1.9683	0.0247	RFXAP	chr13:37394336	Body	chr13:37393213-37394092	South Shore	SMAD9	RFXAP	chr17:15599705
	cg02329670	−1.4981	0.0411	MIR497; MIR195	chr17:6921400	TSS200; TSS1500	chr17:6925400-6927527	North Shelf	ALOX12	LOC103878622	chr16:6672282
	cg18456803	−2.0478	0.0411	ELF1	chr13:41593519	TSS200	NA	NA	NA	WBP4; LOC103879193	chr17:19630529
	cg13030790	−3.5437	0.0411	RUNX1	chr21:36421503	5’UTR; 1stExon	NA	NA	NA	RUNX1	chr3:11500594
	cg24721647	−2.1838	0.0411	ACSL1	chr4:185726836	5’UTR	chr4:185724434-185724647	South Shelf	NA	ACSL1	chr5:173909497
Gene Symbol Filter Probes	cg17890983	−1.9679	0.0374	RFXAP	chr13:37394336	Body	chr13:37393213-37394092	South Shore	SMAD9	RFXAP	chr17:15599705
Gene Symbol Filter Probes	cg04759112	−1.6320	0.0381	CMIP	chr16:81663392	Body; Body	NA	NA	NA	CMIP	chr20:63369656

Provided information includes the 450K array probe ID for each DMP (Probe ID), the log fold change in M values (ΔM), the P values adjusted for multiple testing (P), and additional gene details for each DMP. Results are shown for probes filtered using the alignment criteria and probes filtered using the gene symbol criteria. For all of these DMPs, OA cartilage samples are hypomethylated as compared with healthy cartilage samples. Additional gene details are based on the GRCh37 human reference genome and data available on the UCSC Genome Browser, as well as the Panu_2.0 baboon reference genome. Other OA genes include those genes in close proximity to the DMP that have previously been associated with OA. CRTC1 (Gene ID: 23373) is found about 23kb away from cg0529584, and mechanical stimuli has been shown to promote the nuclear translocation of this gene product into chondrocytes, which in turn regulates downstream gene expression.³⁸ Additionally, CRLF1 (Gene ID: 9244) is found about 53 kb away from cg05295841, is regulated by transforming growth factor–β1 (TGF-β1), and is highly expressed in OA cartilage.³⁹ COMP (Gene ID: 1311) is found about 122 kb from cg05295841, and its level of expression, which is influenced by several polymorphisms, has been associated with OA.⁴⁰ Approximately 25kb away from cg17890983, SMAD9 (Gene ID: 4093) is involved in TGF-β and BMP signaling, which themselves are involved in regulating bone and cartilage development.^41,42 Finally, ALOX12 (Gene ID: 239) is found about 7 kb away from cg02329670, and although polymorphisms in this gene are not associated with bone density, they are associated with fat accumulation,⁴³ which is associated with OA.⁴⁴

Discussion

Here we used the 450K array to identify DNA methylation variation in bone and cartilage tissues from a baboon model of OA. This was done both to determine the effectiveness of this application for baboon DNA and to assess the evolutionary conservation of epigenetic-OA associations in the primate lineage. We show that using the 450K array is feasible in baboon tissues. In silico probe filtering methods^14,15 indicated that 44% of all human probes could be reliably mapped to the baboon genome and contained an informative CpG. This number was lower than expected since previous researchers were able to reliably use these same methods to map 61% of the human probes to the Cynomologus macaque genome,¹⁵ another Old World monkey that is a close phylogenetic relative to baboons. This discrepancy in number may be due to the quality of each nonhuman primates’ genome assembly. While both are well annotated, the average scaffold length (88,649,475) and contig length (86,040) of the macaque genome (Assembly: Macaca_fascicularis_5.0, Accession: GCF_000264685.2) are higher than those (528,927 and 40,262, respectively) of the baboon genome.

Subsequent in silico analyses based on sequence alignment criteria¹⁴ and based on gene symbol criteria¹⁵ retained similar numbers of probes ( Fig. 3 ) that maintained wide and comparable distributions throughout the genome (Table S2). However, only a little more than half of the resulting probes for each filtering technique overlapped with one another (Figure S2). This discrepancy is likely due to the incomplete nature of the baboon genome annotation. More than half of the probes that fit the alignment filter criteria but not the gene symbol criteria (28,699 out of 50,117) are associated with generic gene symbol identifiers (LOC) to indicate the as of yet unknown functions of these regions. Conversely, all of the probes that fit the gene symbol criteria but not the alignment filter criteria have more than 3 mismatches with the baboon genome on average and have a maximum of 9 mismatches with the baboon genome. These high mismatch numbers are a potential concern for proper and accurate probe and baboon DNA hybridization.

Fittingly, after applying the 450K array to measure DNA methylation patterns of genomic material extracted from baboon skeletal tissues, we found that the hybridization efficiency of probes was significantly correlated with the alignment quality of each probe to the baboon genome, and thus, the degree of sequence conservation. The majority of filtered probes for both in silico methods passed quality controls and produced robust signals on the array, indicating that either filtering technique may be appropriate for future research. However, because the filtering method based on the alignment filter criteria retained a larger proportion of successfully hybridized probes than the method based on the gene symbol criteria (Figure S3) and because this method is less influenced by the degree of genome assembly annotation, we recommend that this alignment filter criteria method be preferentially used in subsequent nonhuman primate studies.

This work is an extension of previous work using the 450K array to study DNA methylation patterns in the tissues of nonhuman primates.^14,15 The 450K array is advantageous because it is cost-efficient per sample and simultaneously measures a large number of CpG loci with a broad genomic representation.⁴⁵ Similar to this study, previous researchers have applied the 450K array to measure DNA methylation patterns in great apes,¹⁴ which are closer to humans evolutionarily than baboons, and in macaques,¹⁵ which are comparable in proximity to humans evolutionarily as compared to baboons. In this study, we used a baboon model of OA to assess the evolutionary conservation of epigenetic-OA associations in the primate lineage, and we identified significant DMPs between healthy and OA individuals in cartilage tissues. All these loci showed hypomethylation in OA cartilage samples as compared with healthy cartilage samples. Six DMPs were identified when using the alignment filter criteria, and 2 DMPs were identified when using the gene symbol filter criteria. All together, these eight DMPs are annotated to eight genes that have a variety of functions and are also in close proximity to 5 other genes that have previously been associated with OA ( Table 1 ).

Some of these annotated genes have functions related to skeletal development and maintenance. For instance, RUNX1 (Gene ID: 861), also known as runt related transcription factor 1, is involved in the regulation of bone and cartilage cell development and differentiation.⁴⁶ Additionally, MIR497 (Gene ID: 574456) and MIR195 (Gene ID: 406971) are noncoding microRNAs that are involved in posttranscriptional regulation.⁴⁷ While both of these microRNAs have roles in the development of cancer,^48,49 they also play important regulatory roles in the differentiation of mesenchymal stromal/stem cells into bone related cells.⁵⁰

Other genes have functions related to the immune system which may have proximal roles in the development of OA. In particular, RFXAP (Gene ID: 5994), also known as regulatory factor X–associated protein, codes for a protein that assists in the transcriptional activation of major histocompatibility class II genes which are critical for the development and control of the immune system.⁵¹ Additionally, CMIP (Gene ID: 80790), also known as c-Maf inducing protein, codes for a protein that is involved in the T-cell signaling pathway, and SNPs within this gene have been associated with chronic diseases like diabetes.⁵²

The remaining genes do not have functions related to skeletal phenotypes, which makes their involvement in OA less clear. For example, KLHL26 (Gene ID: 55295), also known as kelch like family member 26, is part of a family of proteins that may be involved in protein ubiquitination.⁵³ Additionally, ACSL1 (Gene ID: 2180), also known as acyl-CoA synthetase long-chain family member 1, codes for a protein that assists in the biosynthesis of lipids and degradation of fatty acids, and SNPs within this gene have been associated with chronic diseases like diabetes.⁵⁴ Finally, ELF1 (Gene ID: 1997), also known as E74 like E26 transformation-specific related transcription factor 1, is an important positive regulator of the Hox cofactor Myeloid ectropic viral integration site 1 (MEIS1) which is involved in developmental processes.⁵⁵

Of all these DMPs and their annotated genes, RUNX1 is the only gene in which differential methylation has previously been associated with OA in humans. Specifically, RUNX1 was found to be differentially methylated in OA and healthy cartilage in humans, with OA cartilage displaying hypomethylation as compared to healthy cartilage.⁵ Additionally, a drug targeting this gene has been proposed as a therapy for the treatment of OA.⁵⁶ As of yet, none of the remaining DMPs and their annotated genes have been identified as candidate methylation loci in human OA studies.^4-10 Nevertheless, some DMPs are located in close proximity to genes that have been previously associated with OA, including CRTC1,³⁸ CRLF1,³⁹ COMP,⁴⁰ SMAD9,^41,42 and ALOX12^43,44 ( Table 1 ). However, the mechanisms by which these DMPs might influence the expression of these distally related genes is unclear.

Overall, these findings suggest that one DNA methylation pattern associated with OA is evolutionarily conserved between humans and baboons. However, further work is necessary to fully elucidate whether the mechanisms contributing to RUNX1 hypomethylation and its downstream effects in primate OA are also conserved, as well. Simultaneously, several DNA methylation patterns associated with OA do not appear to be evolutionarily conserved between humans and baboons. Differences between these 2 species may be due to general speciation events that took place during the evolution of these taxonomic groups, to slight differences in the development or manifestation of OA in these species, or artifacts of the experimental design itself. For instance, the sample size of this study (n = 10) is rather small, and all individuals included were female. The small number of individuals likely reduced our power to detect potentially important OA-related DMPs, and the inclusion of only one sex may have biased our results such that identified OA DMPs are actually female specific DMPs. Thus, to improve the identification of candidate epigenetic alterations that underlie variation in knee OA, a larger sample set that includes both sexes should be considered. Nevertheless, using baboons as a model of OA in this study has begun to clarify the evolutionary conservation of this disorder, and future research in this animal model will help provide insight into the development and progression of OA in order to begin designing preventative and therapeutic agents.¹⁶

In conclusion, we determined that the 450K array can be used to measure genome-wide DNA methylation in baboon tissues and identify significant associations with complex traits. This is the first study to specifically assess DNA methylation in skeletal tissues from a nonhuman primate using this method. Some methylation variation is related to genes that impact skeletal development and maintenance, and this may have direct downstream regulatory and phenotypic effects. Additionally, while some DNA methylation patterns associated with OA in baboons appear to be evolutionarily conserved with humans, others do not. These findings warrant further investigation in a larger and more phylogenetically diverse sample set. Last, the work presented here begins to advance areas of research that incorporate an animal model of disease and an evolutionary perspective of diseases across phylogenies.

Footnotes

Acknowledgments and Funding

We thank members of the Department of Genetics at the Texas Biomedical Research Institute for helpful discussions, including Erin Sybouts, Sophia Johnson, Mel Carless, Kara Peterson, Laura Cox, Kenneth Lange, Jerry Glenn, and Clint Christensen. The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: National Institutes of Health (NIH) Grant P01 HL028972 to AGC; Max and Minnie Tomerlin Voelcker Foundation to LMH; William and Ella Owens Foundation for Biomedical Research to LMH; SNPRC Internship Funds to GH; ASU Chapter of Sigma Xi Grant-in-Aid of Research to GH.

Authors’ Note

Genevieve Housman is currently affiliated with the University of Chicago, and Ellen E. Quillen is currently affiliated with Wake Forest University School of Medicine.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Ethical Approval

Ethical approval was not sought for the present study because the nonhuman primate tissue samples included were opportunistically collected at routine necropsy of these animals. No animals were sacrificed for this study, and no living animals were used in this study.

Animal Welfare

Guidelines for humane animal treatment did not apply to the present study because the nonhuman primate tissue samples included were opportunistically collected at routine necropsy of these animals. No animals were sacrificed for this study, and no living animals were used in this study.

ORCID iD

Genevieve Housman

Supplemental Material

The supplementary material for this article is available online.

References

Cross

Smith

Hoy

Nolte

Ackerman

Fransen

et al . The global burden of hip and knee osteoarthritis: estimates from the Global Burden of Disease 2010 study. Ann Rheum Dis. 2014;73(7):1323-30. doi:10.1136/annrheumdis-2013-204763.

Blagojevic

Jinks

Jeffery

Jordan

KP.

Risk factors for onset of osteoarthritis of the knee in older adults: a systematic review and meta-analysis. Osteoarthritis Cartilage. 2010;18(1):24-33. doi:10.1016/j.joca.2009.08.010.

Macrini

Coan

Levine

Lerma

Saks

Araujo

et al . Reproductive status and sex show strong effects on knee OA in a baboon model. Osteoarthritis Cartilage. 2013;21(6):839-48. doi:10.1016/j.joca.2013.03.003.

Delgado-Calle

Fernández

Sainz

Zarrabeitia

Sañudo

García-Renedo

et al . Genome-wide profiling of bone reveals differentially methylated regions in osteoporosis and osteoarthritis. Arthritis Rheum. 2013;65(1):197-205. doi:10.1002/art.37753.

Fernández-Tajes

Soto-Hermida

Vázquez-Mosquera

Cortés-Pereira

Mosquera

Fernández-Moreno

et al . Genome-wide DNA methylation analysis of articular chondrocytes reveals a cluster of osteoarthritic patients. Ann Rheum Dis. 2014;73(4):668-77. doi:10.1136/annrheumdis-2012-202783.

Jeffries

Donica

Baker

Stevenson

Annan

Humphrey

et al . Genome-wide DNA methylation study identifies significant epigenomic changes in osteoarthritic subchondral bone and similarity to overlying cartilage. Arthritis Rheumatol. 2016;68(6):1403-14. doi:10.1002/art.39555.

Moazedi-Fuerst

Hofner

Gruber

Weinhaeusel

Stradner

Angerer

et al . Epigenetic differences in human cartilage between mild and severe OA. J Orthop Res. 2014;32(12):1636-45. doi:10.1002/jor.22722.

Rushton

Reynard

Barter

Refaie

Rankin

Young

et al . Characterization of the cartilage DNA methylome in knee and hip osteoarthritis: methylation profile of OA cartilage. Arthritis Rheumatol. 2014;66(9):2450-60. doi:10.1002/art.38713.

Alvarez-Garcia

Fisch

Wineinger

Akagi

Saito

Sasho

et al . Increased DNA methylation and reduced expression of transcription factors in human osteoarthritis cartilage: differential DNA methylation in OA cartilage. Arthritis Rheumatol. 2016;68(8):1876-86. doi:10.1002/art.39643.

10.

Aref-Eshghi

Zhang

Liu

Harper

Martin

Furey

et al . Genome-wide DNA methylation study of hip and knee cartilage reveals embryonic organ and skeletal system morphogenesis as major pathways involved in osteoarthritis. BMC Musculoskelet Disord. 2015;16:287. doi:10.1186/s12891-015-0745-5.

11.

Kuyinu

Narayanan

Nair

Laurencin

CT.

Animal models of osteoarthritis: classification, update, and measurement of outcomes. J Orthop Surg. 2016;11:19. doi:10.1186/s13018-016-0346-5.

12.

Cucchiarini

de Girolamo

Filardo

Oliveira

Orth

Pape

et al . Basic science of osteoarthritis. J Exp Orthop. 2016;3(1):22. doi:10.1186/s40634-016-0060-6.

13.

Ameye

Young

MF.

Animal models of osteoarthritis: lessons learned while seekiseeking the “Holy Grail.” Curr Opin Rheumatol. 2006;18(5):537-47.

14.

Hernando-Herraez

Prado-Martinez

Garg

Fernandez-Callejo

Heyn

Hvilsom

et al . Dynamics of DNA methylation in recent human and great ape evolution. PLOS Genet. 2013;9(9):e1003763. doi:10.1371/journal.pgen.1003763.

15.

Ong

Tan

MacIsaac

Mah

Buschdorf

Cheong

et al . Infinium monkeys: infinium 450K array for the Cynomolgus macaque (Macaca fascicularis). G3 (Bethesda). 2014;4(7):1227-34. doi:10.1534/g3.114.010967.

16.

Cox

Comuzzie

Havill

Karere

Spradling

Mahaney

et al . Baboons as a model to study genetics and epigenetics of human disease. ILAR J. 2013;54(2):106-21. doi:10.1093/ilar/ilt038.

17.

Enard

Fassbender

Model

Adorján

Pääbo

Olek

Differences in DNA methylation patterns between humans and chimpanzees. Curr Biol. 2004;14(4):R148-9.

18.

Martin

DIK

Singer

Dhahbi

Mao

Zhang

Schroth

et al . Phyloepigenomic comparison of great apes reveals a correlation between somatic and germline methylation states. Genome Res. 2011;21(12):2049-57. doi:10.1101/gr.122721.111.

19.

Molaro

Hodges

Fang

Song

McCombie

Hannon

et al . Sperm methylation profiles reveal features of epigenetic inheritance and evolution in primates. Cell. 2011;146(6):1029-41. doi:10.1016/j.cell.2011.08.016.

20.

Pai

Bell

Marioni

Pritchard

Gilad

A genome-wide study of DNA methylation patterns and gene expression levels in multiple human and chimpanzee tissues. PLOS Genet. 2011;7(2):e1001316. doi:10.1371/journal.pgen.1001316.

21.

Barnett

Larson

A phenol-chloroform protocol for extracting DNA from ancient samples. Methods Mol Biol. 2012;840:13-9. doi:10.1007/978-1-61779-516-9_2.

22.

Triche

Weisenberger

Van Den Berg

Laird

Siegmund

KD.

Low-level processing of illumina infinium DNA methylation beadarrays. Nucleic Acids Res. 2013;41(7):e90. doi:10.1093/nar/gkt090.

23.

Fortin

Labbe

Lemire

Zanke

Hudson

Fertig

et al . Functional normalization of 450k methylation array data improves replication in large cancer studies. Genome Biol. 2014;15(12):503. doi:10.1186/s13059-014-0503-2.

24.

Aryee

Jaffe

Corrada-Bravo

Ladd-Acosta

Feinberg

Hansen

et al . Minfi: a flexible and comprehensive Bioconductor package for the analysis of Infinium DNA methylation microarrays. Bioinformatics. 2014;30(10):1363-69. doi:10.1093/bioinformatics/btu049.

25.

Huber

Carey

Gentleman

Anders

Carlson

Carvalho

et al . Orchestrating high-throughput genomic analysis with Bioconductor. Nat Methods. 2015;12(2):115-21. doi:10.1038/nmeth.3252.

26.

Zhang

Huang

Jafari

Kibbe

Hou

et al . Comparison of beta-value and M-value methods for quantifying methylation levels by microarray analysis. BMC Bioinformatics. 2010;11:587. doi:10.1186/1471-2105-11-587.

27.

Altschul

Madden

Schäffer

Zhang

Miller

et al . Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 1997;25(17):3389-402.

28.

Chen

Lemire

Choufani

Butcher

Grafodatskaya

Zanke

et al . Discovery of cross-reactive probes and polymorphic CpGs in the Illumina Infinium Humanmethylation450 Microarray. Epigenetics. 2013;8(2):203-9. doi:10.4161/epi.23470.

29.

Maksimovic

Phipson

Oshlack

A cross-package Bioconductor workflow for analysing methylation array data. F1000Research. 2016;5:1281. doi:10.12688/f1000research.8839.1.

30.

Jaffe

Irizarry

RA.

Accounting for cellular heterogeneity is critical in epigenome-wide association studies. Genome Biol. 2014;15(2):R31. doi:10.1186/gb-2014-15-2-r31.

31.

Leek

Johnson

Parker

Jaffe

Storey

JD.

The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics. 2012;28(6):882-3. doi:10.1093/bioinformatics/bts034.

32.

Morris

Beck

Analysis pipelines and packages for Infinium Human Methylation450 BeadChip (450k) data. Methods. 2015;72:3-8. doi:10.1016/j.ymeth.2014.08.011.

33.

Houseman

Molitor

Marsit

CJ.

Reference-free cell mixture adjustments in analysis of DNA methylation data. Bioinformatics. 2014;30(10):1431-9. doi:10.1093/bioinformatics/btu029.

34.

Kaushal

Zhang

Karmaus

WJJ

Wang

JSL

. Which methods to choose to correct cell types in genome-scale blood-derived DNA methylation data? BMC Bioinformatics. 2015;16(Suppl 15):P7. doi:10.1186/1471-2105-16-S15-P7.

35.

Ritchie

Phipson

Law

Shi

et al . limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47. doi:10.1093/nar/gkv007.

36.

Smyth

GK.

Linear models and empirical Bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol. 2004;3:Article3. doi:10.2202/1544-6115.1027.

37.

Benjamini

Hochberg

Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B Methodol. 1995;57(1):289-300.

38.

Ogawa

Kozhemyakina

Hung

Grodzinsky

Lassar

AB.

Mechanical motion promotes expression of Prg4 in articular cartilage via multiple CREB-dependent, fluid flow shear stress-induced signaling pathways. Genes Dev. 2014;28(2):127-39. doi:10.1101/gad.231969.113.

39.

Tsuritani

Takeda

Sakagami

Ishii

Eriksson

Hara

et al . Cytokine receptor-like factor 1 is highly expressed in damaged human knee osteoarthritic cartilage and involved in osteoarthritis downstream of TGF-β. Calcif Tissue Int. 2010;86(1):47-57. doi:10.1007/s00223-009-9311-1.

40.

Ramos

YFM

Metrustry

Arden

Bay-Jensen

Beekman

de Craen

AJM

et al . Meta-analysis identifies loci affecting levels of the potential osteoarthritis biomarkers sCOMP and uCTX-II with genome wide significance. J Med Genet. 2014;51(9):596-604. doi:10.1136/jmedgenet-2014-102478.

41.

Salazar

Gamer

Rosen

BMP signalling in skeletal development, disease and repair. Nat Rev Endocrinol. 2016;12(4):203-21. doi:10.1038/nrendo.2016.12.

42.

Wang

Green

Wang

Deng

Qiao

Peabody

et al . Bone morphogenetic protein (BMP) signaling in development and human diseases. Genes Dis. 2014;1(1):87-105. doi:10.1016/j.gendis.2014.07.005.

43.

Xiao

Zhang

Yue

et al . ALOX12 polymorphisms are associated with fat mass but not peak bone mineral density in Chinese nuclear families. Int J Obes (London). 2011;35(3):378-86. doi:10.1038/ijo.2010.157.

44.

Driban

Lapane

McAlindon

Eaton

CB.

Dietary fat intake and radiographic progression of knee osteoarthritis: data from the osteoarthritis initiative. Arthritis Care Res (Hoboken). 2017;69(3):368-75. doi:10.1002/acr.22952.

45.

Michels

Binder

Dedeurwaerder

Epstein

Greally

Gut

et al . Recommendations for the design and analysis of epigenome-wide association studies. Nat Methods. 2013;10(10):949-55. doi:10.1038/nmeth.2632.

46.

Stein

Lian

van Wijnen

Stein

Montecino

Javed

et al . Runx2 control of organization, assembly and activity of the regulatory machinery for skeletal gene expression. Oncogene. 2004;23(24):4315-29. doi:10.1038/sj.onc.1207676.

47.

Wei

Luo

Sun

Song

et al . microRNA-497 induces cell apoptosis by negatively regulating Bcl-2 protein expression at the posttranscriptional level in human breast cancer. Int J Clin Exp Pathol. 2015;8(7):7729-39.

48.

Zhao

Liu

Chen

Jiang

et al . Analysis of MiR-195 and MiR-497 expression, regulation and role in breast cancer. Clin Cancer Res. 2011;17(7):1722-30. doi:10.1158/1078-0432.CCR-10-1800.

49.

Liu

Chen

microRNA-195 promotes apoptosis and suppresses tumorigenicity of human colorectal cancer cells. Biochem Biophys Res Commun. 2010;400(2):236-40. doi:10.1016/j.bbrc.2010.08.046.

50.

Almeida

Silva

Vasconcelos

Almeida

Caires

Pinto

et al . miR-195 in human primary mesenchymal stromal/stem cells regulates proliferation, osteogenesis and paracrine effect on angiogenesis. Oncotarget. 2016;7(1):7-22. doi:10.18632/oncotarget.6589.

51.

Garvie

Boss

JM.

Assembly of the RFX complex on the MHCII promoter: role of RFXAP and RFXB in relieving autoinhibition of RFX5. Biochim Biophys Acta. 2008;1779(12):797-804. doi:10.1016/j.bbagrm.2008.07.012.

52.

Dastani

Hivert

Timpson

Perry

JRB

Yuan

Scott

et al . Novel loci for adiponectin levels and their influence on type 2 diabetes and metabolic traits: a multi-ethnic meta-analysis of 45 891 individuals. PLoS Genet. 2012;8(3):e1002607. doi:10.1371/journal.pgen.1002607.

53.

Dhanoa

Cogliati

Satish

Bruford

Friedman

JS.

Update on the Kelch-like (KLHL) gene family. Hum Genomics. 2013;7:13. doi:10.1186/1479-7364-7-13.

54.

Manichaikul

Wang

Zhao

Wojczynski

Siebenthall

Stamatoyannopoulos

et al . Genetic association of long-chain acyl-CoA synthetase 1 variants with fasting glucose, diabetes, and subclinical atherosclerosis. J Lipid Res. 2016;57(3):433-42. doi:10.1194/jlr.M064592.

55.

Xiang

Argiropoulos

Lai

Rouhi

Imren

et al . Identification of E74-like factor 1 (ELF1) as a transcriptional regulator of the Hox cofactor MEIS1. Exp Hematol. 2010;38(9):798-8,808.e1-e2. doi:10.1016/j.exphem.2010.06.006.

56.

Yano

Hojo

Ohba

Fukai

Hosaka

Ikeda

et al . A novel disease-modifying osteoarthritis drug candidate targeting Runx1. Ann Rheum Dis. 2013;72(5):748-53. doi:10.1136/annrheumdis-2012-201745.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.21 MB