miRNAexpression profile of retinal pigment epithelial cells under oxidative stress conditions

Deep analysis of regulative mechanisms of transcription and translation in eukaryotes could improve knowledge of many genetic pathologies such as retinitis pigmentosa (RP). New layers of complexity have recently emerged with the discovery that ‘junk’ DNA is transcribed and, among these, miRNAs have assumed a preponderant role. We compared changes in the expression of miRNAs obtained from whole transcriptome analyses, between two groups of retinal pigment epithelium (RPE) cells, one untreated and the other exposed to the oxidant agent oxidized low‐density lipoprotein (oxLDL), examining four time points (1, 2, 4 and 6 h). We found that 23 miRNAs exhibited altered expression in the treated samples, targeting genes involved in several biochemical pathways, many of them associated to RP for the first time, such as those mediated by insulin receptor signaling and son of sevenless. Moreover, five RP causative genes (KLHL7, RDH11,CERKL, AIPL1 and USH1G) emerged as already validated targets of five altered miRNAs (hsa‐miR‐1307, hsa‐miR‐3064, hsa‐miR‐4709, hsa‐miR‐3615 and hsa‐miR‐637), suggesting a tight connection between induced oxidative stress and RP development and progression. This miRNA expression analysis of oxidative stress‐induced RPE cells has discovered new regulative functions of miRNAs in RP that should lead to the discovery of new ways to regulate the etiopathogenesis of RP.

Deep analysis of regulative mechanisms of transcription and translation in eukaryotes could improve knowledge of many genetic pathologies such as retinitis pigmentosa (RP). New layers of complexity have recently emerged with the discovery that 'junk' DNA is transcribed and, among these, miR-NAs have assumed a preponderant role. We compared changes in the expression of miRNAs obtained from whole transcriptome analyses, between two groups of retinal pigment epithelium (RPE) cells, one untreated and the other exposed to the oxidant agent oxidized low-density lipoprotein (oxLDL), examining four time points (1, 2, 4 and 6 h). We found that 23 miRNAs exhibited altered expression in the treated samples, targeting genes involved in several biochemical pathways, many of them associated to RP for the first time, such as those mediated by insulin receptor signaling and son of sevenless. Moreover, five RP causative genes (KLHL7, RDH11, CERKL, AIPL1 and USH1G) emerged as already validated targets of five altered miRNAs (hsa-miR-1307, hsa-miR-3064, hsa-miR-4709, hsa-miR-3615 and hsa-miR-637), suggesting a tight connection between induced oxidative stress and RP development and progression. This miRNA expression analysis of oxidative stress-induced RPE cells has discovered new regulative functions of miRNAs in RP that should lead to the discovery of new ways to regulate the etiopathogenesis of RP.
Analysis of regulative mechanisms of transcription and translation in eukaryotes is a very challenging task. New layers of complexity are daily discovered, such as the preponderant role in regulative functions of miR-NAs. miRNAs represent a group of short non-coding RNAs that induce transcript degradation or translational inhibition of their target mRNAs, acting as post-transcriptional regulators of gene expression [1].
They assume the role of key regulators of several important biological processes, in both physiological and pathological conditions [2], controlling specific pathways by targeting networks of functionally correlated genes. Alterations of miRNA expression, due to mutations in either the miRNA itself or its target genes, could lead to several pathological conditions such as cancers [3], neurodegenerative and genetic diseases [4]. Therefore miRNAs, due to their emerging role as disease biomarkers, might be possible therapeutic targets in human disorders [5].
The retina is the back portion of the eye, photosensitive and able to focus light signals towards the optic nerve first, then towards the brain, after transduction of them into electrical stimuli. This light-sensitive layer of the eye represents the target of a huge number of human inherited pathologies, such as retinitis pigmentosa (RP) [6]. RP is a genetic disease that determines retinal degeneration by inducing a slow and progressive death in photoreceptors and retinal pigment epithelium (RPE) cells [7], leading to loss of ability to transmit to the brain the visual information. The term 'pigmentosa' deals with the characteristic appearance, during the advanced states of the disease, of abnormal areas of pigment in the retina. Much evidence supports the role of miRNAs in normal retinal development and functions [8]. Alterations of miRNA regulation in conditional Dicer mouse mutant eyes reduce and damage normal development of lens, cornea, retina and optic chiasm [9]. Furthermore, post-developmental disruption of miRNA processing in photoreceptors leads to severe functional impairments [10]. Interestingly, the targeted deletion of specific retina-enriched miRNAs has relevant effects on vertebrate eye [11], such as pathogenic roles in human retinal diseases [12]. Nowadays, miRNA expression data come only from the analysis of murine miRNA transcriptome (miRNome) but, due to structural and functional differences between human and mouse retinas, they are not totally useful. Therefore, an improved knowledge of human retina miRNome, especially of patients affected by retinal disease, could lead to a better understanding of the physiopathology of this tissue. In this work, we investigated the complexity of human retina miR-Nome, analyzing data from human RPE cell transcriptomes. RPE represents a single layer of post-mitotic cells, acting as both a selective barrier to and a vegetative regulator of the overlying photoreceptor layer, thereby playing a key role in its maintenance. Due to its specific proteins, RPE helps to renew outer segments by phagocytizing the spent discs of photoreceptor outer segments, regulates the trafficking of nutrients and waste products to and from the retina, protects the outer retina from excessive high-energy light and light-generated reactive oxygen species and maintains retinal homeostasis through the release of diffusible factors. In detail, we compared miRNA expression changes between two group of RPE cells, one exposed to the oxidant agent oxidized low-density lipoprotein (oxLDL) and the other untreated, considering four time points (1, 2, 4 and 6 h) over basal one (time zero). oxLDL was chosen because it has already been tested on several neurodegenerative diseases but, above all, because it was seen that high cholesterol level could be linked to RP development and progression [13]. The main purpose of our work was to discover which miRNAs changed during oxidative stress induction and what their targets are, in order to better understand how reactive oxygen species might lead to RP development.

