# Rattus norvegicus reference genome evaluation for hippocampus RNA-seq data analysis

Christophe Le Priol, Andrée Delahaye-Duriez

<figure><img src="image_1.png" /></figure>

To cite this version:

Christophe Le Priol, Andrée Delahaye-Duriez. Rattus norvegicus reference genome evaluation for hippocampus RNA-seq data analysis. 2023. inserm-04713573

HAL Id: inserm-04713573

https://inserm.hal.science/inserm-04713573v1

Preprint submitted on 29 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Rattus norvegicus reference genome evaluation for hippocampus RNA-seq data analysis

Christophe Le Priol (✉ chris.lepriol@gmail.com)

Université Paris-Cité

Andrée Delahaye-Duriez

Université Paris-Cité

## Research Article

**Keywords:** RNA-seq, reference, genome, annotations, Rattus norvegicus

**Posted Date:** July 10th, 2023

**DOI:** https://doi.org/10.21203/rs.3.rs-3128726/v1

**License:** © This work is licensed under a Creative Commons Attribution 4.0 International License.

Read Full License

**Additional Declarations:** No competing interests reported.
*Rattus norvegicus* reference genome evaluation
for hippocampus RNA-seq data analysis

Christophe Le Priol¹* and Andrée Delahaye-Duriez¹,²,³

¹Inserm UMR 1141 NeuroDiderot, Université de Paris Cité, Paris,
75019, France.

²UFR Santé Médecine Biologie Humaine, Université Sorbonne Paris
Nord, Bobigny, 93000, France.

³Unité fonctionnelle de médecine génomique et génétique clinique,
Hôpital Jean Verdier, AP-HP, Bondy, 93140, France.

*Corresponding author(s). E-mail(s): chris.lepriol@gmail.com;

Contributing authors: andree.delahaye@inserm.fr;

**Abstract**

**Background:** One of the first steps of a usual RNA-seq data analysis work-
flow consists in quantifying gene expression by aligning the sequencing reads
to a reference genome and counting the aligned reads in its annotated regions.
Downstream analysis, such as the identification of differentially expressed genes,
strongly rely on the quality of this process. In addition to the performance of
alignment methods, the choice of reference genome, when several of them are
available, may strongly impact this step.

**Results:** Here, we propose to evaluate the effect of widely used *Rattus norvegicus* reference genomes consisting of Ensembl and RefSeq annotations of the Rnor_6.0 assembly and the ones based on the mRatBN7.2 genome assembly recently published by RefSeq, on a classical differential expression workflow. We re-analyzed published RNA-seq datasets from different hippocampal subregions and revealed that the Ensembl and Refseq reference genome based on the mRatBN7.2 assembly provide an improvement of read mapping statistics. We showed that the RefSeq annotations of this assembly make gene expression quantification and differentially expressed gene identification more reliable thanks to overall longer exon length in comparison with Ensembl annotations. Moreover, we identified specific biologically relevant results using RefSeq annotations of the new genome assembly.

**Conclusion:** Overall, the biological interpretation of the differential expression
analysis of the analyzed datasets may be dramatically impacted by the choice of
reference genome. Therefore, we believe that this choice should be more carefully addressed and that our approach could extend to other tissues and species.

**Keywords:** RNA-seq, reference, genome, annotations, *Rattus norvegicus*

# 1 Introduction

## 1.1 Evaluation of the choice of reference genome on RNA-seq studies

The sequencing of the transcriptome has enabled to assess genome-wide changes in gene expression in a variety of biological contexts. One of the first steps of a usual RNA-seq data analysis workflow is to quantify gene expression by counting the reads mapped to annotated genomic regions. To do so, a reference genome, which consist of both assembly sequences and annotations, is required. Subsequent common analysis, e.g. differential expression analysis or variant calling, depend on the quality of the gene expression quantification. Methods dedicated to each of these analysis have been extensively evaluated, e.g. read mapping [1] or differential expression analysis [2]. On the contrary, the choice of reference genome in genomic studies remains rarely discussed. While tools have been developed to compare genome annotations for many years [3, 4], the impact of the choice of a reference genome on RNA-seq data analysis remains largely underexplored. Yet, a variety of analysis may be impacted, such as differential expression analysis, alternative splicing analysis or variant calling [5, 6]. The variability of entire analysis pipelines, including the use of different reference genomes, were assessed for both bulk and single-cell data [7, 8] but only a few studies specifically focused on this issue. Zhao *et al* have demonstrated that the expression quantification and the differential expression assessment are indeed dramatically affected by the choice of reference human genome for large sets of genes [9]. Besides, they showed that this effect is tissue-dependent. Slabaugh *et al* conducted the same type of analysis with rice genomes, with an additional interest in splicing analysis, and came to the same conclusion [10]. It has been showed that conservative and extensively curated annotations are recommended for accurate gene expression quantification [11] while more complex annotations may negatively impact the identification differentially expressed genes [12] and should be used for exploratory research [13]. When several reference genomes are available for a species, the choice between them may be done arbitrarily or lack strong rationales. To overcome this situation, Chen *et al* quantified the effect of two reference genomes for Chinese hamster ovary cells on gene expression quantification using RNA-seq datasets and proposed a strategy for a reliable identification of differentially expressed genes based on two metrics they introduced: the consistency and the stringency [14].
## 1.2 *Rattus norvegicus* reference genomes

The Normay rat, *Rattus norvegicus* species, is a widely used experimental model in medical and biological research [15]. The Rnor_6.0 assembly, published in 2014, has been the gold standard reference genome sequence in RNA-seq data studies these last few years. Recently, a new genome assembly, mRatBN7.2, was published [16]. This new genome assembly exhibit overall better quality, taking benefits from years of significant improvements in sequencing since the release of the previous genome assembly [17]. In particular, this assembly is composed of fewer and longer scaffolds and its much higher coverage reduces the probability of sequencing errors, which is crucial for variant calling studies (table 1).

<table><caption>**Table 1** Sequence statistics of *Rattus norvegicus* genome assemblies</caption><thead><tr><th>Assembly 1</th><th>Rnor_6.0</th><th>mRatBN7.2</th></tr></thead><tbody><tr><td>Publication date</td><td>2014/07/01<br>rn6, RGSC 6.0</td><td>2020/11/10<br>rn7</td></tr><tr><td>Other names</td><td>RGSC Rnor_6.0<br>RGSC Genome Assembly v6.0</td><td></td></tr><tr><td>GenBank accession</td><td>GCA_000001895.4</td><td>GCA_015227675.2</td></tr><tr><td>RefSeq accession</td><td>GCF_000001895.5</td><td>GCF_015227675.2</td></tr><tr><td>Total sequence length</td><td>2 870 184 193</td><td>2 647 915 728</td></tr><tr><td>Number of scaffolds</td><td>440</td><td>176</td></tr><tr><td>N50</td><td>14 986 627</td><td>135 012 528</td></tr><tr><td>L50</td><td>65</td><td>8</td></tr><tr><td>Coverage</td><td>3x-10.0x</td><td>92.0x</td></tr></tbody></table>

A few months after the assembly, the NCBI published annotations in RefSeq release 108, followed by Ensembl in its release 105. Similarly to assembly metrics, these new annotations exhibit improvements in comparison with the most widely used annotations of the previous assembly in terms of number of genes or exons (table 2).

