**E**ff**ects of A2E-Induced Oxidative Stress on Retinal Epithelial Cells: New Insights on Di**ff**erential Gene Response and Retinal Dystrophies**

**Luigi Donato 1,2,\* , Rosalia D'Angelo 1,2, Simona Alibrandi 1,3, Carmela Rinaldi <sup>1</sup> , Antonina Sidoti 1,2,\* and Concetta Scimone 1,2**


Received: 20 March 2020; Accepted: 8 April 2020; Published: 10 April 2020

**Abstract:** Oxidative stress represents one of the principal inductors of lifestyle-related and genetic diseases. Among them, inherited retinal dystrophies, such as age-related macular degeneration and retinitis pigmentosa, are well known to be susceptible to oxidative stress. To better understand how high reactive oxygen species levels may be involved in retinal dystrophies onset and progression, we performed a whole RNA-Seq experiment. It consisted of a comparison of transcriptomes' profiles among human retinal pigment epithelium cells exposed to the oxidant agent N-retinylidene-N-retinylethanolamine (A2E), considering two time points (3h and 6h) after the basal one. The treatment with A2E determined relevant differences in gene expression and splicing events, involving several new pathways probably related to retinal degeneration. We found 10 different clusters of pathways involving differentially expressed and differentially alternative spliced genes and highlighted the sub- pathways which could depict a more detailed scenario determined by the oxidative-stress-induced condition. In particular, regulation and/or alterations of angiogenesis, extracellular matrix integrity, isoprenoid-mediated reactions, physiological or pathological autophagy, cell-death induction and retinal cell rescue represented the most dysregulated pathways. Our results could represent an important step towards discovery of unclear molecular mechanisms linking oxidative stress and etiopathogenesis of retinal dystrophies.

**Keywords:** RNA-Seq; RPE; Retinitis pigmentosa; A2E

### **1. Introduction**

Oxidative stress, recently defined as "a state where oxidative forces exceed the antioxidant systems due to loss of the balance between them", represents one of the principal inductors of lifestyle-related and genetic diseases [1]. Among them, retinal dystrophies and, in particular, age-related macular degeneration (AMD) and the subgroup of retinitis pigmentosa (RP) are well known to be susceptible to oxidative stress [2,3]. Nowadays AMD is considered the principal cause of visual disability among patients over 50 years [4]. The typical cellular sign of early AMD is represented by drusen, characteristic macular pigmentary deposits, associated with intermediate vision loss [5]. A "dry" and a "wet" form of AMD are currently known, the first more diffused but the second

