Repressor elements provide insights into tissue development and phenotypes in pigs

INTRODUCTION

The domestic pig (Sus scrofa domestica) is an important source of human nutrition and is also highly utilized as a biomedical model for research (Lunney et al., 2021; Yan et al., 2018). The quantity and quality of pork are governed by a complex network of genomic regulatory elements, including both activators and repressors. While recent studies have primarily focused on activator elements, such as enhancers (Wu et al., 2022) and promoters (Jiang et al., 2024; Pan et al., 2021), repressor elements (REs) also play important roles in the modulation of gene expression and, consequently, phenotype manifestation in pigs. For example, a single G-to-A nucleotide transition in intron 3 of the insulin-like growth factor 2 (IGF2) gene has been identified as the causative mutation for increased skeletal muscle mass and decreased backfat thickness in commercial pig breeds (Van Laere et al., 2003). This mutation disrupts the binding site for the repressor protein ZBED6, leading to a three-fold up-regulation of IGF2 expression in the skeletal muscle of pigs (Markljung et al., 2009). Similarly, a 14 bp insertion in albino pigs has been implicated in the creation of a binding site for SOX proteins, which function as transcriptional silencers, leading to the elimination of an isoform of the MITF gene and causing hearing loss in pigs (Chen et al., 2016). Given the crucial role that repressor elements play in gene regulation, a comprehensive genome-wide identification of these elements in pigs would enhance our ability to identify causal mutations underlying diverse phenotypes across pig breeds.

Traditional Chinese philosophy employs the Yin-Yang concept to describe complex dynamic processes. In 1952, Alan Turing proposed the reaction-diffusion system, a hypothetical molecular mechanism that incorporates both an activator and inhibitor (Turing, 1952). Today, reaction-diffusion models have been widely applied to explain a variety of biological phenomena across different levels of organization, such as the stripe patterns on the skin of marine angelfish (Pomacanthus) (Kondo and Asai, 1995), formation of digit patterns during finger development (Raspopovic et al., 2014), and the regulation of gene expression profiles (Farouq et al., 2021). We propose that at the molecular level, gene expression regulation can also be understood through reaction-diffusion systems that involve both activator elements and REs. Consequently, identifying REs involved in gene regulation is essential for understanding these mechanisms. Histone modifications, including H3K9me3, H4K20me3, and H3K27me3, are closely associated with transcriptional repression. H3K9me3 is established by histone methyl transferases (HMTs) and interacts with heterochromatin protein 1 (HP1) to stabilize heterochromatic structures, particular in tandemly repeated chromosomal regions (Hiragami-Hamada et al., 2016). Recent research has shown that H4K20me3 and the enzyme responsible for its formation, SUV420H2, play critical roles in chromatin compaction during cellular quiescence and differentiation, implicating H4K20me3 in the formation of heterochromatin (Nelson et al., 2016). Bivalent promoter regions marked by both H3K27me3 and H3K4me3 inhibit transcription by directly preventing the release of RNA polymerase II from the transcription start site, a process catalyzed and maintained by Polycomb repressive complexes (PRCs) (Zhang et al., 2020). Furthermore, recent studies have utilized H3K27me3 data to map the genomic locations of regulatory elements (REs). (Cai et al., 2021; Huang et al., 2019; Ouyang et al., 2023).

In this study, we analyzed epigenetic and gene expression data from six pig tissues to identify REs involved in gene regulation based on the negative correlation between H3K27me3-chromatin accessibility and neighboring gene expression. Identified REs (silencing promoters, silencing enhancers, and silencers) were then classified into super repressor elements (SREs) and typical repressor elements (TREs). SREs demonstrated a significantly stronger inhibitory effect on target gene expression compared to TREs. Functional enrichment analysis suggested that REs may play a critical role in tissue differentiation by suppressing the expression of genes associated with specific biological functions. Furthermore, our findings revealed a dual regulatory mechanism across different tissues, where both activator and repressor functions converge to modulate the expression of the KLF4 gene. This gene is known to influence several economically important traits in domestic pigs, underscoring the significance of these repressor elements in agricultural genomics.

MATERIALS AND METHODS

ChIP-seq and ATAC-seq data quality control and preprocessing

Raw ChIP-seq and ATAC-seq data from eight tissues of two 6-month-old domestic pigs were downloaded from the SRA database (PRJNA665212; Supplementary Table S1) (Kern et al., 2021). Adapter sequences and low-quality reads were trimmed using fastp (v.0.23.2) (Chen et al., 2018) with default parameters. Trimmed reads were aligned to the reference pig genome (Sscrofa11.1) using BWA mem (v.0.7.17) (Li and Durbin, 2009) and the resulting SAM files were then processed to remove duplicates using the “markdup” function in Sambamba (v.0.8.0) (Tarasov et al., 2015). The BAM files were converted to bw files using the “bamCoverage” function in deepTools (v.3.5.1) (Ramírez et al., 2016) for visualization in the IGV browser (Robinson et al., 2011). Based on a visual inspection of the peak profiles in IGV, samples with significantly poor quality were excluded.