<table><caption>**Table 2** Annotation statistics of *Rattus norvegicus* reference genomes.</caption><thead><tr><th>Annotation</th><th>Ensembl 104</th><th>RefSeq 106</th><th>RefSeq 108</th><th>Ensembl 105</th></tr></thead><tbody><tr><td>Assembly</td><td>Rnor_6.0</td><td>Rnor_6.0</td><td>mRatBN7.2</td><td>mRatBN7.2</td></tr><tr><td>Publication date</td><td>2015/06</td><td>2016/07/27</td><td>2021/01/21</td><td>2021/12</td></tr><tr><td>Genes</td><td>32 883</td><td>39 537</td><td>42 172</td><td>30 560</td></tr><tr><td>Transcripts</td><td>32 945</td><td>56 098</td><td>34 325</td><td>30 565</td></tr><tr><td>Exons</td><td>228 617</td><td>257 370</td><td>263 738</td><td>243 601</td></tr><tr><td>Exon length</td><td>57 272 488</td><td>91 967 979</td><td>108 936 556</td><td>74 911 182</td></tr></tbody></table>

Number of genes and exons and cumulative exon length in the annotations of *Rattus norvegicus* genomes in Ensembl releases 104 and 105 and RefSeq releases 106 and 108. The GTF file of each reference genome was used to calculate these statistics (more details in Methods).
Overall, these better assembly quality and annotation metrics let hope significant refinements in RNA-seq data analysis.

The few studies evaluating the impact of the choice of reference genome on RNA-seq data analysis mostly deal with the human genome. No particular attention has been paid yet to the *Rattus norvegicus* species on this topic. Here, we propose to evaluate the effect of common genome annotations, *i.e.* Ensembl release 104 and RefSeq release 106, of the most widely used *Rattus norvegicus* genome sequence, Rnor_6.0 assembly, and the annotations of the new mRatBN7.2 genome assembly published in Ensembl release 105 and RefSeq release 108, on a classical differential expression workflow based on RNA-seq data. We re-analyzed RNA-seq datasets of published studies composed of samples from different hippocampal subregions, a tissue of interest for neurological diseases such as epilepsy or neurocognitive impairments.

# 2 Results

We evaluated the impact of the choice of reference genomes on gene expression quantification and differential expression analysis using published RNA-seq data from different hippocampal subregions. We made pairwise reference genome comparisons and only considered the pairs composed of two reference genomes from the same database, *i.e.* either Ensembl or RefSeq, and the one composed of the Ensembl and RefSeq genomes based on the mRatBN7.2 assembly.

## 2.1 Hippocampal subregion comparisons in normal condition

First, we analyzed two RNA-seq datasets of different hippocampal subregions in normal condition. Machado *et al* microdissected and sequenced the CA1, CA2, CA3 and DG subregions of non-treated rats in order to characterize the transcriptomic profiles of the different hippocampal subregions in normal condition (GEO Series GSE179101) [18]. Smith *et al* sequenced the hippocampal CA1, CA3 and DG subregions of young and middle-age non-treated rats in order to study age-related changes in spatial pattern separation (GEO Series GSE140420) [19]. We used these two datasets to evaluate the impact of the choice of reference genome in the context of hippocampal subregion comparisons in normal condition. For the latter dataset, we controlled the effects of difference of age in order to compare subregions (see Methods). For both datasets, we quantified gene expression by aligning the sequencing reads against the four *Rattus norvegicus* reference genomes: Ensembl 104 and 105 and RefSeq 106 and 108 and thus obtained a count matrix for each of them.

### 2.1.1 Mapping statistics

The use of annotations of the mRatBN7.2 genome assembly, *i.e.* RefSeq 108 and Ensembl 105 annotations, improves read mapping statistics in comparison with the previous assembly for both Ensembl and RefSeq annotations for the GSE179101 dataset (Figure 1).

In particular, the number of multi-mapped reads is reduced to the benefit of uniquely mapped reads. This disambiguation of mapping location results in an increase
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Read mapping statistics with respect to *Rattus norvegicus* reference genomes for the GSE179101 dataset. The raw sequencing reads of all the samples of GSE179101 dataset were trimmed for adpaters and quality before being aligned against *Rattus norvegicus* reference genomes and assigned to annotations.</figcaption></figure>

of the number of annotated reads, raising from 52.1% to 56.5%. When comparing
annotations of the same genome assembly, either Rnor_6.0 or mRatBN7.2, similar
results are obtained. In particular for the new assembly, Refseq annotations make
the percentage of annotated reads slightly increase from to 55.6% to 57.4%. Over-
all, we observed a similar improvement of mapping statistics using the annotations of
the mRatBN7.2 genome assembly, but for the benefit of Ensembl annotations, when
analyzing the GSE140420 dataset (Supplementary Figure 1).

2.1.2 Consistency and stringency

We used the consistency and stringency metrics introduced by Chen *et al* to measure
how different the expression quantification of the genes expressed in the GSE179101
dataset are for the evaluated pairs of reference genomes [14]. The Figures 2.A to 2.C
show the results of the comparison between RefSeq 108 and Ensembl 105 genomes
while the results of the pairwise comparisons including the Ensembl genomes and
Refseq genomes are showed in Supplementary Figures 2 and 3 respectively.

Matches between gene identifiers used in the two compared expression matrices
were identified since these metrics can be calculated only for genes expressed according
to both reference genomes. While this matching is straigthforward when comparing
two genomes from the same database since the same type of gene identifiers is used
by both genomes, some 1-to-many matches are found when comparing the RefSeq
108 and Ensembl 105 genomes due to unperfect relationships between Ensembl and
RefSeq gene identifiers (Figure 2.A). The majority of the genes were found with a
unique match in the other expression matrix for both genomes, representing 79.68%
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Consistency and stringency of *Rattus norvegicus* reference genomes for the analysis of GSE179101 dataset. A: gene identifier matching, genes with multiple matches and genes whose match is not expressed were gathered in the same category, B: consistency and stringency for all expressed genes in both matrices, C: distribution of mean expression of all the genes expressed in the matrix generated with Ensembl 105 genome, D: gene counts per category defined by consistency and stringency for pairwise genome comparison, E: distribution of cumulative exon length ratio in the annotations of two reference genomes for inconsistent genes.</figcaption></figure>

and 70.27% of the expressed genes in the matrices based on Ensembl 105 and RefSeq 108 genomes respectively.

Consistency and stringency were then calculated for all gene identifier matches between the two expression matrices of all the evaluated genome comparisons. The majority of the genes were found to be consistent, i.e. with less than 5% of difference in read counts over all the samples in average (Figure 2.B and Supplementary Figures 2 and 3). Among these consistent genes, the huge majority of them are stringent, with a significant part of them with identical counts for all the samples, while a few hundreds of them are non stringent. As Chen *et al* showed, most of the consistent and non stringent are very lowly expressed genes (Figure 2.C and Supplementary Figures 2 and 3). More than a third of the expressed genes in both matrices, representing sets of thousands of genes, are found to be inconsistent and possibly non stringent for all genome pairwise comparisons. Overall, the distribution of genes among the four categories defined by
the consistency and stringency metrics are similar for all the evaluated genome pair-wise comparisons (Figure 2.D). Chen *et al* showed that inconsistent genes are due to different numbers of exons and/or different exon lengths in the annotations of these genes. We confirmed this result with an overall trend of higher cumulative exon length in the annotations of the mRatBN7.2 genome assembly in comparison with those of the Rnor_6.0 genome assembly and in RefSeq 108 annotations in comparison with Ensembl 105 annotations (Figure 2.E).

