Differential DNA methylation patterns in human Schlemm's canal endothelial cells with glaucoma.

Purpose: Schlemm's canal (SC) endothelial cells derived from donors with or without glaucoma showed different mechanical properties and gene expression. As an important contributor to the regulation of intraocular pressure (IOP) and pathogenesis of primary open-angle glaucoma (POAG), the heritable key epigenetic changes, methylation may play an important role in the physiologic function of SC cells. This study aims to identify differentially methylated CpG sites (DMSs) in primary cultures of human SC cells with or without glaucoma.
Methods: We examined the methylation pattern of seven strains of primary human cells (two glaucoma and five normal SC cell samples), which were isolated and characterized using established protocols. DNA methylation was profiled using Illumina Human Methylation 450 BeadChip. Raw data were extracted and exported using Illumina GenomeStudio software. After quantile normalization, DNA methylation data were analyzed using R package RnBeads in Bioconductor. DMSs were filtered with p ≤ 1E-5, methylation change ≥ 0.1, and false discovery rate ≤ 0.05. The closest genes and the location of each CpG site were annotated using R package FDb.InfiniumMethylation.hg19. Gene Ontology and pathway analysis was performed using WebGestalt. Selected DMSs were validated using the Zymo qMethyl kit.
Results: We used five non-glaucoma and two glaucomatous SC cell samples to profile genome-wide DNA methylation using Illumina Infinium Methylation BeadChips. Principle component analysis showed the separation between the glaucoma and control samples. After quality control and differential analysis, we identified 298 highly significant DMSs (p ≤ 1E-5). Among them, 221 DMSs were within 1 kb of a nearby gene. Gene Ontology analysis demonstrated significant enrichment in positive regulation of cell migration, negative regulation of endothelial cell proliferation, and stress fiber and actin filament bundles. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis showed enrichment in cell adhesion and gap junctions. Several glaucoma-related genes were identified, including TGFBR3, THBS1, PITX2, DAXX, TBX3, TNXB, ANGPT1, and PLEKHA7. We also examined differentially methylated regions (DMRs) near these CpG sites and identified significant DMRs in TBX3, TNXB1, DAXX, and PITX2. Conclusions: This study represents the first genome-wide DNA methylation profiling in cultured human primary SC cells. The DMSs were enriched in the pathways related to outflow resistance. Several DMRs were validated in glaucoma-associated genes, further suggesting the role of DNA methylation in glaucoma development. This study could provide comprehensive understanding of DNA methylation in glaucoma and its effect on aqueous humor outflow.

Introduction

Glaucoma is a progressive neurodegenerative disorder and the leading cause of irreversible blindness worldwide [1]. Primary open-angle glaucoma (POAG) is the most common subtype with an unspecified cause. POAG is characterized by the progressive loss of retinal ganglion cells, optic nerve cupping, and visual field loss [2]. Risk factors include advanced aging, elevated intraocular pressure (IOP), positive family history, and African ancestry [3]. Among them, elevated IOP is the main risk factor and the only effective target for clinical treatment [4]. The balance between secretion and drainage of aqueous humor in the anterior chamber determines IOP [5]. In most cases, elevated IOP is caused by increased outflow resistance, which is mainly generated by cells in the juxtacanalicular (JCT) region, where trabecular meshwork (TM) and Schlemm’s canal (SC) endothelial cells interact [6,7]. In glaucomatous eyes, the mechanical properties of JCT tissue are altered, becoming more fibrotic [8]. Significantly, cultured glaucomatous SC cells show greater stiffness and similar mechanical properties to glaucomatous JCT tissue [9]. Moreover, expression profiling of cultured SC cells with or without glaucoma revealed significant differential expression of genes involved in extracellular matrix (ECM) remodeling [10]. These corresponding changes suggested a heritable phenotype with cell division, which may be ascribed to epigenetic changes in cellular genome, such as histone modification and DNA methylation.

