Identification of cis-regulatory elements provides insights into tissue-specific gene regulation in the sheep genome

Domain: biology
Published: 2026-03-16
Editorial outcome: major_revision

Abstract
Annotation of regulatory elements is essential for understanding mechanisms underlying gene regulation, particularly tissue-specific regulation in human and animals. Here, we characterize 274,682 enhancers and 25,975 promoters across 24 tissues from an adult female sheep using ChIP-seq, ATAC-seq, CAGE-seq, RRBS, WGBS, and RNA-seq. We identify seven neural development-related genes with over 10 enhancers in brain tissues, highlighting the role of tissue-specific regulation. Cis-regulatory enhancer-promoter combinations provide insights into tissue-specific enhancers, such as the cerebellum-specific enhancer (chr15:57390520-57390685) regulating BDNF, which is expressed in both the cerebellum and cerebral cortex. Comparative analysis of enhancer-promoter combinations in human, mouse, pig, cattle, and sheep reveals ruminant-specific pathways, including pentose catabolism and long-chain fatty acid import regulation. A milk fat yield quantitative trait locus (QTL) identified within an enhancer interacts with the fat metabolism-related gene COMMD1, and a birth weight-associated QTL detected within a cerebellum-specific enhancer regulates XKR4. This study provides a robust framework for exploring cis-regulatory mechanisms and tissue-specific regulation, advancing the functional annotation of the sheep reference genome.

Introduction
The annotation of regulatory elements in the genome and assessment of their function in various tissues are essential for understanding gene regulation and interpreting the potential impact of genomic variants on phenotype. The encyclopedia of DNA elements (ENCODE) project annotated functional sequences on a genome-wide scale in human cells and tissue types1, and subsequently in model organisms2,3. These studies identified similarities and differences in regulatory elements across species. Functional annotation involves identifying and annotating the functional sequences in a genome, which include both expressed sequences (coding and non-coding RNAs) and regulatory sequences. Early research primarily focused on identifying expressed sequences using technologies like ESTs and cDNA libraries4, laying the groundwork for understanding gene expression profiles. The advent of RNA-seq then provided a more comprehensive view of the transcriptome, quantifying gene expression across tissues and developmental stages5,6, thereby enhancing the functional annotation of genomes. In contrast, identifying regulatory sequences on a genome-wide scale has been more limited, largely concentrated on humans and a few model organisms through projects like ENCODE1. These initiatives have created extensive datasets that illuminate the roles of regulatory elements, such as enhancers, promoters, and insulators, in shaping gene expression and driving complex phenotypes7. Both expressed and regulatory sequence annotation are paving the way for a more complete understanding of genome function.

Recently, significant progress has been made in the study of regulatory sequences in animal genomes, particularly in livestock species such as cattle8, pig9, and chicken10 through the project Functional Annotation of Animal Genomes (FAANG). In cattle, research identified associations between weak enhancers, active transcription start sites (TSS), and economic traits during rumen development8. In pigs, multi-omics data were used to annotate chromatin states across 14 tissues, revealing the importance of active promoters and enhancers in complex traits and domestication, with distinct regulatory patterns between Asian and European pigs9. Similarly, researchers annotated 1.57 million regulatory elements in chicken, predicting some enhancer-gene combinations and super-enhancers, providing valuable resources for poultry genetics and genomics research10. These projects not only deepen our understanding of genetic regulation but also fuel comparative and evolutionary studies, as well as enabling the investigation of the genetic control of phenotypes. As an important livestock species, the study of regulatory elements in sheep is equally valuable.

The ovine FAANG project offers a unique opportunity to characterize transcriptional regulatory elements in sheep, advancing our understanding of the biological mechanisms underlying complex traits in this species. Recently, we initiated the characterization of regulatory elements in three tissues (liver, spleen, and cerebellum) using chromatin immunoprecipitation sequencing (ChIP-seq) data to define chromatin states, including active and poised enhancers and promoters7. Other studies have identified regulatory elements in sheep macrophages, cartilage, cecum, pituitary gland, hypothalamus, liver, longissimus dorsi, rumen, spleen, and tail fat tissues11,12. However, current annotations remain limited to a few tissues, especially in the sheep genome. A comprehensive, systematic study identifying regulatory elements across multiple tissues is urgently needed, particularly focusing on enhancers, promoters, and their interactions.

Emerging research increasingly highlights the regulatory impact of enhancer-promoter interactions on gene expression, revealing their dynamic nature and marked tissue- and cell-type specificity13. A common approach for identifying enhancer regulations involves correlating chromatin signals at enhancers with gene transcriptional activity. Studies have utilized correlations between chromatin signals (H3K4me1 and H3K27ac) and gene expression to predict enhancer–promoter interactions. This correlation-based approach has proven an effective strategy for mapping regulatory interactions in research. For instance, Lin et al. identified putative enhancer-gene interactions within the same topologically associated domain based on significant positive correlations between H3K27ac signals at enhancers and gene expression14. Furthermore, transcription signatures, such as bidirectional and unidirectional transcription captured by cap analysis gene expression and deep sequencing (CAGE-seq), have become defining marks of active enhancers and promoters15. The FANTOM5 project utilized CAGE expression profiling to catalog enhancers across diverse human tissues and cell types, providing expression-based support for enhancer-promoter pairing that collectively modulate target gene expression16. Tissue-specific enhancer activity provides compelling examples of how these interactions underpin genomic regulation. Enhancer regions in the human prefrontal cortex selectively activate FOXP217, a gene crucial for neural development and vocal communication. In the mouse, muscle-specific enhancer regions marked by H3K27ac interact with MYF5 during myogenesis18, promoting muscle-specific expression. Similarly, in cattle, muscle-specific enhancers marked by H3K27ac target Myod, regulating muscle growth in skeletal muscle cells19. These reports highlight the pivotal role of enhancer-promoter or enhancer-gene interactions in shaping tissue-specific gene expression, offering valuable insights into the regulatory mechanisms underlying the complexity of genomic regulation.

To meet the need for a comprehensive analysis of regulatory elements across diverse sheep tissues and to elucidate tissue-specific characteristics of enhancer-promoter interactions, we systematically characterized enhancers and promoters across 24 tissues in this study using six high-resolution assays: chromatin immunoprecipitation sequencing (ChIP-seq), assay for transposase-accessible chromatin using sequencing (ATAC-seq), cap analysis gene expression and deep sequencing (CAGE-seq), reduced representation bisulfite sequencing (RRBS), whole-genome bisulfite sequencing (WGBS), and RNA-seq (Supplementary Fig. 1). All tissues were derived from the same individual whose genome was sequenced and assembled de novo, resulting in the current publicly available ovine reference genome ARS_UI_Ramb_v2.020. These datasets facilitated the systematic mapping of enhancers, promoters, and cis-enhancer-promoter combinations throughout the sheep genome. This extensive analysis provides critical insights into tissue-specific gene regulation and the genetic architecture of complex traits in this economically significant species.

Results
Overview of collected data
Twenty-four tissue samples from the animal used to generate the ARS_UI_Rambv2.0 reference genome was collected in this study (Fig. 1A). These samples were used to explore chromatin states and to identify and characterize key regulatory elements, such as enhancers and promoters, through a comprehensive set of 188 high-throughput sequencing datasets (Fig. 1B and Supplementary Data 1). The datasets include chromatin immunoprecipitation sequencing (ChIP-seq) for four histone marks (H3K27ac, H3K27me3, H3K4me1, and H3K4me3) along with control samples, ATAC-seq, CAGE-seq, and RNA-seq. Each contributes a distinct insight into the intricate landscape of regulatory elements for gene regulation. A total of 12.8 billion raw reads were sequenced from 188 sequencing datasets, generating 10.4 billion mapped reads, with an average mapping rate of 80.64% for all datasets (Supplementary Data 1). The average number of reads was 49,493,667; 48,223,579; 57,509.914; 63,299,552; 50,676,069; 4,321,867; 102,757,218; and 72,333,332 for H3K27ac, H3K4me3, H3K4me1, H3K27me3, Control, CAGE, ATAC and RNA-seq data, respectively (Supplementary Data 1).

Characterization of chromatin states
ChIP-seq libraries of H3K27ac, H3K4me3, H3K4me1, H3K27me3, and Control were generated from 24 tissues and sequenced (see Materials and Methods). A total of 8,802,602 peaks were identified in all four histone marks across 24 tissues and annotated on the peaks to determine their genomic features relative to the reference genome’s annotated elements (Figs. 1C and 2 and Supplementary Table 1). The majority of H3K4me1 and H3K27ac peaks were mapped to intron and intergenic regions across all tissues (Fig. 1C and Supplementary Fig. 2A, B), supporting their roles as key indicators of enhancers. Approximately half of the peaks for these two histone marks were located within intron regions (50.99% ± 5.68 for H3K4me1 and 45.07% ± 6.48 for H3K27ac), while peaks near the transcription start site (TSS) within 2 kb were more abundant for H3K27ac (19.37% ± 5.83) than for H3K4me1 (11.77% ± 3.41). In contrast, for H3K4me3, the majority of peaks were found near TSS regions (34.4% ± 11.99), more than in intergenic and intron (Supplementary Fig. 2C), confirming H3K4me3 and H3K27ac as characteristic marks of promoters. H3K27me3 peaks were predominantly located in regions without gene annotation (64.26% ± 2.56), with relatively few mapped near annotated TSS (5.72% ± 0.92) (Supplementary Fig. 2D). Overall, the peak localizations for all four histone marks relative to gene annotations were consistent with known properties of these marks with respect to genes and gene activity, supporting the reliability of the data for annotating promoters and enhancers and their status relative to gene activity.

Chromatin states were further annotated by ChromHMM21, and the ten specific chromatin states were identified by using four histone marks and ATAC-seq datasets (Fig. 1D). The three states (TssA, TssAFlnk and BivFlnk) involving TSS covered ~1.6% of the genome. The active and poised enhancer states, defined by H3K27ac, H3K4me1, and ATAC signals, covered ~17.6%. Repressed polycomb and weak repressed polycomb states defined by H3K27me3 combined to cover ~26.6% (Fig. 1E). These chromatin states represented ~45.9% of the genome with distinguishable modification signals measured with these histones (Fig. 1E and Supplementary Table 2). The chromatin states from sheep were consistent with a previous study of cattle8 and chicken10 chromatin states.

Identification of enhancers and promoters across tissues
To further clarify the specific localization of chromatin states and identify key regulatory elements, particularly enhancers and promoters across tissues, we integrated peaks from CAGE-seq and ATAC-seq data to cross-validate those identified from ChIP-seq (Fig. 2A). By cross-validation of ChIP-seq, CAGE-seq, and ATAC-seq datasets, we identified a total of 424,559 peaks as potential enhancers and 320,425 peaks as potential promoters from 20 tissues, with four tissues excluded due to missing ATAC-seq data (Fig. 2A and Supplementary Fig. 3).

Based on the variation in the peak number of shared peaks across tissues (Supplementary Fig. 3), we further identified the high-confidence peaks identified in at least two tissues for both enhancers and promoters. This resulted in a total of 119,186 and 89,433 peaks, respectively (Fig. 2A and Supplementary Fig. 3), which were subsequently used to infer potential enhancers and promoters specific to each of the 24 tissues. Intersecting the peaks with H3K4me1, H3K27ac, and CAGE data yielded an average of 71,340 potential enhancers per tissue, out of a total of 274,682 regions across all tissues combined (Fig. 2B, Supplementary Fig. 4A, and Supplementary Data 2 and 3). The width of these potential enhancers ranged from 50bp to 1000bp (Supplementary Fig. 4B). For identification of promoter, TSSs identified from RNA-seq and CAGE-seq data, capturing both annotated and unannotated TSSs, were used to facilitate the mapping of promoters (Fig. 2A). A total of 33,801 annotated TSS were derived from the GTF file, and additional 331 unannotated TSS were identified through integrated analysis of RNA-seq and CAGE data (Fig. 2A and Supplementary Fig. 5A). After intersecting the 89,433 peaks with the TSS regions, a total of 25,975 promoters were identified across tissues, comprising 271 unannotated promoters and 25,704 annotated promoters (Fig. 2B, Supplementary Fig. 5B, and Data 2 and 4). The majority of promoters exhibited widths within a 0–500 bp range, with distances to the TSS similarly falling within 0–500 bp (Supplementary Fig. 5C, D).

In the tissues examined, the adrenal cortex exhibited the greatest number of enhancers (135,371) and promoters (15,170), with enhancer variation across tissues notably exceeding that of promoters, while the cerebellum contained a relatively large number of tissue-specific enhancers (4,320) and promoters (323), particularly tissue-specific enhancers, which were markedly enriched compared with other tissues (Fig. 2B and Supplementary Table 3). Furthermore, the distribution patterns of shared enhancers and promoters across tissues differed markedly (Fig. 2C-E and Supplementary Fig. 4C). The shared enhancers numbers declined as the number of tissues increase (Fig. 2C), while promoters show a bimodal distribution, with both tissue-specific and broadly shared promoters more prevalent (Fig. 2D). Across all tissue pairs, the pairwise consistency fraction of promoters exceeded 0.7, substantially higher than that observed for enhancers (0.1) (Fig. 2E and Supplementary Fig. 4C). These findings suggest that promoters are less tissue-specific than enhancers, indicating that different tissues may have unique enhancers, which was consistent with the observation that enhancer variability across tissues is greater than that of promoters, and also supports the relatively high average of enhancers per promoter (Fig. 2B and Supplementary Table 3). Additionally, 271 unannotated promoters were identified across 24 tissues, exhibiting a similar bimodal distribution (Supplementary Fig. 5E and Supplementary Data 4). Interestingly, among these, nine tissue-specific TSS were aligned with annotated gene sequences, of which five were annotated in other species (Bos Taurus, Capra hircus, and Cervus elaphus), implying that they may correspond to unannotated genes in sheep (Supplementary Data 5).

Cis-enhancer and promoter combinations across tissues
To further explore the relationship between gene enhancers and promoters, we identified enhancer-promoter combinations based on correlations of gene expression and enhancer activity. Out of 25,975 genes with identified promoters, 21,507 of them were expressed, and 4468 were not expressed in any tested tissue. We identified a total of 10,197 combined cis-enhancer-promoter regulatory elements involving 4293 genes and 8860 enhancers across 24 tissues (Supplementary Data 6). Among these, 2607 were supported by the identified cis-enhancer and promoter combinations from CTCF or CAGE data (see M&M), involving 1180 genes (Supplementary Data 6). For expressed genes, gene expression abundance (Transcripts per million, TPM) differed significantly across tissues with enhancers compared to those without, with higher expression of genes observed in tissues with enhancers than in those without (Fig. 3A). Furthermore, in tissues where a cis-enhancer was present (colored in Fig. 3B), the corresponding genes were predominantly expressed, whereas non-expressed genes consistently lacked an associated enhancer in these tissues (Fig. 3B and Supplementary Fig. 6). For example, the gene ADORA1 on chromosome 12, which encodes the adenosine A1 receptor with important functions in the brain. Two enhancers Enh1 (Chr12: 682,652-683,139) and Enh2 (Chr12: 737,891-738,081) were paired with ADORA1, both showing a high correlation with gene expression across tissues. This gene was expressed in three brain-related tissues (cerebellum, cerebral cortex, and spinal cord). Enh2 was consistently identified in these tissues (correlation=0.998, R²=0.996). Interestingly, Enh1 associated with ADORA1 was also detected in the bladder (correlation=0.844, R² =0.698) (Fig. 3C). This aligns with previous research highlighting the role of luminal adenosine turnover and A1 receptor activation in regulating bladder function22, providing additional validation for our analysis.

