In silico modeling of epigenetic-induced changes in photoreceptor cis-regulatory elements

Tissue collection

Deidentified post-mortem human donor eyes procured from three individuals with no reported ocular disease (National Disease Research Interchange, Philadelphia, PA) were used to collect ocular tissues for gene-specific quantification of DNA methylation. All donor tissue was collected within 8 h of donor death (Appendix 1). The tissue was shipped from the donor location submerged in saline on wet ice and received within 24 h of donor death (Appendix 1). A scalpel was used to pierce the limbus followed by collection of the cornea using scissors. Whole corneas were flash frozen, ground into a fine powder using a mortar and pestle, and then immediately transferred to nucleic acid extraction buffer and stored at −80 °C. Eyecups were further dissected with scissors by making four radial cuts exposing the retina. Cuts posterior to the ciliary margin were made liberating the retina from the anterior portion of the eye. The vitreous was removed, and the retina was carefully peeled away from the eyecup using fine forceps. Whole retina from each donor was briefly washed in Hank’s Balanced Salt Solution without calcium or magnesium to rinse away contaminating RPE cells, placed in nucleic acid extraction buffer, vigorously vortexed, and stored at −80 °C.

Nucleic acid purification

Genomic DNA was extracted from human ocular tissues using a Qiagen AllPrep Kit (Hilden, Germany) per the manufacturer’s instructions. Briefly, lysates were homogenized using QIAshredder spin columns. Lysate flow-through was then transferred to silica-based spin columns where genomic DNA and RNA were sequentially purified. The quality and quantity of the DNAs were examined using ultraviolet (UV) spectrophotometry.

Bisulfite pyrosequencing

Quantitative analysis of DNA methylation was measured using bisulfite pyrosequencing performed as previously described [14,21]. Briefly, bisulfite conversion was performed on 200 ng genomic DNA using the EZ DNA Methylation-Gold Kit (Zymo, Irvine, CA). Following conversion, 30 μl PCR reactions were performed using 2X JumpStart Taq ReadyMix (Sigma-Aldrich, St. Louis, MO). 5′ biotinylated PCR primers were designed to the 5′ regulatory regions of the target genes using PyroMark Assay Design 2.0 software (Qiagen, Hilden, Germany). Commercially available 5′ biotinylated PCR primers were used to amplify the 5′ promoters of LINE1 retrotransposon repeats (Qiagen). The PCR cycling conditions were 95 °C for 1 min, followed by 45 cycles of 95 °C for 30 s, 50–58 °C for 30 s, and 72 °C for 30 s, with a final extension at 72 °C for 1 min on a Bio-Rad C1000 Touch Thermal Cycler (Bio-Rad, Hercules, CA). Variable PCR annealing temperatures for different primer sets are indicated in Table 1. Biotinylated PCR products were purified and made single stranded to serve as a template in a pyrosequencing reaction using the PyroMark Q24 Vacuum Prep Tool (Qiagen) per the manufacturer’s instructions. A sequencing primer (0.3 μM final) was annealed to the purified single-stranded PCR product, and pyrosequencing reactions were performed using the PyroMark Q24 Pyrosequencing System (Qiagen) per the manufacturer’s instructions. Percent DNA methylation at each CpG dinucleotide in the bisulfite PCR amplicon was determined and averaged between biologic triplicates. Statistical significance between the two sample groups was determined using a one-tailed t test with the significance threshold set at 0.01. All PCR and sequencing primers used in these experiments are shown in Table 1.

Bioinformatics analysis of CRX binding sites

Computationally predicted CBRs in the human genome were generated from CBRs obtained from a previously published CRX chromatin immunoprecipitation study of mouse photoreceptors [19]. Genome coordinates from experimentally validated CBRs in the mouse mm9 genome assembly were downloaded, aligned, and mapped to the human hg19 genome assembly using the LiftOver tool available within the UCSC Genome Browser [22]. These data were then built as a custom track in the UCSC Genome Browser’s hg19 genome assembly viewer and were overlaid with adult human retina RNA-sequencing transcriptome data as determined by Farkas and colleagues [23], vertebrate conservation data [24], and custom tracks created in the UCSC Genome Browser highlighting predicted and experimentally validated CRX binding sites [25].

