Functionally Annotating Regulatory Elements in the Equine Genome Using Histone Mark ChIP-Seq

One of the primary aims of the Functional Annotation of ANimal Genomes (FAANG) initiative is to characterize tissue-specific regulation within animal genomes. To this end, we used chromatin immunoprecipitation followed by sequencing (ChIP-Seq) to map four histone modifications (H3K4me1, H3K4me3, H3K27ac, and H3K27me3) in eight prioritized tissues collected as part of the FAANG equine biobank from two thoroughbred mares. Data were generated according to optimized experimental parameters developed during quality control testing. To ensure that we obtained sufficient ChIP and successful peak-calling, data and peak-calls were assessed using six quality metrics, replicate comparisons, and site-specific evaluations. Tissue specificity was explored by identifying binding motifs within unique active regions, and motifs were further characterized by gene ontology (GO) and protein–protein interaction analyses. The histone marks identified in this study represent some of the first resources for tissue-specific regulation within the equine genome. As such, these publicly available annotation data can be used to advance equine studies investigating health, performance, reproduction, and other traits of economic interest in the horse.


Introduction
In 1992, researchers discovered the first disease mutation in horses, conferring hyperkalemic periodic paralysis (HYPP) in Quarter Horses [1], yet identification of additional equine genetic diseases progressed slowly, with only nine disease-associated variants discovered prior to 2007 [2, 3]. Since the release of the equine reference genome in 2007 [4], 22 additional genes were found to cause or be associated with equine diseases, yet there are at least 200 described genetic disorders for which causal variants are unknown [2,5]. While the majority of characterized equine disease variants are located within coding regions, an increasing amount of research in humans and other animal species suggests that a large number of disease mutations are harbored within regulatory elements and other functional 80 tissues, six fluids, and two cell lines that were collected from two healthy adult thoroughbred mares [26]. The eight tissues evaluated in this study were prioritized for thorough investigation based on (1) continuity with other FAANG species to enable across species comparisons and/or (2) the primary needs of the equine community in terms of health, performance, and reproduction.

