The Functional Map of Ultraconserved Regions in Humans, Mice and Rats

doi:10.21203/rs.3.rs-4837600/v1

Download PDF

Research Article

The Functional Map of Ultraconserved Regions in Humans, Mice and Rats

https://doi.org/10.21203/rs.3.rs-4837600/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

BACKGROUND: Ultraconserved regions (UCRs) encompass 481 DNA segments exceeding 200 base pairs (bp), displaying 100% sequence identity across humans, mice, and rats, indicating profound conservation across taxa and pivotal functional roles in human health and disease. Despite two decades since their discovery, many UCRs remain to be explored owing to incomplete annotation, particularly of newly identified long non-coding RNAs (lncRNAs), and limited data aggregation in large-scale databases. This study offers a comprehensive functional map of 481 UCRs, investigating their genomic and transcriptomic implications: (i) enriching UCR annotation data, including ancestral genomes; (ii) exploring lncRNAs containing T-UCRs across pan-cancers; (iii) elucidating UCR involvement in regulatory elements; and (iv) analyzing population single-nucleotide variations linked to motifs, expression patterns, and diseases.

RESULTS: Our results indicate that, although a high number of protein-coding transcripts with UCRs (1,945 from 2,303), 1,775 contained UCRs outside CDS regions. Focusing on non-coding transcripts, 355 are mapped in 85 lncRNA genes, with 35 of them differentially expressed in at least one TCGA cancer type, seven lncRNAs strongly associated with survival time, and 23 differentially expressed according to single-cell cancer analysis. Additionally, we identified regulatory elements in 373 UCRs (77.5%), and found 353 SNP-UCRs (with at least 1% frequency) with potential regulatory effects, such as motif changes, eQTL potential, and associations with disease/traits. Finally, we identified 4 novel UCRs that had not been previously described.

CONCLUSION: This report compiles and organizes all the above information, providing new insights into the functional mechanisms of UCRs and their potential diagnostic applications.

ultraconserved regions

non-coding RNAs

lncRNAs

ancestral genomes

pan-cancer

variation

Ultraconserved regions (UCRs) are defined as 481 DNA segments longer than 200 bp, with 100% sequence identity in their alignment with humans, mice, and rats [1]. These regions are also highly conserved across disparate taxa, including vertebrates, insects, worms, yeast, and other species, and are essential DNA markers for phylogenomic analyses [1–3]. Approximately 93% of the UCRs were transcribed in at least one human tissue tested, with nearly 79% showing ubiquitous transcription across a wide range of normal human tissues. Because of this characteristic, these regions are called transcribed UCRs (T-UCRs) [4].

The first UCR annotation proposed by Bejerano et al. [1] focused on their overlap with protein-coding genomic regions, with surprisingly only 23.0% of the UCRs overlapping with known mature mRNA sequences [1]. They identified 256 sequences with no evidence of transcription from any matching expressed sequence tag or mRNA, and 111 overlapped the mRNAs of known human protein-coding genes. Additionally, 114 sequences were "inconclusive," where the evidence for transcription was inconclusive. In 2010, Mestdagh and colleagues reorganized all UCR sequences into six categories based on protein-coding information, classifying those regions as intergenic (38.7%), intronic (42.6%), exonic (4.2%), partly exonic (5.0%), and exon-containing (5.6%). Additionally, 3.9% were classified as “multiple” when the genomic annotation varied because of the host gene splice variants [5]. Recently, in 2022, Snetkova and collaborators observed that 371 (77.1%) were non-coding UCRs and 110 (22.9%) were UCRs that overlapped the exons of protein-coding genes. In this study they reviewed the biological importance of these UCRs, mainly in regulating gene expression during embryogenesis [6].

Coding sequences are typically highly conserved and are subject to more stringent negative selection pressures [7]. Nevertheless, more than 75% of the UCRs were in intergenic or intronic, indicating that negative selection pressure can also act on these non-coding sequences[8]. Additionally, among the intronic transcribed UCRs (T-UCRs), almost 58% were detected in the antisense orientation compared with the host gene [4]. These results suggest that many UCRs may be transcribed as long non-coding RNAs (lncRNAs). However, a general annotation of the 481 ultraconserved regions with overlapping lncRNAs was readily unavailable.

In addition to their potential roles in lncRNA function, UCRs are crucial enhancers during development. For example, Visel and colleagues reported that half of the non-exonic UCRs evaluated drove reproducible reporter gene expression in various tissues of developing mouse embryos [9], and this function was reinforced by others [6, 10]. Additionally, highlighting the role of some T-UCRs, deregulated expression has been associated with several pathological conditions, including tumor types and Cancer disease [4, 11–13]. Although the UCRs play an essential role in human health and disease; there is limited information regarding specific T-UCRs, the investigation of annotated lncRNA-T-UCRs largely improves our knowledge of the potential roles of these molecules. Therefore, even two decades after their initial identification, many UCRs remain to be explored. To address this, we constructed a comprehensive functional map of 481 UCRs. however,

In this study, we investigated four significant contributions: (i) we analysed enrichment mining annotation data, including comparisons of the modern genome version and the uncharted ancestral genomes in 15 Neanderthal and ancient humans from different eras; (ii) examination of pan-cancer expression with the support of RNA-Seq and single-cell data analysis, focusing on the T-UCRs in lncRNA genes; (iii) evidence of UCRs in the regulatory elements, including enhancers and microRNA/transcription factor binding sites (TFBS); and (iv) population single-variation analysis, shedding light on their evolutionary dynamics and implications for changes in binding sites, expression, and disease susceptibility. Our research explored the omics landscape of UCRs, shedding light on non-coding regions and deepening our insights into UCR functionality, gene regulation, and disease mechanisms, underscoring their pivotal roles in genomic stability and disease pathogenesis.

2.1. UCR Annotation Enrichment Mining Insights

We analyzed genomic and transcriptomic annotation data using the currently human reference and we compared this to previous UCR annotations, performing the same analysis for recent mouse and rat genomes (Fig. 1, Supplementary Tables S1-S5). Adopting the genomic classification of Bejerano et al. [1], we observed that 126 (26.2%) and 355 (73.8%) UCRs were classified as coding and non-coding, respectively, in the human genome (Fig. 1A). In the mouse and rat genomes, 133 (27.6%) and 348 (72.4%) sequences were classified as coding, respectively, while 82 (17.0%) and 399 (83.0%) were classified as non-coding, respectively (Fig. 1B).

Using the classification of Mestdag et al. [5], for human, we identified 261 (54.5%) UCRs as intronic, four as exonic (0.8%), four as exon-containing (0.8%), five as partly exonic (1.0%), 98 as intergenic (20.4%), and 109 as multiple (22.7%) (Fig. 1C; Supplementary Table S2). In mouse, 223 (44.9%) UCRs were mapped to intronic, 13 (2.7%) to exonic, 10 (2.1%) to exon-containing, 21 (4.3%) to partly exonic, 93 (19.3%) to multiple, and 121 (25.3%) to intergenic regions (Fig. 1C and Supplementary Tables S2 and S4). With respect to the rat genome, our analysis revealed 161 (33.5%) UCRs to an intron, 12 (2.5%) exonic, 14 (2.9%) exon-containing, 16 (3.3%) partly exonic, 28 (5.8%) multiple, and 250 (51.9%) intergenic regions (Fig. 1C; Supplementary Table S2).

To map UCR evidence at transcriptomic levels, we investigate the UCR sequences onto human, mouse and rat transcript annotations. In human, 384 UCRs overlapped with 2,303 transcripts (GENCODE v46). Among these, 1,945 transcripts were mapped to protein-coding T-UCRs (205 protein-coding genes overlapped with 315 T-UCRs), 355 to lncRNA-T-UCR (85 lncRNA genes overlapped with 100 T-UCRs), and three other transcripts were mapped to a pseudogene and two small non-coding RNAs (one miRNA and one miscellaneous RNA (abbrev. miscRNA)), each containing one T-UCR (Fig. 1D). Among the 1,945 protein-coding transcripts annotated, 22 contained UCRs in the coding sequence (CDS) region, 148 overlapped with UCRs in the CDS, and 1,775 contained UCRs outside CDS (UCRs were identified in intronic or exonic non-CDS regions) (Fig. 1E; Supplementary Table S3).

In mouse (GENCODE vM35), 360 UCRs overlapped with 1,125 transcripts, of which 1,022 were protein-coding transcripts (205 protein-coding genes and four protein-coding genes to be experimentally confirmed), 101 were long non-coding transcripts (associated with 46 lncRNA genes), and two were other non-coding RNA (one miRNA and another miscRNA) (Fig. 1D). Of the 1,022 protein-coding transcripts, 14 contained UCRs into the CDS region, 117 overlapped with UCRs in the CDS, and 891 contained UCRs outside the CDS (Fig. 1E; Supplementary Table S4).