responsible for 90% of acute blindness due to AMD [6]. As previously mentioned, the most relevant risk factor associated with AMD etiopathogenesis is represented by high levels of oxidative stress damaging the macula, generally induced by production of advanced glycation end products (AGE) and exposition to environmental factors [7]. The majority of such effects are also exerted by dysregulation of vascular endothelial growth factor (VEGF), impairment proteins and organelles clearance and glial cell dysfunctions [8]. RP consists of a very heterogeneous group of inherited eye disorders characterized by progressive vision loss [9]. RP has an incidence of 1:4000 people worldwide and represents the most prevalent form of photoreceptor-related diseases [10]. Primary symptoms can already occur in childhood or adolescence and, generally, consist of night blindness and gradual reduction of the visual field due to progressive death of rods. Total blindness, instead, is a feature typical of the late stage of the disease, following degeneration of macula photoreceptors [11]. Rods represent about 95% of all photoreceptors, and oxidative metabolism of fatty acids is their main energy source [12]. Main causes of rod death are genetic mutations, and more than 80 RP-causative genes have been already identified (https://sph.uth.edu/RetNet/sum-dis.htm#B-diseases), even if a relevant number of them are still unknown [13]. Conversely, cone degeneration is usually a late event frequently resulting from cytotoxic effects of high oxygen levels in the retina after rod reduction. Thus, oxidative damage is considered the first cause of cone apoptosis and progressive vision loss [14]. Interestingly, AMD and RP can also arise due to mutations in genes expressed in other retinal cell types, such as *MERTK* [15], *RLBP1* [16] and *RPE65* [17] expressed in retinal pigment epithelium (RPE). Originally, only a trophic function was hypothesized for RPE cells. Nowadays, it is well known that RPE is a monolayer of neural-crista-derived pigmented epithelial cells interacting with Bruch's membrane and choriocapillaris on the basolateral side and with the outer segments of the photoreceptors on the apical one [17]. RPE plays many vital roles for photoreceptor cells and the most fascinating certainly is the protection from oxidative stress [18]. Recent studies confirmed a high level of reactive oxygen species (ROS) in RPE, and fatty acids are one of their molecular targets. If oxidized, they may impair transduction pathways and gene expression [19]. Although fatty acid oxidation was already confirmed to cause macular degeneration, oxidative stress mechanism in RP development requires further clarifications [20]. Therefore, to better understand how high ROS levels may lead to retinal dystrophies onset and progression, we performed a comparison of transcriptomes' profiles among human RPE cells exposed to the oxidant agent N-retinylidene-N-retinylethanolamine (A2E). A2E is a toxic bis-retinoid that derives from the condensation and the oxidation of the trans-retinal [21]. Throughout life, A2E and other complex lipids accumulate and form lipofuscin in the RPE, ultimately determining photoreceptor death [22]. Additionally, A2E photo-oxidation is able to generate singlet oxygen, a highly reactive molecule that contributes to the increase of level of toxic metabolites such as epoxides and endoperoxides [23]. Moreover, lipid peroxidation is also responsible for A2E cleavage that releases cytotoxic reactive aldehydes. These reactive species could affect the lipid membranes fluidity and can damage the DNA and the cellular proteins [24]. On these bases, it was hypothesized that lipofuscin, A2E and its oxidized metabolites could accumulate if the cellular antioxidant system is unable to fight the oxidative-damaged lipids in rods and cones. Consequently, accumulation and photo-oxidation of A2E could lead to several retinal dystrophies, like already established for age-related macular degeneration (AMD) [25]. Various studies based on ARPE-19 cell line present pharmacological solution to retinal dystrophies. The same cell line is also used to identify new candidate compounds able to protect RPE against A2E oxidation [26]. As highlighted in the manuscript, the treatment with such a compound determines relevant differences in gene expression and splicing events, involving several new pathways probably related to retinal degeneration.

#### **2. Materials and Methods**

#### *2.1. Cell Culture*

Human RPE-derived Cells (H-RPE—Human Retinal Pigment Epithelial Cells, Clonetics™, Lonza, Walkersville, USA) were grown in T-75 flasks containing RtEGM™ Retinal Pigment Epithelial Cell Growth Medium BulletKit® (Clonetics™, Lonza, Walkersville, MD, USA) with 2% *v*/*v* fetal bovine serum (FBS), 1% of penicillin/streptomycin and incubated at 37 ◦C with 5% CO2. H-RPE cells were then plated into 96-well plates (4 <sup>×</sup> <sup>10</sup><sup>4</sup> cells/well) and cultured for 24 h to reach confluence before treatment. Subsequently, A2E was added to a final concentration of 20 µM for 6 h before rinsing with medium. Control cell groups were incubated without A2E. Confluent cultures were transferred to PBS supplemented with calcium, magnesium and glucose (PBS–CMG) and then exposed to blue light emitted by a tungsten-halogen source (470 <sup>±</sup> 20 nm; 0.4 mW/mm<sup>2</sup> ) for 30 min to induce phototoxicity of A2E and incubated at 37 ◦C. The 1–3 generation of subcultured RPE cells were used in this study.

#### *2.2. MTT Assay*

Cell viability was determined by mitochondrial-dependent reduction of methylthiazolyldiphenyltetrazolium bromide (MTT) (Sigma-Aldrich, St. Louis, MO, USA) to formazan-insoluble crystals. Briefly, 10 µL of 5 mg/mL of MTT in PBS was added to the cells following the A2E treatment. After incubation at 37 ◦C for 2 h, 100 µL of 10% SDS in 0.01 mol/L HCl was added to dissolve the crystals and incubated for 16 h. The absorbance was measured in a Dynatech microplate reader at 570 nm. Results were expressed as percentage of viable cells normalized with control conditions in the absence of A2E.

#### *2.3. Total RNA Sequencing*

Total RNA was extracted by TRIzolTM Reagent (InvitrogenTM, ThermoFisher Scientific, Waltham, MA, USA), following manufacturer's protocol, and quantified at Qubit 2.0 fluorimeter by Qubit® RNA assay kit (InvitrogenTM, 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. The selection of 3 h and 6 h time points was based on previous experiments realized by our research group (unpublished data), confirmed by outcomes from MTT assay in this work. Such results highlighted 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), according to manufacturer's protocols. Sequencing runs were performed on an HiSeq 2500 Sequencer (Illumina, San Diego, CA, USA), using the HiSeq SBS Kit v4 (Illumina, San Diego, CA, USA). The experiment was repeated thrice.

#### *2.4. Quality Validation and Read Mapping*

Sequence reads were generated from RPE-specific cDNA libraries on the Illumina HiSeq 2500 Sequencer. Obtained raw sequences were filtered to remove low-quality reads (average per base Phred score <30) and adaptor sequences. The quality of analyzed data was checked using FastQC (v.0.11.9) (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and QualiMap (v.2.2.1) [27], while trimming was realized by Trimmomatic (v.0.39). Filtered data were then mapped by CLC Genomics Workbench v.20.0 (https://digitalinsights.qiagen.com/products-overview/analysis-andvisualization/qiagen-clc-genomics-workbench/) against the Homo sapiens genome hg38 and the RNA database v.91, on Ensembl database. RNA-seq analysis was conducted using the following settings: quality trim limit = 0.01, ambiguity trim maximum value = 2. Map to annotated reference was as follows: mismatch cost = 2, insertion and deletion costs = 3, minimum length fraction and minimum

similarity fraction = 0.8, maximum number of hits for a read = 10, strand-specific = both, expression value = TPM. Raw data are available upon request.

#### *2.5. Gene Expression and Statistical Analysis*

Mapped reads were quantified by alignment-dependent and alignment-independent methods. The first approach uses the expectation-maximization (EM) algorithm [28] in order to determine expressions even in cases where the majority of 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. The second method has a higher accuracy for the point expression estimation and also allows the user to bootstrap the expression quantification to get an estimate of the technical variability. This approach was applied by the Salmon tool [29] using the transcript fasta files downloaded from GENCODE v.32 (gencode.v32.transcripts.fa). Salmon was run with the following settings for the RNA-seq data: quant –index/index –libType U –unmatedReads /single.fastq –incompatPrior '0.0' –biasSpeedSamp '5 ' –fldMax '1000' –fldMean '250' –fldSD '25' –forgettingFactor '0.65' –maxReadOcc '100' –numBiasSamples '2000000' –numAuxModelSamples '5000000' –numPreAuxModelSamples '5000' –numGibbsSamples '0' –numBootstraps '0' –thinningFactor '16' –sigDigits '3' –vbPrior '1e-05' -o/output. Once obtained, Salmon outputs were imported using tximport R package version 1.10.0 and lengthScaledTPM method [30] to generate read counts and Transcripts Per Million (TPMs). Low expressed transcripts and genes were filtered based on the data mean–variance trend analysis. The expected decreasing trend between data mean and variance was observed when expressed transcripts were determined as which had ≥1 of the 9 samples with count per million reads (CPM) ≥1, which provided an optimal filter of low expression. A gene was expressed if any of its transcripts with the above criteria were expressed. The trimmed mean of M-values (TMM) method was used to normalize the gene and transcript read counts to *log*2-CPM [31]. The principal component analysis (PCA) plot showed that the RNA-seq data did not have distinct batch effects, so it was possible to proceed with downstream analyses.

#### *2.6. DE, DAS and DTU Analysis*

Limma R package was used for differential expression analyses [32]. General linear models were established to compare gene and transcript 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 versus (3 h.treated + 6 h.treated)/2. For differentially expressed (DE) genes/transcripts, the *log*<sup>2</sup> fold change (*L*2*FC*) of gene/transcript abundance were calculated based on contrast groups and significance of expression changes were determined using the t-test [33]. *P*-values of multiple testing were adjusted with BH to correct false discovery rate (FDR) [34]. A gene/transcript was significantly DE in a contrast group if it had adjusted *p*-value <0.01 and *L*2*FC* ≥1. At the alternative splicing level, differential transcript usage (DTU) transcripts were determined by comparing the *L*2*FC* of a transcript to the weighted average of *L*2*FCs* (weights were based on their standard deviation) of all remaining transcripts in the same gene. A transcript was determined as significant DTU if it had adjusted *p*-value <0.01 and ∆PS ≥0.1. For differentially alternative spliced (DAS) genes, *L*2*FC* of each individual transcript was compared to gene level *L*2*FC*, which was calculated as the weighted average of *L*2*FCs* of all transcripts of the gene. Then *p*-values of individual transcript comparisons were summarized to a single gene level *p*-value with F-test. A gene was significantly DAS in a contrast group if it had an adjusted *p*-value <0.01 and any of its transcript had a ∆ Percent Spliced (∆PS) ratio ≥0.1. Finally, time points (0 h, 3 h, 6 h) in groups (untreated, treated) were used for time-series trend analysis. Natural Cubic Spline method with degree of freedom was used to generate time-series trend (Figure S1).

#### *2.7. Gene-Enrichment and Functional Pathway Analysis*

The up- and down-regulated genes were analyzed by the Database for Annotation, Visualization and Integrated Discovery (DAVID) 6.8 [35]. This tool is based on more than 40 annotation categories, including GO terms to protein–protein interactions, from disease associations to gene functional summaries, and many others. In DAVID annotation system, EASE Score, a modified Fisher Exact P-Value, is adopted to measure the gene-enrichment in annotation terms. The EASE score provides a conservative adjustment to the Fisher exact probability that weights significance in favor of themes supported by more genes. The EASE score is calculated by penalizing (removing) one gene within the given category from the list and calculating the resulting Fisher exact probability for that category.

#### *2.8. Selection of Single-Pathway "Master genes" and Selection of Retinitis Pigmentosa Candidate Genes by ToppGene Prioritization*

In order to highlight new candidate genes involved into retinitis pigmentosa, based on oxidative-related candidate pathways, we chose the 15 most altered genes for each one. Firstly, we considered them for each time point; then, we chose the commons in all time points. Subsequently, chosen genes underwent prioritization by ToppGene (https://toppgene.cchmc.org), a web-based tool able to classify a selected group of candidate genes from a large set of genes correlated with a pathology, giving each one a score. The score is based on the intersection of data from various databases of annotations related to cellular and physiological functions, analyzing complex networks shared between genes already known to cause the disease (training genes) and candidate genes (test genes). Training genes were obtained from RetNet online database.

#### *2.9. Data Validation by qRT-PCR.*

Ten most dysregulated mRNAs from candidate genes previously identified were selected and validated by quantitative Real-Time-Polymerase Chain Reaction (qRT-PCR), in order to validate RNA-Seq data. Reverse transcription was carried out according to the manufacturer's protocol of GoScript™ Reverse Transcription System (Promega, Madison, WI, USA). The obtained cDNA was subjected to RT-PCR in the ABI 7500 fast sequence detection system (Applied Biosystems, Foster City, CA, USA), using BRYT-Green-based PCR reaction. PCR amplification was performed in a total reaction mixture of 20 µL containing 20 ng cDNA, 10 µL 2 × GoTaq1qPCR Master Mix (Promega, Madison, WI, USA) and 0.2 µM of each primer. PCR was run with the standard thermal cycle conditions using 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 expression of mRNAs was calculated relative to expression level of endogenous control β-actin, and the relative expression of genes was calculated using the 2−∆∆Ct method [36]. 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 (https://www.ibm.com/analytics/us/en/technology/spss/). Adjusted *p*-values <0.05 were considered statistically significant. The research was approved by the Scientific Ethics Committee of the Azienda Ospedaliera Universitaria—Policlinico "G. Martino" Messina.

#### **3. Results**

#### *3.1. MTT Cell Viability Assay Results*

The MTT cell viability assay showed a relevant and different trend in RPE-treated cells versus control. The addition of A2E to cultures led to a dose-dependent increase in cell death percentage (Figure 1).

*Antioxidants* **2020**, *9*, x FOR PEER REVIEW 6 of 24

**Figure 1.** MTT determination of A2E treatment in retinal pigment epithelium (RPE) cells. Cell death was assessed at considered time points (3 h and 6 h) in A2E treated samples compared to basal untreated group. Results are shown as mean ± standard error of mean (*n* = 3). *p*-value < 0.05. **Figure 1.** MTT determination of A2E treatment in retinal pigment epithelium (RPE) cells. Cell death was assessed at considered time points (3 h and 6 h) in A2E treated samples compared to basal untreated group. Results are shown as mean ± standard error of mean (*n* = 3). *p*-value < 0.05.

#### *3.2. Sequencing Analysis and Mapping Statistics 3.2. Sequencing Analysis and Mapping Statistics*

RNA sequencing carried out on Illumina HiSeq 2500 yielded a total of 96,346,180 quality reads (mean mapping quality = 29) with a percentage of 67.8% uniquely mapped. A total of 16,173 genes and 69,653 transcripts were identified out of 60,609 and 227,462 reference counterparts, respectively, considering the whole human transcriptome. The annotated reference assembly v.32 (GRCh38.p13) was downloaded from GeneCode FTP server (ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode\_human/). All previous mapping statistics were based on average values calculated for all three replicates in each time point. Details are available in RNA sequencing carried out on Illumina HiSeq 2500 yielded a total of 96,346,180 quality reads (mean mapping quality = 29) with a percentage of 67.8% uniquely mapped. A total of 16,173 genes and 69,653 transcripts were identified out of 60,609 and 227,462 reference counterparts, respectively, considering the whole human transcriptome. The annotated reference assembly v.32 (GRCh38.p13) was downloaded from GeneCode FTP server (ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode\_human/). All previous mapping statistics were based on average values calculated for all three replicates in each time point. Details are available in Figure S2.

#### Figure S2. *3.3. Analysis of Gene Expression Profile of RPE Cells*

*3.3. Analysis of Gene Expression Profile of RPE Cells* As previously cited, in our transcriptome study, 16,173 genes (Table S1) and 69,653 transcripts (Table S2) were expressed with log2-CPM values ≥ 1 considered as the average value per sample. We first analyzed differential expression at the gene level (DE) (Figure 2a and Table S3). Considering the stringent criteria we have chosen, we identified 2432 genes that were significantly differentially expressed in response to A2E treatment. Of these, 59.7% resulted as up-regulated, while 40.3% downregulated, showing a decreasing time-dependent trend for up-regulated ones and an increasing trend for down-regulated ones. Afterward, the analysis of transcript-level data allowed us to identify genes that were DAS between the contrast groups (Figure 2b and Table S4). We detected 5119 DAS genes, of which 1101 were also DE genes (regulated by both transcription and AS) and 4108 were regulated by AS only. Therefore, considering all 6540 genes that showed significantly altered levels of differential gene and/or transcript-level expression, 78.3% were differentially alternatively spliced, highlighting a consistent response to the induced oxidative stress. Furthermore, to identify the specific transcripts that characterize a gene as DAS, a DTU analysis was performed (Figure 2b and Table S5). Globally, around 12% (8587) of expressed transcripts were classified as DTU. The next step consisted in the analysis of early response to the induced oxidative stress in comparison with the As previously cited, in our transcriptome study, 16,173 genes (Table S1) and 69,653 transcripts (Table S2) were expressed with log2-CPM values ≥1 considered as the average value per sample. We first analyzed differential expression at the gene level (DE) (Figure 2a and Table S3). Considering the stringent criteria we have chosen, we identified 2432 genes that were significantly differentially expressed in response to A2E treatment. Of these, 59.7% resulted as up-regulated, while 40.3% down-regulated, showing a decreasing time-dependent trend for up-regulated ones and an increasing trend for down-regulated ones. Afterward, the analysis of transcript-level data allowed us to identify genes that were DAS between the contrast groups (Figure 2b and Table S4). We detected 5119 DAS genes, of which 1101 were also DE genes (regulated by both transcription and AS) and 4108 were regulated by AS only. Therefore, considering all 6540 genes that showed significantly altered levels of differential gene and/or transcript-level expression, 78.3% were differentially alternatively spliced, highlighting a consistent response to the induced oxidative stress. Furthermore, to identify the specific transcripts that characterize a gene as DAS, a DTU analysis was performed (Figure 2b and Table S5). Globally, around 12% (8587) of expressed transcripts were classified as DTU. The next step consisted in the analysis of early response to the induced oxidative stress in comparison with the response to late treatment with A2E. About 40% (858) and 55% (2457) of DE and DAS genes, respectively, resulted

as common to both time point observations, while the residual was unique to 3 h and 6 h (Figure 3). Consequently, it is evident how changes in gene-level expression and alternative splicing occurred throughout the whole period, either transiently (occurring after 3 h and returning to initial level after 6 h), occurring later (only after 6 h from treatment) or enduring throughout the whole period. Such results probably suggest the different responses of RPE cells as a growing resistance to oxidative stress stimuli in a time-dependent manner. *Antioxidants* **2020**, *9*, x FOR PEER REVIEW 7 of 24 respectively, resulted as common to both time point observations, while the residual was unique to 3 h and 6 h (Figure 3). Consequently, it is evident how changes in gene-level expression and alternative splicing occurred throughout the whole period, either transiently (occurring after 3 h and returning to initial level after 6 h), occurring later (only after 6 h from treatment) or enduring throughout the whole period. Such results probably suggest the different responses of RPE cells as a growing resistance to oxidative stress stimuli in a time-dependent manner.

**Figure 2.** Summary figure of expressed genes and significant DE, DE + DAS, DE + DTU, DAS and DTU genes from analysis of the considered time points of RPE cell data. (**A**) Number of genes regulated only by transcription (DE), only by alternative splicing (DAS) and by both transcription and alternative splicing (DE + DAS); (**B**) number of transcripts regulated only by transcription (DE), only by alternative splicing (DAS) and by both transcription and alternative splicing (DE + DAS). DTU = Differential Transcript Usage. AS = Alternative Splicing. **Figure 2.** Summary Figure of expressed genes and significant DE, DE + DAS, DE + DTU, DAS and DTU genes from analysis of the considered time points of RPE cell data. (**A**) Number of genes regulated only by transcription (DE), only by alternative splicing (DAS) and by both transcription and alternative splicing (DE + DAS); (**B**) number of transcripts regulated only by transcription (DE), only by alternative splicing (DAS) and by both transcription and alternative splicing (DE + DAS). DTU = Differential Transcript Usage. AS = Alternative Splicing.

*Antioxidants* **2020**, *9*, x FOR PEER REVIEW 8 of 24

**Figure 3***.* Comparison of the gene lists generated during differential expression analyses. Euler diagrams of DE (**A**) and DAS (**B**) genes identified during expression analyses in considered 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 versus (3 h.treated + 6 h.treated)/2. **Figure 3.** Comparison of the gene lists generated during differential expression analyses. Euler diagrams of DE (**A**) and DAS (**B**) genes identified during expression analyses in considered 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 versus (3 h.treated + 6 h.treated)/2.

#### *3.4. DE and DAS Genes Highlighted Different Functionality Patterns 3.4. DE and DAS Genes Highlighted Di*ff*erent Functionality Patterns*

Functional enrichment analyses of DE and DAS genes showed relevant differences, as already supposed by the overlap of only 1011 genes (Figure 2). The most significantly enriched terms for DE genes were linked to vascular events and very heterogeneous responses to oxidative stress. As evidenced by hierarchical clustering of total gene expression levels of DE genes, adaptive, transient and late expression profiles followed the A2E induced stress, and an analog response was highlighted by transcript expression profiles of individual DE genes (Figure 4). Functional annotation of individual DE genes clusters was associated with methylation of RNA (cluster 1), various stress response (cluster 2), microtubule assembly and activity (cluster 3), angiogenesis (cluster 6), lipid biosynthesis (cluster 7), focal adhesions (cluster 8) and respiratory electron transport (cluster 10). For DAS genes, the most enriched functional terms deal with alteration of cellular proliferation and cell death. Hierarchical clustering of the DTU transcripts and expression profiles of individual DAS genes also evidenced heterogeneous response profiles. Functional annotation of DTU clusters of genes revealed enrichment of terms related to misfolded and damaged protein removal (clusters 1 and 3), autophagy (cluster 2), lipid biosynthesis (cluster 6), regulation of DNA damage response (cluster 7), induction of cell death (cluster 8) and translation regulation (cluster 10) (Figure 5). Functional enrichment analyses of DE and DAS genes showed relevant differences, as already supposed by the overlap of only 1011 genes (Figure 2). The most significantly enriched terms for DE genes were linked to vascular events and very heterogeneous responses to oxidative stress. As evidenced by hierarchical clustering of total gene expression levels of DE genes, adaptive, transient and late expression profiles followed the A2E induced stress, and an analog response was highlighted by transcript expression profiles of individual DE genes (Figure 4). Functional annotation of individual DE genes clusters was associated with methylation of RNA (cluster 1), various stress response (cluster 2), microtubule assembly and activity (cluster 3), angiogenesis (cluster 6), lipid biosynthesis (cluster 7), focal adhesions (cluster 8) and respiratory electron transport (cluster 10). For DAS genes, the most enriched functional terms deal with alteration of cellular proliferation and cell death. Hierarchical clustering of the DTU transcripts and expression profiles of individual DAS genes also evidenced heterogeneous response profiles. Functional annotation of DTU clusters of genes revealed enrichment of terms related to misfolded and damaged protein removal (clusters 1 and 3), autophagy (cluster 2), lipid biosynthesis (cluster 6), regulation of DNA damage response (cluster 7), induction of cell death (cluster 8) and translation regulation (cluster 10) (Figure 5).

**Figure 4.** Hierarchical clustering and heatmap of DE genes and Key GO terms**.** DE genes show segregation into 10 coexpressed clusters, of which the main ones related to oxidative stress were highlighted (circled) and linked to GO specific terms. Full results of GO enrichment analyses are shown in Tables S6 and S7. The z-score scale represents mean-subtracted regularized log-transformed transcripts per million (TPMs)*.* **Figure 4.** Hierarchical clustering and heatmap of DE genes and Key GO terms. DE genes show segregation into 10 coexpressed clusters, of which the main ones related to oxidative stress were highlighted (circled) and linked to GO specific terms. Full results of GO enrichment analyses are shown in Tables S6 and S7. The z-score scale represents mean-subtracted regularized log-transformed transcripts per million (TPMs).

*Antioxidants* **2020**, *9*, x FOR PEER REVIEW 10 of 24

**Figure 5***.* Hierarchical clustering and heatmap of DE genes + DTU transcripts from DAS genes and Key GO terms. DE genes and DTU transcripts from DAS genes show segregation into 10 coexpressed clusters, of which the main ones related to oxidative stress were highlighted (circled) and linked to GO specific terms. Full results of GO enrichment analyses are shown in Tables S6 and S7. The z-score scale represents mean-subtracted regularized log-transformed TPMs*.* **Figure 5.** Hierarchical clustering and heatmap of DE genes + DTU transcripts from DAS genes and Key GO terms. DE genes and DTU transcripts from DAS genes show segregation into 10 coexpressed clusters, of which the main ones related to oxidative stress were highlighted (circled) and linked to GO specific terms. Full results of GO enrichment analyses are shown in Tables S6 and S7. The z-score scale represents mean-subtracted regularized log-transformed TPMs.

#### *3.5. Early Cellular Response to Induced Stress Mainly Involves Pre-mRNA Splicing and Glycolysis-Related DE and DAS Genes 3.5. Early Cellular Response to Induced Stress Mainly Involves Pre-mRNA Splicing and Glycolysis-Related DE and DAS Genes*

The applied statistical model used in our analyses permitted us to determine precisely at which time point each DE and DAS gene first showed a significant change, together with the magnitude and course of that change. The top 15 significant genes or transcripts belonging to the considered groups were deeply investigated, allowing the identification of new candidate genes and specific pathways possibly involved in retinal dystrophy etiopathogenesis (Figure 6). After the first 3 h of treatment, only one gene (*HNRNPA3P6*) reached a relevant DE status, and it is involved in cytoplasmic trafficking of RNA and pre-mRNA splicing. Main DAS genes that showed early important differences (*CTSH* and *GPI*), instead, were related to glycolysis, and the same biochemical pathway was the most correlated to initial DTU transcripts expression dysregulation (ENST00000615999.4, ENST00000588991.7, ENST00000586425.2, ENST00000525807.5 and ENST00000550050.5).

**Figure 6.** *Cont.*

**Figure 6.** Volcano plots of significant DE and DAS genes and of DE and DTU transcripts. Volcano plots of significant (adjusted *p*-value < 0.01) DE genes (**A**), DE transcripts (**B**), DAS genes (**C**) and DTU transcripts (**D**). The low expressed genes and transcripts were filtered. The top 15 considered elements with the smallest corrected *p*-values and bigger fold-changes are highlighted, and different colors refer to different contrast groups. DE genes: log2FC vs. −log10(FDR) at gene level; DAS genes: maximum ∆PS of transcript in a gene vs. −log10(FDR) at gene level; DE transcripts: log2FC vs. −log10(FDR) at gene level at transcript level and DTU transcripts: ∆PS vs. −log10(FDR) at transcript level.

#### *3.6. Late RPE Cell Response to A2E Treatment Could Impair Bioenergetic Specific Reactions, Extracellular Matrix Integrity and Neurotransmission-Related DE and DAS Genes*

While the number of significant and relevant DE genes was limited in the early stage of treatment, it grew over time, reaching a climax at 6 h (Figure 6). At this time point, many of the DE genes (*ACTG1*, *CCN2*, *RPL19*, *RPL3*, *P4HB*, *RPS11*, *FTL*, *CAPZB* and *RNA5-8SN2*) (Figures 7 and 8) resulted as linked to new particular pathways, as iron metabolism, plasma lipoproteins assembly and F-actin capping. Furthermore, several over-expressed DE genes at 6 h (*TTC8*, *ARL3*, *REEP6*, *GUCA1B* and *PDE6G*) resulted as associated with alternative splicing of retina-preferred gene transcripts. About DAS genes, instead, the most significant ones (*ACADVL*, *GPI* and *LTBP3*) resulted as correlated to already-involved pathways (e.g., glycolysis). However, more interestingly, the same genes were together with DTU transcripts (ENST00000586425.2 and ENST00000588991.7), but even more interestingly, with other biochemical activities regarding oxidative processes in mitochondria (e.g., fatty acids reactions) and low conductance of potassium channels (as evidenced by DE transcripts ENST00000577650.5, ENST00000451956.1, ENST00000556690.5 and ENST00000551173.5).

*Antioxidants* **2020**, *9*, x FOR PEER REVIEW 14 of 24

**Figure 7***.* Expression profiles of most important DE genes/transcripts; first cluster associated with newly discovered pathways linked to oxidative stress. Detailed gene/transcript expression profiles

and RPS11) linked to newly identified candidate pathways linked to oxidative stress.

**Figure 8.** Expression profiles of most important DE genes/transcripts secondo cluster associated with newly discovered pathways linked to oxidative stress. Detailed gene/transcript expression profiles across the time course of the second cluster of most important DE genes (*P4HB*, *RNA5-8SN2, RPL3* and *RPL19)* linked to newly identified candidate pathways linked to oxidative stress.
