Molecular Vision 2018; 24:153-164
Received 18 August 2017 | Accepted 13 February 2018 | Published 15 February 2018
The first two authors contributed equally to this manuscript.
1School of Optometry, Department of Optometry and Vision Science, University of Alabama at Birmingham, Birmingham, AL; 2School of Medicine, Civitan International Research Center, University of Alabama at Birmingham, AL; 3Department of Chemistry and Biochemistry, Bates College, Lewiston, ME; 4School of Medicine, Department of Pharmacology, Vanderbilt University, TN
Correspondence to: Alecia K. Gross, University of Alabama at Birmingham, 1825 University Blvd, SHEL 906, Birmingham, AL 35294-0010; Phone: (205) 975-8396; FAX; (205) 975-7394; email: firstname.lastname@example.org
Purpose: Epigenetic and transcriptional mechanisms have been shown to contribute to long-lasting functional changes in adult neurons. The purpose of this study was to identify any such modifications in diseased retinal tissues from a mouse model of rhodopsin mutation-associated autosomal dominant retinitis pigmentosa (ADRP), Q344X, relative to age-matched wild-type (WT) controls.
Methods: We performed RNA sequencing (RNA-seq) at poly(A) selected RNA to profile the transcriptional patterns in 3-week-old ADRP mouse model rhodopsin Q344X compared to WT controls. Differentially expressed genes were determined by DESeq2 using the Benjamini & Hochberg p value adjustment and an absolute log2 fold change cutoff. Quantitative western blots were conducted to evaluate protein expression levels of histone H3 phosphorylated at serine 10 and histone H4. qRT-PCR was performed to validate the expression patterns of differentially expressed genes.
Results: We observed significant differential expression in 2151 genes in the retina of Q344X mice compared to WT controls, including downregulation in the potassium channel gene, Kcnv2, and differential expression of histone genes, including the H1 family histone member, H1foo; the H3 histone family 3B, H3f3b; and the histone deacetylase 9, Hdac9. Quantitative western blots revealed statistically significant decreased protein expression of both histone H3 phosphorylated at serine 10 and histone H4 in 3-week-old Q344X retinas. Furthermore, qRT-PCR performed on select differentially expressed genes based on our RNA-seq results revealed matched expression patterns of up or downregulation.
Conclusions: These findings provide evidence that transcriptomic alterations occur in the ADRP mouse model rhodopsin Q344X retina and that these processes may contribute to the dysfunction and neurodegeneration seen in this animal model.
One of the overarching goals in the treatment of blinding diseases is to determine the underlying pathological mechanisms that lead to blindness. Mutations in genes encoding retinal trafficking proteins often manifest as blinding diseases, such as retinitis pigmentosa (RP). RP is a heterogeneous group of hereditary disorders that lead to the progressive loss of retinal function. RP is the most common inherited disease leading to blindness, with a worldwide prevalence of 1 in 3500 people; roughly 25% of these cases are autosomal dominant retinitis pigmentosa (ADRP) [1-3]. To date, more than 25 genes are known to cause ADRP, and over 1000 mutations have been reported in these genes . Transgenic mouse, rat and fly models of ADRP have been constructed based on human rhodopsin mutations, and in mice apoptotic cell death has been shown as the route of retinal degeneration [2,5-7]. Symptoms include reduction in the peripheral visual field, leading to tunnel vision and total blindness. In patients, disease onset typically begins in the early teenage years, and severe visual impairment occurs between the ages of 45 and 60 . Typically, mutations in the C-terminus of the rhodopsin gene cause an earlier onset time course of degeneration in patients . One such mutation, Q344X, is a naturally occurring rhodopsin ADRP mutation that causes a truncation of the protein, effectively removing the trafficking signal VAPA and leading to rhodopsin mislocalization and subsequent cellular apoptosis [2,9]. While it has been shown that mutations in genes encoding retinal trafficking proteins affect processes vital for cell function and overall homeostasis, the underlying pathological mechanisms remain largely unknown.
In recent years, more studies have investigated the genetic and environmental influences associated with the onset and progression of retinal diseases, analyzing how these risk factors contribute to molecular alterations that ultimately lead to pathology [10-12]. Regulated gene expression is essential for proper cellular function and homeostasis. One method cells use to control gene expression is the structural alteration of the chromatin complex consisting of DNA, RNA, and protein [13,14]. Mechanisms that cause modifications of chromatin structure without changing the nucleotide sequence are known as epigenetic mechanisms, which trigger alterations in gene expression [10,15-17]. Histone modifications and DNA methylation are the two most extensively investigated epigenetic marks and contain distinct mechanisms affecting the repression or activation of transcription, which is dependent on factors such as the type of modification, genomic locations, and site specificity [10,14,18]. Recent studies have shown that these epigenetic marks change within the organism throughout its lifetime, finding that both histone subunit composition and DNA methylation are dynamically regulated in cells [16,19-21]. Nevertheless, it is evident that DNA methylation and attendant changes in chromatin structure are capable of self-perpetuation and self-regeneration. Disruption of these mechanisms can contribute to neuron and neural circuit dysfunction.
While it is evident that transcriptional and epigenetic processes are linked to disease and impairment in non-retinal tissue, there is lack of gene profiling studies of RP. Indeed, previous studies have investigated the transcriptome and novel transcripts in the retina [22,23], including the transcriptome from the autosomal-recessive RP model, rd10 . However, to date, the ADRP model rhodopsin Q344X has not been profiled using whole transcriptome techniques. In this study, we analyzed the transcriptomic alterations in the mRNA from 3-week-old mouse retina in the rhodopsin Q344X ADRP model relative to age-matched wild-type (WT) mouse retina to identify the potential targets of epigenetic modifications. This time point was specifically chosen because it has already been shown in this model that outer nuclear layer thickness and outer segment length are maintained and display no structural characteristics of retinal remodeling . Furthermore, to investigate whether epigenetics may have a role in driving and perpetuating persisting changes within rhodopsin Q344X, we examined the levels of protein expression from two histone proteins, histone H3 phosphorylated at serine 10 and histone H4.
Homozygous knock-in mice expressing human rhodopsin Q344X (Q344X mice) have been previously described . The mice were kept on a 12-h light/dark cycle. All procedures were performed using Institutional Animal Care and Use Committee (IACUC)-approved protocols and were conducted in full compliance with the Association for the Assessment and Accreditation of Laboratory Animal Care (AAALAC).
For RNA sequencing (RNA-seq), the retinas were extracted from 3-week-old mice (n=40 Q344X; 20 females, 20 males and n=34 WT; 17 females, 17 males; samples were collected at approximately noon) and placed in stabilization reagent (RNAlater, QIAGEN, Hilden, Germany; Catalog no. 76,104). Total RNA was extracted (RNeasy, QIAGEN; Catalog no. 74,104) and quality controlled (Bioanalyzer, Agilent, Santa Clara, CA). A single WT replicate contained a RIN number <7, which was removed from the study. Poly(A) selection was performed, and the libraries were constructed using the NEBNext Directional RNA Library kit at the HudsonAlpha Institute for Biotechnology. Multiplex sequencing was performed with the Illumina HiSeq 2500 system (HudsonAlpha, v4 sequencing reagents, paired end, 50 bp). All samples contained a minimum of 26 million paired end reads with an average number of 26.4 million reads across all replicates.
The FASTQ files were uploaded to the University of Alabama at Birmingham’s high-performance computer cluster, Cheaha, for bioinformatics analysis. First, the quality and control of the reads were assessed using FastQC, and trimming of the bases with quality scores of less than 20 was performed with Trim_Galore! (v 0.4.4). Following trimming, the reads were aligned with STAR  to the mm10 Ensembl genome (v 2.5.2a. During runMode genomeGenerate, the option sjdbOverhang was set to 49). The alignment resulted in an average of 92.3% of reads that were uniquely mapped. The total number of reads that were uniquely mapped per sample is shown in Appendix 1. Gene-level counts were generated from the binary alignment map (BAM) files using the featureCounts  function in the Rsubread package (v 1.26.1) in R with the Mus_musculus.GRCm38.90.gtf file from the Ensembl database. The options used included: isGTFAnnotationFile=TRUE, useMetaFeatures=TRUE, isPairedEnd=TRUE, requireBothEndsMapped=TRUE, strandSpecific=2, and autosort=TRUE; additional options were kept as default. A summary of the read assignment per sample is given in Appendix 1. Further, the gene body coverage from the BAM files was evaluated with RSeQC (v 2.6.3) and is reported in Appendix 1.
Finally, DESeq2  (v 1.16.1) in R was used to perform count normalization and differential expression analysis. Following count normalization, principle component analysis (PCA) was performed and sample-to-samples distances were computed using the Euclidean distance. The analysis identified a replicate from the Q344 × group as an outlier, which was removed from further analysis (Appendix 2). Therefore, the final data set for the current study includes two Q344X replicates and two WT replicates.
Differential expression was computed in DESeq2 with the application of the Benjamini & Hochberg false discovery rate (FDR) method to adjust the p values from multiple testing. Genes were called differentially expressed genes (DEGs) if they passed a statistical cutoff of FDR p<0.05 and if they contained an absolute log2 fold change (FC) >=1. Functional annotation enrichment analysis was performed in the NIH Database for Annotation, Visualization and Integrated Discovery (DAVID, v 6.8) by submitting all DEGs identified. The Benjamini & Hochberg FDR correction was also applied to determine gene ontology (GO) terms and KEGG pathways with the cutoff of an FDR p<0.05.
The FASTQ files for the current study have been uploaded to NCBI’s Gene Expression Omnibus under accession number GSE102247. Figures, including the heatmaps, scatter plot, volcano plot, and box plots were made in R with the following packages: ggplot2 (v 2.2.1), ggforce (v 0.1.1), pheatmap (v 1.0.8), and gplots (v 3.0.1). The heatmap of the DEG used the Euclidean clustering method shown in the dendrograms.
Retinas from mice (n=5 Q344X, 3 females, 2 males; n=5 WT, 3 females, 2 males; samples were collected at approximately noon) were homogenized in sample application buffer as previously described , separated by standard PAGE, transferred to a nitrocellulose-supported membrane (GVS North America, Sanford, ME; Catalog no. 1,212,590), blocked for 1 h at room temperature with 4% non-fat dry milk in TBST [20 mM Tris-Cl, pH 7.6 and 0.1% (v/v) Tween-20] supplemented with 0.02% sodium azide, and then incubated overnight at 4 °C in the same solution containing primary antibody (anti-phospho-histone H3 (Ser10) cloneE173, Millipore, Billerica, MA; Catalog no. 04–1093; anti-histone H4, Millipore; Catalog no. 7–108; anti-histone H3, Millipore; Catalog no. 06–599). Binding of secondary antibody conjugated to horseradish peroxidase (Invitrogen, Carlsbad, CA; Catalog no. 65–6120) was detected using chemiluminescent reagents by exposure to film. For quantification, individual band intensity was measured using ImageJ. Each band was normalized to its relevant loading control.
Retinas of Q344X and WT mice (n=12 Q344X; 6 females, 6 males and n=12 WT; 6 females, 6 males; samples were collected at approximately noon) were placed in stabilization reagent (RNAlater, QIAGEN; Catalog no. 76,104). Total RNA was extracted (RNeasy, QIAGEN; Catalog no. 74,104) and reverse transcribed using a High-Capacity cDNA Reverse Transcription Kit (Thermo Fischer Scientific, Waltham, MA; Catalog no. 4,368,814). Individual gene assays were purchased from Applied Biosystems, Foster City, CA for each of the RNAs analyzed. ΔΔCt values were generated using Hist1h2be (Mm01166416_s1), Hist2H4 (Mm01952224_u1), Hist3h2a (Mm01701417_s1), Kcnv2 (Mm00807577_m1), Pax7 (Mm01354484_m1), and Prss33 Mm00617657_m1). TaqMan gene assays with Actb (Mm02619580_g1) served as internal standards. qPCR results are shown as the average of three different amplifications of cDNAs that were generated. The ΔΔCT method was applied to determine relative cDNA levels . Unpaired Student t tests were conducted on ΔΔCT values from each genotype to determine their significance.
Statistical analyses for quantitative western blots were performed in JMP software (SAS Institute Inc.). Statistical analyses for qRT-PCR were performed in R (The R Foundation). For all studies where an n is reported, the n represents the number of separate animals. A Student’s unpaired t test was used, and significance was set at p<0.05.
RNA-seq profiling of the ADRP mouse model Q344X relative to WT mice shows that 2151 genes are differentially expressed (FDR p<0.05, absolute log2 FC >=1); 1386 genes were downregulated, and 765 genes were upregulated in Q344X (Figure 1A,B). As expected, RNA-seq indicates the downregulation of the rhodopsin gene, Rho, in Q344X animals (Figure 1C). It is interesting to note that H1foo, a H1 histone family gene that is commonly expressed in early development, was identified among the top upregulated genes in Q344X animals, indicating that epigenetic mechanisms may have a role in RP. Further, additional histone genes were identified, including the histone deacetylase Hdac9, the H3 histone gene H3f3b, and the H4 histone gene Hist2h4 (Figure 1C). A complete list of all histone genes is given in Table 1. In addition, Kcnv2, a voltage-gated potassium channel gene found in high levels in the retina, was found to be downregulated in Q344X. Additional DEGs found that are functionally related to transcriptional regulation include Egr1 (early growth response 1), Taf4b (TATA-box binding protein associated factor 4b), and Med20 (mediator complex subunit 20), suggesting that Q344X may alter specific targets affecting the transcriptional mechanism within this animal model. All genes found to be differentially expressed are represented in Figure 2, Appendix 3, and their normalized counts are shown in Appendix 4.
To identify genes that are functionally related among our DEGs, functional annotation clustering was performed on gene ontology (GO) terms using the Benjamini & Hochberg correction. Of the 2151 DEGs submitted to DAVID, 1927 were identified in the database. The analysis indicates transcriptomic alterations in the Q344x model are linked to the following terms: extracellular matrix (GO:0031012, 92 genes, FDR p=2.7×10−10), metal ion binding (GO:0046872, 394 genes, FDR p=1.6×10−6), visual perception (GO:0007601, 34 genes, FDR p=2.7×10−5), and photoreceptor outer segment (GO:0001750,17 genes, FDR p=2.7×10−3). A KEGG pathway for phototransduction (12 genes, FDR p=1.3×10−4) was also identified. The list of genes within the visual perception and photoreceptor GO terms and the genes within the KEGG pathway are presented in Table 2.
To determine the specificity of transcriptional alterations in the ADRP model, Q344X, the current results were compared to the published list  of DEGs in the autosomal-recessive RP model, rd10. The results indicate that 250 genes were present in both lists, and 196 genes were expressed in the same direction of change (up or downregulated) in both RP models. Of these, 178 genes were downregulated and 18 were upregulated in Q344X and rd10 mice. Therefore, 1955 genes (1901 genes not detected to be DEGs in rd10 + 54 genes with a dissimilar direction of change) in the rhodopsin Q344X mouse retina represent transcriptional alterations specific to Q344X. The genes that were found to be common between both models are indicated in Appendix 3.
To further investigate the role that histones may have in ADRP, we quantified the levels of two well studied histone proteins, histone H3 phosphorylated at serine 10 and histone H4. Quantitative western blots reveal statistically significant (p<0.05) decreased expression of both histone H3 phosphorylated at serine 10 and histone H4 (Figure 3A–D, Appendix 5). These data are associated with the decreased expression found in the RNA-seq results, indicating the proper function of these proteins is slightly inhibited in Q344X retinas.
In addition to Rho, which was confirmed to be downregulated in Q344X animals in our RNA-seq, we performed qPCR on reverse-transcribed RNA (qRT-PCR) isolated from the retinas of 3-week-old Q344X and WT mice. The downregulated genes tested were: potassium channel, subfamily V, member 2, Kcnv2 (p=0.0008); histone cluster 2 H4, Hist2H4 (p=0.032); and histone cluster 3 H2a, Hist3h2a (p=0.017). The upregulated genes tested were: protease, serine 33, Prss33 (p=0.0023) and paired box 7, Pax7 (p=0.034). All genes tested were found to be expressed in the same direction of change (down or upregulated), matching the results found from our RNA-seq data (Appendix 6).
We have shown that transcriptional mechanisms contribute to retinal degeneration in the ADRP mouse model rhodopsin Q344X. Through RNA-seq profiling, we have identified 2151 DEGs linked to ADRP by using a stringent statistical cutoff (FDR p<0.05). Intriguingly, our results show an increase in the RNA levels of the histone subunit H1foo, a subunit typically only expressed early in development. Meanwhile, there was a downregulation in the potassium channel gene, Kcnv2, which has also been shown to be downregulated in the rd10 model . Furthermore, transcriptional regulatory genes were also identified, including the upregulation of Egr1, which has also been linked to RP in the rd10 model and in additional models of retinal disorders, including the retinal degeneration slow mouse model and the retinoschisin knockout model [24,30-32]. Furthermore, the results suggest that we found not only genes related to the visual perception, photoreceptor outer segment, and phototransduction pathway but also a vast majority of genes related to the extracellular matrix, neuron differentiation, and neuron projection [33-37]. Given the abundance of work indicating the presence of epigenetic modifications in neuro-specific genes, the current findings suggest that epigenetic mechanisms may have a role in RP [13,14,21,38-42].
Interestingly, relative to the transcriptomic profile from the autosomal-recessive RP model, rd10, our results show an overlap of 196 genes with the same directional change . This suggests that the overlapping genes are ideal candidates for studying RP across distinct genotypic models. The data also indicate that there is specificity in transcriptional patterns in relation to the heterozygous Q344X knock-in mouse model that may serve as candidates in relation to the autosomal dominant form of RP.
Our data indicate transcripts from several histone family genes were downregulated in Q344X. To characterize the role that histones may have in ADRP at the protein level, we performed quantitative western blots to measure the levels of the proteins from post-translationally modified histone H3 phosphorylated at ser10 and the chromatin core histone H4. The results indicate that both histones are downregulated to a large degree. Post-translational modification, phosphorylation of serine 10 of the N-terminal arm of histone H3, has been shown to be fundamental for mitotic chromosomal segregation and condensation; however, it has also been shown to play a role in regulating transcription [43-45]. Other studies have shown in vivo that several residues within histone H3 and H4 cores are vital for heterochromatin integrity . Therefore, these results strongly support the disruption of the histone core particle in ADRP and that it potentially drives transcriptional changes that contribute to retinal degeneration and its severity.
These data support the idea that the histone subunit composition of the chromatin particle can be modified in response to an ADRP-inducing gene mutation and that chromatin remodeling may be a crucial regulator of lasting functional change in the retina in ADRP. Individual histone variants differ in their capacity to support specific transcriptional modifications, with some likely more associated with increased epigenetic “plasticity” and others more closely allied to epigenetically stable genomic regions . Understanding the molecular mechanisms of transcriptional modifications within retinal degenerations will provide pivotal information for future therapeutic interventions, thus improving the quality of life for patients.
Appendix 3. List of DEG in Q344x model with similarities in rd10 indicated.
Appendix 5. Statistical analyses of normalized protein expression level.
Appendix 6. qRT-PCR on select differentially expressed genes.
This work was supported by the National Eye Institute (R01EY019311, AKG), the Civitan Institute Research Center Precision Medicine Pilot Grant (AKG), and the National Institute of Mental Health (R01MH57014, JDS). We thank Dr. Scott Wilson for qPCR materials and insight, Dr. Jeremy Day for antibodies and helpful discussions, and Evan Boitet for technical support.