Transcriptome Analyses of lncRNAs in A2E-Stressed Retinal Epithelial Cells Unveil Advanced Links between Metabolic Impairments Related to Oxidative Stress and Retinitis Pigmentosa.

Long non-coding RNAs (lncRNAs) are untranslated transcripts which regulate many biological processes. Changes in lncRNA expression pattern are well-known related to various human disorders, such as ocular diseases. Among them, retinitis pigmentosa, one of the most heterogeneous inherited disorder, is strictly related to oxidative stress. However, little is known about regulative aspects able to link oxidative stress to etiopathogenesis of retinitis. Thus, we realized a total RNA-Seq experiment, analyzing human retinal pigment epithelium cells treated by the oxidant agent N-retinylidene-N-retinylethanolamine (A2E), considering three independent experimental groups (untreated control cells, cells treated for 3 h and cells treated for 6 h). Differentially expressed lncRNAs were filtered out, explored with specific tools and databases, and finally subjected to pathway analysis. We detected 3,3’-overlapping ncRNAs, 107 antisense, 24 sense-intronic, four sense-overlapping and 227 lincRNAs very differentially expressed throughout all considered time points. Analyzed lncRNAs could be involved in several biochemical pathways related to compromised response to oxidative stress, carbohydrate and lipid metabolism impairment, melanin biosynthetic process alteration, deficiency in cellular response to amino acid starvation, unbalanced regulation of cofactor metabolic process, all leading to retinal cell death. The explored lncRNAs could play a relevant role in retinitis pigmentosa etiopathogenesis, and seem to be the ideal candidate for novel molecular markers and therapeutic strategies.


Introduction
In recent years, biomedical research and clinical approaches evolved towards the production, translation and implementation of new technologies and practices. Among them, deep sequencing technologies unveiled that the cellular transcriptional world is far more composite than was initially hypothesized, as the most of genomic sequences could be widely transcribed not only into a very heterogeneous spectrum of protein-coding RNAs, but also in probably wider group of molecular and what their targets are, in order to clarify the biochemical pathways altered by oxidative stress and to identify the ideal candidates for novel molecular markers and therapeutic strategies useful for retinitis pigmentosa.

Cell Culture
Human RPE-derived Cells (H-RPE-Human Retinal Pigment Epithelial Cells, Clonetics™, Lonza, Walkersville, MD, USA) were cultivated in T-75 flasks containing RtEGM™ Retinal Pigment Epithelial Cell Growth Medium BulletKit ® (Clonetics™, Lonza, Walkersville, USA) with 2% v/v fetal bovine serum (FBS), 1% of penicillin/streptomycin and incubated at 37 • C with 5% CO 2 . H-RPE cells were, then, aliquoted into 96-well plates (4 × 10 4 cells/well) and grown for 24 h until the confluence prior to the treatment. Afterwards, A2E was added to a final concentration of 20 µM for 3 h and 6 h before rinsing with medium. Control cell groups were incubated without A2E. Reached the confluence, the cultures were shifted to phosphate-buffered saline with calcium, magnesium and glucose (PBS-CMG) and, soon after, exposed for 30 min to blue light produced by a tungsten halogen source (470 ± 20 nm; 0.4 mW/mm 2 ) to induce A2E phototoxicity before being incubated at 37 • C.

MTT Assay
Mitochondrial-dependent reduction of methylthiazolyldiphenyl-tetrazolium bromide (MTT) (Sigma-Aldrich, St. Louis, MO, USA) to formazan insoluble crystals was performed to evaluate RPE cell viability during treatment. Briefly, 10 µL of 5 mg/mL of MTT in PBS were added to the cell cultures following the A2E treatment. Following incubation at 37 • C for 2 h, a volume of 100 µL of 10% SDS in 0.01 mol/L HCl was added with the purpose to dissolve the crystals, before being incubated for 16 h. Afterwards, a Dynatech microplate reader permitted to determine the absorbance at 570 nm. Finally, we obtained the percentage of viable cells normalized with control conditions in the absence of A2E. The number of independent experiments was 3.