Overall, we observed the same results with the GSE140420 dataset (Supplementary Figures 4 to 6), which suggests more complete annotations in the RefSeq 108 reference genome.

### 2.1.3 Differential expression analysis

A very common analysis using RNA-seq data consists in identifying genes whose mean expression change between sets of samples of biological interest in the so-called 'differential expression analysis'. We evaluated the impact of the choice of the reference genome on this type of analysis using the GSE179101 dataset. We identified differentially expressed (DE) genes for all "one subregion vs all others" comparisons. Based on the outcomes of mapping statistics and the results of the consistency and stringency analysis, we focused here on the comparison of the reference genomes based on the mRatBN7.2 assembly, *i.e.* Ensembl 105 and RefSeq 108 genomes (Additional file 2). The Figure 3 shows the comparison of results obtained with the Ensembl 105 and RefSeq 108 genomes for the comparison of CA1 subregion with CA2, CA3 and DG subregions while the results of the other subregion comparisons are shown in Supplementary Figures 7 to 9.

The majority of the DE genes are identified with both Ensembl 105 and RefSeq 108 genomes, representing 89.00% and 84.35% of the entire sets of DE genes respectively (Figure 3.A). Nevertheless, hundreds of DE genes are only found according to only one of the two reference genomes. Among these specific DE genes, most of them do not have a corresponding gene identifier, while the others are not identified as differentially expressed when using the other reference genome. We calculated the consistency and stringency for all the DE genes expressed in both matrices and found that most of the common DE genes are consistent and stringent while DE genes specific to one of the two reference genomes are mostly not annotated in the other reference genome annotation or inconsistent (Figure 3.B). We calculated the cumulative exon length of inconsistent common and specific genes and found a similar result to the one obtained with all the genes of these categories in the entire matrices: the cumulative exon length is higher in RefSeq 108 reference genome in comparison with Ensembl 105 reference genome (Figure 3.C). We observed similar results in the differential expression analysis for all other subregion comparisons (Supplementary Figures 7 to 9).

To assess the impact of the reference genome on the biological interpretation of the differential expression analysis outcomes, we performed a Gene Ontology (GO) term enrichment analysis among the sets of down-regulated and upregulated genes obtained with the Ensembl 105 and RefSeq 108 genomes for all subregion comparisons (Additional file 3). The differences in the differential expression analysis results generate GO
<figure><img src="image_5.png" /><figcaption>**Fig. 3** Impact of *Rattus norvegicus* reference genomes on the comparison of CA1 subregion to other subregions. Samples of the CA1, CA2, CA3 and DG subregions from the GSE179101 dataset were used for the differential expression analysis. A: gene identifier matching, genes with multiple matches and genes whose match is not expressed were gathered in the same category, B: differentially expressed (DE) gene counts per category defined by consistency and stringency, C: distribution of cumulative exon length ratio in the annotations of the two reference genomes for inconsistent DE genes, p-values of one-sided Wilcoxon signed-rank tests for paired difference in exon length where the alternative is a greater exon length in RefSeq 108 annotations are displayed per specificity and for inconsistent stringent genes (left) and inconsistent non stringent genes (right). D: enriched Gene Ontology term counts among down-regulated and upregulated genes.</figcaption></figure>

terms specifically found enriched for one of the reference genome (Figure 3.D). Even if most of the enriched GO terms are found regardless the reference genome used in the differential expression analysis (from 64.7% to 90.3% of the overall sets of enriched GO terms), dozens of enriched GO terms are found specifically to a single reference genome for some subregion comparisons (Supplementary Figures 7 to 9). Among these specifically enriched GO terms, we found that most of them are identified using the RefSeq 108 geneome. Overall, this significant amount of specific enriched GO terms suggests that the choice of the reference genome may have a dramatic impact on the biological interpretation of gene expression studies of rat hippocampal tissues. Similar results were observed for all the subregion comparisons enabled by the GSE140420 dataset (Additional files 4 and 5 and Supplementary Figures 10 to 12),
which confirms that the more complete annotations in RefSeq 108 reference genome impact the outcome of differential expression analysis for hippocampal subregion comparisons.

We explored into more details specifically enriched GO terms depending on the reference genome used for the differential expression analysis and focused on genome-specific differentially expressed genes related to these GO terms (Additional file 6). The Figure 4 shows the exon length log ratio of these genes and the p-value obtained with the other reference genome stratified by dysregulation, i.e. upregulation or down-regulation, and by genome specificity for the comparison of the CA1 subregion to the other three subregions.

<figure><img src="image_6.png" /><figcaption>**Fig. 4** Genome-specific differentially expressed genes related to specifically enriched Gene Ontology terms. Exon length log ratio and p-value obtained with the other reference genome for the comparison of CA1 subregion to all other subregions using samples from the GSE179101 dataset. Results are stratified by dysregulation and by genome specificity: specifically enriched GO terms identified among upregulated genes with RefSeq 108 genome (A), upregulated genes with Ensembl 105 genome (B), downregulated genes with RefSeq 108 genome (C), downregulated genes with Ensembl 105 genome (D). P-values of one-sided Wilcoxon signed-rank tests for paired difference in exon length where the alternative is a greater exon length or a lower exon length in RefSeq 108 annotations are displayed on top or bottom of the panels respectively for differentially expressed genes specifically identified with the two different genomes.</figcaption></figure>
Results obtained for the other subregion comparisons are showed in Supplementary Figures 13 to 15. For most combinations of dysregulation and genome specificity, differentially expressed genes specifically identified by one of the two compared genomes have longer exon length in RefSeq 108 in comparison with the annotations of the Ensembl 105 genome. Besides, the p-values obtained for the differential expression analysis with Ensembl 105 genome for differentially expressed genes specifically identified with RefSeq 108 these genes tend to be more noticeably nonsignificant. On the contrary, the p-values obtained for the differential expression analysis with the RefSeq 108 genome for the differentially expressed genes specifically identified using the Ensembl 105 genome are more close to the significance threshold set to 0.05. Since the exons of these genes are not longer in Ensembl 105 genome, this difference in differential expression outcome may be the consequence of a more stringent multiple testing p-value correction due to a higher number of genes (Figure 3.A and Supplementary Figures 7 to 9) rather than due to a difference in gene annotations. Taken together, these results shows that the difference in the annotations of the two genomes is more detrimental for the Ensembl 105 genome than for the RefSeq 108 genome. Overall, we observed the same trends in the results obtained with the GSE140420 dataset (Additional file 7 and Supplementary Figures 16 to 18).