Homology modeling of human CRX binding site

The UniProt accession number O43186 was used for modeling and assembly of human CRX [26]. The DNA binding domain of CRX, consisting of amino acids 39–98, was generated using SWISS-MODEL [27]. Models of the CRX DNA binding domains were energy minimized in Yet Another Scientific Artificial Reality Application (YASARA) before the docking experiments [28,29].

Molecular modeling of the rhodopsin enhancer region and RHO promoter region

The human rhodopsin enhancer region (RER) and RHO (OMIM 180380; Gene ID: 6010) promoter sequences were obtained from UCSC Genome Browser and three-dimensional (3D) B-DNA structures generated using 3D-DART [30]. The 39-bp DNA construct for the RER consisted of the following bases: 5′-ACC TCA TTA GCG TTG GGC ATA ATC ACC AGG CCA AGC GCC-3′. The RHO promoter region consisted of the following 37 bps: 5′-TCT GCA GCG GGG ATT AAT ATG ATT ATG AAC ACC CC-3′. The DNA sequences were numbered in the 5′ to 3′ direction, continuing with the complementary strands. CpG methylation regions previously identified in this study were generated through YASARA. Methyl groups were built onto the C-5 position of the cytosine ring by using the “Swap>Atom” function. Molecular dynamics modeling was conducted using the AMBER14 force field in YASARA. The cutoff for electrostatic interactions was set at 10.5 Å. The cell boundaries were periodic and filled with 0.3% magnesium chloride and water molecules at pH 7.0. The simulation was run at 298.1 K. The save interval was every 0.1 ns over the 100 ns of the simulation. All other parameters remained in the default setting.

Root mean square deviation (RMSD), Root mean square fluctuation (RMSF), and groove widths were calculated in YASARA. Groove widths were calculated by measuring the distance between the phosphate of one nucleotide to the phosphate of a nucleotide on opposite face of the groove. For the CRX binding sequence calculations, groove widths were localized to nucleotides 20–25 for the RER and 15–19 for the RHO promoter. Calculations and statistics were performed in R v3.4.0, and data were plotted using ggplot2, ggridges v0.4.0, and cowplot v0.9.2.9900 [31,32].

CRX docking

CRX was bound to the RER and RHO sequences through High Ambiguity-Driven protein-protein DOCKing (HADDOCK) to determine the interaction mechanism [33,34]. The previously modeled CRX DNA binding domain structure consisting of residues 39–98 was used in the docking simulations with the energy minimized average structures from the promoter and RER simulations. Active residues for CRX include amino acids 40–46, 63, 69, 82, 84–85, and 88–93 [26]. Active residues for the DNA molecules were bases 20–25 for the RER and 15–19 for the RHO promoter. CRX was then docked to the sequences using the Prediction interface on the HADDOCK server. We then used YASARA to determine the interactions occurring between the residues and the consensus sequence.

DNA methylation is inversely correlated with gene expression of photoreceptor-specific genes in human ocular tissues