Total RNA Sequencing
Total RNA was extracted by TRIzol TM Reagent (Invitrogen TM , ThermoFisher Scientific, Waltham, MA, USA), following manufacturer's protocol and quantified at Qubit 2.0 fluorimeter by Qubit ® RNA assay kit (Invitrogen TM , ThermoFisher Scientific, Waltham, MA, USA). The RNA-seq samples consisted of 3 factor groups, represented by Human RPE cells before the treatment with A2E and at 2 following different time points of 3 h and 6 h, respectively. For each group 3 biological replicates were considered, for a total of 9 samples. 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), according to manufacturer's protocols. Sequencing runs were performed on an HiSeq 2500 Sequencer (Illumina), using the HiSeq SBS Kit v4 (Illumina).

Gene Expression Quantification and Normalization
Mapped reads were quantified by alignment-dependent expectation-maximization (EM) algorithm [28], particularly useful when most reads map equally well to multiple genes or transcripts. Once the algorithm has converged, every non-uniquely mapping read was assigned randomly to a particular transcript according to the abundances of the transcripts within the same mapping. The transcript per million reads (TPM) values were, then, computed from the counts assigned to each transcript, after normalization by the trimmed mean of M-values (TMM) method [29].

Long Non-Coding RNAs Alignment-Free Algorithms of Analysis
Previously described approaches are alignment-based, either applying multiple alignments to compute the conservation score or pairwise homology search for similar proteins. However, in order to avoid the limitations of alignment-based methods, the most reliable alignment-free algorithms have been applied. Among them, Coding-Potential Assessment Tool (CPAT) [43] and lncScore [44] discriminated non-coding transcripts from protein-coding by logistic regression model as machine learning classifier adopting sequence-based features such as open-reading frame (ORF) length, ORF coverage, ORF size, GC content, hexamer and Fickett scores. CNCI catalogued lncRNAs analyzing adjoining nucleotide triplets (ANT) to detect most-like CDS (MLCDS) regions in each transcript [45]. PLEK uses SVM classifier from LIBSVM package to calibrate k-mer frequencies of a sequence and sliding-window approach as features for classification [46]. Finally, FEELnc [47] accurately annotated lncRNAs using a Random Forest model trained with general features such as multi k-mer frequencies (from k = 1 to 12) and relaxed open reading frames.

Specific Circular RNAs Detection Pipelines
To detect CircRNAs in RNA-seq data an analysis of sequence reads spanning the back-splice junctions produced during circRNAs biogenesis is needed. Back-splice reads align to the genome in chiastic order, so circRNA detection from RNA-seq reads needs specific methods for non-collinear read mapping and analysis. To obtain a reliable result on circRNAs detection, several different methods were applied. The first, a "pseudo-reference-based" tool called KNIFE [48], directly constructed all the putative out-of-order exon-exon junction sequences from gene annotation data before alignment. Next, CIRCexplorer [49] and UROBORUS [50] followed the "fragmented-based" approach, which identified backsplicing junctions from the mapping information of a multiple-split read's alignment to the genome. Finally, CIRI [51] exploited its own algorithm based on paired chiastic clipping (PCC) signals detection, with an important reduction of potential false positives. Instructions provided in each tool manual were followed, filtering circRNAs with ≥2 back-spliced junction reads.