Materials and methods
This study was approved by the Ethics Committee of Azienda Policlinico Universitario 'G. Martino' Messina.
Total RNA sequencing RNA was extracted by TRIzol TM reagent (Invitrogen, Thermo Fisher Scientific, Waltham, MA, USA), following manufacturer's protocol, and quantified with a Qubit 2.0 fluorimeter by QubitÒ RNA assay kit (Thermo Fisher Scientific). Expression analysis was realized comparing human RPE cells treated with 100 lgÁmL À1 of oxLDL and untreated ones, both at the moment of treatment and for four different time points (1, 2, 4 and 6 h). Libraries were generated using 1 lg of total RNA and the Ion Total RNA-Seq Kit v2 (Thermo Fisher Scientific), then purified by Dyn-abeadsÒ magnetic beads and quantified with a Qubit 2.0 fluorimeter with Qubit dsDNA HS Assay Kit. An appropriate library amount was used for clonal amplification performed with the Ion PI TM Template OT2 200 Kit v2 (Thermo Fisher Scientific) on the Ion One Touch TM 2 System; template-positive Ion Sphere Particles were enriched with the Ion One Touch TM Enrichment System. Sequencing runs were performed on an Ion Proton TM Sequencer (Ion Torrent technology; Thermo Fisher Scientific), using the Ion PI TM Sequencing 200 Kit v2 and the Ion PI TM Chip Kit v2 (Thermo Fisher Scientific). The experiment was repeated three times. were filtered to remove low quality reads (average per base Phred score < 28). Furthermore, the reads containing adaptor sequences and low-quality sequences (reads presenting ambiguous bases denoted as 'N') were also trimmed from the raw data. The quality of analyzed data was checked using FASTQC (v.0.11.5) and QUALIMAP (v.2.2.1) software. The filtered data were then analyzed by CLC GENOMICS WORK-BENCH v.10.0.1 (Qiagen Aarhus, Denmark; https://www.qia genbioinformatics.com/products/clc-genomics-workbench/) using Homo sapiens genome hg19 and Ensembl RNA database v.74 as references. RNA-Seq analysis was conducted using the following settings: quality trim limit = 0.01, ambiguity trim maximum value = 2. Map to annotated reference: minimum length fraction and minimum similarity fraction = 0.8, maximum number of hits/read = 2, type of organism = eukaryote, paired settings =default.

Small RNA analysis
The applied approach counted the different types of small RNAs in the data and compared them with databases of miRNAs or other small RNAs. Once whole RNA-Seq data were imported, the small RNAs were extracted and counted, in order to create a small RNA sample that could be used for further steps. Sequences were filtered based on length (reads below 15 bp and above 55 bp were discarded) and on minimum sampling count (set at 1). The aligned and selected reads were grouped on the sequence of the mature miRNAs allowing up to two mismatches within the exact length of the reference mature sequence (i.e. excluding trimming or extension variants). Subsequently, the number of reads mapping on each mature miRNA was counted and then normalized using either the trimmed mean of Mvalues (TMM) method [14] or reads per million (CPM). Finally, the small RNA sample produced when counting the tags was enriched by comparing the tag sequences with the annotation resources miRBase (v21) and Ensembl noncoding RNA database (v74).