The contribution of genetic defects to POAG has been recognized for decades [11,12]. Several causal genes, such as myocilin (Gene ID 4653, OMIM 601652), optineurin (Gene ID 10133, OMIM 602432), and TBK1 (Gene ID 29110, OMIM 604834), have been identified [13-15]; however, mutations in these Mendelian genes account for ≤ 10% POAG cases worldwide [16,17]. Additionally, large-cohort genome-wide association studies (GWASs) have identified many associated single nucleotide polymorphisms (SNPs) associated with POAG and IOP [18-22]. Currently, POAG is considered to be a complex disease affected by genetic and environmental factors [17,23]. DNA methylation, the most well-studied epigenetic regulators in response to environmental factors, may inhibit gene transcription by covalent addition of methyl groups to DNA using DNA methyltransferase (DNMT), and may reverse the inhibition by cell division or removal of methyl groups through tet methylcytosine dioxygenase (TET) enzymes [24]. The methyl group is susceptibly added to CpG sites, where cytosine is followed by a guanine. The effect of DNA methylation is dependent on the density of CpG sites, especially high CG content regions named CpG islands, and their locations within a gene [25]. Although the mechanism of environmental factors in DNA methylation remains unclear, especially the gene-specific DNA methylation patterns, numerous studies support the susceptibility of DNA methylation to environmental risk factors in relation to eye diseases [26,27]. However, little is known about the effect of DNA methylation in the outflow facility glaucoma pathogenesis. In this study, we aimed to identify differentially methylated DNA loci in glaucomatous SC endothelial cells using genome-wide methylation profiling. The identification of such DNA methylation patterns holds promise to provide more information about how epigenetic factors play a role in glaucoma pathogenesis.

Methods

Isolation and culture of SC cells

Human SC endothelial cells were isolated from postmortem eyes provided by Miracles in Sight (Winston-Salem, NC), Midwest Eye Bank (Ann Arbor, MI), National Disease Research Interchange (Philadelphia, PA), or Life Legacy (Tucson, AZ) within 36 h of death with enucleation occurring ≤ 6 h after death. Isolation and culture of primary SC cells were processed according to the established techniques [28]. All strains of primary SC cells were characterized using five inclusion criteria: (1) the expression of vascular endothelial cadherin, (2) a net trans-endothelial electrical resistance ≥10 ohms·cm2, (3) lack of myocilin induction by dexamethasone, (4) monolayers that formed non-overlapping, linear-arranged morphology, and (5) contact-inhibited growth. SC cells were isolated from two donors with glaucoma and five donors without glaucoma. The determination of glaucoma was based on a combination of information provided by the eye banks: a history of ocular hypertension or glaucoma or the presence of glaucoma eye drops on the patient’s medication list. All the cultured SC cells were harvested between passages 3 and 5, as cell stiffness remains consistent through passage 6 [9].

DNA extraction

Genomic DNA was extracted using the AllPrep DNA Mini spin column kit (Qiagen, Valencia, CA) according to the manufacturer’s instructions. Briefly, cultured cells were harvested and lysed using provided buffer RLT. The lysate was collected and loaded to the spin column. After purification, 50 µl of RNase-free water was added to elute DNA from the spin column. The DNA concentration was measured using the TECAN Infinite M200PRO microplate reader (Tecan, Mannedorf, Switzerland).

Identification of differentially methylated CpG sites

DNA samples were labeled and hybridized to the Illumina Infinium HumanMethylation450 BeadChips and scanned by the Illumina iScan System following the manufacturer’s protocol. Briefly, 1 µg of high-quality DNA samples were treated with bisulfite to convert unmethylated cytosines to uracils while methylated cytosines were protected and remained cytosine. Predesigned probes for each CpG site, flagged with special fluorescence, were used to determine whether the base at a given locus was converted or not converted, leading to a profile of the DNA methylation status at this locus. Methylation status was calculated as the β value: the ratio of the signal from a methylated probe relative to the sum of the methylated and unmethylated probes. This value ranges from 0 (unmethylated) to 1 (fully methylated). This platform enabled sensitive and reproducible methylation profiling with quantitative measurements. Quality controls from the Illumina BeadChips were checked for experimental procedures, including bisulfite conversion, staining, hybridization, target removal, specificity, negative, and non-polymorphism. The internal bisulfite conversion control Infinium I and II probes indicated the successful conversion upon initial bisulfite treatment of the DNA samples. The experiment was performed by the Genotyping Core Facility at University of Miami (Coral Gables, FL).