Among specifically enriched GO terms depending on the reference genome used for the differential expression analysis identified for GSE179101 and GSE140420 datasets (Additional files 6 and 7), we focused on those with genome-specific differentially expressed genes having big differences in exon lengths and differential expression p-values and which are relevant in the context of hippocampal tissue analysis. Among them, we noted terms related to neurogenesis, e.g. "positive regulation of neurogenesis", downregulated in CA1 subregion, axons ("axon guidance" and "axon development") upregulated in CA1 subregion and calcium ion ("regulation of calcium ion transport" and "intracellular calcium ion homeostasis"), regulation of kinase activity ("positive regulation of MAP kinase activity" and "positive regulation of protein serine/threonine kinase activity"), "positive regulation of nervous system development" and "neuron projection development" downregulated in DG subregion. All of these terms were specifically enriched when using the RefSeq 108 genome.

## 2.2 Treated vs control condition comparison

After analyzing datasets of different hippocampal subregions in normal condition, we analyzed RNA-seq samples from a study with a classical "treated vs control" experimental design. O'Leary et al sequenced the rat dorsal hippocampal CA1 subregion in an early life seizure induced by kainic acid model in order to study neurocognitive deficits (GEO Series GSE136913) [20]. We reproduced the same analysis workflow as for the analysis of the previous datasets representing normal conditions and observed similar overall outcomes. Mapping statistics were improved when using the reference genomes based on the mRatBN7.2 assembly with a slightly higher amount of annotated reads when using the Ensembl 105 annotations in comparison with RefSeq 108 annotations (Supplementary Figure 19). Most of expressed genes which are present in both annotations have consistent read counts (Supplementary Figure 20). A third of the commonly expressed genes across the two annotations have inconsistent read
counts and most of these genes have longer exon length in the RefSeq 108 genes. The
results of the differential expression analysis confirms those obtained with datasets
for normal conditions even though the sets of differential expressed genes according
to each reference genome are much smaller in this study (Additional file 8 and Sup-
plementary Figure 21). Most of the differentially expressed genes for the “treated *vs*
control” comparison are identified by both reference genomes, differentially expressed
specific to one of the two reference genomes have longer exon length in the annotations
of RefSeq 108 genome. Regarding the functional analysis of these sets of differentially
expressed genes, a higher amount of GO terms are commonly identified (92.0% and
77.8% of the GO terms for downregulated and upregulated genes respectively) but
most of the specifically enriched GO terms are still identified using the RefSeq 108
genome (Additional file 9). Differentially expressed genes related to these GO terms
and specifically identified using the RefSeq 108 genome have longer exon length and
the p-values obtained with the Ensembl 105 genome for several of them are not close
to the significance threshold (Supplementary Figure 22). On the contrary, differen-
tially expressed genes related to these GO terms and specifically identified using the
Ensembl 105 genome do not have longer exons and the p-values obtained with the
RefSeq 108 genome are closed to 0.05. Similarly to the GSE179101 and GSE140420
datasets, we noticed relevant biological processes in the context of this study among
enriched GO terms specifically identified using the RefSeq 108, *e.g.* synaptic transmis-
sion (“modulation of chemical synaptic transmission”, “regulation of transmembrane
transport”), calcium transport (“positive regulation of calcium ion-dependent exocy-
tosis”), “generation of neurons” or “regulation of Ras protein signal transduction”
(Additional file 10).

3 Discussion

We showed the impact that the choice of reference genome may have on RNA-seq
data analysis. In particular, the outcomes of a differential expression analysis, either
in normal condition or in the context of a treated *vs* control comparison study, may
dramatically change depending on this choice. Some relevant dysregulated biological
processes may indeed be missed. Our results revealed that the RefSeq 108 genome is
the most complete reference genome of *Rattus norvegicus* species for the analysis of
hippocampal tissues to date.

## 3.1 Use of consistency and stringency to identify the reference genome to use in RNA-seq studies

3.1.1 Generalities

Initially, Chen *et al* introduced the consistency and stringency metrics to evaluate
the impact of two closely related reference genomes of *Cricetulus griseus* species (Chi-
nese hamster) on the analysis of RNA-seq data of ovary cells [14]. In this context, the
two metrics were used to reveal discrepancies in gene definition between the two com-
pared reference genomes. In particular, they questioned the reliability of differentially
expressed genes found to be inconsistent. Here, in the context of RNA-seq studies of
*Rattus norvegicus* hippocampal tissues, where several reference genomes are available, these two metrics revealed the higher completeness of the RefSeq 108 genome. First, the annotation of reads increased when using the reference genomes composed of either the Ensembl or RefSeq annotations of the mRatBN7.2 genome assembly thanks to significant improvements in sequencing technologies since the publication of the previous assembly. Second, we showed an overall higher cumulative exon length of inconsistent genes when comparing these two new reference genomes, highlighting more complete annotations in the RefSeq 108 genome.

### 3.1.2 Differential expression analysis results

While the consistency and stringency metrics are valuable measures to emphasize the differences in gene expression quantification due to the use of different reference genomes, some differential expression analysis results may appear confusing. For example, it is expected that most of the genes that are commonly found to be differentially expressed by using two different reference genomes would be consistent and stringent. Although this is the case for most common differentially expressed genes in our analysis (Figure 3 and Supplementary Figures 7 to 12 and 21), a significant part of these genes are inconsistent and possibly non stringent (from 27.6 to 36.5%). As Chen et al showed, the differences of annotations in the two compared reference genomes for these genes do not have consequences on the differential expression analysis outcome. Another surprising results deal with differentially expressed genes specifically identified when using one of the two reference genomes that are nevertheless found to be consistent and stringent. These results are explained by slight differences of annotations in the two compared reference genome and the effect of the FDR-controlling procedure which made the differential expression test p-value higher than 0.05 only for one of the two reference genomes, as it can be seen for some genes in Figure 4 and Supplementary Figures 13 to 18 and 22). Overall, all the examples discussed here are not the most frequent. Indeed, the categories according to the consistency and stringency metrics are mostly in agreement with the differential expression analysis results: most of the commonly identified differentially expressed genes are consistent and stringent and most of the differentially expressed genes specifically identified by only one reference genome are inconsistent.

## 3.2 RefSeq 108 genome as the new gold standard reference genome for *Rattus norvegicus*

### 3.2.1 Biological interpretation of differential expression analysis

When comparing the Ensembl and RefSeq annotations of the new mRatBN7.2 genome assembly, we showed that most of the enriched GO terms specifically identified with one of the two reference genomes are found using the RefSeq 108 geneome. Besides, we found that several of them are relevant in the context of the respective hippocampal tissue studies. We identified that “positive regulation of neurogenesis” GO term is enriched among downregulated genes when comparing samples from the CA1 subregion with those from other hippocampal subregions, which is in agreement with the fact that the DG subregion is a brain area where neurogenesis occurs [21, 22]. We
found that several GO terms related to kinase activity ("positive regulation of MAP
kinase activity", 'positive regulation of protein serine/threonine kinase activity") are
enriched among downregulated genes in the DG subregion in comparison with other
hippocampal subregions. A mitogen-activated protein kinase cascade is indeed known
to be activated in the CA1 and CA2 subregions but not in the CA3 and DG subregions
[23]. We also identified that "regulation of calcium ion transport" and "intracellular
calcium ion homeostasis" are specifically enriched among downregulated genes in the
DG subregion when using the RefSeq 108 genome. Calcium ion has a major role in the
synaptic transmission [24] and it was shown than the DG subregion has lower synapse
densities than the CA1 and CA3 subregions [25].