In rat (ENSEMBLE 111), we identified 224 T-UCRs overlapping 935 transcripts, with 826 protein-coding transcripts (118 protein-coding genes), 50 lncRNA transcripts (associated with 17 lncRNA genes), 58 miscRNA transcripts (58 genes), and one miRNA (Fig. 1D and Supplementary Table S5). Among the protein-coding transcripts, two contained UCRs in the CDS, 54 overlapped with UCRs in the CDS, and 770 contained UCRs outside the CDS (Fig. 1E). This unique miRNA gene is associated with the uc.420/miR-3064 gene in the three species.

2.2. Conservation of UCRs in Ancient Human Genomes

We also investigated the presence of UCRs in five Homo neanderthalensis genomes (one male and four females) and ten Homo sapiens genomes (nine males and one female) from various periods, including the Palaeolithic era to more recent times, and present-day tribal populations, such as those from the Khoisan Bantu people from Kalahari desert (Fig. 2A) (Supplementary Table S6).

We obtained all of 481 UCRs with 100% coverage and identity in these genomas, although some had regions with sequences containing few bases represented as ‘N’ without the exact nucleotide defined (Supplementary Table S7). Furthermore, both modern and ancient human genomes contained an additional set of novel UCRs with high conservation levels (sequences of at least 200 base pairs, minimum alignment coverage of 75%, and identity threshold of 80% for the UCR set) in distinct chromosomal regions, except in the Mezmaiskaya genome (Fig. 2B and Supplementary Table S8). We found four novel UCRs align perfectly in the human, mouse, and rat genomes and follow Bejerano’s rules for UCR definition (Fig. 2C; Table 1 and Supplementary Table S9), and we denominated them ‘extra’ UCRs (exUCRs as abbreviations and “euc.” for the code sequence). In particular, we observed the presence of euc.397 (chr18) in the analyzed genomes. We identified euc.129 (chr13) and euc.425 (chr16) in all genomes, except for ZKU and Bantu, respectively. Finally, euc.175 (chr10) was present in these genomes, except in the ZKU and Scandinavian genomes(Fig. 2C).

Table 1

ExUCR vs UCR coordinates, and overlapping gennomic annotations for human, mouse and rat.
genome	exUCR	exUCR coordinates	exUCR gene overlapping	exUCR gene overlapping class	UCR	UCR coordinates	UCR gene overlapping	UCR gene overlapping class
human	euc.129	chr13:97356685–97356888	MBNL2 (exonic)	protein-coding	uc.129	chr3:152446598–152446809	MBNL1 (exonic)	protein-coding
mouse	euc.129	chr14:120633013–120633216	MBNL2 (exonic)	protein-coding	uc.129	chr3:60521998–60522209	MBNL1 (exonic)	protein-coding
rat	euc.129	chr15:97508498–97508701	EFCAB12-ARSJ (intergenic)	NA (intergenic)	uc.129	chr2:144798866–144799077	MBNL1 (exonic)	protein-coding
human	euc.175	chr10:129893328–129893555	EBF3 (intronic)	protein-coding	uc.175	chr5:158914830–158915079	EBF1 (intronic)	protein-coding
mouse	euc.175	chr7:136852289 136852516	EBF3 (intronic)	protein-coding	uc.175	chr11:44691125–44691374	EBF1 (intronic)	protein-coding
rat	euc.175	chr1:192052056 192052283	EBF3 (intronic)	protein-coding	uc.175	chr10:29284049–29284298	EBF1 (intronic)	protein-coding
human	euc.397	chr18:25285204–25285497	ZNF521 (intronic)	protein-coding	uc.397	chr16:49701937–49702247	ZNF423 (intronic)	protein-coding
mouse	euc.397	chr18:14029595–14029888	ZFP521 (intronic)	protein-coding	uc.397	chr8:88563296–88563607	ZFP423 (intronic)	protein-coding
rat	euc.397	chr18:4997387–4997680	ZFP521 (intronic)	protein-coding	uc.397	chr19:19231570–19231880	ZFP423 (intronic)	protein-coding
human	euc.425	chr16:49701982–49702237	ZNF423 (intronic)	protein-coding	uc.425	chr18:25285243–25285567	ZNF521 (intronic)	protein-coding
mouse	euc.425	chr8:88563341–88563596	ZNF423 (intronic)	protein-coding	uc.425	chr18:14029634–14029958	ZNF521 (intronic)	protein-coding
rat	euc.425	chr19:19231580–19231835	SYNPOL2L-MYO5B (intergenic)	NA (intergenic)	uc.425	chr18:4997426–4997750	ZNF521 (intronic)	protein-coding

Owing to the lack of comprehensive annotations for the Neanderthal and ancient human genomes, we utilized the human GENCODE v46 reference annotation to describe the exUCRs and compared these annotations with those from mouse and rat. All exUCRs overlapped with protein-coding genes in human and mouse annotations, whereas euc.129 and euc.425 were intergenic in rat (Table 1 and Supplementary Tables S10 and S11). Euc.129 is located between EFCAB12 gene (upstream) and ARSJ gene (downstream), whereas euc.425 is located between the SYNPOL2L gene (upstream) and the MYO5B gene (downstream). Interestingly, all extra UCRs were located in protein-coding genes from the same family from the related UCR, even maintaining overlap within the region (exonic or intronic) of the gene (Table 1). This finding suggests that the ultraconserved sequence is essential for the structure of these proteins or for transcript regulation.

2.3. The T-UCR in LncRNAs: Insights in Pan-Cancer Analyses

There is limited information regarding T-UCRs (Supplementary Table S12). The investigation of annotated lncRNA-T-UCRs improves our knowledge of the potential functional roles in specific pathologic conditions, such as Cancer disease. We performed differential expression analysis of lncRNA-T-UCRs via RNA-seq data (bulk and single-cell) from public cancer databases, including The Cancer Genome Atlas Program (TCGA) and the Curated Cancer Cell Atlas (3CA).

Differential expression of lncRNA-T-UCRs in cancer

We selected datasets from the TCGA database using sample size as the only criterion (at least 30 samples per condition, tumor, and non-tumoral). Ten cancer datasets were obtained: breast invasive carcinoma (BRCA), colon adenocarcinoma (COAD), head and neck squamous cell carcinoma (HNSC), kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), prostate adenocarcinoma (PRAD), thyroid carcinoma (THCA), and uterine corpus endometrial carcinoma (UCEC). The number of non-tumoral and tumoral samples for each cancer is detailed in Supplementary Table S13.

This analysis aimed to detect lncRNAs containing at least one T-UCR transcript isoform. Thus, to identify differentially expressed lncRNA genes, we used the 85 lncRNAs identified in our transcript mapping analysis (i.e., 85 lncRNAs mapped to 100 T-UCRs). However, the RNA-Seq data did not include data for two lncRNAs (ENSG00000288692 and ENSG00000289413). Therefore, we performed tumor and single-cell differential expression analyses of 83 lncRNAs.

Using the Wilcoxon rank-sum test to calculate differentially expressed genes and filtering the results by p-values (≤ 0.05), we detected differentially expressed (DE) lncRNAs in each cancer type (Supplementary Figure S1 and Supplementary Table S14). Considering the threshold of + 2 log2foldChange, the number of genes upregulated in cancers ranges from three (HNSC and PRAD) to 11 (LUSC). With respect to the downregulated lncRNAs (log₂foldchange ≤ -2), we identified only one lncRNA-T-UCR in LUAD, THCA, and UCEC cancer types and four in KIRC. BRCA and PRAD cancers did not present downregulated lncRNA-T-UCR (Supplementary Table S14).

When we compared the expression of lncRNA genes across all tumor types, we identified 35 lncRNAs that were differentially expressed in at least one cancer type (Supplementary Table S15). Using TPM values, we constructed a heatmap using these 35 DE lncRNA-T-UCRs for each cancer type (Supplementary Figure S2). None of the lncRNA-T-UCRs were differentially expressed in any cancer type. However, we detected two genes, FEZF1-AS1/uc.232 and DLX6-AS1/uc.220/uc.221, which were upregulated in the seven cancer types (Fig. 3A). The lncRNA FEZF1-AS1/uc.232 (ENSG00000230316.7) was upregulated in BRCA, COAD, HSNC, LUAD, LUSC, PRAD, and UCEC (fold-change ranging from 2.3 in UCEC to 8,82 in COAD); whereas DLX6-AS1/uc.220/uc.221 was upregulated in BRCA, COAD, HSNC, KIRC, LUAD, LUSC, and UCEC (fold-change ranging from 1.07 in PRAD to 6,7 in LUSC).

