We downloaded H3K27ac and H3K4me3 ChIP-seq and RNA-seq data for six pig tissues (Figures 1A,B), i.e., adipose, skeletal muscle, brain cortex, liver, lung, and spleen, publicly. For human and mouse, we downloaded the available H3K27ac and H3K4me3 ChIP-seq and RNA-seq data for tissues such as adipose, brain (forebrain or brain cortex), muscle, liver, lung, and spleen (Figures 1A,B). The individual accession IDs are listed in Supplementary Table S1.

For pig ChIP-seq data, the fastq raw files were trimmed by using Trim_galore (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) with default parameters in single-end mode. Thus, reads with low quality (Q < 20) and shorter than 20 bp were removed. Then, high-quality reads were aligned to the pig reference genome (susScr11) using bowtie2 (Langmead and Salzberg, 2012) with default parameters. Reads with low-mapping quality scores (MAPQ <30), unmapped, and duplicated were filtered (-F 1796) using samtools (Li et al., 2009). The filtered reads were used for downstream analyses. For human and mouse ChIP-seq data, uniquely mapped reads were directly recovered from ENCODE, which were aligned to the human and mouse genomes (hg38 and mm10), respectively. H3K27ac and H3K4me3 peaks were identified using macs2 (Zhang et al., 2008) with parameters “--keep-dup = auto”, while “-g” was set to “2.7e9”, “hs”, and “mm” for pig, human, and mouse tissues, respectively. For tissues with two replicates, their peaks overlapped at least 50% reciprocally and were used for downstream analyses using bedtools (Quinlan and Hall, 2010) with parameters “intersect -f 0.5 -r”. For tissues without replicates, their top 20,000 highest signals regions were used for downstream analyses. For analyses on H3K4me3 peak signal intensity, similar to Benayoun et al. (2014), pile-up of ChIP-seq reads mapped within the significant macs2 peak was computed as the coverage in the specified intervals using function “coverage” of bedtools with default settings and normalized to the input signal in the same interval. Then, H3K4me3 peak signal was normalized to peak breadth to obtain tags per base pair of peak.

