Multiple Differentially Methylated Regions Specific to Keratoconus Explain Known Keratoconus Linkage Loci

Abstract

Purpose: Keratoconus (KTCN) is a complex eye disorder resulting in loss of visual function. Its development is affected by genetic and environmental components. The aim of this study was to unravel the role of epigenetic factors in KTCN.

Methods: To verify if DNA methylation may play a role in KTCN development, reduced representation bisulfite sequencing of five KTCN and five non-KTCN human corneas was performed.

Results: Multiple KTCN-specific differentially methylated regions were detected and many of them overlap previously identified KTCN linkage loci (3p14.3, 5q35.2, 13q32.3, 15q24.1, and 20p13) and chromosome arms that have been linked to KTCN (2q, 4q, 5p, 9p, 14q, and 17q). Reanalysis of the previously described RNA sequencing dataset of 25 KTCN and 25 non-KTCN human corneas revealed that 12 genes downregulated in KTCN and 6 upregulated genes overlapped or were located in the near vicinity of the identified differentially methylated regions. Particularly interesting were the DNA methylation changes in WNT3 and WNT5A encoding Wnt ligands, as they provide a potential explanation for the Wnt signaling pathway dysregulation observed in KTCN.

Conclusions: We presented the results of data analysis from the first study of DNA methylation changes in human KTCN corneas compared to non-KTCN samples. We were able to identify genomic regions with distinct patterns of DNA hypo- and hypermethylation and link them to previously found KTCN susceptibility loci as well as transcriptomic disruption of Wnt signaling pathway observed in KTCN.

Keratoconus (KTCN; OMIM 148300) is a degenerative eye disorder characterized by progressive thinning of the cornea, leading to its conical shape and impairment of visual function.¹ The estimated prevalence of KTCN in the general population is between 1 in 2000¹ and 1 in 375 individuals.² KTCN is a complex, heterogeneous disease with multifactorial etiology.³ Its development is affected by environmental and behavioral factors, such as eye rubbing and wearing contact lenses, but also has a strong genetic component indicated by familial inheritance, discordance between dizygotic twins, and association with other known genetic disorders.⁴ Multiple loci and gene variants have been linked to KTCN, but they tend to be specific to particular families or populations.⁵ 

Epigenetic modifications, such as DNA methylation and histone marks, have been suggested to play vital role in shaping complex phenotypes, including disease etiology.⁶ These modifications have been shown to be both highly plastic, responding to environmental conditions, and able to alter gene expression in heritable manner.⁷ As a result, they offer a potential solution to the problem of “missing heritability,” explaining why identified DNA variants account for only a fraction of observed heritability of many complex traits.⁷ 

Large-scale epigenetic studies have never been implemented in KTCN research. However, methylation status in KTCN was studied in the TIMP3 promoter. TIMP3 encodes tissue inhibitor of metalloproteinase 3 and, due to its overexpression leading to apoptosis in corneal stromal cell, is proposed as one of the candidate genes for KTCN.⁸ Methylation-specific quantitative PCR performed in three keratoconic and three healthy corneal tissues excluded promoter methylation in the TIMP3, suggesting that factors other than epigenetic factors affect TIMP3 expression in KTCN.⁹ 

Herein, we present the results of a DNA methylation study performed in human corneas derived from five KTCN and five non-KTCN individuals by using reduced representation bisulfite sequencing (RRBS).¹⁰ We were able to identify multiple regions that show distinct patterns of DNA hypo- and hypermethylation in KTCN. The majority of these differentially methylated regions (DMRs) were located within annotated gene bodies, including genes coding for two Wnt ligands that we previously discussed in the context of KTCN.¹¹ We also found that many of detected DMRs overlap previously described KTCN susceptibility regions.⁵ 

Methods

Patients

All individuals underwent a complete ophthalmic evaluation in the Department of Ophthalmology, Medical University of Warsaw, Poland. The detailed ascertainment and examination processes in KTCN patients have been previously described.^11,12 Briefly, the diagnosis of KTCN was based on both corneal topography and pachymetry assessment, a best-corrected visual acuity testing, and a slit-lamp examination. Patients were considered as having KTCN if they had at least one clinical sign, including corneal thinning, Vogt's striae, or Fleischer rings in slit-lamp biomicroscopy concurrent with the changes in topography and pachymetry assessment (Supplementary Table S1). 

The study protocol was approved by the Institutional Review Board at Poznan University of Medical Sciences. All individuals provided informed consent after the possible consequences of the study were explained, in accordance with the Declaration of Helsinki. 