In addition, processed peak files were obtained from the GEO database (GSE158430). For biological replicate samples from the same tissue, the peaks of the two replicates were first merged using BEDTools (Quinlan and Hall, 2010). The intersection of the merged peaks was then identified, retaining only those with overlapping regions.

RNA-seq data processing

Raw sequencing data were trimmed using Trimmomatic (v.0.39) (Bolger et al., 2014) and then aligned to the reference pig genome (Sscrofa11.1) using Hisat2 (v.2.2.0) (Pertea et al., 2016). Gene expression quantification was performed using StringTie (v.2.2.0) (Pertea et al., 2016), with the gtf file (Sus_scrofa.Sscrofa11.1.104.gtf) from Ensembl used for annotation. All available RNA-seq data were obtained for each tissue, using at least 40 samples for each tissue to ensure the reliability of the results.

RE identification

RE identification was conducted in accordance with previously proposed methods (Huang et al., 2019). (1) Overlapping regions between the chromatin accessible regions (ATAC-seq) and inhibitory regions (H3K27me3) were defined as target regions (H3K27me3-OCRs). (2) Expression levels of target genes were obtained from the samples of the various tissues. For target regions overlapping with genes, the host gene was considered as the target gene. For target regions located in intergenic regions, upstream and downstream genes of the target regions were defined as target genes. (3) The correlation corr(G,h) between H3K27me3-OCRs (h, h=1 and 0 representing the presence and absence of a region from a given tissue, respectively) and the expression level of target genes (G) was calculated using the following formula. Target regions with −0.05<corr(G,h)<0 were identified as REs. The Wilcoxon rank-sum test was used to measure significant differences between h=1 and h=0 in different tissues.

SRE identification

The “super” module in Homer (v.4.11) (Heinz et al., 2010) was used to identify broad peaks marked by H3K27me3, with the maximum distance between flanking regions of two single pre-stitched peaks set to 4kb (Ouyang et al., 2023). H3K27me3 peak clusters with signal slopes larger than 1 were identified as H3K27me3-rich regions, as reported in previous research (Ouyang et al., 2023). REs located within the H3K27me3-rich regions were defined as SREs, while the remainder were defined as TREs.

Peak annotation

Peak annotation was conducted using ChIPseeker (v.1.36.0) (Yu et al., 2015), with annotation information extracted from Sus_scrofa.Sscrofa11.1.104.gtf using GenomicFeatures (v.1.52.1) (Lawrence et al., 2013).

Motif enrichment analysis

A comprehensive collection of 1 284 motifs was retrieved from the CIS-BP database (Weirauch et al., 2014) and employed for motif enrichment analysis. Motifs enriched in RE regions were identified using the “enrich_motifs” function from the R package “universalmotif”, employing the default algorithm. After adjusting the P-value using the “FDR” method, significant enrichment was declared if the P-adjust value was less than 0.05. Heatmap visualizations were generated using the R packages “pheatmap” and “ggplot2”.

Functional enrichment analysis

ClusterProfiler (v.4.8.1) (Wu et al., 2021) was used to perform functional enrichment analysis on the gene sets of the target genes. Annotation information was sourced from org.Ss.eg.db (v.3.17.0) for Gene Ontology (GO) annotation, with only biological process (BP) GO terms retained. The remaining parameters were set to default.

Gene set enrichment analysis (GSEA)

GSEA (Subramanian et al., 2005) was performed using mouse transcriptomics data (Mukhtar et al., 2020). Analysis was performed using the “GSEA” module in ClusterProfiler (v.4.8.1) (Wu et al., 2021) with default parameters. Gene sets were derived from the mouse orthologs, projected from the target genes of the pig REs. A ranked gene list was generated by sorting genes based on their log₂(FC) values derived from differential expression analyses across different developmental stages of the neuronal lineage.

Hi-C data processing

Hi-C data from two Bama pig adipose tissues were downloaded from the SRA database (PRJNA733023; Supplementary Table S1) (Jin et al., 2023). After trimming using fastp (Chen et al., 2018), the Hi-C raw sequencing data were aligned to the pig reference genome (Sscrofa11.1) with BWA mem (Li and Durbin, 2009) using the parameters “-SP -A1 -B4 -E50 -L0”, with all other parameters set to default. The aligned SAM files were converted to pair files using the pairtools (v1.0.2) (Open2c et al., 2023) “parse2” tool, with additional settings “--drop-sam --drop-seq --add-columns mapq –flip”. Sorting was performed using “sort”. The “dedup” tool was used to remove duplicates. Filtered pairs were loaded into a cool matrix using the cooler (v.0.9.1) (Abdennur and Mirny, 2020) “cload pairs” tool. The pairs were also loaded into the “hic” matrix using the “pre” module in juicertools (v.1.22.01) (Durand et al., 2016), while the “hiccups” module was used for loop calling.

