Multi-omics analysis in human retina uncovers ultraconserved cis -regulatory elements at rare eye disease loci

  1. Victor Lopez Soriano
  2. Alfredo Dueñas Rey
  3. Rajarshi Mukherjee
  4. Genomics England Research Consortium
  5. Frauke Coppieters
  6. Miriam Bauwens
  7. Andy Willaert
  8. Elfride De Baere

Reviewed by

Read the full article See related articles

Listed in

Log in to save this article

Abstract

Cross-species genome comparisons have revealed a substantial number of ultraconserved non-coding elements (UCNEs). Several of these elements have proved to be essential tissue- and cell type-specific cis -regulators of developmental gene expression. Here, we characterized a set of UCNEs as candidate CREs (cCREs) during retinal development and evaluated the contribution of their genomic variation to rare eye diseases, for which pathogenic non-coding variants are emerging. Integration of bulk and single-cell retinal multi-omics data revealed 594 genes under potential cis -regulatory control of UCNEs, of which 45 are implicated in rare eye disease. Mining of candidate cis -regulatory UCNEs in WGS data derived from the rare eye disease cohort of Genomics England revealed 178 ultrarare variants within 84 UCNEs associated with 29 disease genes. Overall, we provide a comprehensive annotation of ultraconserved non-coding regions acting as cCREs during retinal development which can be targets of non-coding variation underlying rare eye diseases.

