Transcriptome Sequencing Reveals Tgf-β-Mediated Noncoding RNA Regulatory Mechanisms Involved in DNA Damage in the 661W Photoreceptor Cell Line

Transforming growth factor β (Tgf-β), a pleiotropic cytokine, can enhance DNA repair in various cells, including cancer cells and neurons. The noncoding regulatory system plays an important role in Tgf-β-mediated biological activities, whereas few studies have explored its role in DNA damage and repair. In this study, we suggested that Tgf-β improved while its inhibitor LSKL impaired DNA repair and cell viability in UV-irradiated 661W cells. Moreover, RNA-seq was carried out, and a total of 106 differentially expressed (DE)-mRNAs and 7 DE-lncRNAs were identified between UV/LSKL and UV/ctrl 661W cells. Gene ontology and Reactome analysis confirmed that the DE-mRNAs were enriched in multiple DNA damaged- and repair-related biological functions and pathways. We then constructed a ceRNA network that included 3 lncRNAs, 19 miRNAs, and 29 mRNAs with a bioinformatics prediction. Through RT-qPCR and further functional verification, 2 Tgf-β-mediated ceRNA axes (Gm20559-miR-361-5p-Oas2/Gbp7) were further identified. Gm20559 knockout or miR-361-5p mimics markedly impaired DNA repair and cell viability in UV-irradiated 661W cells, which confirms the bioinformatics results. In summary, this study revealed that Tgf-β could reduce DNA damage in 661W cells, provided a Tgf-β-associated ceRNA network for DNA damage and repair, and suggested that the molecular signatures may be useful candidates as targets of treatment for photoreceptor pathology.


Introduction
Photoreceptor cells (PRs) are principal light-sensitive, terminally differentiated neurons in the outer retina. Malfunction or loss of PRs is involved in multiple blindness-causing and intractable diseases, such as retinitis pigmentosa (RP) [1] and age-related macular degeneration (AMD) [2]. Notably, due to their constant exposure to light and metabolically active characteristics, PRs are rather susceptible to oxidative stress and DNA damage, which act as important factors of pathogenesis of the aforementioned diseases [3][4][5]. As such, to improve vision retention and disease prognosis, we sought to determine the mechanism of DNA damage/repair regulation in PRs to promote proper DNA repair and subsequent cell survival.
Neurons are highly susceptible to DNA damage, including single-strand breaks (SSBs) and double-strand breaks (DSBs). SSBs are among the most common lesions and can arise directly from attack or indirectly during the repair of base and nucleotide damage [6]. Meanwhile, DSBs are routinely formed in neurons under normal conditions [7]. γ-H2AX is at the center of a repair signaling cascade initiated by DSBs, and its locus formation is now
661W cells were divided into 5 groups: (1) UV group, in which a control vehicle was added before UV exposure (30 J/m 2 , 254 nm, irradiated for 40 s); (2) Tgf-β group, in which Tgf-β (Sinobiological, Beijing, China, final concentration at 1 ng/mL [19]) was added 1 h prior to UV exposure; (3) LSKL group, in which LSKL (MCE, HY-P0299, final concentration at 10 µg/mL) was added 24 h prior to UV exposure; (4) Tgf-β + LSKL, in which a pretreatment with Tgf-β for 1 h and LSKL for 24 h followed by UV exposure was added; and (5) Sham, in which there was no UV exposure. For UV exposure, the cell medium was removed and replaced with PBS, and after 40 s of UV radiation, the cells were reincubated with a fresh cell culture medium and allowed to repair in the dark.