We did not identify any downregulated lncRNA-T-UCRs across all cancer types. However, we observed downregulated three lncRNAs in three cancers: MECOM-AS1/uc.136 (HNSC, KIRC, KIRP), MEF2C-AS1/uc.167 (COAD, LUSC, and UCEC), and ENSG00000257277.2/uc.70, an lncRNA without a gene symbol defined and antisense to the coding ARHGAP15 gene (differentially expressed in COAD, LUAD, LUSC). Furthermore, four lncRNAs were up-regulated in one or two cancer types and downregulated in another (SAMMSON/uc.116, HOXB-AS3/uc.415/uc.416/uc.417, LINC01505/uc.266, and ERICH2-DT/uc.95). Finally, regarding lncRNAs expressed exclusively in one type of cancer, 10 were upregulated and seven were downregulated (Fig. 3A).

The impact of lncRNA-T-UCRs on patient prognosis

To analyze the impact of lncRNA-T-UCRs on patient prognosis via TCGA survival data, we investigated the expression of all 83 human lncRNA-T-UCRs in the 10 TCGA cancer types selected, and 19 lncRNAs were associated with survival in the Cox proportional-hazard analysis, seven of which also exhibited significance in the Kaplan-Meier model (Supplementary Figure S3 and Supplementary Table S16).

Focusing on lncRNA-T-UCRs that were significantly associated with survival in both analyses, we identified LINC00583/uc.250, FEZF1-AS1/uc.233, FOXG1-AS1/uc.361, DPH6-DT/uc.382, PKN2-AS1/uc.31, AQP4-AS1/uc.427, and ENSG00000270087.5/uc.286. Interestingly, the lncRNA-UCR analysis revealed a cancer-specific clinical impact, since none were shared among more than one cancer type. This result emphasizes the tissue specificity of these molecules, which is a biologically prominent characteristic of these molecules, highlighting their use as prognostic markers.

lncRNA-T-UCR single-cell RNA-seq (scRNA-seq) data analysis

To identify the 83 lncRNA-T-UCRs expressed within the tumor microenvironment, we obtained publicly available scRNA-seq data from the Curated Cancer Cell Atlas database. Among the lncRNA-T-UCRs, 35 and 23 were differentially expressed according to single-cell analysis (Fig. 3B). We discovered that NR2F1-AS1/uc.169 was the most highly expressed lncRNA in the fibroblasts of all the tumor types analyzed. NR2F1-AS1/uc.169 was also expressed in pericytes from BRCA and COAD patients and in malignant cells from LUAD and COAD patients (Fig. 3B).

Additionally, IQCH-AS1/uc.389 was highly expressed in epithelial cells from HNSC and mast cells from BRCA and LUAD. MEF2C-AS1/uc.167 was highly expressed in B cells from HNSC, BRCA, and COAD. Specifically, in LUAD, FEZF1-AS1/uc.232 was highly expressed in malignant cells, HOXB-AS3/uc.415 to uc.417 was highly expressed in fibroblasts, and AQP4-AS1/uc.427 was highly expressed in epithelial cells. LINC01117/uc.109 was highly expressed in the stromal and malignant cells of BRCA and COAD. LINC01505/uc.266 was highly expressed in NK cells of KIRC, whereas SAMMSON/uc.116 was expressed in the B cells of KIRC. SATB1-AS1/uc.113 was highly expressed in immune cells, and PKN2-AS1/uc.31 was highly expressed in malignant and epithelial cells in HNSC (Fig. 3B).

2.4. UCRs in Regulatory Elements

We performed enrichment mining for the presence of regulatory elements within UCRs via the ORegAnno tool [14]. We identified regulatory elements in 373 UCRs (77.5%), comprising 234 (62.7%) instances as regulatory regions (RR) (including promoters, enhancers, etc.), 73 (19.6%) as transcription factor binding sites (TFBS), 64 (17.1%) as combinations of TFBS and RR, one as a TFBS and miRNA binding site (SMARCA4 and miR-429, respectively), and one as a miRNA binding site (miR-214-3p) in the human genome (Supplementary Table S17). Additionally, 308 UCRs were 100% inserted into RRs and TFBS. Among the 45 distinct TFBSs, SMARCA4 was the binding site most representative of 62 UCRs, followed by EGR1, FOXA1, CTCF, and CEBPB, which contained 20, 13, 12, and 11 UCRs, respectively. Despite SMARCA4, 51 UCRs were completely inserted into the TFBS sequence, and another 11 UCRs partially overlapped. Among the UCRs related to EGR1, FOXA1, CTFC, and CEBPB TFs, 18, 5, 5, and 4 UCRs, respectively, were contained in these TFBSs.

The 19 UCRs exhibited agglomerative numbers of TFBS: uc.213 and uc.246 related to six TFBS; uc.88, uc.329, uc.374, and uc.392 with five TFs; uc.95, uc.162, uc.393, and uc.417 with four TFBS; and uc.74, uc.189, uc.232, uc.277, uc.398, uc.409, uc.412, uc.416, and uc.456 with three TFBS. Finally, in 26 UCRs, we identified two TFBS associated with the same UCR, and a single TFBS in the remaining 93 UCRs.

According to the VISTA Enhancer database, 295 of the 373 (79.1%) regulatory regions within UCRs were identified as human enhancers [15]. For four UCRs (uc.29, uc.347, uc.350 and uc.355) we report the presence of two simultaneous enhancers, totaling 299 enhancers. Among these enhancers, almost 80% (240/299) were experimentally validated and linked with the respective UCR, according to previous studies [9, 10]. In addition, 152 and 147 enhancers were classified as positively and negatively regulating transcription, respectively (Supplementary Table S17).

2.5. Personalized Genomics: Pan-SNP Data Analysis

SNP-UCR compendium identified in human databases and pangenomes

The availability of complete human genomes in various databases allowed us to examine the presence of single nucleotide polymorphisms (SNPs) in UCR sequences and investigate their potential functional effects. In our study, we identified 353 SNPs (with at least 1% frequency) that mapped to 223 UCRs in the Genome Aggregation Database (gnomAD) [16], 1000 Genomes Project (1Kgenomes) and Archive of Brazilian Mutations (ABraOM) [17] (see Supplementary Table S18).

We found that 17 (4.8%) of these SNPs are shared across all three databases (gnomAD, 1kgenomes, and ABraOM), 301 (85.3%) were shared between gnomAD and 1k genomes but were not present in ABraOM, 33 (9.3%) are unique to gnomAD, and two (0.6%) are unique to 1kgenomes (Fig. 4A). The predominant SNP types were single nucleotide variants (SNVs) (332/353), followed by insertion or deletion markers (INDELs) (17/353), deletions (3/353), and insertions (1/353). Most of the variants were located in introns of the coding gene, including 146 (41%) SNPs. Additionally, 13 SNPs were found in exonic-coding regions (of them, six were synonymous and seven were nonsynonymous SNVs), seven in the 3’UTR, three in the 5’UTR of protein-coding genes, and six in exonic and 58 in intronic regions of lncRNA genes. A total of 120 SNPs were identified in the intergenic regions.

Concerning the ethnicity identified in these databases, 36 SNPs were present in African, European, American, East, and South Asian populations, whereas 63, 48, 3, 19, and 20 SNPs were unique to African, European, American, East, and South Asian populations, respectively. In addition, 164 SNPs were distributed among at least two ethnicities, and 17 SNPs were shared with the Brazilian population (Fig. 4A).

UCR-SNPs in regulatory regions associated with diseases

To understand the effects of some SNPs on the regulation of nearby genes and those associated with diseases, we assessed whether UCR-SNPs mapped to a sequence of a transcription factor binding motif (called a motif SNP), which SNPs affects gene expression (called eQTL SNPs), are associated with a complex disease or trait through a Genome-Wide Association Study (GWAS SNPs) [18]. Using the Super Enhancer database (SEdb v2.0), we identified 48 SNPs with potential regulatory effects in 43 UCRs (Supplementary Tables S19 and S20). 21 UCRs containing exclusively motif SNPs; 12 UCRs with eQTL and motif SNPs; 11 with eQTL SNPs; one with motif, GWAS, and eQTL SNPs; one containing eQTL and GWAS SNPs; and two UCRs containing motif and GWAS SNPs (Supplementary Table S19 and S20). For UCRs with GWAS SNPs, we identified that rs12981/uc.268 (RC3H2 gene, 3’UTR variant) was associated with coronary artery calcification [19], rs17291131/uc.211 (SKAP2 gene, intron variant) was associated with type 2 diabetes [20], rs2056117/uc.140 (intergenic variant) was associated with systemic lupus erythematosus [21], and rs11190870/uc.302 (intergenic variant) was associated with scoliosis [22–24].

