Comparative Analysis of Promoters and Enhancers in the Pituitary Glands of the Bama Xiang and Large White Pigs

The epigenetic regulation of gene expression is implicated in complex diseases in humans and various phenotypes in other species. There has been little exploration of regulatory elements in the pig. Here, we performed chromatin immunoprecipitation coupled with high-throughput sequencing (ChIP-Seq) to profile histone H3 lysine 4 trimethylation (H3K4me3) and histone H3 lysine 27 acetylation (H3K27ac) in the pituitary gland of adult Bama Xiang and Large White pigs, which have divergent evolutionary histories and large phenotypic differences. We identified a total of 65,044 non-redundant regulatory regions, including 23,680 H3K4me3 peaks and 61,791 H3K27ac peaks (12,318 proximal and 49,473 distal), augmenting the catalog of pituitary regulatory elements in pigs. We found 793 H3K4me3 and 3,602 H3K27ac peaks that show differential activity between the two breeds, overlapping with genes involved in the Notch signaling pathway, response to growth hormone (GH), thyroid hormone signaling pathway, and immune system, and enriched for binding motifs of transcription factors (TFs), including JunB, ATF3, FRA1, and BATF. We further identified 2,025 non-redundant super enhancers from H3K27ac ChIP-seq data, among which 302 were shared in all samples of cover genes enriched for biological processes related to pituitary function. This study generated a valuable dataset of H3K4me3 and H3K27ac regions in porcine pituitary glands and revealed H3K4me3 and H3K27ac peaks with differential activity between Bama Xiang and Large White pigs.

Introduction

Cis-regulatory elements play fundamental roles in the regulation of gene expression. Recent advances in ChIP-Seq, e.g., targeting H3 lysine 27 acetylation (H3K27ac) and H3 lysine 4 trimethylation (H3K4me3) (Barski et al., 2007), have enabled a dramatic increase in the annotations of regulatory regions, including enhancers and promoters in organisms including human, mice, and Drosophila melanogaster (Gerstein et al., 2010; Ernst et al., 2011; Shen et al., 2012). This has led to extensive data resources for understanding the functions of regulatory elements and molecular mechanisms underlying regulatory variants associated with diseases (Corradin and Scacheri, 2014; Zhao et al., 2018). Meanwhile, a number of studies have annotated the regulatory elements in farm animals, such as pigs (Zhao et al., 2021), cows, sheep, and chickens, coordinated by the Functional Annotation of Animal Genomes (FAANG) consortium (Andersson et al., 2015; Kern et al., 2021).

The pig is an important livestock species, serving as both a major food source and an excellent animal model for biomedical research in humans (Vodicka et al., 2005; Meurens et al., 2012). European and Chinese pigs have a divergent evolutionary history of ∼0.6 million years and were domesticated in Europe and China independently ∼9,000 years ago (Larson et al., 2005; Frantz et al., 2013). The Bama Xiang pig is a typical indigenous breed from south China, characterized by its small body size (∼60 kg at 300 days old) (Gong et al., 2019) and two-end black coat color. The Large White pig is a European breed characterized by its large and long body size (∼100 kg at 168 days old) with excellent hams, a white coat color, and a high level of reproduction; it serves as an important maternal breed for the international pork production industry (Rubin et al., 2012). The large phenotypic differences between the Bama Xiang and the Large White make them valuable animal models for comparative genetic and epigenetic studies. The pituitary gland is a master endocrine organ with an important function in regulating the growth, metabolic status, immunoregulation, reproduction, nervous system, and lactation of animals (Brinkmeier et al., 2009; Kiezun et al., 2014; Shan et al., 2014). Thus, it may play an important role in the formation of phenotypic differences between the two breeds. A previous study compared the pituitary gland transcriptome between the Bama Xiang and the Large White (Rubin et al., 2012). However, the landscapes of promoter and enhancer activity in the pituitary glands of the two breeds have not been reported.

This study profiles and characterizes the distribution of H3K4me3 and H3K27ac in the pituitary glands of two Bama Xiang and two Large White adult pigs and identifies the differential activity of H3K4me3 and H3K27ac peaks between the two breeds.

Materials and Methods

Samples

We carefully dissected the pituitary glands of two biological replicates (full siblings) of 150-day-old LW (one male, one female) and 132-day-old female BMX pigs, respectively, following the standardized sample collection protocols of the FAANG Project¹ (Andersson et al., 2015). The ChIP sequencing data related to these samples were deposited at the Gene Expression Omnibus with the accession code GSE178380².

Chromatin Immunoprecipitation Sequencing

Chromatin immunoprecipitation was performed using a SimpleChIP^® Plus Enzymatic Chromatin IP Kit (Magnetic Beads, 9005) with 500 μg chromatin and 3 μg antibody H3K4me3 (active motif, 39,915) and 5 μg antibody H3K27ac (active motif, 39,133). The detailed protocols were obtained from https://www.encodeproject.org/about/experiment-guidelines/ and https://www.animalgenome.org/community/FAANG. The dissected tissues were treated with 37% formaldehyde to cross-link the proteins covalently with DNA. This was followed by cell disruption and sonication to shear the chromatin to a target size of 100–300 bp (Landt et al., 2012; Shen et al., 2012). Then, the modified histone with its bound DNA was enriched using the corresponding antibody. After protein removal, the DNA was purified, and real-time quantitative polymerase chain reaction (qPCR) was performed. ChIP and input (control sample) library construction and sequencing procedures were carried out according to the Illumina protocols, with minor modifications (Illumina, San Diego, CA, United States).

