The validation of Short Interspersed Nuclear Elements (SINEs) as a RT-qPCR normalization strategy in a rodent model for temporal lobe epilepsy

Background In gene expression studies via RT-qPCR many conclusions are inferred by using reference genes. However, it is generally known that also reference genes could be differentially expressed between various tissue types, experimental conditions and animal models. An increasing amount of studies have been performed to validate the stability of reference genes. In this study, two rodent-specific Short Interspersed Nuclear Elements (SINEs), which are located throughout the transcriptome, were validated and assessed against nine reference genes in a model of Temporal Lobe Epilepsy (TLE). Two different brain regions (i.e. hippocampus and cortex) and two different disease stages (i.e. acute phase and chronic phase) of the systemic kainic acid rat model for TLE were analyzed by performing expression analyses with the geNorm and NormFinder algorithms. Finally, we performed a rank aggregation analysis and validated the reference genes and the rodent-specific SINEs (i.e. B elements) individually via Gfap gene expression. Results GeNorm ranked Hprt1, Pgk1 and Ywhaz as the most stable genes in the acute phase, while Gusb and B2m were ranked as the most unstable, being significantly upregulated. The two B elements were ranked as most stable for both brain regions in the chronic phase by geNorm. In contrast, NormFinder ranked the B1 element only once as second best in cortical tissue for the chronic phase. Interestingly, using only one of the two algorithms would have led to skewed conclusions. Finally, the rank aggregation method indicated the use of the B1 element as the best option to normalize target genes, independent of the disease progression and brain region. This result was supported by the expression profile of Gfap. Conclusion In this study, we demonstrate the potential of implementing SINEs -notably the B1 element- as a stable normalization factor in a rodent model of TLE, independent of brain region or disease progression.


Introduction
The reverse transcription quantitative polymerase chain reaction (RT-qPCR) is considered as combined with a unique 120 bp sequence. A neuronally expressed BC1 gene is believed to be the origin of the ID elements [24]. In contrast, the B4 element resembles a fusion of a B1 and an ID element. The B elements, as well as the Alu repeats, are dispersed throughout the transcriptome, with the highest density in the intronic regions and at 16 kb upstream to the start position of specific protein-coding transcripts. In contrast, they are under-represented in exonic regions [25][26][27].
Previous RT-qPCR studies tried to obtain a better understanding of the molecular mechanisms of epilepsy, without or with a limited reference gene validation study [28][29][30][31]. In the present study, we have evaluated the stability of nine well-accepted reference genes and two SINEs (B1 and B2 elements) in the systemic kainic acid (KA) post status epilepticus (SE) rat model of Temporal Lobe Epilepsy (TLE). The expression stability of the transcripts was determined in the hippocampus and the cortex during the acute phase (day 10 post-SE) and chronic phase (day 80 and 120 post-SE) of disease progression (i.e. epileptogenesis). The RT-qPCR results were analyzed with geNorm [32], NormFinder [33] and rank aggregation [34] algorithms. As a validation strategy, normalization using the geNorm and NormFinder results was compared to normalization using the B elements or every single reference gene in the context of evaluation of the relative expression of glial fibrillary acidic protein (Gfap) mRNA. To our knowledge, this is the first large-scale RT-qPCR study which validated the robustness of B elements in a rodent model.

Animals
Male Sprague-Dawley rats (Harlan Laboratories B.V., the Netherlands) weighing 242.6 ± 15.1 g (9 weeks) were treated according to the guidelines of the European Communities Council Directive (2010/63/EU). The Animal Experimental Ethical Committee of Ghent University Hospital (ECD 16/06) approved the study protocol. The animals were conventionally housed in a temperature-controlled (20-23˚C) and humidity-controlled (50%) environment under a 12h/12h light/dark cycle, where food and water intake was ad libitum.
The rats were randomly distributed for KA administration (n = 19) and controls (n = 12). In order to induce SE, rats received 2-6 times KA (5 mg/kg; Tocris Bioscience, Bristol, UK) by intraperitoneal (i.p.) injections according to the protocol of Hellier et al., 1998. Seizure activity of all rats was continuously monitored visually and electrographically. The KA treatment was repeated hourly until the animals displayed a stable self-sustained SE for �3 h (i.e., >10 behavioral seizures per hour). Animals that exhibited excessive motor or excessive lethargic behavior were no longer injected with KA to avoid mortality [35,36]. The rats treated with KA were sampled after 10 (T10, n = 6), 80 (T80, n = 7) or 120 days (T120, n = 6). The control rats were treated with vehicle (saline, i.p.) and were sampled after 10 (acute phase control, n = 6) or 120 days (chronic phase control, n = 6). The animals were anesthetized via 5% isoflurane (Affygility Solutions, Broomfield, CO, USA) followed by decapitation and hippocampi and cortices were dissected on an ice-cold plate, the left and right parts were stored separately at -80˚C until RNA isolation.