Article activity feed

  1. Version published to 10.1101/2023.06.08.23291052v2 on medRxiv
  2. Review Commons

    Note: This rebuttal was posted by the corresponding author to Review Commons. Content has not been altered except for formatting.

    Learn more at Review Commons

    Reply to the reviewers

    We would like to thank all reviewers for their highly valuable comments, reviewing the article, and suggesting changes to improve its overall structure, clarity, and comprehension. Please find below our point-by-point responses to each reviewer’s comments. Lines, figures, and tables we refer to in the responses correspond to the clean copy of the revised manuscript.

    Reviewer #1:

    1. In the introduction the authors list 4 databases of ultra-conserved non-coding elements, and choose to work with the UCNEbase that detects UCNE regions by comparing human and chicken, as they state that it is one of the most comprehensive resources. Can the authors comment on the extent of overlap of UCNE regions identified in this database compared to other databases. It would also be helpful to add an explanation in the Introduction on the advantage of comparing the human genome to the chicken genome, e.g., is chicken the most distant vertebrate with a high quality genome, or another reason? And can the authors comment on the extent of conservation of the UCNEs compared to other species - in the UCNEbase paper, orthologues for the UCNEs are identified in 18 vertebrate species including reptiles, amphibians and fish.

    We have addressed these comments and incorporated them in the introduction (lines 46-47; 49-52). As stated, it is not straightforward to evaluate to which extent UCNE regions overlap with those collected in other databases due to the different scopes and methods of these resources. We clarified why we selected this database, which was precisely based on the comments mentioned by the reviewer, i.e., sufficient evolutionary distance to identify functional regions confidently and high quality genome assemblies. Regarding the extent to which UCNEs are conserved in other species, the UCNEdatabase indeed provides additional information with respect to UCNE orthologues in other vertebrate species, including reptiles, amphibians, and fish. As we consider that this comparison is beyond the scope of the present study, it was not included in the main manuscript. However, to address this interesting point raised by the reviewer, we evaluated the proportion of UCNEs found in different species with respect to those annotated for human-chicken. As show in the figure below, UCNEs are conserved to a larger extent up to Xenopus tropicalis, for which most of the UCNEs annotated for human-chicken have corresponding orthologues in this species. Although to a minor extent, UCNEs are also conserved across more distant species (e.g., fish), for which approximately half of the UCNEs annotated for human-chicken have orthologues.

    1. The authors briefly describe how potential target genes were assigned to UCNEs in the Method section (section begins on page 21, line 407 and on lines 379-381). They note that they used the Genomic Regions Enrichment of Annotations Tool (GREAT). Can the authors provide additional details on what this method does and also provide a high level sentence in the Results section on page 8, lines 144-146, on how target genes were assigned to UCNEs, as well as in the Discussion on page 16, line 289. Based on the Methods description on lines 409-414, a UCNE was associated with any gene within 1Mb of the UNCE that was expressed in retina and whose curated regulatory domains overlapped the UNCE. Is this correct? How were the regulatory domains curated (line 411-412)? Can the authors please clarify this important point.

    Additional details about the GREAT algorithm have been included in the Methods section (lines 410-424) as well as a high-level sentences in the Results (lines 140-141) and Discussion (lines 276-280). It is correct that a UCNE was associated with any gene within 1Mb of the UCNE that was expressed in retina and whose curated regulatory domains overlapped the UCNE. As stated now (lines 420-424), these regulatory domains are supported by experimental evidence demonstrating that a gene is directly regulated by an element located beyond of its putative regulatory domain. We specified these domains for the utilized version of GREAT as well.

    1. There are other approaches for linking putative transcriptionally active genomic regions with target genes in specific tissues or cell types through correlation analysis of ATAC-seq peaks (chromatin accessibility) and gene expression in RNA-seq data. A commonly used peak-to-gene linkage method is implemented in ArchR (Granja et al, 2021, PMID: 33633365). The authors note that they used scATAC-seq gene scores and scRNA-seq gene expression to further characterize the proposed target genes of UCNEs. It is not clear what the authors mean by "scATAC-seq gene scores". Please define "scATAC-seq gene score" and add a reference. Can the authors compare the target gene mapping to UCNEs using GREAT and gene expression filtering to that using peak-to-gene linkage based on retina tissue or single cell ATAC-seq and RNA-seq data?

    We have defined explicitly what scATAC-seq gene scores are in the Methods section (lines 436-437) along with its reference. We have also addressed this important point and compared the overlap between the set of target genes predicted by GREAT and those assigned by the peak-to-gene linkage method implemented by ArchR. Details of this analysis, its results, interpretation are included in the Methods (lines 441-446), Results (lines 158-163; Supplementary Table 4), and Discussion (lines 280-282) sections.

    1. For gene expression filtering (line 419) the authors quantify transcript expression of retina from FASTQ samples using the Kallisto method, and then note that they did the filtering on gene expression levels (TPM<0.5). Please add details on how you went from transcript expression levels to gene expression levels for the filtering, or was the filtering performed at the transcript level?

    We have added details on how transcript-level quantification estimates were summarized at gene level, for which the filtering was performed (lines 429-430).

    1. The authors use the words "active UCNE", first mentioned in the Results on line 144. Can the authors define what they mean by "active UCNE". What information/evidence is used to ascertain that a UCNE is indeed active. Overlap of a UCNE with a chromatin accessibility region from ATAC-seq or DNAase-Seq would only suggest that the UCNE may be active. Intersection with enhancer activity measured with in vivo enhancer reporter assays in transgenic mice from the VISTA enhancer browser provides stronger evidence of transcriptional activity. The authors might want to distinguish between putatively actively and active based on the functional support.

    We thank the reviewer for this relevant comment to address the nuances of defining active UCNEs. The reviewer’s assumptions are correct and hence these terms were clarified throughout the entire text. The term functional is now only used when referring to UCNEs for which there is functional support (e.g., PAX6-associated UCNE in line 193) .

    1. The authors assessed the significance of depletion of common variation (MAF>1%) in the UCNE regions compared to a background of randomly selected genomic regions. In generating the random distribution of regions, did the authors match on the distribution of distances of the UNCEs to the TSS of genes in the randomly selected regions? This may be a confounder. Also, in the legend of Figure 3, lines 191-192, it is stated: "Variant population frequencies within putative retinal UCNEs normalized to a background composed of randomly selected sequences (see Methods).", but we did not find a description of this analysis in the Methods section.

    Evaluating the potential confounding effect of the genomic background was indeed a very important point. We have now incorporated the details showing the suitability of such background well as a detail description of how such background was generated (lines 479-481; Figure S1). Additionally, to further support our analysis demonstrating the depletion of common variation within UCNEs, we have included an evaluation of the distribution of genome-wide residual variation intolerance score (gwRVIS) values (PMID: 33686085) compared to this background of randomly selected genomic regions in a human reference cohort (lines 173-178; 487-495; Fig. 3C).

    1. In regards to intersecting UCNEs with epigenetic marks that detect active or repressed enhancers in retina, the authors used data from Aldiri et al 2017 that measured epigenetic changes during retinal development. Did this dataset contain epigenetic measurements in adult retina? The authors might want to consider using the epigenetic marks/ChIP-seq data from adult human retina in Cherry et al. PNAS 2020 (PMID: 32265282)

    We have incorporated the adult-stage data suggested by the reviewer to provide a more comprehensive characterization. Details about the integration of this dataset as well as the results and their corresponding interpretation have been incorporated in the Methods (lines 372-374), Results (lines 115-117; Supplementary Table 2), and Discussion (lines 271-273) sections accordingly.

    1. With respect to the examination of rare variants that may be associated with rare eye disease in retina active UCNEs, for the interpretation of the results, it would be helpful to get more information on the distribution of rare variants found in UCNEs associated in this study to known IRD genes in all affected individuals in families, if this information is available in the 100,000 Genomes Project.

    Although it is indeed a relevant point, this information cannot be retrieved in the 100,000 Genomes Project. As it is a restricted research environment, we are only allowed to query sequencing data corresponding to participants enrolled within the framework of our specific sub-domain, namely “Hearing and Sight”, and therefore evaluating the distribution of rare variants in all affected individuals is not feasible.

    1. In the Methods section on lines 450-451, the authors mention that they performed variant screening of retinal disease genes, referencing the Genomics England Retinal Disorder panel and Martin et al., 2019. Can the authors add to the Methods and Results sections how many retinal disease genes were initially tested. Also, to get a sense of the specificity of the overlap of rare variants in the 100k Genome Project cohort with UNCE regions, it would be informative to show a distribution of the number of rare variants <0.5% that passed the filtering in gnomAD per eye disease gene before the overlap with UNCEs.

    We specified the number of retinal genes that were tested in the Method section (line 471). In addition, as suggested by the reviewer, we generated allele frequency distributions for all variants retrieved within a selection of 25 disease-gene associated loci and their corresponding UCNEs in order to assess the specificity of the overlap between rare variants and UCNE regions (lines 181-182, 496-501; Figure S2).

    1. The authors found "an ultrarare SNV (chr11:31968001T>C) within a candidate cis-regulatory UCNE located ~150kb upstream of PAX6. This variant was found in four individuals of a family segregating autosomal dominant foveal abnormalities". They tested the functional effect of this element with a reporter assay in zebrafish and found that the UNCE affects expression in the eye, forebrain, and neural tube. It would add further value if the authors were to test the effect of this SNV in the UCNE on the reporter expression pattern, using CRISPR/cas9?

    That is a very relevant point. We have tested the effect of this SNV in the UCNE on the reported expression pattern using the same experimental setup that we used for testing the wild-type construct, namely transgenic enhancer zebrafish assays. However, we could not obtain conclusive results, most likely due to the limitations posed by testing these regions outside their native genomic context. Therefore, additional experimental work (e.g., CRISPR-based) should be performed to address this question. This is, however, beyond the scope of the present study, for which the main focus was the identification and functional annotation of ultraconserved cCREs. We have incorporated the details, results, and interpretation of the assays performed mutant construct in the Methods (lines 525-527; Supplementary Table 12), Results (lines 235-238; Supplementary Table 10), and Discussion (lines 350-353) sections.

    1. The authors found rare variants in UCNEs linked to 45 IRD genes. Can the authors provide additional information on the functional genomic annotations of these UCNEs and distance to the target genes. The UCNEs were characterized with respect to their genic features in the original paper (UCNEbase, Dimitrieva et al., NAR, 2013), e.g., intergenic, intronic and 3'/5' UTR. Also, it would be useful for clinical applications to provide the start and end positions of the UCNEs that contain the rare variants associated with their 45 IRD genes in Supplementary Table 6.

    Additional functional genomic annotations, genic features following those of the original UCNE paper, and the distance to the TSS of these 45 target disease-associated genes have been incorporated in (new) Supplementary Table 5. The start and end positions of the UCNEs that contain the rare variants have also been indicated in new Supplementary Table 7.

    1. A total of 724 target genes were assigned to 1,487 UCNEs that displayed candidate cis-regulatory activity. Given the interest in using UCNEs to help identify potential pathogenic mutations that lead to IRDs, can the authors note in the Results section how many of the 724 target genes are IRDs.

    We thank the reviewer for this important point. From the total of 724, a total of 27 genes are annotated as IRD genes, of which (interestingly) 23 were kept as found to be expressed in the retina. This has been clarified in the Results section (line 166-168).

    1. In the Discussion on page 15 line 259, can the authors clarify if variation found in UCNEs were only associated with rare disease or also with common diseases.

    We have clarified that variation found in UCNEs has only been associated with rare diseases (line 247).

    Minor edits:

    1. In abstract, the authors might consider changing the words "rare eye diseases" on lines 20 and 22 to "rare retina degeneration diseases", and on lines 88-89.

    We thank the reviewer for this comment. However, we consider that rare eye diseases is a more suitable term for our purpose as it includes diseases primarily characterized by stationary and non-progressive phenotypes such as North Carolina Macular Dystrophy and fovea hypoplasia.

    1. In the Introduction on line 49, there seems to be a typo in the number of UCNE regions reported. 4,135 UCNE regions is supposed to be 4,351, based on the original paper (https://academic.oup.com/nar/article/41/D1/D101/1057253).

    We have corrected this typo accordingly.

    1. In the introduction on lines 75-76, these references: Lyu et al., Cell Reports 2021, PMID: 34788628, and Zhang et al., Trends in Genetics 2023, PMID: 37423870, could be added to the following sentence to provide additional: "This cellular complexity is the result of spatiotemporally controlled gene expression programs during retinal development”.

    We have now included these relevant references.

    1. On lines 77 and 84, I would write IRD as plural: IRDs.

    This has been amended in the new version.

    1. In introduction on lines 89-90, it can be further added that you provide experimental support for an ultra rare SNV in a cis regulatory UCNE affecting PAX6.

    We have explicitly stated that we provided functional evidence for this UCNE.

    1. On line 98, the authors refer to Figure 1A when noting that the integration of UCNEs with multi-omics data in human retina allows to evaluate the regulatory capacity of UCNEs across the major developmental stages of human retina. However panel A in figure 1 does not seem to show this point. It shows the comparison of elements across species. Please make appropriate changes to the main text and figure legend.

    We have made the appropriate changes and located the reference to this figure in a more relevant part of the text (line 87).

    1. Please explain what the names appended to the gene symbols in the first column "UCNE ID" in Supplementary Tables 1 and 2 refer to.

    We have clarified what these refer to.

    1. On line 145, can the authors clarify what they mean by "active gene expression in the retina". Is this just another way of referring to genes found to be expressed in retina? If so, it might be clearer just to write: "We annotated the identified active UCNEs to assign them potential target genes and thus assess their association with genes expressed in the retina"

    We indeed meant genes found to be expressed in retina. As this phrasing might not be completely clear, we have now changed to the wording suggested by the reviewer.

    1. One line 156, I would write "regulation of transcription" as listed in the gene ontology terms in Figure 2C, instead of "regulation of gene expression". The authors might want to add this to the Discussion. Can the authors include the full gene set enrichment results from Enrichr in a supplementary table at Padj<0.05 since only the top gene sets are shown in Fig. 2C (at P<1E-13)?

    We changed the term to “regulation of transcription” to keep the nomenclature consistent to that of Figure 2C. We have also provided a full gene set enrichment from Enrichr as well (Supplementary Table 3).

    1. On page 12, line 214, what does "EH38E1530321" Stand for? It seems to specify a distal enhancer-like signatures in bipolar neurons, but I couldn't find this ID in any database.

    This refers to the identifier of ENCODE:

    https://screen-v2.wenglab.org/search/?q=EH38E1530321&assembly=GRCh38

    Additionally, when mentioning a specific UCNE, VISTA enhancer, or ENCODE cCRE (as in this case), we have explicitly included its corresponding identifier.

    1. In the Methods section on lines 391-392, can the authors give some high-level explanation of the unconstrained integration method: "Single-nucleus RNA-seq of the same tissue and timepoints (GSE183684) were integrated using the unconstrained integration method". Also, can they comment on how retinal cell class identities were assigned (line 393). Was it based on known markers or on previous identification of cell classes and highly variable genes between clusters?

    We have included a high-level explanation of the unconstrained method in the Methods section (lines 387-392). We also clarified that the assignment of cell class identities was based on known markers (line 394).

    1. In the integration of UCNEs with bulk and single cell ATAC-seq and Dnase hypersentitivity regions, can the authors note in the Methods section (lines 400-404) what peak width was used to test for overlap with the UCNEs.

    We have specified the peak widths that were used for the overlap with UCNEs (lines 397 and 403).

    1. On line 436, the word 'and' is missing between "(SNVs, and indels < 50bp)" and "large structural variants".

    This has been corrected.

    1. On lines 443-444, please provide references to the computational tools listed. Please note if default settings/parameters were used.

    We have specified that default parameters were used in the analysis (line 464).

    1. In the following sentence in the Methods section on lines 447-449, it is not clear in the Results section how this was used in the flow of the analysis, and how many cases showed such a similarity in phenotype: "For each candidate variant, we compared the similarities between the participant phenotype (HPO terms) and the ones known for its target gene through literature search and clinical assessment by the recruiting clinician when possible." Can the authors add more detail to the Results section.

    As the evaluation of the candidate variants was essentially performed on a case-by-case basis, we opted to include a rather general description of the workflow, which indeed included a comparison of the reported phenotypes with those associated with the putative target gene. An example of such comparison has been included in the Results (lines 186-187) section regarding the cases for which a NCMD-like phenotype was suspected.

    1. It would be helpful to have a table that describes the different omics datasets used in the paper, with some basic annotations (tissue type, sample size, reference).

    This has now been incorporated in Supplementary Table 11.

    1. Can the authors add references to their sentence in the Discussion on page 17 lines 299-301: "As has been shown before, the phenotype caused by a coding mutation of a developmental gene can be different from the phenotype caused by a mutation in a CRE controlling spatiotemporal expression of this gene."

    We clarified that this phrase referred to the case of PRDM13, for which bi-allelic coding pathogenic variants are linked to hypogonadotropic hypogonadism and perinatal brainstem dysfunction in combination with cerebellar hypoplasia (Whittaker et al., 2021), while non-coding variants within its regulatory regions are associated with NCMD.

    Reviewer #2:

    Minor discretionary suggestions for improving the presentation:

    1. Wherever a specific UCNE, Vista enhancer or ENCODE cCRE is mentioned, the element should be identified by name or accession code: For instance (Iine 212): "this variant is located within a UCNE (PAX6_Veronica) that is catalogued as a cCRE in ENCODE (EH38E1530321)". UCNE names are particularly important, because they are systematically used as identifiers in the supplemental Tables and thus would enable the reader to easily find additional information about the element mentioned in the main text.

    We have now explicitly included all corresponding identifiers throughout the text.

    1. I also recommend inclusion of the UCNE, Vista and ENCODE cCRE tracks in all genome browser screen shots. The UCNE track is currently included only in Figure 1. Vista and ENCODE cCRE tracks are missing in all browser views.

    We have now included UCNE, VISTA, and ENCODE cCREs tracks in the main genome browser figures (Figures 1 and 4).

    1. Supplementary Table 6: It would be useful to indicate for each variant, the type of ophthalmological disorder (Table S5, column C) it is associated with.

    We agree this is a relevant point. However, due to limitations in the (bulk) export of clinical information from the protected Research Environment of Genomics England, inclusion of this type of information is not feasible.

    1. Fig S2 and supplementary Table 3 are not referred to in the main Text.

    We have corrected this and updated the figure and table accordingly.

    1. Supplementary Table 8: The Table caption should be expanded.

    The contents of each column should be explained. For instance, column F: what means Homo_sapiens|M01298_1.94d|Zoo_01|2337? Where does this information come from, what data and software resources were used?

    We have expanded the caption of this table to clarify this output, which is derived from the QBiC-Pred tool, a software used for predicting quantitative TF binding changes due to nucleotide variants.

    1. Line 401 probable typo: 103-105 days (103-125?)

    Indeed, this typo has now been corrected.

    Reviewer #3:

    1) Given that UCNE only accounts for a small fraction of gene regulatory elements, this study is likely with low sensitivity in terms of identifying potential regulatory mutations. Although one would expect that variants in UCNE are more likely to be pathogenic, it is hard to extrapolate from the results to assess the contribution of gene regulatory variant to the disease.

    We agree that restricting our analysis to these particular regions is one of the limitations of the study, as stated in the Discussion section (line 304-311). However, one of our main aims was to provide a strategy to reduce the search space for pathogenic variants with a potential regulatory effect. Given the substantial body of literature supporting a regulatory role for these regions and, particularly, the availability of already-existing functional data, we considered that this set of regions could represent a suitable target for such analysis. Indeed, the features evaluated, and the methods presented in this study could be extrapolated and applied in other settings involving other candidate regulatory regions and/or tissues of interest, and their associated disease-phenotypes, for which, in any case, the overall contribution of regulatory pathogenic variation to disease might vary greatly.

    1. I am wondering how many UCNE overlaps with open chromatin regions specific to the fetal retina and how many UCNE overlaps with adult only. Are UCNEs enriched for developmental genes? If so, how many patients are due to developmental defect?

    We have now integrated into our analysis epigenetic measurements in adult retina, in particular the candidate cis-regulatory elements nominated by Cherry et al. PNAS 2020 (PMID: 32265282) based on accessible chromatin and enrichment for active enhancer-related histone modifications in adult human retina. Details about the integration of this dataset as well as the results and their corresponding interpretation have been incorporated in the Methods (lines 372-374), Results (115-117; Supplementary Table 2), and Discussion (lines 271-273) sections accordingly. In particular, out of the 111 UCNEs identified to display the active enhancer mark H3K27ac, 33 were found to maintain this signature at adult stage. Regarding the specific question from the reviewer, the majority of UCNEs overlapping with open chromatin regions are specific to the fetal retina (1,346), with only 7 UCNEs overlapping with open chromatin regions exclusively in adult state. This indeed further supports the major role of the identified candidate cis-regulatory UCNEs in the regulation of developmental gene expression programs, which was already suggested by the Gene Ontology enrichment analysis performed using the set of UCNE target genes as input (Figure 2C; new Supplementary Table 3). As far as the number of patients with development defects that were included in this study, these included: corneal abnormalities (n=62), Leber congenital amaurosis (n=142), ocular coloboma (n=111), developmental foveal and macular dystrophy (n=230), developmental glaucoma (n=94), anophthalmia or microphthalmia (n=120).

    1. I am wondering if the 431 ultrarare variants found in the UCNEs is higher than expected. This can be tested by comparing patients without visual disorders.

    Although it is indeed a relevant point, retrieving sequencing data from patients without visual disorders is not feasible for us. As it is a restricted research environment, we are only allowed to query sequencing data corresponding to participants enrolled within the framework of our specific sub-domain, namely “Hearing and Sight”, and therefore evaluating additional patients from other sub-domains is not doable. Based on previous studies and our observations, common variants are precisely the ones depleted within UCNEs, while ultrarare variation seems to occur at levels comparable to those observed elsewhere in the genome. Therefore, it is reasonable to speculate that this amount of ultrarare variants is not higher than expected as compared to patients without visual disorders. To further demonstrate the high intolerance of UCNEs to common variation, we have included an evaluation of the distribution of genome-wide residual variation intolerance score (gwRVIS) values compared to a set of randomly selected genomic regions in a human reference cohort (lines 173-178; 487-495; Fig. 3C). Additionally, to address this question further, we have also generated allele frequency distributions for all variants retrieved within a selection of 25 disease-gene associated loci and their corresponding UCNEs in order to assess the specificity of the overlap between rare variants and UCNE regions (lines 181-182, 496-501; Figure S2).

    1. It seems that the ultrarare variants listed in sup table 6 are more abundant in a small number of genes. Is this due to the number/size of UCNEs is larger in these genes?

    Indeed, the clustering of UCNEs in genomic regions containing genes coding for transcription factors and developmental regulators (e.g., OTX2, PAX6, ZEB2) seems to be one of their intrinsic properties, hence the observed enrichment for a small number of genes. One reason can be that these neighboring UCNEs cooperate to achieve higher degrees of tissue-specific regulatory accuracy needed for these genes.

    1. The variant in the putative Pax6 gene regulatory element is intriguing. It would be much more informative if the enhancer with and without the variant is tested in parallel in fish.

    That is a very relevant point. We have now tested the effect of this SNV in the UCNE on the reported expression pattern using the same experimental setup that we used for testing the wild-type construct, namely transgenic enhancer zebrafish assays. However, we could not obtain conclusive results, most likely due to the limitations posed by testing these regions outside their native genomic context. We have incorporated the details, results, and interpretation of the assays performed mutant construct in the Methods (lines 525-527; Supplementary Table 12), Results (lines 235-238; Supplementary Table 10), and Discussion (lines 350-353) sections.

    1. (optional) it would be quite interesting to check the phenotype in fish or mice with the element repressed via technique such as CRISPRi.

    Indeed, we fully agree that CRISPR-based techniques would be the ideal experimental approaches to further validate the functionality of the PAX6-associated UCNE and the identified variant in their native genomic context. Conducting these detailed and focused mechanistic studies is, however, beyond the scope of the present work, for which the main focus was the identification and functional annotation of ultraconserved cCREs.

  3. Review Commons

    Note: This preprint has been reviewed by subject experts for Review Commons. Content has not been altered except for formatting.

    Learn more at Review Commons

    Referee #3

    Evidence, reproducibility and clarity

    In the report, the authors combined bulk and single cell multi-omics data from the retina to identify ultraconserved noncoding elements (UCNE) that might function as regulatory elements based on chromatin openness, DNAse sensitive regions, and histone marks. The candidate UCNEs are intersected with whole genome sequencing data from 3.220 patients with rare eye disease to identify potential rare variants that might affect the activity of UNCE. The goal of the project is intriguing.

    My comments are listed below:

    1. Given that UCNE only accounts for a small fraction of gene regulatory elements, this study is likely with low sensitivity in terms of identifying potential regulatory mutations. Although one would expect that variants in UCNE are more likely to be pathogenic, it is hard to extrapolate from the results to assess the contribution of gene regulatory variant to the disease.
    2. I am wondering how many UCNE overlaps with open chromatin regions specific to the fetal retina and how many UCNE overlaps with adult only. Are UCNEs enriched for developmental genes? If so, how many patients are due to developmental defect?
    3. I am wondering if the 431 ultrarare variants found in the UCNEs is higher than expected. This can be tested by comparing patients without visual disorders.
    4. It seems that the ultrarare variants listed in sup table 6 are more abundant in a small number of genes. Is this due to the number/size of UCNEs is larger in these genes?
    5. The variant in the putative Pax6 gene regulatory element is intriguing. It would be much more informative if the enhancer with and without the variant is tested in parallel in fish.
    6. (optional) it would be quite interesting to check the phenotype in fish or mice with the element repressed via technique such as CRISPRi.

    Significance

    Given that UCNE only accounts for a small fraction of gene regulatory elements, this study is likely with low sensitivity in terms of identifying potential regulatory mutations. Although one would expect that variants in UCNE are more likely to be pathogenic, it is hard to extrapolate from the results to assess the contribution of gene regulatory variant to the disease.

  4. Review Commons

    Note: This preprint has been reviewed by subject experts for Review Commons. Content has not been altered except for formatting.

    Learn more at Review Commons

    Referee #2

    Evidence, reproducibility and clarity

    Starting with a comprehensive analysis of multi-omics data, the authors identify a subset of ultra-conserved non-coding elements (UCNEs) likely to play a role in retinal development. Restricting subsequent analyses to these genomic regions, they identify ultra-rare mutations associated with eye disease, using data from the 100k genome project. They then follow-up on one newly discovered, disease associated UCNE and a presumably causal mutation within this UCNE. The follow-up experiments involve fine mapping of the spatiotemporal expression pattern of a UCNE-driven reporter gene in zebra fish, as well as a re-examination of the disease phenotype in carriers of the mutation.

    The computational pipeline used for prioritization of UCNEs is sound. The evidence supporting the claim that the identified ultra-rare SNV is causal is highly convincing. This study also constitutes proof of concept for a novel methodology to search for causal disease associated mutations in the nowadays still under-investigated non-coding part of the genome.

    The paper is clearly written. Methods are described in enough detail to allow for reproduction of the results. Overall, this study is of high scientific quality, self-contained, and complete. Publication in a peer-reviewed journal should not be delayed by additional, perhaps interesting but non-essential experiments.

    However, if the authors intend to undertake similar studies in the future, I would recommend to carry out the reporter-gene assays in zebrafish with both the wild-type and the mutant version of the UCNE. Comparison of the spatiotemporal expression patterns of the two alleles could provide valuable insights into the mechanism of action of the deleterious mutation under investigation.

    Minor discretionary suggestions for improving the presentation:

    Wherever a specific UCNE, Vista enhancer or ENCODE cCRE is mentioned, the element should be identified by name or accession code: For instance (Iine 212):

    "this variant is located within a UCNE (PAX6_Veronica) that is catalogued as a cCRE in ENCODE (EH38E1530321)"

    UCNE names are particularly important, because they are systematically used as identifiers in the supplemental Tables and thus would enable the reader to easily find additional information about the element mentioned in the main text.

    I also recommend inclusion of the UCNE, Vista and ENCODE cCRE tracks in all genome browser screen shots. The UCNE track is currently included only in Figure 1. Vista and ENCODE cCRE tracks are missing in all browser views.

    Supplementary Table 6: It would be useful to indicate for each variant, the type of ophthalmological disorder (Table S5, column C) it is associated with.

    Fig S2 and supplementary Table 3 are not referred to in the main Text.

    Supplementary Table 8: The Table caption should be expanded. The contents of each column should be explained. For instance, column F: what means Homo_sapiens|M01298_1.94d|Zoo_01|2337? Where does this information come from, what data and software resources were used?

    Line 401 probable typo: 103-105 days (103-125?)

    Significance

    This work is of interest to different research communities: Biomedical researchers working on neural disorders, human geneticists engaged in GWAS studies, computational biologists trying to make sense out of omics data, molecular biologists exploring the "dark matter" of the genome, and finally the small community tackling the enigma of UCNEs. As with many omics papers, the most valuable parts of this study are in the supplemental tables, in particular tables 1,4, and 6. It can be hoped that some prospective readers will follow up on the leads presented in these tables.

    The detailed computational and experimental characterization of a likely causal ultra-rare disease associated mutation may serve as a guiding and motivating example for medical geneticists working on other syndromes.

    Back to UCNEs: They are enigmatic entities, which so far have largely resisted molecular and physiological characterization. It took 10 years to finally uncover a phenotype in ko mice missing one or several UCNEs, after the surprising initial observation that such mice were viable and fertile. The difficulties in studying the function of UCNEs may be due to their conjectured pleiotropic activity in different cell types at different developmental stages, their apparent cooperative interactions with many other control elements (limiting the power of reporter gene assays with single elements), and their putative involvement in morphogenetic processes (minimizing the relevance of epigenetic data collected from cell lines). In view of these considerations, I consider UCNE research starting from human disease phenotypes more promising, than ab initio approaches using reverse genetics in model organisms.

    The impact of this paper is potentially very high. Note the following statement in the paper:

    "For each instance for which only the UCNE variant remained as candidate, we placed a clinical collaboration request with Genomics England."

    We thus can expect more exiting stories from the same team. The strategy and computational pipeline introduced here are of course applicable to other congenital diseases, and it can reasonably be hoped that researchers inspired by this study will apply components of the methodology in other contexts. The prioritization of UCNEs in studying the "dark matter" of the genome and the "missing heredity" would likely lead to new insights into the function of these enigmatic elements and the reasons for their extreme conservation.

    My background: I'm a bioinformatician with first training in molecular genetics. My research focus is on gene regulation: promoters, enhancers, transcription factor binding sites. I also made some contributions to the UCNE field, having co-developed UNCEbase with Slavica Dimitrieva. On the other hand, I don't claim to be an expert in medical genetics, and more specifically, I know very little about eye diseases and retinal development.

  5. Review Commons

    Note: This preprint has been reviewed by subject experts for Review Commons. Content has not been altered except for formatting.

    Learn more at Review Commons

    Referee #1

    Evidence, reproducibility and clarity

    In the manuscript entitled "Multi-omics analysis in human retina uncovers ultraconserved cis-regulatory elements at rare eye disease loci", Soriano et al. explore the role of ultraconserved noncoding elements (UCNEs) as candidate cis-regulators (cCREs) in developing and adult human retina by integrating UCNEs with multi-omic data from fetal and adult human retina. They further examine the potential contribution of UCNEs to rare eye diseases by testing for rare variants that fall in UCNEs that are potentially active in retina in patients with inherited retina degeneration diseases (IRDs), using whole genome sequencing (WGS) data from Genomics England. This work is based on the assumption that genomic regions under strong evolutionary constrain are good predictors of important functionality. The authors use 4,135 UCNEs identified in Dimitrieva and Bucher, 2013, defined as sequences with more than 95% identity between human and chicken with >200bp in length. By integrating bulk tissue RNA-Seq, single cell (sc)RNA-Seq, DNAase-Seq, scATAC-Seq and ChIP-Seq from retinal tissues from different developmental stages and adulthood with the 4,135 UCNEs, they identified one third of UCNEs (1,487) that may be associated with various retinal development stages. Using bioinformatic approaches they identify potential target genes for 724 UCNEs that may be regulated in cis. Based on ATAC-Seq, the authors show that 834 UCNEs fall in open chromatin regions, and of these 111 have active histone markers, supporting an active functional role for the identified UCNEs. The authors demonstrate an application of this approach in identifying potential pathogenic noncoding variants for rare eye diseases. Using WGS from 100,000 Genomes project the authors identify 45 UCNEs that were bioinformatically associated with mendelian IRD genes that contain rare variants, proposing potential solutions for unsolved cases. They demonstrate the utility of their analysis by highlighting an ultra-rare SNP they identified in a UCNE candidate CRE located upstream of PAX6 that belongs to a family of genes with foveal abnormalities. They confirm the effect of the UCNE 150k upstream of PAX6 on gene expression in the eye, forebrain and neuronal tube using an enhancer reporter assay in Zebrafish.

    The authors demonstrate a strategy for narrowing down the genomic space to putative functional noncoding regions involved in development and that may contribute to rare disease. This work that integrates various genomic modalities is of broad interest and can be applied to other diseases and tissues. There are several points that should be addressed to assess the robustness of the results and strengthen the conclusions of the paper, and in some places to better clarify the analyses conducted.

    Major edits:

    1. In the introduction the authors list 4 databases of ultra-conserved non-coding elements, and choose to work with the UCNEbase that detects UCNE regions by comparing human and chicken, as they state that it is one of the most comprehensive resources. Can the authors comment on the extent of overlap of UCNE regions identified in this database compared to other databases. It would also be helpful to add an explanation in the Introduction on the advantage of comparing the human genome to the chicken genome, e.g., is chicken the most distant vertebrate with a high quality genome, or another reason? And can the authors comment on the extent of conservation of the UCNEs compared to other species - in the UCNEbase paper, orthologues for the UCNEs are identified in 18 vertebrate species including reptiles, amphibians and fish.
    2. The authors briefly describe how potential target genes were assigned to UCNEs in the Method section (section begins on page 21, line 407 and on lines 379-381). They note that they used the Genomic Regions Enrichment of Annotations Tool (GREAT). Can the authors provide additional details on what this method does and also provide a high level sentence in the Results section on page 8, lines 144-146, on how target genes were assigned to UCNEs, as well as in the Discussion on page 16, line 289. Based on the Methods description on lines 409-414, a UCNE was associated with any gene within 1Mb of the UNCE that was expressed in retina and whose curated regulatory domains overlapped the UNCE. Is this correct? How were the regulatory domains curated (line 411-412)? Can the authors please clarify this important point.

    There are other approaches for linking putative transcriptionally active genomic regions with target genes in specific tissues or cell types through correlation analysis of ATAC-seq peaks (chromatin accessibility) and gene expression in RNA-seq data. A commonly used peak-to-gene linkage method is implemented in ArchR (Granja et al, 2021, PMID: 33633365). The authors note that they used scATAC-seq gene scores and scRNA-seq gene expression to further characterize the proposed target genes of UCNEs. It is not clear what the authors mean by "scATAC-seq gene scores". Please define "scATAC-seq gene score" and add a reference. Can the authors compare the target gene mapping to UCNEs using GREAT and gene expression filtering to that using peak-to-gene linkage based on retina tissue or single cell ATAC-seq and RNA-seq data?

    For gene expression filtering (line 419) the authors quantify transcript expression of retina from FASTQ samples using the Kallisto method, and then note that they did the filtering on gene expression levels (TPM<0.5). Please add details on how you went from transcript expression levels to gene expression levels for the filtering, or was the filtering performed at the transcript level?

    1. The authors use the words "active UCNE", first mentioned in the Results on line 144. Can the authors define what they mean by "active UCNE". What information/evidence is used to ascertain that a UCNE is indeed active. Overlap of a UCNE with a chromatin accessibility region from ATAC-seq or DNAase-Seq would only suggest that the UCNE may be active. Intersection with enhancer activity measured with in vivo enhancer reporter assays in transgenic mice from the VISTA enhancer browser provides stronger evidence of transcriptional activity. The authors might want to distinguish between putatively actively and active based on the functional support.
    2. The authors assessed the significance of depletion of common variation (MAF>1%) in the UCNE regions compared to a background of randomly selected genomic regions. In generating the random distribution of regions, did the authors match on the distribution of distances of the UNCEs to the TSS of genes in the randomly selected regions? This may be a confounder.

    Also, in the legend of Figure 3, lines 191-192, it is stated: "Variant population frequencies within putative retinal UCNEs normalized to a background composed of randomly selected sequences (see Methods).", but we did not find a description of this analysis in the Methods section.

    1. In regards to intersecting UCNEs with epigenetic marks that detect active or repressed enhancers in retina, the authors used data from Aldiri et al 217 that measured epigenetic changes during retinal development. Did this dataset contain epigenetic measurements in adult retina? The authors might want to consider using the epigenetic marks/ChIP-seq data from adult human retina in Cherry et al. PNAS 2020 (PMID: 32265282)
    2. With respect to the examination of rare variants that may be associated with rare eye disease in retina active UCNEs, for the interpretation of the results, it would be helpful to get more information on the distribution of rare variants found in UCNEs associated in this study to known IRD genes in all affected individuals in families, if this information is available in the 100,000 Genomes Project.
    3. In the Methods section on lines 450-451, the authors mention that they performed variant screening of retinal disease genes, referencing the Genomics England Retinal Disorder panel and Martin et al., 2019. Can the authors add to the Methods and Results sections how many retinal disease genes were initially tested. Also, to get a sense of the specificity of the overlap of rare variants in the 100k Genome Project cohort with UNCE regions, it would be informative to show a distribution of the number of rare variants <0.5% that passed the filtering in gnomAD per eye disease gene before the overlap with UNCEs.
    4. The authors found "an ultrarare SNV (chr11:31968001T>C) within a candidate cis-regulatory UCNE located ~150kb upstream of PAX6. This variant was found in four individuals of a family segregating autosomal dominant foveal abnormalities". They tested the functional effect of this element with a reporter assay in zebrafish and found that the UNCE affects expression in the eye, forebrain, and neural tube. It would add further value if the authors were to test the effect of this SNV in the UCNE on the reporter expression pattern, using CRISPR/cas9?
    5. The authors found rare variants in UCNEs linked to 45 IRD genes. Can the authors provide additional information on the functional genomic annotations of these UCNEs and distance to the target genes. The UCNEs were characterized with respect to their genic features in the original paper (UCNEbase, Dimitrieva et al., NAR, 2013), e.g., intergenic, intronic and 3'/5' UTR. Also, it would be useful for clinical applications to provide the start and end positions of the UCNEs that contain the rare variants associated with their 45 IRD genes in Supplementary Table 6.
    6. A total of 724 target genes were assigned to 1,487 UCNEs that displayed candidate cis-regulatory activity. Given the interest in using UNCEs to help identify potential pathogenic mutations that lead to IRDs, can the authors note in the Results section how many of the 724 target genes are IRDs.
    7. In the Discussion on page 15 line 259, can the authors clarify if variation found in UNCEs were only associated with rare disease or also with common diseases.

    Minor edits:

    1. In abstract, the authors might consider changing the words "rare eye diseases" on lines 20 and 22 to "rare retina degeneration diseases", and on lines 88-89.
    2. In the Introduction on line 49, there seems to be a typo in the number of UCNE regions reported. 4,135 UCNE regions is supposed to be 4,351, based on the original paper (https://academic.oup.com/nar/article/41/D1/D101/1057253).
    3. In the introduction on lines 75-76, these references: Lyu et al., Cell Reports 2021, PMID: 34788628, and Zhang et al., Trends in Genetics 2023, PMID: 37423870, could be added to the following sentence to provide additional: "This cellular complexity is the result of spatiotemporally controlled gene expression programs during retinal development,"
    4. On lines 77 and 84, I would write IRD as plural: IRDs.
    5. In introduction on lines 89-90, it can be further added that you provide experimental support for an ultra rare SNV in a cis regulatory UCNE affecting PAX6.
    6. On line 98, the authors refer to Figure 1A when noting that the integration of UCNEs with multi-omics data in human retina allows to evaluate the regulatory capacity of UCNEs across the major developmental stages of human retina. However panel A in figure 1 does not seem to show this point. It shows the comparison of elements across species. Please make appropriate changes to the main text and figure legend.
    7. Please explain what the names appended to the gene symbols in the first column "UCNE ID" in Supplementary Tables 1 and 2 refer to.
    8. On line 145, can the authors clarify what they mean by "active gene expression in the retina". Is this just another way of referring to genes found to be expressed in retina? If so, it might be clearer just to write: "We annotated the identified active UCNEs to assign them potential target genes and thus assess their association with genes expressed in the retina"
    9. One line 156, I would write "regulation of transcription" as listed in the gene ontology terms in Figure 2C, instead of "regulation of gene expression". The authors might want to add this to the Discussion. Can the authors include the full gene set enrichment results from Enrichr in a supplementary table at Padj<0.05 since only the top gene sets are shown in Fig. 2C (at P<1E-13)?
    10. On page 12, line 214, what does "EH38E1530321" Stand for? It seems to specify a distal enhancer-like signatures in bipolar neurons, but I couldn't find this ID in any database.
    11. In the Methods section on lines 391-392, can the authors give some high-level explanation of the unconstrained integration method: "Single-nucleus RNA-seq of the same tissue and timepoints (GSE183684) were integrated using the unconstrained integration method". Also, can they comment on how retinal cell class identities were assigned (line 393). Was it based on known markers or on previous identification of cell classes and highly variable genes between clusters?
    12. In the integration of UCNEs with bulk and single cell ATAC-seq and DNase hypersentitivity regions, can the authors note in the Methods section (lines 400-404) what peak width was used to test for overlap with the UCNEs.
    13. On line 436, the word 'and' is missing between "(SNVs, and indels < 50bp)" and "large structural variants".
    14. On lines 443-444, please provide references to the computational tools listed. Please note if default settings/parameters were used.
    15. In the following sentence in the Methods section on lines 447-449, it is not clear in the Results section how this was used in the flow of the analysis, and how many cases showed such a similarity in phenotype: "For each candidate variant, we compared the similarities between the participant phenotype (HPO terms) and the ones known for its target gene through literature search and clinical assessment by the recruiting clinician when possible." Can the authors add more detail to the Results section.
    16. It would be helpful to have a table that describes the different omics datasets used in the paper, with some basic annotations (tissue type, sample size, reference)
    17. Can the authors add references to their sentence in the Discussion on page 17 lines 299-301: "As has been shown before, the phenotype caused by a coding mutation of a developmental gene can be different from the phenotype caused by a mutation in a CRE controlling spatiotemporal expression of this gene."

    Significance

    General assessment: This is the first study to our knowledge to integrate ultraconserved noncoding elements (UNCEs) with a range of multi-omic data to identify UNCEs that may be active and their target genes in a specific tissue, in this case human retina. They also experimentally test one of the UNCE predicted to affect the expression of a given gene in retina (and potentially other tissues) in Zebrafish and confirm its transcriptional activity in the eye, to provide some functional support to their strategy. This study also demonstrates the potential value of inspecting UCNEs in prioritizing pathogenic mutations for rare disease. One limitation is the lack of statistical significance assessment of the potential causal effect of a rare variant found in a UCNE on the expression of its predicted target gene, which is a rare eye disease gene, since the linkage of the UCNE to the disease gene was performed based on bioinformatic analysis of multi-omic data. Experimental testing of the effect of some of these mutations on their target gene expression could provide additional support.

    Advance: This study addresses an important unsolved problem in the field of human genetics and rare diseases, namely the challenge of identifying pathogenic mutations that lead to rare Mendelian diseases, in particular, in noncoding regions in unsolved cases. This is the first study to consider UCNEs together with tissue or cell type specific expression, epigenetics and chromatin accessibility in detecting pathogenic mutations for retina degeneration diseases. See for example Ellingford et al., Genome Medicine 2023 (PMID: 35850704). Their demonstration of rare variants in UCNE associated with inherited retinal degeneration diseases in patients with IRD in the 100k Genomics Project suggests that the role of UCNEs in disease should be further investigated and functionally tested. This work could have important clinical implications, and also proposes a strategy for integrating UCNEs with multi-omic and genomic data.

    Audience: This work will be of interest to the human genetics and genomics community, in particular to researchers interested in uncovering the genetic basis and causal mechanisms of rare diseases, and to scientists interested in clinical applications of genetic variation. This work will also be of interest to scientists interested more broadly in understanding the regulatory effects of the noncoding regions in the genome.

  6. Version published to 10.1101/2023.06.08.23291052v1 on medRxiv