Probe and intensity data for each sample were imported into R using package RnBeads in Bioconductor [29], followed by preprocessing to filter probes outside the CpG context, SNP probes, sex chromosome probes, probes without intensity value, and probes with a low standard deviation. Then preprocessed data were normalized using the subset quantile normalization (SWAN) function. Data from duplicated samples were verified to have similar β results and averaged for further analysis. The Rnbeads.hg19 package was used to annotate the location of each probe in chromosomes and CpG islands. Principle component analysis (PCA) was performed to show genome-wide methylation variation between each biologic sample. Next, the glaucoma samples were grouped and compared to the set of control samples to identify differentially methylated CpG sites (DMSs) using the limma package. The most significant DMSs were filtered as p ≤ 1E-5, false discovery rate (FDR) ≤ 0.05, and absolute difference in β values (|Δβ|) greater than 0.1. Because there were only two glaucoma SC samples, we assumed that the control and glaucoma groups had similar variations. The FDb.InfiniumMethylation.hg19 package [30] and Infinium Human-Methylation450K v1.2 product files were imported to annotate all CpG sites with nearest genes, distance to the gene, and location within the gene.

Gene Ontology and pathway analysis

First, the list of annotated genes within 1 kb of significant CpG sites (p ≤ 1E-5, |Δβ| ≥ 0.1) was uploaded into WebGestalt (WEB-based Gene SeT AnaLysis Toolkit) for Gene Ontology analysis, categorized by biologic process, molecular function, and cellular component. WebGestalt exhibited a network view in each category to reveal the relationship of each class and the location of proteins produced by target genes. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was also performed with WebGestalt to explore the underlying biology that may contribute in the development of glaucoma.

Differentially methylated regions

The OneStep qMethyl Kit (Zymoresearch, Irvine, CA) was used to examine differentially methylated CpG regions containing selected candidate CpG sites. To perform quantitative real time polymerase chain reaction (RT-PCR) with methylation sensitive restriction enzymes (MSREs), we designed primers that included two or more MSRE targeting sites in the amplicon flanking targeted CpG sites (Table 1). The methylated sites in the region were protected from digestion of restriction enzymes and failed RT–PCR amplification. Then the percentage of methylation was estimated using the Ct between the test group (with restriction enzymes) and the reference group (without restriction enzymes). Statistical analysis was performed with GraphPad 8 using an unpaired t test with Welch’s correction (p ≤ 0.05).

Table 1

Selected primers for DMR validation using Zymo OneStep qMethyl Kit.

Probe ID	Gene Symbol	Gene Description	Forward Primer	Reverse Primer	Amplicon Size, bp
cg13078798	TGFBR3	Transforming growth factor, β receptor III	TGA GGA TGC ATT TGG ATG AG	AAG GCC TGT AGG CTG CTG TA	200
cg16055869	CXCL5	Chemokine (C-X-C motif) ligand 5	CTC TAC CGA GAG GGA CGT TG	AGG AGG AGC ATC TCC CAG AG	171
cg21951797	TNXB	Tenascin XB	GGA ACA CAA ATG ACC CGT AGA	AGC CCA CGT TTT TCT AGT GA	345
cg07139946			GGC TCA GTC AGA CCA GGA GA	AGG GCC AGT TTG ACT CCT TT	255
cg17982478	PLEKHA7	Pleckstrin homology domain contaning, family A member 7	CCC GAA CGC TCA GTA AAC C	CCA GCT GAA AAC TTG GGA AA	273
cg03133735	PITX2	Paired-like homeodomain 2	TAC TAT GCG TTG CCG ATT CC	AAA GAC CCC TGC TCC AAA AT	268
cg01733176			GTG TGG AAG GAG CTG GAC AT	ATT GCT TGC TTT GCT TGG AC	295
cg26500914 cg24498636	DAXX	Death-domain associated protein	GTA CCC CAT CCA CAC CTC AC	GGG CTG AGT GCT CTG ACT TT	263
cg05431670			GTC ACA GAG TTT CCG CCT TC	GTA GAA GCA CCG GGT GAA AA	268
cg16277169	TBX3	T-Box 3	GAA ACT GGA CGA AAG GTG GA	GGG CCA ACA GTT CTT CAA C	300
cg18161956			TCA GCA GCG AAA AGG TGA G	CAG AGA GGC TAA GGG GCT TT	300
cg09413529			TGC CCG TTG AAG AAC TGT TG	CAC CAT CTC GTC CAG CAC T	227
cg11246938			AGT GCT GGA CGA GAT GGT G	GAG AGC AAA GAG GAG CAT GG	201
cg09053536			GAA TTC AGT TTC GGG GAA CA	TTT GGC AAC TGA GGA GCA AT	322
cg13661397	NET1	Neuroepithelial cell transforming 1	GAT GCG CTC AGG AGT TAA GG	CAG CGA TCA GCC AAT CAG T	130