For pig RNA-seq data files, adapters were trimmed similar to ChIP-seq data but with paired-end mode (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/). Similar to Kern et al. (2021), reads were aligned to pig genome using STAR (Dobin et al., 2013) (--outFilterMultimapNmax 20 --alignSJoverhangMin 8 --alignSJDBoverhangMin 1 --outFilterMismatchNmax 999 --alignIntronMin 20). Meanwhile, reads with low-mapping quality scores (MAPQ <30) were filtered using samtools (Li et al., 2009). For human and mouse RNA-seq, reads that uniquely aligned to hg38 and mm10 genome were directly recovered from ENCODE. Gene counts were calculated based on bedtools (Quinlan and Hall, 2010) with function “multicov”. The gtf annotation files of pig (Suscra 11), human (hg38), and mouse (mm10) (Ensemble genes annotation) were downloaded from UCSC table browser (http://www.genome.ucsc.edu/cgi-bin/hgTables). Gene expressions were normalized to reads per million mapped reads per kilobase (rpkm) with house script.

SEs were identified based on H3K27ac enriched peaks using ROSE (https://bitbucket.org/young_computation/rose) with default parameters. Broad H3K4me3 domains were defined as the top 5% broadest H3K4me3 enriched peaks (Benayoun et al., 2014). A SE or BD was assigned to the gene with the closest transcription start site using R package ChIPseeker (Yu et al., 2015). GO analyses were performed using R package ClusterProfiler (Yu et al., 2012). GO terms with p values <0.001 were considered significantly enriched. The calculation of tissue specificity of SEs and BDs and their typical counterparts was based on the function “mergePeaks” in Homer (Heinz et al., 2010), and parameter “-d” was set to “2000”, “4,000”, and “8,000”, respectively.

For the orthologous SE and BD comparisons (Figure 1E), SEs and BDs for different tissues in pig (sucScr11) or mouse (mm10) were converted into human (hg38) coordinates by using UCSC LiftOver tools (https://genome.ucsc.edu/cgi-bin/hgLiftOver), the minimum match of which was set to 0.2. The identification of functionally conserved SEs and BDs was based on the function “mergePeaks” in Homer (Heinz et al., 2010) using default parameters. Sequence conservation scores were calculated based on the vertebrate conservation phastCons tracks (phastCons100way) from UCSC table browser.

We first investigated the characteristic features of SEs in six pig tissues, i.e., adipose, skeletal muscle, brain cortex (brain), liver, lung, and spleen. For the comparisons, we also assessed the characteristic features of SEs in six human and mouse tissues, i.e., adipose, skeletal muscle, brain, liver, lung, and spleen. We identified a median of 946, 642, and 616 SEs for pig, human, and mouse, respectively (Figure 2A). Pig SEs (in median, 3,241 bp) were longer than TEs (in median, 847 bp) in each tissue (Wilcoxon rank-sum test, p < 0.001) (Figure 2B and Supplementary Table S2). The H3K27ac signal levels of pig SEs (in median, 13,149 rpm) were higher than TEs (in median, 946 rpm) in each tissue (Wilcoxon rank-sum test, p < 0.001) (Figure 2C and Supplementary Table S2). Furthermore, expression levels of pig SE-associated genes (in median, 8.6 rpkm) were higher than TE-associated genes (in median, 2.4 rpkm) (Wilcoxon rank-sum test, p < 0.001) (Figure 2D and Supplementary Table S2). The results suggest that the characteristic features of SEs in pig are similar to human and mouse in our study and previous studies (Hnisz et al., 2013; Whyte et al., 2013).

We then examined the characteristic features of BDs in six pig tissues and compared them with six human and mouse tissues. We identified a median of 1,582, 1,042, and 942 BDs for pig, human, and mouse (Figure 2E). Similar to human and mouse, pig BDs (in median, 2,134 bp) were longer than narrow peaks (NDs) (in median, 814 bp) in each tissue (Wilcoxon rank-sum test, p < 0.001) (Figure 2F and Supplementary Table S3). In addition, the median peak signal intensity (tags per bp) of BDs was higher than NDs in pig, human, and mouse tissues (Figure 2G and Supplementary Table S3). Peak signal intensity was measured as the normalized number of mapped ChIP reads per base pair for each peak determined by macs2 (Benayoun et al., 2014). Values were normalized to the sequencing depth of the ChIP sample and peak breadth. Meanwhile, the levels of BD-associated genes (in median, 6.9 rpkm) were higher than NDs (in median, 3.0 rpkm) (Wilcoxon rank-sum test, p < 0.001) (Figure 2H and Supplementary Table S3). This indicates that there are no consistent relationships between H3K4me3 signal levels and width.

To compare the tissue specificity between SEs and TEs, all elements were unified to the size of 2 kb according to peak median location. Across six pig tissues, about 57–83% SEs and about 32–66% TEs were tissue specific (paired Wilcoxon test, p < 0.05) (Figure 3A). Consistently, across six human tissues, about 64–88% SEs and about 32–66% TEs were tissue specific (Figure 3B) (paired Wilcoxon test, p < 0.05). Across six mouse tissues, about 79–90% SEs and about 53–77% TEs were tissue specific (paired Wilcoxon test, p < 0.05) (Figure 3C). When elements were unified to 4 and 8 kb, it also showed higher tissue specificity in SEs than in TEs (paired Wilcoxon test, p < 0.05) (Supplementary Table S4). This indicates that SEs exhibit higher tissue specificity than TEs in pig, human, and mouse.

When BDs and NDs were unified to 2 kb, about 14–31% BDs and about 4–24% NDs were tissue specific (paired Wilcoxon test, p < 0.05) (Figure 3D). At the same time, about 26–41% BDs and 1–24% NDs were tissue specific for six human tissues (paired Wilcoxon test, p < 0.05) (Figure 3E). In addition, about 15–43% BDs and 3–14% were tissue specific for six mouse tissues (paired Wilcoxon test, p < 0.05) (Figure 3F). When BDs and NDs were unified to 4 and 8 kb, the frequency of the tissue-specific BDs were also higher than NDs (Supplementary Table S4) (the paired Wilcoxon test, p < 0.05). This suggests that tissue specificity of BDs is higher than NDs in three species.

To assess whether SEs and BDs were associated with tissue identity, we classified the top 20 most significant enriched GO terms of their associated genes into seven categories (C1–C7), respectively. Each category represented tissue-specific function (Figure 4A and Supplementary Tables S4, S5). In brief, C1 corresponds to neuron, axon, and synapse-related processes. C2 contained muscle, actin, and myosin-related processes. C3 consisted of sterol metabolic, organic anion transport, insulin-related processes. C4 was composed of immune-related processes. C5 represented processes such as response to stimulus. C6 contained processes such as cell migration, adhesion, and assembly. C7 contained tube and vasculature-related processes. The details of enriched terms in each category are provided in Supplementary Table S5. For pig cerebellum, brain cortex, and hypothalamus, they showed similar enriched terms such as neuron, nervous, axon, and synapse-related processes (Supplementary Table S5). For the comparisons, brain cortex was selected at random for the SE and BD comparison with human and mouse brains.

We then used these categories (C1–C7) to explore whether SEs were involved in tissue identity–related processes in each tissue (Figures 4A–C and Supplementary Table S5). Except for pig lung, in five out of six tissues among three species, GO annotations showed SE-associated genes were enriched terms relating to tissue identity. For example, in brain, they were related to neuron-related processes (C1); in muscle, they were associated with muscle and myosin-related processes (C2); in liver, they were correlated with small molecule catabolic processes (C3); in spleen, they were correlated with immune-related processes (C4); and in adipose, they were enriched in cell mobility and cell migration–related processes (C6) (Figure 4A). In pig lung, the top 20 terms were largely distributed in immune-related processes (C4), while in human and mouse lung, they were mainly distributed in cell mobility, cell migration (C6), and tube and vasculature-related processes (C7). Moreover, in mouse lung, SE-associated genes were related to respiratory system/tube development and lung development. Our results indicate that SEs are preferable to define tissue identity in most tissues of three species except for pig lung.

At the same time, we investigated whether BDs were related to tissue identity in each tissue based on the seven categories (C1–C7) (Figures 4D–F). BDs were well coupled with tissue identity in three out of six pig tissues, i.e., brain, liver, and adipose. For human BDs, they had connection with tissue identity in four human tissues, i.e., brain, muscle, spleen, and adipose. For mouse BDs, they were also capable of defining four tissues, i.e., brain, muscle, spleen, and adipose. This demonstrates BDs have the ability to define tissue identity in majority of tissues in pig, human, and mouse, but show a weaker determination ability than SEs.

To investigate the functional conservation of SEs among species, we analyzed the distribution of orthologous SEs in six tissues among pig, human, and mouse. Across six tissues, about 5–12% (55–182) pig functionally conserved SEs were identified because they had orthologous SEs in human and mouse (Figure 5A, Supplementary Figure S1 and Supplementary Table S6). Then, we estimated the sequence conservation of these functionally conserved SEs (Figure 5B). In pig brain, muscle, and spleen, phastCons scores in functionally conserved SEs were significantly higher than un-conserved ones, while the contrast pattern was found in pig adipose, liver, and lung (Wilcoxon rank-sum test, p < 0.001) (Figure 5B). Thus, this implied that functionally conserved SEs were not sequence conserved. Later, we explored the possible biological roles of the pig functionally conserved SEs. We also classified the top 20 most significant terms into seven categories (C1–C7) for six tissues. In four out of six pig tissues, including adipose, muscle, liver, and spleen, these SEs are able to control tissue identity–related genes (Figure 5C and Supplementary Table S5). The result indicates that pig functionally conserved SEs are able to define tissue identity in majority of tissues.

Likewise, to explore the functional conservation of BDs, we analyzed the distribution of BD-associated orthologous genes in six tissues among pig, human, and mouse. About 8–16% (99–309) of pig functionally conserved BDs were determined as they kept orthologous BDs in human and mouse (Figure 5D, Supplementary Figure S2, and Supplementary Table S6). Then, we estimated the sequence conservation of these functionally conserved BDs (Figure 5E). In all tissues, pig functionally conserved BDs showed lower phastCons scores than their un-conserved BDs (Wilcoxon rank-sum test, p < 0.001) (Figure 5E), while in brain and muscle, the patterns were reversed. This suggested that functionally conserved BDs were not sequence conserved in all detected pig tissues, which was similar to pig functionally conserved SEs. Later, we investigated the functional features of the pig functionally conserved BDs by using the seven categories (C1–C7) to summarize the top 20 most significant terms in each tissue. In contrast to SEs, conserved BDs were well relevant to tissue identity in brain and liver (Figure 5F and Supplementary Table S5). This may indicate that the ability of functionally conserved BDs in defining tissue identity is weaker than SEs.