Trimmed clean reads were mapped to the pig reference genome Sscrofa 11.1³ using the Burrows-Wheeler Aligner (Abuin et al., 2015), allowing two mismatches and removing duplicated ChIP-seq reads using SAMtools. After that, Model-based Analysis for ChIP-Seq (MACS) version 2.1.0 peak caller was applied to infer the histone modification regions, i.e., the peaks (Zhang et al., 2008) with the following parameters: -t Sample_ac_R1_sorted.bam -c Sample_input_R1_sorted.bam –broad -g 2.48e9 –broad-cutoff 0.1 -n Sample_ac. The peak regions identified in individual sample were integrated by intersecting all peaks across datasets, with less than 1 kb distance between summits, using bedtools version 2.27.0 (Quinlan and Hall, 2010). We calculated the Spearman correlation coefficients of the BAM files of various ChIP-seq samples using the deepTools bamCorrelate module (Ramirez et al., 2016), and TSS enrichment analysis was performed using deepTools bamCoverge module. The library complexity was calculated using bedtools bamtobed module.

RNA Sequencing

We downloaded RNA sequencing data (only BAM files are available) of the pituitary glands of three Bama Xiang pigs and three Large White pigs from gsa.big.ac.cn under the GSA number CRA000876⁴ (Yang et al., 2018). We extracted fastq files from the BAM files using bedtools (bedtools bamtofastq -i sample.bam -fq sample.fq) and remapped the fastq files to Sus scrofa 11.1 using STAR-2.5.3a (Dobin et al., 2013). The transcripts were assembled and merged with Stringtie with reference to Ensembl GTF (98.111), which were further merged with Ensembl GTF (98.111) to obtain a customized GTF file as the reference for the gene expression quantification. The expression levels of the genes were quantified using FeatureCounts 1.5.3 (Liao et al., 2014). After removing the mitochondrial genes, the RPKM algorithm from DESeq2 (Love et al., 2014) was used to normalize the expression values for each sample.

Annotations of Putative Promoters and Enhancers Across Pig Genomic Features

The GAT program⁵ was used to investigate H3K4me3 and H3K27ac enrichment across genomic features with the following parameters: gat-run.py –segment-file = segments.bed.gz –workspace-file = workspace.bed.gz –annotation-file = annotations.bed.gz. The bedtools output was visualized with Integrative Genomics Viewer version 2.4, a genomic dataset viewer that allows for the visualization of genomic features. We downloaded whole genome GERP score from Ensembl Database and computed the GERP score of a peak by averaging the GERP scores of bases within that peak. GC-content of a peak was calculated using EMBOSS geecee software (Rice et al., 2000).

Quantification of the H3K4me3 and H3K27ac Peaks and Identification of Peaks With Differential Activity Between the Two Breeds

The H3K4me3 and H3K27ac peaks were merged using the merge command in bedtools (Quinlan and Hall, 2010). The read depths in a peak region were calculated using the SAMtools bedcov utility (version 1.2) and then normalized to the RPKM value by DESeq2 R package (Love et al., 2014), where the normalized value is defined as peak activity. We further used the DESeq2 program to identify peaks to show the differential activity between the two breeds. Regions with P of less than 0.05 and fold changes greater than 2 were assigned as peaks with differential activity.

Annotation and Gene Ontology Analysis of the Gene Set

Gene functions were annotated using DAVID⁶ and ClueGO (Bindea et al., 2009), and the background gene sets for enrichment analyses were all genes from the human genome. The P values for the enrichment of GO terms were corrected using the Benjamini-Hochberg approach. HOMER⁷ was adopted to investigate the enrichment of binding motif of transcription factors (TFs) in the corresponding set of peaks.

Identification of Super Enhancer

The ROSE algorithm (Whyte et al., 2013) was used to identify super enhancers in all H3K27ac samples. The definition of super enhancer activity was the same as that for all peaks of H3K27ac and H3K4me3.

Results and Discussion

Mapping Promoters and Enhancers in the Pituitary Glands of BMX and LW Pigs

We profiled H3K4me3 and H3K27ac in the pituitary glands of two Bama Xiang (132 days old) and two Large White pigs (150 days old) using chromatin immune-precipitation, followed by high-throughput sequencing (ChIP-Seq) of single end 50 bp reads (Supplementary Figure 1A and Table 1). The average number of reads uniquely mapped to the reference genome was 25.5 M, with a range from 22.3 M to 26.8 M (Table 1). Clustering analyses performed on the BAM files based on the read depths of a set of 65,044 merged peaks first grouped the samples into two clusters represented by the two histone markers; within both clusters, the samples from the same breed were grouped together (Supplementary Figure 1B), suggesting that between-breed differences in the profiles of both H3K4me3 and H3K27ac were greater than those within breeds. PCA analyses yielded similar results (Supplementary Figure 2). We calculated the library complexity of the ChIP-seq data and obtained average Non-Redundant Fraction (NRF), PCR Bottlenecking Coefficient 1 (PBC1), and PCR Bottlenecking Coefficient 2 (PBC2) values of 0.97, 0.98, and 44.61, respectively, for the investigated samples, suggesting that the library constructions were reasonably good (Supplementary Table 1).

TABLE 1

Table 1. Summaries of ChIP-Seq experiments.