Regarding the analysis of the treated *vs* control comparison in the context of epilepsy, we also identified relevant GO terms among those specifically enriched among upregulated genes using the RefSeq 108 genome. In particular, we identified some GO terms related to the synaptic transmission. Dysfunction in the regulation of synaptic transmission plays a major role in epileptogenesis [26], especially in the CA1 subregion [27] which is the biological material analyzed in the GSE136913 dataset. In this context, similarly to the hippocampal subregion comparisons, finding biological processes related to calcium ion ("positive regulation of calcium ion-dependent exocytosis") also makes sense. Besides, we found the GO terms "generation of neurons" or "regulation of Ras protein signal transduction" among those specifically identified with the RefSeq 108 among upregulated genes while neurogenesis is aberrant in epilepsy [28] and recent studies proved the involvement of Ras pathway in epileptogenesis [29].

All these biological processes specifically identified when using the RefSeq 108 genome highlight that the differences in annotations in comparison with the Ensembl 105 genome are beneficial to the RefSeq 108 in the context of the different hippocampal studies analyzed here.

### 3.2.2 RefSeq 108 genome: the new gold standard reference genome for *Rattus norvegicus* hippocampal studies

Usually, RT-qPCR data are used as ground truth to validate RNA-seq data analysis
results. A positive correlation between the expression quantifications of genes obtained
with these two techniques is indeed required to confirm the results obtained with the
latter [11, 12]. Such data accompagnying the samples which compose the RNA-seq
datasets analyzed here are missing. Nevertheless, we consider the identical conclu-
sions of the reference genome evaluation using multiple RNA-seq datasets in different
contexts, i.e. improved read mapping statistics, longer exon length of inconsistent
genes and relevant specific differential expression analysis results, as solid evidences
toward the reference genome to choose. Taken together, these results make the RefSeq
108 genome the new gold standard reference genome for *Rattus norvegicus* RNA-seq
studies of hippocampal tissues.

### 3.2.3 Extension to other tissues

We believe that this approach should be reproduced for RNA-seq studies of any tissue when several reference genomes are available for the concerned species. Here, we
focused on the analysis of RNA-seq datasets from *Rattus norvegicus* hippocampal tissues. As Zhao *et al* showed for the human genome, the impact of the reference genome on gene expression quantification based on RNA-seq data can significantly vary depending on the analyzed tissues [9]. Thus, our approach should be extended to other tissues to eventually make the RefSeq 108 genome a new gold standard reference genome for all *Rattus norvegicus* RNA-seq studies.

### 3.2.4 RefSeq 108 genome is still unperfect

While the RefSeq 108 genome enables to reach a higher proportion of mapped reads and increases the proportion of annotated reads, a significant amount of them are mapped to genomic regions without any annotation (from 57.4 to 75.5% of reads aligned to annotated genomic regions over the three analyzed datasets). This gap may reveal that the annotations of this new reference genome still remain incomplete for the analysis of hippocampus RNA-seq data.

Waiting for the publication of new annotations of the mRatBN7.2 assembly or improvements in the sequencing of new genome assemblies, strategies have been developed to overcome the incompleteness of reference genomes, either dealing with assembly sequences or annotations. Morillon *et al* recommend to use a catalogue of *de novo* assemblies to reflect the transcript diversity [30]. A common strategy consists in taking benefits from data available for closely related species of the species of interest of a study. Davidson *et al* built a superTranscriptome by combining genome-guided and *de novo* assemblies for the species of interest and closely related species to solve the incompleteness of genome sequences or annotations [31]. Torres-Oliva *et al* reciprocally re-annotated the genomes of several *Drosophila* species in order to enable inter-specific comparative analysis [32]. Another strategy consists in building assemblies from large RNA-seq datasets and using classical annotation pipelines, e.g. BLASTx or Pfam. Ji *et al* implemented this approach to complete the *Rattus norvegicus* transcriptome by aligning the sequencing reads of 320 RNA-seq samples from 11 different organs and assembling them into transcripts [33]. They obtained a transcriptome with much more annotations than the version of RefSeq and Ensembl they compared to. Nevertheless, the assembly used for these annotations was the Rnor_6.0 assembly. A similar work with the mRatBN7.2 assembly would be of great interest. Similar approaches were applied to improve the annotations of the domestic pig and the zebrafish genomes [34, 35].

Overall, all these strategies should be considered in the context of gene expression quantification to improve the quality of RNA-seq data analysis, especially for tissues whose annotations are known to be incomplete in the available reference genomes of the species of interest. Nevertheless, these strategies should be carefully applied since errors and contamination in draft assemblies may lead to errors in annotations which possibly may be propagated across species [36].

## 4 Conclusion

The choice of reference genome in RNA-seq studies is an underexplored issue that may have dramatic effects on common downstream analysis, such as differential expression
analysis. Here, we focused on the analysis of RNA-seq datasets from *Rattus norvegicus*
hippocampal tissues and revealed that the reference genome in the release 108 of Ref-
Seq provides significant improvements in mapping statistics and more confident gene
expression quantification and differentially expressed gene identification. We propose
that our approach, based on the consistency and stringency metrics, should be repro-
duced for RNA-seq studies of any tissue to validate this genome as an overall new gold
standard reference genome for *Rattus norvegicus* RNA-seq studies. More generally, we
believe that the choice of reference genome in RNA-seq studies should be carefully
addressed when several reference genomes are available for a particular species.

# 5 Materials and Methods

## 5.1 RNA-seq datasets

All three RNA-seq datasets analyzed in our study were downloaded from the NCBI Gene
Expression Omnibus.

The GSE179101 dataset [18] contrains the RNA-seq data of the hippocampal subre-
gions CA1, CA2, CA3 and DG laser-microdissected from four rat brains. The extracted
RNA of the 16 samples was sequenced using the Illumina HiSeq 2500 platform.

The GSE140420 dataset [19] contains the RNA-seq data of the hippocampal subre-
gions CA1, CA3 and DG coming from young (5 months) and middle-age (12 months)
rats. It is composed of 84 samples coming from the three hippocampal subregions of
12 young rats and 16 middle-aged rats. Two sequencing platforms (Ion Torrent Pro-
ton and S5) were used to generate the entire RNA-seq dataset.

The GSE136913 dataset [20] contains the RNA-seq data of an early life seizure induced
by kainic acid model. It includes the sequencing data of the dorsal hippocampal CA1
subregion of 16 rats: 8 rats treated with kainate and 8 non treated rats. The control
group refers to the non treated rats while the treated group refers to rats with a seizure
induced by kainic acid. The rats belong to 8 litters composed of a control and a treated
rat in order to control the potential impact of litter variability on results. RNA was
sequenced in two rounds of 8 samples, coming from 4 control and 4 treated rats, using
two different sequencing platforms (Illumina HiSeq 2500 and 4000). For each round,
the samples were pooled over two lanes of sequencer, resulting in a RNA-seq dataset
composed of 32 sequencing samples.

## 5.2 *Rattus norvegicus* reference genomes