RNA extraction and cDNA synthesis
Total RNA was isolated from left hippocampi and cortices with the RNeasy Plus Universal Mini Kit (Qiagen, Hilden, Germany), according to the manufacturer's instructions. The absorbance at 230 nm, 260 nm, 280 nm and total RNA concentration were measured with the Nano-Drop 1000 Spectrophotometer (ThermoFisher Scientific, Wilmington, DE, USA). The 260/ 280 nm ratio was used as RNA purity measure and all samples had values over 2.0. The integrity of the RNA was assessed by the analysis of the ratio of 28S to 18S rRNAs after agarose gel electrophoresis. Total RNA samples (15 μg) were treated with the Heat&Run gDNA removal kit (ArticZymes, Tromsø, Norway) to avoid amplification of genomic DNA during the reverse transcription.
The concentrations and purities were verified again, as described above. Only 1.5 μg of the total RNA was converted to cDNA using Oligo(dT)18 primers (Cell Signaling, Danvers, MA, USA) and M-MLV Reverse Transcriptase (Promega, Madison, WI, USA), according to the manufacturer's protocol. The absence of DNA contamination was controlled for each individual sample by omitting the reverse transcriptase (RT) from the procedure (i.e. minus RT samples). The minus RT and the plus RT samples of the control and the KA-treated samples of a similar disease progression stage (i.e. acute or chronic phase) were reverse transcribed simultaneously. After 1:20 dilution of the samples, a 10 μl fraction of every cDNA was pooled together in order to determine the PCR efficiency for every primer set. All samples were stored at -80˚C until further analysis.

Primer design and real-time PCR
Nine commonly used reference genes were chosen based on their differences in cellular or signaling pathways ( Table 1). The primers for the nine reference genes and Gfap were found in literature [9,14]. The primers for the B1 and B2 elements were designed based on the available consensus sequences [20,37]. The optimal primers were chosen from the regions with the highest consensus sequences and specificity was controlled in silico through the UCSC genome browser (Santa Cruz, CA, USA). All primers were synthesized by Sigma Aldrich (St. Louis, MO, USA).
Real-time PCRs were performed in white 96-well plates using the LightCycler 480 Instrument (Roche Diagnostics, Mannheim, Germany). In a reaction volume of 15 μl, 5 μl of the synthesized cDNA sample was pipetted. Thereby, the synthesized cDNA was diluted 1:40 for the reference genes and Gfap and 1:8000 for the B elements, to 3.75 ng and 0.019 ng of total cDNA respectively. The primer concentrations were 250 nM. Finally, 7.5 μl iQ SYBR Green Supermix (Bio-Rad Laboratories, Hercules, CA, USA) and 1.75 μl DEPC H 2 O was added to the reaction volume. The two-step amplification protocol of the real-time PCRs was as follows: initial denaturation at 95˚C for 5 minutes and 40 cycles of 95˚C for 15 seconds and 60˚C for 45 seconds. A light signal was acquired at the end of each cycle at 60˚C. The specificity of product formation was confirmed by a melting curve analysis (55˚C to 95˚C in increments of 0.11˚C/s). The Cq values were determined in technical triplicates with the LightCycler software v1.1.5 (Roche Applied Science, Mannheim, Germany) according to the second derivative method.
The PCR efficiencies were determined by technical duplications of a three-or five-point serial dilution of pooled cDNA, B elements or reference genes respectively. All the PCR efficiencies of the primer sets were between 90 and 110% (Table 2), demonstrating the robustness and reproducibility of the performed RT-qPCR assay.