On average, we found 19,296 H3K4me3 and 43,733 H3K27ac peaks per sample (Table 1). The average/median peak lengths were 3,421/2,429 bp and 3,835/1,792 bp for the H3K4me3 and H3K27ac markers, respectively (Table 1), comparable to a previous report (Villar et al., 2015) (Supplementary Figure 3). Next, we merged all H3K4me3 and H3K27ac peaks in each individual, generating 23,680 H3K4me3 (21,396 in the BMX and the 19,553 in LW) and 61,791 H3K27ac (55,585 in the BMX and 48,333 in the LW) non-redundant peaks (Supplementary Figure 1A). The 61,791 H3K27ac peaks were further grouped into 12,138 proximal peaks that intersected with a ± 1 kb region from the transcription start site (TSS) of any transcript provided by the Ensembl database (release 98) and the remaining 49,473 distal peaks. These regions constituted a valuable resource of regulatory elements in the pituitary of the pigs. Next, we merged all H3K4me3 and H3K27ac peaks and identified 16,099 regions that overlapped with both histone modification markers, then observed significant positive correlations between the quantitative activities of the two histone modification markers across the 16,099 regions in all four individuals (r = 0.28∼0.47, P < 2.2 × 10^–16) (Figures 1A,B), suggesting that H3K4me3 and H3K27ac coordinately regulated gene expression. In addition, the peak region of H3K4me3 modifications had a higher GC content and extent of sequence conservation across species (GERP score) compared to the H3K27ac peak regions (Figures 1C,D).

FIGURE 1

Figure 1. The intersection and conservation of H3K27ac and H3K4me3 peaks. (A) The intersection of H3K27ac and H3K4me3 peaks. (B) The correlation of 16,099 overlapped peaks between two histone modifications in four individuals. (C,D) Boxplots comparing the GC contents and GERP scores of H3K27ac and H3K4me3 peaks.

We assumed that the genes close to the peak with high activity may play an important function in the pituitary gland. At the 16,099 regions, we identified 100 peaks ranked within the top 5% of both histone markers in all four individuals based on their activity (see section ‘‘Materials and Methods’’). The 100 peaks overlapped with 114 genes involved in the intracellular estrogen receptor signaling pathway, positive regulation of blood pressure, and neuroepithelial cell differentiation. Among these genes, BCL6, CDKN1B, ESRP1, PITX2, and SIX6 genes had functions related to the pituitary based on the annotation information from the DAVID bioinformatics resources⁸. For example, BCL6 is involved in the reverse coordination of endogenous impulses released by pituitary growth hormone (GH) by interacting with STAT5 (Meyer et al., 2009), and CDKN1B is associated with adreno-cortico-tropic hormone hypersecretion and the formation of pituitary neoplasms (Tichomirowa et al., 2012).

Genomic Annotations of Putative Promoters and Enhancers

Next, we examined the overlap of genomic features with the H3K4me3 and H3K27ac peaks. The H3K4me3 and H3K27ac proximal peaks were mostly located in the intron (61.36 and 66.96%) and intergenic (24.61 and 20.57%) regions, while the H3K27ac distal peaks tended to be located in intron (71.88%) and intergenic regions (23.81%) (Figure 2A). The H3K4me3 (14.03%) and H3K27ac proximal (12.47%) peaks had stronger enrichment in the transcription initiation region than the H3K27ac distal peaks (4.31%), agreeing with the fact that the promoters tend to be located near the TSS of genes, while enhancers are usually distal to the target genes (Heintzman and Ren, 2009; Pennacchio et al., 2013) (Figure 2A). This suggests that intronic sequences close to the TSS regions could play important regulatory roles. In addition, both histone modifications were enriched at 3 kb near the gene TSS, and the H3K4me3 ChIP-Seq reads showed a stronger enrichment at the TSS region than H3K27ac ChIP-Seq reads (Figure 2B). The expression levels of the genes showed significant positive correlations with the H3K4me3 and H3K27ac activity around their TSS regions (e.g., ± 5 kb) across the genome, although ChIP-Seq and RNA-Seq data (Yang et al., 2018) were generated from different samples. The correlations with gene expression were stronger for H3K4me3 than for H3K27ac. Moreover, the strengths of the correlations of the gene expressions with the ChIP-Seq peak activities decayed when greater regions centered on the TSS were examined (Supplementary Figure 4).

FIGURE 2

Figure 2. Enrichment of H3K4me3 and H3K27ac peaks across genomic features and transcription start sites. (A) Distribution of H3K4me3 and H3K27ac peaks across genome. (B) Heatmaps depicting normalized ChIP-seq signal (H3K27ac and H3K4me3) at 3 kb near the TSS, sorted by signal intensity.

Diversity of Promoter and Enhancer Peaks Among Individuals and Breeds

Next, we quantified the peak activity in each sample based on read depths in the corresponding peak region normalized by peak length, using the R package DESeq2 (see section “Materials and Methods”). We evaluated the correlation among samples in terms of genome-wide H3K4me3 and H3K27ac peak activity. The mean correlations based on H3K4me3 peak activity (0.81) were higher than those based on the H3K27ac activity (0.42), suggesting that the activity of the promoters was more conservative than those of enhancers.

We then sought to identify H3K4me3 and H3K27ac peaks that showed differential activity between the two breeds. At the threshold of fold change ≥ 2 & P ≤ 0.05, we identified 793 H3K4me3 peaks (509 and 284 showed higher activity in BMX and LW, respectively) (Figure 3A), and 3,602 H3K27ac peaks showed differential activity between the two breeds (1,614 and 1,918 showed higher activity in BMX and LW, respectively) (Figure 3B). We displayed the activity of the 20 H3K4me3 and 20 H3K27ac peaks that show the most significant differential activity among two breeds (Figures 3C,D).