In addition, we identified eight UCRs with SNPs in more than 1% of the global population that were correlated with disease/trait according to the GWAS catalogue [25] (Supplementary Table S21). rs1861100/uc.53 contains an intron variant in LINC01122 and was related to body mass index (obesity) [26–28]). Additionally, rs142252570/uc.380, an intergenic variant related to the LINC02304/LINC02325 genes, was associated with menarche (age at onset) [29]. The SNP rs142252570/uc.380, an intergenic variant related with the RSRC1 gene, was associated with body mass index and drinks per week [30, 31]. The SNP rs73174306/uc.136, an intron variant in MECOM and MECOM-AS1 genes, was related to blood glucose levels and pancreatic hormone levels [32–36]. The SNP rs1538101/uc.252, an intron variant in the BNC2 gene, was related to velopharyngeal dysfunction [37]. The SNP rs11190870/uc.302, an intergenic variant reported with the LINC01514/LBX1 genes, was related to scoliosis [22, 24, 38, 39]. The SNP rs182041872/uc.315, an intergenic variant reported with the ACADSB/HMX3 genes, was related to antidepressant treatment resistance [33].

Investigation of SNPs in UCR sequences for ancient humans, Neanderthals, and human pangenomes (Fig. 4B and Supplementary Table S22) revealed six altered SNPs with a prevalence of > 50%. Among these SNPs, five were located in intronic regions correlated with GWAS/SEdb with impact disease. Rs2682406/uc.133 (Fig. 4B in the upper triangle) and rs1861100/uc.53 (Fig. 4B in the lower triangle) were related to breast cancer and obesity, respectively [40, 41]. Rs1538101/uc.252, mapped to the protein-coding gene BNC2, was associated with velopharyngeal dysfunction and rs11896224/uc.82, an intron variant of the lncRNA LINC01876 was also associated with velopharyngeal dysfunction. Finally, two SNVs were found to be intergenic: rs11190870/uc.302 (associated with scoliosis) and rs7092999/uc.295 with no disease association.

Some regions related to the 481 UCRs initially described by Bejerano and colleagues in 2004 [1] have been described as important in cellular processes and are associated with diseases [12]. However, despite two decades since the first description, the role of most regions has yet to be discovered. Previous UCR annotations have focused predominantly on DNA coding to provide functional evidence. In this study, we innovatively enriched the annotation of UCR data at both the genomic and transcriptomic levels by incorporating the latest human, mouse, and rat genome versions. Our analysis reinforces the presence of most UCRs in the non-coding and outside CDS regions. Despite this, the transcripts associated with these UCRs are mostly annotated in the intronic regions of protein-coding transcripts, with a quarter annotated as lncRNA genes on the basis of biotype annotation. To the best of our knowledge, this is the first study to focus on UCRs mapped to lncRNAs, providing accessible and organized information such as gene symbols and IDs.

Comparing human, mouse, and rat annotations is an interesting way to observe human and mouse similarities. For example, the percentages of intergenic UCRs in human and mouse is approximately 20 and 25%, respectively; in contrast, 52% of UCRs were found in intergenic positions in rat. This difference was also verified by examining the transcripts that included UCR. The observed differences may result from the limited information available for the rat genome because the rat transcript data are significantly less complete than those available for the mouse, particularly in terms of transcribed and exon regions [42].

With respect to transcript analysis, a unique miRNA gene, MIR3064, was wholly contained in the uc.420 genomic region of the three organisms. MiR-3064 is related to diverse diseases [43–47]. In cancer, miR-3064 is targeted by several lncRNAs/circRNAs associated with important pathways such as PI3K/Akt [48], hTERT [49], glycolysis [50], and Src, and plays crucial tumor suppressive roles [48, 51–59].

Given the extreme conservation of UCRs in mammals, they are expected to be preserved in the ancestral hominids. However, this has not been previously confirmed. Considering the analysis of the 481 UCRs within ancestral genomes, including Neanderthal and ancient human genomes from different eras, as expected, UCR sequences in hominid genomes were also 100% conserved. Notably, we identified four novel UCRs in humans (modern and ancestral genomes), mouse, and rat that had not been previously described. Interestingly, these four novel UCRs (exUCRs) are protein-coding genes from the same family as the UCR related genes, suggesting a potential motif associated with the conserved sequence in these molecules. euc.129 is a novel UCR related to uc.129, which maps to the MBNL1 gene. In contrast, euc.129 is located in the MBNL2 gene, which plays a role in MBNL1. Similarly, the new UCR exuc.175 overlaps the EBF3 gene while uc.175 the EBF1; exuc.397 overlaps the ZNF521 gene while the uc.397 the ZNF423; and the exuc.425 overlaps the ZNF423 gene, while uc.425 the ZNF521. MBNL and ZNF are both zinc finger proteins, and EBF are DNA-binding proteins, all of which act as transcription factors [60].

In human, our approach enabled us to associate 85 lncRNA-T-UCRs with previously unrelated functions, demonstrating its potential for analysis in public databases through Ensemble codes and gene symbols. This is the first study to analyze the lncRNA-T-UCR panel in a pancancer context. These lncRNAs play essential roles in cancer [61, 62]; however, they are not previously recognized as containing T-UCRs. For example, the uc.232 sequence is almost entirely in exon 2 of FEZF1-AS1-202 and in the opposite strand in exon 1 of FEZF1-203 (with no protein associated with the retained intron). A unique study focused on uc.232 showed that the T-UCR was differentially expressed in TCGA gastric cancer samples compared with their non-tumoral counterparts, but no expression of this T-UCR was observed when a different cohort of GC samples was utilized [63]. Here, we reinforced the importance of FEZF1-AS1 in multiple tumor types and highlighted the relationship between this lncRNA and uc.232.

In the survival analysis, we identified seven lncRNA-T-UCRs in cancer cohorts: FEZF1-AS1/uc.233, LINC00583/uc.250, FOXG1-AS1/uc.361, DPH6-DT/uc.382, PKN2-AS1/uc.31, AQP4-AS1/uc.427, and ENSG00000270087.5/uc.286. Two of these genes were previously associated with prognosis in cancer patients, including the widespread deregulation of FEZF1-AS1, which is associated with worse patient prognosis in 13 cancer types [64]. Additionally, PKN2-AS1 was previously identified as being related to the prognosis of patients with bladder cancer [65]. The expression of the other five lncRNAs is associated with survival time for the first time herein.

To the best of our knowledge, this is the first study to analyze T-UCR expression via single-cell data. Data from 35 lncRNA-T-UCRs were obtained, and NR2F1-AS1/uc.169 was the most highly expressed lncRNA in the fibroblasts of all the tumor types analyzed. NR2F1-AS1/uc.169 was also expressed in pericytes from BRCA and COAD patients and in malignant cells from LUAD and COAD patients. The uc.169 sequence overlapped with the coding gene NR2F1 and its antisense transcript NR2F1-AS1. NR2F1 encodes a nuclear hormone receptor that functions as a transcription factor. This protein is critical in the developmental process of the brain[66]. In contrast, antisense transcripts have mostly been studied in the context of cancer [67]. NR2F1-AS1 plays an oncogenic role in several tumor types and is associated with proliferation, migration, and drug resistance [68]. In colorectal and cervical cancers, this lncRNA is downregulated and associated with a tumor-suppressive role[67–69].

Single-cell analysis of NR2F1-AS1 expression in different cell types is underexplored. This lncRNA is usually highly expressed in cancer samples, but in our study, this lncRNA was explicitly found in fibroblast cells for all tumor types analyzed. Cancer-associated fibroblasts (CAFs) are part of the tumor stroma. Although not malignant, fibroblasts may play a role in cancer progression, such as supporting tumor growth, epithelial-mesenchymal transition/metastasis, and therapy resistance [70]. Fibroblasts may also be associated with angiogenesis, as the interaction between tumors and stromal cells can increase vascular endothelial growth factor (VEGF) [71]. Furthermore, NR2F1-AS1/uc.169 was expressed in pericytes in BRCA and COAD. Pericytes are fibroblast-like cells that wrap around endothelial cells in arterioles, capillaries, and venules and are strongly associated with angiogenesis [72]. Our findings underscore the importance of studying lncRNA-T-UCRs in the context of single cells, suggesting more specific functions than previously known.

By organizing data on regulatory elements within the UCRs of the human genome, we found evidence of regulatory elements in 373 UCRs (77.5% of all UCRs). These regulatory elements include enhancers, transcription factors, and miRNA binding sites. The two most abundant TFs in the UCR sequences are SMARCA4 and EGR1. Both TFs are essential for the regulation of critical networks under both physiological and pathological conditions [73, 74], are deregulated during cancer development, and may deregulate the regulation of transcripts related to UCRs.