The Fasta files of the Rnor_6.0 and mRatBN7.2 genome assemblies and their
corresponding RefSeq annotations were downloaded from the NCBI FTP site in
the RefSeq releases 106 and 108 subdirectories respectively (https://ftp.ncbi.nlm.
nih.gov/genomes/all/annotation_releases/10116/106/GCF_000001895.5_Rnor_6.0/
and https://ftp.ncbi.nlm.nih.gov/genomes/all/annotation_releases/10116/108/
GCF_015227675.2_mRatBN7.2/). Ensembl annotations of the two genome assemblies
were downloaded from the Ensembl FTP site in releases 104 and 105 subdirecties
respectively (http://ftp.ensembl.org/pub/release-104/gtf/rattus_norvegicus/ and
http://ftp.ensembl.org/pub/release-105/gtf/rattus_norvegicus/).
## 5.3 Expression quantification

Trim Galore was used to trim raw sequencing reads for adapters and quality with default parameters [37]. The trimmed and filtered reads were then aligned against the reference genomes using STAR [38] and counted using htseq-count [39]. After mapping, Illumina paired-end reads were counted using htseq-count -s option set to 'reverse'.

## 5.4 Annotation file processing

Matches between gene identifiers used in Ensembl and RefSeq annotations were found using the biomaRt R package [40]. The 'gene', 'transcript' and 'exon' keys were used for gene, transcript and exon annotations respectively. The number of features of an annotation type and the cumulative length of an annotation type were calculated after merging overlapping genomic ranges with the reduce() function of the GenomicRanges R package for each annotation file [41]. Gene exon lengths according to two different annotations were compared using a one-sided Wilcoxon signed-rank test where the alternative is a lower exon length according to one of the two compared annotations.

## 5.5 Consistency and stringency

Consistency and stringency are two metrics based on the ratio $R_{ij}$ of read counts for a gene $i$ derived from the mapping of sequencing reads of a sample $j$ against two different reference genomes [14]:

$$R_{ij} = \log_2 \left( \frac{C_{ij}^1 + n}{C_{ij}^2 + n} \right)$$

where,

* $C_{ij}^1$ refers to the read count obtained with a first reference genome and $C_{ij}^2$ refers to the read count obtained with a second reference genome;

* $n$ is an integer value added to all the read counts to avoid null values ($n = 1$).

Contrary to Chen et al, we set $n$ to 1 instead of 8 to minimize the impact of the addition of this value, especially for low read counts.

The consistency is defined by the absolute value of the ratios over all the samples which compose the RNA-seq dataset while the stringency is defined by the standard deviation:

$$\text{consistency}_i = - \log_{10} (|\text{mean}(R_{ij})|)$$

$$\text{stringency}_i = - \log_{10} (\text{sd}(R_{ij}))$$

The minimum non-null absolute value of all $R_{ij}$ was added to null values to avoid infinite values after log-transformation. A threshold corresponding to a difference of 5% in read counts between the two compared reference genome was used to define categories according the consistency and stringency metrics.
## 5.6 Differential expression analysis

Lowly expressed genes were filtered out using a threshold of 5 counts per million mapped reads and read counts were normalized using the TMM method [42] prior to differential expression identification.

Differentially expressed genes were identified using the edgeR package [43]. Specific designs were used to take benefits from the metadata of each RNA-seq dataset. For the GSE136913 dataset, we used a condition factor of interest with two levels to make comparisons between treated and control samples. Variations due to litters and sequencing platforms were handled by supplying the litter to which rats belong and the sequencers used to generate the RNA-seq samples as blocking factors in the generalized linear model (GLM). For the GSE140420 dataset, the age of the rats was included as a blocking factor in the GLM to account for age-specific effects. One *vs* all other subregions comparisons were made for both GSE179101 and GSE140420 datasets. Differentially expressed genes were identified using a fold-change threshold of 1 and thanks to a quasi-likelihood F-test p-value below 0.05 after correction by the Benjamini-Hochberg procedure [44] for the coefficient in the GLM corresponding to the comparison of interest.

To perform Gene Ontology term over-representation analysis, Ensembl and RefSeq identifiers of differentially expressed genes were converted into Entrez identifiers using both the biomaRt [40] and the org.Rn.eg.db [45] R packages. The org.Rn.eg.db and biomaRt R packages were used as the main conversion source for RefSeq and Ensembl identifiers respectively due to their respective higher performances. In case of disagreement between the two packages, biomaRt Entrez identifier was kept. Enriched Gene Ontology terms were identified among differentially expressed genes thanks to a one-sided hypergeometric test p-value below 0.05 after correction by the Benjamini-Yekutieli procedure [46].

**Supplementary information.** Additional file 1:

* filename: supplementary_figures.pdf

* file format: PDF (.pdf extension)

* title: Supplementary figures

* description: File containing all supplementary figures.

Additional file 2:

* filename: additional_file_2

* file format: Microsoft Excel 2007-365 (.xlsx extension)

* title: Differentially expressed genes identified for all "one subregion *vs* all others" comparisons using Ensembl 105 and RefSeq 108 genomes for the GSE179101 dataset.

* description: One sheet per reference genome and subregion comparison.

Additional file 3:

* filename: additional_file_3

* file format: Microsoft Excel 2007-365 (.xlsx extension)
* title: Enriched GO terms among differentially expressed genes identified for all "one subregion vs all others" comparisons using Ensembl 105 and RefSeq 108 genomes for the GSE179101 dataset.

* description: One sheet per reference genome, dysregulation, i.e. upregulation or downregulation, and subregion comparison.

Additional file 4:

* filename: additional_file_4

* file format: Microsoft Excel 2007-365 (.xlsx extension)

* title: Differentially expressed genes identified for all "one subregion vs all others" comparisons using Ensembl 105 and RefSeq 108 genomes for the GSE140420 dataset.

* description: One sheet per reference genome and subregion comparison.

Additional file 5:

* filename: additional_file_5

* file format: Microsoft Excel 2007-365 (.xlsx extension)

* title: Enriched GO terms among differentially expressed genes identified for all "one subregion vs all others" comparisons using Ensembl 105 and RefSeq 108 genomes for the GSE140420 dataset.

* description: One sheet per reference genome, dysregulation, i.e. upregulation or downregulation, and subregion comparison.

Additional file 6:

* filename: additional_file_6

* file format: Microsoft Excel 2007-365 (.xlsx extension)

* title: Genome-specific differentially expressed genes related to specifically enriched GO terms depending on the reference genome used for the differential expression analysis for the GSE179101 dataset.

* description: One sheet per reference genome, dysregulation, i.e. upregulation or downregulation, and subregion comparison.

Additional file 7:

* filename: additional_file_7

* file format: Microsoft Excel 2007-365 (.xlsx extension)

* title: Genome-specific differentially expressed genes related to specifically enriched GO terms depending on the reference genome used for the differential expression analysis for the GSE140420 dataset.

* description: One sheet per reference genome, dysregulation, i.e. upregulation or downregulation, and subregion comparison.

Additional file 8:

* filename: additional_file_8

* file format: Microsoft Excel 2007-365 (.xlsx extension)

* title: Differentially expressed genes identified for the "treated vs control samples" comparison using Ensembl 105 and RefSeq 108 genomes for the GSE136913 dataset
* description: One sheet per reference genome.

Additional file 9:

* filename: additional_file_9

* file format: Microsoft Excel 2007-365 (.xlsex extension)

* title: Enriched GO terms among differentially expressed genes identified for the "treated vs control samples" comparison using Ensembl 105 and RefSeq 108 genomes for the GSE136913 dataset.

* description: One sheet per reference genome and dysregulation, i.e. upregulation or downregulation.

Additional file 10:

* filename: additional_file_10

* file format: Microsoft Excel 2007-365 (.xlsex extension)

* title: Genome-specific differentially expressed genes related to specifically enriched GO terms depending on the reference genome used for the differential expression analysis for the GSE136913 dataset.

* description: One sheet per reference genome and dysregulation, i.e. upregulation or downregulation.

## Declarations

Some journals require declarations to be submitted in a standardised format. Please check the Instructions for Authors of the journal to which you are submitting to see if you need to complete this section. If yes, your manuscript must contain the following sections under the heading 'Declarations':

* Funding: This work was supported by the Institut National pour la Santé et la Recherche Médicale (Inserm, France), Université Paris Cité, Université Sorbonne Paris Nord, the French National Research Agency (#ANR-18-CE17-0009-01)(C.L.P., A.D.-D).

* Conflict of interest/Competing interests: The authors declare that they have no competing interests.

* Ethics approval: Not applicable.

* Consent to participate: Not applicable.

* Consent for publication: Not applicable.

* Availability of data and materials: The datasets analysed during the current study are available in NCBI Gene Expression Omnibus under the accession numbers GSE136913 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE136913), GSE140420 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE140420) and GSE179101 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE179101). The Fasta files of the Rnor_6.0 and mRatBN7.2 genome assemblies and their corresponding RefSeq annotations were downloaded from the NCBI FTP site in the RefSeq releases 106 and 108 subdirectories respectively (https://ftp.ncbi.nlm.nih.gov/genomes/all/annotation_releases/10116/106/GCF_000001895.5_Rnor_6.0/ and https://ftp.ncbi.nlm.nih.gov/genomes/all/
annotation_releases/10116/108/GCF_015227675.2_mRatBN7.2/). Ensembl annota-
tions of the two genome assemblies were downloaded from the Ensembl FTP
site in releases 104 and 105 subdirecties respectively (http://ftp.ensembl.org/pub/
release-104/gtf/rattus_norvegicus/ and http://ftp.ensembl.org/pub/release-105/
gtf/rattus_norvegicus/).