Promoter and enhancer hubs reveal two predominant interaction patterns in gene transcription regulation, where enhancer hubs involve multiple promoters connected to a single enhancer, and promoter hubs involve multiple enhancers connected to a single promoter. In both hubs, the most common interaction type is a single enhancer paired with a single promoter (Supplementary Fig. 7). Approximately 52% of promoters were associated with a single enhancer (Supplementary Fig. 7A), while 89% of enhancers associated with a unique promoter (Supplementary Fig. 7B). These interaction patterns observed across all tissues were consistent with the trends in each tissue-specific hub (Fig. 3D, E and Supplementary Data 7). Notably, the vast majority of genes (99%) in promoter hubs interacted with fewer than ten enhancers (Fig. 3D), with most enhancers regulating only one gene, a smaller subset regulating two to three genes, and very few regulating more than three (Fig. 3E). Interestingly, brain and adrenal-related tissues exhibited similar hub patterns within each cluster (Fig. 3F), with more than 10 enhancers consistently identified for seven genes (DNAJC6, FMNL2, NAV3, AMER2, NCAM1, OTUD7A, and LRRTM3) in both brain-related tissues (cerebellum and cerebral cortex). These genes play critical roles in neural development and function: DNAJC6 encodes an HSP40 chaperone protein that regulates ATPase activity and metabolism; FMNL2 governs cell migration, size, and shape changes; NAV3 is essential for normal nerve development, with both over-expression and under-expression linked to neural deficits; AMER2 functions as a tumor suppressor via the WNT pathway; NCAM1 is involved in cell adhesion, development, and synaptic regulation; OTUD7A is crucial for neural development; and LRRTM3 likely supports normal neuronal function and has been implicated in Alzheimer’s disease. The presence of multiple enhancers regulating these genes highlights their critical importance in maintaining proper neural development and function. This intricate modulation by enhancers ensures precise gene expression during critical windows of development while minimizing the risk of dysregulation-associated pathologies. The observed greater diversity and variability in enhancers compared to promoters aligns with known biological processes, reflecting the dynamic and tissue-specific roles of enhancers in gene regulation.

Tissue-specific enhancer-promoter combinations and functional annotation
To investigate the functional roles of enhancers and promoters, we identified and compared the genes linked to tissue-specific promoters and enhancers from enhancer-promoter combinations. Notably, the genes associated with tissue-specific promoters were all included in the genes of enhancers from tissue-specific enhancer-promoter combinations (Fig. 4A and Supplementary Data 8). The tissue-specific enhancer-promoter combinations provide additional insights, as they incorporate interactions of enhancers and promoters that are not captured when considering tissue-specific gene expression alone. While previous analyses have focused on promoter-gene specificity derived from gene annotations, the inclusion of the unrecognized regulatory roles of enhancers from enhancer-promoter combinations uncovers insights into gene regulation.

To validate the contribution of enhancers from combinations, we performed GO enrichment analysis on genes associated with tissue-specific enhancer genes and compared with these to promoter-associated genes. GO terms reflected the biological functions of the respective tissues, with particularly strong associations observed in brain and heart tissues (Fig. 4B and Supplementary Data 9). Among the 24 tissues analyzed, genes associated with both enhancers and promoters displayed overlapping GO terms in 14 tissues. However, in the remaining ten tissues, GO analysis using promoter-enhancer combinations yielded stronger associations, including lower p values, a greater number of GO terms, and higher tissue correlation, compared to results obtained using promoter-only information (Supplementary Data 10). Furthermore, we performed GO analysis specifically for the 151 genes uniquely associated with enhancers in the cerebellum (Fig. 4A), five significant biological processes were identified and related to synaptic transmission and neuronal differentiation (Supplementary Table 4). Notably, 11 key genes (e.g., BDNF, NLGN3, and NRCAM) were shared across multiple pathways (Supplementary Fig. 8), highlighting their central roles in neuronal development and synaptic function, and reinforcing the biological relevance of the enhancer-specific gene set.

Moreover, the cerebellum-specific enhancer-promoter combination involving the BDNF gene further exemplifies this regulatory relationship (Fig. 4C and Supplementary Data 11). Although BDNF is expressed in both the cerebellum and cerebral cortex, its cerebellum-specific expression was not detectable using promoter-gene specificity alone. By incorporating enhancer data, we identified a cerebellum-specific enhancer (Enh1, chr15: 57,390,520-57,390,685) active exclusively in the cerebellum, whereas a second enhancer (Enh2, chr15: 57,150,237-57,150,878) was shared between the cerebellum and cerebral cortex. This discovery yields a reasonable hypothesis that the cerebellum-specific enhancer may be the key to tissue-specific BDNF expression relevant to the distinct cerebellar phenotypes when compared to the cerebrum. Differential H3K4me1 signals in Integrative Genomics Viewer (IGV) from ChIP-seq data further validated tissue-specific enhancer activity (Fig. 4C). RNA expression of exon and read junctions from RNA-seq revealed three BDNF isoforms in the cerebellum, of which two were exclusive to this tissue (XM_015101092.3 and XM_042233303.1), compared to one isoform in the cortex. These findings suggest a regulatory association between enhancer Enh1 (chr15: 57,390,520-57,390,685) and cerebellum-specific isoforms of BDNF (Fig. 4C).

To further elucidate enhancer function, we analyzed tissue-specific enhancer sequences to identify associated transcription factors (TFs). A total of 81 TFs were identified (Fig. 4D and Supplementary Data 12), many of which are known regulators of the gene expression of corresponding tissues. For instance, Nkx2.2, identified in the cerebral cortex, is critical for oligodendrocyte differentiation and CNS development, and Atoh1 from cerebellum tissue plays a pivotal role in neurogenesis. Forkhead box (Fox), prominently associated with cardiac tissues, cardiovascular development, cell differentiation, and stress response. Similarly, SOX15 from SOX family TFs contributes to the development and structural integrity of cardiac tissue (Fig. 4D). These results highlight the importance of enhancers in tissue-specific regulation and emphasize their role in defining enhancer-promoter interactions. These discoveries establish opportunities for continued hypothesis-based research that will elucidate the nuances of tissue-specific phenotypes established through tissue-specific gene expression regulation.

To account for organ-level specificity, we extended our analysis to identify enhancers associated with broadly tissue-specific combinations within key organs, including the brain, heart, adrenal gland, stomach, and colon. By analyzing tissues derived from the same organ, we aimed to explore how different tissues within an organ utilized distinct enhancers to regulate the expression of shared genes, thereby highlighting tissue-specific regulatory mechanisms within an organ. This expanded analysis revealed that although tissues from the same organ often share target genes regulated by enhancers, different tissues frequently employ distinct enhancers for these genes, reflecting tissue-specific regulatory diversity within the organ (Fig. 4E and Supplementary Data 13). For instance, in the brain, enhancers specific to the cerebellum and cerebral cortex were found to regulate common neural genes, such as BDNF, but the specific enhancers used differed between these two tissues, as previously discussed. Similarly, in both the stomach and colon, distinct enhancers targeting shared genes involved in epithelial barrier maintenance and digestive enzyme production emphasized tissue-level specialization within the gastrointestinal system. These findings demonstrate that, even within the same organ, the regulatory landscape is finely tuned to the unique functional requirements of individual tissues and likely at very specific timing and life stages. The identification of such enhancer combinations not only underscore the complexity of gene regulation at the organ level but also provides valuable insights into the mechanisms driving tissue specificity within broader organ systems. This enhanced understanding of intra-organ regulatory dynamics offers a critical foundation for future studies on tissue- and organ-specific gene regulation.

Tissue-specific enhancer-promoter combination and methylation
To explore the specific impact of enhancer-promoter combinations on gene expression across tissues, we collected DNA methylation data from whole-genome bisulfite sequencing (WGBS) and reduced representation bisulfite sequencing (RRBS) across 22 tissues (Fig. 1B) and analyzed the association between DNA methylation within enhancer-promoter combinations and the RNA expression of related genes. Results from both WGBS (all five tissues) and RRBS (21 tissues) revealed that methylation levels in enhancer regions were negatively correlated with the expression of their cis-regulatory genes (Fig. 5A and Supplementary Figs. 9A and 10), highlighting the distinct regulatory roles of enhancers in gene expression. For promoter regions, RRBS methylation levels did not exhibit a clear correlation with gene expression (Supplementary Fig. 11). However, methylation patterns derived from WGBS (five tissues) were consistently negative (Fig. 5B and Supplementary Fig. 9B), aligning with previous studies reporting that promoter methylation decreases gene expression23. The findings indicated that the repressive effect of enhancer methylation on gene expression is consistent with that of promoter methylation. This study demonstrates the suppressive impact of enhancer methylation on gene expression in sheep, which aligns with known biological mechanisms. Methylation likely significantly influences the interaction between enhancers and promoters, thus regulating gene activation. This commonality in methylation effects on both enhancers and promoters is biologically plausible and aligns with findings from previous studies24, highlighting the role of methylation in impairing enhancer-promoter interactions. Furthermore, RRBS data revealed that enhancer methylation exhibits a stronger repressive effect on gene expression than promoter methylation, consistent with findings in human cancer studies25.

Cross-validation of cis-enhancer-promoter combinations across species
To further substantiate the functionality of the identified enhancers and promoters in sheep, we conducted a comparative analysis of four selected tissues (cerebellum, cortex, lung, and muscle) across five species (Fig. 6A). Relative to human, we identified an average of 8914 promoters in sheep tissues, compared to 7736 in cattle, 7963 in pig, 5313 in mouse, and 21,792 in human. Enhancers were identified more than promoters, averaging 42,904 in sheep, 17,808 in cattle, 18,770 in pigs, 123,066 in mouse, and 280,892 in humans (Supplementary Data 14).

Interspecies comparison of enhancers revealed that 61%, 68%, 68%, and 72% of sheep enhancers in cerebellum (69,399), cortex (42,467), lung (46,371), and muscle (13,378), respectively, were conserved in at least one other species. (Fig. 6B and Supplementary Fig. 12A-C). Promoters exhibited even higher conservation rates across species, with 86, 89, 88, and 92% of promoters in the same tissues conserved in other species (Fig. 6B and Supplementary Fig. 12D-F). Moreover, on average, approximately 17% of promoters (17.4% in cerebellum, 17.7% in cortex, 16.5% in lung, and 18.3% in muscle) were conserved across all five species, substantially higher than enhancers 0.7% (0.6, 0.6, 0.7, and 1% in the respective tissues) (Fig. 6B and Supplementary Fig. 12). Additionally, sheep exhibited more conservation with cattle than pig, reflecting the closer evolutionary relationship between these ruminant species. (Fig. 6B and Supplementary Fig. 12). Notably, enhancers showed higher overlap observed for human and mouse may result from their more comprehensive annotation and the larger catalog of regulatory elements available for these species. Promoters displayed consistently higher conservation rates across all five species compared to enhancers, highlighting a fundamental difference in conservation patterns between the two types of regulatory elements (Fig. 6C, Supplementary Fig. 13, and Supplementary Data 15), which is consistent with previous observations in mammalian species26. These cross-species comparisons underscore the greater variability of enhancers relative to promoters, emphasizing the enhancer regulatory specificity and functional diversity in sheep.

We further focused on tissue-specific cis-enhancer-promoter combinations, we identified 596, 107, 117, and 23 unique combinations in cerebellum, cortex, lung, and muscle, respectively (Fig. 6D). Functional enrichment analysis of these cis-regulatory elements revealed strong associations with tissue-specific biological functions: cerebellum-related combinations were significantly enriched in synaptic organization, cortex combinations in MAPK cascade regulation, lung combinations in smooth muscle apoptosis regulation, and muscle combinations in G protein-coupled opioid receptor signaling pathways (Fig. 6E and Supplementary Data 16). Furthermore, the cerebellum and cortex tissues exhibited the highest gene overlap (178 genes), with functional terms pointing to brain-specific roles, including nervous system development and synaptic transmission (Fig. 6E and Supplementary Data 16). Additionally, 69 gene-associated enhancer-promoter combinations conserved exclusively in ruminants (sheep and cattle) (Supplementary Data 17), or for which the GO enrichment was relevant to ruminant physiology, such as pentose catabolic processes and regulation of long-chain fatty acid import (Fig. 6E and Supplementary Data 16). These findings validate the robustness of the identified enhancers and cis-enhancer-promoter combinations, which effectively reflect tissue-specific regulatory functions and species-specific biological processes.