Finally, we investigated frequent SNPs in UCR sequences, including thousands (~ 128.700) of human whole-genome sequences. Previously, a broad analysis and utilizing ultraconserved regions that include comparisons with other species revealed that some polymorphisms have also been identified, however, the vast majority of these polymorphisms are rare events [75]. In the present study, we focused on 481 UCRs and in the most frequent SNPs (> 1%), identifying 353 SNPs 25 of which were associated with eQTL effects. For example, rs34384113/uc.109 is in the intron portion of LINC01117 with an eQTL effect in HOXD10 and LINC01117 expression; rs117486161/uc.220, in the intron of the DLX6-AS1 gene has an eQTL effect on DLX6-AS1 expression; and rs56805315/uc.417 mapped at the intron of HOXB-AS3 and 5’UTR of the HOXB6 gene has an eQTL effect on AC091133.1, CDK5RAP3, and HOXB6 expression. Interestingly, the lncRNAs LINC01117/uc.109, DLX6-AS1/uc.220, and HOXB-AS3/uc.417 were dysregulated in cancer, as highlighted in our expression analysis.

Our in-depth and comprehensive functional map revealed that almost all UCRs provide tips for understanding their functions. These complete maps of UCRs provide new insights into these regions, accelerating the understanding of functional mechanisms and their potential use as biomarkers of these related regions.

Data sources

UCR data acquisition

The human UCR coordinates were obtained from Bejerano et al. (2004). The sequences of the UCRs were extracted from the human genome via coordinates (obtained from GENCODE release 45 to GRCh38 (hg38).p14)). We used bedtools[76] version 2.31.0 with the getfasta algorithm. General command: bedtools getfasta [OPTIONS] -fi < input FASTA> -bed < BED/GFF/VCF>

The genomic annotation data