* Code availability: The source code and data used to produce the results and analyses presented in this manuscript are available at https://github.com/lepriolc/R_norvegicus_ref_genome_evalution and we have archived our code on Zenodo (DOI: 10.5281/zenodo.8099337).

* Authors' contributions: Conceptualization: CLP, Methodology: CLP, Software: CLP, Validation: CLP, Formal analysis: CLP, Investigation: CLP, Resources: ADD, Writing - Original Draft: CLP, Writing - Review & Editing: ADD, Visualization: CLP, Funding acquisition: ADD

If any of the sections are not relevant to your manuscript, please include the heading
and write 'Not applicable' for that section.

Editorial Policies for:

Springer journals and proceedings: https://www.springer.com/gp/editorial-policies

Nature Portfolio journals: https://www.nature.com/nature-research/editorial-policies

Scientific Reports: https://www.nature.com/srep/journal-policies/editorial-policies

BMC journals: https://www.biomedcentral.com/getpublished/editorial-policies

References

[1] Srivastava, A., Malik, L., Sarkar, H., Zakeri, M., Almodaresi, F., Soneson, C., Love, M.I., Kingsford, C., Patro, R.: Alignment and mapping methodology influence transcript abundance estimation. Genome Biol 21(1), 239 (2020)

[2] Soneson, C., Delorenzi, M.: A comparison of methods for differential expression analysis of RNA-seq data. BMC Bioinformatics 14, 91 (2013)

[3] Richardson, J.E.: fjoin: simple and efficient computation of feature overlaps. J Comput Biol 13(8), 1457-1464 (2006)

[4] Standage, D.S., Brendel, V.P.: ParsEval: parallel comparison and analysis of gene structure annotations. BMC Bioinformatics 13, 187 (2012)

[5] Frankish, A., Uszczynska, B., Ritchie, G.R., Gonzalez, J.M., Pervouchine, D., Petryszak, R., Mudge, J.M., Fonseca, N., Brazma, A., Guigo, R., Harrow, J.: Comparison of GENCODE and RefSeq gene annotation and the impact of reference geneset on variant effect prediction. BMC Genomics 16 Suppl 8, 2 (2015)
[6] Li, H., Dawood, M., Khayat, M.M., Farek, J.R., Jhangiani, S.N., Khan, Z.M., Mitani, T., Coban-Akdemir, Z., Lupski, J.R., Venner, E., Posey, J.E., Sabo, A., Gibbs, R.A.: Exome variant discrepancies due to reference-genome differences. Am J Hum Genet 108(7), 1239–1250 (2021)

[7] Arora, S., Pattwell, S.S., Holland, E.C., Bolouri, H.: Variability in estimated gene expression among commonly used RNA-seq pipelines. Sci Rep 10(1), 2734 (2020)

[8] Bruning, R.S., Tombor, L., Schulz, M.H., Dimmeler, S., John, D.: Comparative analysis of common alignment tools for single-cell RNA sequencing. Gigascience 11 (2022)

[9] Zhao, S., Zhang, B.: A comprehensive evaluation of ensembl, RefSeq, and UCSC annotations in the context of RNA-seq read mapping and gene quantification. BMC Genomics 16, 97 (2015)

[10] Slabaugh, E., Desai, J.S., Sartor, R.C., Lawas, L.M.F., Jagadish, S.V.K., Doherty, C.J.: Analysis of differential gene expression and alternative splicing is significantly influenced by choice of reference genome. RNA 25(6), 669–684 (2019)

[11] Chisanga, D., Liao, Y., Shi, W.: Impact of gene annotation choice on the quantification of RNA-seq data. BMC Bioinformatics 23(1), 107 (2022)

[12] Hamaguchi, Y., Zeng, C., Hamada, M.: Impact of human gene annotations on RNA-seq differential expression analysis. BMC Genomics 22(1), 730 (2021)

[13] Wu, P.Y., Phan, J.H., Wang, M.D.: Assessing the impact of human genome annotation choice on RNA-seq expression estimates. BMC Bioinformatics 14 Suppl 11, 8 (2013)

[14] Chen, C., Le, H., Goudar, C.T.: Evaluation of two public genome references for chinese hamster ovary cells in the context of rna-seq based gene expression analysis. Biotechnol Bioeng 114(7), 1603–1613 (2017)

[15] Smith, J.R., Bolton, E.R., Dwinell, M.R.: The Rat: A Model Used in Biomedical Research. Methods Mol Biol 2018, 1–41 (2019)

[16] Howe, K., Dwinell, M., Shimoyama, M., Corton, C., Betteridge, E., Dove, A., Quail, M.A., Smith, M., Saba, L., Williams, R.W., Chen, H., Kwitek, A.E., McCarthy, S.A., Uliano-Silva, M., Chow, W., Tracey, A., Torrance, J., Sims, Y., Challis, R., Threlfall, J., Blaxter, M.: The genome sequence of the Norway rat, Rattus norvegicus Berkenhout 1769. Wellcome Open Res 6, 118 (2021)

