Integrative RNA-seq and ATAC-seq Analysis of Retinitis Pigmentosa Caused by PDE6 Mutation


 Purpose: Inhibition of PARP1 could relieve PDE6 mutation-induced Retinitis pigmentosa (RP). However, the mechanism related with PARP1 overexpression in the RP has not been clarified. We attempted to explore the potential regulatory mechanism related with PARP1 underlying RP. Methods: ATAC-seq and RNA-seq were preformed for retina tissues of C3H and rd1 mice. The differential expressed genes (DEGs) were identified, followed by PARP1-DEG coexpression network, and PPI network construction. GO-BP and pathway enrichment of DEGs of interest were performed by clusterprofiler software. The overlapped genes that might play regulatory role in PARP1 expression were mined by integrated analysis of RNA-seq and ATAC-seq data. Results: Total 1061 DEGs were identified between C3H and rd1 group. Co-expression network was constructed with 313 PARP1-gene coexpression pairs. The down-regulated DEGs were closely related with visual perception, light stimulus related biological process, while the up-regulated DEGs were significantly enriched in phototransduction and PPAR signaling pathway. PPI network was constructed with 202 nodes and 375 edges, which was clustered into 3 modules. Module 1 genes were closely related with detection of light stimulus, visual perception related biological process and phototransduction pathway (involved with Gnat1/Guca1b/Gnat2/Sag/Pde6g). By integrated analysis of the RNA-seq and ATAC-seq, the overlapped up-regulated genes were Asxl3 and Nyap2, while the down-regulated genes were Tmem136 and Susd3. Conclusion: Gnat1 may play a key role in RP development by interacting with PARP1. Susd3 may play a regulatory role in PARP1 expression and affect RP formation.


Introduction
Retinitis pigmentosa (RP) is a genetic retinal neurodegenerative disease, characterized by the pathogenic gene mutation and progressive apoptosis of retina photoreceptor cells. RP is one of the most common causes of blindness and severe visual decline [1]. It is estimated that the incidence rate of RP is about 1/4000 worldwide [2] and it mostly common affects people aged 20 to 60 years [3]. Despite the advances in the treatment measures for RP, such as gene therapy [4], stem cell therapy [5], there is no cure for RP currently. RP gene mutation is responsible for majority of RP cases, which is associated with retinal molecular pathway defects, leading to the nal photoreceptor cell apoptosis. Thus, targeting the prevention of photoreceptor cell apoptosis may shed light on the effective therapy for RP.
Many efforts have been made to the deep understanding of the photoreceptor cell apoptosis underlying RP. It is reported that PDE (Phosphodiesterases) and PARP1 (poly-ADP-ribose polymerase) play a key role in the apoptosis process of photoreceptor cell [6]. PDE plays a regulatory role in signal transduction. The PDE defects-induced by RP gene mutation could cause the accumulation of cGMP, which posed toxicity to photoreceptor cell [7]. PDE6 mutation is responsible for 4-8% RP in human homologous animal models [7,8]. In addition, PARP1 exerts various biological functions in cell apoptosis and gene expression regulation [9]. PARP-1 is involved in the chromatin function regulation by interacting with transcriptional factor (TF), DNA methyltransferas 1, and histone modifying enzymes [10,11]. Previous evidence has revealed that PARP1 is proposed as the therapeutic target for intervening retinal photoreceptor function and survival in RP induced by PDE6A mutation [12], indicating the causal relationship between PARP1 dysregulation and PDE6 gene mutations.
The previous study for the progression of retinal degeneration in RP mice model with PDE6A mutation suggested that the apoptotic photoreceptors cells caused by PDE6 mutation peaked at postnatal 15 day (P15) [13]. RP caused by PDE6 dysfunction is primarily studied in rd1 mice with Pde6b gene mutations [14]. Therefore, in this study, we performed RNA-seq and ATAC-seq for retina tissues from PDE6B mutant mice (Rd1) and C3H wild type (wt) mice.
The differential accessible peaks and differential expressed genes were analyzed. We aimed to mine the regulatory mechanism underlying PARP1 overexpression and photoreceptors cells apoptosis.

Animals and sampling
The C3H wild type (wt) mice during the rst 31 postnatal days and age-paired Rd1 (Pde6b rd1) mice with Pde6b mutation were presented by Tubingen university. All the speci c pathogen free (SPF) animals were maintained in the IVC cages with food and water freely. The mice were acclimatized to experimental environment under a 12:12 light/dark cycle at 22°C. The C3H mice and Rd1 mice at P15 were sacri ced, respectively. Then, the mice retinal tissues were extracted and the single cells were isolated according the previous description [12]. Our study was approved by the Ethics Committee of Yunnan University (No. YNUCARE20210024) and Yunnan University a liated Hospital (No. 2021045). All the animal experimental procedures conformed to Animal Care guidelines.