Previous reports in the literature have demonstrated an inverse correlation between DNA methylation and gene expression of photoreceptor-specific genes [11,12,35]. These studies have been limited, however, to murine model organisms or immortalized human cell lines and have yet to be characterized in primary human tissues. To determine whether this trend is also observed in human ocular tissues, corneas and retinas were collected from three sex-matched postmortem eyes procured from human donors 75 years or older (Appendix 1). DNAs extracted from these tissues were used for quantitative bisulfite pyrosequencing analysis of DNA methylation on 5′ regulatory regions upstream of the phototransduction genes RHO (Figure 1A) and PDE6B (Gene ID: 5158 OMIM: 180072; Figure 1B), the embryonic eye field transcription factor PAX6 (Gene ID: 5080; OMIM: 607108), and the multicopy long interspersed nuclear element (LINE1) retrotransposon repeats. The rod-specific genes RHO and PDE6B have been previously shown to have cell-specific patterns of gene expression in murine and human rod photoreceptors [11,12,36]. Previous human retina RNA-sequencing transcriptome data demonstrate that these rod-specific genes are highly transcribed in the adult human retina (Figure 1) [23]. Bisulfite pyrosequencing analysis of conserved 5′ regulatory regions upstream of RHO and PDE6B demonstrate lower levels of DNA methylation in retinal tissue relative to the cornea, a tissue in which neither gene is expressed (Figure 2A,B). Levels of global DNA methylation between retina and corneal tissues were determined to be similar based on measurement of the constitutively methylated and silent LINE1 retrotransposon repeats, as well as the unmethylated embryonic eye field transcription factor PAX6 (Figure 2C,D). Collectively, these data demonstrate an inverse correlation between DNA methylation and transcription of RHO and PDE6B in human ocular tissues. These findings are consistent with previous observations of epigenetic regulation of phototransduction genes in the mouse retina, as well as immortalized cell lines derived from human retinal tissue.

Differentially methylated regions upstream of photoreceptor-specific genes correspond to CRX binding sites

The epigenetic analysis of the human eye in this study demonstrates tissue-specific patterns of DNA methylation in 5′ regulatory sequences of photoreceptor-specific genes that are inversely correlated with mRNA expression. These and previous observations in other mammalian retinal model systems suggest a functional role for DNA methylation in repressing transcription at these loci. However, a mechanism for this repression remains uncharacterized. A commonality between RHO, PDE6B, and many other photoreceptor-specific genes is that they are transcriptionally regulated by the homeodomain transcription factor CRX [19,36,37]. This observation led us to question whether DNA methylation plays a role in regulating temporal and cell-specific binding of CRX to cis-regulatory regions within photoreceptor genomes. To test this hypothesis, we used a computational approach to align previously determined genome-wide CBRs in the rod-rich wild-type mouse retina to the human hg19 2009 genome assembly (Figure 1; Predicted CBRs). To further assess the functionality of these presumptive regulatory regions, predicted CBRs were searched for sequences containing CRX binding motifs in the human PDE6B upstream promoter region (Figure 1) or experimentally validated CRX binding sites in the human RHO locus (Figure 3). These analyses demonstrated that differentially methylated regions (DMRs) identified in this study are adjacent to experimentally validated CRX binding sites in the well-characterized RHO 5′ promoter and RER, as well as predicted CRX binding sites in the PDE6B 5′ promoter. Given these results, we predict that differential methylation of these CpG sites may play a prominent role in modulating CRX binding to cis-regulatory elements upstream of RHO, PDE6B, and other CRX-regulated genes.

The CRX DNA binding domain interacts with DNA grooves but not known methylation sites