Quantitative trait locus (QTL) analysis

QTL data for pigs were obtained from pigQTLdb (Hu et al., 2022). Fine QTLs were defined as QTLs containing a maximum of five genes.

Dual-luciferase reporter assays

REs were amplified using polymerase chain reaction (PCR) from Duroc pig genomic DNA (primers listed in Supplementary Table S2). The amplified genomic fragments were cloned into the KpnI and MluI sites of the pGL3-RPL32 promoter plasmid (Firefly luciferase (fLUC) driven by the human RPL32 gene promoter) and verified by Sanger sequencing. Individual constructs were tested by co-transfecting 2×10⁵ cells with the respective fLUC construct (0.3 pmol) and a Renilla luciferase (rLUC) control plasmid (pRL-TK plasmid promoter was replaced by the pig ACTB promoter; 0.06 pmol) using jetPRIME® transfection reagent (Polyplus, cat. no. 101000046, France). Luciferase reporter gene activity, along with normalized fLUC to rLUC activity, was measured using the Vazyme Dual Luciferase Reporter Assay Kit. (Vazyme, cat. no. DL101-01, China). Luciferase activity of the constructs was compared to the negative control (no amplified fragment cloned into the pGL3-RPL32 promoter plasmid) using a t-test, with P<0.05 indicating a significant difference.

RESULTS

Data summary

We compiled and integrated sequencing data from 66 epigenomes and 4745 transcriptomes across eight major tissues (cortex, cerebellum, hypothalamus, liver, lung, spleen, adipose tissue, and muscle), derived from three distinct embryonic layers. (Supplementary Table S1). The datasets included three key histone modifications (H3K4me3, H3K27ac, and H3K27me3), chromatin accessibility data from ATAC-seq for eight tissues, Hi-C data from two adipose tissues, and RNA-seq transcriptome data from at least 40 samples per tissue. However, the H3K27me3 data for two tissues (liver and spleen) exhibited a low signal-to-noise ratio, as determined through a review of their aligned IGV peak plots (Supplementary Figure S1). Consequently, only data from the remaining six tissues (cortex, cerebellum, hypothalamus, lung, adipose, and muscle) were used for subsequent analyses (Figure 1A).

Figure  1.  Identification of repressor elements (REs) using H3K27me3-OCR signal and gene expression data

A: Six sampled tissues derived from three embryonic layers. Cortex, cerebellum, and hypothalamus originate from the ectoderm, muscle and fat from the mesoderm, and lung from the endoderm. B: RE identification based on correlation of cross-tissue profiles of gene expression and H3K27me3-OCR activity. H3K27me3-OCRs were classified into three groups based on their correlation with target gene expression: positively correlated (Pos cor), negatively correlated (Neg cor), and uncorrelated (None cor). Negatively correlated H3K27me3-OCRs were defined as REs. C: Statistical analysis of identified REs in each tissue. UpSet plot shows shared and unique REs among different tissue groups, with colors corresponding to embryonic layers as indicated in Panel A.

DownLoad: Full-Size Img PowerPoint

RE identification using H3K27me3-OCR signal and gene expression data

While H3K27me3 is commonly associated with transcriptional repression, it is not a uniquely specific marker and often co-exists with activating histone modifications (Huang et al., 2019). Therefore, relying solely on H3K27me3 is insufficient for accurate identification of silencer elements. Following previously described methods (Huang et al., 2019), we initially focused on open chromatin regions (OCRs), which are accessible for activity region, and assessed whether they exhibited a significant negative correlation between the presence of H3K27me3 markers and transcription levels across different tissues (Wilcoxon rank-sum test, P<0.05; Figure 1B). This approach enables the identification of negatively correlated elements, including silencing promoters, silencing enhancers, and silencers, all categorized as REs. Therefore, to identify REs we first used peak calling results from H3K27me3 and ATAC-seq data sources from the above 6 tissues. To enhance reliability, peaks shared between two biological replicates were merged, yielding an average of 97 325 OCRs and 51 744 H3K27me3 regions per tissue. Subsequently, 54 920 non-redundant regions (H3K27me3-OCRs) were identified by isolating OCRs marked by H3K27me3 across the six tissues (Supplementary Table S3).