Materials and Methods
Due to the limited number of previous histone ChIP-Seq experiments across equine tissues, we first performed quality control (QC) testing to determine appropriate experimental parameters for each tissue and mark. Sequencing, rather than targeted qPCR, was used to assess the quality of the QC data due to the limited knowledge of appropriate tissue-specific control genes in the horse. All ChIP-Seq experiments, including QC, were conducted by Diagenode ChIP-Seq Profiling Service (Diagenode, Cat# G02010000, Liège, Belgium), and a complete summary of the final protocols used for all tissues can be accessed at ftp://ftp.faang.ebi.ac.uk/ftp/protocols/assays/. Quality control involved chromatin extraction, ChIP, library preparation, and sequencing of one training sample from each of the eight tissues investigated. Quality control training samples were obtained from previously banked tissues collected at the University of California, Davis. Bioinformatic analysis was performed on the QC data in order to calculate library complexity and ChIP enrichment metrics and to call peaks to evaluate the genomic distribution of detected marks. Library complexity metrics included non-redundant fraction (NRF) and PCR bottleneck coefficients 1 and 2 (PBC1 and PBC2). Metrics for ChIP enrichment included normalized strand cross-correlation coefficient (NSC), relative strand cross-correlation coefficient (RSC), and Jensen-Shannon distance (JSD). Non-redundant fraction, PBC1 and 2, NSC, and RSC are all standardized metrics of the ENCODE project [27] and were compared against ENCODE standards to determine the efficacy of the ChIP protocols. Jensen-Shannon distance is a common statistic used to compare two distributions that can be applied to ChIP datasets using deepTools version 2.4.3 [28], and a threshold was determined by agreement among FAANG collaborators. Tissues of interest (TOI) used in the final experiments that generated combined peak-calls were collected from two thoroughbred mares (referred to in this manuscript as ECA_UCD_AH1 for SAMEA104728862 and ECA_UCD_AH2 for SAMEA104728877) as part of the FAANG equine biobank [26], and all protocols for this work were approved by the University of California, Davis Institutional Animal Care and Use Committee (Protocol #19037). Laboratory procedures that varied based on tissue are summarized in Table 1.

Chromatin Extraction
Chromatin was extracted from adipose using the True MicroChIP Kit (Diagenode, Cat# C01010130, Liège, Belgium) and from the other seven tissues following the iDeal ChIP-Seq Kit for histones (Diagenode, Cat# C01010059, Liège, Belgium) with the modifications or specifications described in this paper. The starting amount varied depending on tissue, such that those with extensive extracellular matrices and/or low ratio of nuclei to cellular matter required larger amounts of starting material compared to those tissues that homogenized easily. Samples were homogenized either by douncing (liver) or grinding with the Tissue Lyser II (Qiagen, Germany) at 25 strokes/minute for a length of time that varied by tissue (Table 1).
In order to reach the desired fragment length (approximately 200 bp), chromatin was sheared with the Bioruptor Pico (Diagenode, Cat# B01060001, Liège, Belgium) combined with the Bioruptor ® Water cooler for 8-12 cycles of 30 s with 30 s rest between cycles. The temperature was maintained at 10 • C for adipose and 4 • C for all other tissues during shearing. The number of cycles varied based on tissue (Table 1), and the chromatin quality was assessed using the Fragment Analyzer (Aligent, Santa Clara, CA, USA).

Immunoprecipitation
Immunoprecipitation (IP) of the four histone marks, along with a negative IP control (IgG), was performed on tissue-specific amounts of chromatin using the IP-Star Compact Automated System (Diagenode, Cat# B03000002, Liège, Belgium). The antibodies used were all previously validated by Diagenode, and antibody concentrations were determined during QC for every tissue and histone mark combination (Table S1). An aliquot of chromatin from each tissue was set aside for the input to characterize sequencing background and identify true ChIP enrichment.

Sequencing
Libraries were prepared using the IP-Star ® Compact Automated System (Diagenode, Cat# B03000002, Liège, Belgium) and the MicroPlex Library Preparation Kit v2 (Diagenode, Cat# C05010013, Liège, Belgium) for the input and four ChIPs per tissue. Libraries were amplified prior to sequencing for at least 10 cycles, and additional cycles were performed as needed to reach a concentration of 3-10 nM. Using Agencourt AMPure XP (Bechman Coulter, Brea, CA, USA), libraries were purified and fragments were size-selected for approximately 200 bp. Libraries were sequenced as 50 bp, single-end reads on the HiSeq 4000 platform (Illumina, San Diego, CA, USA) to generate approximately 55-80 M raw reads for H3K27me3 and 30-50 M raw reads for the other marks and inputs.

Data Processing
A complete summary of the bioinformatic workflow can be accessed at ftp://ftp.faang.ebi.ac. uk/ftp/protocols/analysis/, and bioinformatic parameters that varied by mark are summarized in Table 2. Reads were trimmed using Trim-Galore version 0.4.0 [29,30] under the default parameters and aligned to EquCab3.0 [31] with BWA-MEM version 0.7.9a [32], such that split hits were marked as secondary alignments. Alignments were converted to BAM file format, processed, and filtered using SAMtools version 1.9 [33]. For strict quality filtering, reads were removed if they did not map, had secondary alignments, failed platform/vendor quality tests, were identified as optical duplicates, or had an alignment quality score lower than 30. PCR duplicates were marked with PicardTools version 2.7.1 [34] and then removed with SAMtools. For peak-calling, MACS2 version 2.1.1.20160309 [35] was used to call peaks for all marks, and SICERpy version 0.1.1 [36], which is a wrapper for executing SICER [37], was also used to call peaks for the broad mark, H3K27me3. Combining peak-calls involved identifying overlapping regions of enrichment in both biological replicates where at least one replicate was significantly enriched based on a set of significance thresholds that varied by mark (Table 2). Additionally, enrichment tracks (bigWig files) were generated using deepTools version 2.4.3, which subtracted background characterized by the input and then combined enrichment from both replicates. Gap Size n/a n/a n/a 4 Window Size n/a n/a n/a 200 Genome Fraction n/a n/a n/a 0.63 1 Mark was treated as broader than other narrow marks due to being categorized previously as broad by ENCODE [19]. 2 SICERpy was only used to call peaks for the broad mark, H3K27me3.

Data Analysis
As with the QC data, the datasets from the eight TOI were assessed by calculating library complexity and ChIP enrichment metrics, as well as evaluating the genomic distribution of detected marks. Identity between peaks called for ECA_UCD_AH1 and ECA_UCD_AH2 were compared to assess the similarity of the biological replicates using the Jaccard Index [38], also known as the Jaccard Similarity Coefficient, from BEDtools version 2.27.1 [39]. Unique peaks were defined as a peak for a given mark that does not have any overlap with peaks from the same mark in the other prioritized tissues. BEDtools version 2.27.1 was used to identify unique peaks, as well as calculate the percent of the genome covered by peaks. Graphs were generated using ggplot2 with R software version 3.4.3 [40,41]. To characterize the average peak topology in relation to gene annotations and calculate normalized enrichment patterns, we used deepTools version 2.4.3. RNA-Seq data from the two FAANG replicates allowed for site-specific validation of the histone peaks (ERR2584116, ERR2584168, ERR2584153, ERR2584205, ERR2584194, ERR2584142, ERR2584135, ERR2584187, ERR2584195, ERR2584143, ERR2584197, ERR2584145, ERR2584144, ERR2584196, ERR2584152, and ERR2584204). Analysis of Motif Enrichment (AME) from the MEME Suite version 5.0.5 [42] was used to identify known transcription factor binding motifs within peaks based on the JASPAR 2018 vertebrate database [43], and Biological Process Gene Ontologies (GOs) from Swiss Prot were used to perform a GO term analysis [44]. Novel motifs were characterized with DREME and MEME, and each of the novel motifs was manually investigated to identify additional known motifs that were not included in the JASPAR database. The Integrated Genome Viewer [45] was used to visualize peak-calls in conjunction with the Ensembl annotation (release 95) for the EquCab3.0 reference assembly [46]. String version 11.0 was used to perform a protein-protein interaction analysis on the transcription factors implicated in each tissue based on the enriched motifs identified in the tissue-specific active enhancer elements [47].

Results
Quality control testing was performed to determine tissue-specific laboratory parameters such as antibody concentration and shearing time ( Table 1 and Table S1) by comparing the quality metrics and peak-calls to ENCODE and FAANG standards (Table S2). The raw and processed data are available on https://data.faang.org/home under the study accession PRJEB35307. The processed files include read alignments to EquCab3.0 and peak-calls for each biological replicate, as well as the combined peak-calls.

Assessing Data Quality
Using the Jaccard Index to compare replicates for all of the marks, we found the highest identity between the two replicates of the same tissue ( Figure 1), with the exception of the brain replicates for H3K4me1 and the ovary replicates for H3K27me3. For the brain replicates, ECA_UCD_AH2 had 65,327 peaks compared to 143,328 peaks for ECA_UCD_AH1 (Table 3). In addition to a lower peak number, two of the library complexity quality scores for the ECA_UCD_AH2 brain replicate were lower than the ENCODE quality thresholds (Table 4). In terms of enrichment, however, the cross-correlation metrics for this sample were both above the established thresholds and indicate that the data had sufficient ChIP enrichment. Indeed, we were able to call 95,918 combined peaks, which is consistent with the range of H3K4me1 values from the other tissues. For the ovary replicates, two of the three library complexity metrics were below ENCODE standards, yet all of the ChIP enrichment metrics were consistent with those for the broad mark in other tissues, indicating that we had sufficient ChIP enrichment despite lower library complexity. When comparing SICER and MACS2 peak-calls for the broad mark, the number of combined peaks (8479 and 40,825, respectively) and the percentage of the genome covered (2.1% and 1.2%, respectively) for ovary were all consistent with the same measures for the H3K27me3 peaks from the other tissues evaluated in this study.
Genes 2019, 10, x FOR PEER REVIEW 6 of 22 65,327 peaks compared to 143,328 peaks for ECA_UCD_AH1 (Table 3). In addition to a lower peak number, two of the library complexity quality scores for the ECA_UCD_AH2 brain replicate were lower than the ENCODE quality thresholds (Table 4). In terms of enrichment, however, the cross-correlation metrics for this sample were both above the established thresholds and indicate that the data had sufficient ChIP enrichment. Indeed, we were able to call 95,918 combined peaks, which is consistent with the range of H3K4me1 values from the other tissues. For the ovary replicates, two of the three library complexity metrics were below ENCODE standards, yet all of the ChIP enrichment metrics were consistent with those for the broad mark in other tissues, indicating that we had sufficient ChIP enrichment despite lower library complexity. When comparing SICER and MACS2 peak-calls for the broad mark, the number of combined peaks (8479 and 40,825, respectively) and the percentage of the genome covered (2.1% and 1.2%, respectively) for ovary were all consistent with the same measures for the H3K27me3 peaks from the other tissues evaluated in this study.   While demonstrating more variation than the other marks, H3K27me3 had peaks called by SICER that were more similar between replicates of the same tissue compared to the H3K27me3 peaks called by MACS2 software ( Figure S1). Peak-calls from the two biological replicates were combined by identifying regions of overlapping enrichment in which at least one replicate had significant enrichment based on a q-value that differed by mark (Table 2). For H3K4me1, H3K4me3, and H3K27ac, the number of combined peaks ranged from 93-121 K, 26-29 K, and 64-88 K, respectively. The number of combined peaks called for the broad mark was lower than the three activating marks for both MACS2 and SICER (24-68 K and 7-11 K, respectively). Although the combined peak numbers for the MACS2-H3K27me3 datasets were more similar to the ENCODE equivalent than the number of peaks called by SICER (Table S3), the SICER-H3K27me3 combined peaks covered a larger portion of the genome, similar to that for the other marks (Table 3). Files for H3K27me3 peaks called by MACS2 and SICER are both publicly available for every tissue.
As a proof-of-principle, we investigated a small number of regions near well-characterized genes to compare the histone peaks with RNA-seq data generated from the same tissues. We found consistent activating marks across all tissues for a widely expressed gene, ACTB ( Figure 2A). Conversely, liver was the only tissue enriched for a set of active histone marks near the transcription start site (TSS) of the liver-specific gene CYP2E1 ( Figure 2B) [48]. Table 4. Quality metrics for assessing library complexity and ChIP enrichment. Thresholds for NRF, PBC1, PBC2, NSC, and RSC represent those developed by the Encyclopedia of DNA Elements (ENCODE) [27]. JSD threshold was established among members of the Functional Annotation of ANimal Genomes (FAANG) consortium.  As a proof-of-principle, we investigated a small number of regions near well-characterized genes to compare the histone peaks with RNA-seq data generated from the same tissues. We found consistent activating marks across all tissues for a widely expressed gene, ACTB ( Figure 2A). Conversely, liver was the only tissue enriched for a set of active histone marks near the transcription start site (TSS) of the liver-specific gene CYP2E1 ( Figure 2B) [48]. Proof-of-principle investigating house-keeping gene, ACTB, and liver-specific gene, CYP2E1, for appropriate regulatory elements. For adipose, brain, and liver tissue, combined peaks are displayed for H3K4me1 (aqua), H3K4me3 (light blue), H3K27ac (dark blue), H3K27me3 from SICER (orange), and mRNA expression (purple). (A) ACTB is a housekeeping gene that is commonly expressed for many tissues. (B) CYP2E1 is a liver enzyme, which displays tissue-specific expression. Note the presence of the H3K27me3 repressive mark (orange) within adipose and brain samples.

Characterizing Tissue-Specific Features
Brain tissue had the highest percentage of unique peaks, defined as peaks that were only found in that tissue, for H3K27ac (31%) and H3K27me3 (20%), while liver had the highest percentage of Figure 2. Proof-of-principle investigating house-keeping gene, ACTB, and liver-specific gene, CYP2E1, for appropriate regulatory elements. For adipose, brain, and liver tissue, combined peaks are displayed for H3K4me1 (aqua), H3K4me3 (light blue), H3K27ac (dark blue), H3K27me3 from SICER (orange), and mRNA expression (purple). (A) ACTB is a housekeeping gene that is commonly expressed for many tissues. (B) CYP2E1 is a liver enzyme, which displays tissue-specific expression. Note the presence of the H3K27me3 repressive mark (orange) within adipose and brain samples.

Characterizing Tissue-Specific Features
Brain tissue had the highest percentage of unique peaks, defined as peaks that were only found in that tissue, for H3K27ac (31%) and H3K27me3 (20%), while liver had the highest percentage of unique peaks for H3K4me1 (32%) and H3K4me3 (16%), along with the second highest for H3K27ac (26%) and H3K27me3 (14%) (Figure 3). Lamina tissue also had a high percentage of unique peaks for the three activating marks with 24, 10, and 26% for H3K4me1, H3K4me3, and H3K27ac, respectively. unique peaks for H3K4me1 (32%) and H3K4me3 (16%), along with the second highest for H3K27ac (26%) and H3K27me3 (14%) (Figure 3). Lamina tissue also had a high percentage of unique peaks for the three activating marks with 24, 10, and 26% for H3K4me1, H3K4me3, and H3K27ac, respectively. In addition to characterized genes, we also investigated a small number of genomic regions with putative tissue-specific functions in liver and muscle. For liver, a potential tissue-specific regulatory element was identified in the 59th intron of PKHD1 (Ensembl Transcript ID: ENSECAT00000024985.1; Figure 4), a gene which has been previously associated with liver fibrosis [49]. Similarly, when considering a genomic region associated with racing ability [50], peaks for H3K4me3 in both muscle tissues were discovered at the start of a predicted lncRNA from Ensembl genebuild [51], indicating that this uncharacterized gene may be particularly informative for the function of contractile tissue ( Figure 5).
Across all tissues, H3K4me1 was enriched around the TSS, with a decrease in enrichment at the actual annotated start site that created a bimodal distribution across the average gene body ( Figure  6A). Additionally, H3K4me1 maintained a moderate level of enrichment throughout the gene body, as well as 3 Kb up-and downstream, as expected for distal and proximal enhancer elements. Alternatively, H3K4me3 ( Figure 6B) and H3K27ac ( Figure 6C) had peaks of enrichment at the TSS, although H3K27ac also showed enrichment, to a lesser extent, just upstream of the TSS. H3K27me3 had lower enrichment than the other three marks overall, but the average enrichment was essentially constant throughout the gene body, as well as 3 Kb up-and downstream ( Figure 6D). While patterns of enrichment for each mark were highly consistent between tissues, the enrichment of H3K27me3 for lamina was lower at the TSS compared to the other tissues. In addition to characterized genes, we also investigated a small number of genomic regions with putative tissue-specific functions in liver and muscle. For liver, a potential tissue-specific regulatory element was identified in the 59th intron of PKHD1 (Ensembl Transcript ID: ENSECAT00000024985.1; Figure 4), a gene which has been previously associated with liver fibrosis [49]. Similarly, when considering a genomic region associated with racing ability [50], peaks for H3K4me3 in both muscle tissues were discovered at the start of a predicted lncRNA from Ensembl genebuild [51], indicating that this uncharacterized gene may be particularly informative for the function of contractile tissue ( Figure 5).  Across all tissues, H3K4me1 was enriched around the TSS, with a decrease in enrichment at the actual annotated start site that created a bimodal distribution across the average gene body ( Figure 6A). Additionally, H3K4me1 maintained a moderate level of enrichment throughout the gene body, as well as 3 Kb up-and downstream, as expected for distal and proximal enhancer elements. Alternatively, H3K4me3 ( Figure 6B) and H3K27ac ( Figure 6C) had peaks of enrichment at the TSS, although H3K27ac also showed enrichment, to a lesser extent, just upstream of the TSS. H3K27me3 had lower enrichment than the other three marks overall, but the average enrichment was essentially constant throughout the gene body, as well as 3 Kb up-and downstream ( Figure 6D). While patterns of enrichment for each mark were highly consistent between tissues, the enrichment of H3K27me3 for lamina was lower at the TSS compared to the other tissues.

Identifying Motifs and Biological Process GO Terms
We identified between 16 and 61 transcription factor binding motifs in the unique active regions for each tissue, and the full results from the GO term analysis can be found in Tables S4-S19. While a large proportion of the identified transcription factors were still shared between tissues despite being identified in tissue-specific active regions, all of the tissues except for adipose and lung had at least one uniquely detected transcription factor binding site. A uniquely detected binding site was defined as a motif that was only identified in the tissue-specific active regions for a single tissue. The motifs were ranked based on significance of enrichment after multiple testing correction, and the top five enriched and identified motifs for each tissue are listed in Table 5. Ovary was the only tissue that had uniquely detected transcription factor motifs in this top five list, and the most significant motif was associated with FOXO3, which is a transcription factor (TF) characterized by 39 GO terms, including several for tissue-specific functions such as oocyte maturation (GO:0001556) and ovulation (GO:0001542). In fact, upon closer inspection, the FOXO3 motif was found within a tissue-specific regulatory element near a gene with recognized roles in mammalian reproduction, NR5A1 (Figure 7) [52]. Using a network analysis for TFs implicated in each tissue, we found that six networks contained EP300 as a central node, although these networks did not link every TF for a given tissue ( Figure 8A). Interestingly, MYC was the central node for brain, liver, and skeletal muscle ( Figure 8B), and TP53 was the central node for lung ( Figure 8C). Table 5. Top five significantly enriched transcription factor binding motifs identified in tissue-specific active enhancers. Tissue specificity of the active enhancers was defined by overlap of H3K4me1 and H3K27ac in the same tissue and no overlap of these marks in this region in any other tissues, and tissue specificity of the binding motifs was then defined by detection of an enriched binding motif in only one tissue.

Rank
Motif

Identifying Motifs and Biological Process GO Terms
We identified between 16 and 61 transcription factor binding motifs in the unique active regions for each tissue, and the full results from the GO term analysis can be found in Tables S4-S19. While a large proportion of the identified transcription factors were still shared between tissues despite being identified in tissue-specific active regions, all of the tissues except for adipose and lung had at least one uniquely detected transcription factor binding site. A uniquely detected binding site was defined as a motif that was only identified in the tissue-specific active regions for a single tissue. The motifs were ranked based on significance of enrichment after multiple testing correction, and the top five enriched and identified motifs for each tissue are listed in Table 5. Ovary was the only tissue that had uniquely detected transcription factor motifs in this top five list, and the most significant motif was associated with FOXO3, which is a transcription factor (TF) characterized by 39 GO terms, including several for tissue-specific functions such as oocyte maturation (GO:0001556) and ovulation (GO:0001542). In fact, upon closer inspection, the FOXO3 motif was found within a tissue-specific regulatory element near a gene with recognized roles in mammalian reproduction, NR5A1 (Figure 7) [52]. Using a network analysis for TFs implicated in each tissue, we found that six networks contained EP300 as a central node, although these networks did not link every TF for a given tissue ( Figure 8A). Interestingly, MYC was the central node for brain, liver, and skeletal muscle ( Figure 8B), and TP53 was the central node for lung ( Figure 8C). Table 5. Top five significantly enriched transcription factor binding motifs identified in tissue-specific active enhancers. Tissue specificity of the active enhancers was defined by overlap of H3K4me1 and H3K27ac in the same tissue and no overlap of these marks in this region in any other tissues, and tissue specificity of the binding motifs was then defined by detection of an enriched binding motif in only one tissue.

Rank
Motif

Discussion
As part of the FAANG consortium, we mapped more than 1 million putative regulatory sites across the equine genome, which will contribute significantly to our understanding of genome regulation in horses, as well as regulatory differences across species. To ensure that we obtained high quality data for the horse, we first performed QC experiments to optimize tissue-specific laboratory protocols (Table S2). Adipose tissue presented a distinct challenge for chromatin extraction due to the low number of nuclei and large amount of cellular material, including lipid deposits. Efforts to obtain high quality data were successful, as quality metrics for this tissue surpassed ENCODE standards or scored within the range of the other tissues for a given mark. From this work, we suspect that other difficult tissues (i.e., those with low nuclei density, extensive and persistent extra-or intracellular material, etc.) will also require similar alternate approaches for chromatin extraction, such as additional starting material or specialized kits.
To ensure accurate and relevant peak-calling, we employed MACS2 software for calling H3K4me1, H3K4me3, and H3K27ac peaks and both MACS2 and SICER for H3K27me3 peaks. While several attempts have been made to develop a bench-marking method for ranking ChIP peak-callers, there is no established gold standard for selecting a particular application [53]. Steinhauser et al., for example, found that SICER was better able to detect true-positive differential regions (DRs) when H3K36me3 data were down-sampled to as low as 10%, while MACS2 was only effective at detecting true-positive DRs at 60% or higher [54]. SICER, however, had a higher number of false-positive DRs compared to MACS2, indicating that it may be sacrificing some specificity to obtain higher sensitivity, while MACS2 does the opposite. Based upon generating peaks with higher identity between replicates ( Figure S1) and higher proportions of the genome covered (Table 3), we found SICER software to be more consistent for calling broad peak-calls compared to MACS2, which is

Discussion
As part of the FAANG consortium, we mapped more than 1 million putative regulatory sites across the equine genome, which will contribute significantly to our understanding of genome regulation in horses, as well as regulatory differences across species. To ensure that we obtained high quality data for the horse, we first performed QC experiments to optimize tissue-specific laboratory protocols (Table S2). Adipose tissue presented a distinct challenge for chromatin extraction due to the low number of nuclei and large amount of cellular material, including lipid deposits. Efforts to obtain high quality data were successful, as quality metrics for this tissue surpassed ENCODE standards or scored within the range of the other tissues for a given mark. From this work, we suspect that other difficult tissues (i.e., those with low nuclei density, extensive and persistent extra-or intracellular material, etc.) will also require similar alternate approaches for chromatin extraction, such as additional starting material or specialized kits.
To ensure accurate and relevant peak-calling, we employed MACS2 software for calling H3K4me1, H3K4me3, and H3K27ac peaks and both MACS2 and SICER for H3K27me3 peaks. While several attempts have been made to develop a bench-marking method for ranking ChIP peak-callers, there is no established gold standard for selecting a particular application [53]. Steinhauser et al., for example, found that SICER was better able to detect true-positive differential regions (DRs) when H3K36me3 data were down-sampled to as low as 10%, while MACS2 was only effective at detecting true-positive DRs at 60% or higher [54]. SICER, however, had a higher number of false-positive DRs compared to MACS2, indicating that it may be sacrificing some specificity to obtain higher sensitivity, while MACS2 does the opposite. Based upon generating peaks with higher identity between replicates ( Figure S1) and higher proportions of the genome covered (Table 3), we found SICER software to be more consistent for calling broad peak-calls compared to MACS2, which is consistent with previous research [55]. For that reason, we continued our investigations using the SICER-called H3K27me3 peaks. Given that peaks from both callers are supported in the literature for other species studied, both sets of peak-calls are available at https://data.faang.org/dataset/PRJEB35307.
Using the Jaccard Similarity Index to compare significant peak-calls, we determined that biological replicates were highly similar for each tissue with the exception of the brain replicates for H3K4me1 and the ovary replicates for H3K27me3 (Figure 1). While these low identity scores could be the result of underlying biological differences between the samples of brain or ovary tissue, this is unlikely given that the differences between replicates were only found in one mark for each of the two tissues. Using a bioinformatic method that relied on identifying overlapping enrichment rather than the intersection of significant peaks ensured accurate combined peak calling despite lower similarity scores. Combined peak-calls were assessed for significant enrichment surrounding genes with known expression patterns. In particular, the three active marks (H3K4me1, H3K4me3, and H3K27ac) had consistent peaks between tissues near ACTB ( Figure 2A) and tissue-specific peaks near the liver enzyme CYP2E1 (Figure 2B), which support the validity of our peak-calling methods. Additionally, these peak-calls were all consistent with RNA expression from the same tissues for ECA_UCD_AH1 and ECA_UCD_AH2, further supporting the validity of the enrichment patterns detected for each histone mark.
Comparing replicated peaks between tissues, we found that liver, brain, and lamina had the highest percentage of peaks that were unique to only that tissue ( Figure 3). Liver is known to have hundreds of distinct biological functions [48]. Since many of the functions are entirely unique to liver, a high proportion of unique active regulatory elements is consistent with the specialization of this tissue. Additionally, the mammalian brain is thought to have hundreds of different neuronal cell types within the cortex to coordinate numerous neurologic functions simultaneously, making it one of the most complex tissues of the body [56]. Therefore, these results are consistent with the expectation that there is a high degree of regulatory specificity to control numerous coordinated functions within complex tissues.
Interestingly, lamina was the other tissue with a high percentage of unique peaks for H3K4me1, H3K4me3, and H3K27ac ( Figure 3). Due to the role of lamina in disease, this finding may be particularly impactful for research into the physiological changes associated with laminitis, a syndrome that was established as a priority for equine research by the American Association of Equine Practitioners [57,58]. Additionally, when comparing the distribution of marks across the average gene body, we found that lamina tissue had an unusually low level of the repressive mark, H3K27me3, at the TSS ( Figure 6). While all of the tissues had a portion of genes that appear to have a dip in enrichment directly at the TSS, lamina had the most extreme decrease. Perhaps the dip in enrichment for lamina is an indication of increased levels of expression across more of the genome and, when combined with the high percentage of unique peaks, could suggest that hoof lamina may perform additional biological processes that are currently uncharacterized. In order to understand all of the molecular functions of lamina tissue and its role in laminitis, more work is needed to further annotate and functionally characterize the extent of the cellular processes within healthy and diseased tissue. Without further validation, however, we cannot exclude that this decrease in enrichment may be a technical artifact such as decreased detection of H3K27me3 in lamina tissue due to low cell numbers and excess extracellular matrix.
To further characterize tissue-specific regulatory regions, we identified numerous transcription factor binding motifs within the unique active elements for each tissue (Table 5). Despite investigating elements that were only enriched in one tissue, we still identified many transcription factor motifs that were shared between tissues (Table 4 and Table S19) and, based on a GO term analysis, many of these factors had numerous associated biological processes (Tables S4-S18). Our network analysis identified EP300 as a central secondary factor needed to connect many implicated TFs ( Figure 8A). This gene encodes a histone acetyltransferase that is known to interact with many TFs by protein-protein interactions rather than DNA binding [59,60]. Given that the TF motifs were all identified within active elements based on the presence of H3K27ac, finding a connection with EP300 supports our ability to detect relevant peaks for this histone modification. Interestingly, there were three tissues (brain, liver, and muscle) that had MYC rather than EP300 as the central node for their TF networks ( Figure 8B). We found that binding motifs for MYC were enriched in all three tissues, which is consistent with its role as a TF for many common cellular processes, including growth and regulation of the cell cycle [61]. Similarly, TP53 was the central node in the protein interaction network for lung tissue ( Figure 8C), which is consistent with its role as a key regulator of the cell cycle when functional and as a major tumorigenesis factor in lung cancer when dysfunctional [62].
By developing tissue-specific maps of the H3K4me1, H3K4me3, H3K27ac, and H3K27me3 histone modifications, we aimed to build upon our current knowledge about genomic regulation in the horse and provide new resources for further advancing research on health, performance, and reproduction traits. In particular, incorporating histone ChIP-seq data with current annotations is expected to lead to the discovery of more tissue-specific functional variants that are outside of protein-coding regions. For example, a genome-wide investigation of liver fibrosis in the Swiss Franches-Montagnes breed led to the identification of an associated region containing a promising candidate gene, PKHD1, which was also associated with kidney and liver disease in humans [49]. Sequencing the coding region of this gene identified two synonymous coding variants strongly associated with disease; however, their causal role remains unclear. Comparing peaks from multiple tissues in that region, there are also several liver-specific regulatory peaks within intron 59 of PKHD1 (Ensembl Transcript ID: ENSECAT00000024985.1) that may affect expression of this or another hepatic gene nearby ( Figure 4) and represent another avenue to explore for the molecular mechanism underlying liver fibrosis.
Tissue-specific histone modifications have also been implicated in the complex traits of performance [63] and, throughout modern horse domestication, performance abilities such as speed have been a major factor for selective breeding [64]. For example, recent research comparing three trotting breeds in Scandinavia (North Swedish Draught, Norwegian-Swedish Coldblooded trotter, and Standardbred trotter) found a 684 Kb genomic region associated with trotting racing ability containing numerous annotated and unannotated genes [50]. Using a protein-coding annotation, only one gene in the region, TRIM37, was identified as a potential candidate due to its associations with growth phenotypes in humans. Utilizing the locations of histone mark enrichment could help to identify other candidates such as tissue-specific regulatory regions or the genes that they regulate. In particular, there are H3K4me3 peaks at the start of a novel lncRNA predicted by Ensembl genebuild [51] within the region of interest for racing ability ( Figure 5A). While uncharacterized, LOC111775680 has a set of unique promoters found in skeletal and cardiac muscle with corresponding transcript data (as determined by RNAseq of the same tissues from ECA_UCD_AH1 and ECA_UCD_AH2), suggesting that the lncRNA may play a role in muscle physiology ( Figure 5B). Additional work is needed to further characterize LOC111775680 and investigate its potential role in muscle tissue and any corresponding effects on racing ability.
In addition to health and performance, these annotation data can also be used to identify important genomic regulatory regions for the reproductive system. Previously, homology with humans and mice was utilized to create a panel of candidate regions for determining the genetic cause of gonadal dysgenesis disorders across many breeds of horses [65], yet the panel was focused on coding variants and those within untranslated regions that may affect RNA and protein synthesis. By evaluating histone marks, we found that there are additional regulatory regions relevant to the reproductive system that were missed by the homology-based approaches. For example, NR5A1, a gene that is implicated in many dysgenesis cases [52,66], has a regulatory region as identified by H3K4me1 and H3K27ac histone modifications in equine ovary tissue that is characteristic of an active enhancer ( Figure 7A). Upon further inspection, we identified two FOXO3 binding sites within the region enriched with peak-calls ( Figure 7B,C). FOXO3 is a transcription factor with known roles in oocyte maturation [GO:0001556] and ovulation [GO:0001542] that was identified as an ovary-specific TF based on our motif analysis (Table S18).
Along with the regulatory regions identified in this study, the equine FAANG consortium is characterizing the full RNA profile of more than 50 tissues and cell types from the FAANG equine biobank supported by a research community initiative [26] and the DNA methylation profiles for the TOI. Moreover, the equine FAANG consortium is currently mapping a major insulator protein known as CTCF, with the goal of generating tissue-specific chromatin state predictions for the eight TOI. CTCF plays a key role in defining chromatin looping and, therefore, topologically associated domains and gene-enhancer interactions when combined with histone ChIP-Seq [67,68]. The same panel of genomic investigations are also being conducted in two stallions, so that males and females are represented in the annotations. We anticipate that the integration and utilization of these functional annotation datasets by the equine genomics community will lead to the identification of causal, non-coding variants underlying many traits of interest for equine medicine, performance, and reproduction.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4425/11/1/3/s1: Figure S1. Jaccard Index to measure similarity of peaks called for H3K27me3 peaks based on either MACS2 or SICER peak-calling; Table S1. Concentrations and catalog numbers for antibodies used during ChIP; Table S2. Quality metrics, peak-calling summary, and comparison to Encyclopedia of DNA Elements (ENCODE) data from quality control (QC) testing; Table S3. Comparing peak numbers to ENCODE data; Table S4. Motif and gene ontology (GO) term analysis results for adipose-specific active enhancers; Table S5. Motif and GO term analysis results for adipose-specific active promoters; Table S6. Motif and GO term analysis results for brain-specific active enhancers; Table S7. Motif and GO term analysis results for brain-specific active promoters; Table S8. Motif and GO term analysis results for heart-specific active enhancers; Table S9. Motif and GO term analysis results for heart-specific active promoters; Table S10. Motif and GO term analysis results for lamina-specific active enhancers; Table S11. Motif and GO term analysis results for lamina-specific active promoters; Table S12. Motif and GO term analysis results for liver-specific active enhancers; Table S13. Motif and GO term analysis results for liver-specific active promoters; Table S14. Motif and GO term analysis results for lung-specific active enhancers; Table S15. Motif and GO term analysis results for lung-specific active promoters; Table S16. Motif and GO term analysis results for muscle-specific active enhancers; Table S17. Motif and GO term analysis results for muscle-specific active promoters; Table S18. Motif and GO term analysis results for ovary-specific active enhancers; Table S19. Summary of the GO term analysis results for unique active elements from all tissues. Funding: Funding for experimental materials, and other resources was provided by the Grayson Jockey Club Foundation, USDA NRSP-8 equine species coordinator funds, and the Center for Equine Health with funds provided by the State of California pari-mutuel fund and contributions by private donors. Support for N.B.K. was also provided by the Grayson Jockey Club Foundation, Morris Animal Foundation (D16-EQ-028), and a fellowship from the Center for Equine Health at the University of California-Davis. Support for C.J.F. was provided by the National Institutes of Health (NIH) (K01OD015134 and L40 TR001136).