DNA methylation is known to alter the three-dimensional structure of double-stranded DNA and interactions with DNA binding proteins, such as transcription factors [38,39]. Computational modeling data suggest that the presence of a bulky methyl group in 5mC results in widening of the major groove and a concomitant narrowing of the minor groove [39]. Further, it has been experimentally validated that 5mC DNA bases biochemically mimic thymine bases due to a similar hydrophobic interaction with major groove edge methyl groups [40]. More recent evidence has demonstrated a structural basis for affinity modulation of human transcription factors to methylated DNA via interactions of hydrophobic amino acid residues with methyl groups [38]. Taken together, these observations suggest a mechanism for DNA methylation to narrow minor groove width and increase available hydrophobic interactions of cis-regulatory elements as a reversible mode of modulating CRX binding. However, a lack of CRX structural information, as well as the interactions of CRX with DNA, has hindered development of a specific mechanism. To test the hypothesis that methylated DNA sequences are less favorable binding partners with CRX than their unmethylated counterparts, we generated a homology model of the human CRX DNA binding domain and models of the human RHO promoter and RER. The CRX binding domain model showed the expected three-helix bundle characteristic of a helix-turn-helix homeobox protein (Figure 4). In the promoter and RER, the modeled DNA duplexes spanned one experimentally validated CRX binding site, as well as one (promoter) or two (RER) CpG cytosine substrates for DNA methylation. We then docked the model of CRX to the DNA models using HADDOCK (Figure 4A,B) [34]. The CRX binding site within the promoter is within the minor groove, and the C-terminal helix of CRX fits into the groove making contact with the DNA backbone. Glutamic acid 80, lysine 88, and arginine 90 are known disease-causing mutations within CRX, and in our model of CRX bound to DNA, we saw that K88 and R90 make electrostatic interactions with DNA [41-43]. E80 makes contact with R69 and Q84, which contact the minor groove of DNA, suggesting a structural role for this amino acid (Figure 4C,D). Thus, we felt confident that this CRX model could provide insight into how CRX interacts with DNA and is regulated by methylation. We next repeated the docking to the CRX binding sequence in the RER. The best interaction score produced by HADDOCK placed CRX in the major groove side of the CRX binding sequence. In addition to the previous backbone interactions noted for the promoter, CRX makes interactions with the DNA bases via N89 and K93 (Figure 4C,D). Given that there currently is no structural information on CRX or these DNA sequences beyond the models we present here, we cannot exclude either binding mode. What is clear is that CRX likely interacts with a single DNA groove and does not make direct contact with methylation sites at this locus. These data suggest a mechanism of inhibition based on direct interaction with the CpG methyl sites is unlikely.

CpG methylation enhances DNA flexibility

The data showing that CpG methylation negatively regulates CRX binding raises the question of the molecular mechanism for binding inhibition. In addition to steric effects, CpG methylation have been proposed to change DNA dynamics and change the local structure of DNA; however, which of these scenarios applies to CRX was not clear [38,39,44]. To address this question of the inhibition mechanism, we created a second set of DNA models that were methylated at the known sites, and then, we equilibrated the models using molecular dynamics. We then analyzed the global and local effects of CpG methylation on the dynamics and structure of the DNA sequences and the CRX binding sites. Groove width statistics for these simulations are summarized in Table 2. Although the mean and median values for the groove widths are similar, global analysis of the major and minor grooves shows that the two structures are more dynamic (Table 2 and Figure 5). The increase in flexibility is indicated by the increased variance and standard deviations of the data measurements and visually by the wider distribution of values in the density plots (Table 2 and Figure 5). Analysis of the ensemble volume ratio shows that the structures of the methylated promoter have a larger volume, characteristic of higher dynamics. The effect of methylation is more apparent in the groove width measurement for the CRX binding sequence bases, which show wider major grooves and narrower minor grooves (Figure 5). These data show that the binding site is part of the dynamic region of DNA suggested by the global values calculated from the simulation. Moreover, the observations indicate that methylation does not alter the structure of DNA or the CRX binding site within the RHO promoter but changes the dynamics of the structure.

We then performed molecular dynamics simulations and analyses using the human RER sequence, which is methylated at sites on either side of the CRX binding motif. Methylation of the RER induced increased bending of the RER sequence in simulations relative to the unmodified sequence (Figure 6). Further analysis of the groove widths in the RER sequence showed that with methylation the major groove widens while the minor groove narrows relative to the unmodified form (Figure 6 and Table 2). The standard deviation and variance for the major groove also increase, indicative of increased dynamics. The ensemble volume ratio increased by about 25% suggesting drastic changes in the dynamics of the RER upon methylation. The effect of methylation is especially apparent in the density plots showing data for the CRX binding site (Figure 6). The bases in the major groove that include the CRX binding show a population of wider measurements that is not observed in the unmodified simulations. There are notable changes in the minor groove measurements upon methylation as well. Collectively, these data show that the presence of methylation alters the dynamics of DNA and the groove widths, possibly preventing efficient binding by CRX to the promoter and RER sequence.