To identify REs, we assigned the two nearest neighboring genes or the host genes as target genes for each H3K27me3-OCR region. For each H3K27me3-OCR locus, the six tissues were divided into two groups based on the presence or absence of the H3K27me3-OCR signal in the relevant tissue (Figure 1B; existing/non-existing: ●/○; Supplementary Table S3). In total, 23 301 H3K27me3-OCR loci were identified as REs, with target gene expression levels significantly lower in tissues where the H3K27me3-OCR signal was present compared to those where it was absent (Wilcoxon rank-sum test; P<0.05; Supplementary Table S4). Among the tissues examined, the cortex exhibited the highest number of REs (8 393), while the hypothalamus had the fewest (4 097 REs) (Figure 1C). These findings suggest that REs are widespread across pig tissues.

Next, we investigated the tissue-specific characteristics of the REs. The cortex had the greatest number of tissue-specific REs (3 481), while the hypothalamus had the fewest (943) (Figure 1C). Interestingly, over half of the identified REs in the cerebellum were tissue-specific (50.4%) (Figure 1C). The number of REs per gene was calculated for each tissue, which indicated that the cerebellum had the highest RE/gene ratio (1.75:1), even when intergenic and overlapping REs were counted separately (Supplementary Table S5). This observation suggests that genes in the cerebellum may require more precise RE regulation. Although the underlying molecular mechanisms remain unclear, these findings suggest that the regulation of gene expression in the cerebellum differs from that in other tissues. Furthermore, we examined the sharing of REs among tissues derived from the same embryonic layer. Apart from the endoderm, where only one tissue was analyzed, the two mesoderm-derived tissues shared 263 REs, and the three ectoderm-derived tissues shared 283 REs (Figure 1C). This suggests that some REs may be acquired during embryonic layer differentiation, with a greater number of tissue-specific REs emerging during tissue formation.

Stronger inhibitory effects of SREs identified through H3K27me3 clusters on gene expression

Previous studies have demonstrated that H3K27me3 peaks, similar to H3K27ac, are enriched in specific regions (Cai et al., 2021; Pang et al., 2023). Furthermore, knockout of H3K27me3-rich regions has been shown to alter cell identity and influence the growth of xenografted tumors (Cai et al., 2021). Building on this understanding, we utilized a strategy analogous to super-enhancer identification to identify H3K27me3-rich regions (Figure 2A) (Hnisz et al., 2013). In total, 4 477 H3K27me3-rich regions were identified across six tissues, characterized by clusters of H3K27me3 peaks within 4 kb regions of genomic DNA (Methods; Supplementary Table S6). Interestingly, more than half of these H3K27me3-rich regions (2 608) did not contain any above-identified REs, potentially due to the exclusion of chromatin accessibility and gene expression data in the identification process. The remaining H3K27me3-rich regions (1 869) contained an average of 1.62 REs, suggesting that some REs may function synergistically within these regions (Supplementary Table S7).

Figure  2.  Identification and characterization of super repressor elements (SREs)

A: Identification of H3K27me3-rich regions, similar to the identification of super-enhancer elements. B: Statistical analysis of identified SREs in each tissue. UpSet Plot shows shared and unique SREs among different tissue groups. C: Comparison of expression levels among target genes of SREs, typical repressor elements (TREs), and positively correlated H3K27me3-OCRs (Pos cor) (using adipose tissue as an example). SREs showed stronger inhibitory effects on gene expression. P-values were calculated using Wilcoxon rank-sum test. D, E: Distribution (D) and distance from transcription start site (TSS) (E) of SREs and TREs in pig genome.

DownLoad: Full-Size Img PowerPoint

Subsequently, the identified REs were categorized into SREs and TREs based on their overlap with H3K27me3-rich regions, resulting in the identification of 2 034 non-redundant SREs and 22 223 non-redundant TREs (Supplementary Table S8). We examined the expression of target genes with positive correlation with H3K27me3-OCR, SREs or TREs in each tissue. The expression levels of target genes showed lowest in SREs, followed by TREs, and highest in positive correlation H3K27me3-OCRs (Figure 2C, Supplementary Figure S2). These findings support the reliability of our SRE results.

Recent perspectives have suggested that silencer elements can exert their repressive effects on promoters either at short-range or long-range (Pang and Snyder, 2020; Zhang et al., 2022). To explore this further, we determined the distance between the transcriptional start site (TSS) of associated genes and SREs or TREs. Results showed that approximately 50% of both SREs and TREs were located in distal intergenic regions, more than 3 kb from the TSS (Figure 2D, E). Notably, SREs were more frequently located within gene promoter regions compared to TREs (Figure 2D). Furthermore, REs were distributed across a wide range of distances from the associated TSSs (0–1 kb, 1–3 kb, 3–5 kb, 5–10 kb, 10–100 kb, and >100 kb), with most being situated within the 10–100 kb range (Figure 2E).