[17] Jong, T.V., Chen, H., Brashear, W.A., Kochan, K.J., Hillhouse, A.E., Zhu, Y., Dhande, I.S., Hudson, E.A., Sumlut, M.H., Smith, M.L., Kalbfleisch, T.S., Doris, P.A.: mRatBN7.2: familiar and unfamiliar features of a new rat genome reference assembly. Physiol Genomics 54(7), 251–260 (2022)
[18] Machado, J.P.D., Athie, M.C.P., Matos, A.H.B., Lopes-Cendes, I., Vieira, A.S.: The transcriptome of rat hippocampal subfields. IBRO Neurosci Rep **13**, 322–329 (2022)

[19] Smith, G., Rani, A., Kumar, A., Barter, J., Foster, T.C.: Hippocampal Subregion Transcriptomic Profiles Reflect Strategy Selection during Cognitive Aging. J Neurosci **40**(25), 4888–4899 (2020)

[20] O'Leary, H., Vanderlinden, L., Southard, L., Castano, A., Saba, L.M., Benke, T.A.: Transcriptome analysis of rat dorsal hippocampal CA1 after an early life seizure induced by kainic acid. Epilepsy Res **161**, 106283 (2020)

[21] Ninkovic, J., Mori, T., tz, M.: Distinct modes of neuron addition in adult mouse neurogenesis. J Neurosci **27**(40), 10906–10911 (2007)

[22] Abbott, L.C., Nigussie, F.: Adult neurogenesis in the mammalian dentate gyrus. Anat Histol Embryol **49**(1), 3–16 (2020)

[23] Blum, S., Moore, A.N., Adams, F., Dash, P.K.: A mitogen-activated protein kinase cascade in the CA1/CA2 subfield of the dorsal hippocampus is essential for long-term spatial memory. J Neurosci **19**(9), 3535–3544 (1999)

[24] Mateos-Aparicio, P., guez-Moreno, A.: Calcium Dynamics and Synaptic Plasticity. Adv Exp Med Biol **1131**, 965–984 (2020)

[25] Santuy, A., s-Roca, L., guez, J.R., lez-Soriano, J., Zhu, F., Qiu, Z., Grant, S.G.N., DeFelipe, J., Merchan-Perez, A.: Estimation of the number of synapses in the hippocampus and brain-wide by volume electron microscopy and genetic labeling. Sci Rep **10**(1), 14014 (2020)

[26] Casillas-Espinosa, P.M., Powell, K.L., O'Brien, T.J.: Regulators of synaptic transmission: roles in the pathogenesis and treatment of epilepsy. Epilepsia **53 Suppl 9**, 41–58 (2012)

[27] Owen, B., Bichler, E., Benveniste, M.: Excitatory synaptic transmission in hippocampal area CA1 is enhanced then reduced as chronic epilepsy progresses. Neurobiol Dis **154**, 105343 (2021)

[28] Jessberger, S., Parent, J.M.: Epilepsy and Adult Neurogenesis. Cold Spring Harb Perspect Biol **7**(12) (2015)

[29] Sran, S., Bedrosian, T.A.: RAS pathway: The new frontier of brain mosaicism in epilepsy. Neurobiol Dis **180**, 106074 (2023)

[30] Morillon, A., Gautheret, D.: Bridging the gap between reference and real transcriptomes. Genome Biol **20**(1), 112 (2019)

[31] Davidson, N.M., Oshlack, A.: Necklace: combining reference and assembled
transcriptomes for more comprehensive RNA-Seq analysis. Gigascience **7**(5) (2018)

[32] Torres-Oliva, M., Almudi, I., McGregor, A.P., Posnien, N.: A robust (re-)annotation approach to generate unbiased mapping references for RNA-seq-based analyses of differential expression across closely related species. BMC Genomics **17**, 392 (2016)

[33] Ji, X., Li, P., Fuscoe, J.C., Chen, G., Xiao, W., Shi, L., Ning, B., Liu, Z., Hong, H., Wu, J., Liu, J., Guo, L., Kreil, D.P., Łabaj, P.P., Zhong, L., Bao, W., Huang, Y., He, J., Zhao, Y., Tong, W., Shi, T.: A comprehensive rat transcriptome built from large scale RNA-seq-based annotation. Nucleic Acids Res **48**(15), 8320–8331 (2020)

[34] Beiki, H., Liu, H., Huang, J., Manchanda, N., Nonneman, D., Smith, T.P.L., Reecy, J.M., Tuggle, C.K.: Improved annotation of the domestic pig genome through integration of Iso-Seq and RNA-seq data. BMC Genomics **20**(1), 344 (2019)

[35] Lawson, N.D., Li, R., Shin, M., Grosse, A., Yukselen, O., Stone, O.A., Kucukural, A., Zhu, L.: An improved zebrafish transcriptome annotation for sensitive and comprehensive detection of cell type-specific genes. Elife **9** (2020)

[36] Salzberg, S.L.: Next-generation genome annotation: we still struggle to get it right. Genome Biol **20**(1), 92 (2019)

[37] Krueger, F.: TrimGalore, 0.6.10, (2023). https://doi.org/10.5281/zenodo.7598955

[38] Dobin, A., Davis, C.A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., Batut, P., Chaisson, M., Gingeras, T.R.: STAR: ultrafast universal RNA-seq aligner. Bioinformatics **29**(1), 15–21 (2013)

[39] Anders, S., Pyl, P.T., Huber, W.: HTSeq–a Python framework to work with high-throughput sequencing data. Bioinformatics **31**(2), 166–169 (2015)

[40] Drost, H.G., Paszkowski, J.: Biomartr: genomic data retrieval with R. Bioinformatics **33**(8), 1216–1217 (2017)

[41] Lawrence, M., Huber, W., s, H., Aboyoun, P., Carlson, M., Gentleman, R., Morgan, M.T., Carey, V.J.: Software for computing and annotating genomic ranges. PLoS Comput Biol **9**(8), 1003118 (2013)

[42] Robinson, M.D., Oshlack, A.: A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol **11**(3), 25 (2010)

[43] Robinson, M.D., McCarthy, D.J., Smyth, G.K.: edgeR: a Bioconductor package
for differential expression analysis of digital gene expression data. Bioinformatics **26**(1), 139–140 (2010)

[44] Benjamini, Y., Hochberg, Y.: Controlling the false discovery rate: A practical and powerful approach to multiple testing. Journal of the Royal Statistical Society. Series B (Methodological) **57**(1), 289–300 (1995)

[45] Carlson, M.: org.Rn.eg.db: Genome Wide Annotation for Rat, 3.10.0, (2019)

[46] Benjamini, Y., Yekutieli, D.: The control of the false discovery rate in multiple testing under dependency. Ann. Statist. **29**(4), 1165–1188 (2001)
Supplementary Files

This is a list of supplementary files associated with this preprint. Click to download.

* additionalfile1.pdf

• additionalfile2.xlsx

• additionalfile3.xlsx

• additionalfile4.xlsx

• additionalfile5.xlsx

• additionalfile6.xlsx

• additionalfile7.xlsx

• additionalfile8.xlsx

• additionalfile9.xlsx

* additionalfile10.xlsx