ATAC sequencing and processing
For ATAT-seq, 50,000 single cells were prepared and collected by centrifugation at 500 g for 5 min at 4°C. Cell nuclei extraction was achieved with cold lysis buffer followed by transposase reaction with Nextera kit. After puri ed, DNA was ampli ed by PCR analysis in a 50 µl reaction mixture with the reaction condition of 1 cycle of 72°C for 5 min, 98°C for 30 s, followed by 5 cycles of 98°C for 10s, 63°C fro 30s and 2°C for 2 min. After puri ed, the generated DNA libraries was pooled and subjected to Illumina HiSeq. The quality of raw data was controlled by FastQC. The pair-end reads were mapped to reference genome sequences by the bwa program. Peak calling was performed by MACS2 software with q value <0.05. The differential peaks (open chromatin regions) between C3H and Rd1 samples were analyzed by DEGseq with the cutoff value of |log2FC| > 1 and P value < 0.05. The differential Peak adjacent genes were annotated by annotatePeaks.pl tool in HOMER.

RNA sequencing and processing
Total RNA was extracted from retinal tissues and the quality control was achieved by Nanodrop and Agilent 2100.
After rRNA was removed, RNA samples were fragmented by fragmentation buffer. The rst strand of cDNA was synthesized with random hexamers, followed by the second strand synthesis with dNTPs, DNA polymerase and RNase H. The cDNAs were subjected to end repair, A tail addition, and the second strand (containing U) degradation. The cDNA library was constructed by PCR ampli cation and sequenced based on HiSeq platform. The raw data were ltered by removing the reads with adapter, >10% N bases and low quality. The clean reads were mapped to reference genome based on modi ed BWT algorithm by Hisat2 software. The gene expression levels were analyzed by reads count located at genomic regions or exon regions of genes according to FPKM method with the application of HTSeq software. The differentially expressed genes (DEGs) between C3H and Rd1 samples were screened by DESeq2 software with cutoff value of logFC > 0.585 and adjust P value < 0.05.

Coexpression of PARP1 with DEGs
The expression correlation between PARP1 and DEGs was estimated by Pearson correlation coe cient (r). The absolute value of r ranged from 0 to 1. When the absolute of r was closed to 1, the most signi cant correlation between the PARP1 and DEG pair was considered. The coexpression gene pairs were collected with |r| >0.9 and p value <0.05, followed by coexpression network construction.

Protein-protein interaction (PPI) network and module analysis
The protein interaction pairs were retrieved from STRING [15]  Integrated analysis of RNA-seq and ATAC-seq data The analysis of ATAC-seq data unravels transcription factor binding sites. To explore the regulatory role of TFs in gene expression, ATAC-seq and RNA-seq data were integrated. The down-regulated genes obtained by RNA-seq data were overlapped with those around Peaks with depressed reads signals analyzed by ATAC-seq data, otherwise the up-regulated genes were overlapped. In the present study, the differential Peak adjacent genes analyzed by ATACseq data were overlapped with nodes in co-expression network. The overlaps were subjected to functional enrichment analysis.

Identi cation of DEGs
The genes that may be perturbed by Pde6b mutation and play a key role in the development of photoreceptor cells in the retina were identi ed by RNA-seqbn n . Total 1061 genes were found to be differentially expressed in Rd1 samples compared with C3H ones at P15, including 340 up-regulated genes and 721 down-regulated ones. The gene differential expression was illustrated by heatmap, which showed differential expression pro le of DEGs between Rd1 and C3H samples ( Figure 1A). The volcano plot represented differential expressed genes based on the de ned p value and fold change cutoff values ( Figure 1B).

Coexpression network of PARP1 with DEGs
Based on Pearson correlation analysis, 313 PARP1-DEG coexpressed pairs were obtained, including 102 pairs with positive correlation and 211 ones with negative expression. The coexpression network was illustrated in Figure 2, which was comprised with 102 up-regulated genes and 211 down-regulated genes. Function enrichment analysis was performed for up-and down-regulated genes, respectively. Results showed that the up-regulated genes were closely related with tube formation related biological process, PPAR signaling pathway and B cell receptor signaling pathway ( Table 1). The down-regulated genes were signi cantly enriched in visual perception and sensory perception of light stimulus related biological process. The most signi cant pathway for down-regulated genes was phototransduction (mmu04744) ( Table 1).

PPI network and the functional modules
The PPI network was constructed with 202 nodes connected with 375 edges ( Figure 3A). Three signi cant network modules were divided by MCODE with PPI score ≥ 5, including Module1(score=8.222, 10 nodes and 37 edges), Module 2(score=6.333, 7 nodes and 19 edges) and Module 3 (7 nodes and 18 edges) ( Figure 3B). After topological property analysis, the node degrees of modular genes were obtained ( Table 2 Table 3).