Material Collection and DNA Extraction

The corneal tissues were obtained from five KTCN patients and five non-KTCN individuals undergoing a penetrating keratoplasty procedure. The non-KTCN corneas, used as controls, were collected from patients who were referred for corneal transplantation for reasons other than KTCN, including bullous keratopathy, corneal scarring, and perforations (Supplementary Table S1). The corneas were preserved in the Eusol-C solution (Alchimia, Ponte San Nicolo PD, Italy) directly after corneal button excision during the surgical procedure, and stored at 4°C until DNA extraction. Genomic DNA samples from the corneas were extracted using the Cells and Tissue DNA Isolation Kit (Norgen Biotek, Thorold, ON, Canada) according to the manufacturer's protocol. 

RRBS Library Preparation and Sequencing

RRBS libraries were performed using the NEXTflex Bisulfite Library Preparation Kit (Bioo Scientific, Austin, TX, USA) according to manufacturer's instructions. In brief, 500 ng of high-molecular-weight genomic DNA was used for MspI digestion (New England Biolabs, Ipswich, MA, USA) followed by DNA purification using Agencourt AMPure XP magnetic beads (Beckman Coulter, Beverly, MA, USA). Purified fragmented DNA samples with 5′-CG-3′ overhangs were subject into end repair, adenylation, and adaptor ligation reactions. Bisulfite conversion was performed using the EZ DNA Lightning Kit (Zymo Research, Irvine, CA, USA), according to the manufacturer's instructions, and it was followed by PCR amplification. The PCR products were purified using Agencourt AMPure XP magnetic beads (Beckman Coulter). Prior to sequencing, RRBS libraries were quantified using Qubit instrument (ThermoFisher Scientific, San Jose, CA, USA) and qualified using Agilent 2100 Bioanalyzer High Sensitivity chips (Agilent Technologies, Santa Clara, CA, USA). Paired-end sequencing (2 × 100 bp) was performed on the HiSeq 1500 platform (Illumina, Scientific, San Diego, CA, USA). 

RRBS Data Analysis