Gene expression and statistical analysis
The original expression values were log2 transformed and normalized in order to ensure that samples are comparable and assumptions on the data for analysis are met [15]. In order to highlight the miRNAs with different level of expression between untreated and treated samples, and for the four considered time points, we categorized them into two groups, based on count ratios (fold-change): (a) upregulated (fold change > 1); (b) down-regulated (0 < fold change < 1). Furthermore, because the fold changes are linear, for any value smaller than 1 (i.e. for down-regulation), we chose to replace the value by its negative reciprocal value, in order to make the variation more noticeable (for instance, 2-fold downregulation is indicated by a value of -2 instead of 0.5). Due to the few biological replicates available for each of the experimental groups studied (only three replicates for each considered time point), but with many features to be studied simultaneously (miRNAs in a whole transcriptome), we applied the empirical analysis of DGE (EDGE) statistical algorithm, which implements the 'exact test' for two-group comparisons developed by Robinson and Smyth [16]. The test is based on the assumption that the count data follow a negative binomial distribution, which in contrast to the Poisson distribution has the characteristic that it allows for a non-constant meanvariance relationship. The 'exact test' of Robinson and Smyth is similar to Fisher's exact test, but also accounts for over-dispersion caused by biological variability. The miRNAs uniquely identified in the RPE cells with at least five unique gene reads, greater than one-fold (up-regulated) or lower than one-fold (down-regulated) changes in expression based on the ratio of expression values, and with Bonferroni-adjusted P values lower than 0.05 were selected for functional categorization of differentially expressed miR-NAs.

Validation of miRNAs by qRT-PCR
To confirm the transcriptome results, the 23 analyzed miRNA were validated by qRT-PCR. Complementary DNA synthesis from miRNA samples was performed using the miScript II RT Kit and HiSpec Buffer (Qiagen). The obtained cDNA was subjected to the RT-PCR in the ABI 7500 fast sequence detection system (Applied Biosystems, Thermo Fisher Scientific), using the BRYT-Green based PCR reaction. PCR amplification was performed in a total reaction mixture of 20 lL, containing 20 ng cDNA, 10 lL 29 GoTaq1qPCR Master Mix (Promega, Madison, WI, USA) and 0.2 lM of each primer. The 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 run in triplicate, considering all selected time points (1, 2, 4 and 6 h), and the average threshold cycle (C t ) was calculated for each replicate. The expression of miRNAs was calculated relative to the expression level of the endogenous control U6, and the relative expression of a gene was calculated using the 2 ÀDDCt method. The correlation of the fold-change of the gene expression ratios between qRT-PCR and RNA-Seq was checked by linear regression analysis in SPSS Statistics v24.0 software (IBM Corp., Armonk, NY, USA).

Sequencing analysis and mapping statistics
RNA sequencing carried out with Ion Torrent yielded a total of 11 214 300 quality reads (mean mapping quality = 32.92) with mean read length of 155.03 bp. All reads were previously aligned to GRCh37/hg19 reference assembly, and then to known human miRNAs (miRBase v21) [21] and GRCh37 non-coding RNAs [22]. About 71 500 small RNAs were founded in all samples, 69 158 of which were annotated and about 2341 unannotated. In details, 23 miRNAs, each one considered as a 'group' made of known precursors and/or mature miRNAs resulting from mapping (Table 1), showed expression alterations at the analyzed time points (Table S1). All previous mapping statistics are based on average values calculated for all three replicates in each time point. Detailed information on RNA-Seq statistics are available in Table 2.

Expression analysis
The predominant length of the resulting mature miRNAs was 21 bp. We identified a total of 55 mature miRNAs with an average expression level of three reads across all considered treated and untreated RPE cells cultures ( Table 1). The variability was low across samples, with an expected higher trend for miRNAs expressed at lower levels. As highlighted in Table 3 Table 1. Altered miRNAs with ranking. The RNA-Seq analysis highlighted the 23 grouped miRNAs (mature, precursors and precursor variants) with expression alterations, ranked on their abundance (based on read count). As seen in table, the 10 top-ranked miRNAs accounted for almost 80% of the total count, and the top five for 60%.