ATAC-seq and integrated analysis
Based on ATAC-seq, we obtained 1394 peaks up-regulated signals and 891 differential peaks with down-regulated signal. The differential accessible peaks covered genes were overlapped with DEGs in co-expression network. The Venn diagram depicted that the overlapped up-regulated genes were Asxl3 and Nyap2 ( Figure 4A, while the downregulated overlaps were Tmem136 and Susd3 ( Figure 4B).

Discussion
RP is one of the leading cause of inherited blindness, which has affected 1/3000 to 1/5000 individuals globally.
PDE6B mutation is responsible for 5-8% RP patients, which affects rod-speci c cyclic guanosine monophosphate (cGMP) level and results in rod degeneration [21,22]. Recent evidence showed that PARP1 was a promising therapeutic target for PDE6 mutation-induced RP. However, the regulatory mechanism of PARP1 in photoreceptor cell apoptosis underlying RP has not been clari ed. In this study, we attempted to explore the potential mechanism related with photoreceptor cell apoptosis in PDE6 mutation-induced RP with ATAC-seq and RNA-seq.
Our results of RNA-seq showed that there were 1061 DEGs between Rd1 and C3H mice at P15. The expression pro le of DEGs could distinguish the samples from Rd1 and C3H clearly, which revealed that DEGs identi ed were signi cant and could play key roles in the pathogenesis of RP. In order to explore the mechanism related with PARP1 in RP, PARP1-DEG coexpression network was constructed. Functional analysis showed that the downregulated genes with coexpression correlation with PARP1 were closely related with visual perception, sensory perception of light stimulus related biological process and phototransduction pathway. The manifestations of RP were characterized by loss of depth perception, night blindness, and loss of central vision [23]. The variety of eye disorders in RP patients was related with rod-opsin protein disruption, in parallel with the defects of phototransduction cascade of central nervous [24]. Our ndings were consistent with the reports mentioned above and revealed that the visual perception and phototransduction pathway were dysregulated in RP. Additionally, our data showed that PPAR signaling pathway was a signi cant pathway enriched by up-regulated genes. It is reported that PPAR plays a regulatory role in the transcription of genes involved in intracellular signaling pathways [25].
PPAR is involved in lipid metabolism and homeostasis [26]. Accumulating evidence indicates that abnormal lipid metabolism is associated with RP development [26,27]. In this study, the PPAR pathway was found to be enhanced by up-regulated genes. Thus, we speculated that PARP1 could be associated with abnormal lipid metabolism induced photoreceptor cell apoptosis in RP.
Besides, the PPI network indicated that the Module 1 genes were also enriched in phototransduction pathway, which was in agree with GO-BP analysis of DEGs in coexpression network. Gnat 1 was a down-regulated gene in retinal tissues of Rd1 mice with highest degree of 19 in Module 1 PPI network. As described in previous study, Gnat1 encodes the α-subunit of transducin in rod cells. Gnat1 variants were previously reported to be associated with congenital stationary night-blindness [28]. New evidence showed that Gnat1 mutation was responsible for the slow and late onset of retinal degeneration [29]. Similarly, our GO-BP enrichment analysis showed that Gnat1 was closely related with visual perception and light stimulus related biological process. Gnat1 showed negative expression correlation with PARP1 (r=-0.90 and p = 0.013). All these above suggested that the down-regulation of Gnat1 was correlated with PARP1 overexpression, which caused the dysfunction of visual perception, light stimulus, and phototransduction, leading to the onset of RP.
Furthermore, the RNA-seq and ATAC-seq integrated analysis in our study revealed that Susd3 (sushi domaincontaining protein 3) was a down-regulated overlap. Susd3 is cell surface protein and has been suggested as the potential biomarker for breast cancer [30]. Besides, Susd3 was found to be one of the genes encoding protein at the surface of young rod photoreceptors [31]. Susd2 was identi ed to be specially expressed for photoreceptors in the retina. However, the signi cant role of Susd3 in RP onset has not been clari ed. ATAC-seq is a simple protocol for open chromatin analysis, which provides possibility of discovering TF binding sites. Susd3 was the gene covered by differential peaks by ATAC-seq, which indicated that Susd3 could exert TF function. Susd3 showed coexpression correlation with PARP1. Currently, little is known about whether Susd3 elicited transcriptional regulation on PARP1 and affected PDE6-mutation induced RP, which should be veri ed in the following studies.
In summary, the down-regulated genes with coexpression correlation of PARP1 were signi cantly enriched in visual perception, sensory perception of light stimulus related biological process and phototransduction pathway. The PPAR signaling pathway was enhanced by up-regulated genes. Gnat1 could play a key role in visual perception, light stimulus, and phototransduction of RP by interacting with PARP1. Susd3 might play a regulatory role in PARP1 expression and affect RP formation. Our study may provide new insights in understanding the mechanism of PDE6mutation induced RP.   Protein-protein interaction network and its modules A, PPI network of DEGs were constructed with 202 nodes and 375 edges. B, PPI network was clustered into 3 modules, such as module-1(score=8.222), module-2(score=6.333) and module-3(score=6). Red triangle represents as PARP1, yellow circle indicates up-regulated genes, green prismatic indicates down-regulated genes. The node size is based on degree value.