Discovery of GLO1 New Related Genes and Pathways by RNA-Seq on A2E-Stressed Retinal Epithelial Cells Could Improve Knowledge on Retinitis Pigmentosa

Endogenous antioxidants protect cells from reactive oxygen species (ROS)-related deleterious effects, and an imbalance in the oxidant/antioxidant systems generates oxidative stress. Glyoxalase 1 (GLO1) is a ubiquitous cellular enzyme involved in detoxification of methylglyoxal (MG), a cytotoxic byproduct of glycolysis whose excess can produce oxidative stress. In retinitis pigmentosa, one of the most diffuse cause of blindness, oxidative damage leads to photoreceptor death. To clarify the role of GLO1 in retinitis pigmentosa onset and progression, we treated human retinal pigment epithelium cells by the oxidant agent A2E. Transcriptome profiles between treated and untreated cells were performed by RNA-Seq, considering two time points (3 and 6 h), after the basal one. The exposure to A2E highlighted significant expression differences and splicing events in 370 GLO1 first-neighbor genes, and 23 of them emerged from pathway clustered analysis as main candidates to be associated with retinitis pigmentosa. Such a hypothesis was corroborated by the involvement of previously analyzed genes in specific cellular activities related to oxidative stress, such as glyoxylate and dicarboxylate metabolism, glycolysis, axo-dendritic transport, lipoprotein activity and metabolism, SUMOylation and retrograde transport at the trans-Golgi network. Our findings could be the starting point to explore unclear molecular mechanisms involved in retinitis pigmentosa etiopathogenesis.


Introduction
In recent years, extensive evidence has associated oxidative stress to the pathophysiology of many human diseases, and this phenomenon has now become a key focus of translational research. Oxidative stress was defined as an unbalanced level of reactive oxygen species (ROS) and intrinsic antioxidant defenses [1]. ROS act as second messengers in numerous intracellular signaling cascades, helping to maintain cellular homeostasis [2]. However, their excess can cause undiscriminating damage to biological molecules, inducing loss of function and even cell death [3]. In addition to ROS, high levels of advanced glycation end products (AGEs), heterogeneous molecules produced from nonenzymatic

Cell Culture
Human Retinal Pigment Epithelial Cells (H-RPE, Clonetics™, Lonza, Walkersville, MD, USA) were cultivated, and then grown for 24 h to reach confluence, as previously described [24]. Subsequently, a group of cells was treated with A2E 20 µM for 3 and 6 h before rinsing with medium, while a control group was incubated without the oxidant agent. Finally, confluent cultures were transferred to PBS-CMG and, then, subjected to blue light delivered by a tungsten halogen source (470 ± 20 nm; 0.4 mW/mm 2 ) for 30 min, in order to induce phototoxicity of A2E, and incubated at 37 • C for 24 h. Each group of cells had three biological replicates.

MTT Assay
The mitochondrial-dependent reduction of methylthiazolyldiphenyl-tetrazolium bromide (MTT) (Sigma-Aldrich, St. Louis, MO, USA) to formazan insoluble crystals was realized to analyze cell viability, following a protocol already described [24]. Finally, a Dynatech microplate reader permitted evaluation of the absorbance at 570 nm, and results were expressed as a percentage of viable cells compared with control conditions in the absence of A2E. Multiple t-tests were performed for statistical comparison (p-value < 0.05), considering 3 independent experiments, each one characterized by 3 replicates.

Total RNA Isolation and RNA-Seq Profiling
Total RNA was extracted, checked for degradation and contamination, and quantified as previously reported [24]. The RNA-seq samples were divided in 3 factor groups, consisting of human RPE cells before A2E treatment and at the different time points of 3 and 6 h, respectively. Each group was biologically replicated three times, for a total of 9 samples. Both 3 and 6 h time points were chosen on the basis of previous experiments realized by our research group (unpublished data), confirmed by results obtained from the MTT assay in this work. Such outcomes showed that in wider time intervals, the death rate of oxidative stressed cells might be so high as to invalidate the following expression analysis. Libraries were generated using 1 µg of total RNA by the TruSeq Stranded Total RNA Sample Prep Kit with Ribo-Zero H/M/R (Illumina, San Diego, CA, USA), following the manufacturer's protocols. The last step involved the sequencing of the libraries on an HiSeq 2500 Sequencer (Illumina, San Diego, CA, USA), using the HiSeq SBS Kit v4 (Illumina, San Diego, CA, USA).