FIGURE 3

Figure 3. Differential activity analysis on the H3K4me3 and H3K27ac peaks between two breeds and representative validation of differential peaks exhibiting distinct activity of the H3K27ac signal. (A,B) Volcano plot showing the differential activity analysis of H3K4me3 and H3K27ac peaks between BMX and LW breeds. (C,D) Unsupervised hierarchical clustering of the top differential H3K4me3 and H3K27ac peaks between BMX and LW samples. (E,F) The activity track of H3K27ac regions of (Chr13: 75916467–75920943) and (Chr2: 25338342–25347879) with increased hyper-acetylation in the BMX pituitary was related to genes of EPHB1 and FJX1 respectively. (G,H) The activity track of H3K27ac regions of (Chr7: 34554890–34556263) and (Chr3: 17096021–17100322) with higher H3K27ac enrichment in LW was related with KCNK5 and TGFB1I1, respectively.

Candidate Genes Associated With Breed-Specific Regulatory Elements

We performed gene ontology enrichment analysis on the genes that overlapped with the H3K4me3 peak regions that showed breed-specific activity (334 for BMX and 169 for LW). These genes were enriched in melanosome organization (P = 0.004, AP3B1, HPS4, and ZEB2), regulation of T cell differentiation (P = 0.018, ADAM8, GTF2H3, and SART1), and activated T cell proliferation (P = 0.016, GPAM, RIPK3, and SATB1) (Supplementary Table 2). We displayed two peaks that displayed higher activity in Large White and Bama Xiang pigs. These included a peak (Chr15:111509199–111511820) located at the intron of PTH2R, a gene that is expressed in nervous system and plays a role in modulating pituitary functions (Della Penna et al., 2003) (Supplementary Figures 5A,C), and a peak in intron of ZEB2 (Supplementary Figures 5B,D), which encodes a zinc finger protein that is involved in the regulation of pituitary GH (Kumar et al., 2010). According to published RNA-Seq data, this gene also shows higher expression in the pituitary of Bama Xiang pigs than in Large White pigs (Yang et al., 2018).

For H3K27ac peaks, we found that 951 and 757 genes overlapped with BMX- and LW-specific peak regions, respectively. These genes were enriched in the regulation of growth (P = 0.0003, ABL1, ACVR1, and BGLAP), response to GH (P = 0.03, GHR, GHRL, and PFDN1), the thyroid hormone signaling pathway (P = 0.0003, ACTG1, MAPK1, and TP53), and the Notch signaling pathway (P = 0.01, ADAM10, CDH6, and TSPAN5) (Supplementary Table 3). These results match the function of the pituitary, which secretes GH and thyroid stimulating hormone that regulates the thyroid hormone (Brent, 2012). The Notch signaling pathway regulates endocrine cell specifications in the anterior pituitary (Dutta et al., 2008). We provide examples indicating that the peak (Chr13:75916467–75920943) shows higher H3K27ac activity in BMX, which is located close to EPHB1 gene (Figure 3E) that plays an important role in the regulation of synapse formation and maturation, migration of neural progenitors, and development of the immune organs (Wei et al., 2017), and locates in the gonadotropes of the pituitary gland, which plays an important role in various niches (Yoshida et al., 2017). Another peak Chr2:25338342–25347879, with higher H3K27ac activity in the BMX, was adjacent to the FJX1 gene (Figure 3F), which is expressed in telencephalon neurons and has an important function in neuron growth and the differentiation of limbs (Ashery-Padan et al., 1999) and as a direct target of Notch signaling (Rock et al., 2005). Among the H3K27ac peaks that show higher activity in LW than in BMX, we found that Chr7:34554890–34556363 is close to KCNK5, which participates in the regulation of platelet size and maturity (Christiansen et al., 2017) (Figure 3G); likewise, Chr3:17096021–17100322 is near TGFB1I1, which is involved in male sexual differentiation, cell proliferation, and vascular development (Wang et al., 2011) (Figure 3H) and the regulation of follicle-stimulating hormone in pituitary gonadotropes (Gore et al., 2005). The pituitary regulates the growth and sexual maturation of individuals by secreting GH and gonadotropin (Girard et al., 1999; Themmen and Huhtaniemi, 2000; Kumar, 2016).

Enrichment of Transcription Factor-Binding Motifs in Breed-Specific H3K27ac Regions

As TF plays an essential role in triggering epigenetic reprogramming (Sindhu et al., 2012), next we investigated the enrichment of TF binding motif in peaks that showed differential activity in H3K27ac between the two breeds using the HOMER program⁹ (Heinz et al., 2010). In total, we found 31 significantly enriched TF-binding motifs (q ≤ 0.05) in LW-specific peak regions (Supplementary Table 4), while no significant enrichment was observed in BMX-specific peaks. Although the existence of TF-binding motifs does not guarantee binding of TFs, of these enriched TFs, ATF3 (Hunt et al., 2012) has been linked to brain development and the nervous system; BATF (Williams et al., 2003), JUN (Garcia et al., 2009), and JunB (Szremska et al., 2003) may be related to the immune system, and ESRRG is a sequence motif of a susceptibility gene for obesity-related traits (Dong et al., 2016).

Diversity of Super Enhancer Among Individuals and Breeds