We further analyzed the regulatory relationship between gene expression activity and chromatin interactions associated with REs (SREs or TREs) at the loop level using Hi-C data. Focusing on adipose tissue as an example (Jin et al., 2023), 62.8% of SREs and 61.3% of TREs were located within the same chromatin loop as their target genes (Figure 3A). Interestingly, genes within loops containing REs in adipose tissue did not exhibit significantly lower expression levels compared to those with loops lacking REs (non-negatively correlated OCRs) (Figure 3B). These findings suggest that only a limited number of genes within a loop were regulated by REs. Moreover, genes within loops containing TREs had lower expression levels compared to those within loops containing SREs (Figure 3B). Based on the expression levels of genes adjacent to SREs and TREs (Figure 2C), the distances between REs (SREs or TREs) and the TSSs of their associated genes (Figure 2E), the expression levels of genes within loops containing SREs or TREs (Figure 3B), and the minimal effect of the small number of extreme values in the distribution, we infer that SREs exert strong and tight regulation over a small number of target genes within a specific genomic region, while TREs exert broader but weaker regulatory influence over a wider range of target genes (Figure 3C). This aligns with the observation that H3K27me3-rich regions (MRRs) are characterized by dense chromatin interactions connecting target genes (Cai et al., 2021).

Figure  3.  Super repressor elements (SREs) and typical repressor elements (TREs) exhibit different regulatory patterns

A: Statistical analysis of whether repressor elements (REs) and their target genes were located within the same loop (using adipose tissue as an example). B: Comparison of gene expression levels among genes located in the same loop with SREs, TREs, or non-negatively correlated OCRs (No-neg Cor) (using adipose tissue as an example). P-values were calculated using Wilcoxon rank-sum test. C: Differential regulatory mechanisms of SREs and TREs. SREs exhibited a “tight” model, while TREs exhibited a “loose” model. In a loop, the impact of TREs on gene expression reduction was weaker but more widespread. D: Comparison of proportion of transcription factors (TFs) in SRE target genes and TRE target genes. P-values were calculated using the Wilcoxon rank-sum test. E: Enrichment analysis for motifs of repressor TFs in TREs across tissues. Colors in heatmap represent value of −log₁₀(P-value). ns: Not significant. F: Visualization of normalized signals for ATAC-seq, H3K27me3, and H3K27ac within 100 kb upstream and downstream of SRE, TRE, and positively correlated H3K27me3-OCR (Pos cor) centers (using adipose tissue as an example).

DownLoad: Full-Size Img PowerPoint

Previous studies have indicated that transcription factor genes require exceptionally precise regulatory mechanisms (Niwa et al., 2000). Thus, we examined whether the SRE and TRE target genes were transcription factors (TFs). Surprisingly, approximately 25% of SRE target genes were found to encode TFs, a significantly higher proportion than observed for TREs (Figure 3D). This finding is consistent with the proposed model of tight regulation by SREs versus looser regulation by TREs (Figure 3C).

To further investigate the potential involvement of TFs in the formation and function of REs, we performed motif enrichment analyses across all six tissues, using random genomic regions as controls, similar to the approach of previous study (Jayavelu et al., 2020). Results revealed a consistent and significant enrichment of motifs for TFs with known repressive activity (ATF3, ZFX, REST, FOSTL2, and RFX family TFs (RFX2, RFX3, and RFX5)) across all TREs in the six tissues studied (Figure 3E; Supplementary Table S9). Furthermore, the ZFX motif was significantly enriched in both SREs and TREs across all tissues (Supplementary Figure S3 and Table S9). Motifs for other repressive TFs, such as PTF1A, BACH2, ETS1, and the TEAD family (TEAD2 and TEAD4), were also significantly enriched in one or more tissues (Figure 3E; Supplementary Figure S3 and Table S9). These findings provide evidence that the identified REs may exert their repressive functions through the regulation of TFs.

Finally, the chromatin landscape surrounding SREs and TREs was analyzed by examining ATAC-seq, H3K27me3, and H3K27ac signals within 100 kb upstream and downstream of the centers of these elements (Figure 3F; Supplementary Figure S4). As expected, SREs exhibited stronger H3K27me3 signals compared to TREs and positively correlated H3K27me3-OCRs. Notably, both SREs and TREs, as well as positively correlated H3K27me3-OCRs from all tissues, displayed high levels of H3K27ac modification at their centers, consistent with previous studies (Cai et al., 2021; Jayavelu et al., 2020). These results suggest that some of the identified REs may act as different functional REs (silencing promoters, silencing enhancers, and silencers) or activator elements (enhancers and promoters) in different cell types.

Potential role of REs in tissue differentiation via inhibition of genes with specific biological functions

The binding of repressor element 1-silencing transcription factor (REST) to silencers, which recruits repressive cofactors and chromatin modifiers to suppress neuronal gene transcription in non-neuronal cells (Chong et al., 1995; Schoenherr and Anderson, 1995), served as a basis for examining the functions of target genes regulated by REs.