Primary Retinal Neurons Culture and Treatment
Postnatal 1-3 days old C57/BL6 mice were sacrificed with an intraperitoneal injection of 4% chloral hydrate (Sigma-Aldrich, St. Louis, MO, USA). The eyes were then resected, and the retinas were separated. The separated retinas were incubated in a solution containing 0.125% trypsin for 15 min at 37 • C to dissociate the cells. Then, the trypsinization was terminated by adding DMEM (Gibco, Carlsbad, CA, USA) supplemented with 10% FBS (Gibco, Carlsbad, CA, USA). The cells were seeded on plates precoated with 0.01% poly-D-lysine (Sigma-Aldrich, St. Louis, MO, USA) at approximately 1 × 10 6 cells/mL. After cultured for 24 h, the culture media was changed to neurobasal media.
Primary retinal cells were divided into 2 groups: (1) UV group, in which a control vehicle was added 24 h before UV exposure (30 J/m2, 254 nm, irradiated for 30 s) and (2) LSKL group, in which LSKL (MCE, HY-P0299, final concentration at 10 µg/mL) was added 24 h prior to UV exposure. For UV exposure, the cell medium was removed and replaced with PBS, and after 30 s of UV radiation, the cells were reincubated with a fresh cell culture medium and allowed to repair in the dark.

Immunofluorescence Assay
Cells were seeded on coverslips in six-well plates and treated with different procedures. At the 2 h mark after UV radiation, cells on coverslips were fixed with 4% paraformaldehyde (Leagene, Beijing, China) for 15 min at room temperature, permeabilized with 0.1% Triton X-100 (Sigma, St Louis, MO, USA) for 10 min, and blocked with 5% BSA for 1 h. Afterwards, 661W cells were incubated with anti-phospho-H2AX (Ser139) antibody (clone JBW301, Millipore, Billerica, MA, USA) diluted in 5% BSA (1:100) overnight at 4 • C. Primary retinal neurons were incubated with anti-phospho-H2AX (Ser139) antibody (1:100) and MAP2 (Boster, Wuhan, China) antibody (1:100) overnight at 4 • C. Cells were then washed three times for 10 min each in PBS. The 661W cells were incubated with Alexa Fluor 555 anti-rabbit lgG (Cell Signaling Technology, Inc.) diluted in 5% BSA (1:500) for 1 h at room temperature, and the primary retinal neurons were incubated with both Alexa Fluor 555 anti-rabbit lgG (1:500) and Alexa Fluor 488 anti-mouse lgG (1:500). The cells were then washed again three times for 10 min each with PBS. Nuclei was visualized with DAPI staining and mounted on slides. Images were obtained using a Zeiss Axio Imager.Z2 microscope. Cells with more than 10 foci were recognized as γ-H2AX positive. The percentages of γ-H2AXpositive cells among DAPI-positive 661W cells or MAP2-positive primary retinal neurons were determined in samples from at least three independent experiments. Images from 8 randomly selected fields were used for analysis in ImageJ software. At least 3000 cells were quantified.

Comet Assay
The alkaline comet assay was performed using the comet assay reagent kits (Trevigen, Gaithersburg, MD, USA, 4250-050-K) according to the manufacturer's protocol. Briefly, 661W cells were washed with ice-cold PBS and suspended in cold PBS 30 min after UV radiation. Then, the cells were combined at 1 × 10 5 /mL with molten LMAgarose (at 37 • C) at a ratio of 1:10, and they were pipetted 50 µL to the CometSlide TM . The slides were placed at 4 • C in the dark for 30 min, and then the slides were incubated in 4 • C Lysis Solution overnight. The slides were then immersed in freshly prepared Alkaline Unwinding Solution for 20 min at room temperature. Then, the electrophoresis was performed with 4 • C Alkaline Electrophoresis Solution at 25 Volt for 20 min. The slides were immersed twice in dH2O for 5 min each and in 70% ethanol for 5 min. After the slides were placed at 37 • C for 15 min, 100 µL diluted Gel-Red (Beyotime Biotechnology, Shanghai, China) was added onto the gel area and stained for 30 min. Then, images were obtained using a Zeiss Axio Imager.Z2 microscope. The tail length and olive tail moment were calculated using comet assay software.