Super-enhancer elements typically include multiple enhancer elements and represent major drivers of transcriptional activation (Lovén et al., 2013). In this study, we identified 1,192, 1,299, 1,017, and 859 super enhancers in four individuals, respectively, using the ROSE algorithm (Whyte et al., 2013). After merging the peaks identified in the four individuals using bedtools, we obtained 2,025 non-redundant peaks. Among these, 302 were shared in all individuals under study (Figure 4A). The genes covered by the 302 super enhancers were enriched in the regulation of the Wnt signaling pathway (P = 0.004), central nervous system neuron development (P = 0.004), and the thyroid hormone signaling pathway (P = 0.01) (Figure 4B and Supplementary Table 5). These pathways match the function of the pituitary. Furthermore, using the DESeq2 program, we identified one and three super enhancers that showed higher activity in Bama Xiang and Large White pigs, respectively (Figure 4C). Chr15:109448708–109526513 and Chr13:143955853–144005651 (Figures 4D,E) are representative examples that display higher activity in Large White and Bama Xiang pigs, and the region of the two super enhancers covers the genes GPR1, EEF1B2, NDUFS1, and LSAMP; these can be used as potential targets for neuropsychiatric disorders (Innos et al., 2013).

FIGURE 4

Figure 4. Analysis of super enhancers. (A) The distribution of super enhancers in the four individuals. (B) Enriched Gene Ontology terms for genes covered by the 302 super enhancers shared cross all samples. (C) Volcano plot showing the differential activity analysis of super enhancer between two breeds. (D,E) The tracks of H3K27ac activity and gene in the two representative super enhancers that show differential activity between two breeds for reference (Chr15:109448708–109526513 and Chr13:143955853–144005651).

Conclusion

To the best of our knowledge, we present the first comprehensive exploration of regulatory elements in the pituitary glands of pigs using ChIP-Seq target H3K4me3 and H3K27ac, two histone modifications with important roles in regulating gene expression. We identified 65,044 non-redundant cis-regulatory sequences, including 23,680 in the H3K4me3 region (putative promoters) and 61,791 in the histone H3K27ac region, which include 12,318 proximal peaks (putative promoters) and 49,473 distal peaks (putative enhancers). It is worth noting that there are 86 H3K4me3 breed-specific regions that overlap with breed-specific regions of H3K27ac. This study enlarges the database of regulatory elements in the pituitary glands of pigs. The identified peaks are valuable for annotating the putative functional mutations identified in genome wide association studies of complex traits in pigs, as well as cross-species comparative studies.

Data Availability Statement