Differential lncRNAs Expression and Statistical Analysis
Differential expression analysis of lncRNAs was realized using Limma R package [52]. General linear models were applied to compare lncRNA expression changes at the different conditions of experimental design, setting the contrast groups as 0 h.untreated versus 3 h.treated, 0 h.untreated versus 6 h.treated, 3 h.treated versus 6 h.treated, [(0 h.untreated + 3 h.treated)/2] -[(0 h.untreated + 6 h.treated)/2]. The latter came from multiple group mean comparison, that allowed to indirectly capture the differences in expression level determined by the whole period of treatment, hereafter called "Due to Time" effect. For differentially expressed lncRNAs, the log 2 fold change (log 2 FC) of lncRNA abundance was calculated based on contrast groups and significance of expression changes was determined using the t-test [53]. p-values of multiple testing were, then, adjusted with Benjamini & Hochberg (BH) to correct false discovery rate (FDR) [54]. The lncRNAs uniquely identified in the RPE cells, showing 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) expression changes based on expression values ratio, and with BH-adjusted p-values lower than 0.05, were selected for functional classification of differentially expressed lncRNAs.

lncRNAs Validation by qRT-PCR
We selected the ten most dysregulated lncRNAs, obtained from RNA-seq data, to be validated by quantitative Real-Time polymerase chain reaction (qRT-PCR). Reverse transcription was performed following the manufacturer's protocol of GoScript™ Reverse Transcription System (Promega, Madison, WI, USA). The retrotranscribed cDNA was subjected to the qRT-PCR in the ABI 7500 fast sequence detection system (Applied Biosystems, Foster City, CA, USA), using the BRYT-Green based PCR reaction. PCR amplification was performed in a total reaction mixture of 20 µL, consisting of 20 ng cDNA, 10 µL 2 × GoTaq1qPCR Master Mix (Promega, USA) and 0.2 µM of each primer. The reaction was completed with the standard thermal cycle conditions applying the two-step qRT-PCR method: an initial denaturation at 95 • C for 30 s, followed by 40 cycles of 30 s at 95 • C and 30 s at 60 • C. Each reaction was replicated six times, considering all analyzed time points (18 samples), and the average threshold cycle (Ct) was calculated for each replicate. The lncRNAs expression was normalized to the expression level of most stable reference lncRNA, after being identified as combination of Delta Ct [55], GeNorm [56], NormFinder [56] and BestKeeper [57] algorithms results. The relative gene expression was, then, calculated using the 2 -∆∆Ct method, and the results were shown as the mean ± SEM (Standard Error of Mean). Statistical significance was determined by analysis of variance between groups (ANOVA), followed by Bonferroni post-hoc test. Finally, a linear regression analysis was performed to check the correlation of the FC of the gene expression ratios between qRT-PCR and RNA-Seq. The statistical analyses were all performed using IBM SPSS 26.0 software (IBM Corp, Armonk, NY, USA) [58]. Adjusted p-values < 0.05 were considered as statistically significant. The research was approved by the Scientific Ethics Committee of the Azienda Ospedaliera Universitaria-Policlinico "G. Martino" Messina.

lncRNA Host and Target Genes Pathway Analysis
GO term enrichment analysis for the most altered lncRNA host and target genes was performed using the ClueGO (v. 2.5.6) (INSERM, Paris, France) [59] and CluePedia (v. 1.5.6) (INSERM, Paris, France) [60] plugins in Cytoscape (v. 3.7.2) (National Institute of General Medical Sciences, Bethesda, MD, USA) [61]. Default parameters were used, except for Min#/% Genes = 40 and K-score threshold = 0.7, due to the huge number of nodes and edges deriving from analysis. Finally, only GO terms with p < 0.01 were selected.

Pathway Analysis of microRNA Targeting to Most Altered RPE Expressed lncRNAs
To evaluate the crosstalk between lncRNAs and miRNAs, we exploited the computational resource DIANA TOOLS mirPath v.3.0 (University of Thessaly, Thessaly, Greece) [62]. This bioinformatic web-server platform permitted us to obtain accurate statistics for possible pathways involving miRNAs targeting most dysregulated lncRNAs.

MTT Cell Viability Assay Showed an Exposure Time-Related Increased Death
The MTT cell viability assay highlighted a significant and different trend in RPE treated cells versus control. The addition of A2E to cultures led to a decreased percentage of cell viability, with the lowest peak at 6 h after treatment ( Figure 1).