RNA-seq Analysis and Differential Expression Analysis
The 661W cells were pretreated with 10 µg/mL LSKL for 24 h and irradiated with UV (30 J/m 2 , 254 nm) for 40 s. The total RNA was extracted from the 661W cells with Tri Reagent (Sigma-Aldrich, St. Louis, MO, USA, T9424) 1 h after UV radiation. Sequencing was performed on an Illumina NovaSeq 6000 with paired-end 150 bp using a TruSeq SBS Kit v4-HS (Illumina, Inc, San Diego, CA, USA).
For gene differential expression analysis, we used R software and the R package biomaRT to annotate the obtained RNA-seq data (HT-seq counts) [20,21]. Afterward, differential mRNA and lncRNA expression was analyzed between the two groups using the R package edgeR with library sizes normalized by the median count ratios across transcripts [22,23]. The p-values were calculated based on the Poisson distribution and corrected with the Benjamini-Hochberg (BH) correction. A false discovery rate (FDR) < 0.05 and |log2FoldChange| > 1 were set as the criteria for DE-genes screening.
Volcano plots were generated with the R package ggplot2 to visualize gene expression differences [24], and a heatmap was also produced to illustrate the gene expression profile of DE-genes with the pheatmap package of R. The expression data matrix was row-normalized for each gene prior to the application of average linkage clustering.