Analysis of genes repressed by both SREs and TREs across various tissues revealed a notable enrichment in terms related to “cell fate commitment”, “organ development”, and “embryonic morphogenesis” (Figure 4A; Supplementary Figure S5A and Table S10), suggesting involvement in organ development and cell differentiation. Surprisingly, even in neuronal tissues such as the cerebellum and hypothalamus, there was a pronounced enrichment in functions associated with nervous system differentiation, such as “generation of neurons”, “neuron differentiation”, and “neurogenesis” (Figure 4A; Supplementary Figure S5A and Table S10).

Figure  4.  Biological functions of repressor elements (REs) in tissue development and cell lineage differentiation

A: Bubble plot displaying top results of Gene Ontology (GO) functional enrichment analysis of super repressor element (SRE) target genes in various tissues. No significant GO enrichment was found in the cortex. Bottom numbers display number of other significant GO enrichment terms, which are not shown. B: Gene set enrichment analysis (GSEA) of REs from neural tissues (cerebellum, cortex, and hypothalamus). Gene sets used in GSEA were derived from mouse orthologs projected from target genes of pig REs. Ranked gene lists were sorted based on log₂(FC) values from differential expression analyses of mouse basal progenitor over neural stem cells. GSEA of “newborn neuron vs basal progenitor” and “newborn neuron vs neural stem cell” are shown in Supplementary Figure S5. Normalized enrichment score (NES) and false discovery rate (FDR) were calculated using the “GSEA” module in ClusterProfiler. C: GSEA of target genes of REST-associated REs. Ranked gene lists were sorted based on log₂(FC) values from differential expression analyses at different developmental stages in a murine neuronal lineage. D: GO analysis using merged target genes of SREs in nervous system tissues (cortex, cerebellum, and hypothalamus). Gene-concept network displaying major GO terms and their core genes. E: Visualization of histone modifications (H3K4me3, H3K27ac, and H3K27me3) in REs of MFSD2A in six tissues.

DownLoad: Full-Size Img PowerPoint

These findings suggest that the developmental stage-specific repression of genes may play a crucial role in tissue development and cell lineage differentiation. To further explore the functions of the identified REs, transcriptomic data from various developmental stages of a murine neuronal lineage (neural stem cells (NSCs), basal progenitors (BPs), and newborn neurons (NBNs)) (Mukhtar et al., 2020) were analyzed using GESA. Notably, RE target genes in the cerebellum, cortex, and hypothalamus began to be repressed during the transition from NSCs to BPs (Figure 4B). In contrast, during the transition from BPs to NBNs, a small subset of these RE target genes were reactivated (Supplementary Figure S5B, C). This pattern suggests that RE-mediated repression is initiated during the differentiation of stem cells into progenitors and is sustained through maturation into fully differentiated cells. These results underscore the importance of gene silencing in tissue development, similar to gene activation.

Many studies have noted the critical role of REST in non-neuronal cells (Chong et al., 1995; Schoenherr and Anderson, 1995). Interestingly, REST motifs were widely enriched within REs across all six tissues examined (Figure 3E; Supplementary Figure S3). Thus, we further investigated the biological functions of REST-enriched RE-associated genes in the context of cell lineage differentiation. GESA revealed that REST-associated genes were significantly up-regulated during the differentiation of neural progenitor cells into neurons (Figure 4C, Supplementary Table S9). Integrating this with the known functions of REST from earlier studies (Chong et al., 1995; Schoenherr and Anderson, 1995), it is evident that REST-associated genes play pivotal roles in neuronal tissue development. Conversely, the repression of these genes is essential for the proper development of non-neuronal tissues.

To further clarify the role of REs in nervous system development, GO enrichment analysis was conducted on the combined target genes of REs in nervous system tissues. Results showed that the silenced genes were predominantly associated with non-nervous system functions, including “renal system development”, “regulation of ossification”, and “embryonic morphogenesis” (Figure 4D). Additionally, four key regulatory genes (OSR2, SIX2, SFRP1, and TBX18) were consistently implicated in these enriched pathways (Figure 4D).

We next analyzed genes activated in the nervous system and repressed in the non-nervous system. Interestingly, the major facilitator superfamily domain-containing protein 2 lysolipid transporter A (MFSD2A) gene, known as a major transporter for docosahexaenoic acid (DHA) imported across the blood-brain barrier (BBB) into the brain (Nguyen et al., 2014), was identified. The promoter region of MFSD2A exhibited tissue-specific modifications, with strong H3K27ac and H3K4me3 modifications observed in the cortex, cerebellum, and hypothalamus, but weaker H3K4me3 and stronger H3K27me3 modifications observed in adipose, muscle, and lung tissue (Figure 4E). This pattern suggests that this region acts as an activator (enhancer and promoter) in nervous system tissues but as a repressor in non-nervous system tissues (Figure 4E).