Lnc-RNAs Validation by qRT-PCR
In order to validate the reliability of the RNA-Seq outputs, 10 among the most dysregulated lncRNAs (AC105052.4, LINC00968, AL645940.1, PSMG3-AS1, RDH10-AS1, SAMD12-AS1, GABPB1-AS1, NEAT1, AC068580.3, AC013451.2) were analyzed by qRT-PCR, and the obtained expression profiles resulted very similar to the transcriptome analysis profile (Figure 3 and Table S4). Primer characteristics are listed in Table 1. Results outputted from 2 -ΔΔCt method foresaw the normalization based on control group and the best stable lncRNAs GABPB1-AS1, arisen from geometric mean of individual rankings of the four applied algorithms ( Figure S1). The analysis of variance (ANOVA) method, conducted to compare the means among multiple groups, highlighted high significance, also resisted to conservative Bonferroni correction (p-values < 0.01). Linear regression analysis showed a

Lnc-RNAs Validation by qRT-PCR
In order to validate the reliability of the RNA-Seq outputs, 10 among the most dysregulated lncRNAs (AC105052.4, LINC00968, AL645940.1, PSMG3-AS1, RDH10-AS1, SAMD12-AS1, GABPB1-AS1, NEAT1, AC068580.3, AC013451.2) were analyzed by qRT-PCR, and the obtained expression profiles resulted very similar to the transcriptome analysis profile (Figure 3 and Table S4). Primer characteristics are listed in Table 1. Results outputted from 2 -∆∆Ct method foresaw the normalization based on control group and the best stable lncRNAs GABPB1-AS1, arisen from geometric mean of individual rankings of the four applied algorithms ( Figure S1). The analysis of variance (ANOVA) method, conducted to compare the means among multiple groups, highlighted high significance, also resisted to conservative Bonferroni correction (p-values < 0.01). Linear regression analysis showed a significantly positive correlation between gene expression ratios of qRT-PCR and RNA-Seq for each evaluated time point (R 2 = 0.99, r = 0.98), confirming our transcriptomic data validity.

Pathway Analysis of Selected lncRNAs Target and Host Genes Shed Light on Metabolic Pathways Impaired by Induced Oxidative Stress
Most dysregulated lncRNAs were firstly analyzed by RNA Interactome Database (Table S5), and found interactors (including query lncRNAs) were, then, pathway enriched by Cytoscape and its plugins ClueGO and CluePedia. Significant associations emerged between lncRNAs and pathways related to cell cycle alterations, along with cellular response to chemical stress and oxygen levels. However, the most interesting results came from the strong link shown by differentially expressed lncRNAs and metabolism impairments. As evidenced by GO Biological Process and Reactome databases in particular, a huge number of interactors showed high significance (corrected p-values near zero) in relationship to their involvement in nucleotides and nucleic acids metabolism ("Metabolism of nucleotides", 41 interactors; "DNA metabolic process", 380 interactors; "Positive regulation of DNA metabolic process", 109 interactors; "Metabolism of RNA", 276 interactors; "Regulation of mRNA metabolic process", 190 interactors; "Negative regulation of mRNA metabolic process", 44 interactors). Additionally, other metabolism-related pathways evidenced strong associations (corrected p-values near zero) with carbohydrate metabolism ("Glucose metabolism", 40 interactors; "C-type lectin receptor signaling pathway", 44 interactors; "Insulin signaling pathway", 55 interactors; "Insulin resistance", 47 interactors) and lipid metabolism ("Adipogenesis", 69 interactors, "Regulation of lipid metabolism by PPARalpha", 53 interactors, "PPARA activates gene expression", 52 interactors). Furthermore, a novel interesting output of pathway enrichment analysis is related to metabolism of cellular amide ("Regulation of cellular amide metabolic process", p-value = 0.000, 241 interactors; "Positive regulation of cellular amide metabolic process", p-value = 0.000, 92 interactors; "Negative regulation of cellular amide metabolic process", p-value = 0.000, 87 interactors). Finally, two individual regulative pathways, "ABC-family proteins mediated transport" (p-value = 0.000, 44 interactors) and "TP53 Regulates Metabolic Genes" (p-value = 0.000, 41 interactors), could be also involved in metabolism alterations induced by oxidative stress. Detailed results of most altered pathways and sub-pathways are available in Figures 4 and 5 and Table S6. The ClueGO plugin of Cytoscape permitted to find a first rich cluster of overrepresented GO processes and a network of connected GO terms was elaborated. Each node represents a GO biological process, and the colors refer to the GO group. Sixteen GO groups are present in the network, one representing GO biological process per group is named in the figure. The edges reflect the relationships between the terms based on the similarity of their associated genes.