Functional Enrichment and Protein-Protein Interaction Analysis
A GO enrichment analysis was carried out using the R clusterProfiler package in the mouse annotation database org.Mm.eg.db [25]. A pathway enrichment analysis based on the Reactome database was conducted with the R package ReactomePA. The p-values were calculated based on a hypergeometric distribution using the Benjamini-Hochberg (BH) method. The pathways with a corrected p-value < 0.05 and q-value < 0.2 were considered significantly enriched. The results were visualized using custom plotting functions, with DNA damage-related GO and Reactome results displayed.
Differentially expressed mRNAs (DE-mRNAs) were used to construct a protein-protein interaction (PPI) network with the STRING database (http://www.string-db.org, accessed on 10 February 2021) [26]. The minimum required interaction score was set as high confidence (>0.9), and the active interaction sources were chosen as follows: Textmining, Experiment, Database, and co-expression. The results were downloaded in the TSV format and imported into Cytoscape software for visualization and further analysis of the functional protein-protein interaction networks [27]. The cytoHubba plugin was used to analyze hub genes, and maximal clique centrality (MCC) was used for calculating [28].

Reverse Transcription-Polymerase Chain Reaction (RT-qPCR)
The total RNA was isolated with a TRIzol reagent (Invitrogen, Carlsbad, CA, USA) 1 h after UV radiation. The cDNA was synthesized from the total RNA with The Prime-ScriptTM RTPCR kit (Takara, Dalian, China) to detect mRNA and lncRNA. Mature miRNAs were converted into cDNA using reverse transcription performed by the RevertAid First Strand cDNA Synthesis Kit (Thermo Scientific, Waltham, MA, USA) with miRNA-specific primers. The RT-qPCR reaction was performed in triplicate. Expression levels were measured with RT-qPCR with the SYBR Green system (Roche, Indianapolis, IN, USA). β-Actin served as a reference control. Each sample was analyzed in triplicate. The primer pairs used in this article are listed in Supplementary Table S1.

Cell Viability Assayed by Cell Counting Kit-8 (CCK8) Assay
The cell viability was assessed with a CCK8 assay (Biosharp, Guangzhou, China) according to the manufacturer's protocol. Briefly, 5 × 10 3 cells were seeded into a 96 well plate. A 10 µL CCK8 reagent was added to each well 6 h after UV radiation. The cells were then incubated for 2 h at 37 • C. After incubation, the optical density (OD) value at 450 nm was measured. A cell-free well with the same volume of DMEM was used as a blank. Cell viability was determined by the ratio of the OD value of a treated culture to that of an untreated control.

siRNAs and miRNA Mimics
The siRNAs and miRNA mimics were provided by Tsingke Biological Technology company (TsingKe, Beijing, China). Details of the sequences are shown in Supplementary  Table S2.

Transfection
Cells were transfected with 50 nM siRNA or miRNA mimics for 48 h using Lipofectamine 3000 (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's protocol. The 661W cells were 70% confluent at the time of transfection and transfected in Opti-MEM (Gibco, Carlsbad, CA, USA). The serum-free medium was changed to DMEM (Gibco, Carlsbad, CA, USA) containing 10% fetal bovine serum 6 h after transfection.

Statistical Analysis
The data is expressed as the mean ± SD. The differences among mean values were evaluated with a two-tailed Student's t-test (for 2 groups) and analysis of variance (ANOVA; for >2 groups). All calculations and statistical tests were performed with SPSS (IBM SPSS statistics 26). For all analyses, p < 0.05 was considered significant.
To further confirm the effect of Tgf-β inhibition on DNA damage and repair in 661W cells, we performed an alkaline comet assay, which is an effective and reliable method for DNA damage detection in single cells. The 661W cells were pretreated with LSKL (10 µg/mL) or a control vehicle for 24 h before UV radiation (30 s). It was demonstrated by our previous study that a radiation dose of 30 J/m 2 (254 nm) did not cause immediate cell death or apoptosis of the 661W cells while most cells showed apoptosis or death 12 h after radiation [6]. Therefore, the cells were harvested 30 min after UV radiation to ensure that the results of the comet assay represented the degree of DNA damage instead of the level of cell apoptosis. As shown in Figure 1D-F, the tail length (UV, 22.65 ± 2.53; UV + LSKL, 29.66 ± 4.14, * p = 0.017) and olive tail moment (UV, 3.85 ± 0.34; UV + LSKL, 5.50 ± 0.74, ** p = 0.006) of the comets were significantly increased in LSKL-treated groups. These results support that Tgf-β inhibition by LSKL could impair DNA stability in UV-radiated 661W cells.

Enrichment Analysis and PPI Network Establishment Revealed Mechanisms and Pathways Associated with LSKL Prompting DNA Damage
The DE-mRNAs were further analyzed with GO and pathway enrichment. When restricting to DNA damage and repair-related functions, the following GO terms were identified: nucleotidyl transferase activity, GTP binding, GTPase activity, and regulation of nuclease activity ( Figure 3A). The DE-mRNAs were also significantly enriched in some DNA damage-and repair-related pathways in the Reactome, such as posttranslational protein phosphorylation, SUMOylation of DNA damage response and repair proteins, DNA damage bypass, termination of translesion DNA synthesis, nitric oxide-stimulated guanylate cyclase, and cGMP effects ( Figure 3B). Taken together, these results showed that the DE-mRNAs we obtained are closely related to the DNA damage and repair process.
DNA damage-and repair-related pathways in the Reactome, such as posttrans protein phosphorylation, SUMOylation of DNA damage response and repair p DNA damage bypass, termination of translesion DNA synthesis, nitric oxide-stim guanylate cyclase, and cGMP effects ( Figure 3B). Taken together, these results that the DE-mRNAs we obtained are closely related to the DNA damage and rep cess. To further investigate the interactions among the proteins encoded by the id DE-mRNAs, we conducted a PPI analysis. A PPI network with 28 nodes (27 dow lated and 1 upregulated) and 140 edges was obtained at the highest confidence le ( Figure 3C). The node size was mapped according to |log2FoldChange| values, lines between the nodes were mapped by interaction scores. As shown in Figure  lines between blue nodes are thicker, whereas the red nodes are only connected To further investigate the interactions among the proteins encoded by the identified DE-mRNAs, we conducted a PPI analysis. A PPI network with 28 nodes (27 downregulated and 1 upregulated) and 140 edges was obtained at the highest confidence level (0.9) ( Figure 3C). The node size was mapped according to |log2FoldChange| values, and the lines between the nodes were mapped by interaction scores. As shown in Figure 3C, the lines between blue nodes are thicker, whereas the red nodes are only connected to one blue node. These results suggest that the proteins of downregulated genes have closer and more complex interactions with each other.
The PPI network was further analyzed to identify important central genes. According to the calculation of the node connectivity with cytoHubba, the top 10 genes were selected as hub genes (see Table 2) and visualized in Figure 3D. Our results show that all the central genes are downregulated genes. Among them, Isg15, Usp18, Oasl2, and Oasl1 were all enriched in terms related to DNA damage. Usp18 1806 4 Irf7 1705 5 Rsad2 Oasl2 1561 7 Ifi44 848 8 Oasl1 720 9 Isg15 248 10 Stat1 134

Establishment and Verification of a lncRNA-miRNA-mRNA ceRNA Network
To further explore the potential noncoding regulatory mechanisms underlying LSKLmediated DNA damage, we established a ceRNA network through in silico prediction and experimental validation. First, we used an online database to predict the target miRNAs of the hub genes, the DE-mRNAs enriched in DNA damage/repair-related GO terms, and all DE-lncRNAs. The predicted results were processed in R, including intersection and filtration, and a final ceRNA network containing 3 lncRNAs, 19 miRNAs, and 29 mRNAs was obtained. We then separated the ceRNA network into up-and downregulated subgroups based on the expression profiles of DE-lncRNAs. The upregulated ceRNA network (up_ceRNA network) consisted of 2 lncRNAs (Gm20559 and A530040E14Rik), 16 miRNAs, and 28 mRNAs while the downregulated ceRNA network (down_ceRNA network) consisted of 1 lncRNA (Gm15564), 3 miRNAs, and 1 mRNA ( Figure 4A).

Discussion
In the present study, we found that Tgf-β attenuated DNA damage in PRs upon UV radiation while its inhibitor LSKL promoted DNA damage. Moreover, multiple downstream DNA damage/repair-related mRNAs and pathways of Tgf-β signaling were identified through transcriptome sequencing and bioinformatics analysis. Furthermore, 2 downstream lncRNA-miRNA-mRNA axes (Gm20559-miR-361-5p-Gbp7/Oas2) of Tgf-β signaling were demonstrated through the ceRNA network prediction followed by stepwise experimental verification, which revealed the noncoding mechanisms underlying the Tgfβ-regulated DNA damage and repair process in 661W cells. Taken together, our study confirmed the protective role of Tgf-β against DNA damage in PRs and further revealed its underlying mechanisms.
In our study, we first demonstrated that Tgf-β could attenuate DNA damage and improve cell viability of UV-radiated 661W cells ( Figure 1A-F), which is supported by previous studies [11,12,15,[32][33][34]. The general neuroprotective role of Tgf-β against various harmful stimuli has been demonstrated [35][36][37][38][39]. In terms of the role of Tgf-β and photoreceptor damage in ocular diseases, a recent study showed that Tgf-β signaling is vital for photoreceptor viability, while deletion of Tgf-β signaling resulted in thinner outer nuclear layer and accelerated retinal degeneration [38]. More importantly, Tgf-β exerts its neuroprotective role through direct impact on photoreceptors [35], which is also consistent with our study. The capability of maintaining DNA stability could partially accounts for the neuroprotective role of Tgf-β. For example, Tgf-β could induce Gadd45b expression, enhance DNA repair, and therefore protect retinal ganglion cells [11]. Moreover, we previously found that Tgf-β1 was significantly downregulated upon UV treatment in the 661W cells, and the inhibition of Tgf-β1 enhanced UV-induced DNA damage [6]. Although studies showed that Tgf-β could protect neurons through maintaining DNA stability, some studies have reached different conclusions. For example, Hicks et al. reported that Tgf-β potentiated ethanol-induced DNA damage in neural stem cells [32]. Moreover, Tgf-β pathway inhibition promotes DNA repair and improves hematopoietic stem cell (HSC) survival in Fanconi anemia [40]. Notably, inhibition of Tgf-β signaling results in elevated homologous recombination (HR) repair with a concomitant decrease in nonhomologous end-joining (NHEJ). As such, we speculate that differences in cell types (proliferative/nonproliferative) and forms of DNA damage might account for the discrepancies between different studies.
Moreover, the inhibition of Tgf-β signaling affected a variety of biological functions and pathways associated with DNA damage and repair, such as nucleotidyltransferase activity, regulation of nuclease activity, posttranscriptional modification of DNA damage response, and repair proteins and translesion DNA synthesis (TLS) (Figure 3A,B). Among these, nucleotidyltransferases are a family of enzymes that can add nucleotides to substrates such as DNA, and they are involved in DNA damage and repair [41]. For example, Tang et al. reported that DNA polymerase β (Pol β) could act as a nucleotidyltransferase and is involved in the repair pathway for single nucleotide base excision repair (BER) [42]. Second, nucleases are various enzymes that promote the hydrolysis of nucleic acids; these activities are key for most DNA repair systems, including BER, NHEJ, and HR [43]. Last, TLS is a DNA damage tolerance mechanism allowing cells to bypass DNA lesions, thus avoiding the collapse of replication forks [44]. Changes in these pathways further confirmed that Tgf-β is involved in DNA repair in PRs. Furthermore, we identified novel and important downstream molecular targets of Tgf-β signaling effectively through the combination of RNA-seq and bioinformatic analysis, which is one of the strengths of this study. Isg15, Usp18, Oasl2, and Oasl1 were identified as hub genes enriched in terms related to DNA damage and repair, which were all downregulated upon LSKL treatment ( Figure 3C,D). Isg15 is a ubiquitin-like protein (UBL) which has been recognized as an important regulator of genome stability through covalent modification of-or noncovalent interaction with-key proteins involved in the DNA damage response (DDR) [44]. Associated pathways include p53 signaling, replication fork acceleration, and TLS. Usp18 is also a key modulator of TLS. It acts as an Isg15-specific protease responsible for the removal of Isg15 from substrates and is responsible for resumption of DNA replication [45]. Downregulation of these two genes implies reduced TLS capability, resulting in unsuccessful DNA repair. Moreover, Oasl1 and Oasl2 are associated with the regulation of nuclease activity. Accordingly, they have 2 -5 oligoadenylate synthase (OAS)-like qualities, such as activation by interferon (IFN) and binding to dsRNA [46]. However, few studies have reported the relationship between Oasl and DNA damage and repair, which requires further study.
To comprehensively depict the interactions of both protein-coding and noncoding RNAs, we then established a ceRNA network through prediction and stepwise experimental validation ( Figure 4A-D). Several novel molecular targets associated with retinal DNA damage and their interactions were identified based on the ceRNA network. The combination of bioinformatic analysis and experimental validation improved the reliability of these results. We found that Gm20559 knockdown caused the release and upregulation of miR-361-5p, which then reduced the expression of its target mRNAs Oas2 and Gbp7 ( Figure 4E,F). Most importantly, our study indicated that Gm20559/miR-361-5p could stabilize/impair DNA stability in PRs, respectively ( Figure 5A-F), which is indirectly supported by previous studies. For example, it was demonstrated that miR-361-5p was associated with PARP1, ubiquitin protein ligase E3 component N-recognin 5 (UBR5) and ataxiatelangiectasia mutated interactor (ATMIN) in tumors, which are enriched in DNA damage and repair [47,48]. Of note is the fact that miR-361-5p upregulation could reduce UBR5 expression [47], therefore inhibiting ATMIN ubiquitination, attenuating the restoration of ATM and impairing DNA repair [49]. Among the other ceRNA components, Gm20559 is a lncRNA with 3 exons located on chromosome 6. Limited studies have shown that Gm20559 is associated with neuronal diseases, such as Alzheimer's Disease [50] and rabies virus infection of the mouse brain [51], while the biofunction of Gm20559 is largely unknown. This is the first study demonstrating that Gm20559 could regulate DNA repair in retinal 661W cells. Moreover, Gbp7 belongs to the guanylate-binding protein family, which is part of an IFN-γ-inducible guanosine triphosphatase (GTPase) superfamily [52]. GTPases are actively involved in various DNA repair processes [53,54], which might account for the regulatory role of Gbp7 in DNA repair. Similar to Gbp7, Oas2 could also be induced by IFN and serve a role in host defense against infection [55]. Previous studies supported the assertion that the upregulation of Oas2 could enhance resistance to DNA damage and improve cell viability in cancers [56,57]. However, few studies have discussed the function of Oas2 in neuronal pathologies. As such, further studies are required to explore the roles of Gm20559, Gbp7, and Oas2 in DNA damage and repair.
Of note is that inhibition of Tgf-β signaling also accelerated UV-induced DNA damage in primary retinal neurons, and similar expression profiles of Gm20559 and miR-361-5p were also identified. In addition, the general neuroprotective role of Tgf-β was supported by multiple studies as mentioned above. A great portion of downstream hub genes and ceRNA components of Tgf-β signaling were also shown to modulate DNA damage and repair in other tissues. Taken together, these findings and studies supported the general nature of our results. However, there are some limitations in this study. For example, this study was only performed in cell culture, and some experiments were performed on only one cell line. Thus, further studies are required to validate the effects of Tgf-β and its