Results

Clinical phenotypes

This study examined human primary SC endothelial cells from seven donors, two with glaucoma and five without glaucoma (Table 2). The average age was approximately 56 years for the controls, while the average age was 78 years for the patients with glaucoma. Although age was included as a cofactor using the R package limma in Bioconductor [31-35] to reduce the influence of varying age in the analysis, PCA showed variation of methylation pattern in samples with aging, as well as a clear separation between the control group and the glaucoma group (Figure 1).

Table 2

Clinical characteristics of donors to derive the primary endothelial cells.

Group	Cell Strain No.	Gender	Age (years)
Glaucoma	SC57 g	Male	78
	SC64 g	Male	78
Control	SC71	Male	44
	SC73	Female	37
	SC76	Female	59
	SC78	Male	77
	SC80	NA	62

Figure 1

Principal component analysis of SC cells using genome-wide methylation data. F: female. M: male. NA: no gender information. Blue: non-glaucoma control SC samples. Red: glaucomatous SC samples.

Differential DNA methylation analysis

The Illumina internal control probes for bisulfite conversion indicated the successful bisulfite conversion of all DNA samples. After quality control, an average of 485,577 genomic loci were detected with methylation status in each sample. There were 383,976 identified methylation sites in the control and glaucomatous SC cells. The distribution of identified methylation status across these sites was binomial, meaning with two peaks nearing 0.1 or 0.9 of the β value, instead of one peak around 0.5 (Figure 2A).

Figure 2

Differential methylation in human primary cultured SC cells. A: Binomial distribution of the Schlemm’s canal (SC) β values. Two peaks were found near 0.1 and 0.9. B: A Manhattan plot of differential methylation changes in single sites. The red line represents the threshold of a p value 5E-8 while the blue line represents the threshold of a p value 1E-5. C: The heat map shows different methylation patterns in SC cells with or without glaucoma.