Figure 5.
ClueGO enrichment analysis of a second cluster of overrepresented GO terms. The ClueGO plugin of Cytoscape permitted to find a second cluster of overrepresented GO processes and a network of connected GO terms was elaborated. Each node represents a GO biological process, and the colors refer to the GO group. Sixty-four GO groups are present in the network, one representing GO biological process per group is named in the figure. The edges reflect the relationships between the terms based on the similarity of their associated genes.

Discussion
The relevance of lncRNA genes has been established by their nearness to developmental regulators in the genome, enrichment of tissue-specific and developmental stage-specific expression patterns, and recurrent association with genetic traits [64]. LncRNAs regulate a huge number of molecular and cellular processes, such as epigenetic modifications, RNA splicing, mRNA decay, mRNA translation and molecular scaffold for structural/functional complexes [65]. Together with microRNAs and other small non-coding RNAs, expression alteration of lncRNAs could underlie pathogenesis of a broad range of human diseases [66]. Aside from their well-established roles in cancer [67], lncRNAs are crucial to cardiovascular [68], neurological [69], respiratory [70] and neurodegenerative [71] diseases. Among them, retinitis pigmentosa, an eye-related group of pathologies characterized by very heterogeneous genotypes but quite overlapping phenotypes, shows unusually complex molecular genetic causes, most of which are still unknown [72]. Recently, many experiments showed the important involvement of lncRNAs in retinal disorders, especially regulating angiogenesis, photoreceptor maturation, cell cycle in photoreceptor progenitor cells, apoptosis and cell viability, retinal vessel dysfunction, endothelial cell proliferation, vulnerability of the optic nerve, Muller cell proliferation and neurodegeneration [73]. Although the non-coding scenario of retinal dystrophies is now clearer than before, many obscure sides remain, especially about retinitis pigmentosa specific pathways and the connection between them.
Using recent methods of deep sequencing, we carried out an interesting analysis of the whole transcriptome of RPE cells during a follow-up of two time points (3 h and 6 h) after treatment with A2E. Oxidative stress is considered one of the most significant biochemical pathways involved in RP etiopathogenesis, especially targeting the high metabolic demand of RPE cells. Metabolic dysfunctions (e.g., ceramide synthesis [74]) and impairment of high energy processes as life-long light illumination or physiological phagocytosis could determine pathological modifications in blood-retinal barrier (BRB) and in other extracellular matrix molecules, leading to photoreceptor outer segments (POSs) processing inhibition and RPE cells apoptosis [75,76].
Thus, a deep knowledge of metabolism regulation of RPE cells, along with the identification of links between impairments in metabolic demand and downstream pathways, could shed new light on etiopathogenesis mechanisms of RP.
The realized lncRNA differential expression analysis highlighted a significant impairment in several pathways related to RPE metabolism, starting from nucleic acids. It is well known that lncRNAs are able to be involved in DNA metabolic processes, exerting epigenetic functions by acting as scaffolds for chromatin-modifying complexes [77]. The global up-regulation of DNA metabolic process (involving more than 500 elements between analyzed lncRNAs and their potential interactors) evidenced through observational time-points could determine an increased rDNA silencing, due to chromatin-associated lncRNAs tethered to specific genome sites by forming RNA-DNA hybrids acting in cis, via ATP-remodeling complex NoRC. Such effects could alter cellular growth and synthesis of ribosomes, leading to RPE cell death [78]. Interestingly, two of dysregulated lncRNAs, NEAT1 and MALAT1, are known to be involved in the constitution of paraspeckles nuclear bodies, involved in the trapping of adenosine-to-inosine edited RNA and in retention of serine/arginine (SR) splicing factors [79,80]. MALAT1, whose expression resulted mainly down-regulated after induced oxidative stress, could play its activity on cell survival, especially Muller and ganglion cells (RGCs). Reduced MALAT1 expression may impair neurotrophic factor-and stress-related gene expression by PI3K/Akt, MAPK, Ca2+/CaMK and cAMP/PKA pathways, all converging on the CREB family of leucine-zipper transcription factors [81]. NEAT1, slightly up-regulated during treatment period, could protect retinal Muller cells from apoptosis, also regulating splicing events involving genes related to cell cycle and cell death [82]. The opposite expression trends by MALAT1 and NEAT1 might represent a critical metabolic scenario with a final attempt by impaired cells to survive. Due to the wide spectrum of functions realized by lncRNAs, it is very probable that many of them could arise upon DNA damage, as evidenced by many of detected dysregulated lncRNAs involved in in double-strand break repair, DNA ionizing radiation (IR)-damage and cellular response via ATR, DNA damage and integrity checkpoints, signal transduction in response to DNA damage, also regulated by p53 class mediators. Two of the most interesting lncRNAs probably induced by DNA damage are MNX1-AS1 and MIR31HG, both over-expressed throughout observed time points. These lncRNAs interact with Cyclin D1 (CCND1) mRNA, whose encoding gene is already known to transcribe a specific lncRNA in DNA damage condition, and that exists on chromatin both as RNA:DNA hybrids and ssRNA, acting as transcription repressor [83]. Furthermore, various lncRNAs could exert their role outside the nucleus, hybridizing with the 3 untranslated regions (3 UTRs) of mRNAs to regulate their stability in the cytoplasm and/or interacting with miRNAs and RNA decay factors [84]. We found at least 40 lncRNAs, the most of which are antisense lncRNAs, acting as miRNA sponges, and they were mostly down-regulated. However, 26 antisense lncRNAs showed an heterogenous trend: 8 resulted over-expressed and 18 down-expressed. More interestingly, among them, four are the lncRNAs that highlighted the highest fold-change during the observed time points. The two most up-regulated (AC118658.1 and AC145207.2) sponge miRNAs are involved in transcriptional regulation, RNA transport and FoxO/AMPK signaling pathways, probably trying to positive regulate these processes, against the possible transcriptional and mitotic misregulation, along with extracellular matrix alterations, determined by miRNAs interacting with the other two most down-regulated (AL110504.1 and AC078909.1).
As expected, the induced oxidative stress determined changes also in carbohydrate metabolism, particularly in glucose metabolism. About 40 interactors between lncRNAs and their targets/host genes resulted implicated in bioenergetic reactions related to glucose high demand in RPE. In particular, BDNF-AS and TUG1, down-expressed and over-expressed respectively, are already known to affects cell viability by regulating glycolysis [85]. A dysregulation of these two lncRNAs is linked to high-glucose induced (DGI) apoptosis both in RPE and in neuron cells, like other retinal cytotypes [86,87]. Thus, the reduced expression of BDNF-AS and the increased expression of TUG1 could be interpreted as an attempt by RPE cells to protect themselves from apoptosis. Moreover, TUG1 is also an interactor of ARF1 and AKT1, both involved in insulin signaling pathway [88]. Very curiously, other two deregulated lncRNAs, CRNDE (down-expressed) and CYTOR (over-expressed), interact with the already cited ARF1 and AKT1, and with CAP1 and ACACA, also involved in insulin-related pathways [89,90]. Recently, it was experimentally proved that insulin signaling in the RPE may provide a paracrine signal to the retina for maintenance of photoreceptor function and survival, even if partially able to induce oxidative stress [91]. Additionally, high glucose level influences the synthesis of IGF-1, PEDF, advanced glycosylation end (AGE) products and their receptors (RAGE). In particular, it was proved that AGE-RAGE system, activated by glucose-stimulated RPE cells, could determine oxidative stress and inflammatory reactions, ultimately inducing retinal visual cycle impairments [92]. Thus, the combination of both CRNDE and CYTOR dysregulation could enforce the possible protection of RPE cells from DGI by the reduction of glucose intracellular intake and metabolism, as already discussed in relationship with colorectal cancer. Moreover, eight miRNAs, interacting with 10 found dysregulated lncRNAs, could be involved in lacto-and neolacto-glycosphingolipid biosynthesis [93]. Such analyzed lncRNAs showed a global up-regulation, probably silencing sponged miRNAs related to ceramide-ER Stress-AMPK Signaling axis, with the final consequence to induce RPE cell apoptosis.
One of the most investigated metabolic side of the retina regards lipids, mostly phospholipids and vitamin A derivatives. RPE expresses factors required for lipoprotein intake, cholesterol and lipoprotein synthesis and secretion, as well as reverse cholesterol transport (RCT) [94]. Phagocytosis and degradation of the POSs are the principal lipid source for the RPE [95], as well as processing and recycling back to the photoreceptors of vitamin A and of DHA [96,97]. Moreover, excess cholesterol is loaded on HDLs by ABC transporters and, then, removed, or it is also secreted in VLDL-like lipoproteins basolaterally into Bruch's membrane where it forms drusen [98]. Since HDLs cycle cholesterol between the RPE and photoreceptors, HDLs might accumulate the cholesterol coming out from rod outer segments in the subretinal space, as result of the compromised lipid transport following RPE impairment [99]. In order to perform these tasks, RPE cells present fine-regulated intracellular signaling pathways linking phagocytosis to the transcription of genes involved in lipid metabolism and homeostasis. Such pathways are mainly transcriptionally regulated by the peroxisome proliferator-activated receptor (PPAR) family of nuclear hormone receptors, which act as ligand-dependent transcription factors. This family also includes receptors for retinoids, vitamin D, and thyroid and steroid hormones [100]. The alfa subtype, in particular, stimulates fatty acid catabolism in several tissues known for their high rates of fatty acid oxidation (e.g., liver and eye) [101]. Our experiment highlighted a prevalent down-regulation of identified lncRNAs possibly involved in fatty acids biosynthesis and metabolism. In particular, the down-expression of AC007283.1 (interacting with UPF1 and ELAVL1) and AC012442.2 (interacting with HNRNPC and CSTF2T), along with the over-expression AC089983.1 (the antisense of TXNRD1, interacting with CRNDE) could alter the gene expression and lipid metabolism regulation by PPAR-alfa. Additionally, several down-regulated lncRNAs, such as the antisense AC004943.2 and the sense-overlapping AC007036.3, showed the possibility to interact with various miRNAs involved in fatty acids metabolism and biosynthesis, as well as in thyroid hormone synthesis, depicting a probable silencing scenario which could impair the integrity and functionality of lipidic retinal structures. Very interesting, 44 interactors between identified dysregulated lncRNAs and their targets/host genes resulted involved in ABC-family proteins mediated transport, probably impair the excess cholesterol removal, leading to cell death, especially following the 6h-treatment time point.
Finally, the protein metabolism also showed the involvement of numerous identified lncRNAs. It is well known how misfolding of a huge number of proteins related to retinal survival and vision process, like rhodopsin, can result in disruptions of cellular protein homeostasis [102]. Moreover, it is high probable that the interaction between several chaperones could be altered by oxidant agents, due to activation/deactivation switches that induce certain chaperone-encoding genes or co-chaperones, and repress other ones [103]. We detected two clusters made of dysregulated lncRNAs and their interactors/host genes, involved in positive (92 interactors) and in negative regulation of cellular amide metabolism (87 interactors), respectively. Among them, particularly interesting resulted the up-regulated PTEN-induced putative kinase protein 1 (PINK1) antisense RNA, and the down-regulated FMRP Translational Regulator 1 (FMR1)-IT1 sense intronic and vimentin (VIM) antisense 1 RNA. The first is well known to regulate mitochondrial damage, promote mitophagy and protect cells from death and apoptosis, especially during high glucose mediated regulation of RPE [90]. Thus, the over-expression of PINK1-AS could represent an oxidative stress induced signal reflecting the apoptotic status of treated RPE cells. In the meantime, the down-expression of both FMR1-IT1 and VIM-AS1, related to synaptogenesis, intracellular trafficking and cellular stability [104,105], could be seen as a final attempt of RPE cells to boost the production and sorting of vital proteins towards the most essential cellular districts. Additionally, nine miRNAs interacting with several down-regulated lncRNAs resulted involved in lysine degradation pathway. Several studies have documented N ε -Carboxy methyl lysine (N ε -CML), one of the most prevalent AGE, as a biomarker of diabetic retinopathy, linked to disruption of retinal photoreceptor external limiting membrane and ellipsoid zone [106]. Therefore, the reduced sponge activity by the previously cited lncRNAs could block lysine degradation, probably trying to arrest the production of AGEs, towards the final purpose to avoid the total cell death.