Biological insights of cis-regulatory elements with genetic variation in sheep
A total of 15,816,440 SNPs were identified from whole-genome sequencing data of 355 sheep, which were downloaded from the iSheep database (https://ngdc.cncb.ac.cn/isheep/download), of which 8% were located within exon, promoter, and enhancer regions (Fig. 7A). The number of SNPs in intergenic (without enhancer and promoter regions) and introns were significantly higher than in regulatory elements (promoters and enhancers) and exon regions, and SNP density was higher in enhancers and promoters compared to exon regions (Fig. 7B). To further explore the biological functions of SNPs within enhancer and promoter regions, we retrieved sheep QTL from the sheep QTLdb (https://www.animalgenome.org/QTLdb). A total of 38 QTLs colocalized with 17 genes in promoter regions and 215 QTLs overlapping 214 enhancers were identified (Supplementary Data 18). Among these, one QTL for milk fat yield colocalized with an enhancer interacting with the COMMD1 gene. COMMD1 is expressed in at least 21 tissues, with the same enhancer consistently identified across all of them (Fig. 7C). The correlation between this enhancer and COMMD1 gene expression was strong (r=0.90, P=1.17 x 10-8), suggesting a potential regulatory relationship. COMMD1 is a compelling candidate gene given its established role in copper metabolism, which has been linked to lipid metabolism and milk fat composition in ruminants27,28.

To further assess the reliability and functional relevance of the identified cis-regulatory elements associated genes, we analyzed 23,056,026 SNPs from a recent genome-wide association study (GWAS) of birth weight (BW) in 3007 sheep29. Among these, 28,736 SNPs were significantly associated with BW (P≤0.001). Intersecting these variants with our enhancer–promoter interactions across 24 tissues identified 18 SNPs located within 17 enhancers active in the cerebellum and/or cerebral cortex (Supplementary Data 19). Notably, seven cerebellum-specific enhancers were predicted to regulate genes including XKR4, GALNT9, RET, TRNAC-ACA-97, TRNAG-CCC-197, TRNAC-GCA-185, LOC114110132, and LOC121817986. Of particular interest, the most significant SNP (P=2.23×10-7), located within a cerebellum-specific enhancer, is predicted to regulate XKR4. In the same GWAS study29, XKR4 had previously been associated with BW based on an intragenic SNP. Our analysis demonstrates a potential regulatory mechanism by which this, and an additional variant within a regulatory enhancer controlling XKR4, may exert regulatory control on XKR4, thereby providing additional evidence that variation within regulatory elements can influence phenotypic traits in sheep. In addition, GALNT9 encodes a brain-specific glycosyltransferase with high expression in the human cerebellum30, and RET is known to regulate cell growth, differentiation, and nervous system development31. Further, we identified another significant SNP (P=7.47×10-5) within an enhancer active in the cerebellum and spinal cord that interacts with TSPAN5 (Supplementary Data 19), a gene highly expressed in the brain during postnatal development in mice32.

To obtain more direct evidence of enhancer-associated gene expression regulation, we analyzed 181 published significant eQTLs from the same study, we identified one variant significantly associated with LOC121819880 expression (P=4.69×10⁻⁸), which resides within its interacting enhancer (chr6:35907531-35908224) (Supplementary Data 20). This observation provides additional evidence that sequence variants within enhancers can directly influence gene expression, reinforcing the accuracy and biological validity of our enhancer annotations. Collectively, these findings demonstrate how incorporating enhancer information enhances the biological interpretation of gene regulation in sheep and provides insights into ruminant-specific functional mechanisms. More functional validation is needed to confirm the causal roles of these enhancer-associated variants.

Methods
Tissue collection
Tissues were collected from a healthy 6-year-old Rambouillet ewe (Benz 2616, USMARC ID 200935900), whose genome was sequenced and assembled de novo to produce the ovine reference genome ARS_UI_Ramb_v2.020. Ethics approval for all animal procedures were obtained by the Utah State University Institutional Animal Care and Use Committee on April 29, 2016, at the Utah Veterinary Diagnostic Laboratory. Postmortem tissue samples were immediately transferred to cryogenic vials and snap frozen in liquid nitrogen, with subsequent transfer to −80°C freezers for storage. Further details regarding the tissue collection strategy are available through the FAANG Data Coordination Center’s data portal (https://data.faang.org/home), and the metadata are recorded according to FAANG guidelines under BioSample accession SAMEG32960715,43.

Chromatin immunoprecipitation with sequencing (ChIP-seq)
Chromatin immunoprecipitation (ChIP) was performed for 24 collected tissues from Benz 2616 (Supplementary Table 1). The SimpleChIP Plus Enzymatic Chromatin IP Kit (Cell Signaling Technologies, cat. #9005, Danvers, MA, USA) was used for immunoprecipitation according to the manufacturer’s instructions44,45. Antibodies specific for histone modifications H3K4me3 (Abcam, cat. #ab8580), H3K27ac (Abcam, cat. #ab4729), H3K4me1 (Abcam, cat. #ab8895), and H3K27me3 (Abcam, cat. #ab6002) were used for immunoprecipitation. Briefly, tissue was crosslinked with 37% formaldehyde, disaggregated with Dounce homogenization, and chromatin was sheared with micrococcal nuclease (MNase) by incubating on a thermomixer for 20 min at 200 rpm and 37°C. Following nuclear membrane lysis and isolation of sheared chromatin by centrifugation at 15,000×g for 1 min, chromatin was incubated with 1µg of antibody overnight (16 h) at 4°C on a HulaMixer (Thermo Fisher Scientific, USA). An input control (20µl from each immunoprecipitation) sample was removed from sheared chromatin prior to incubation with antibodies and stored at −20°C until purification. Protein-G-coated magnetic beads were used to isolate chromatin bound to specific antibodies, and purification of enriched chromatin was performed with the DNA Purification Buffers and Spin Columns Kit following manufacturer’s instructions (Cell Signaling Technologies, cat. #14209, Danvers, MA, USA).

Purified DNA was quantified using a Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific, cat. #Q32854, Waltham, MA, USA), and size and integrity was verified with a Fragment Analyzer (Agilent, Santa Clara, CA, USA). Libraries were prepared as 75 bp paired end reads with the TruSeq ChIP Library Preparation Kit (Illumina, Inc., cat. #IP-202-1012, San Diego, CA, USA) following manufacturer’s instructions. Libraries for H3K4me3, H3K27ac, and Input were sequenced to achieve at least 20 million uniquely mapped reads, and H3K4me1 and H3K27me3 libraries were sequenced to achieve at least 45 million uniquely mapped reads to comply with ENCODE standards for ChIP-seq1.

Quality control of sequence data was performed with FastQC v0.11.9 and MultiQC v1.1446, and the peaks of H3K4me3, H3K27ac, H3K4me1, and H3K27me3 were called by the chipseq model of Nextflow v24.04.247 with Bowtie2 aligner v2.4.2 and MACS2 peak caller v2.2.7.1 with “-macs_pavalue 0.005”48,49, for narrow histone marks H3K4me3 and H3K27ac with additional parameter “--narrow_peak”. The reference genome for alignment was from NCBI (ARS-UI_Ramb_v2.0, GenBank accession GCF_016772045.1)20, which was also used for all other assays, and the annotation was generated by the standard NCBI annotation pipeline v9.0. The genomic distribution of H3K4me1, H3K27me3, H3K4me3, and H3K27ac peaks across feature regions, including promoters (TSS±2kb), exons, introns, and intergenic regions, was analyzed using the R packages GenomicFeatures v1.50.450 and ChIPseeker v1.34.151.

Cap analysis gene expression sequencing (CAGE-seq)
CAGE-seq data from 24 tissues was generated from ref. 15, and the libraries were prepared for each sample from a starting quantity of 5 μg of DNase-treated total RNA. Random primers were used to ensure conversion of all 5’-cap-trapping RNAs according to ref. 15. The full protocol is available via the FAANG Data Coordination Centre (https://data.faang.org/protocol/experiments/cageProtocol-CAGE-seq-TSS). CAGE-seq data were used to identify the transcription start site (TSS) for each gene. The analysis pipeline and metrics are available on Bitbucket (https://bitbucket.org/msalavat/cagewrap_public/src/master). Briefly, raw CAGE reads were trimmed using tagDust252 and mapped to the ARS-UI Ramb v2.0 using bowtie2 v2.4.2 (--very-sensitive mode)49. The mapped BAM files were converted to strand-specific bedGraphs using BEDTools v2.29.053 and subsequently converted to bigWig format using the UCSC BedGraphToBigWig tool54. The clustering analysis was carried out in R using CAGEfightR v1.12.055, applying “unidirectional” and “bidirectional” clustering algorithms to identify putative TSSs and TSS-Enhancer regions, respectively. The interaction between TSSs and enhancers were further determined by the findLinks function in CAGEfightR v.1.12.055.

Assay for transposase-accessible chromatin with sequencing (ATAC-seq)
The assay for transposase accessible chromatin (ATAC) was performed on 20 tissues as a service provided by the University of California, San Diego Center for Epigenomics. Libraries were sequenced to a depth exceeding 90 million reads (Supplementary Data 1). Quality control and peaks calling parameters for ATAC-seq were consistent with those used for ChIP-seq, utilizing the atacseq model of Nextflow V24.04.247.

RNA sequencing
RNA was isolated from snap-frozen tissue by cryopulverizing tissue and placing it in 2 mL of Trizol. The tissue was then homogenized in the Trizol and split into 1 mL aliquots in 1.5 mL microcentrifuge. Samples were then centrifuged for 5 min at 14,000×g and 4°C. The aqueous (clear) layer was transferred to a new tube, avoiding any fat that may have risen to the top of the liquid. Then, 400 µL of chloroform was added to each tube, shaken vigorously, and incubated at room temperature for 15 min. The sample was centrifuged for 15 min and 14,000×g and 4°C, followed by transfer of the aqueous layer into another 1.5 mL microcentrifuge tube. About 500 µL of isopropanol was then added to each tube and incubated for 10 min at room temperature, followed by centrifugation at 14,000×g and 4°C for 30 min. The supernatant was discarded, and the pellet was washed with 75% ethanol. The purified RNA was resuspended in nuclease-free water and quantified with a Nanodrop, as well as quality evaluated with a Fragment Analyzer. RNA was then treated with DNase I (1 µL of 2000 U/mL per 10 µg RNA) and cleaned up with the Rneasy MinElute spin column kit (Qiagen).

RNA-seq libraries were prepared using the Illumina TruSeq Stranded mRNA Library Preparation Kit (Illumina, Inc.) according to the manufacturer’s instructions. Sequencing was performed on an Illumina NextSeq 500, generating at least 65 million reads per sample (Supplementary Data 1). This sequencing depth ensures that fragments of most expressed genes are covered at least 30-fold, providing sufficient confidence for expression quantification. The reads were aligned to the ARS-UI Ramb v2.0 reference genome using STAR v2.7.8a56. Gene expression levels were quantified as transcripts per million (TPM) using StringTie v2.2.057, and genes with TPM ≥ 0.5 were classified as expressed only for the quantitative analysis of enhancer-promoter expression correlations, while no expression filtering was applied in subsequent qualitative comparisons.

Chromatin state characterization
Chromatin states for the 24 tissues were characterized using the identified ChIP-seq and ATAC-seq peaks and analyzed with the Hidden Markov Model in ChromHMM21, which integrates overlapping histone modification signals and chromatin accessibility to define chromatin states, including active, poised, repressed promoter and enhancer and repressor. The optimal number of chromatin states was examined by implementing ChromHMM with 2 through 16 chromatin states in the LearnModel function, followed by the median Pearson correlation output by the CompareModels function, where each tissue plateaued and was tightly correlated with the model with the greatest number of states37,58. The percentage of the genome occupied by each chromatin state was calculated and averaged across tissues, providing a comprehensive view of chromatin landscape variability.

Identification of enhancers and promoters
Histone modifications have been widely used to characterize enhancers and promoters10,41. Specifically, H3K4me3 and H3K27ac are commonly associated with promoters, while H3K4me1 and H3K27ac indicate enhancers. To achieve a more reliable and accurate identification of enhancers and promoters, we developed a comprehensive strategy integrating four distinct high-throughput omics datasets: ChIP-seq, CAGE-seq, ATAC-seq, and RNA-seq. The 5’ end of the TSS and enhancer RNA were characterized by unidirectional and bidirectional sequences (CAGE-unid and CAGE-id) from the CAGE-seq data, respectively. Chromatin accessibility was mapped using ATAC-seq data. The TSS and active gene expression were profiled by RNA-seq (Fig. 2A).

The specific processes of this strategy are:
I. Definition of candidate enhancer or promoter regions for each tissue:
Peak_enh = P_H3K4me1 ∪ P_H3K27ac (1)
Peak_prom = P_H3K4me3 ∪ P_H3K27ac (2)
Where Peak_enh and Peak_prom are candidate peak regions for enhancer and promoter, respectively, and P. represents peaks from the respective histone marks from ChIP-seq data.
II. Identification of confident potential enhancers or promoters:
To refine candidate regions, additional datasets were integrated. For each tissue (T_i), the potential enhancers and promoters were defined as follows:
P_enh_Ti = (Peak_enh_Ti ∩ P_CAGE_bid_Ti) ∪ (P_H3K4me1_Ti ∩ P_H3K27ac_Ti) ∪ (P_H3K4me1_Ti ∩ P_ATAC_Ti) (3)
P_prom_Ti = (Peak_prom_Ti ∩ P_CAGE_unid_Ti) ∪ (P_H3K4me3_Ti ∩ P_H3K27ac_Ti) ∪ (P_H3K4me3_Ti ∩ P_ATAC_Ti) (4)
Given ATAC-seq data availability for only 20 tissues, i = 1, 2, …, 20. The union of confident peaks across all tissues yielded:
P_enh_all = ∪_{i=0}^{24} P_enh_Ti, P_prom_all = ∪_{i=0}^{24} P_prom_Ti (5)
We used P_enh_all and P_prom_all as the potential enhancer and promoter regions for filtering across all 24 tissues. For a given i-th tissue (T_i), the any potential enhancer (Enh_ij) and promoter (Prom_ij) regions can be represented as follows:
Enh_ij ⊂ P_enh_all, Prom_ij ⊂ P_prom_all (6)
Where j is the j-th enhancer or promoter from tissue i (T_i), respectively. To refine potential enhancer and promoter regions, we applied two criteria:
1 The region must be identified in at least two tissues (∑_{i=1}^{24} t(Enh_ij) ≥ 2), where t(.) is an indicator function:
t(x) = 1 if x ⊂ P_enh_all or P_prom_all, 0 otherwise (7)
2 The region must have a p value (Pval_H3K4me1 for enhancers, Pval_H3K4me3 for promoters) below 0.001.
The final confident enhancer and promoter regions for each tissue were defined as:
Enh_Ti = {Enh_ij ⊂ P_enh_all | ∑_{i=1}^{24} t(Enh_ij) ≥ 2 or Pval_H3K4me1_ij ≤ 0.001} (8)
Prom_Ti = {Prom_ij ⊂ P_prom_all | ∑_{i=1}^{24} t(Prom_ij) ≥ 2 or Pval_H3K4me3_ij ≤ 0.001} (9)
Taken together, the potential enhancer and promoter regions were defined using integrative evidence from histone modification profiles (H3K4me1, H3K27ac, and H3K4me3), CAGE bidirectional and unidirectional transcription signals (CAGE_bid and CAGE_unid), and ATAC-seq chromatin accessibility data (ATAC)(Supplementary Table 6). Promoters were further refined using transcription start site (TSS) information (±2kb of TSS). TSS were derived from two sources: (1) annotated genes in the reference genome, and (2) unannotated genes from RNA-seq and CAGE-seq. identified unannotated TSSs were extracted from StringTie-assembled transcripts with TPM ≥ 10, overlapping with CAGE-seq TSS regions, and detected in at least two tissues.
III. Consistent genomic coordinates for enhancers and promoters:
To unify the genomic coordinates of enhancers and promoters across tissues, overlapping regions in Enh_Ti or Prom_Ti from multiple tissues were consolidated using a stepwise median approach. For k overlapping tissues (k ≤ 24), with start positions {s_1, s_2, …, s_k} and end positions {e_1, e_2, …, e_k}, the consolidated region [S_m, E_m] was calculated iteratively:
S_m = s_(k+1)/2 if k is odd, (s_(k/2) + s_(k/2+1))/2 if k is even, E_m = e_(k+1)/2 if k is odd, (e_(k/2) + e_(k/2+1))/2 if k is even (10)
Regions outside the interval [S_m, E_m] were iteratively removed, and the process was repeated until no points remained outside. The resulting unified enhancer or promoter regions were { [S_m,i, E_m,i] }_{i=1}^k, where k is the number of unique enhancer or promoter regions.
This approach ensured a consistent reference genome annotation for enhancers and promoters, integrating signals from multiple tissues into a unified framework.

Cis-regulatory elements of enhancers and promoters on genes
To investigate cis-regulatory interactions between enhancers and promoters, all enhancers within ±1 Mbp of transcription start sites (TSS) were identified based on their genomic coordinates. Potential enhancer-promoter combinations were predicted using three complementary approaches. First, enhancer-gene combinations were identified by calculating the Spearman’s correlation (r ≥ 0.8) between the tissue-specific activity profile of each enhancer (presence or absence across 24 tissues) and the corresponding gene expression (TPM) across the same tissue. Second, CTCF ChIP-seq data from cerebellum and liver tissues, collected in a previous study7, were used to identify CTCF-mediated chromatin loops. CTCF peaks were overlapped with CTCF motifs identified using FIMO to determine CTCF orientation59,60. Using this data, loops were identified following the pipeline from ref. 61. Loop containing single promoter and enhancers from the same cluster were used as candidate enhancer-promoter cis-regulatory combinations and validated across all 24 tissues. Third, enhancer-promoter combinations were identified based on promoter and enhancer links from CAGE data using the CAGEfightR v1.12.055. These potential combinations from the three sources were combined, and those with an FDR < 0.05 (Benjamini-Hochberg correction) were retained as confident enhancer-promoter interactions. For each tissue, TPM values were compared between expressed genes with identified enhancers and those without, to evaluate the regulatory impact of enhancers. Additionally, enhancer-promoter interactions were categorized into two primary patterns: enhancer hubs (based on the number of promoters associated with a single enhancer) and promoter hubs (based on the number of enhancers associated with a single promoter). This comprehensive integration of multiple datasets provides a robust framework for elucidating tissue-specific cis-regulatory mechanisms.

Tissue-specific enhancer and promoter
To investigate tissue-specific regulatory elements, we identified and compared genes linked to tissue-specific enhancers and promoters derived from enhancer-promoter combinations. Gene Ontology (GO) enrichment analysis was conducted on genes associated with tissue-specific enhancers and compared to those linked to tissue-specific promoters using Enrichr-KG (https://maayanlab.cloud/enrichr-kg)62. Additionally, tissue-specific enhancer sequences were analyzed to identify associated transcription factors (TFs) using the HOMER annotation tool (http://homer.ucsd.edu/homer) with the annotatePeaks.pl function63. This approach provides insight into the functional roles and regulatory networks associated with tissue-specific enhancers and promoters.

Comparisons of enhancer and promoter across five species
To analyze species-specific regulatory elements in sheep, we compared ChIP-seq data with four other mammalian species, including cattle, pig, mouse, and human. Histone modification data from four tissues (cerebellum, cerebral cortex, lung, and muscle) were used for this comparison. Data for cattle64, pig9, mouse (ENCSR738IWZ), and human (ENCSR000FYQ) were sourced from publicly available datasets. Promoter regions in these species were defined as regions overlapping a flanked TSS region (upstream 2.5 kb and downstream 1 kb)65. And the annotated enhancer regions of cattle, pig and mouse were downloaded from (GIGA)nDB66, while human enhancer regions were obtained from the EnhancerAtlas 2.0 database67.

Whole-genome alignment chain files for human (hg38), mouse (mm10), pig (SusScr11), and cattle (bosTau9) were downloaded from the University of California, Santa Cruz (UCSC) database. The chain files processed using reciprocal best chain methods from the UCSC Genome Wiki (http://genomewiki.ucsc.edu/index.php). Genomic coordinates were converted between the sheep reference and other mammal assemblies using the UCSC Liftover tool (https://genome.sph.umich.edu/wiki/LiftOver) with a threshold of 0.8 sequence identity. Human enhancer and promoter regions were used as a reference, and annotations for sheep, cattle, pig, and mouse were lifted over to human coordinates for direct comparison. Shared enhancers and promoters from both sheep and cattle were identified as the ruminant-specific regulatory elements. The KEGG pathway analysis of genes associated with ruminant-specific enhancers and promoters were performed using KOBAS 3.068. Pathways with p value <0.05 were considered significant.

Biological insights into enhancers and promoters
Whole-genome sequencing (WGS)-based SNPs of 355 sheep were obtained from the iSheep database (https://ngdc.cncb.ac.cn/isheep/download), which contains 15,816,440 SNPs mapped to ARS_UI_Ramb_v2.0. The genomic distribution of these SNPs was analyzed across exonic, intronic, promoter, and enhancer regions to identify their functional implications. To explore the biological functions of SNPs within enhancer and promoter regions, we retrieved sheep QTL data from the sheep QTLdb (Release 57, Aug 26, 2025; https://www.animalgenome.org/QTLdb), which includes 8,288 QTL records curated from 293 publications covering 182 base traits and 266 trait variants, providing valuable insights into the functional significance of regulatory elements in sheep. Additionally, birth weight-associated SNPs (28,736 variants with P value ≤ 0.001) identified from a genome-wide association study (GWAS) of 3007 sheep were analyzed, together with 181 eQTL loci, to further validate the reliability and biological relevance of the identified regulatory elements29.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors present a comprehensive collection of enhancers and promoters for sheep. This will undoubtedly be an important resource for studies aimed at bridging the genome to phenome gap in this species. The study is thorough and well described and the authors already provide a few interesting examples of the power of this new data set to do so. Nevertheless, I was not completely convinced about the story for the QTL for milk fat that is presented in lines 437-443. While the authors show that there is an association between this QTL and an enhancer interacting with the COMMD1 gene, they do not present further evidence that indeed differences in expression of this gene are causing this QTL effect. I think the authors should tone down this claim and indicate that this is an interesting strong candidate or provide further evidence that indeed the expression of this gene is responsible for the observed QTL.\n\nIn fact, I thought the section about biological insights of cis-regulatory elements (lines 427-449) the weakest section of the manuscript. The authors mention in lines 428 of this section (and in line 842 of the methods), that they identified 577,834 SNPs. It is not clear how these SNPs were identified. They refer to the HD_SNP data, but no further description or reference is provided for this data. I was also very surprised with the results presented in figure 7A. Normally one would expect a higher SNP density in intergenic regions and introns compared to functional sequences like exons, enhancers and promoters. It makes me wonder whether the HD_SNP data represents the SNPs that are present on the HD genotyping array for sheep. However, in that case these results are not very meaningful as this collection of SNPs will be highly biased and not provide any insight in the actual variation observed in sheep genomes.\n\nI think the results described in lines 229-241 represent “circular reasoning”. The authors have identified “enhancer-promoter combinations based on correlations of gene expression and enhancer activity” (line 230). The observation that “in tissues where a cis-enhancer was present, the corresponding genes were predominantly expressed, whereas non-expressed genes consistently lacked an associated enhancer in these tissues” is a direct result of the definition to identify such enhancers in the first place.\n\nIt was not clear what exactly is presented in figure 3F. Do the numbers on the X-axis refer to number of enhancers and promoters in the promoter and enhancer hubs respectively? What than is represented by the heatmap (range of 0-1000)?\n\nThe GO analyses described in lines 292-297 is done for the 281 genes from promoters and the 432 genes from the enhancers. It would be interested to also see the results of a GO analysis for the 151 genes that were specific for the enhancers.\n\nIn figure 4C the structure of the BDFN gene is shown as well as 3 different transcripts for this gene in the brain. The gene shows an exon 4 that is not present in any of these three transcripts. Are there other transcripts (in other tissues) that do contain this exon? If so, this could briefly be mentioned for clarity.\n\nThe Supplementary tables are not labeled which sometimes made it difficult to know for sure if I was looking at the correct table. E.g. Table S16 I assumed is data set 17 (file 572519_0_data_set_10118522_spclhz.xlsx). I assume this table presents the sheep enhancers and promoters for the four selected tissues used for cross validation with other species and that the first 3 columns present the location in the human genome? This also is not described anywhere in this table. There are several other tables for which it is not very clearly described what exactly is presented. Please add some further information.\n\nLines 331-337 and figure 4E. The authors describe tissue specific enhancers for cerebellum and for the cortex. However, in figure 4E only cerebellum specific enhancers are shown. The authors also should indicate what the numbers 5, 6, 9 represent (I assumed the number of enhancers that are not all shown).\n\nThe conclusion section is rather repetitive and repeating information already presented twice (in the results and in the discussion sections). Please shorten the conclusion section and make more concise.\n\nMinor comments:\nLine 46: Insert “into” between “insights” and “tissue-specific”\nLines 241-242: I suggest changing this sentence to: “For example, the gene ADORA1 on chromosome 12, which encodes the adenosine A1 receptor with important functions in the brain.”\nIn line 250 the authors introduce the terms “promoter and enhancer hubs”. However, how these are defined is only described in the Material and Methods section in lines 802-804. I think it would help the reader if that were done in this section.\nThe sentence in lines 297-298 sounds rather awkward: “However, in rest 10 tissues, using promoter-enhancer combinations for GO analysis yielded greater associations,”\nThe sentence “red lines in RNA-seq track are the joint reads number spread different exons.” also sounds rather awkward. Furthermore “red” should be with a capital “Red”.\nLine 362-363: Change “within enhancers and promoters from combinations” to “within enhancer-promoter combinations”\nLines 400-403. Although the authors refer to earlier findings for this observation in the discussion (Villar 2015, reference 31), I think they already should do that here.\nFigure 6C. Consider adding a supplementary table showing the details for all these individual promoters.\nLines 509-510: Suggest changing “By providing unified genomic coordinates of reference genome” to “By providing unified reference genomic coordinates”", "reviewer_2": "This article presents comprehensive multi-tissue and multi-omics functional annotation data for sheep, providing a significant data foundation for understanding sheep evolution and advancing sheep breeding. However, given the quality of the data currently presented in the article, I have reservations regarding the reliability of the results proposed by the authors. If these major concerns cannot be resolved, I do not support the publication of this article in Nature Communications.\n\nMajor concerns:\n1. The data presented in the article are indeed remarkable in terms of the number of tissues and the diversity of data types. However, it appears that none of the data includes replicates, not even technical replicates. This is particularly concerning for epigenomic data, which heavily relies on the quality of library preparation experiments. Therefore, the inclusion of replicates is essential.\n2. Although the Fraction of Reads in Peaks (FRiP) for ChIP-seq data is not directly comparable across different proteins, the FRiP values within the same type of ChIP-seq data in this study exhibit significant variability. For instance, the FRiP for H3K27ac data ranges from 2% to 22%, with a median of only 6.93%, indicating overall low data quality and substantial variability. Additionally, fluctuations in data quality are also observed in the reads mapping ratio metric. As an example, for the Tonsil tissue, the reads mapping ratio for the control is merely 27.67%, whereas the mapping ratios for other ChIP-seq data from Tonsil are above 80%. This discrepancy likely affects the number of peaks and FRiP values for all ChIP-seq data from Tonsil. Overall, these variations are more likely attributable to the quality of library preparation rather than sample specificity. Such significant variability undermines the biological relevance of the statistical results presented in Figure 2 and related sections. The quality of the data is concerning.\n3. Recently, a similar study on the sheep epigenomic landscape has been published (https://doi.org/10.1002/imt2.254). The primary innovation of the current study, in comparison, lies in the generation of data for brain related tissues. The authors have devoted considerable effort to describing cerebellum-specific enhancers and their associated genes. However, the related biological insights remain insufficient. A compelling narrative could emerge from comparing sheep-specific brain-related regulatory elements with previous data on sheep domestication. Such an analysis might provide a more intriguing and insightful story.\n\nMinor suggestion.\n1. The authors identified candidate regulatory enhancers within a 1 Mbp range upstream and downstream of transcription start sites (TSS). This distance is excessively large and likely exceeds the boundaries of topologically associating domains (TADs). It is recommended to utilize TADs for the identification of candidate regulatory enhancers, as this approach would provide a more accurate and biologically relevant framework for understanding gene regulation.\n2. In line 756 ‘represent’ should be ‘represented’\n3. In line 790, please clarify the meaning of ‘enhancer activity’.", "reviewer_3": "This paper by Xie et al presents a very useful data set on the functional annotation of promoters and enhancers in sheep by analysing a range of techniques on 24 tissues from a single ewe, including ChIPSeq, ATAC Seq, CAGE Seq, methylation studies by WGBS and RRBS, and transcript analysis. These data are highly complementary for characterising promoters and enhancers and their functionality in different tissues. We would also have liked to have HiC data to reinforce the study of the enhancer-promoter link, but it is likely that the samples required for this were not available. This would have strengthened the relevance of the promoter-enhancer link, since you are taking into account distances of up to 2 Mb.\nEven if we can always expect more, globally I find this paper very complete and well written. For me it deserves publication in Nature Communications.\n\nGeneral Comments :\nI would have appreciated it if the results presented had focused a little more on the functions of interest in ruminants. If you could add something to tell us a story about this, either specific to sheep or common to sheep and cattle, I would find that interesting.\nYou have a dataset for which you have the complete genomic sequence of the individual being analysed. You could use the RNAseq data to see if you have allele specific expression ... and if that allele specific expression can be linked with a possible polymorphism in an enhancer or promoter region of that same gene, with information by other techniques, showing that this promoter or enhancer is non functional in the studied tissue. It would be great.\nI was lost by the numbering of several additional files tables. Can you check numbering when you cite a supplementary data file please. It would facilitate to name the excel file tab with the number of the corresponding supplementary table number\nTo make it easier to understand the contribution of each of the techniques, could you propose a figure (possibly in supplementary material) showing the respective contributions of ChIPSeq with the 4 histone marks, ATAC Seq, CAGE-Seq, RRBS and WGBS, as well as RNA Seq? This would better highlight the respective contribution of each of these different techniques.\nSee how you can make your data available more easily than only in supplementary files\n\nDetailed comments :\nIt would be a good idea to state in the introduction, line 133, that the genome of the animal you have used has been sequenced and assembled to produce a de novo reference genome.\nFigure 1B :\nThe color scale does not clearly show the number of tissues analyzed for each of the techniques used.\nDifferentiating WGBS and RRBS on 2 different tracks\nI wonder if a 2-colour scale wouldn't be more visual and informative. The same color scale could also be used for all the applications used.\nFigure 1C & figure S1 :\nIt would be desirable for the 4 categories TSS, exons, introns, intergenic to be in the same order ... and with similar colors\nIn Figure S1, the color of the green is not the same for the graph showing the percentage and that showing the average over the 24 tissues. It is not so important, but it would be nicer.\nFigure 1E :\nThe scale, which is neither log nor standard, does not appear to be the best. The SD value for the poised enhancer, for example, is not presented on the same scale for the minimum and maximum values.\nFor a better understanding, perhaps write also the % of the genome occupied by each of the states.\nLine 179 : it seems that the figure of 27.1% corresponds to the sum of the repressed polycomb and weak repressed polycomb states, and not simply to repressed polycomb.\nFig 2C :\nLegend of y axis : number of enhancer   number of enhancers\nLines 206-208 : These numbers might be added on fig 2A\nFig S3B, S4C & S4D : On y axis, it is numbers rather than frequency\nA question, seeing fig 2C & D: Would a heatmap showing the number of enhancers and promoters for a pair of tissues give anything interesting to see in terms of the use of these promoters?\nIn Figure 2B, rather than simply showing the number of enhancers and promoters identified in each tissue, would it be interesting for each tissue to distinguish the number of promoters specific to that tissue and promoters observed in more than one tissue? The number of promoters and enhancers could thus be represented, highlighting their specific or more or less ubiquitous characteristics.\nLines 228-229 : Of the 9 genes indicated in this table, 4 show a comparison with a sheep gene. Is it 9 genes or 5 genes? In any case, Ovis Aries is not mentioned in the species indicated on lines 228-229.\nLine 239 : in any tissue   in any tested tissue.\nLine 243: Since you are using TPM for the first time, you should define it.\nFigure 3A & 3B:\nAre you sure that TPM reaches 10^5, which means that for some tissues 10-20 % transcripts are from the same gene ?\nAs you considered that a gene is considered as expressed if TPM is > 0.5, it means that quite all genes without enhancer are in fact NOT EXPRESSED. Am I wrong ?\nOn Figure 3B, we see a significant number of colored dot for genes with a TPM <0.5, and thus considered as non expressed, but you write line 247 that “these non expressed genes lacked an associated enhancer”. So you should not have colored dots below a TPM of 0.5. Am I wrong ?\nLine 260 : 51.52% of genes, or 51.52 % of promoters ?\nLines 260 & 261 : 52% not 51.52%. 89% not 88.53% (give fewer numbers)\nFigure 3F:\nI wonder if a 2 colors scale can be more visual and informative\nEven if you used a log scale, numbers are not in log2 as indicated in the title, as you indicate 2^5 and 2^10. So remove “(log2)”.\nFigure 4A:\nIt is not clear for me. I understand that 281 + 151 = 432, but I do not understand.\nWhat means “specific genes from promoters”\nTo what correspond these 151 genes in light orange ?\nFig 4C :\nYou indicate in the text a single isoform in the cortex, but you indicate 2 and 4 in the RNASeq track, suggesting 2 different isoforms exists.\nFig 5 :\nI wonder if using a log scale for the y axis is relevant (moreover you show 0 on a log scale). Have you tested a standard scale? It might better put in evidence the slope of the correlation.\nLine 382 : There is no fig 5C & 5D, you have 2 figures on 5A and 2 on 5B. Alternatively change the figure and its legend.\nLine 404 : I did not find table S16. I should say that the supplementary tables are not properly numbered, at least in the version I downloaded for review. For a simpler reading, the tab title of Excel files should indicate the table number. I did not also find the legend of supplementary tables.\nLine 407 : the value of 24467 enhancers in the cortex do not correspond the value on fig S10A, which is 42467. So in the text 24467   42467. With this value, the indicated percent is correct.\nLine 406 & 410 : use fewer numbers to show the percentages\nLines 412-413 : If you intend to say that sheep is closer to cattle than to pig, you should say that on fig 6B, sheep is even closer to human than to cattle, which is funny ! Is it due to a really larger number of promoters put in evidence in human (as suggested on fig S11), which may induce that a larger number of promoters identified in sheep can also be found in human ?\nFor figure S11, I wonder if it wouldn't be more informative simply to have (or in addition to) a table with the 5 species in rows and columns, and at the intersection of the rows and columns the percentage of promoters (or enhancers) of the species in that column found in the species in the row. This would better highlight the proximity between species :\nNote : the array will not be symmetrical because the number of promoters found in the two species under consideration is not the same.\nLine 414 : you may add the percentage of enhancers and promoters found in all 5 species. Based on your data it is on average for the 4 tissues around 0.8 % for enhancers and 17% for promoters.\nFig 6D : it would be easier to have the tissue in the same order in the text and in the figure, whereas we have cerebellum, cortex, lung and muscle in the text, but cerebellum, lung, muscle and cortex on the figure.\nLines 431-433: I think that in the discussion part, you indicate 69 pairs of enhancer-promoter conserved in cattle and sheep and not conserved in other species. If I am right, this number can be indicated also in this result part..\nLine 445-446 : You did not explain here or in M&M how these SNPs have been identified. This greatly reduces the scope of your conclusions, and in particular the frequency of presence of the polymorphism in the various compartments (exon, intergenic region, enhancers, promoters, etc.) which may be influenced by the way in which the polymorphism was identified. In addition, was the polymorphism identified following sequencing of the ewe, or of a large number of individuals from different breeds?\nFig 7a : If you used density expressed in % on the legend of y axis, replace 0.00022 by 0.022 on the figure\nLine 473 : such as; enhancers   such as enhancers\nLines 578-580 : could you check whether this number of 69 conserved enhancer-promoter pairs is higher than the number you would expect at random between two species taken at random from the 6 species considered.\nLigne 594 : For the QTL, I find the word « associated » too strong, as it suggests causality. Presently it is only a colocalisation. You also say : “this study successfully identified 38 QTL associated with …”. I would say in the other direction that you demonstrate the presence of promoters or enhancers in regions where QTL have been identified.\nM&M :\nIn tissue collection, you may remember that the genome of this ewe has been de novo sequenced. Indicate also the ref of the publication (29).\nLine 734 : As you reported the number of RNA seq reads per tissue of 65-125 Millions, you might add in the text, that it means that for most tissues, fragments of genes considered as expressed are found at least 30 fold.\nLines 766 & 767, if you use _Ti underscript to indicate the tissue, you should do the same for the 4 types PH3K4me1 , PH3K27ac, PCAGE-bid, PATAC.\nLine 773 : represent   represented ?\nLine 795 :\nI am not sure I understand how you calculate Sm & Em. You use Sn, but n is the last tissue. May be you should use an other letter than n. In fact Sn+1 does not exist.\nBe careful you have large s and e for n odd, and small s and e for n even.\nIf the formula (s_(n/2)+s_(n/2+1))/2 is correct, it could be more readable writing (s_(n/2)+s_((n+2)/2))/2 or\n(s_(n/2)+s_((n/2)+1))/2\nLine 804 you indicate ± 1 Mbp of TSS, but in line 545 you indicates “2 Mbp range of promoter”\nLine 859 : can you indicate the source of this collection of SNP ?", "reviewer_1_attached": "The authors present a comprehensive collection of enhancers and promoters for sheep. This will\nundoubtedly be an important resource for studies aimed at bridging the genome to phenome gap in this\nspecies. The study is thorough and well described and the authors already provide a few interesting\nexamples of the power of this new data set to do so. Nevertheless, I was not completely convinced\nabout the story for the QTL for milk fat that is presented in lines 437-443. While the authors show that\nthere is an association between this QTL and an enhancer interacting with the COMMD1 gene, they do\nnot present further evidence that indeed differences in expression of this gene are causing this QTL\neffect. I think the authors should tone down this claim and indicate that this is an interesting strong\ncandidate or provide further evidence that indeed the expression of this gene is responsible for the\nobserved QTL."}, "author_response": "Response: Thank you very much for this valuable suggestion. We agree that our current evidence\nhighlights a strong colocalization between the milk-fat QTL, the candidate enhancer, and COMMD1\nexpression, but does not definitively prove that differential expression of COMMD1 is responsible for\nthe observed phenotype. In the revised manuscript, we have therefore toned down our wording:\ninstead of stating that the enhancer-gene relationship is “confirmed,” we now describe it as\n“suggestive,” and we refer to COMMD1 as a compelling candidate gene rather than a proven causal\ngene. We have also added supporting literature linking COMMD1 to copper and lipid metabolism in\nmammals27,28 (Li et al., 2020; Wang et al., 2012), which further strengthens the biological plausibility\nof its role in milk-fat regulation, while acknowledging that functional validation will be required to\nestablish causality.\n\nAdditionally, to obtain more direct evidence of enhancer-associated gene expression regulation, we\nanalyzed 23,056,026 SNPs from a recent genome-wide association study (GWAS) of birth weight (BW)\nin 3,007 sheep29, and 181 181 published significant eQTLs from the same study. Intersecting these\nvariants with our enhancer-promoter interactions across 24 tissues, we identified the most significant\nSNP (P = 2.23 × 10⁻⁷) was located within a cerebellum-specific enhancer predicted to regulate XKR4. In\nthe same GWAS study29, XKR4 had previously been associated with BW based on an intragenic SNP.\nOur analysis not only confirms this association but also pinpoints an additional variant within a\nregulatory enhancer controlling XKR4, thereby providing new evidence that variation within regulatory\nelements can influence phenotypic traits in sheep. In eQTL result, we identified one variant\nsignificantly associated with LOC121819880 expression (P = 4.69 × 10⁻⁸), which resides within its\ninteracting enhancer (chr6: 35907531-35908224). This observation provides additional evidence that\nsequence variants within enhancers can directly influence gene expression, reinforcing the accuracy\nand biological validity of our enhancer annotations\nWe have revised the relevant text in the Results, Discussion, and M&M sections of revised manuscript.\n\nIn fact, I thought the section about biological insights of cis-regulatory elements (lines 427-449) the\nweakest section of the manuscript. The authors mention in lines 428 of this section (and in line 842 of\nthe methods), that they identified 577,834 SNPs. It is not clear how these SNPs were identified. They\nrefer to the HD_SNP data, but no further description or reference is provided for this data. I was also\nvery surprised with the results presented in figure 7A. Normally one would expect a higher SNP density\nin intergenic regions and introns compared to functional sequences like exons, enhancers and\npromoters. It makes me wonder whether the HD_SNP data represents the SNPs that are present on the\nHD genotyping array for sheep. However, in that case these results are not very meaningful as this\ncollection of SNPs will be highly biased and not provide any insight in the actual variation observed in\nsheep genomes.\nResponse: We thank the reviewer for this constructive comment. In the revised analysis, we replaced\nthe previously used HD_SNP dataset with whole-genome sequencing (WGS)-based SNPs obtained\nfrom the iSheep database (https://ngdc.cncb.ac.cn/isheep/download), which contains 15,816,440\nSNPs identified from WGS sequenced sheep genomes rather than array-based SNPs. Using this WGS-\nderived dataset, we repeated the SNP genomic feature distribution and density analysis (Figure 7A\nand B) and now observe the expected enrichment of SNPs in intergenic and intronic regions compared\nto exons, enhancers, and promoters. We have updated both the M&M and Results sections to clearly\ndescribe the source of SNPs.\n\nI think the results described in lines 229-241 represent “circular reasoning”. The authors have identified\n“enhancer-promoter combinations based on correlations of gene expression and enhancer activity” (line\n230). The observation that “in tissues where a cis-enhancer was present, the corresponding genes were\npredominantly expressed, whereas non-expressed genes consistently lacked an associated enhancer in\nthese tissues” is a direct result of the definition to identify such enhancers in the first place.\nResponse: We sincerely thank the reviewer for this thoughtful comment. We would like to clarify that\nthe identification of enhancers and promoters in our study was not solely based on correlations with\ngene expression. We employed an integrative strategy to identify enhancer–promoter links. Instead,\nwe applied an integrative, multi-layered strategy to define enhancer-promoter interactions.\nSpecifically, the ±1 Mbp window around transcription start sites (TSSs) served as an initial search\nrange to capture potential distal regulatory elements. The final enhancer-promoter interactions were\nestablished by integrating three independent lines of evidence: (i) correlation between enhancer\nactivity and gene expression across 24 tissues; (ii) Chromatin loop information inferred from CTCF-\nmediated interactions; and (iii) CAGE-based enhancer-promoter associations. Only interactions\nsupported by these orthogonal criteria and passing an FDR < 0.05 were retained, ensuring high\nconfidence in their biological relevance. Thus, our enhancer–promoter relationships were not defined\nby expression correlation alone, but rather by converging evidence from multiple omics layers.\nFurthermore, for the quantitative correlation analysis, we considered genes with TPM ≥ 0.5 as\nexpressed. However, in subsequent qualitative comparisons of gene expression with and without\nidentified enhancers across different tissues (e.g., Figures 3A and 3B), we deliberately did not apply\nexpression-based filtering, to avoid potential bias introduced by predefined thresholds.\n\nIt was not clear what exactly is presented in figure 3F. Do the numbers on the X-axis refer to number of\nenhancers and promoters in the promoter and enhancer hubs respectively? What than is represented by\nthe heatmap (range of 0-1000)?\nResponse: Thank you very much for this valuable comment. We have clarified the presentation of\nFigure 3F. The X-axis is divided into two parts: promoter hubs (E1–E20, E>20) and enhancer hubs (P1–\nP5, P>5). For the promoter hub (E1–E20, E>20), the labels indicate the number of enhancers\nconnected to each promoter hub. For example, “E1” refers to promoter hubs that are linked to exactly\none enhancer, “E2” refers to promoter hubs with two enhancers, and so forth, while “E>20” indicates\nhubs that are connected with more than 20 enhancers. For the enhancer hub (P1–P5, P>5), the labels\nindicate the number of promoters connected to each enhancer hub. For instance, “P1” means\nenhancer hubs that are linked to a single promoter, “P2” means enhancer hubs linked to two\npromoters, and “P>5” represents enhancer hubs that connect to more than five promoters.\nThe heatmap color scale (log2 transformed in the figure) represents the counts of hubs in each tissue\nthat fall into the corresponding category. For example, the value at “E1” for “CerebralCortex” indicates\nthe number of promoter hubs in cerebral cortex that are linked to exactly one enhancer. Similarly, the\nvalue at “P3” for a given tissue indicates the number of enhancer hubs in that tissue connected to\nthree promoters.\nTo improve clarity, we have revised the figure legend in the revised manuscript:\n“F Distribution of enhancer and promoter numbers in pair hubs across tissues (The X-axis shows the\nnumber of enhancers (1, 2, 3, …, 20, and >20) per promoter hub and the number of promoters (1, 2, 3,\n4, 5, and >5) per enhancer hub. The heatmap represents the counts of promoter hubs and enhancer\nhubs in each tissue with the indicated number of enhancer or promoter connections).”\n\nThe GO analyses described in lines 292-297 is done for the 281 genes from promoters and the 432 genes\nfrom the enhancers. It would be interested to also see the results of a GO analysis for the 151 genes that\nwere specific for the enhancers.\nResponse: We thank the reviewer for this constructive suggestion. Following the recommendation, we\nperformed GO enrichment analysis specifically for the 151 genes uniquely associated with enhancers\nin the cerebellum. This analysis revealed five significantly enriched biological processes: positive\nregulation of synaptic transmission, synapse assembly, positive regulation of neuron differentiation,\nmodulation of chemical synaptic transmission, and regulation of neuron differentiation\n(Supplementary Table S14). Importantly, we identified 11 key genes that were supported by two or\nmore of these processes. Among them, NLGN3, BDNF, and NRCAM appeared in three pathways\n(Supplementary Figure S8), while RIMS2, RIMS1, SYT1, SLC8A2, NEUROD1, NEUROD2, BRINP1, and\nVWC2 were shared by two pathways. These genes represent core neuronal regulators with pleiotropic\nroles in synaptic transmission, neuronal differentiation, and circuit assembly. Their recurrence across\nmultiple enriched categories underscores that enhancer-specific genes in the cerebellum converge on\nfundamental mechanisms of neuronal development and synaptic function, thereby reinforcing both\nthe robustness and biological relevance of our findings. We have included these in the revised\nmanuscript.\nFigure. S8. Network of GO enrichment analysis of enhancer-specific associated 151 genes in the\ncerebellum. (pink nodes: five significantly GO biological processes; and green nodes: the 11 genes)\n\nIn figure 4C the structure of the BDFN gene is shown as well as 3 different transcripts for this gene in the\nbrain. The gene shows an exon 4 that is not present in any of these three transcripts. Are there other\ntranscripts (in other tissues) that do contain this exon? If so, this could briefly be mentioned for clarity.\nResponse: We thank the reviewer for this helpful comment. In Figure 4C, we illustrated three\nrepresentative brain-specific transcripts of the BDNF gene, which do not include exon 4. We would like\nto clarify that exon 4 is indeed annotated in other BDNF transcripts expressed in non-brain tissues. For\nexample, in the NCBI database (release 104), transcripts XM_042233304.1 contain exon 4 and are\nreported to be expressed with RNA-seq reads\n(https://www.ncbi.nlm.nih.gov/projects/sviewer/?id=NC_056068.1&tkey=5syE4enu4uvs5ejnw9L369G\nZi4Ssl5evu4SbsZIFvZPol_lH0CeiT6QUlzyNCP44jVw&assm_context=GCF_016772045.1&app_context=g\nenome&v=57132748:57206019&c=null&select=null&slim=0).\nWe have added a brief note in the Results section to clarify this point and avoid potential confusion.\n\nThe Supplementary tables are not labeled which sometimes made it difficult to know for sure if I was\nlooking at the correct table. E.g. Table S16 I assumed is data set 17 (file\n572519_0_data_set_10118522_spclhz.xlsx). I assume this table presents the sheep enhancers and\npromoters for the four selected tissues used for cross validation with other species and that the first 3\ncolumns present the location in the human genome? This also is not described anywhere in this table.\nThere are several other tables for which it is not very clearly described what exactly is presented. Please\nadd some further information.\nResponse: We thank the reviewer for pointing this out. We have revised all supplementary materials\n(Supplementary Table S1-S27 and Supplementary Figure S1-S13) by adding explicit legends and\ndescriptions for each table and figure. In particular, we now state what dataset each table represents,\nthe biological content of the data, and the meaning of each column (including whether genomic\ncoordinates correspond to the sheep or human genome).\n\nLines 331-337 and figure 4E. The authors describe tissue specific enhancers for cerebellum and for the\ncortex. However, in figure 4E only cerebellum specific enhancers are shown. The authors also should\nindicate what the numbers 5, 6, 9 represent (I assumed the number of enhancers that are not all shown).\nResponse: We thank the reviewer for this valuable comment. In Figure 4E, we presented only a subset\nof tissue-specific enhancers due to space limitations. The numbers 5, 6, and 9 indicate the number of\nenhancers identified for a given gene in the corresponding tissue (as illustrated across adrenal gland\nand brain). For clarity, we have revised the figure legend to explain this notation.\nIn addition, we have now included two cortex-specific genes and their enhancers in Figure 4E. The\ncomplete list of tissue-specific enhancers across all examined tissues is provided in Supplementary\nTable S17.\n\nThe conclusion section is rather repetitive and repeating information already presented twice (in the\nresults and in the discussion sections). Please shorten the conclusion section and make more concise.\nResponse: Thank you for this helpful suggestion. In the revised manuscript, we have substantially\nshortened and streamlined the conclusion to avoid repetition and to present the key take-home\nmessages more concisely.\n\nMinor comments:\nLine 46: Insert “into” between “insights” and “tissue-specific”\nResponse: We have inserted this into the abstract section of the revised manuscript.\nLines 241-242: I suggest changing this sentence to: “For example, the gene ADORA1 on chromosome 12,\nwhich encodes the adenosine A1 receptor with important functions in the brain.”\nResponse: We have revised this as suggested.\nIn line 250 the authors introduce the terms “promoter and enhancer hubs”. However, how these are\ndefined is only described in the Material and Methods section in lines 802-804. I think it would help the\nreader if that were done in this section.\nResponse: We have added a brief definition in the main text of revised manuscript.\nThe sentence in lines 297-298 sounds rather awkward: “However, in rest 10 tissues, using promoter-\nenhancer combinations for GO analysis yielded greater associations,”\nResponse: We have revised this in the revised manuscript.\nThe sentence “red lines in RNA-seq track are the joint reads number spread different exons.” also sounds\nrather awkward. Furthermore “red” should be with a capital “Red”.\nResponse: We have revised this in the revised manuscript.\nLine 362-363: Change “within enhancers and promoters from combinations” to “within enhancer-\npromoter combinations”\nResponse: We have revised this as suggested.\nLines 400-403. Although the authors refer to earlier findings for this observation in the discussion (Villar\n2015, reference 31), I think they already should do that here.\nResponse: We have included this reference in main text as suggested.\nFigure 6C. Consider adding a supplementary table showing the details for all these individual promoters.\nResponse: We have added the supplementary Table S19.\nLines 509-510: Suggest changing “By providing unified genomic coordinates of reference genome” to “By\nproviding unified reference genomic coordinates”\nResponse: We have revised this in the revised manuscript.\n\nReviewer #2 (Remarks to the Author):\nThis article presents comprehensive multi-tissue and multi-omics functional annotation data for sheep,\nproviding a significant data foundation for understanding sheep evolution and advancing sheep\nbreeding. However, given the quality of the data currently presented in the article, I have reservations\nregarding the reliability of the results proposed by the authors. If these major concerns cannot be\nresolved, I do not support the publication of this article in Nature Communications.\nResponse: We sincerely thank the reviewer for carefully evaluating our work and for raising this\nimportant concern about data reliability. We have taken this comment very seriously and have made\nadditional efforts to strengthen the robustness of our analyses. Specifically:\n(1) Reproducibility across replicates: We have performed technical replicates from 16 tissues (all\ntissue with enough remaining samples), re-examined and confirmed that the consistency of\noriginal and replicates data is high, showing strong concordance with the results presented in\nthe original submission. This provides further confidence in the reliability of our data.\n(2) Integration of multiple independent omics datasets: Our analysis strategy was designed to\nminimize potential biases by integrating multiple independent assays, including ChIP-seq,\nCAGE-seq, ATAC-seq, and RNA-seq. In particular, candidate enhancers and promoters were\ndefined only when supported by at least two independent tissues and validated by at least\ntwo different omics datasets (as described in the Methods). This stringent approach ensures\nthat our functional element annotations are robust and reliable.\n(3) Biological plausibility: Importantly, the patterns we observed for tissue-specific enhancers and\nthe results of pathway analyses are consistent with known biological functions of the\ncorresponding tissues. For instance, XKR4, linked with cerebellum-specific enhancers, has\nbeen implicated in neuronal development and cerebellar function. This concordance provides\nan additional layer of validation and indicates that our findings are not only statistically\nsupported but also biologically meaningful. These results are consistent with known biological\nfunctions and thus support the robustness of our integrative strategy.\nWe greatly appreciate the reviewer’s concerns, which have prompted us to carry out additional\nvalidations and to clarify our methodology more explicitly in the revised manuscript. We hope that\nthese efforts will help address the reviewer’s concerns and demonstrate that the results we present\nare both reliable and valuable to the community.\n\nMajor concerns:\n1. The data presented in the article are indeed remarkable in terms of the number of tissues and the\ndiversity of data types. However, it appears that none of the data includes replicates, not even technical\nreplicates. This is particularly concerning for epigenomic data, which heavily relies on the quality of\nlibrary preparation experiments. Therefore, the inclusion of replicates is essential.\nResponse: We sincerely thank the reviewer for this important comment regarding the inclusion of\nreplicates. To directly address this concern, we generated replicates from 16 tissues derived from the\nsame individual, with newly obtained ChIP-seq data at higher sequencing depth (Supplementary Table\nS25).\nWe carefully evaluated data quality using both the standard Nextflow ChIP-seq pipeline (https://nf-\nco.re/chipseq) and an independent assessment with MarkDuplicates from Picard v2.26.10. The\nevaluation included sequencing depth, mapping rate, and duplicate rate. Overall, the replicate data\nshowed higher sequencing depth and good quality. In most tissues, duplicate rates were below 20%,\nwhile in a few tissues the rate reached ~50%; nevertheless, the number of uniquely mapped reads\nafter duplicate removal was still above 20M in all cases, which is considered sufficient for reliable\ndownstream analysis.\nIn addition, we examined the consistency of enhancer and promoter annotations by comparing the\nregulatory elements identified in the original data with peaks obtained from replicate ChIP-seq.\nSpecifically, enhancers overlapped strongly with replicate H3K4me1 and H3K27ac peaks, while\npromoters showed concordance with replicate H3K4me3 and H3K27ac peaks. Across tissues, on\naverage, 62% of enhancers and 82% of promoters identified in the original data were supported by\nthe corresponding replicate datasets, and in some tissues, this overlap reached over 90%\n(Supplementary Table S26).\nTaken together, these additional analyses confirm that (i) the replicate data are of high quality, (ii) the\nregulatory elements identified show strong consistency with replicates, and (iii) our integrative multi-\nomics strategy yields biologically reliable enhancer and promoter annotations. We have added these\nresults to a detailed discussion in the revised manuscript.\n\n2. Although the Fraction of Reads in Peaks (FRiP) for ChIP-seq data is not directly comparable across\ndifferent proteins, the FRiP values within the same type of ChIP-seq data in this study exhibit significant\nvariability. For instance, the FRiP for H3K27ac data ranges from 2% to 22%, with a median of only 6.93%,\nindicating overall low data quality and substantial variability. Additionally, fluctuations in data quality are\nalso observed in the reads mapping ratio metric. As an example, for the Tonsil tissue, the reads mapping\nratio for the control is merely 27.67%, whereas the mapping ratios for other ChIP-seq data from Tonsil\nare above 80%. This discrepancy likely affects the number of peaks and FRiP values for all ChIP-seq data\nfrom Tonsil. Overall, these variations are more likely attributable to the quality of library preparation\nrather than sample specificity. Such significant variability undermines the biological relevance of the\nstatistical results presented in Figure 2 and related sections. The quality of the data is concerning.\nResponse: We thank the reviewer for the valuable comments on data quality. In response to the\nconcerns regarding the variability of FRiP values and the mapping ratio of the Tonsil samples, we\nperformed additional analyses and provide further clarification below.\nFirst, regarding the FRiP values, we observed that except for three tissues: AdrenalCortex (22.10%),\nCerebellum (21.70%), and LymphNodeMesenteric (17.90%), all other tissues show FRiP values ranging\nfrom 2% to 13%. This indicates that the overall variability is relatively modest and within the\ncommonly reported range in the literature. For the tissues with higher FRiP values, we found that their\npeak numbers were substantially higher than in other tissues (exceeding 200,000 peaks), which\nreasonably explains the higher proportion of reads in peaks.\nSecond, concerning the low mapping ratio observed in the Tonsil control and H3K27me3 samples, we\nfurther evaluated data quality by considering both total sequencing depth and the number of\ndeduplicated reads used for peak calling. The Tonsil samples had approximately 80M total reads, and\neven with a mapping ratio of 27%, more than 20M deduplicated reads remained available for peak\ncalling, which is sufficient to ensure reliable analysis. Moreover, the corresponding CAGE-seq, ATAC-\nseq, and RNA-seq datasets from Tonsil exhibited good quality, suggesting that the main results are\nrobust.\nThird, and perhaps most important, our strategy for enhancer and promoter identification was\nspecifically designed to minimize the influence of data variability from any single dataset. For\nenhancers, for example: (1) we used the union of H3K4me1 and H3K27ac peaks in each tissue as the\ninitial candidate regions, thereby eliminating the reliance on any single histone mark; (2) these regions\nwere further intersected with CAGE-seq and ATAC-seq peaks from the same tissue to obtain higher-\nconfidence candidates; (3) we combined peaks across all tissues into a unified candidate enhancer set,\nrequiring that each enhancer must be independently identified in at least two tissues to reduce the\nimpact of data quality from a single sample; and (4) we standardized genomic boundaries across the\nreference genome. This integrative framework relies not only on ChIP-seq but also on CAGE-seq and\nATAC-seq datasets, and does not incorporate H3K27me3 data. Therefore, the lower mapping ratio\nobserved in a subset of Tonsil samples does not introduce systematic bias into the final enhancer and\npromoter calls.\nIn addition, we performed both direct and indirect validation of the Tonsil data quality: (1) Direct\nevidence: Based on 20M deduplicated reads, the Tonsil samples yielded 99,509 peaks with a FRiP\nvalue of 6.29%, which is comparable to other tissues, suggesting that higher duplicate levels did not\nsystematically bias peak identification. (2) Indirect evidence: In SpinalCord tissue, which had a similar\nsituation (original control with only 20M deduplicated reads), we compared the original control to a\nhigher-depth replicate control (94M deduplicated reads). The overlap between the peaks from\nreplicate and original enhancers was 68.51%, and promoter concordance reached 93.07%. This\nsupports the conclusion that even with 20M deduplicated reads, reliable results can still be obtained.\nIn summary, while we acknowledge that some variability exists across tissues and samples, our\nintegrative strategy combining multiple omics datasets and cross-tissue validation effectively\nminimizes the impact of individual data quality differences. We are therefore confident that our\nresults remain robust and biologically meaningful.\n\n3. Recently, a similar study on the sheep epigenomic landscape has been published\n(https://doi.org/10.1002/imt2.254). The primary innovation of the current study, in comparison, lies in\nthe generation of data for brain related tissues. The authors have devoted considerable effort to\ndescribing cerebellum-specific enhancers and their associated genes. However, the related biological\ninsights remain insufficient. A compelling narrative could emerge from comparing sheep-specific brain-\nrelated regulatory elements with previous data on sheep domestication. Such an analysis might provide\na more intriguing and insightful story.\nResponse: We sincerely thank the reviewer for this constructive and valuable suggestion, which\nindeed helps enhance the biological insights of our study. The recent sheep epigenomic landscape by\nDeyin Zhang et al. provides an important resource across nine tissues. In this study, we integrate\nmulti-omics data from 24 tissues, including CAGE-seq, and propose a new strategy for enhancer and\npromoter identification and enhancer-promoter interactions. These 24 tissues broadly cover key ovine\ntissues, notably including brain-related tissues cerebellum and cerebral cortex.\nTo further elucidate the contribution of brain-specific enhancers and their target genes to sheep\ndomestication, we analyzed 23,056,026 SNPs from a birth weight GWAS of 3,007 sheep whole-\ngenome sequences data29 (https://doi.org/10.1186/s40104-025-01251-4). Of these, 28,736 SNPs were\nassociated with birth weight at P-value ≤ 0.001. Intersecting these variants with our enhancer-\npromoter interaction links across 24 tissues identified 18 SNPs within 17 enhancers active in\ncerebellum and/or cerebral cortex. Among them, 7 cerebellum-specific enhancers serve as regulatory\nelements interacting with genes, including XKR4, GALNT9, RET, TRNAC-ACA-97, TRNAG-CCC-197,\nTRNAC-GCA-185, LOC114110132, and LOC121817986. Of particular interest, the most significant SNP\n(P = 2.23 × 10⁻⁷) located within a cerebellum-specific enhancer identified to regulate XKR4 gene. In\nthe same GWAS study of sheep29, XKR4 was previously associated with BW based on intragenic SNP.\nOur analysis not only confirms the BW association at XKR4 but also pinpoints an additional variant\nlocated within a regulatory enhancer controlling XKR4, providing a novel association between variant\nwithin regulatory element and phenotype. In addition, GALNT9 encodes a brain-specific enzyme with\nhigh expression in the human cerebellum30 (https://pubmed.ncbi.nlm.nih.gov/10978536), and RET is\nrelated to cell growth, differentiation, and nervous system development31\n(https://www.tandfonline.com/doi/full/10.1080/07853890701646256). We also highlight another\nsignificant SNP (P = 7.47 × 10⁻⁵) within an enhancer active in cerebellum and spinal cord that interacts\nwith TSPAN5 gene; notably, TSPAN5 is highly expressed in the brain during postnatal development in\nmouse32 (https://www.sciencedirect.com/science/article/pii/S2211124719312392).\nThese findings provide direct evidence linking brain-active regulatory elements to a domestication-\nrelated trait (BW) and further support the reliability and accuracy of our enhancer identification\nstrategy and enhancer–promoter interactions. We have included these results in the revised\nmanuscript.\n\nMinor suggestion.\n1. The authors identified candidate regulatory enhancers within a 1 Mbp range upstream and\ndownstream of transcription start sites (TSS). This distance is excessively large and likely exceeds the\nboundaries of topologically associating domains (TADs). It is recommended to utilize TADs for the\nidentification of candidate regulatory enhancers, as this approach would provide a more accurate and\nbiologically relevant framework for understanding gene regulation.\nResponse: We sincerely thank the reviewer for this insightful suggestion. We fully agree that TADs\nprovide a biologically meaningful framework for constraining enhancer–promoter interactions. In our\nstudy, the ±1 Mb window was applied only as an initial search space to capture potential distal\nregulatory elements rather than as a fixed interaction boundary. The final enhancer–promoter links\nwere stringently defined based on multi-omics evidence, including correlation analyses, CTCF-\nmediated chromatin loops, and CAGE-based transcriptional associations, all under FDR control.\nAlthough a Hi-C dataset from sheep liver is available from our previous work20, its tissue coverage and\nsequencing depth are insufficient to delineate high-confidence TAD boundaries across all 24 tissues\nanalyzed in the present study. Therefore, we adopted an integrative multi-omics strategy that\nemphasizes both spatial and functional reliability. We anticipate that as high-resolution Hi-C data\nbecome available across additional tissues, our framework can be further refined and enhanced to\ncapture TAD-level organization more precisely.\n2. In line 756 ‘represent’ should be ‘represented’\nResponse: Thank you, we have corrected “represent” to “represented” in the revised manuscript.\n3. In line 790, please clarify the meaning of ‘enhancer activity’.\nResponse: We thank the reviewer for this helpful comment. In this study, “enhancer activity” refers to\nthe tissue-specific activation status of candidate enhancers across the 24 profiled tissues. For each\nenhancer, we determined its activity in each tissue based on enhancer identification, thereby\ngenerating an activity profile across tissues. This profile was then correlated with the corresponding\ngene expression profile to infer enhancer-promoter relationships. We have revised the text to clarify\nthis definition in the revised manuscript.\n\nReviewer #3 (Remarks to the Author):\nThis paper by Xie et al presents a very useful data set on the functional annotation of promoters and\nenhancers in sheep by analysing a range of techniques on 24 tissues from a single ewe, including\nChIPSeq, ATAC Seq, CAGE Seq, methylation studies by WGBS and RRBS, and transcript analysis. These\ndata are highly complementary for characterising promoters and enhancers and their functionality in\ndifferent tissues. We would also have liked to have HiC data to reinforce the study of the enhancer-\npromoter link, but it is likely that the samples required for this were not available. This would have\nstrengthened the relevance of the promoter-enhancer link, since you are taking into account distances of\nup to 2 Mb.\nEven if we can always expect more, globally I find this paper very complete and well written. For me it\ndeserves publication in Nature Communications.\nResponse: We sincerely thank the reviewer for the encouraging and constructive comments. We agree\nthat Hi-C data would further strengthen the analysis of enhancer–promoter interactions. In our study,\nthe individual used for multi-omics profiling is the same animal as that used to generate the reference\ngenome20, and we indeed have Hi-C data available for this sample. However, this dataset represents\nonly one tissue with ~30× coverage (https://www.ncbi.nlm.nih.gov/sra/?term=SRX3399085), and thus\nis insufficient to support a comprehensive, tissue-resolved analysis across all 24 profiled tissues.\nTo address this limitation, we employed an integrative strategy to identify enhancer-promoter links.\nThe ±1 Mbp window around transcription start sites was used only as an initial search space to\ncapture potential distal regulatory elements. The final enhancer–promoter interactions were defined\nnot by genomic distance alone, but by integrating three complementary lines of evidence: (i)\ncorrelation between enhancer activity and gene expression across 24 tissues, (ii) CTCF-mediated\nchromatin loop information, and (iii) promoter-enhancer associations inferred from CAGE data. Only\ninteractions supported by these orthogonal criteria and passing an FDR < 0.05 were retained, ensuring\nhigh confidence in their biological relevance. We agree that future availability of high-resolution Hi-C\ndata across additional tissues will further refine this framework.\n\nGeneral Comments :\nI would have appreciated it if the results presented had focused a little more on the functions of interest\nin ruminants. If you could add something to tell us a story about this, either specific to sheep or common\nto sheep and cattle, I would find that interesting.\nResponse: We thank the reviewer for this insightful suggestion. To further highlight the biological\nrelevance and functional insights of the identified enhancers and promoters in ruminants, particularly\nin sheep, we integrated our enhancer-promoter interactions with significant GWAS and eQTL variants\nfrom a recent study of 3,007 sheep29 (https://doi.org/10.1186/s40104-025-01251-4). The most\nsignificant SNP (P = 2.23 × 10⁻⁷) was located within a cerebellum-specific enhancer predicted to\nregulate the XKR4 gene. In the same GWAS study, XKR4 had previously been associated with body\nweight (BW) based on an intragenic SNP. Our analysis not only confirms this BW association but also\nidentifies an additional variant within a regulatory enhancer controlling XKR4, providing novel\nevidence that variation within regulatory elements contributes to phenotypic differences in sheep.\nMoreover, among 181 previously published significant eQTLs, we identified one variant significantly\nassociated with LOC121819880 expression (P = 4.69 × 10⁻⁸), which resides within its interacting\nenhancer (chr6:35907531–35908224). This finding further supports that sequence variants within\nenhancers can directly influence gene expression, reinforcing the accuracy of our enhancer\nannotation. Collectively, these examples demonstrate how incorporating enhancer information\nenhances the biological interpretation of gene regulation in sheep and provides new insights into\nruminant-specific functional mechanisms.\n\nYou have a dataset for which you have the complete genomic sequence of the individual being analysed.\nYou could use the RNAseq data to see if you have allele specific expression ... and if that allele specific\nexpression can be linked with a possible polymorphism in an enhancer or promoter region of that same\ngene, with information by other techniques, showing that this promoter or enhancer is non functional in\nthe studied tissue. It would be great.\nResponse: We sincerely thank the reviewer for this insightful suggestion. Indeed, integrating enhancer\npolymorphisms with allele-specific expression (ASE) would be the most direct and ideal approach to\nvalidate enhancer-gene interactions. However, all 24 tissues in our study were collected from the\nsame individual sheep, which also serves as the reference genome for the current Ovis aries\nARS_UI_Ramb v2.0 assembly. As a result, variants detected from the ChIP-Seq, CAGE, and ATAC-Seq\ndata are likely insufficient to capture genome-wide genetic variation. Therefore, it is not feasible to\ndirectly assess enhancer or promoter polymorphisms within this individual or to correlate them with\nASE derived from the RNA-seq data.\nNevertheless, this concept can be indirectly validated at the population level through expression\nquantitative trait loci (eQTL) analyses. By associating SNP polymorphisms with gene expression across\na population, eQTL mapping allows us to identify variants strongly linked to transcriptional changes.\nWhen such SNPs are located within enhancer or promoter regions, they provide strong evidence that\nsequence variation within cis-regulatory elements influences gene expression.\nTo this end, we examined 181 previously reported eQTLs identified in the sheep population and\nintegrated them with our enhancer–promoter interaction data. Notably, we detected one variant\nsignificantly associated with LOC121819880 expression (P = 4.69 × 10⁻⁸), which resides within its\ninteracting enhancer (CM028709.1:35907531–35908224). This observation provides additional\nevidence that sequence variants within enhancers can directly modulate gene expression, thereby\nreinforcing the accuracy and biological relevance of our enhancer annotations.\n\nI was lost by the numbering of several additional files tables. Can you check numbering when you cite a\nsupplementary data file please. It would facilitate to name the excel file tab with the number of the\ncorresponding supplementary table number\nResponse: We appreciate the reviewer’s helpful comment. We have carefully checked and corrected\nthe numbering of all supplementary tables and renamed each Excel tab according to its corresponding\nsupplementary table number.\n\nTo make it easier to understand the contribution of each of the techniques, could you propose a figure\n(possibly in supplementary material) showing the respective contributions of ChIPSeq with the 4 histone\nmarks, ATAC Seq, CAGE-Seq, RRBS and WGBS, as well as RNA Seq? This would better highlight the\nrespective contribution of each of these different techniques.\nResponse: Thank you for this helpful comment. To better illustrate the respective contributions of\neach technique, we have added a new figure in the revised supplementary materials (Supplementary\nFig. S1), which summarizes the roles of ChIP-seq (H3K27ac, H3K4me3, H3K4me1, and H3K27me3),\nATAC-seq, CAGE-seq, RRBS/WGBS, and RNA-seq in our analysis. The Supplementary Fig. S1 have been\nincluded to the revised manuscript.\nFigure S1. Summary of the respective contributions of ChIP-seq, ATAC-seq, CAGE-seq, RRBS/WGBS, and RNA-seq\nin this study.\n\nSee how you can make your data available more easily than only in supplementary files\nResponse: Thank you for this valuable suggestion. The data of this study have been deposited in the\nFigshare repository (https://figshare.com/account/mycontent/projects). The data are currently under\nrestricted access and will be publicly available upon publication of the manuscript.\n\nDetailed comments :\nIt would be a good idea to state in the introduction, line 133, that the genome of the animal you have\nused has been sequenced and assembled to produce a de novo reference genome.\nResponse: Thank you for this helpful suggestion. We have added a statement in the Introduction to\nindicate that the sheep used in this study is the same individual whose genome was sequenced and\nassembled to produce the current ovine reference genome (ARS_UI_Ramb_v2.0).\n\nFigure 1B :\nThe color scale does not clearly show the number of tissues analyzed for each of the techniques used.\nDifferentiating WGBS and RRBS on 2 different tracks\nI wonder if a 2-colour scale wouldn't be more visual and informative. The same color scale could also be\nused for all the applications used.\nResponse: We have revised Figure 1B by differentiating WGBS and RRBS on separate tracks and using\na clearer two-color scale.\n\nFigure 1C & figure S1 :\nIt would be desirable for the 4 categories TSS, exons, introns, intergenic to be in the same order ... and\nwith similar colors\nIn Figure S1, the color of the green is not the same for the graph showing the percentage and that\nshowing the average over the 24 tissues. It is not so important, but it would be nicer.\nResponse: In the revised manuscript, the previous Figure S1 has been renumbered as Figure S2. We\nhave revised Figure S2 to ensure consistent order and color schemes across all categories in Figure 1C,\nand unified the color tones as suggested.\n\nFigure 1E :\nThe scale, which is neither log nor standard, does not appear to be the best. The SD value for the poised\nenhancer, for example, is not presented on the same scale for the minimum and maximum values.\nResponse: We thank the reviewer for this valuable comment. In the revised Figure 1E, the x-axis has\nbeen unified to a log₂ scale, and both the mean and standard deviation (SD) values are clearly\nindicated. The updated Figure 1E has been included in the revised manuscript.\nFor a better understanding, perhaps write also the % of the genome occupied by each of the states.\nLine 179 : it seems that the figure of 27.1% corresponds to the sum of the repressed polycomb and weak\nrepressed polycomb states, and not simply to repressed polycomb.\nResponse: Thank you for this helpful comment. The value of 27.1% indeed corresponds to the\ncombined proportion of the repressed polycomb and weak repressed polycomb states, and we have\nupdated this in the revised manuscript.\n\nFig 2C :\nLegend of y axis : number of enhancer  number of enhancers\nLines 206-208 : These numbers might be added on fig 2A\nResponse: We have revised these in Fig.2.\n\nFig S3B, S4C & S4D : On y axis, it is numbers rather than frequency\nResponse: We have revised these in the corresponding figures.\n\nA question, seeing fig 2C & D: Would a heatmap showing the number of enhancers and promoters for a\npair of tissues give anything interesting to see in terms of the use of these promoters?\nResponse: We appreciate the reviewer’s helpful suggestion. We performed an heatmap to visualize\nthe shared intersections of enhancers and promoters between pairs of tissues (Fig. 2E and\nSupplementary Fig. S4C). The results revealed that enhancers exhibit greater tissue diversity, whereas\npromoters are more conserved across tissues. This analysis provides additional insight into the tissue-\nspecific regulatory diversity.\n\nIn Figure 2B, rather than simply showing the number of enhancers and promoters identified in each\ntissue, would it be interesting for each tissue to distinguish the number of promoters specific to that\ntissue and promoters observed in more than one tissue? The number of promoters and enhancers could\nthus be represented, highlighting their specific or more or less ubiquitous characteristics.\nResponse: We thank the reviewer for this helpful suggestion. We have updated Fig. 2B to include\nadditional bars representing tissue-unique enhancers and promoters. Furthermore, we added a new\nsupplementary table (Table S7) listing the numbers of tissue-specific and shared enhancers and\npromoters across all 24 tissues. These additions provide a clearer overview of both unique and\nubiquitous regulatory elements among tissues.\n\nLines 228-229 : Of the 9 genes indicated in this table, 4 show a comparison with a sheep gene. Is it 9\ngenes or 5 genes? In any case, Ovis Aries is not mentioned in the species indicated on lines 228-229.\nResponse: We thank the reviewer for the helpful comment. The sentence has been revised for\nclarification. Specifically, nine tissue-specific novel TSSs were aligned with annotated gene sequences,\nof which five were annotated in other species (Bos taurus, Capra hircus, and Cervus elaphus)\n\nLine 239 : in any tissue  in any tested tissue.\nResponse: We have revised this in the revised manuscript.\n\nLine 243: Since you are using TPM for the first time, you should define it.\nResponse: We have included this in the revised manuscript.\n\nFigure 3A & 3B:\nAre you sure that TPM reaches 10^5, which means that for some tissues 10-20 % transcripts are from\nthe same gene ?\nResponse: We thank the reviewer for the thoughtful comment. We have rechecked the data, and the\nmaximum TPM is 17,770.03. The apparent higher values were a visual artifact caused by the log₁₀\nscale on the y-axis. We have changed it to a log₂ scale in the revised Fig. 3A.\n\nAs you considered that a gene is considered as expressed if TPM is > 0.5, it means that quite all genes\nwithout enhancer are in fact NOT EXPRESSED. Am I wrong ?\nResponse: We thank the reviewer for the insightful comment. Genes with TPM ≥ 0.5 were classified\nas expressed only in the quantitative model analyzing enhancer-promoter expression correlations.\nHowever, in subsequent qualitative comparisons of gene expression with and without identified\nenhancers across different tissues (e.g., Figures 3A and 3B), no filtering was applied to avoid potential\nbias introduced by predefined thresholds.\n\nOn Figure 3B, we see a significant number of colored dot for genes with a TPM <0.5, and thus considered\nas non expressed, but you write line 247 that “these non expressed genes lacked an associated\nenhancer”. So you should not have colored dots below a TPM of 0.5. Am I wrong ?\nResponse: We thank the reviewer for the valuable comment. The threshold of TPM ≥ 0.5 was only\napplied in the quantitative model for enhancer-promoter correlation analysis. In Figure 3B, no\nexpression threshold was applied, and all genes were shown to illustrate the overall distribution of\nenhancer-associated and non-enhancer-associated genes.\n\nLine 260 : 51.52% of genes, or 51.52 % of promoters ?\nLines 260 & 261 : 52% not 51.52%. 89% not 88.53% (give fewer numbers)\nResponse: Thank you, we have revised these in the revised manuscript.\n\nFigure 3F:\nI wonder if a 2 colors scale can be more visual and informative\nEven if you used a log scale, numbers are not in log2 as indicated in the title, as you indicate 2^5 and\n2^10. So remove “(log2)”.\nResponse: We thank the reviewer for the helpful suggestions. The heatmap represents the number of\npromoters associated with each enhancer hub (and the number of enhancers in each promoter hub),\nwhere values are either zero or greater than zero. Therefore, a two-color scale from white to yellow\nwas used to better reflect these differences across tissues. In addition, “(log₂)” has been removed in\nthe revised figure.\n\nFigure 4A:\nIt is not clear for me. I understand that 281 + 151 = 432, but I do not understand.\nWhat means “specific genes from promoters”\nTo what correspond these 151 genes in light orange ?\nResponse: We thank the reviewer for the helpful comment. In Figure 4A, 281 genes were identified\nbased on tissue-specific promoters alone, while 432 genes were obtained from tissue-specific\nenhancer-promoter combinations. Among them, 151 genes (in light orange) were specifically\nidentified from enhancer-promoter combinations but not from promoter-based analysis. These genes\nlikely represent cases where the same promoter is active in multiple tissues, but different enhancers\nregulate the gene in a tissue-specific manner, thereby revealing additional enhancer-driven tissue\nspecificity.\nThe tissue-specific enhancer-promoter combinations provide additional insights, as they incorporate\ninteractions of enhancers and promoters that are not captured when considering tissue-specific gene\nexpression alone. While previous analyses have focused on promoter-gene specificity derived from\ngene annotations, the inclusion of the previously unrecognized regulatory roles of enhancers from\nenhancer-promoter combinations uncovers new insights into gene regulation.\nIn addition, we performed GO enrichment analysis specifically for the 151 genes uniquely associated\nwith enhancers in the cerebellum. This analysis revealed five significantly enriched biological\nprocesses: positive regulation of synaptic transmission, synapse assembly, positive regulation of\nneuron differentiation, modulation of chemical synaptic transmission, and regulation of neuron\ndifferentiation (Supplementary Table S14). Importantly, we identified 11 key genes that were\nsupported by two or more of these processes. Among them, NLGN3, BDNF, and NRCAM appeared in\nthree pathways (Supplementary Figure S8), while RIMS2, RIMS1, SYT1, SLC8A2, NEUROD1, NEUROD2,\nBRINP1, and VWC2 were shared by two pathways. These genes represent core neuronal regulators\nwith pleiotropic roles in synaptic transmission, neuronal differentiation, and circuit assembly. Their\nrecurrence across multiple enriched categories underscores that enhancer-specific genes in the\ncerebellum converge on fundamental mechanisms of neuronal development and synaptic function,\nthereby reinforcing both the robustness and biological relevance of our findings. We have included\nthese in the revised manuscript.\n\nFig 4C :\nYou indicate in the text a single isoform in the cortex, but you indicate 2 and 4 in the RNASeq track,\nsuggesting 2 different isoforms exists.\nResponse: We thank the reviewer for the helpful comment. The numbers “2” and “4” in the RNA-seq\ntrack represent junction reads rather than distinct isoforms. Given the small number of supporting\nreads and the annotation information indicating that both junctions point to the same exon 1, these\nsignals likely represent technical variation rather than alternative splicing. So we considered them to\nbelong to a single isoform.\n\nFig 5 :\nI wonder if using a log scale for the y axis is relevant (moreover you show 0 on a log scale). Have you\ntested a standard scale? It might better put in evidence the slope of the correlation.\nResponse: We thank the reviewer for the helpful suggestion. We tested a standard (linear) scale and\nexamined an enlarged view of the 0-100 TPM range; however, the log₁₀ scale provided a clearer\nrepresentation of the overall trend. In the revised figure, the y-axis represents true TPM values, and\nTPM values of zero were converted back to their original values (not log-transformed) to avoid\ndistortion near the lower range.\nStandard scale of enhancer-gene expression and methylation level in lung tissue:\nEnlarged view of the 0-100 TPM range:\n\nLine 382 : There is no fig 5C & 5D, you have 2 figures on 5A and 2 on 5B. Alternatively change the figure\nand its legend.\nResponse: Thank you, we have revised this.\n\nLine 404 : I did not find table S16. I should say that the supplementary tables are not properly numbered,\nat least in the version I downloaded for review. For a simpler reading, the tab title of Excel files should\nindicate the table number. I did not also find the legend of supplementary tables.\nResponse: We thank the reviewer for the comment. The supplementary tables have been reordered\nand renamed according to the suggestion.\n\nLine 407 : the value of 24467 enhancers in the cortex do not correspond the value on fig S10A, which is\n42467. So in the text 24467  42467. With this value, the indicated percent is correct.\nResponse: Thank you, we have corrected this.\n\nLine 406 & 410 : use fewer numbers to show the percentages\nResponse: Thank you, we have corrected this.\n\nLines 412-413 : If you intend to say that sheep is closer to cattle than to pig, you should say that on fig\n6B, sheep is even closer to human than to cattle, which is funny ! Is it due to a really larger number of\npromoters put in evidence in human (as suggested on fig S11), which may induce that a larger number of\npromoters identified in sheep can also be found in human ?\nResponse: We thank the reviewer for this interesting observation. Indeed, in Figure 6B, sheep appears\nslightly closer to human than to cattle. This pattern is likely influenced by the larger number of\npromoters identified in human, as shown in Figure S12, which increases the number of shared\npromoters when comparing with other species. We have clarified this point in the revised manuscript.\n\nFor figure S11, I wonder if it wouldn't be more informative simply to have (or in addition to) a table with\nthe 5 species in rows and columns, and at the intersection of the rows and columns the percentage of\npromoters (or enhancers) of the species in that column found in the species in the row. This would\nbetter highlight the proximity between species :\nNote : the array will not be symmetrical because the number of promoters found in the two species\nunder consideration is not the same.\nResponse: Thank you for the helpful suggestion. We have added an additional supplementary\nTableS19 in the revised manuscript.\n\nLine 414 : you may add the percentage of enhancers and promoters found in all 5 species. Based on your\ndata it is on average for the 4 tissues around 0.8 % for enhancers and 17% for promoters.\nResponse: Thank you, we have included these in the revised manuscript.\n\nFig 6D : it would be easier to have the tissue in the same order in the text and in the figure, whereas we\nhave cerebellum, cortex, lung and muscle in the text, but cerebellum, lung, muscle and cortex on the\nfigure.\nResponse: Thank you, we have revised the Fig6D.\n\nLines 431-433: I think that in the discussion part, you indicate 69 pairs of enhancer-promoter conserved\nin cattle and sheep and not conserved in other species. If I am right, this number can be indicated also in\nthis result part..\nResponse: Thank you for the helpful suggestion. We have added the corresponding information and\none supplementary table (Table S21) in the Results section of revised manuscript.\n\nLine 445-446 : You did not explain here or in M&M how these SNPs have been identified. This greatly\nreduces the scope of your conclusions, and in particular the frequency of presence of the polymorphism\nin the various compartments (exon, intergenic region, enhancers, promoters, etc.) which may be\ninfluenced by the way in which the polymorphism was identified. In addition, was the polymorphism\nidentified following sequencing of the ewe, or of a large number of individuals from different breeds?\nResponse: Thank you, we have included these in the M&M section of the revised manuscript.\n\nFig 7a : If you used density expressed in % on the legend of y axis, replace 0.00022 by 0.022 on the figure\nResponse: Thank you, we have corrected this in the new figure.\n\nLine 473 : such as; enhancers  such as enhancers\nResponse: Thank you, we have corrected this.\n\nLines 578-580 : could you check whether this number of 69 conserved enhancer-promoter pairs is higher\nthan the number you would expect at random between two species taken at random from the 6 species\nconsidered.\nResponse: We thank the reviewer for this constructive suggestion. The 69 conserved pairs refer to 38\nenhancers and 31 promoters showing reciprocal pairing in ruminants, rather than complete\nenhancer–gene matches. Thus, a direct comparison with randomly paired species would not be fully\nappropriate in this context.\n\nLine 594 : For the QTL, I find the word « associated » too strong, as it suggests causality. Presently it is\na colocalization. You also say : “this study successfully identified 38 QTL associated with …”. I would\nsay in the other direction that you demonstrate the presence of promoters or enhancers in regions\nwhere QTL have been identified.\nResponse: We appreciate the reviewer’s valuable comment. We agree that our findings represent\ncolocalization rather than causality. In the revised manuscript, we have replaced “associated” with\n“colocalized”, and we rephrased the sentence “this study successfully identified 38 QTL associated\nwith …” to “this study identified 38 QTLs colocalized in promoter regions of 17 genes and 215 QTLs\noverlapping to 214 enhancers…” This wording better reflects the current evidence and avoids\noverinterpretation. We have revised the relevant text in both the Results and Discussion sections of\nthe revised manuscript.\n\nM&M :\nIn tissue collection, you may remember that the genome of this ewe has been de novo sequenced.\nIndicate also the ref of the publication (29).\nResponse: Thank you, we have included this.\n\nLine 734 : As you reported the number of RNA seq reads per tissue of 65-125 Millions, you might add in\nthe text, that it means that for most tissues, fragments of genes considered as expressed are found at\nleast 30 fold.\nResponse: We thank the reviewer for the helpful suggestion. We have added a clarification in the text\nindicating that this sequencing depth corresponds to an average coverage of at least 30-fold for most\nexpressed genes.\n\nLines 766 & 767, if you use _Ti underscript to indicate the tissue, you should do the same for the 4 types\nPH3K4me1 , PH3K27ac, PCAGE-bid, PATAC.\nResponse: Thank you, we have revised these in the revised manuscript.\n\nLine 773 : represent  represented ?\nResponse: Thank you, we have corrected this.\n\nLine 795 :\nI am not sure I understand how you calculate Sm & Em. You use Sn, but n is the last tissue. May be you\nshould use an other letter than n. In fact Sn+1 does not exist.\nBe careful you have large s and e for n odd, and small s and e for n even.\nIf the formula (s_(n/2)+s_(n/2+1))/2 is correct, it could be more readable writing (s_(n/2)+s_((n+2)/2))/2\nor (s_(n/2)+s_((n/2)+1))/2\nResponse: Thank you, we have corrected this.\n\nLine 804 you indicate ± 1 Mbp of TSS, but in line 545 you indicates “2 Mbp range of promoter”\nResponse: Thank you, we have corrected this.\n\nLine 859 : can you indicate the source of this collection of SNP ?\nResponse: Thank you, we have revised this."}, {"round": 2, "reviewer_comments": {"reviewer_1": "In their revised version of the manuscript, the authors have adequately addressed all my previous concerns which has resulted in a much clearer and stronger manuscript. As I already stated in my previous report, this comprehensive collection of enhancers and promoters will be an important resource for studies aimed at bridging the genome to phenome gap in sheep.\nThe added descriptions of the figures and supplementary tables have greatly improved the readability of the manuscript. I have no further comments for the revised version and support publication in its current form.", "reviewer_2": "The authors have extensively added new data and analyses, and largely addressed my concerns. I have no more question.", "reviewer_3": "The new version of this paper takes into account the comments. For me it deserves publication in Nature communication.\nI have a single remaining comment concerning the notation which could be clearer.\nIn line 836-837, Penh_Ti is not a potential enhancer, but a peak corresponding to a potential enhancer. I mean that you distinguished Penh_Ti and EnhTj. Psomething is not something , it is a peak revealing something (for instance a potential enhancer).\nThe notation is ever confusing on lines 831-832 and 838-839. On line 838 for instance you indicate:\nI guess that the one I marked as 2, corresponds to what you define in line 831. Thus, you should find a way to differently indicate the “candidate enhancer” (and promoter) defined in lines 831 & 832, and the “potential enhancer” (and promoter) defined in lines 838 & 839."}, "author_response": "Response: We sincerely thank the reviewer for the positive and supportive comments.\nResponse: We sincerely thank the reviewer for the thoughtful comments.\nResponse: We sincerely thank the reviewer for the valuable and thoughtful comments. As noted by the reviewer, the definition marked as 2 corresponds to the candidate enhancer defined in lines 831–832 (Penh = P(H3K4me1) ∪ P(H3K27ac)). To clearly distinguish candidate regulatory regions from potential enhancers and promoters, we have revised the notation and now denote candidate peak regions in revised manuscript as Peak and Peak . And the P and P in lines 831 and 832 enh_Ti prom_Ti enh prom were also revised as Peak and Peak , respectively. enh prom\nWe greatly appreciate this helpful suggestion.\nIn lines 844-845 you indicate that i is for a given tissue, but you do not define to what corresponds j as used in Enhij. I imagine that j is for the j-th enhancer on tissue i ?\nResponse: We thank the reviewer for pointing out this lack of clarity. Yes, j denotes the j-th enhancer (or promoter) from tissue i. We have included this in the revised manuscript by the following definition: “where j is the j-th enhancer or promoter from tissue i (Ti), respectively”.\nIn line 848, you write . But t(.) is only defined in line 855, too far from its first use.\nResponse: We thank the reviewer for this helpful suggestion. We have moved the definition of t(·) to immediately follow its first occurrence in the equation, and this change has been implemented in the revised manuscript.\nIn fact for enhancers, you use 5 different notations: Penh, Penh_Ti, Penh_all, Enhij, EnhTi , and it is a bit complicated to follow precisely these 5 different kinds. Please make that part as clear as possible with the adequate notation.\nI am wondering, if it would not be simpler to add a table telling which kind of signals you need to identify a possible enhancer or promoter. As far as I see, it can be something like:\nResponse: We thank the reviewer for this helpful suggestion. We have summarized the integration of multi-omics data to define potential enhancers and promoters in Supplementary Table 6, and the relevant text has been included in the revised manuscript.\nI paste a PDF file easier to read.\nResponse: We greatly appreciate the reviewer’s careful review and the annotated PDF, which was very helpful."}]

Decision letter