Gene expression stability analysis and reference gene validation
The stability of the reference genes and B elements was assessed by geNorm (https://genorm. cmgg.be) and NormFinder (https://moma.dk/normfinder-software) [32,33]. In order to perform the geNorm analysis, the Cq values were imported into the qBaseplus software version 3.0 (Biogazelle, Ghent, Belgium) [38]. This algorithm calculates the gene stability measure M, which is the average pairwise variation of a particular gene with all other genes. The genes are ranked according to the calculated M value from the least stable (highest M value) to the most stable (lowest M value). The NormFinder algorithm estimates the stability through a model-based approach by comparing the overall gene expression variation and the gene variation between sample subgroups. This results in a stability value and ranks the genes based on their stable expression. The Cq values used in both algorithms were efficiency adjusted.
All the ranked gene lists obtained by the geNorm and NormFinder analyses were used in the rank aggregation R package called "RankAggreg" [34]. The weighted rank aggregation was applied to determine the most stable reference genes across the different disease stages, brain regions and algorithms. Thereby, we used the Cross-Entropy (CE) Monte Carlo algorithm, which generates random lists. The starting list was converged towards the best optimal list through an iteration procedure with the weighted Spearman's foot-rule distance, resulting in a consensus list of ranks. According to literature, patients and different animal models of TLE show significant upregulation in Gfap mRNA expressions [10, 14,[39][40][41]. For that reason, Gfap mRNA expression was chosen to evaluate the impact of the different kind of normalization strategies. The Gfap expression was evaluated with the most stable combination of references genes according to the two algorithms or to the eleven references independently, using the 2 ΔΔCT method [38].
The Cq values of each individual reference gene and the Gfap expression values, which were determined by different normalization methods, in the hippocampus and cortex were imported into GraphPad Prism version 6.00 for Windows (GraphPad Software, La Jolla, CA, USA). The various groups were statistically analyzed with the Mann-Whitney two-tailed U test (i.e. acute phase) or the nonparametric ANOVA by ranks of Kruskal-Wallis test followed by the Dunn's multiple comparisons post-hoc test (i.e. chronic phase). Differences were considered statistically significant at P � 0.05.

Transcription profiles
The mean Cq values represent a different intragroup variation profile of the reference genes in the hippocampus compared to the cortex (Fig 1). The genes B2m and Gusb had a significantly higher expression in the TLE model than in the controls at the acute phase (acute phase control vs. T10), with the exception of the Gusb cortical expression. The chronic phase shows a more disperse pattern, where the mRNA expression of Hprt1, Pgk1, Ywhaz, Actb and Rpl13a differs significantly between the control and TLE rats in the hippocampus, while Gusb, B2m, Rpl13a, Tbp, Ywhaz and the B2 element differ significantly between groups in the cortex. The significant differences in mean Cq values may indicate the instability of a reference gene or B element. The Cq value profile displayed a wide range of mRNA expression levels, where Gusb has high Cq values (30.42), thus a low mRNA expression, and Gapdh (18.37), Actb (19.13) and B2 element (19.41) were highly expressed.

Stability determination
The expression stability of the reference genes and the B elements was assessed by two algorithms. The geNorm algorithm labels the most stable genes with the lowest M value. This value is defined as the mean standard deviation of the logarithmically transformed expression values of the compared genes, which is calculated by the average variation among pairs of genes through the comparison of a control gene with other genes. First, the algorithm selects a pair of genes, which have the lowest M value, the additional genes are ranked based on the highest degree of compatibility with the other and with a geometric mean of the first pair [38]. The order of stability of the reference genes and B elements varied between the acute and chronic phase in the hippocampus (Fig 2A and 2C).
Nowadays, to perform correct normalization of a target gene (e.g. Gfap), multiple reference genes are recommended [4]. The qBaseplus software package also performs a pairwise variation analysis. This analysis calculates the variation between successive pairs, which is expressed as the normalization factor (NF). Hereby, the influence of a subsequent reference gene added to the previous pair is determined. The comparison is made with NF n for n number of genes to NF n+1 , which contain the same set of genes with an additional gene having a consecutive higher M value. If the added gene to the NF has a significant impact, the variation V n/n+1 between two NF factors is high. Supported by the data presented in the original publication, a cutoff value of 0.15 for V n/n+1 was determined [32]. When the value is below the cutoff, the addition of an extra reference gene is not necessary for the normalization of a target gene. Using two reference genes in the acute (Pgk1/Hprt1) or chronic (B1 and B2 element) model is sufficient, with values of 0.060 and 0.069 respectively (Fig 2B and 2D). The validation of SINEs for RT-qPCR normalization in a rodent model of TLE An equivalent geNorm analysis was performed in the cortex of the rat model for TLE (Fig  3). The algorithm determined that the Hprt1 and Ywhaz reference genes were sufficient for the normalization (V 2/3 = 0.032) in the acute phase. Similar as in the hippocampus, the B elements were also the most stably expressed in the chronic phase. The V value (0.123) is below the defined cutoff value (0.15), so according to the geNorm algorithm, the use of only the two B elements as reference genes is adequate for normalization.
The NormFinder algorithm uses a model-based estimation of variance approach to determine a stability value. This ANOVA model-based approach takes into account the average influence of a gene within the group and the individual impact of a group compared to the other groups, called the intra-and intergroup variations, respectively. The final value depends The validation of SINEs for RT-qPCR normalization in a rodent model of TLE on the number of candidate genes and samples analyzed [33]. The obtained results of the NormFinder analysis of this study are shown in Table 3. In the acute phase, the two most stable genes in the hippocampus are Rpl13a (0.164) and Actb (0.185), while Gusb (0.611) and B2m (0.693) are the two most unstable genes. In the chronic phase the three most stable genes were Rpl13a (0.129), Gapdh (0.143) and the B1 element (0.159). The most unstable genes were B2m (0.314) and Tbp (0.345), which differed slightly from the geNorm analysis. In the acute phase, the most stable reference genes in the cortex are Actb (0.112) and Gapdh (0.145). As before, the two most unstable genes were Gusb (0.356) and B2m (0.508). In the chronic phase, B2m (0.196) and B1 element (0.240) have the least diversity, while Rpl13a (0.359) and Ywhaz (0.437) turn out to be unreliable as reference genes. The calculated intra-and intergroup variations per gene are presented in S1 Table.

Rank aggregation
Not surprisingly, the different algorithms yielded a different ranking of the genes and B elements. In order to obtain a final consensus list of ranks, weighted rank aggregation was conducted. The R script for the weighted rank aggregation, via a CE Monte Carlo algorithm with the weighted Spearman's foot-rule distance, is presented in S1 File. The weight of the M values of the geNorm algorithm and the stability values of the NormFinder algorithm were taken into account in this iterative process [34]. As presented in Fig 4, the B1 element is the most stably expressed, independent of brain region and disease stage.

Validation of reference genes and B elements
Finally, the B elements and reference genes used in this study were validated to normalize the expression analysis of Gfap. From literature, it is known that expression of this gene is  The validation of SINEs for RT-qPCR normalization in a rodent model of TLE upregulated in several epilepsy models as well as in human patients with TLE [10, 14,[39][40][41]. The best pairs of genes resulting from the geNorm and NormFinder analysis and all the individually expressed genes and B elements were used to evaluate the pattern of Gfap in the hippocampus and the cortex, as depicted in Fig 5A and 5B, respectively. Normalized Gfap expression was very similar for both B elements as reference genes, in the hippocampus, as well as in the cortex. Moreover, normalization using the B elements was almost equivalent to normalization using the geNorm and NormFinder method. Only the chronic phase control versus the T120 in the NormFinder differed slightly compared to the B1 or B2 element normalization in the hippocampus. Even the use of unstably expressed genes as Gusb and B2m still allowed detection of the significant upregulation of Gfap mRNA expression in the two brain regions, in both the acute and chronic phase. Nevertheless, the upregulation of Gfap normalized by these reference genes differed in significance level compared to the normalization of Gfap with geNorm and NormFinder. Moreover, these reference genes were significantly upregulated in the acute phase of the TLE model, independent of brain region (S1 Fig). Remarkably, in the chronic phase of the TLE model normalization using Ywhaz didn't yield significant Gfap upregulation in the cortex.

Discussion
Currently, RT-qPCR is the 'gold standard' to determine variation in mRNA expression levels in biological materials between experimental conditions. Hereby, an accurate normalization strategy is key to infer correct conclusions from generated data. Predominantly, the expression of one or more endogenous 'reference' genes is used to correct for technical variations in the RT-qPCR protocol. Using multiple reference genes is recommended and the optimal number applied should be experimentally determined [4]. Importantly, some experimental procedures also influence reference gene expression, rendering specific genes unsuitable for normalization of target gene expression. Applying only a single reference gene could lead to erroneous conclusions in that case [32].
Our study is the first to validate the use of SINEs as a normalization strategy, comparing these to nine widely used and accepted reference genes as a normalization method in a rodent model. Previously, SINEs were validated in zebrafish tissues, human blood samples and human cell lines [15][16][17][18]. In order to perform the comparison, we selected nine reference genes from different functional classes, for example genes involved in major histocompatibility complex (B2m), transcription (Tbp) or metabolism (Gapdh), to avoid co-regulation ( Table 1). The use of ribosomal RNA (rRNA) as a normalizer is discouraged, because the rRNA (e.g. 18S rRNA) is not polyadenylated, transcribed by a different RNA polymerase and has a different function in the cell than mRNA. Although around 90% of the total RNA concentration contains rRNA, the mRNA/rRNA ratio can fluctuate depending on the experimental conditions [42,43]. Therefore, the rRNA was not included in our research. In contrast, the SINEs (i.e. B elements in the Rattus norvegicus) are part of the mRNA and are dispersed abundantly throughout the transcriptome.
The implementation as a normalization factor in RT-qPCR of another class of retrotransposons, the Long Interspersed Nuclear Elements (LINE-1), is discouraged as well. A difference in transcription of these autonomous retrotransposons cannot only be induced by various chemical and biological stressors, but also has been observed in cancers, neurodegenerative disorders and autoimmune diseases [44][45][46][47][48][49]. This would be due to hypomethylation of CpG islands and chromatin rearrangements [50][51][52]. In addition, a qPCR analysis of human brain tissue showed a high variation in LINE-1 copy number expression between individuals, especially in the hippocampus, and differential expression between neurons [53,54].
The majority of the B elements are either located in upstream regions of specific transcripts, or in their intronic regions [25][26][27]. Hence, PCR amplification of a specific repeat from a B element will multiply the same transcript of various genes. Predictably, the detection of many transcripts simultaneously will be less influenced by an unstably expressed gene. As a consequence, it can be expected that the detected signal will be less prone to variation during various biological conditions, than when a single gene is detected in RT-qPCR.
As a proof-of-concept we used a rodent model of TLE and sampled two different brain regions at various disease stages (i.e. acute phase and chronic phase). KA will induce SE, followed by a latency period, and finally spontaneous recurrent seizures occur [55,56]. The glutamate analog, KA, is a neuroexcitotoxic agent that acts via kainate receptors. Different pathologic changes will be manifested in the rodent, such as behavioral changes (e.g. occurrence of wet-dog shakes), electrophysiological alterations, induction of oxidative stress, The validation of SINEs for RT-qPCR normalization in a rodent model of TLE astrogliosis and aberrant mossy fiber sprouting. Primarily, the CA3 region of the hippocampus has been shown as the most susceptible region. To a lesser extent the entorhinal-and piriform cortex are targeted by KA [57][58][59]. Several studies have proposed reference genes for intrahippocampal KA and pilocarpine (PILO) induced models of TLE and human epileptic brain tissue [10, 14,60,61]. Because the expression of a reference gene may be unstable in different (pathological) models, this study evaluated the suitability of nine specifically chosen reference genes, as well as B elements, to normalize gene expression in the rodent model of TLE via systemic KA administration [14]. Where most studies analyze only the hippocampus, we also included the cortex, which is known to be affected as well.
An extra added value of our analysis is the application of both the geNorm and NormFinder algorithm to determine the stability of the B elements, unlike previous studies, which validated SINEs and inferred their conclusions only from one algorithm (i.e. geNorm) [16][17][18]. As demonstrated by our data, there is a clear difference between the rankings generated by the two algorithms ( Table 3, Figs 2 and 3). Where a target gene (e.g. Gfap) uses reference genes (e.g. Gapdh) as a point of reference to infer differences in expression over various experimental conditions, a reference gene itself does not have a point of reference. In order to select the most stable reference gene(s) over the experimental conditions, the two algorithms employ different mathematical schemes [32,33]. To merge the results of both algorithms, we applied a rank aggregation strategy, resulting in one ranking scheme, independent of biological variables [34]. From this analysis, we could conclude that the B1 element, but not the B2 element, functions as the most stable reference in our hands. Others have also reported that not every repeat element is stably expressed under different conditions [17]. Upregulation of mouse B2 elements under heat shock conditions has been reported. Hereby, the B2 element binds to RNA polymerase II, acting as a transcriptional repressor of protein-coding genes to prevent translation, thus reducing the number of misfolded proteins under hyperthermia [62][63][64]. An opposite effect was seen in the hippocampus under stress conditions, wherein it was hypothesized that proteins should be translated in the hippocampus to retain memories of successful escapes and danger cues for future situations. The above demonstrates the 'modulability' of the B2 element, possibly explaining why in our model its expression was somewhat less stable than that of the B1 element [65].
Expression of the astrocyte-specific cytoskeleton protein Gfap was evaluated, relative to every single reference gene or B element independently, or relative to the results of geNorm or NormFinder. It has been well established for several models of TLE that Gfap expression is highly upregulated during astrogliosis in human and animal epileptogenic hippocampi [10, 14,[39][40][41]. The expression pattern, when normalizing with the B elements, was similar to that obtained when using geNorm or NormFinder, which is an indication of a specific upregulated pattern. Remarkably, B2m and Gusb were significantly upregulated in both brain regions of the acute phase. In a PILO model, when these two genes were used in normalization, a significant upregulation of Gfap expression was not detectable [14]. Clearly, the Gusb and B2m genes are unsuitable as references genes in the acute phase, since both were significantly upregulated (S1 Fig). The rank aggregation also indicated these two genes as most unstable (Fig 4). Notable is Gapdh, which comes out as the second-best normalizer, while in a KA intrahippocampal model this gene was found to be an inadequate reference [10]. Furthermore, related models mimicking the same disease could differ significantly in reference gene stabilities, since Gusb was an unstable gene in the systemic PILO model, but stable in the intrahippocampal PILO model [14].
In the end, we believe that the reference genes are less affected in the acute phase, except for B2m and Gusb which were significantly upregulated. The systemic KA administration does influence the global expression of some well-accepted reference genes (e.g. Yhwaz and Hprt1) and their stability decreases with progression of the disease, when aberrant mossy fiber sprouting is more prominent. The B1 element remains more or less stable in M values or stability values independent of brain region or disease stage. This study demonstrates the potential of implementing of the B1 element as a new reference for correct normalization of target genes in the systemic KA induced rat model of TLE. We are currently investigating whether its implementation as a normalization strategy can be expanded to other disease models or rodent species.

Conclusion
Accurate normalization is essential in RT-qPCR studies to infer correct conclusions from the generated data. In this study, we investigated the stability of nine well-accepted reference genes and the transcripts of two SINEs (i.e. B1 and B2 element). The stability of the B elements and the performance as a normalizer were validated in a TLE rodent model induced by the systemic administration of KA. The exact same expression data resulted in different ranks between the geNorm and NormFinder algorithms. It is advisable to implement more than one algorithm in reference studies, or skewed conclusions may be made. The weighted rank aggregation generated a consensus list from both algorithms, but, importantly, this list was independent of brain region (i.e. hippocampus or cortex) and disease stage. Overall, the expression of the B1 element appears to be most stable, while B2m and Gusb were rather unreliable, even significantly upregulated in the acute phase of the TLE model in both evaluated brain regions. The validation of the new normalization strategy using Gfap, showed that the B1 element was comparable with the best pair of reference genes generated by geNorm and NormFinder. Thus, the B1 element can be implemented for the normalization of genes in this rodent model. Its general application in rodents should be studied further and confirmed by other pathology models and tissues.