The human genomic annotation (GFF/GTF formats) was obtained from GENCODE[77] release 45 to the GRCh38(hg38).p14 genome (available at https://www.gencodegenes.org/human/). For mouse genomic annotation, we used the GENCODE release M34 for the GRcm39/mm39 genome (available at https://www.gencodegenes.org/mouse). Finally, the rat genomic annotation was obtained from RefSeq 106 (ENSEMBLE release 111) for the mRatBN7.2 rn7.2 genome (available at https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_015227675.2).

UCR genome remapping

For human, we remapped the 481 UCRs for NCBI37(hg16) genomic coordinates (BED format) via the liftOver tool (https://genome.ucsc.edu/cgi-bin/hgLiftOver) to human GRCh37(hg19), mouse NCBI34(mm6) and rat Baylor 3.1(rn3). The coordinates were converted from GRCh37(hg19) to GRCh38(hg38).p14 for reference to the human genome, NCBI34(mm6) to GRCm39/mm39 for mouse and Baylor 3.1(rn3) to mRatBN7.2 rn7 for rat via the NCBI-remapping tool (https://www.ncbi.nlm.nih.gov/genome/tools/remap).

UCR overlapping functions

The UCR coordinates (BED format) were overlapped with respective genomic annotation releases in GFF/GTF (GENCODE and RefSeq annotations) or BED formats (mirRNA, regulatory, and SNP/SNVs annotations) and extracted via the R (v4.3.2) script built with the foverlaps function from package data.table (v1.14.2) available at https://www.rdocumentation.org/packages/data.table/versions/1.14.2.

Identification of UCRs in the miRNA annotation data

To detect which UCRs overlapped with miRNAs and mirtrons, we obtained these features from miRBase release 22.1 [78], MirGeneDB 2.1 [79], and MirtronDB 1.1 [80].

Neanderthal and ancient human genomes

Neanderthal genomes were obtained from direct genome projects in BAM format: Altai (https://www.eva.mpg.de/genetics/genome-projects/neandertal/) [81], Chagyrskaya 8 [82] (https://www.eva.mpg.de/genetics/genome-projects/chagyrskaya-neandertal/), Mezmaskaya 3 (PRJNA765125), Vindija 19.33 [83] (http://cdna.eva.mpg.de/neandertal/Vindija/bam/Pruefer_etal_2017/Vindija33.19/), and Denisovan 8 [84] (https://www.eva.mpg.de/genetics/genome-projects/denisova/). The ancient humans in the ENA database: Ust'-Ishim (PRJEB6622), Yana (PRJEB29700), AHUR_2064 (PRJEB29074), Sunghir Burial 3 (SIII) (PRJEB22592), Zlatý kůň (ZKU) (PRJEB39040), Sumidouro 5 (PRJEB29074) Scandinavian hunter-gatherers (PRJEB32786), Ötzi (Tyrolean Iceman) (PRJEB56570), and Ayayema A460 (PRJEB29074). The Kalahari genome was obtained from RAW data (FASTq, paired-end sequences) from the work of Schuster et al. 2010 work [85]. More information about the samples is provided in Supplementary Table S6.

Kalahari genome assembly

All genomic FASTAs were obtained in the BAM format via the consensus algorithm (http://www.htslib.org/doc/samtools-consensus.html) from SAMtools [86] (http://www.htslib.org/). The Kalahari genome was assembled chromosome by chromosome and aligned with GRCh38(hg38) as a reference assembly. The BAM format was also converted to FASTA via the aforementioned method.

Human pangenome project data

From the human pangenome project [87] (https://humanpangenome.org/) we selected the 47 genome haplotypes (primary assemblies) data from the NCBI bioproject (PRJNA730822) in genomic FASTA format, which included the following: 04 Southern Han Chinese (SHC), 07 Puerto Rican in Puerto Rico (PUR), 03 Colombian in Medellin, Colombia (CLM), 02 Esan in Nigeria (ESN), 07 African ancestries from Barbados (Caribbean) (ACB), 04 Peruvian in Lima, Peru (PEL), 01 Vietnamese Kinh In Ho Chi Minh City, Vietnam (KHV), 08 Mandinka in Gambia - Western Division (MWG), 04 Mende in Sierra Leone (MSL), 01 Punjabi in Lahore, Pakistan (PLP), 02 Yoruban in Ibadan, Nigeria (YIN), 02 Kinyawa, Kenya (KWK), 02 Ashkenazim Son HG002 (AIS), 01 Han Chinese Son/HG005 (HCS). More information about the samples is available in Supplementary Table S6.

Identification of ‘extra’ UCRs (exUCRs)

For this analysis, we used the Neanderthal and human ancient genomes previously described and the pangenome genome and human genome references GRCh38(hg38).p14 and T2T-CHM13(hs1) in genomic FASTA format (T2T-CHM13(hs1) reference has complete the sequence of the Y chromosome[88] and was included to mapping the UCRs). First, each genome was converted into the BLAST database format via makeblastdb -in genome.fasta -dbtype nucl -out genomedb.DB. We subsequently performed similarity searches via the BLASTn algorithm of the 481 UCRs against the genomes via the following command line: blastn -query UCRs.fasta -db genome.DB -out ucr-genomeDB.txt -evalue 0.000001 -max_target_seqs 200 -outfmt 6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore qlen slen sstrand qcovs qcovhsp.

We then filtered the results on the basis of size of the aligned region (at least 200 bp), percent identity (at least 75%), and query coverage (at least 80%). Finally, we consider those results located in different coordinates in relation to the reference sequence. That is, a region was considered an “extra UCR” if its coordinates were different from the original UCR in the reference genome.

Literature review of UCRs and UCR overlapped genes

The electronic search was performed up to June 2024 in the PubMed and Google Scholar databases. We have included original studies, reviews, and chapters written in English. In the PubMed database, the terms “uc.x” or “gene name” were searched, and the title and abstracts were screened. In Google Scholar, the terms “ultraconserved region” + “uc.x” were used. Duplicates, theses, preprints, and unrelated texts were excluded from the study.

Cancer Expression Data Sources from TCGA

For gene expression analysis, we obtained RNA-Seq datasets for different cancer types from the TCGA database. We selected cancer datasets containing at least 30 samples per condition (tumor and normal tissue). To obtain the gene read counts and TPM values for each cancer type, we used the TCGAbiolinks package from R [89]. For the differential expression analyses, we selected genes annotated as lncRNAs that contained at least one UCR in some of its isoforms.

Differentialexpression (DE) analysis and visualization of results

For differential expression analysis, we used the Wilcoxon rank-sum test, as described in [90]. For this analysis, we used the gene read counts of the 83 lncRNA-T-UCRs previously selected. As a first step, we filtered genes that had very low counts via the filterByExpr function from the egdeR package. After that, we normalized the gene counts via the trimmed mean of M values (TMM) method. We subsequently calculated each gene’s counts-per-million (CPM) values and used it as input in the wilcox.test function in R to calculate the p-value. Finally, we set a p-value cut-off on the basis of an FDR threshold using the Benjamini & Hochberg method. To visualize the differentially expressed lncRNA genes across the 10 cancer types, we constructed Vulcan plots for each of them via the ggplot2 package in R. In addition, we selected genes that were differentially expressed in at least one type of cancer. Using this list, we constructed heatmaps for each cancer type using the TPM values obtained from TCGA database. The heatmaps were constructed via the function Heatmap of the ComplexHeatmap package from R. To analyze the differentially expressed lncRNA-T-UCRs among the 10 cancer types, we constructed a bar plot via package of ggplot2 from R [91].

Cancer survival analysis from TCGA impact measurement

To analyze the impact of lncRNA-T-UCRs on patients prognosis via TCGA survival data, we employed both the Cox proportional-hazards model and the Kaplan-Meier (KM) method via the survival R package [92]. TCGA survival data were downloaded via XenaBrowser (https://xenabrowser.net/datapages/) and only patients used for differential expression analysis were included. For the Cox model, Bonferroni correction was applied post-hoc for multiple testing adjustments. When a significant result was found with the Cox model, the same lncRNA-T-UCR was further analyzed via the KM model, which groups patients on the basis of median expression value (high and low). Graphical visualizations were generated via the survminer R package (https://cran.r-hub.io/web/packages/survminer/index.html).

Single-cell RNA-sequencing (scRNA-seq) analysis

We analyzed publicly available scRNA-seq data to determine the transcriptional profiles of 83 lncRNA-T-UCRs mapped to UTR and annotated in TCGA. BRCA (GSE148673, GSE161529, and E-MTAB-8107), KIRC (EGAS00001002325, GSE159115, phs002065.v1.p1, nc9bc8dn4m.1), LUAD (GSE131907, GSE123904, HRA000154, and E-MTAB-6149), PRAD (GSE141445), HNSC (nasopharyngeal cancer; GSE150430), and COAD (EGAS00001003779) scRNA-seq data were downloaded from the Curated Cancer Cell Atlas [93] (https://www.weizmann.ac.il/sites/3CA/) database. We selected 162 untreated primary samples belonging to the 10X platform whose count data were available. The selected datasets and samples are described in Supplementary Table S23. THCA and UCEC were excluded from this analysis because of the lack of data availability and our inclusion criteria. The count data were processed via the Seurat (v4.0.6) package in R (v4.2.2) [94]. KIRC, BRCA, and LUAD had more than one available dataset and were integrated via the IntegrateData function [95]. The identities of the cell types were based on the expression of canonical markers. In total, the RNA sequencing data from 596,400 single cells were analyzed. The heatmaps were generated and analyzed via the Morpheus online tool (https://software.broadinstitute.org/morpheus/)

Regulatory element annotation and enrichment sources

To identify possible regulatory evidence in UCRs, we obtained transcription factors (TFs) and enhancers from the human genome (hg38) from the Open Regulatory Annotation database (ORegAnnoDB) version 3.0 [14]. These datasets are also accessible through the UCSC Genome Browser https://genome.ucsc.edu/cgi-bin/hgTrackUi?hgsid=686342163_2it3aVMQVoXWn0wuCjkNOVX39wxy&c=chr1&g=oreganno). The Enhancers dataset is enriched with the VISTA Enhancer [15], Visel et al. (2008) and Snetkova et al. (2021) publications [9, 10].

SNP/SNV data acquisition

SNP identification and annotation were performed via dbSNP version 156 for GRCh38 VCF as a reference [96]. Tabix (HTSlib v1.13) [97] extracted variants in UCRs, resulting in a filtered VCF file. To add variant frequency information, the VCF file was annotated via Annovar (v2020-06-08) [98] with the following formats and versions provided by the tool: ABraOM (v20181204) [99], 1,000 Genomes (v20150824) [100], and gnomAD 4.0 (v20231127) [16]. Variants were filtered to retain those with a frequency equal to or greater than 1% in at least one subpopulation from the analyzed databases (ABraOM; 1,000 Genomes: African, Admixed American, East Asian, European, South Asian; gnomAD: African, Amish, Latino/Admixed American, Ashkenazi Jewish, East Asian, Finnish, Non-Finnish European, Middle Eastern, South Asian, Other, XX, XY). Additionally, to find SNPs in UCRs in ancient humans and pangenome data we aligned all the UCRs with the MUSCLE algorithm [101]. Here, we selected the common SNPs from the Super Enhancer database (SEdb) release v2.0 [18] (available at https://bio.liclab.net/sedb/download.php) to search for SNP motifs, SNP GWASs, and SNP eQTLs.

Author Contribution

B.T.L.N. wrote the main manuscript text and prepared all figures and tables. L.S.O. contributed to drafting the manuscript and participated in the data analysis, data curation and validation of RNA-seq analyses. A.C.R. and J.C.O. contributed to the literature review and manuscript editing. C.M. assisted in data analysis and interpretation of survival patient prognostic. D.F.G. provided critical revisions and contributed to the drafting of the manuscript. A.H.U. helped with methodology development and data curation of SNVs and SNPs. F.P. provided oversight throughout the research project. V.L.S.C., S.S.C. and A.P.S. participated in the data analysis, data curation and validation of scRNA-seq analyses. G.A.C., R.F.C., A.R.P. and J.C.O. provided expert advice and critical feedback on the study's findings. A.R.P. and J.C.O. supervised the project, reviewed the manuscript and contributed to designing the study. All authors reviewed and approved the final manuscript.

Acknowledgement

The results published here are in whole or in part based upon data generated by the TCGA Research Network: https://www.cancer.gov/tcga. This study is supported by Araucaria Foundation - Fundação Araucária in NAPI Bioinformática (# PDI 66/2021) and CNPq (# 440412/2022-6).

Bejerano G, Pheasant M, Makunin I, Stephen S, Kent WJ, Mattick JS, Haussler D. Ultraconserved elements in the human genome. Science. 2004;304:1321–5.
Stephen S, Pheasant M, Makunin IV, Mattick JS. Large-scale appearance of ultraconserved elements in tetrapod genomes and slowdown of the molecular clock. Mol Biol Evol. 2008;25:402–8.
Carter JK, Kimball RT, Funk ER, Kane NC, Schield DR, Spellman GM, Safran RJ. Estimating phylogenies from genomes: A beginners review of commonly used genomic data in vertebrate phylogenomics. J Hered. 2023;114:1–13.
Calin GA, Liu C-G, Ferracin M, et al. Ultraconserved regions encoding ncRNAs are altered in human leukemias and carcinomas. Cancer Cell. 2007;12:215–29.
Mestdagh P, Fredlund E, Pattyn F, et al. An integrative genomics screen uncovers ncRNA T-UCR functions in neuroblastoma tumours. Oncogene. 2010;29:3583–92.
Snetkova V, Pennacchio LA, Visel A, Dickel DE. Perfect and imperfect views of ultraconserved sequences. Nat Rev Genet. 2022;23:182–94.
Prabh N, Rödelsperger C. Are orphan genes protein-coding, prediction artifacts, or non-coding RNAs? BMC Bioinformatics. 2016;17:226.
de Oliveira JC. Transcribed Ultraconserved Regions: New regulators in cancer signaling and potential biomarkers. Genet Mol Biol. 2023;46:e20220125.
Visel A, Prabhakar S, Akiyama JA, Shoukry M, Lewis KD, Holt A, Plajzer-Frick I, Afzal V, Rubin EM, Pennacchio LA. Ultraconservation identifies a small subset of extremely constrained developmental enhancers. Nat Genet. 2008;40:158–60.
Snetkova V, Ypsilanti AR, Akiyama JA, et al. Ultraconserved enhancer function does not require perfect sequence conservation. Nat Genet. 2021;53:521–8.
Fabris L, Calin GA. (2017) Chapter Four - Understanding the Genomic Ultraconservations: T-UCRs and Cancer. In: Galluzzi L, Vitale I, editors International Review of Cell and Molecular Biology. Academic Press, pp 159–172.
Pereira Zambalde E, Mathias C, Rodrigues AC, de Souza Fonseca Ribeiro EM, Fiori Gradia D, Calin GA, de Carvalho J. Highlighting transcribed ultraconserved regions in human diseases. Wiley Interdiscip Rev RNA. 2020;11:e1567.
Gibert MK Jr, Sarkar A, Chagari B, Cells et al. https://doi.org/10.3390/cells11101684
Lesurf R, Cotto KC, Wang G, Griffith M, Kasaian K, Jones SJM, Montgomery SB, Griffith OL, Open Regulatory Annotation Consortium. ORegAnno 3.0: a community-driven resource for curated regulatory annotation. Nucleic Acids Res. 2016;44:D126–32.
Visel A, Minovitsky S, Dubchak I, Pennacchio LA. VISTA Enhancer Browser–a database of tissue-specific human enhancers. Nucleic Acids Res. 2007;35:D88–92.
Chen S, Francioli LC, Goodrich JK, et al. A genomic mutational constraint map using variation in 76,156 human genomes. Nature. 2024;625:92–100.
Naslavsky MS, Scliar MO, Yamamoto GL, et al. Whole-genome sequencing of 1,171 elderly admixed individuals from São Paulo, Brazil. Nat Commun. 2022;13:1004.
Wang Y, Song C, Zhao J, et al. SEdb 2.0: a comprehensive super-enhancer database of human and mouse. Nucleic Acids Res. 2023;51:D280–90.
Ferguson JF, Matthews GJ, Townsend RR, et al. Candidate gene association study of coronary artery calcification in chronic kidney disease: findings from the CRIC study (Chronic Renal Insufficiency Cohort). J Am Coll Cardiol. 2013;62:789–98.
Diabetes Genetics Initiative of Broad Institute of Harvard and MIT, Lund University, and Novartis Institutes of BioMedical Research, Saxena R, Voight BF et al. (2007) Genome-wide association analysis identifies loci for type 2 diabetes and triglyceride levels. Science 316:1331–1336.
Chung SA, Brown EE, Williams AH, et al. Lupus nephritis susceptibility loci in women with systemic lupus erythematosus. J Am Soc Nephrol. 2014;25:2859–70.
Takahashi Y, Kou I, Takahashi A, et al. A genome-wide association study identifies common variants near LBX1 associated with adolescent idiopathic scoliosis. Nat Genet. 2011;43:1237–40.
Miyake A, Kou I, Takahashi Y, et al. Identification of a susceptibility locus for severe adolescent idiopathic scoliosis on chromosome 17q24.3. PLoS ONE. 2013;8:e72802.
Ogura Y, Kou I, Miura S, et al. A Functional SNP in BNC2 Is Associated with Adolescent Idiopathic Scoliosis. Am J Hum Genet. 2015;97:337–42.
Sollis E, Mosaku A, Abid A, et al. The NHGRI-EBI GWAS Catalog: knowledgebase and deposition resource. Nucleic Acids Res. 2023;51:D977–85.
Zhu Z, Guo Y, Shi H, et al. Shared genetic and experimental links between obesity-related traits and asthma subtypes in UK Biobank. J Allergy Clin Immunol. 2020;145:537–49.
Koskeridis F, Evangelou E, Said S, Boyle JJ, Elliott P, Dehghan A, Tzoulaki I. Pleiotropic genetic architecture and novel loci for C-reactive protein levels. Nat Commun. 2022;13:6939.
Huang J, Huffman JE, Huang Y, et al. Genomics and phenomics of body mass index reveals a complex disease network. Nat Commun. 2022;13:7973.
Horikoshi M, Day FR, Akiyama M, et al. Elucidating the genetic architecture of reproductive ageing in the Japanese population. Nat Commun. 2018;9:1977.
Pulit SL, Stoneman C, Morris AP, et al. Meta-analysis of genome-wide association studies for body fat distribution in 694 649 individuals of European ancestry. Hum Mol Genet. 2019;28:166–74.
Saunders GRB, Wang X, Chen F, et al. Genetic diversity fuels gene discovery for tobacco and alcohol use. Nature. 2022;612:720–4.
Sinnott-Armstrong N, Tanigawa Y, Amar D, et al. Genetics of 35 blood and urine biomarkers in the UK Biobank. Nat Genet. 2021;53:185–94.
Sakaue S, Kanai M, Tanigawa Y, et al. A cross-population atlas of genetic associations for 220 human phenotypes. Nat Genet. 2021;53:1415–24.
Sun BB, Maranville JC, Peters JE, et al. Genomic atlas of the human plasma proteome. Nature. 2018;558:73–9.
Lagou V, Jiang L, Ulrich A, et al. GWAS of random glucose in 476,326 individuals provide insights into diabetes pathophysiology, complications and treatment stratification. Nat Genet. 2023;55:1448–61.
Pietzner M, Wheeler E, Carrasco-Zanini J, et al. Mapping the proteo-genomic convergence of human diseases. Science. 2021;374:eabj1541.
Chernus J, Roosenboom J, Ford M, et al. GWAS reveals loci associated with velopharyngeal dysfunction. Sci Rep. 2018;8:8470.
Khanshour AM, Kou I, Fan Y, et al. Genome-wide meta-analysis and replication studies in multiple ethnicities identify novel adolescent idiopathic scoliosis susceptibility loci. Hum Mol Genet. 2018;27:3986–98.
Kou I, Otomo N, Takeda K, et al. Genome-wide association study identifies 14 previously unreported susceptibility loci for adolescent idiopathic scoliosis in Japanese. Nat Commun. 2019;10:3685.
Shen H, Lu C, Jiang Y, et al. Genetic variants in ultraconserved elements and risk of breast cancer in Chinese population. Breast Cancer Res Treat. 2011;128:855–61.
Yang R, Frank B, Hemminki K, et al. SNPs in ultraconserved elements and familial breast cancer risk. Carcinogenesis. 2008;29:351–5.
Ji X, Li P, Fuscoe JC, et al. A comprehensive rat transcriptome built from large scale RNA-seq-based annotation. Nucleic Acids Res. 2020;48:8320–31.
Xu C, Wang Z, Liu YJ, Duan K, Guan J. Harnessing GMNP-loaded BMSC-derived EVs to target miR-3064-5p via MEG3 overexpression: Implications for diabetic osteoporosis therapy in rats. Cell Signal. 2024;118:111055.
Yang W, Tu H, Tang K, Huang H, Ou S, Wu J. MiR-3064 in Epicardial Adipose-Derived Exosomes Targets Neuronatin to Regulate Adipogenic Differentiation of Epicardial Adipose Stem Cells. Front Cardiovasc Med. 2021;8:709079.
Huang M, Li X, Li G. Mesenchyme homeobox 1 mediated-promotion of osteoblastic differentiation is negatively regulated by mir-3064-5p. Differentiation. 2021;120:19–27.
Grosu Ș-A, Dobre M, Milanesi E, Hinescu ME. (2023) Blood-Based MicroRNAs in Psychotic Disorders-A Systematic Review. Biomedicines. https://doi.org/10.3390/biomedicines11092536
Patel RB, Bajpai AK, Thirumurugan K. Differential Expression of MicroRNAs and Predicted Drug Target in Amyotrophic Lateral Sclerosis. J Mol Neurosci. 2023;73:375–90.
Luo Z, Hao S, Yuan J, Zhu K, Liu S, Zhang J, Yao L. Long non-coding RNA LINC00958 promotes colorectal cancer progression by enhancing the expression of LEM domain containing 1 via microRNA miR-3064-5p. Bioengineered. 2021;12:8100–15.
Bai L, Wang H, Wang A-H, Zhang L-Y, Bai J. MicroRNA-532 and microRNA-3064 inhibit cell proliferation and invasion by acting as direct regulators of human telomerase reverse transcriptase in ovarian cancer. PLoS ONE. 2017;12:e0173912.
He R, Zhang FH, Shen N. LncRNA FEZF1-AS1 enhances epithelial-mesenchymal transition (EMT) through suppressing E-cadherin and regulating WNT pathway in non-small cell lung cancer (NSCLC). Biomed Pharmacother. 2017;95:331–8.
Zhang P, Ha M, Li L, Huang X, Liu C. MicroRNA-3064-5p sponged by MALAT1 suppresses angiogenesis in human hepatocellular carcinoma by targeting the FOXA1/CD24/Src pathway. FASEB J. 2020;34:66–81.
Xiao E, Zhang D, Zhan W, Yin H, Ma L, Wei J, Kang Y, Mao Z. circNFIX facilitates hepatocellular carcinoma progression by targeting miR-3064-5p/HMGA2 to enhance glutaminolysis. Am J Transl Res. 2021;13:8697–710.
Wei M, Chen Y, Du W. LncRNA LINC00858 enhances cervical cancer cell growth through miR-3064-5p/ VMA21 axis. Cancer Biomark. 2021;32:479–89.
Shih C-H, Chuang L-L, Tsai M-H, Chen L-H, Chuang EY, Lu T-P, Lai L-C. Hypoxia-Induced MALAT1 Promotes the Proliferation and Migration of Breast Cancer Cells by Sponging MiR-3064-5p. Front Oncol. 2021;11:658151.
Wang S, Ping M, Song B, Guo Y, Li Y, Jia J. Exosomal CircPRRX1 Enhances Doxorubicin Resistance in Gastric Cancer by Regulating MiR-3064-5p/PTPN14 Signaling. Yonsei Med J. 2020;61:750–61.
Yan J, Jia Y, Chen H, Chen W, Zhou X. Long non-coding RNA PXN-AS1 suppresses pancreatic cancer progression by acting as a competing endogenous RNA of miR-3064 to upregulate PIP4K2B expression. J Exp Clin Cancer Res. 2019;38:390.
Khalilian S, Mohajer Z, Khazeei Tabari MA, Ghobadinezhad F, Ghafouri-Fard S. circGFRA1: A circular RNA with important roles in human carcinogenesis. Pathol Res Pract. 2023;248:154588.
Meng M, Wu Y-C. (2022) LMX1B Activated Circular RNA GFRA1 Modulates the Tumorigenic Properties and Immune Escape of Prostate Cancer. J Immunol Res 2022:7375879.
Ji X, Lv C, Huang J, Dong W, Sun W, Zhang H. ALKBH5-induced circular RNA NRIP1 promotes glycolysis in thyroid cancer cells by targeting PKM2. Cancer Sci. 2023;114:2318–34.
Brown GR, Hem V, Katz KS, et al. Gene: a gene-centered information resource at NCBI. Nucleic Acids Res. 2015;43:D36–42.
Shi C, Sun L, Song Y. FEZF1-AS1: a novel vital oncogenic lncRNA in multiple human malignancies. Biosci Rep. 2019. https://doi.org/10.1042/BSR20191202.
Hu C, Liu K, Wang B, Xu W, Lin Y, Yuan C. DLX6-AS1: An Indispensable Cancer-related Long Non-coding RNA. Curr Pharm Des. 2021;27:1211–8.
Khalafiyan A, Emadi-Baygi M, Wolfien M, Salehzadeh-Yazdi A, Nikpour P. Construction of a three-component regulatory network of transcribed ultraconserved regions for the identification of prognostic biomarkers in gastric cancer. J Cell Biochem. 2023;124:396–408.
Zhou Y, Xu S, Xia H, Gao Z, Huang R, Tang E, Jiang X. Long noncoding RNA FEZF1-AS1 in human cancers. Clin Chim Acta. 2019;497:20–6.
Yang FL, Hong K, Zhao GJ, Liu C, Song YM, Ma LL. [Construction of prognostic model and identification of prognostic biomarkers based on the expression of long non-coding RNA in bladder cancer via bioinformatics]. Beijing Da Xue Xue Bao. 2019;51:615–22.
Tocco C, Bertacchi M, Studer M. Structural and Functional Aspects of the Neurodevelopmental Gene NR2F1: From Animal Models to Human Pathology. Front Mol Neurosci. 2021;14:767965.
Hu J, Peng F, Qiu X, Yang J, Li J, Shen C, Yuan C. NR2F1-AS1: A Functional Long Noncoding RNA in Tumorigenesis. Curr Med Chem. 2023;30:4266–76.
Ghafouri-Fard S, Khoshbakht T, Hussen BM, Baniahmad A, Taheri M, Samsami M. A review on the role of NR2F1-AS1 in the development of cancer. Pathol Res Pract. 2022;240:154210.
Luo D, Liu Y, Yuan S, Bi X, Yang Y, Zhu H, Li Z, Ji L, Yu X. The emerging role of NR2F1-AS1 in the tumorigenesis and progression of human cancer. Pathol Res Pract. 2022;235:153938.
Tao L, Huang G, Song H, Chen Y, Chen L. Cancer associated fibroblasts: An essential role in the tumor microenvironment. Oncol Lett. 2017;14:2611–20.
Gomes FG, Nedel F, Alves AM, Nör JE, Tarquinio SBC. Tumor angiogenesis and lymphangiogenesis: tumor/endothelial crosstalk and cellular/microenvironmental signaling mechanisms. Life Sci. 2013;92:101–7.
Jiang Z, Zhou J, Li L, Liao S, He J, Zhou S, Zhou Y. Pericytes in the tumor microenvironment. Cancer Lett. 2023;556:216074.
Dreier MR, Walia J, de la Serna IL. (2024) Targeting SWI/SNF Complexes in Cancer: Pharmacological Approaches and Implications. Epigenomes. https://doi.org/10.3390/epigenomes8010007
Wang B, Guo H, Yu H, Chen Y, Xu H, Zhao G. The Role of the Transcription Factor EGR1 in Cancer. Front Oncol. 2021;11:642547.
Habic A, Mattick JS, Calin GA, Krese R, Konc J, Kunej T. Genetic Variations of Ultraconserved Elements in the Human Genome. OMICS. 2019;23:549–59.
Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.
Frankish A, Carbonell-Sala S, Diekhans M, et al. GENCODE: reference annotation for the human and mouse genomes in 2023. Nucleic Acids Res. 2023;51:D942–9.
Kozomara A, Birgaoanu M, Griffiths-Jones S. miRBase: from microRNA sequences to function. Nucleic Acids Res. 2019;47:D155–62.
Fromm B, Høye E, Domanska D, et al. MirGeneDB 2.1: toward a complete sampling of all major animal phyla. Nucleic Acids Res. 2022;50:D204–10.
Da Fonseca BHR, Domingues DS, Paschoal AR. mirtronDB: a mirtron knowledge base. Bioinformatics. 2019;35:3873–4.
Prüfer K, Racimo F, Patterson N, et al. The complete genome sequence of a Neanderthal from the Altai Mountains. Nature. 2014;505:43–9.
Mafessoni F, Grote S, de Filippo C, et al. A high-coverage Neandertal genome from Chagyrskaya Cave. Proc Natl Acad Sci U S A. 2020;117:15132–6.
Prüfer K, de Filippo C, Grote S, et al. A high-coverage Neandertal genome from Vindija Cave in Croatia. Science. 2017;358:655–8.
Meyer M, Kircher M, Gansauge M-T, et al. A high-coverage genome sequence from an archaic Denisovan individual. Science. 2012;338:222–6.
Schuster SC, Miller W, Ratan A, et al. Complete Khoisan and Bantu genomes from southern Africa. Nature. 2010;463:943–7.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R, 1000 Genome Project Data Processing Subgroup. The Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009;25:2078–9.
Liao W-W, Asri M, Ebler J, et al. A draft human pangenome reference. Nature. 2023;617:312–24.
Rhie A, Nurk S, Cechova M, et al. The complete sequence of a human Y chromosome. Nature. 2023;621:344–54.
Colaprico A, Silva TC, Olsen C, et al. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res. 2016;44:e71.
Li Y, Ge X, Peng F, Li W, Li JJ. Exaggerated false positives by popular differential expression methods when analyzing human population samples. Genome Biol. 2022;23:79.
Wickham H. ggplot2: Elegant Graphics for Data Analysis. Springer; 2016.
Borgan Ø, Patricia M, grambsch. Springer-Verlag, New York, 2000. No. Of pages: Xiii + 350. Price: $69.95. ISBN 0‐387‐98784‐3. Stat Med 20:2053–2054.
Gavish A, Tyler M, Greenwald AC, et al. Hallmarks of transcriptional intratumour heterogeneity across a thousand tumours. Nature. 2023;618:598–606.
Satija R, Farrell JA, Gennert D, Schier AF, Regev A. Spatial reconstruction of single-cell gene expression data. Nat Biotechnol. 2015;33:495–502.
Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM 3rd, Hao Y, Stoeckius M, Smibert P, Satija R. Comprehensive Integration of Single-Cell Data. Cell. 2019;177:1888–e190221.
Sherry ST, Ward MH, Kholodov M, Baker J, Phan L, Smigielski EM, Sirotkin K. dbSNP: the NCBI database of genetic variation. Nucleic Acids Res. 2001;29:308–11.
Bonfield JK, Marshall J, Danecek P, Li H, Ohan V, Whitwham A, Keane T, Davies RM. (2021) HTSlib: C library for reading/writing high-throughput sequencing data. Gigascience. https://doi.org/10.1093/gigascience/giab007
Wang K, Li M, Hakonarson H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res. 2010;38:e164.
Naslavsky MS, Yamamoto GL, de Almeida TF, et al. Exomic variants of an elderly cohort of Brazilians in the ABraOM database. Hum Mutat. 2017;38:751–63.
1000 Genomes Project Consortium, Auton A, Brooks LD, et al. A global reference for human genetic variation. Nature. 2015;526:68–74.
Edgar RC. MUSCLE: a multiple sequence alignment method with reduced time and space complexity. BMC Bioinformatics. 2004;5:113.

No competing interests reported.

Download PDF

Version 1

posted

You are reading this latest preprint version

The Functional Map of Ultraconserved Regions in Humans, Mice and Rats

Status:

Version 1

Abstract

Figures

1. INTRODUCTION

2. RESULTS

2.1. UCR Annotation Enrichment Mining Insights

2.2. Conservation of UCRs in Ancient Human Genomes

2.3. The T-UCR in LncRNAs: Insights in Pan-Cancer Analyses

2.4. UCRs in Regulatory Elements

2.5. Personalized Genomics: Pan-SNP Data Analysis

3. DISCUSSION

4. METHODS

Declarations

Author Contribution

Acknowledgement

References

Additional Declarations

Supplementary Files

Status:

Version 1