The glaucoma samples were grouped (n = 2) and compared to the set of control samples (n = 5) to identify DMSs using the limma package. For greater accuracy of differential methylation analysis, the significance level (-logP value) of each DMS was plotted across all chromosomes in the Manhattan plots (Figure 2B). The highly statistically significant DMSs (p ≤ 1E-5, FDR < 0.05) were then filtered with significant changes ((|Δβ| ≥ 0.1) in the glaucomatous SC cells to identify 298 significant DMSs, including 106 increased DMSs and 192 decreased DMSs (Appendix 1). These sites were then used to generate a heat map to visualize the difference in methylation patterns between the two groups (Figure 2C).

To identify glaucoma-related genes related with epigenetic alterations, we annotated all significant DMSs with their relative distance to the nearest genes. The majority of DMSs are located in the body of genes; 15% are in the untranslated region (UTR), and 13% within 1,500 bases of the transcription start site (TSS; Figure 3). According to previous studies, gene expression has a greater chance of being inhibited through methylation presented at or near the TSS, including the promoter region [25]. To identify methylation that potentially contributes to the altered expression of related genes, we filtered the distance to the nearest genes at 1,000 base pairs for greater chance of regulation. In this circumstance, we identified 221 significant DMSs located near or within 208 genes. Many of these genes are related to glaucoma, such as TGFBR3 (Gene ID 7049, OMIM 600742), LAMA3 (Gene ID 3909, OMIM 600805), PITX2 (Gene ID 5308, OMIM 601542), TNXB (Gene ID 7148, OMIM 600985), ANGPT1 (Gene ID 284, OMIM 601667), and PLEKHA7 (Gene ID 144100, OMIM 612686)[36-41]. Gene Ontology analysis indicated that these 208 genes were significantly enriched in functions related with positive regulation of cell migration, negative regulation of endothelial cell proliferation, and stress fiber and actin filament bundles (Appendix 2). KEGG pathway analysis showed significant enrichment in cell adhesion and gap junctions (Appendix 3). These results suggest that glaucoma SC cells may provide greater resistance to aqueous humor outflow with inheritable gene expression alteration.

Figure 3

The percentage of DMS location related to the nearest genes. UTR: untranslated region. TSS: transcript start sites. Intergenic: no annotation information.

Differentially methylated region in glaucoma-related genes

The differential methylation of one CpG site has less chance to turn over gene expression, yet differential methylation of a CpG region (DMR) has a larger possibility [25]. Based on the distance of the genes, the relationship with the CpG islands, and the function of the nearest gene, we selected ten glaucoma-related genes to examine potential DMRs, including previously reported glaucoma associated genes CXCL5 (Gene ID 6374, OMIM 600324), TNXB (Gene ID 7148, OMIM 600985), ANGPT1 (Gene ID 284, OMIM 601667), PLEKHA7 (Gene ID 144100, OMIM 612686), THBS1 (Gene ID 7057, OMIM 188060), PITX2 (Gene ID 5308, OMIM 601542), DAXX (Gene ID 1616, OMIM 603186), and TBX3 (Gene ID 6926, OMIM 601621), and two genes that were differentially expressed in glaucoma SC cells, TGFBR3 (Gene ID 7049, OMIM 600742) and NET1 (Gene ID 10276, OMIM 606450)[10]. Distinct from bilsulfite sequencing, the Zymo OneStep qMethyl kit can measure the methylation level of a region with a small amount of DNA sample (Zymoresearch Cat No.D5310). Using MSRE-based RT–PCR, we examined DMRs containing one or multiple significant DMSs in these candidate genes (Table 3). Corresponding to the microarray data, significant DMRs showed the same direction of methylation status alteration and mainly located at the gene with multiple significant DMSs, such as TBX3, which has five significant DMSs and four significant DMRs flanking these CpG sites, suggesting a greater chance of gene expression inhibition. This result indicated different methylation levels between the DMSs and the DMRs.

Table 3

Selected genes with differential methylation and relative location of these CpG sites in the genes.

Probe ID	Gene	Description	∆β in array	p in array	∆β in validation	p in validation	Location in gene
cg13078798	TGFBR3	Transforming growth factor, β receptor III	0.54	2.17E-07	−0.08	0.376	Body
cg16055869	CXCL5	Chemokine (C-X-C motif) ligand 5	0.24	4.58E-06	0.08	0.138	TSS200
cg21951797	TNXB	Tenascin XB	0.47	7.13E-06	0.06	0.002	Body
cg07139946			−0.18	1.42E-06	0.17	0.16
cg15724328	ANGPT1	Angiopoietin 1	−0.31	5.76E-06	N/A*		Body
cg17982478	PLEKHA7	Pleckstrin homology domain contaning, family A member 7	0.36	5.35E-07	0.005	0.42	TSS1500
cg04827020	THBS1	Thrombospondin 1	0.27	1.74E-06	N/A*		Body
cg16933229	PITX2	Paired–like homeodomain 2	0.6	3.72E-07	N/A*
cg03133735			0.6	1.17E-07	0.02	0.001
cg01733176			0.38	1.34E-08	0.1	0.06
cg26500914	DAXX	Death-domain associated protein	0.5	4.11E-06	0.32	0.08	Body
cg24498636			0.48	2.26E-06	0.32	0.08
cg05431670			0.28	2.37E-09	0.07	0.001	TSS200
cg16277169	TBX3	T-Box 3	−0.32	1.10E-07	−0.56	0.04	Body
cg18161956			−0.44	4.25E-08	−0.14	0.734
cg09413529			−0.49	7.55E-06	−0.82	0.001
cg11246938			−0.39	5.36E-08	−0.19	0.051
cg09053536			−0.51	7.33E-06	−0.89	0.001
cg13661397	NET1	Neuroepithelial cell transforming 1	0.47	8.37E-10	−0.00043	0.861	TSS200

* Not available due to shortage of sites for methylation sensitive restriction enzyme. TSS200: 200 base pairs upstream of the transcription start site; TSS1500: 1500 base pairs upstream of the transcription start site.

Discussion

To our knowledge, this is the first genome-wide DNA methylation profiling of primary human SC endothelial cells derived from patients with or without glaucoma. The glaucomatous SC cells displayed a significantly different methylation pattern than control SC cells. The majority of DMSs are located at the body of genes. Significant DMSs were enriched in genes related with stress fiber and actin filament bundles, cell adhesion and gap junctions, which may contribute to the outflow resistance. Using a microarray and the Zymo OneStep qMethyl kit, we identified several DMSs and DMRs in glaucoma-related genes, showing the potential contribution of epigenetic regulation in glaucoma development.

Despite being the most studied epigenetic regulators related to environmental factors, the role of DNA methylation is understudied in human ocular diseases, especially those related with the anterior chamber. Due to the relatively large number of TM cells and easy access in the eye, it is easier to isolate and culture primary TM cells than SC cells [42]. A detection of DNMT1 (Gene ID 1786, OMIM 126375)transcripts in TM tissues and primary TM cells showed consistent DNMT expression in vitro [43]. The global DNA methylation was increased in age-matched glaucoma TM cells, accompanied with increased expression of TGFβ1 (Gene ID 7040, OMIM 190180) [44]. The treatment of normal TM cells with TGFβ1 leads to increased expression of DNMT1 (Gene ID 1786, OMIM 126375) and decreased RASAL1 (Gene ID 8437, OMIM 604118) expression, similar to those in glaucoma TM cells [44]. Reduced TGFβ1 promoter methylation and enhanced global DNA methylation were also examined in glaucomatous lamina cribrosa cells [45], further suggesting the association of alteration in DNA methylation and development of glaucoma. However, there was no significant difference between global DNA methylation level of glaucoma and normal SC cells in this study, which may be ascribed to the loss of global methylation in aging samples of both groups [46].

The deep location and limited number of cells in the SC lead to an extremely low successful rate of primary human SC cell culture using postmortem human donor tissues [28]. Moreover, the late onset nature of POAG raises greater difficulty in populating cells derived from older donors with glaucoma. The small sample size and age variation limited the statistical analysis of profiling, so we focused on the methylation of the previously reported glaucoma-related genes. Integrated with these SC expression data [10], TGFBR3 and NET1 showed potential correlation of DMSs with the opposite direction of expression change. The regions around the significant DMSs of TGFBR3 and NET1 were not significantly different. The variations between the DMS and DMR validation were probably a result of probe-based microarray that CpG sites were differentially methylated, but the level of the region was not significantly different. In the whole view, there were 40 DMSs located in the body of TGFBR3, but only four DMSs (10%) with a p value of less than 0.05. The comparison of multiple probes in a consecutive genetic region showed the potential failure in validation of DMRs. However, genes with multiple significant DMSs have several validated DMRs, including PLEKHA7, PITX2, DAXX, and TBX3 (Table 3). These results indicated the complexity of methylation regulation by altering at several locations, but not the whole gene. There was no significant expression change in correspondence with these DMRs, which may be due to different set of samples used for methylation and expression profiling or other expression-regulating mechanisms such as histone modifications. Among these identified genes with DMRs, PLEKHA7 and PITX2 were associated with primary glaucoma [41,47]; while TBX3 was crucial in retina development and DAXX was a key component mediating retinal cell death [48,49]. The different methylation level on these genes may propose the existence of inheritable gene malfunction instead of mutation in patients with glaucoma.

This study was limited by several factors, especially the small sample size. Age variation may have contributed to the differences in methylation. Additionally, due to the limited amount of DNA for methylation profiling, we used Illumina BeadChips instead of DNA methylation sequencing, so we may have missed some DMSs. Thus, we increased the significance of DMSs, targeted on the glaucoma-related genes, and validated DMRs containing DMSs. It will be critical to derive more human primary cells, which was limited by the decreasing availability and high cost of human donor eyes. We believed that this study was able to provide comprehensive understanding of glaucoma methylation in outflow facility. By updating our knowledge of the mechanism of methylation regulation, we can integrate these human data with genetics, expression, and functions.