Subsequently, we examined the functional enrichment of genes associated with tissue-specific SREs and TREs (Figure 5A, B). Results showed a significant enrichment of ion transmembrane transport-related functions for genes associated with SREs and TREs in adipose tissue (Figure 5A, B). This suggests that the silencing of genes related to ion transmembrane transport plays a crucial role in maintaining adipocyte function (Figure 5A, B). Taken together, these results suggest that the repression of specific functional genes is critical at various developmental stages and in specific tissues.

Figure  5.  Functional analysis of tissue-specific repressor elements (REs)

A, B: Bubble plot displaying top results of Gene Ontology (GO) functional enrichment analysis for target genes of tissue-specific super repressor elements (SREs) (A) and typical repressor elements (TREs) (B) in various tissues. Bottom numbers are number of other significant GO enrichment terms, which are not shown. Transport GO terms enriched in adipose tissues are highlighted by the red boxes. SRE target genes from hypothalamus and cortex, and TRE target genes from cerebellum and cortex did not show significant GO enrichment.

DownLoad: Full-Size Img PowerPoint

Conversion of activator elements and REs regulates phenotypic traits in pigs

Recent research has demonstrated that certain genomic regions can function as either activators or repressors, depending on the cellular context, associated with histone modifications such as H3K27ac and H3K27me3 (Bandara et al., 2021; Gisselbrecht et al., 2020; Ogbourne and Antalis, 1998). Identifying these bifunctional elements across different tissues is crucial for unraveling the complex molecular mechanisms by which QTLs affect traits. In this study, 5 794 bifunctional elements were identified across six tissues (Supplementary Table S11), 381 of which overlapped with fine QTLs (Supplementary Table S11) (Hu et al., 2022). These QTLs were associated with economically important traits such as “Average daily gain”, “Average backfat thickness”, and “Intramuscular fat content” (Supplementary Table S11).

Two bifunctional elements (chr1:248 588 115–248 590 169, chr1:248 606 955–248 611 087, Sscrofa11.1) were identified as activators in adipose tissue but as repressors in the cortex, lung, hypothalamus, cerebellum, and muscle. These elements are located near the induced pluripotent stem cell (iPSCs)-related TF Krüppel-like transcription factor 4 (KLF4), which is within the pig “drip loss” QTL region (Figure 6A). Analysis of chromatin signatures for active and repressive histone modifications (H3K4me3, H3K27ac, and H3K27me3) (Figure 6A) revealed that the histone modification patterns of these two bifunctional elements across different tissues did not entirely match the expression levels of KLF4 (Figure 6A, B). This discrepancy suggests the presence of additional regulatory elements influencing KLF4 expression. Further investigation identified 32 additional REs surrounding KLF4 in the six tissues (Figure 6A; Supplementary Table S2). Regression analysis demonstrated a negative correlation between the number of REs and the expression level of KLF4 (Figure 6B), suggesting that KLF4 expression may be regulated through the additive effects of multiple regulatory elements.

Figure  6.  Activation-repression transition of KLF4 regulatory elements between different tissues regulates economically important traits

A: Visualization of activation and repression of histone modifications, repressor elements (REs), and economic trait QTLs around KLF4. B: Negative correlation between expression level of KLF4 (ln-transformed TPM) and number of identified REs in six tissues. Regression equation, R², and P-value were calculated using the “Line_Regression” module in Hiplot.

DownLoad: Full-Size Img PowerPoint

To validate the potential dual functions of the identified repressor regions, experiments were conducted using three cell lines (PIEC, ST, and PK15 cells). We successfully cloned 20 of the 32 identified REs (Supplementary Table S2). In PIEC cells, the majority of these elements demonstrated enhancer function, whereas in ST and PK15 cells, only a small fraction exhibited such activity. Comparative between different cell lines revealed that most of these elements functioned as silencing enhancers in ST and PK15 cells (Figure 7). This finding suggests that some of these elements possess dual roles as activators and repressors, with the potential to switch between these functions depending on the tissue type. This functional flexibility could influence the regulation of economically important phenotypes in pigs. Consequently, relying on H3K27me3, ATAC-seq, and RNA-seq data is insufficient for the accurate identification of silencers, as it does not adequately exclude between silencing enhancers and promoters.

Figure  7.  Dual-luciferase reporter analysis for identified repressor elements around KLF4 gene

Experiments were performed in PIEC, ST, and PK15 cells. Relative luciferase activity was defined as the normalized ratio (ratio of Firefly luciferase and Renilla luciferase (fLuc/rLuc)) by activity of the promoter (P). P-values were calculated using a t-test (^***: P≤0.001; ^**: P≤0.05; ^*: P≤0.1). Elements significantly different to the promoter are highlighted in blue.

DownLoad: Full-Size Img PowerPoint

DISCUSSION

Over the past decade, genome-wide association studies (GWAS) have been applied to investigate the associations between genetic variations and traits of interest across species (Wei et al., 2023; Yang et al., 2021; Yengo et al., 2022). In human disease research, GWAS have revealed that over 88% of phenotype-associated variations are located outside protein-coding regions, predominantly within gene regulatory regions (Hindorff et al., 2009).