Rank
ID Sequence  (Fig. 1). We ranked the expressed miRNAs based on their abundance ( Table 1). The 10 top-ranked miRNAs accounted for almost 80% of the total count and the top five for 60%.

qRT-PCR verification
To validate the authenticity and reproducibility of the RNA-Seq results, the 23 analyzed miRNAs were selected for qRT-PCR analysis, and the obtained expression profiles were similar to the results of transcriptome analysis (data not shown). The linear regression analysis showed a significantly positive correlation of the relationship between gene expression ratios of qRT-PCR and RNA-Seq for all selected time points (Fig. 2), confirming our transcriptomic data validity.

miRNA target identification and pathway analysis
There were 1402 genes that were experimentally validated targets of 23 miRNAs found altered in miR TarBase.

Discussion
Retinitis pigmentosa is an ocular disease with very heterogeneous phenotypes with unusually complex molecular genetic causes [23]. Such an intricate picture, primarily determined by locus and allelic heterogeneity [24], is worsened by the actual lack of knowledge on all possible causative genes and their function. Moreover, little is known about involvement of regulative non-coding RNAs [25], which are already considered the most promising targets of experimental therapies [26]. Using high-throughput sequencing technologies, we analyzed the whole transcriptome of RPE cells treated with oxLDL during a follow-up of four time points (1, 2, 4 and 6 h) after exposure, and compared the results with untreated cells. Due to the high coverage of our sequencing experiment, we were able to detect miRNAs and sequence variants that had a low level of expression. Furthermore, the parallel analysis of three replicates for each selected group for each time point, along with a specific statistical test, allowed us to obtain reliable data, overcoming possible biasrelated variability in miRNA expression levels and nucleotide sequences. Oxidative stress plays a critical role in the etiopathogenesis of RP [27], leading to pathobiological changes including outer blood-retina barrier breakdown and senescence of RPE cells [28]. RPE cells are very susceptible to oxidative stress, due to high metabolic activity, including physiological phagocytosis and life-long light illumination [29]. Therefore, as high cholesterol could be linked to RP development and progression [30], oxLDL could represent the source of pathobiological changes of RPE cells [31], such as outer blood-retina barrier dysfunction [32], inhibition of processing of photoreceptor outer segments by RPE [33], expression of transforming growth factor-b 2 [34], synthesis alterations of extracellular matrix components [35], increasing RPE apoptosis [36] and senescence changes [31]. Moreover, several miRNAs were already associated to retinal development and function in vertebrates, such as miR-216a and miR-23a [37,38], although at very low levels.

Conclusions
We analyzed whole transcriptomes of two group of RPE cells, treated with oxLDL and untreated, comparing miRNA expression changes at four selected time points (1, 2, 4 and 6 h) from time zero. We found that 23 grouped miRNAs exhibited expression alterations in treated samples, targeting genes involved in several biochemical pathways. Most of these, such as fatty acid metabolism and the ubiquitin proteasome pathway, are already known to be directly involved in retinal degenerations. Several others, instead, might be associated for the first time to RP etiopathogenesis, such as IRS-and SOS-mediated pathways. Moreover, five RP causative genes (KLHL7, RDH11, CERKL, AIPL1 and USH1G) emerged as already validated targets of five altered miRNAs (hsa-miR-1307, hsa-miR-3064, hsa-miR-4709, hsa-miR-3615 and hsa-miR-637), suggesting a tight connection between induced oxidative stress and RP development and progression, thanks to the important junction ring represented by regulative functions of miRNAs. Nevertheless, many other important aspects have to be investigated, such as variant miRNA targets and isomiR, along with predicted target validation. Additionally, we have to underline that predicted miRNA targets resulted from in silico analyses and, even if they are based on statistical significance algorithms and literature data, they will need to be experimentally validated. Furthermore, a deeper transcriptome sequencing on a larger number of samples could permit us to increase the number of detected miRNAs, improving knowledge of regulative functions of these small RNAs and RP.

Author contributions
LD planned experiments and wrote the manuscript; PB contributed reagents; CS performed experiments; CR performed manuscript supervision; RD analyzed data; AS wrote and drafted the manuscript.

Supporting information
Additional Supporting Information may be found online in the supporting information tab for this article: Table S1. miRNAs precursors and precursors variants coming from RNA-Seq analysis. Table S2. mirPath KEGG and GO analysis. Table S3. ClueGO detailed pathway analysis of miR-NAs target genes from miRTarBase and microT databases.