Data Analysis
Obtained raw reads were quality trimmed and then mapped against the hg38 reference genome and Ensembl RNA database v.99 (EMBL-EBI, Hinxton, Cambridgeshire, UK) by Qiagen CLC Genomics Workbench v.20.0.2 (Qiagen, Hilden, Germany) [25], following criteria adopted by Donato et al. [24]. Aligned reads were quantified by a mapping-dependent expectation-maximization (EM) algorithm [26], and the transcript per million reads (TPM) values were then computed after normalization by the trimmed mean of M-values (TMM) method [27]. Differential expression analysis was realized using the Limma R package [28], setting the contrast groups as 0 (untreated) versus 3 h (treated), 0 (untreated) versus 6 h (treated), 3 versus 6 h (both treated) and ((0 h.untreated+3 h.treated)/2)-((0 h.untreated+6 h.treated)/2). The latter resulted from a multiple group mean comparison, that permitted the indirect detection of the differences in expression level induced by the whole period of treatment, hereafter called "Due to Time" effect. For differentially expressed genes, the log 2 fold change (log 2 FC) of their abundance was calculated based on contrast groups, and significance of expression changes was determined using the t-test [29]. Benjamini and Hochberg (BH) post hoc test was then applied to correct false discovery rate (FDR) on p-values obtained by multiple testing [30]. The genes uniquely identified in the RPE cells, presenting at least 3 unique gene reads, lower than two-fold (log 2 FC < −1, down-regulated) or greater than two-fold (log 2 FC > 1, up-regulated) changes in expression, and with BH-adjusted p-values lower than 0.05, were chosen for functional classification.

Quantitative RT-PCR (qRT-PCR) Validation
To assess the reliability of RNA-Seq data, ten of most dysregulated genes from the whole RNA-Seq analysis were selected to be validated by quantitative Real-Time polymerase chain reaction (qRT-PCR). Reverse transcription was performed with a GoScript™ Reverse Transcription System (Promega, Madison, WI, USA), according to manufacturer's protocol. The qRT-PCR was, then, applied on obtained cDNA in the ABI 7500 fast sequence detection system (Applied Biosystems, Foster City, CA, USA), using the BRYT-Green based PCR reaction, as previously performed [24]. Each reaction was replicate 6 times, considering all analyzed time points (18 samples), and the average threshold cycle (Ct) was calculated for each replicate. Gene expression was normalized to the expression level of most stable reference gene, identified as combination of GeNorm [31], BestKeeper [32], Delta Ct [33] and NormFinder [34] algorithm results. The relative gene expression was then estimated using the 2 −∆∆Ct method, and the results were shown as the mean ± SEM (Standard Error of Mean). The analysis of variance between groups (ANOVA), corrected by Bonferroni post hoc test, permitted the assessment of the statistical significance. Lastly, a linear regression analysis was carried out to check the correlation of fold change ratios between RNA-Seq and qRT-PCR. The whole statistical analyses were executed using IBM SPSS 26.0 software (https://www.ibm.com/analytics/us/en/technology/spss/). Corrected p-values < 0.05 were considered as statistically significant.