To gain a better understanding of the molecular consequences of mutations within repressive regions and their association with phenotypic traits, we systematically identified 23 301 REs across six different tissues in pigs. Analysis of gene expression levels surrounding the identified REs revealed significantly reduced expression, confirming the reliability of our findings. While our methodology is effective, it is not without limitations, including the influence of tissue-specific factors and specific histone modifications analyzed. The development of novel linear models for the systematic identification of REs is therefore essential, although devising such models necessitates a comprehensive understanding of the diverse factors that regulate gene expression. At present, distinguishing between silencing promoters, silencing enhancers, and silencers remains challenging when relying solely on ATAC-seq, H3K27me3, and RNA-seq data. Current studies often assume that all identified REs function as silencers (Cai et al., 2021; Huang et al., 2019; Ouyang et al., 2023). However, we argue that new analytical approaches need to be developed to accurately differentiate between these various types of regulatory elements.

To exclude silencing promoters and enhancers from REs, it is essential to refine the identification pipeline based on their fundamental property of reducing gene expression. The first step involves identifying genomic regions that regulate expression levels using an expression quantitative trait loci (eQTL) approach (Zhong et al., 2023). Next, the activity of each genotype within these regions should be evaluated using a deep learning model trained with STARR-seq data. Silencers can then be identified based on the combination of their activity values and the expression levels associated with each genotype (Figure 8). Given the significant biological role of silencers, developing a precise, genome-wide method for their identification using this strategy is crucial.

Figure  8.  Schematic for an approach to identify silencers

Expression quantitative trait loci (eQTL) analysis was performed in pig tissues and cells. Based on STARR-seq data, a deep learning model was constructed to predict the sequence activity for each allele from eQTL analysis. Silencers were identified by comparing sequence activity and expression level between alleles.

DownLoad: Full-Size Img PowerPoint

One of the earliest identified silencers is the HMRE element in brewer’s yeast, responsible for suppressing gene expression at the HML and HMR loci (Brand et al., 1985), crucial for mating type determination (Brand et al., 1985). In fruit flies, the VRE silencer mediates ventral gene suppression in response to dorsal morphogen gradients, initiating the differentiation of mesoderm and neuroectoderm during embryogenesis (Jiang et al., 1993). Silencers associated with the CD4 gene influence T cell lineage development (Sawada et al., 1994). Although research on inhibitory elements is less extensive compared to studies on activators, their role in organ development is likely just as important. In the present study, we observed widespread inhibition of neuro-related genes across various tissues, with RE target genes in neuronal tissues suppressed during the differentiation from stem cells to progenitors, a state maintained through full maturation. This is an intriguing finding, and further exploration of inhibitory factors in organ development is warranted. Elucidating the molecular mechanisms underlying this phenomenon is essential for a comprehensive understanding of developmental and differentiation processes.

The use of dual-functional sites that can switch between activator and repressor roles appears to be a common mechanism by which REs exert their functions (Bandara et al., 2021; Huang and Ovcharenko, 2022). The same sequence may function as a RE in certain tissues and an activator in others, thereby influencing gene regulation in a tissue-specific manner. This phenomenon may offer a partial explanation for the question, “Why we are so different?”. Consequently, future epigenetic studies should aim to simultaneously measure both promotive histone modifications (such as H3K27ac and H3K9ac) and inhibitory histone modifications (such as H3K27me3 and H3K9me3) across multiple tissues, even in the single cell level (Xiao et al., 2024), to further dissect the complexities of these dual-functional elements. In earlier research, we identified imprinted genes in pigs that exhibit parent-of-origin-specific differential expression (Wu et al., 2020). Whether these genes are regulated by REs remains an open question that urgently needs to be answered. Moreover, while we have proposed a model to explain the impact of REs on KLF4 expression, further verification through detailed animal experiments is required to substantiate these findings.

CONCLUSIONS

In this study, we identified 2 034 SREs and 22 223 TREs within the pig genome across six tissues. SREs exhibited tight regulation of target genes, while TREs exerted a broader but loose influence. Through the integration of transcriptome data from neural cells at different differentiation stages, we found that RE-mediated repression of target genes began during the transition from neural stem cells to basal progenitors and persisted throughout their maturation into fully differentiated cells. Furthermore, our findings showed that the silencing of genes involved in ion transmembrane transport may play an important role in maintaining the function of adipose tissues, with specific genes inhibited in the tissue differentiation process. We also observed that multiple regulatory elements worked synergistically, with additive effects, to regulate the expression of KLF4, thereby influencing economically important phenotypes in pigs. This RE map of the pig genome provides valuable insights for future research on REs and holds significant potential for advancing pig breeding strategies.