Conclusions
Nucleic acids are vital for cell survival, and hence for life. RNA plays innumerable roles in the cell-from gene expression modulation to regulation of protein translation or to control the architecture of whole chromosomes. LncRNAs represent a particular class of molecules characterized by the lack of protein-coding potential but with fundamental features including unique regulatory mechanisms, cis-regulatory activities, alternative forms of biogenesis and functional structured RNA domains. With modern high-throughput RNA-sequencing and advanced epigenomic technologies, the discovery rate of new lncRNA genes is rapidly increasing. We realized a whole transcriptome analysis to evaluate the lncRNA expression changes following the treatment with A2E in RPE cells. About 600 lncRNAs showed significant expression alterations in treated samples, involving target/host genes related to biochemical pathways associated to all major fields of cellular metabolism. Many specific sub-pathways were never linked to mechanisms that underlie retinal degeneration and related pathologies like retinitis pigmentosa. Despite this, a deeper transcriptome sequencing on a wider number of samples could permit to increase the number of detected lncRNAs, further clarifying regulative aspects of these non-coding RNAs in RP etiopathogenesis. Additionally, performing functional assays, (e.g., RNA pull-down, electrophoretic mobility shift assay (EMSA), RNA structure mapping, crosslinking immunoprecipitation (CLIP), phylogenetic lineage tracing, ribosome profiling, etc.) will permit to experimentally confirm the interaction between detected lncRNAs and their interactors. Exploiting innovative techniques, we will surely discover even more intriguing and exclusive features and functions of lncRNAs, improving their link to retinal disease etiopathogenesis.
Supplementary Materials: The following are available online at http://www.mdpi.com/2076-3921/9/4/318/s1, Figure S1: qRT-PCR selected lncRNAs stability results. The funnel chart shows the ranking of analyzed lncRNAs stability, obtained from the geometric mean of rankings values computed by Delta CT, GeNorm, NormFinder and BestKeeper algorithms. The lower stability value corresponds to more stably expressed gene., Table S1: RNA-Seq analysis on the gene level, Table S2: RNA-Seq differential expression analysis of the most dysregulated lncRNAs, Table S3: Annotation enrichment of most dysregulated lncRNAs, Table S4: Fold change expression of selected lncRNAs after A2E treatment, calculated by ∆∆Ct method, Table S5: RNAInteractor Database analysis of lncRNAs, Table S6: ClueGo pathway analysis of most dysregulated lncRNAs, Table S7: mirPath analysis of miRNAs sponged by most dysregulated lncRNAs.