Pathway Analysis
GO term enrichment analysis for the most dysregulated genes was performed using the GeneMANIA (v. 3 [37]. GeneMANIA has been set to find the top 20 related genes and at most attributes using automatic weighting. ClusterMaker2 performed a clustering based on BestNeighbor Filter, set with the Proportion of node edges in cluster = 0.5. ClueGO options have been set as follow: CLINVAR, GO (Biological Process, Cellular Component, Molecular Function and Immune System Process), INTERPRO, KEGG, REACTOME (Pathways and Reactions), WIKIPATHWAYS and CORUM 3.0 as selected ontologies; GO Tree Interval Min Level = 3 and Max Level = 8; GO Term/Pathway Selection Min # Genes = 3 and % Genes = 4.000; GO Term/Pathway Network Connectivity (Kappa Score) = 0.4; Statistics Options set on Enrichment/Depletion (Two-Sided hypergeometric test), with pV correction = Bonferroni step-down. CluePedia was used following default settings. Finally, only GO terms with p < 0.01 were selected.

A2E Treatment Highlighted a Significant Negative Effect on RPE Cell Viability
The MTT assay showed a relevant impact of A2E treatment on RPE cell survival in a time-dependent manner. In contrast to the control group, the viability of treated RPE cells was significantly decreased (p < 0.05), especially after 6 h from treatment ( Figure 1).

A2E Treatment Highlighted GLO1 Down-Regulation, as Well as the Most of Its Related Genes
The RNA sequencing globally generated about 100 million quality reads (mean mapping quality=29) and with a percentage of ~70% uniquely mapped. A total of 58,243 differentially expressed genes (DEGs) were detected out of 59,661 reference counterparts, considering the whole human genome annotations. All previous mapping statistics were based on average values calculated for all three replicates in each time point. Among all detected DEGs, 24,465 showed expression alterations in evaluated time points ( Figure S1), and 370 of them were highly correlated to GLO1, as an output of GENEMANIA analysis. In detail, 244 were down-regulated, with the lowest expression value reached by PIP5K1B (log2FC = −2.207, p-value = 0.000), and 126 were up-regulated, with the highest value shown by ARID5A (log2FC = 1.918, p-value = 0.000) ( Table S1). The subsequent clustering by ClusterMaker2 shed light on a particular group of genes, the best neighbors of GLO1, showing the same expression trend ( Figure 2). Of this cluster consisted of 23 genes, 20 were globally down-expressed and only 3 were over-expressed. In detail, the lowest expression value was shown precisely by GLO1 (log2FC = −1.088, p-value = 0.000), while the highest expression change resulted from ANKH (log2FC = 0.564, p-value = 8.656E-05) ( Table 1).

A2E Treatment Highlighted GLO1 Down-Regulation, as Well as the Most of Its Related Genes
The RNA sequencing globally generated about 100 million quality reads (mean mapping quality = 29) and with a percentage of~70% uniquely mapped. A total of 58,243 differentially expressed genes (DEGs) were detected out of 59,661 reference counterparts, considering the whole human genome annotations. All previous mapping statistics were based on average values calculated for all three replicates in each time point. Among all detected DEGs, 24,465 showed expression alterations in evaluated time points ( Figure S1), and 370 of them were highly correlated to GLO1, as an output of GENEMANIA analysis. In detail, 244 were down-regulated, with the lowest expression value reached by PIP5K1B (log 2 FC = −2.207, p-value = 0.000), and 126 were up-regulated, with the highest value shown by ARID5A (log 2 FC = 1.918, p-value = 0.000) ( Table S1). The subsequent clustering by ClusterMaker2 shed light on a particular group of genes, the best neighbors of GLO1, showing the same expression trend ( Figure 2). Of this cluster consisted of 23 genes, 20 were globally down-expressed and only 3 were over-expressed. In detail, the lowest expression value was shown precisely by GLO1 (log 2 FC = −1.088, p-value = 0.000), while the highest expression change resulted from ANKH (log 2 FC = 0.564, p-value = 8.656E-05) ( Table 1).

qRT-PCR Validation
To confirm the reliability of DEGs identified by deep sequencing, a total of ten genes were chosen among most the dysregulated and pathway-related to GLO1, for confirmation in a biologically independent experiment using qRT-PCR. qRT-PCR analysis showed that the expression trends of the chosen genes matched with those observed by the RNA-Seq (Table 2 and Figure 3), although there were some little differences in the degree of the changes. The analysis of variance (ANOVA) method, conducted to compare the means among multiple groups, highlighted high significance (p-values < 0.05).

qRT-PCR Validation
To confirm the reliability of DEGs identified by deep sequencing, a total of ten genes were chosen among most the dysregulated and pathway-related to GLO1, for confirmation in a biologically independent experiment using qRT-PCR. qRT-PCR analysis showed that the expression trends of the chosen genes matched with those observed by the RNA-Seq (Table 2 and Figure 3), although there were some little differences in the degree of the changes. The analysis of variance (ANOVA) method, conducted to compare the means among multiple groups, highlighted high significance (p-values < 0.05).

Positive regulation of mTOR signaling SIK3
Nuclear import of proteins IPO7

Catalysis of mucin-type oligosaccharides GALNT10
PPi transport regulation ANKH

Discussion
The development of eye-related dystrophies is based on both environmental and genetic features. A common pathogenic factor is the accumulation of AGEs which are the toxic byproducts of the nonenzymatic reaction of proteins, lipids and nucleic acids with reducing sugars [38]. Advanced glycation (or glycosylation) occurs in all cytotypes and regards the post-translational modification of glucose-derived dicarbonyl compounds and amino residues present in proteins, lipids and DNA [39]. In the retina, AGEs could promote vascular dysfunction by altering intra-and extracellular protein structure and by increasing inflammation and oxidative stress [40]. The retinal AGE deposition could determine an upregulation of vascular endothelial growth factor (VEGF) and a downregulation of pigment epithelium-derived factor (PEDF), as well as a serious disruption of the inner blood-retinal barrier (iBRB) [41]. Additionally, boosted advanced lipoxidation end-products (ALEs) accumulation was also detected in the outer retina. This portion is mainly rich in polyunsaturated fatty acids that are highly susceptible to lipid peroxidation, further contributing to retinal aging and diseases [42].
Intracellular AGEs are able to induce post-translational modifications of regulatory proteins, especially those in the proteasome, and can directly bind and impair mitochondrial proteins involved in the electron transport chain, inhibiting oxidative ATP production [43]. Intracellular AGE precursors such as methylglyoxal (MG) and glyoxal (GO) can also modify and inhibit the function of important enzymes such as Glyceraldehyde-3-Phosphate Dehydrogenase (GAPDH) and glyoxalase 1 (GLO1) [44]. Glyoxalase 1, in particular, is part of the endogenous detoxification system that converts toxic AGE-forming dicarbonyls to less reactive products [45]. In detail, GLO1 metabolizes MG and prevents MG-induced damage. An excess of MG inactivates antioxidant enzymes such as glutathione peroxidase and SOD enzymes, impairing degradation of MG and determining a positive feedback loop [46]. The activity of GLO1 gradually declines during aging induces the accumulation of AGE in tissues [47]. GLO1 overexpression was already associated with protective effects against neuroglial and vascular lesions [48], and its mutations could be involved in cancer onset [49], as well as in several genetic pathologies like cerebral cavernous malformations (CCMs) [50] and retinitis pigmentosa [11]. Even if it is known that GLO1 defects could play a role in retinitis pigmentosa, contributing to the accumulation of AGEs in the retina, very little is known about molecular mechanisms and pathways by which GLO1 determines its deleterious effects.
In order to achieve this purpose, we applied next generation sequencing (NGS) technologies to investigate the whole transcriptome of RPE cells during a follow-up of two time points (3 and 6 h) after exposure to activated oxidant compound A2E. Oxidative stress is currently recognized as one of the most relevant biochemical pathways involved in RP etiopathogenesis, especially targeting the high metabolic demand of RPE cells. In addition to the down-expression of GLO1, we found several GLO1-related genes dysregulated in treated cells. The analysis of these genes shed light on candidate pathways, shared with GLO1, probably responsible for the cellular damage phenotype.
Among these pathways, the one related to cytoskeleton dynamics and RhoGTPases activation cascades needs to be mentioned. Cytoplasmic Activator of Transcription and Developmental Regulator AUTS2 activates the Rho family small GTPase Rac1, a key coordinator of actin polymerization and microtubule dynamics, controlling neuronal migration and neurite extension. The nuclear form of AUTS2 acts as a transcriptional activator of several target genes, in complex with the polycomb complex 1 (PRC1) [51]. The other two genes encoding for proteins acting as RhoGTPase regulators are Rho GTPase Activating Protein 21 (ARHGAP21) and Protein Tyrosine Phosphatase Non-Receptor Type 13 (PTPN13), both involved in the down-regulation of cell migration and proliferation, cell polarity, cell adhesion, stress fiber formation and cell differentiation. ARHGAP21 also plays a relevant role in Golgi regulation and positioning, intracellular trafficking and glucose homeostasis [52]. PTPN13, instead, is a tyrosine phosphatase which mediates phosphoinositide 3-kinase (PI3K) signaling trough dephosphorylation of Phosphoinositide-3-Kinase Regulatory Subunit 2 (PIK3R2) [53] and regulates negative apoptotic signaling [54], as the antagonist of RhoGTPase A, SLIT-ROBO Rho GTPase Activating Protein 1 (SRGAP1) [55]. The latter is also involved in the modulation of contractility during epithelial junction maturation [56], a role also played by Formin Like 2 (FMNL2), probably by regulating actin polymerization and organization of the cytoskeleton [57]. Thus, the global down-expression of these genes (the only up-regulated gene is AUTS2) could indicate the alteration of actin filament structure and activity, leading to impairment in cell polarity and adhesion, with the final result of an increase in RPE apoptosis.
Moreover, FMNL2 expression reduction, already known to be involved in glaucoma [58], causes huge Golgi dispersal, malformations of vesicular organelles and defective anterograde transport from the Golgi to plasma membrane [57]. The same consequences could occur with down-expression of Ubiquitin C (UBC) [59], Myosin XVIIIA (MYO18A) and Epidermal Growth Factor Receptor Pathway Substrate 15 (EPS15), and with over-expression of ANKH Inorganic Pyrophosphate Transport Regulator (ANKH). In particular, MYO18A is also involved in intracellular transport processes, in retrograde treadmilling of actin and in its transport from focal adhesions to the leading edge [60]. ANKH regulates the transmembrane efflux of ATP and the trans-Golgi network trafficking, as well as endocytosis [61]. The latter process is also modulated by EPS15, whose encoded protein is involved in clathrin-coated pit maturation, including invagination or budding [62]. Furthermore, EPS15 is involved in cell growth regulation, synaptic vesicle recycling and recruitment of alpha-adaptin [63]. Thus, the dysregulation of these genes could affect the normal vesicular trafficking of RPE cells, fundamental for many retinal processes, such as POS renewal and visual cycle intermediate regeneration, as well as avoiding AGE accumulation.
In order to reduce the deposition of such compounds, the ubiquitin-proteasome system represents an intracellular system that has to be fully functional. Three down-regulated GLO1-related genes, Ring Finger and FYVE Like Domain Containing E3 Ubiquitin Protein Ligase (RFFL), F-Box and WD Repeat Domain Containing 2 (FBXW2) and Cullin Associated and Neddylation Dissociated 1 (CAND1), play an important role inside this system. RFFL is an E3 ubiquitin-protein ligase that directly regulates cell migration through the mTORC2 complex and negatively regulates cell death downstream of death domain receptors in the extrinsic pathway of apoptosis [64,65]. FBXW2 is also an E3 ligase, which promotes the ubiquitylation and degradation of ß-catenin, up-regulated in the WNT/ß-catenin pathway as ca onsequence of inflammatory processes. FBXW2 functions as a substrate recognition receptor in the SCF (Skp1/Cullin/F-box) E3 ubiquitin ligase complex, involved in ubiquitin-dependent proteasomal degradation of cyclins and key metabolite enzymes [66], also promoted by CAND1. Additionally, CAND1 is able to regulate the tubular endoplasmic reticulum network, remodeling through elongation and retraction of tubules [67]. Thus, the down-expression of these three genes could impair the ubiquitin-proteasome system activity, favoring the accumulation of misfolded proteins and AGEs, inducing ER stress.
Interestingly, both FBXW2 and CAND1 play a role in glycolytic metabolism. Through the WNT/ß-catenin pathway, FBXW2 normally activates aerobic glycolysis (also called the Warburg effect), triggering the glycolytic enzymes Glut/PDK1/LDH-A/MCT-1 via PI3K/Akt/HIF-1α. This phenomenon allows the retinal cells to metabolize glucose and protect them against oxidative damage [68]. A similar metabolic effect is exerted by CAND1 by controlling the abundance of the glycolysis-promoting enzyme 6-Phosphofructo-2-Kinase/Fructose-2,6-Biphosphatase 3 (PFKFB3). This is a substrate for the SCF complex, acting during the S phase of cell cycle. Thus, the presence of PFKFB3 is tightly controlled to ensure the up-regulation of glycolysis at a specific point of the G1 phase [69]. Moreover, PFKFB3 maintains a high level of glycolytic intermediates (e.g., for synthesis of nonessential amino acids or hexosamine), whose regulation highlighted the importance of anaplerosis by glycolysis to overcome the restriction point of the G1 phase [70]. This scenario suggests that the down-expression of both FBXW2 and CAND1 could impair the glycolytic metabolism in RPE cells, altering the following cellular respiration process and inducing ROS production and cell death.
A weakening of glucose metabolism might also be determined by down-expression of Glutamine-Fructose-6-Phosphate Transaminase 1 (GFPT1), the first and rate-limiting enzyme of the hexosamine biosynthetic pathway (HBP), involved in ubiquitous glycosylation processes [71]. HBP activation results in the synthesis of UDP-N-acetylglucosamine (UDP-GlcNAc), a donor of N-acetylglucosamine (GlcNAc) for O-linked and N-linked glycosylation of a wide range of proteins, especially ones involved in signal transduction and the TGF-ß pathway [72]. The reduction of glycosylation seems to be compensated by over-expression of Polypeptide N-Acetylgalactosaminyltransferase 10 (GALNT10) which, however, only catalyzes the initial reaction in O-linked oligosaccharide biosynthesis, such as mucin-type oligosaccharides [73]. The opposite trend of GFPT1 and GALNT10 expression probably suggests a serious induced ER stress and cell-cell and cell-extracellular matrix attachment impairment, due to a decrease of GFPT1-related N-glycosylation, balanced by an increased O-glycosylation by GALNT10, which could be interpreted as a final attempt to guarantee mitochondrial respiration and redox homeostasis, as already demonstrated in retinal aging [74]. A defect in energy metabolism might also depend on dysregulation of SIK Family Kinase 3 (SIK3), which encodes for a specific kinase that positively regulates mTOR and CREB signaling, as well as cholesterol biosynthesis by coupling with retinoid metabolism and melanogenesis [75]. Thus, the observed down-expression of SIK3 could decrease mitochondrial respiration and up-regulate removal of dysfunctional cellular components via autophagy, compromising cellular antioxidant mechanisms [76,77].
Intriguingly, a cluster of down-expressed GLO1-related genes is involved, directly or indirectly, in the regulation of translation machinery and cell survival. Decrease of Importin 7 (IPO7) could trigger p53 activation and p53-dependent growth arrest, also determining ribosomal biogenesis stress and nucleolar morphology changes [78]. Reduced levels of Mitochondrial Ribosomal Protein S33 (MRPS33) might damage mitochondrial protein synthesis [79]. MORC Family CW-Type Zinc Finger 4 (MORC4) and Microcephalin 1 (MCPH1) deregulation promotes apoptosis and arrests DNA damage repair [80,81]. Nuclear Factor I A (NFIA) down-regulation impairs mitotic exit and cell differentiation [82]. Cap Binding Complex Dependent Translation Initiation Factor (CTIF) down-expression alters mRNA and protein quality control, compromising the crosstalk between translation and the aggresome-autophagy pathway [83].

Conclusions
We report data obtained from whole transcriptome analysis performed on RPE cells, after exposure to A2E. Particularly, we focused on fold changes of genes sharing pathways with GLO1, comparing differential gene expression among control group cells and cultures at two different time points (3 and 6 h). We identified 22 GLO1-related genes that change their expression, and these are involved in a complex network of biochemical mechanisms that might be associated to RP onset and progression. Such pathways regarded microtubules and actin assembly, ubiquitin-proteasome activity, RE and Golgi integrity, vesicular trafficking, transcriptional and translational regulation, glycolytic metabolism control and glycosylation modifications.
However, obtained results could become more significant when an in vivo experiment will confirm what was observed in RPE cells, but with the wider scenario of the whole retina in a model organism. Moreover, other experimental procedures, such as immunoblots to validate obtained results at the protein level, will further clarify the relationship correlation among GLO1 and these related genes here discussed.
Despite the aforementioned limitations, our results could represent an important step towards clarification of new GLO1-related molecular mechanisms behind RP etiopathogenesis.
Supplementary Materials: The following are available online at http://www.mdpi.com/2076-3921/9/5/416/s1, Figure S1: Volcano plots based on p-values and fold-changes produced by Limma, Table S1: RNA-Seq differential expression analysis of the most dysregulated genes related to GLO1, Table S2: ClueGo and CluePedia pathway analysis of identified GLO1-related genes, Table S3: Pathways analysis of the 22 GLO1 most related neighbor genes.

Conflicts of Interest:
The authors declare no conflict of interest.