The ChIP-Seq data and processed data of BMX and LW individuals analyzed for this study can be found at Gene Expression Omnibus (GEO) with accession code: GSE178380 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE178380). The RNA sequencing data of pituitary glands downloaded from three Bama Xiang pigs and three Large White pigs from gsa.big.ac.cn under the GSA number of CRA000876 (https://ngdc.cncb.ac.cn/gsa/browse/CRA000876) (Yang et al., 2018).

Ethics Statement

The animal study was reviewed and approved by Committee on Animal Biosafety of Jiangxi Agricultural University. Written informed consent was obtained from the owners for the participation of their animals in this study.

Author Contributions

WL and BY conceived and designed the experiments. ZiZ, ZeZ, and TJ performed the experiments. ZiZ, BY, YZ, TJ, and ZL analyzed the data. YZ, BY, and ZiZ wrote the manuscript. All authors read and approved the final manuscript.

Funding

BY was supported by the National Natural Science Foundation of China (31760657).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank Prof. Lusheng Huang for his contribution in the experimental design and manuscript revision. We thank sample collectors in State Key Laboratory of Pig Genetic Improvement and Production Technology, Jiangxi Agricultural University for their great contributions.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.697994/full#supplementary-material

Footnotes

1. ^ https://www.faang.org
2. ^ https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE178380
3. ^ http://hgdownload.cse.ucsc.edu/goldenPath/susScr11/bigZips/
4. ^ https://ngdc.cncb.ac.cn/gsa/browse/CRA000876
5. ^ https://gat.readthedocs.io/en/Latest/tutorialGenomicAnnotation.html
6. ^ https://david.ncifcrf.gov/home.jsp
7. ^ http://homer.ucsd.edu/homer/motif/motifDatabase.html
8. ^ https://david.ncifcrf.gov/
9. ^ http://homer.ucsd.edu/homer/motif/motifDatabase.html

References

Abuin, J. M., Pichel, J. C., Pena, T. F., and Amigo, J. (2015). BigBWA: approaching the Burrows-Wheeler aligner to Big Data technologies. Bioinformatics 31, 4003–4005. doi: 10.1093/bioinformatics/btv506

PubMed Abstract | CrossRef Full Text | Google Scholar

Andersson, L., Archibald, A. L., Bottema, C. D., Brauning, R., Burgess, S. C., Burt, D. W., et al. (2015). Coordinated international action to accelerate genome-to-phenome with FAANG, the Functional Annotation of Animal Genomes project. Genome Biol. 16:57. doi: 10.1186/s13059-015-0622-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Ashery-Padan, R., Alvarez-Bolado, G., Klamt, B., Gessler, M., and Gruss, P. (1999). Fjx1, the murine homologue of the Drosophila four-jointed gene, codes for a putative secreted protein expressed in restricted domains of the developing and adult brain. Mech. Dev. 80, 213–217.

Google Scholar

Barski, A., Cuddapah, S., Cui, K., Roh, T.-Y., Schones, D. E., Wang, Z., et al. (2007). High-Resolution Profiling of Histone Methylations in the Human Genome. Cell 129, 823–837. doi: 10.1016/j.cell.2007.05.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Bindea, G., Mlecnik, B., Hackl, H., Charoentong, P., Tosolini, M., Kirilovsky, A., et al. (2009). ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics 25, 1091–1093. doi: 10.1093/bioinformatics/btp101

PubMed Abstract | CrossRef Full Text | Google Scholar

Brent, G. A. (2012). Mechanisms of thyroid hormone action. J. Clin. Invest. 122, 3035–3043. doi: 10.1172/JCI60047

PubMed Abstract | CrossRef Full Text | Google Scholar

Brinkmeier, M. L., Davis, S. W., Carninci, P., MacDonald, J. W., Kawai, J., Ghosh, D., et al. (2009). Discovery of transcriptional regulators and signaling pathways in the developing pituitary gland by bioinformatic and genomic approaches. Genomics 93, 449–460. doi: 10.1016/j.ygeno.2008.11.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Christiansen, M. K., Larsen, S. B., Nyegaard, M., Neergaard-Petersen, S., Wurtz, M., Grove, E. L., et al. (2017). The SH2B3 and KCNK5 loci may be implicated in regulation of platelet count, volume, and maturity. Thromb. Res. 158, 86–92. doi: 10.1016/j.thromres.2017.08.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Corradin, O., and Scacheri, P. C. (2014). Enhancer variants- evaluating functions in common disease. Genome Med. 6:85. doi: 10.1186/s13073-014-0085-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Della Penna, K., Kinose, F., Sun, H., Koblan, K. S., and Wang, H. (2003). Tuberoinfundibular peptide of 39 residues (TIP39): molecular structure and activity for parathyroid hormone 2 receptor. Neuropharmacology 44, 141–153. doi: 10.1016/s0028-3908(02)00335-0

CrossRef Full Text | Google Scholar

Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., et al. (2013). STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21. doi: 10.1093/bioinformatics/bts635

PubMed Abstract | CrossRef Full Text | Google Scholar

Dong, S. S., Guo, Y., Zhu, D. L., Chen, X. F., Wu, X. M., Shen, H., et al. (2016). Epigenomic elements analyses for promoters identify ESRRG as a new susceptibility gene for obesity-related traits. Int. J. Obes. 40, 1170–1176. doi: 10.1038/ijo.2016.44

PubMed Abstract | CrossRef Full Text | Google Scholar

Dutta, S., Dietrich, J. E., Westerfield, M., and Varga, Z. M. (2008). Notch signaling regulates endocrine cell specification in the zebrafish anterior pituitary. Dev. Biol. 319, 248–257. doi: 10.1016/j.ydbio.2008.04.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Ernst, J., Kheradpour, P., Mikkelsen, T. S., Shoresh, N., Ward, L. D., Epstein, C. B., et al. (2011). Mapping and analysis of chromatin state dynamics in nine human cell types. Nature 473, 43–49. doi: 10.1038/nature09906

PubMed Abstract | CrossRef Full Text | Google Scholar

Frantz, L. A., Schraiber, J. G., Madsen, O., Megens, H.-J., Bosse, M., Paudel, Y., et al. (2013). Genome sequencing reveals fine scale diversification and reticulation history during speciation in Sus. Genome Biol. 14:R107. doi: 10.1186/gb-2013-14-9-r107

PubMed Abstract | CrossRef Full Text | Google Scholar

Garcia, C. A., Wang, H., Benakanakere, M. R., Barrett, E., Kinane, D. F., and Martin, M. (2009). c-jun controls the ability of IL-12 to induce IL-10 production from human memory CD4+ T cells. J. Immunol. 183, 4475–4482. doi: 10.4049/jimmunol.0901283

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerstein, M. B., Lu, Z. J., Van Nostrand, E. L., Cheng, C., Arshinoff, B. I., Liu, T., et al. (2010). Integrative analysis of the Caenorhabditis elegans genome by the modENCODE project. Science 330, 1775–1787. doi: 10.1126/science.1196914

PubMed Abstract | CrossRef Full Text | Google Scholar

Girard, N., Boulanger, L., Denis, S., and Gaudreau, P. (1999). Differential in vivo regulation of the pituitary growth hormone-releasing hormone (GHRH) receptor by GHRH in young and aged rats. Endocrinology 140, 2836–2842. doi: 10.1210/endo.140.6.6760

PubMed Abstract | CrossRef Full Text | Google Scholar

Gong, H., Xiao, S., Li, W., Huang, T., Huang, X., Yan, G., et al. (2019). Unravelling the genetic loci for growth and carcass traits in Chinese Bamaxiang pigs based on a 1.4 million SNP array. J. Anim. Breed Genet. 136, 3–14. doi: 10.1111/jbg.12365

PubMed Abstract | CrossRef Full Text | Google Scholar

Gore, A. J., Philips, D. P., Miller, W. L., and Bernard, D. J. (2005). Differential regulation of follicle stimulating hormone by activin A and TGFB1 in murine gonadotropes. Reprod. Biol. Endocrinol. 3:73. doi: 10.1186/1477-7827-3-73

PubMed Abstract | CrossRef Full Text | Google Scholar

Heintzman, N. D., and Ren, B. (2009). Finding distal regulatory elements in the human genome. Genet. Dev. 19, 541–549. doi: 10.1016/j.gde.2009.09.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Heinz, S., Benner, C., Spann, N., Bertolino, E., Lin, Y. C., Laslo, P., et al. (2010). Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol. Cell 38, 576–589. doi: 10.1016/j.molcel.2010.05.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Hunt, D., Raivich, G., and Anderson, P. N. (2012). Activating transcription factor 3 and the nervous system. Front. Mol. Neurosci. 5:7. doi: 10.3389/fnmol.2012.00007

PubMed Abstract | CrossRef Full Text | Google Scholar

Innos, J., Koido, K., Philips, M.-A., and Vasar, E. (2013). Limbic system associated membrane protein as a potential target for neuropsychiatric disorders. Front. Pharmacol. 4:32. doi: 10.3389/fphar.2013.00032

PubMed Abstract | CrossRef Full Text | Google Scholar

Kern, C., Wang, Y., Xu, X., Pan, Z., Halstead, M., Chanthavixay, G., et al. (2021). Functional annotations of three domestic animal genomes provide vital resources for comparative and agricultural research. Nat. Commun. 12:1821. doi: 10.1038/s41467-021-22100-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Kiezun, M., Smolinska, N., Maleszka, A., Dobrzyn, K., Szeszko, K., and Kaminski, T. (2014). Adiponectin expression in the porcine pituitary during the estrous cycle and its effect on LH and FSH secretion. Am. J. Physiol. Endocrinol. Metab. 307, E1038–E1046. doi: 10.1152/ajpendo.00299.2014

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumar, P. A., Kotlyarevska, K., Dejkhmaron, P., Reddy, G. R., Lu, C., Bhojani, M. S., et al. (2010). Growth hormone (GH)-dependent expression of a natural antisense transcript induces zinc finger E-box-binding homeobox 2 (ZEB2) in the glomerular podocyte: a novel action of gh with implications for the pathogenesis of diabetic nephropathy. J. Biol. Chem. 285, 31148–31156. doi: 10.1074/jbc.M110.132332

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumar, T. R. (2016). Mouse Models for the Study of Synthesis, Secretion, and Action of Pituitary Gonadotropins. Prog. Mol. Biol. Transl. Sci. 143, 49–84. doi: 10.1016/bs.pmbts.2016.08.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Landt, S. G., Marinov, G. K., Kundaje, A., Kheradpour, P., Pauli, F., Batzoglou, S., et al. (2012). ChIP-seq guidelines and practices of the ENCODE and modENCODE consortia. Genome Res. 22, 1813–1831. doi: 10.1101/gr.136184.111

PubMed Abstract | CrossRef Full Text | Google Scholar

Larson, G., Dobney, K., Albarella, U., Fang, M., Matisoo-Smith, E., Robins, J., et al. (2005). Worldwide phylogeography of wild boar reveals multiple centers of pig domestication. Science 307, 1618–1621. doi: 10.1126/science.1106927

PubMed Abstract | CrossRef Full Text | Google Scholar

Liao, Y., Smyth, G. K., and Shi, W. (2014). featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30, 923–930. doi: 10.1093/bioinformatics/btt656

PubMed Abstract | CrossRef Full Text | Google Scholar

Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15:550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Lovén, J., Hoke Heather, A., Lin Charles, Y., Lau, A., Orlando David, A., Vakoc Christopher, R., et al. (2013). Selective Inhibition of Tumor Oncogenes by Disruption of Super-Enhancers. Cell 153, 320–334. doi: 10.1016/j.cell.2013.03.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Meurens, F., Summerfield, A., Nauwynck, H., Saif, L., and Gerdts, V. (2012). The pig: a model for human infectious diseases. Trends Microbiol. 20, 50–57. doi: 10.1016/j.tim.2011.11.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Meyer, R. D., Laz, E. V., Su, T., and Waxman, D. J. (2009). Male-specific hepatic Bcl6: growth hormone-induced block of transcription elongation in females and binding to target genes inversely coordinated with STAT5. Mol. Endocrinol. 23, 1914–1926. doi: 10.1210/me.2009-0242

PubMed Abstract | CrossRef Full Text | Google Scholar

Pennacchio, L. A., Bickmore, W., Dean, A., Nobrega, M. A., and Bejerano, G. (2013). Enhancers: five essential questions. Nat. Rev. Genet. 14, 288–295. doi: 10.1038/nrg3458

PubMed Abstract | CrossRef Full Text | Google Scholar

Quinlan, A. R., and Hall, I. M. (2010). BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26, 841–842. doi: 10.1093/bioinformatics/btq033

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramirez, F., Ryan, D. P., Gruning, B., Bhardwaj, V., Kilpert, F., Richter, A. S., et al. (2016). deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Res. 44, W160–W165. doi: 10.1093/nar/gkw257

PubMed Abstract | CrossRef Full Text | Google Scholar

Rice, P., Longden, I., and Bleasby, A. (2000). EMBOSS: the European Molecular Biology Open Software Suite. Trends Genet. 16, 276–277.

Google Scholar

Rock, R., Heinrich, A. C., Schumacher, N., and Gessler, M. (2005). Fjx1: a notch-inducible secreted ligand with specific binding sites in developing mouse embryos and adult brain. Dev. Dyn. 234, 602–612. doi: 10.1002/dvdy.20553

PubMed Abstract | CrossRef Full Text | Google Scholar

Rubin, C. J., Megens, H. J., Martinez Barrio, A., Maqbool, K., Sayyab, S., Schwochow, D., et al. (2012). Strong signatures of selection in the domestic pig genome. Proc. Natl. Acad. Sci. 109, 19529–19536. doi: 10.1073/pnas.1217149109

PubMed Abstract | CrossRef Full Text | Google Scholar

Shan, L., Wu, Q., Li, Y., Shang, H., Guo, K., Wu, J., et al. (2014). Transcriptome profiling identifies differentially expressed genes in postnatal developing pituitary gland of miniature pig. DNA Res. 21, 207–216. doi: 10.1093/dnares/dst051

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, Y., Yue, F., McCleary, D. F., Ye, Z., Edsall, L., Kuan, S., et al. (2012). A map of the cis-regulatory sequences in the mouse genome. Nature 488, 116–120. doi: 10.1038/nature11243

PubMed Abstract | CrossRef Full Text | Google Scholar

Sindhu, C., Samavarchi-Tehrani, P., and Meissner, A. (2012). Transcription factor-mediated epigenetic reprogramming. J. Biol. Chem. 287, 30922–30931. doi: 10.1074/jbc.R111.319046

PubMed Abstract | CrossRef Full Text | Google Scholar

Szremska, A. P., Kenner, L., Weisz, E., Ott, R. G., Passegue, E., Artwohl, M., et al. (2003). JunB inhibits proliferation and transformation in B-lymphoid cells. Blood 102, 4159–4165. doi: 10.1182/blood-2003-03-0915

PubMed Abstract | CrossRef Full Text | Google Scholar

Themmen, A. P. N., and Huhtaniemi, I. T. (2000). Mutations of gonadotropins and gonadotropin receptors: elucidating the physiology and pathophysiology of pituitary-gonadal function. Endocr. Rev. 21, 551–583. doi: 10.1210/edrv.21.5.0409

PubMed Abstract | CrossRef Full Text | Google Scholar

Tichomirowa, M. A., Lee, M., Barlier, A., Daly, A. F., Marinoni, I., Jaffrain-Rea, M. L., et al. (2012). Cyclin-dependent kinase inhibitor 1B (CDKN1B) gene variants in AIP mutation-negative familial isolated pituitary adenoma kindreds. Endocr. Relat. Cancer 19, 233–241. doi: 10.1530/ERC-11-0362

PubMed Abstract | CrossRef Full Text | Google Scholar

Villar, D., Berthelot, C., Aldridge, S., Rayner, T. F., Lukk, M., Pignatelli, M., et al. (2015). Enhancer evolution across 20 mammalian species. Cell 160, 554–566. doi: 10.1016/j.cell.2015.01.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Vodicka, P., Smetana, K. Jr., Dvorankova, B., Emerick, T., Xu, Y. Z., Ourednik, J., et al. (2005). The miniature pig as an animal model in biomedical research. Ann. N. Y. Acad. Sci. 1049, 161–171. doi: 10.1196/annals.1334.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., Hu, G., Betts, C., Harmon, E. Y., Keller, R. S., Water, L. V. D., et al. (2011). Transforming growth factor-β1-induced transcript 1 protein, a novel marker for smooth muscle contractile phenotype, is regulated by serum response factor/myocardin protein. J. Biol. Chem. 286, 41589–41599,

Google Scholar

Wei, W., Wang, H., and Ji, S. (2017). Paradoxes of the EphB1 receptor in malignant brain tumors. Cancer Cell Int. 17:21. doi: 10.1186/s12935-017-0384-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Whyte, W. A., Orlando, D. A., Hnisz, D., Abraham, B. J., Lin, C. Y., Kagey, M. H., et al. (2013). Master Transcription Factors and Mediator Establish Super-Enhancers at Key Cell Identity Genes. Cell 153, 307–319. doi: 10.1016/j.cell.2013.03.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Williams, K. L., Zullo, A. J., Kaplan, M. H., Brutkiewicz, R. R., Deppmann, C. D., Vinson, C., et al. (2003). BATF transgenic mice reveal a role for activator protein-1 in NKT cell development. J. Immunol. 170, 2417–2426. doi: 10.4049/jimmunol.170.5.2417

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Y., Adeola, A. C., Xie, H. B., and Zhang, Y. P. (2018). Genomic and transcriptomic analyses reveal selection of genes for puberty in Bama Xiang pigs. Zool. Res. 39, 424–430. doi: 10.24272/j.issn.2095-8137.2018.068

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoshida, S., Kato, T., Kanno, N., Nishimura, N., Nishihara, H., Horiguchi, K., et al. (2017). Cell type-specific localization of Ephs pairing with ephrin-B2 in the rat postnatal pituitary gland. Cell Tissue Res. 370, 99–112. doi: 10.1007/s00441-017-2646-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Liu, T., Meyer, C. A., Eeckhoute, J., Johnson, D. S., and Bernstein, B. E. (2008). Model-based Analysis of ChIP-Seq (MACS). Genome Biol. 9:R137. doi: 10.1186/gb-2008-9-9-r137

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, J., Cheng, F., Jia, P., Cox, N., Denny, J. C., and Zhao, Z. (2018). An integrative functional genomics framework for effective identification of novel regulatory variants in genome-phenome studies. Genome Med. 10:7. doi: 10.1186/s13073-018-0513-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, Y., Hou, Y., Xu, Y., Luan, Y., Zhou, H., Qi, X., et al. (2021). A compendium and comparative epigenomics analysis of cis-regulatory elements in the pig genome. Nat. Commun. 12:2217. doi: 10.1038/s41467-021-22448-x

PubMed Abstract | CrossRef Full Text | Google Scholar