Animal-eRNAdb: a comprehensive animal enhancer RNA database.

Enhancer RNAs (eRNAs) are a class of non-coding RNAs transcribed from enhancers. As the markers of active enhancers, eRNAs play important roles in gene regulation and are associated with various complex traits and characteristics. With increasing attention to eRNAs, numerous eRNAs have been identified in different human tissues. However, the expression landscape, regulatory network and potential functions of eRNAs in animals have not been fully elucidated. Here, we systematically characterized 185 177 eRNAs from 5085 samples across 10 species by mapping the RNA sequencing data to the regions of known enhancers. To explore their potential functions based on evolutionary conservation, we investigated the sequence similarity of eRNAs among multiple species. In addition, we identified the possible associations between eRNAs and transcription factors (TFs) or nearby genes to decipher their possible regulators and target genes, as well as characterized trait-related eRNAs to explore their potential functions in biological processes. Based on these findings, we further developed Animal-eRNAdb (http://gong_lab.hzau.edu.cn/Animal-eRNAdb/), a user-friendly database for data searching, browsing and downloading. With the comprehensive characterization of eRNAs in various tissues of different species, Animal-eRNAdb may greatly facilitate the exploration of functions and mechanisms of eRNAs.

INTRODUCTION

Enhancers are a type of distal regulatory genomic element, which could couple with promoters to establish an enhancer–promoter loop to determine spatiotemporal and quantitative transcription of genes in response to developmental or environmental cues (1,2). Recent advances in transcriptomic and (epi)genomic research have revealed that active enhancers can open local chromatin structure and expose the DNA motifs to attract transcription factors (TFs); these TFs further recruit RNA polymerases to generate non-coding RNAs (ncRNAs), which are defined as enhancer RNAs (eRNAs) (3). In early studies, eRNAs were only regarded as a by-product of enhancer transcription without specific functions (4). Traditionally, the identification of enhancer activity mainly depends on the binding of TFs, histone modifications in the enhancer region and chromatin accessibility (5). With increasing attention to eRNAs, more studies show that the production of eRNAs is a widespread phenomenon and eRNAs are highly related to enhancer activity (6–8). For example, inactive enhancers often show lower expression levels of PRO-cap and CAGE signals compared with active enhancers (6), and eRNAs can promote the formation and dynamic stabilization of the enhancer–promoter loop (7,9,10). Accumulating evidence shows that the expression levels of eRNAs are associated with multiple traits and characteristics (11). For example, the eRNA OLMALINC can influence body weight by regulating the gene stearoyl-coenzyme A desaturase (12). Furthermore, eRNAs are regulated by TFs and play an important role in regulating gene expression (8,13,14). For example, estrogen receptor 1 (ESR1) induces thousands of eRNAs to maintain transcriptional circuitry in breast cancer (9), and an eRNA transcribed from a distal regulatory MyoD enhancer can mediate cohesin recruitment and promote myogenin gene expression during myogenic differentiation (15).

Due to the importance of eRNAs, numerous eRNAs have been identified across various human tissues (16), and several human eRNA databases have been developed. Andersson et al. systematically annotated ∼65 000 eRNAs by the FANTOM project in ∼400 human tissues and cell types using the cap analysis of gene expression (CAGE-seq) technique targeting the molecules with 5’cap (17). Besides humans, eRNAs also play important roles in other eukaryotic animals. For example, antisense oligonucleotides targeting hypoxia-inducible eRNA (HERNA) can protect mice from stress-induced pathological hypertrophy (18). However, the expression landscape of eRNAs in animals has not been fully elucidated and few databases have been developed for animal eRNAs. In addition, many enhancers are highly constrained and some are even conserved across long evolutionary distances (19). The function of enhancers is usually explored based on conservation (20). Hence, further research on conserved eRNAs can help to build a better understanding of the potential functions of enhancers and eRNAs across different species. In summary, it is highly valuable to comprehensively investigate the expression landscape and conservation of eRNAs to uncover the mechanisms underlying the regulation of gene expression and the phenotypes of animals.

With the development of high-throughput sequencing technology, several techniques can be used to identify eRNAs, such as global nuclear run-on sequencing (GRO-seq) (21), cap analysis of gene expression (CAGE-seq) (22), precision nuclear run-on sequencing (PRO-seq) (23), chromatin immunoprecipitation sequencing (ChIP-seq) (24) and routine RNA-seq. The rapid development of RNA-seq technology in the past few decades has contributed to the accumulation of large amounts of data, which will facilitate the convenient and efficient studies of gene regulation, transcript structure and ncRNAs. Although the sensitivity of RNA-seq in detecting eRNAs is lower than that of other techniques, the relatively abundant data and low costs make it possible to study animal eRNAs at the genome-wide level through RNA sequencing (2). In humans, a large number of detectable eRNAs have been identified and a high-resolution map of eRNA loci has been generated using RNA-seq data respectively from the Genotype-Tissue Expression (GTEx) project (11) and The Cancer Genome Atlas (TCGA) (3).

In this study, we collected vast amounts of RNA-seq data from 10 species, including chimpanzee, rhesus, mouse, rat, sheep, chicken, clawed frog, zebrafish, fruitfly and worm from public databases, as well as their annotations from SEA 3.0 (25) and EnhancerAtlas 2.0 (26). By integrating these datasets, we systematically characterized the eRNA profiles in 5085 samples of the 10 species. In addition, we analyzed the correlations between eRNAs and traits/TFs/genes to find the possible trait-related eRNAs as well as the putative eRNA regulators and target genes. To uncover their potential evolutionary conservation, we investigated the sequence similarity of eRNAs among multiple species using blastn. Finally, we developed Animal-eRNAdb (http://gong_lab.hzau.edu.cn/Animal-eRNAdb/), a user-friendly database for the browsing, searching and downloading of eRNA-related information.

MATERIALS AND METHODS

Collection and processing of data and identification of eRNAs

Annotations of enhancers were collected from SEA 3.0 (http://sea.edbc.org/) and EnhancerAtlas 2.0 (http://www.enhanceratlas.org/indexv2.php), and then adjusted using LiftOver (27) according to the corresponding genome version of the species (28) (Figure 1). We defined ± 3 kb around the middle loci of the enhancer as the eRNA region (11,16,29). To avoid the potential influence of known transcripts, we discarded the eRNAs which overlapped with known annotations including protein-coding RNAs and ncRNAs (e.g. lncRNA, pseudogene and snoRNA) and whose length was <6001 bp. The RNA-seq data of animal samples were firstly downloaded from the Sequence Read Archive (SRA, https://www.ncbi.nlm.nih.gov/sra) of the National Center for Biotechnology Information (30,31). Then, we extracted certain species with eRNA annotation files. As a result, a total of 5085 samples of 10 species were selected for further study. The raw RNA-seq data were downloaded, converted into standard fastq, subjected to quality control using FastQC (version: v0.11.8), cleaned with Trim Galore (version: 0.6.4_dev), and then aligned to the corresponding reference genome using HISAT2 (32). Subsequently, we captured those RNA-seq reads mapping on the regions of eRNAs by SAMtools (version: 1.11) (33) and calculated the read counts per gene using FeatureCounts (version: v2.0.1) (34). Then, we used Reads Per Million (RPM) to normalize the eRNA expression (35) and used transcripts per million (TPM) to normalize the gene expression (36). Detectable eRNAs were defined as those with average expression values >1 (RPM ≥ 1) in at least one tissue in one bioproject. All codes and scripts used in this study are available upon reasonable request.

Figure 1.

Flow chart of Animal-eRNAdb.

Identification of trait-related eRNAs

The information of traits, including gender and developmental stage, was also downloaded from SRA. Then, we calculated the association between the expression of individual eRNAs and each trait across tissues for each species (11,37). We used Student’s t-test to assess the statistical difference of eRNAs between male and female samples and defined |fold change (FC)| ≥ 1.5 and false discovery rate (FDR) < 0.05 as statistical significance. Considering the differences between embryo and postnatal samples, we divided the tissues into three types: tissues with both embryo and postnatal samples, tissues with only embryo samples and tissues with only postnatal samples. For tissues with only embryo or postnatal samples, we used Spearman’s correlation to assess the association between each eRNA expression and the developmental stage if the developmental index is numerical variables and defined |Rho| ≥ 0.3 and FDR < 0.05 as statistical significance. For tissues categorized by developmental stage, the statistical difference of eRNAs was evaluated by Student’s t-test for dichotomous variables (|FC| ≥ 1.5 and FDR < 0.05) and the analysis of variance (ANOVA) test for variables presenting more than two categories (FDR < 0.05) (11). For tissues with both embryo and postnatal samples, we first used Student’s t-test to identify the differentially expressed eRNAs between the embryo and postnatal samples. Then, the same methods were used to identify differentially expressed eRNAs for embryo and postnatal samples, respectively.

Identification of putative regulators and target genes of eRNAs

Putative regulators of eRNAs were defined as TFs that have significant associations with eRNA expression. Annotations of TFs were collected from AnimalTFDB (http://bioinfo.life.hust.edu.cn/AnimalTFDB/) (38), and the expression data of TFs were extracted from the above total gene expression matrix. Putative regulators of eRNAs were identified based on the co-expression between eRNAs and TFs across tissues for each species. We used Spearman’s correlation to evaluate the co-expression and defined |Rho| ≥ 0.3 and FDR < 0.05 as statistical significance.

Putative target genes of eRNAs were defined as genes that are relatively close to eRNAs (distance ≤ 1 MB) and have significant co-expression with eRNAs (Spearman’s correlation |Rho| ≥ 0.3 and FDR < 0.05) in each tissue of species. We removed those eRNAs-target pairs in which eRNAs are located in the intronic region of the target genes following previous studies (39,40).

Identification of sequence similarity of eRNAs

Considering the conservation of enhancers, we evaluated the similarity degree to uncover the potential evolutionary conservation of eRNAs in different species. First, we extracted the sequences of each eRNA per species from genome.fa using Bedtools (version: v2.29.2) (41) and downloaded human eRNA sequences from HeRA (11). Subsequently, we calculated the similarity of each eRNA in a given species to all eRNAs in other species through blastn and defined similarity ≥ 0.5 and expectation value (E-value) ≤ 1e-5 as statistical significance.

IMPLEMENTATION

Animal-eRNAdb (http://gong_lab.hzau.edu.cn/Animal-eRNAdb/) was built based on the THINKPHP framework (http://www.thinkphp.cn/) and Bootstrap 4 (https://getbootstrap.com/), which runs on the Apache 2 web server (https://httpd.apache.org/) with MySQL (https://www.mysql.com/) as its database engine and R (https://www.r-project.org/) for graph drawing. Animal-eRNAdb is available online without registration and optimized for Chrome (recommended), Internet Explorer, Opera, Firefox, Windows Edge and macOS Safari.

DATABASE CONTENT AND USAGE

Samples in Animal-eRNAdb

In total, 5085 samples across 187 tissues of 10 species were included in Animal-eRNAdb, ranging from 199 samples in zebrafish to 901 samples in rat (Table 1), and from one tissue in frog to 73 tissues in sheep. The detailed information, including the number of samples per species, reference genome versions and the number of eRNAs, is available on the ‘Document’ page. As a user-friendly data portal, Animal-eRNAdb displays eRNA-related information across a large number of tissues of various species.

Table 1.

eRNAs identified in Animal-eRNAdb

			eRNAs
Species	No. of samples	No. of tissues	No. of total eRNAs	eRNAs in ≥ 2 tissues	eRNAs in ≥ 3 tissues
Gallus gallus (Chicken)	656	15	24 067	0.62	0.43
Pan troglodytes (Chimp)	262	23	11 995	0.73	0.60
Drosophila melanogaster (Fruitfly)	774	3	1807	0.53	0.21
Xenopus tropicalis (Clawed frog)	284	1	1081	-	-
Mus musculus (Mouse)	703	35	72 956	0.47	0.33
Rattus norvegicus (Rat)	901	16	31 871	0.71	0.57
Macaca mulatta (Rhesus)	257	13	12 913	0.53	0.38
Ovis aries (Sheep)	730	73	5697	0.77	0.65
Caenorhabditis elegans (Worm)	319	1	505	-	-
Danio rerio (Zebrafish)	199	7	22 285	0.52	0.36
Summation	5,085	187	185 177	0.57	0.42
Minimun	199	1	505	0.47	0.21
Maximun	901	73	72 956	0.77	0.65
Median	488	14	12 454	0.57	0.41

eRNA expression landscape in Animal-eRNAdb

We identified a total of 185 177 eRNAs in these species, ranging from 505 in worm to 72 956 in mouse at species level, and from 276 in the psoas major muscle of sheep to 43 690 in the testis of mouse at the tissue level. Many eRNAs are expressed in multiple tissues, which is consistent with previous studies (3,16,42). For example, a median of 57% eRNAs were expressed in more than one tissue, with the maximum of 77% in sheep, and a median of 41% eRNAs were expressed in three or more tissues, with the maximum of 65% in sheep (Table 1).

Trait-related eRNAs, regulators and target genes of eRNAs in Animal-eRNAdb

We identified a total of 100 723 trait-related eRNAs in these species (from 99 in worm to 42 358 in mouse). We found a total of 157 500 eRNAs associated with 12 232 TFs (from 505 eRNAs associated with 740 TFs in worm to 58 249 eRNAs associated with 1574 TFs in mouse). We also identified a total 135 501 eRNAs related to 151 673 target genes (from 505 eRNAs related to 13,943 target genes in worm to 48 674 eRNAs related to 22 664 target genes in mouse). More detailed information is presented in Table 2.

Table 2.

Data summary of eRNAs

		Regulators		Target genes		Sequence similarity
Species	No. of trait-related eRNAs	eRNAs	Regulators	eRNAs	Target genes	eRNAs	eRNAs in ≥ 2 species	eRNAs in ≥ 3 species
Gallus gallus (Chicken)	16 337	22 989	956	21 771	20 630	1275	0.34	0.16
Pan troglodytes (Chimp)	1822	5082	1552	4526	8708	11 537	0.95	0.74
Drosophila melanogaster (Fruitfly)	230	1785	632	1634	9358	16	0.94	0.81
Xenopus tropicalis (Clawed frog)	839	1081	649	1074	9665	91	0.41	0.22
Mus musculus (Mouse)	42 358	58 249	1574	48 674	22 664	67 651	0.44	0.24
Rattus norvegicus (Rat)	27 837	31 866	1443	29 968	22 228	26 640	0.27	0.11
Macaca mulatta (Rhesus)	5749	9241	1213	8046	17 153	12 334	0.94	0.69
Ovis aries (Sheep)	3397	5158	1267	4809	13 237	2471	0.47	0.24
Caenorhabditis elegans (Worm)	99	505	740	505	13 943	4	0.25	0.25
Danio rerio (Zebrafish)	2055	21 544	2206	14 494	14 087	2455	0.70	0.68
Summation	100 723	157 500	12 232	135 501	151 673	124 474	0.51	0.31
Minimun	99	505	632	505	8708	4	0.25	0.11
Maximum	42 358	58 249	2206	48 674	22 664	67 651	0.95	0.81
Median	2726	7200	1240	6428	14 015	2463	0.46	0.25

Sequence similarity of eRNAs in Animal-eRNAdb

We identified a total of 124 474 eRNAs with sequence similarities in these species, ranging from four in worm to 67 651 in mouse. Besides, we found that some eRNAs were similar in multiple species. For example, a median of 46% eRNAs have similar eRNAs in more than one species and a median of 25% eRNAs have similar eRNAs in three or more species (Table 2). For the latter, we further counted the number of eRNAs in a given species sharing at least one common related trait with those in other species, and as a result, a total of 17 333 eRNAs were found in different species. For example, chr10:3541913–3547913 in chicken, chr7:54337787–54343787 in rhesus and chr8:60860550–60866550 in rat were all associated with the developmental stage. These eRNAs with sequence similarities will provide a new scope for studying the evolution of species and functions of the eRNAs.

Web interface

Animal-eRNAdb provides a user-friendly interface. Six main modules, including ‘eRNA expression’, ‘trait-related eRNA’, ‘eRNA regulator’, ‘eRNA target gene’, ‘sequence similarity’ and ‘download’ (Figure 2A), are provided for users. Several search/selection boxes are designed on each page, including species selection box, tissue selection box, trait selection box and eRNA search box. In the eRNA search box, users can query a unique eRNA by entering an eRNA ID, or query all eRNAs located in a certain region by entering the genomic region, or query all eRNAs located at ± 1 MB around the start site of the selected gene by entering a gene symbol/Ensembl ID.

Figure 2.

Overview of Animal-eRNAdb. (A) Main functions of Animal-eRNAdb, including the ‘Expression’, ‘Trait’, ‘Regulator’, ‘Target gene’, ‘Sequence’ and ‘Download’ modules. (B) A table of queried eRNAs in the ‘Expression’ module. (C) The expression graph of the queried eRNA chr1:3817931–3823931 across tissues in PRJEB11710 of mouse. (D) Differential expression of the eRNA chr13:15113803–15119803 between female and male in liver in PRJNA471297 of mouse. (E) Significant correlation of the eRNA chr13:15113803–15119803 with the developmental stage of embryo in liver in PRJEB26869 of mouse. (F) Co-expression of the eRNA chr12:112711986–112717986 and putative regulator Zbtb7c (ENSMUSG00000044646) in forebrain in PRJEB26869 of mouse. (G) Co-expression of the eRNA chr15:85681973–85687973 and putative target gene Celsr1 (ENSMUSG00000016028) in testis in PRJEB27404 of mouse. (H) A table of queried eRNAs in the ‘Sequence’ module.

On the ‘eRNA expression’ page, users can query the eRNA expression in specific tissues of a given species. A table with columns of species, bioproject, tissue, eRNA ID, expression and plotAll of the queried eRNAs will be provided (Figure 2B). Users can view the eRNA expression across tissues of a given species and the selected tissue is marked as red in the graph (Figure 2C) by clicking the ‘PlotAll’ button.

On the ‘trait-related eRNA’ page, users can query the eRNAs associated with the developmental stage and gender in specific tissues of a given species. A table with the columns of species, bioproject, tissue, eRNA ID, trait, FC/F/Rho, FDR and plot of the queried eRNAs will be presented. An association diagram will be shown when clicking ‘Plot’ (Figure 2D–E). Users can click the ‘Download’ button to download the queried data or click the ‘?’ button for more information.

On the ‘eRNA regulator’ page, users can search for TFs significantly associated with eRNAs in specific tissues of a given species. Two tables will be shown: one provides the species, tissue, eRNA ID, number of correlated TFs or number of correlated eRNAs and detail, and the other presents the bioproject, tissue, eRNA ID, TF ID, TF symbol, Rho, FDR and plot. Users can view the co-expression correlation diagram between the eRNA and the TF by clicking the ‘Plot’ button (Figure 2F).

On the ‘eRNA target gene’ page, users can search for target genes (within 1 MB) that may be regulated by eRNAs in specific tissues of a given species. A table with the columns of species, bioproject, tissue, eRNA ID, gene ID, gene symbol, gene start, gene end, distance, Rho, FDR and detail of the queried eRNAs will be provided, and the correlation diagram between the eRNA and the target gene will be displayed upon the clicking of the ‘Detail’ button (Figure 2G).

On the ‘sequence similarity’ page, the users can browse eRNAs with sequence similarities in multiple species. A table with columns including species, eRNA ID, seq, match species, match eRNA ID, seq, identify and E-value of the queried eRNAs will be shown (Figure 2H). Users can download the sequence of the eRNA by clicking the ‘Seq’ button and download the list of all eRNAs by clicking the ‘Download’ button.

In Animal-eRNAdb, the users can comprehensively investigate one eRNA through different modules (Supplementary Figure S1). On the ‘Download’ page, users can obtain free access to the main datasets of specific tissues for each species. The ‘Document’ page provides the sample information, reference genome versions, eRNA summary, pipeline of database construction and some other information. Besides, Animal-eRNAdb welcomes any feedback with the email address provided on the ‘Contact’ page.

SUMMARY AND FUTURE DIRECTIONS

Recent advances in experimental techniques and available computing power have led to an exponential growth of biological data of animals besides humans. Using these public resources, many animal-related databases such as AnimalTFDB 3.0 (38), AnimalQTLdb (43) and Animal-imputeDB (44) have been constructed. At present, some progress has been made in the research on human eRNAs, such as TCeA (3) and eRic (16). However, there has been limited research on the mechanisms and functions of eRNAs in other animals. In this study, we developed Animal-eRNAdb by collecting public available data, which provides comprehensive information of eRNAs in different tissues of 10 species. To the best of our knowledge, Animal-eRNAdb is the first and most comprehensive animal eRNA database. In this version of Animal-eRNAdb, by using the data of 5085 samples, we systematically quantified the expression of eRNAs, and identified the trait-related eRNAs, putative eRNA regulators, putative eRNA target genes and eRNAs with sequence similarities across different tissues in various species, which are expected to greatly expand our knowledge of eRNAs in evolution and phenotype. However, considering that only routine RNA-seq data were used and annotations of enhancers were insufficient, the database probably misses many potential eRNAs. In the future, we will further collect available enhancer annotations and sequence data including routine RNA-seq data and nascent RNA-seq data (e.g., GRO-seq and PRO-seq) to characterize eRNAs and update the database. With a comprehensive characterization of eRNAs in various tissues across different species, we believe that Animal-eRNAdb will be a valuable resource for understanding the functions and mechanisms of eRNAs across tissues of multiple species.

DATA AVAILABILITY

Animal-eRNAdb is freely available to the public without registration or login requirements (http://gong_lab.hzau.edu.cn/Animal-eRNAdb/).

Click here for additional data file.