Illumina adapters and poor-quality regions (mean Phred quality, <5) were trimmed from short reads by using the BBDuk2 program from the BBTools suite (http://jgi.doe.gov/data-and-tools/bbtools/). Bismark¹³ was used to align filtered read pairs to the human reference genome (GRCh38) and calculate DNA methylation percentage at each CpG site. The first five bases from the 5′ end of all reads were excluded from methylation calling based on the methylation bias observed in the M-bias plots from Bismark reports. The methylKit package was used to create sample correlation and principal component analysis (PCA) plots (Supplementary Figs. S1 and S2).¹⁴ The DSS package¹⁵ was used to detect DMRs between KTCN and non-KTCN samples by using recommended settings for RRBS experiments. The DSS package performs individual differential methylation tests for each CpG-context cytosine and merges the adjacent significant results (<50 bp distance, minimum 3 CpGs, and minimum 10 percentage point difference in methylation) into DMRs. The raw P value threshold of 0.001 was selected, as suggested by DSS documentation for RRBS experiments. BEDTools¹⁶ was used to find the closest gene for each identified DMR. The karyotype plot showing found DMRs and KTCN susceptibility loci (Fig. 1) was created using the karyoploteR package.¹⁷ KTCN linkage loci overlapped by DMRs were visualized (Fig. 2) using the Gviz package.¹⁸ The plot showing the overlap between DMRs and different types of annotated genomic features (Supplementary Fig. S3) was prepared using the ChIPseeker package.¹⁹ 

Figure 1

 

Differential methylation regions (DMRs) between KTCN and non-KTCN samples. The karyotype plot shows hypermethylated (red) and hypomethylated (green) DMRs located on each chromosome (centromeres are colored red). KTCN linkage regions (blue) and chromosomal arms linked to KTCN (aquamarine) are shown below chromosomes and were taken from a paper by Karolak and Gajecka.⁵

Figure 2

 

KTCN linkage loci overlapped by differential methylation regions (DMRs). The plots show the detailed view of the KTCN linkage loci 3p14.3 (A), 5q35.2 (B), 13q32.3 (C), 15q24.1 (D), and 20p13 (E) overlapped by DMRs. Hypermethylated and hypomethylated DMRs are colored red and green, respectively. Only protein-coding genes are shown.

Figure 2

 

Continued.

RNA Sequencing (RNA-Seq) Data Reanalysis

To compare the results of the differential methylation analysis with changes on the transcriptome level, the previously described RNA–Seq dataset of 25 KTCN and 25 non-KTCN human corneas was used.¹¹ Adapter trimming and poor-quality region removal were performed in the identical way to the RRBS data analysis. Short reads matching human rRNA sequences were additionally removed using BBDuk2. Salmon²⁰ was used to estimate the expression values for genes and transcripts from GENCODE 25 (Ensembl 87) annotations while accounting for sequence-specific and GC content bias. Differential expression analysis was performed using a previously published protocol²¹ using the limma package.²² Genes were determined to be differentially expressed based on the 0.01 false discovery rate threshold value and 1.5-fold change cutoff. 

Results

All of the analyzed samples had a similar number of cytosines covered by at least 10 reads (2–3 million). Based on the methylation percentage of 457,621 sites that had this coverage in all samples, sample correlation plot (Supplementary Fig. S1) and sample PCA plot (Supplementary Fig. S2) were prepared. 

Overall, 112 DMRs between KTCN and non-KTCN samples were detected. Of these, 43 regions were hypermethylated in KTCN corneas, while 69 showed signs of hypomethylation (Fig. 1). Seven DMRs (Fig. 2) overlapped known KTCN linkage loci, namely, 3p14.3,²³ 5q35.2,²⁴ 13q32.3,²⁵ 15q24.1,²⁶ and 20p13,²⁷ whereas 15 DMRs were located on six chromosome arms, namely, 2q, 4q, 5p, 9p, 14q, and 17q, that have been linked to KTCN.^5,28 The detailed information regarding each region is presented in the Table. None of the identified DMR regions overlapped known genome-wide association study loci.^29–35 

Table

 

A List of Annotated Differential Methylation Regions Detected in This Study

 

The majority of DMRs (65%, 73 regions) overlapped known genes (the detailed overlap status between DMRs and different categories of genomic annotations is presented in Supplementary Fig. S3). For the remaining 39 regions, the nearest annotated gene was located within a distance between 93 bp and over 360 kb (Table). The comparison of the set of DMR-related genes with the results of differential expression analysis (Supplementary Table S2) based on the previously described corneal RNA-Seq data¹¹ showed that 18 genes were common to both datasets. Among the genes overlapping recognized DMRs, 12 genes (IQGAP2, SYNJ2, CYP1B1, MYO1G, WNT5A, PARVB, MGLL, CDC25B, PSG3, FHL2, CAMK1D, and THEMIS) showed downregulated expression in KTCN corneas, whereas 6 genes (WNT3, RB1, AC098617.1, RPS6KA2, PELI2, and PLXNA4) were upregulated in KTCN tissues (Table). 

Discussion

Despite many years of research, the genetic factors responsible for the development of KTCN remain elusive. A number of genomic regions and sequence variants have been associated with KTCN but turned out to be specific to particular families or populations and failed to explain the cause of the disease in the general population.⁵ A potential solution to this problem of missing heritability has been recently suggested, as epigenetic modifications, including DNA methylation and histone modifications, were shown to affect phenotype in heritable manner.^6,7 To verify if this is the case for KTCN, we performed the comprehensive analysis of DNA methylation in human KTCN and non-KTCN corneas using the RRBS approach.¹⁰ 

Large-scale epigenetic studies have never been implemented in KTCN research. However, methylation status in KTCN was studied in TIMP3 promoter. TIMP3 encodes tissue inhibitor of metalloproteinase 3 and due to its overexpression leading to apoptosis in corneal stromal cell is proposed as one of the candidate genes for KTCN. Methylation-specific quantitative PCR performed in three keratoconic and three healthy corneal tissues excluded promoter methylation in the TIMP3, suggesting that other than epigenetic factors affect TIMP3 expression in KTCN. 

In this study, RRBS was used to analyze DNA methylation in human corneas from five KTCN and five non-KTCN individuals. We identified 112 DMRs with distinct patterns of DNA hypermethylation (43 regions) or hypomethylation (69 regions) in KTCN and assigned the nearest gene to each locus. Interestingly, the majority (65%) of the DMRs were located within annotated gene bodies and not in the intergenic space, as it would be expected for regulatory regions (e.g., promoters) commonly associated with CpG islands studied in RRBS experiments. Many of the identified DMRs overlapped previously described KTCN susceptibility loci. Seven DMRs were located on the known KTCN linkage regions 3p14.3,²³ 5q35.2,²⁴ 13q32.3,²⁵ 15q24.1,²⁶ and 20p13,²⁷ which further verifies the previously described linkage between these chromosomal loci and KTCN. Moreover, 15 DMRs were located on the six chromosome arms 2q, 4q, 5p, 9p, 14q, and 17q, previously reported as possibly linked to KTCN.^5,28 These linkage regions contain hundreds of genes and our analysis may allow for significant narrowing of the genomic region of interest and reduce the list of putative KTCN genes from these chromosome arms. None of the identified DMRs region overlapped known genome-wide association study loci, including HGF, RAB3GAP1, LOX, RXRA-COL5A1, MPDZ-NF1B, and BANP-ZNF469.^29–35 

The next step was comparing the RRBS results with reanalyzed data from the published RNA-Seq study previously performed in human KTCN and non-KTCN corneas. The most interesting finding was that two DMRs located on chr3:55,474,327-55,474,562 and chr17:46,769,843-46,769,996 overlapped the exons of Wnt family member genes WNT5A and WNT3, respectively. These genes encode important ligands of the Wnt signaling pathway that was previously found to be downregulated in KTCN corneas.¹¹ Wnt signaling is involved in regulating the proliferation of corneal epithelial stem cells,³⁶ which suggests that the dysregulation of genes involved in this biological pathway is not incidental and may play a role in KTCN development. Moreover, the first DMR located on chr3:55,474,327-55,474,562 also overlaps the 3p14-q13 KTCN susceptibility locus detected by linkage analysis in an Italian family.²³ COL8A1 has been first suggested as candidate gene located within this locus, but no pathogenic variants of COL8A1 have been found in KTCN patients.^23,37 On the other hand, the second region chr17:46,769,843-46,769,996 is located on the 17q chromosome arm that shows linkage to KTCN in the Hispanic population.³⁰ 

Although both regions are hypermethylated in KTCN, WNT5A is transcriptionally downregulated, whereas WNT3 is transcriptionally upregulated in KTCN corneas.¹¹ This discrepancy can be explained by the fact that DNA methylation of the gene body can affect alternative splicing and indirectly gene expression by using multiple different mechanisms.³⁸ 

Other differentially expressed genes overlapping or located in the near vicinity of the DMRs include IQGAP2, SYNJ2, CYP1B1, MYO1G, PARVB, MGLL, CDC25B, PSG3, FHL2, CAMK1D, THEMIS (downregulated), RB1, AC098617.1, RPS6KA2, PELI, and PLXNA4 (upregulated). 

Apart from providing some explanation behind the transcriptomic changes observed in KTCN, DNA methylation can also be potentially linked to other important KTCN-related phenomena. Increased levels of interleukin-6 and tumor necrosis factor–α have been repeatedly found in the tears of KTCN patients.³⁹ Their levels in the tear fluid are also increased by eye rubbing, a known KTCN risk factor.³⁹ Both proteins have also been shown to be able to induce changes in DNA methylation that affect gene expression.^40,41 In this way, DNA methylation provides us with the unified framework linking hereditary and environmental factors with transcriptomic and phenotypic effects specific to KTCN. 

There was high heterogeneity of the control group used in this study, as it was not possible to obtain “healthy” corneas from living individuals for comparative analyses. The small sample size disabled a reliable statistical analysis, and due to this limitation, the general conclusions could not be made. Additionally, the lower sample correlation (Supplementary Fig. S1) and higher sample variability on PCA plot (Supplementary Fig. S2) observed in KTCN samples may suggest that the methylation dysregulation in KTCN results in biological heterogeneity exceeding that of the control group. However, the results of this pilot study have given an additional insight into the role of methylation in KTCN development and showed the possible indications for further KTCN research. 

In summary, we presented herein the results of data analysis from the first study of DNA methylation changes in human KTCN corneas compared to non-KTCN samples. We were able to identify genomic regions with distinct patterns of DNA hypo- and hypermethylation and link them to previously found KTCN susceptibility loci as well as transcriptomic disruption of Wnt signaling pathway observed in KTCN. These findings provide novel insights into the KTCN etiology, allowing us to elucidate some of the common mechanisms behind the development of this complex and heterogenic disease. 

Acknowledgments

Supported by National Science Centre in Poland, Grant 2012/05/E/NZ5/02127 (to MG). 

Disclosure: M. Kabza, None; J.A. Karolak, None; M. Rydzanicz, None; M. Udziela, None; P. Gasperowicz, None; R. Ploski, None; J.P. Szaflik, None; M. Gajecka, None