Genome-Wide Bovine H3K27me3 Modifications and the Regulatory Effects on Genes Expressions in Peripheral Blood Lymphocytes

Background Gene expression of lymphocytes was found to be influenced by histone methylation in mammals and trimethylation of lysine 27 on histone H3 (H3K27me3) normally represses genes expressions. Peripheral blood lymphocytes are the main source of somatic cells in the milk of dairy cows that vary frequently in response to the infection or injury of mammary gland and number of parities. Methods The genome-wide status of H3K27me3 modifications on blood lymphocytes in lactating Holsteins was performed via ChIP-Seq approach. Combined with digital gene expression (DGE) technique, the regulation effects of H3K27me3 on genes expressions were analyzed. Results The ChIP-seq results showed that the peaks of H3K27me3 in cows lymphocytes were mainly enriched in the regions of up20K (∼50%), down20K (∼30%) and intron (∼28%) of the genes. Only ∼3% peaks were enriched in exon regions. Moreover, the highest H3K27me3 modification levels were mainly around the 2 Kb upstream of transcriptional start sites (TSS) of the genes. Using conjoint analysis with DGE data, we found that H3K27me3 marks tended to repress target genes expressions throughout whole gene regions especially acting on the promoter region. A total of 53 differential expressed genes were detected in third parity cows compared to first parity, and the 25 down-regulated genes (PSEN2 etc.) were negatively correlated with H3K27me3 levels on up2Kb to up1Kb of the genes, while the up-regulated genes were not showed in this relationship. Conclusions The first blueprint of bovine H3K27me3 marks that mediates gene silencing was generated. H3K27me3 plays its repressed role mainly in the regulatory region in bovine lymphocytes. The up2Kb to up1Kb region of the down-regulated genes in third parity cows could be potential target of H3K27me3 regulation. Further studies are warranted to understand the regulation mechanisms of H3K27me3 on somatic cell count increases and milk losses in latter parities of cows.


Introduction
Dairy cows provide a mass of milk and protein for the human beings world-wide. High quality milk is required for human health. For dairy cattle, ten or more lactations are possible and the highest milk production of dairy cows should be at the fifth parity from physiological perspective. However, the average herd life of Holsteins is under three lactations at present [1]. Studies reported that over 90% of all cows are culled due to reasons of infertility, mastitis, lameness and lower milk production [2,3,4]. Although somatic cells are normal constituent of milk, high somatic cell counts (SCC) in the milk have been shown to be significantly related with milk quality, milk losses and mastitis risk [5,6].
Peripheral blood lymphocytes are the main source of somatic cells in the milk of dairy cows in response to the infection or injury of the mammary gland [7]. Lymphocytes include three major types of cells, natural killer cells, T cells and B cells. Natural killer cells are a part of the innate immune system, whereas, T cells and B cells are the major components of the adaptive immune response [8]. Proper lymphocyte development and function are the result of multiple cell lineage choices and cellular transitions with the hematopoietic system. Numerous studies are working on understanding the contribution of transcriptional and epigenetic regulation during lymphopoiesis at both physiological and pathological levels in recent years [9,10,11,12,13]. The complex relationships between transcription level and the epigenetic machinery of lymphocytes identity and development are partially revealed [12,14,15]. Posttranslational histone modifications are the major part of epigenetics, which include acetylation, methylation, phosphorylation and ubiquitination [16]. Acetylation is the most frequently studied histone modification and is commonly localized to active chromatin. However, histone methylation plays a dual role as it is important in both activation and repression, depending on where the methyl group is modified. Trimethylation of lysine four on histone H3 (H3K4me3) is associated with transcriptional activity [17]. By contrast, trimethylaiton of lysine 27 on histone H3 (H3K27me3) is related with inactivation of genes expressions or located in heterochromatin [18,19]. Studies reported that H3K27me3 are relatively stable through lymphocyte divisions and for preserving cell identity, and are therefore considered as potential marks for transmitting the epigenetic information [20,21].
In bovine, the remodeling of H3K27me3 was investigated during development of in vitro fertilization (IVF) and somatic cell nuclear transfer (SCNT) bovine embryos [22,23]. The global dynamics of H3K27me3 during bovine oocyte maturation and preimplantation development were detected by semiquantitative immunofluorescence [24]. However, genome-wide H3K27me3 modifications and regulatory function are still unclear in bovine lymphocytes. The objectives of the study were to map the genomewide H3K27me3 landscape of bovine via high-throughput chromatin immunoprecipitations coupled with sequencing (ChIP-Seq) approach, and to declare the regulatory effects of H3K27me3 modification on genes expressions in peripheral blood lymphocytes by means of digital gene expression (DGE) techniques. The potential functions of H3K27me3 modification on the parity-associated differentially expressed genes were also analyzed and discussed.

Genome-wide Maps of H3K27me3 Modifications in Cows Lymphocytes
To reveal epigenetic features of lymphocytes in bovine peripheral blood, we generated global maps of H3K27me3 modifications via the ChIP-seq approach. A total of 9.7 million short reads from each sample were aligned to the bovine reference genome (Btau4.0). We utilized the MACS model-based algorithm for identification of significant clustering of windows with enriched ChIP signals, which were called ''peaks'' of H3K27me3 in the study. The scaling analysis indicated that the sequence reads were sufficient to identify most regions associated with the H3K27me3 modifications in the bovine genome [25] (see Figure S1 for details). To understand the functional consequences of H3K27me3 modification, gene expression in bovine lymphocytes were analyzed by the DGE assay. Therefore, we made a combined analysis for ChIP-seq and DGE to gain insights into the H3K27me3 regulation pattern in bovine lymphocytes ( Figure 1).
Distribution of H3K27me3 reads. To obtain an overall picture of H3K27me3 distribution, all reads of ChIP-seq were aligned to the bovine genome and only uniquely matching reads were retained after filtering dirty reads. Mapping rate was 96% and unique mapping rate was 77%, which laid down a basis for further analysis (Figure 2A). To study H3K27me3 distribution in genome region, we divided the bovine genome into five kinds of regions-up20K (20 kb upstream of transcription start site (TSS)), exon, intron, down20K (20 kb downstream of transcription end site (TES)) and intergenic regions-on the basis of annotation of ''known genes'' from BioMart Query System of ENSEMBL Btau4.0 database developed jointly by the OICR and the EBI. The proportion of each region of the total genome was indicated ( Figure 2B). As shown in Figure 2B, more reads were identified in intergenic regions (,83%) and few reads in exon regions (,0.35%). However, their proportions of the total genome were different with intergenic of 78.7% and exon regions of 0.71%. Reads located in up20K and down20K were 8.4% and 6.7%, but their proportions of the total genome were the same (7.4%). For intron region, approximately reads of 11% were distributed and this region accounted for the entire genome of 13%. In view of this point, abundance of reads located in these regions was drawn in Figure 2C. It showed that the abundance of up20K was highest and that of exon region was lowest. To visualize the distribution trends of H3K27me3 in the genic regions, a composite profile of H3K27me3 for all known genes was generated, spanning their gene bodies and extending 20 kb upstream and 20 kb downstream ( Figure 2D). It is notable that H3K27me3 signals levels were high on upstream20Kb of TSS. Moreover, it appeared that H3K27me3 distribution decreased dramatically within gene bodies, in accordance with the observation that less of the H3K27me3 reads were located in exon regions ( Figure 2B). These results were consistent with previous observations that H3K27me3 associates extensively with proximal promoters in human T cells, as well as human and murine embryonic stem cells [25,26,27].
Distribution of H3K27me3 peaks. As the final ChIP DNA fragments contain immunoprecipitated target DNA only at minor fraction (reportedly ,1%) within a large proportion of background input DNA [28,29], many researchers seek to identify locations of ChIP target DNA in which sequence reads cluster into many enrichment regions, wherein overlapping reads appear as peaks. In the present study, MACS1.4.0 software was used to identify enriched ChIP regions of H3K27me3. In total, 831, 607, 982 and 703 peaks were found, while the known genes related to peaks were 285, 195, 311 and 165 in the four cows, respectively ( Figure 3A). To study the pattern of epigenetic modifications in different regions of genes, we also calculated distribution percentages of peaks identified in four kinds of genic regions. The result showed that most peaks enriched in gene regions of up20K (,50%), down20K (,30%) and intron (,28%), whereas only 3% peaks were enriched in exon regions ( Figure 3B), suggesting that H3K27me3 modification in bovine lymphocytes may play a major role in regulatory region of genes. It is worth mentioning that a higher percentage of enrichment in up20K of genes was marked with H3K27me3 modification compared to in down20K region, consistent with the observation that more H3K27me3 enrichments were located 5 Kb upstream of TSS than in 5 Kb downstream of TES in murine T cell subset [25]. In summary, our data indicated that high H3K27me3 modification levels of bovine lymphocytes were mainly in regulatory regions of genes.

Correlation of H3K27me3 Modifications with Gene Expression
To reveal the functional consequences of H3K27me3 on target genes, we generated expression profiles for the four individuals' lymphocytes using DGE tag profiling technology, which uses Illumina high-throughput sequencing to generate a large number of tags of length 21 bp. For the four individuals, 3.66610 6 , 3.52610 6 , 3.55610 6 and 3.64610 6 raw tags were generated, and the numbers of unambiguous tag-mapped genes were 6862, 6894, 6764 and 6923, respectively.
Correlation of H3K27me3 modifications near TSS with gene expression. Gene promoters near TSS contain critical regulatory elements necessary for transcription. Above genome-wide analyses have revealed that H3K27me3 modification was enriched in the promoter regions of the target genes. To correlate this modification with gene transcription, the bovine genes whose expression levels were determined by the DGE assay were separated into four groups of around 700 genes in each group according to their expression levels. The H3K27me3 tag numbers in each region were calculated and normalized near the TSS for four sets genes corresponding to highly expressed, two types of intermediately expressed (medium and low) and silent genes ( Figure 4, Table S1). As expected, H3K27me3 signals were correlated with gene repression (Figure 4), which is in agreement with previous studies in human T cell and murine embryonic stem cells [30,31]. Obviously, H3K27me3 levels were higher at silent genes than at active genes. Intriguingly, H3K27me3 levels were elevated surrounding the TSSs and peaked at upstream 2 K of TSS for the silent genes sets (dotted line), though were not significant for the other three sets. In order to confirm it, 113 silent genes were selected randomly to analyze their H3K27me3 profiles in all four individuals. The result showed that H3K27me3 levels indeed rise surrounding of the upsteam2K of TSS ( Figure S5D), suggesting that H3K27me3 repressed gene expression by acting on the surrounding of upstream2K. A significant dip in the signal was observed at the TSS (Figure 4), which may correlate with lower nucleosome density at this region [32]. Moreover, it is evident that the modification levels of up10K relative to TSS were higher than of downstream10K regions, which conform to the result of Figure 2C and suggest that H3K27me3 occurs mainly in gene promoter region.
Correlation of H3K27me3 modifications across transcribed regions with gene expression. To reveal H3K27me3 features in transcribed regions of genes, we plotted their average modification levels in their transcribed regions and extending 20 kb upstream and 20 kb downstream for four sets genes. The combined analysis with H3K27me3 levels and genes expression profile showed that H3K27me3 marks tended to repress expressions of target genes across the whole gene regions ( Figure 5, Table S2). Moreover, H3K27me3 were mainly enriched in upstream and downstream regions of the target genes, which is consistent with our above results and similar to H3K27me3 modifications in human T cell [20].
H3K27me3 modifications characteristics in cytokines and transcription factors genes. Different effectors of lymphocytes are characterized by lineage-specific expression of cytokine genes, and transcription factors play a role in this lineage commitment [33]. To further characterize the regulatory effects of H3K27me3 on the marked genes, we combined H3K27me3 maps with DGE data to identify groups of genes associated with cytokines and transcription factors. Two hierarchical clusters were generated: the first cluster contained all cytokines and their related genes (total 138 genes, see Table S3 for details; Figure 6A); the second cluster had all transcription factors and their related genes (total 178 genes, see Table S4 for details; Figure 6B). These two clusters revealed that most genes had high expression but low H3K27me3 levels, conversely, a small number of genes had very low expression and their H3K27me3 modifications were high. The results demonstrated that H3K27me3 was associated with inactive genes acting on the promoter regions. Notably, 15.2% cytokine genes and 6.2% TF genes were not negatively regulated by their H3K27me3 modifications (Table 1). For example, high amounts of CD40, NLRP3 and TFEC mRNA showed high levels of H3K27me3 modifications at their promoters, while IL17RC and MYT1L had low mRNA expressions and low H3K27me3 modifications levels. Therefore, H3K27me3 regulates genes expressions but could not permanently repress their expressions [25].

H3K27me3 Modifications Effects on the Parity-specific Differentially Expressed Genes
In order to reveal whether the bovine parity-related genes correlate with H3K27me3 modifications, differentially expressed genes were screened and analyzed between the first and the third parity cows. A total of 53 differentially expressed genes were found (FDR #0.01 and | log2Ratio |$1), in which 25 genes were high expressed in the first parity and low expressed in the third parity ( Figure 7A), and the other 28 genes were conversely ( Figure 7C). To explore the relationship between H3K27me3 modification and bovine parity-related genes, we calculated H3K27me3 levels in the regions of up5K to TSS of the parity-specific differentially expressed genes. As shown in Figure 7B, the H3K27me3 modification levels of the thirdparity-specific down regulated genes were significantly higher in the third parity cows than of the first parity cows, and the remarkable modification region was only located in the up2K-up1K region relative to TSS (P,0.05, t-test) ( Figure 7B).  The mapping result of H3K27me3 reads. Raw reads were generated by Solexa sequencing, and then clean reads were obtained after filtering dirty reads. All clean reads were mapped to the bovine reference genome and only uniquely matching reads were retained to use for subsequent analysis. Mapping rate is the ratio of mapped reads to clean reads and unique mapping rate is the ratio of unique mapped reads to clean reads. (B) Distribution of H3K27me3 reads among different genomic regions. The bovine genome was divided into five kinds of regions: 20 kb upstream of transcription start site (TSS), exon, intron, 20 kb downstream of transcription end site (TES) and intergenic regions. The histogram described the percentage of unique mapped reads among five genomic regions and the proportion of each region of the total genome. (C) Abundance of H3K27me3 reads among different genomic regions. The percentages of reads distribution were normalized to the abundance values. (D) Coverage depth of H3K27me3 reads among genic regions. For each gene, the tag numbers detected in every 5% of the gene-body region and every 1 kb outside of the gene-body region were summed to obtain methylation levels. These numbers were then normalized by the total number of base pairs in each region [20]. doi:10.1371/journal.pone.0039094.g002  However, the up-regulated genes in the third parity cows did not showed any relationship with H3K27me3 level ( Figure 7D).
To explore the function significance of DGE differentially expressed genes regulated by H3K27me3 modification in the third parity cows, we further conducted GO analyses in DAVID database for the 53 differentially expressed genes. The result indicated that HEXIM2, PSEN2 and PTX3, PSAP genes of downregulated in the 3rd parity cows were involved in the GO term of transcription regulator, immune response and reproductive process, respectively ( Figure S6). Moreover, the pathway analysis for the 25 down-regulated genes was also performed. The result showed that down-regulated genes mainly participated in cancerand disease-related pathways ( Figure S7).

Real-time PCR Validation
To assess the accuracy of the ChIP-seq mapping result, bovine CD4 and IL10 cytokine genes were used to confirm their H3K27me3 enrichment profile using ChIP -quantitative PCR (ChIP-qPCR) approach [34]. We arbitrarily chose five enriched regions from the H3K27me3 maps for the two genes (4 sites on the gene body region of CD4 and 1 site on the promoter region of IL10). Relative enrichment was quantified for each site with real-time PCR reactions using 0.5 ng ChIP DNA or 0.5 ng input DNA for each sample and normalized by the negative controls [34]. For the four sites on CD4 gene ( Figure 8A), the relative enrichments was mostly consistent with the profiles observed in ChIP-seq ( Figure 8B and 8C), so did for IL10 gene ( Figure 8E and 8F).
In order to confirm the results of DGE and explore the relationship between H3K27me3 enrichment and gene expression, the expression levels of CD4 and IL10 genes were also detected with reverse transcription -quantitative PCR (RT-qPCR) and standardized with three housekeeping genes. The results showed that the expressions of CD4 and IL10 were mostly consistent with the data in the DGE ( Figure 8D and 8G, except CD4 expression level in C2 individual in DGE) and were negatively related with the level of H3K27me3. Consequently, the ChIP-seq and DGE data were reliably and efficiently used for the analysis.

Discussion
Genome-wide measurements for studying the interactions between histone and DNA as well as their regulation on transcriptional activity were increased rapidly by high-throughput DNA sequencing techniques [35,36]. In the present study, we first described the genome-wide H3K27me3 profiling of peripheral blood lymphocytes in dairy cows. The results indicated that H3K27me3 modification were associated with gene silencing in bovine lymphocytes. Five H3K27me3 modified   regions were identified, which included enrichments in the promoter, at the transcriptional start site, at the transcriptional end site, across the gene body region and 3' region. All regions present in the four Holstein cows depicting that H3K27me3 were enriched mainly in upstream and downstream regulation regions of genes. It was notable that differences of H3K27me3 were observed between the first and the third parity cows. The results provided new insights into gene expression regulated by H3K27me3 modification in bovine peripheral blood lymphocytes.
Lymphocytes are a key part of the immune system and involved in innate and adaptive immunity. The bovine genome-wide patterns of H3K27me3 in peripheral blood lymphocytes are in agreement with previous findings of CD4+ T cells in human [20], and also provided novel information on lymphocytes epigenetic regulation. As previous studies reported, higher H3K27me3 signals in bovine lymphocytes were detected at silent promoters but decreased levels were present at active promoters. In particular, the upstream2Kb region relative to the TSS showed significant correlation with gene silencing, while the upstream region from 2 Kb to 10 Kb relative to TSS of the silent genes as well as the highly, moderately and lowly expressed genes did not present this feature ( Figure S5A, S5B, S5C and S5D). The possible reason was that the region of up2K-up10K of the genes might be overlapped with intergenic regions that did not show much modification changes of H3K27me3 [20]. The results indicated that the H3K27me3 modification on 2 Kb upstream of the TSS have a significant regulatory effect on genes suppression in bovine lymphocytes.
Somatic cell count (SCC) is an indicator of milk quality and udder health in dairy herds. General agreement relies on the values of SCCs that less than 100,000 cells/ml for healthy cows, when SCC greater than 300,000/ml is indicated a problem with subclinical or clinical mastitis and inferior milk quality [6]. It was found that SCC increases with increased parity that could be due to increased risk of infectious pathogens entering the  udder [37,38,39]. Considering the SCC of the third-parity cows (SCC.400,000 cells/ml) significantly higher than the first-parity cows (SCC ,150,000 cells/ml) [40], we compared the differences of H3K27me3 modification and gene expression levels between the third-parity cows and the first-parity cows. As expected, the 25 decreased expression genes in the third parity cows were suppressed by the H3K27me3 modification, which were mainly located in the up2K-up1K region relative to TSS of the genes. Among the genes, PSEN2 (presenilin-2) is extremely correlated with Alzheimer's disease (AD) risk. Presenilin has been expressed all through the mammal development process and the loss of its function causes neuron apoptosis and impairs long term memory and cognitive deficits in AD patients. For PSEN1 and 2 double knock-out mice, loss of presenilins function leads to up-regulation of inflammatory markers in the cerebral cortex [41]. Thus the potential function of silenced PSEN2 on increased SCC companion with the thirdparity cows is worth further study. Unexpected, the down expressed genes involved in immune response in the first parity (but up expressed in the third parity, such as IFI6 and IFI47) cannot be found the regulation effect of H3K27me3 ( Figure 7C and 7D), which might be regulated by other epigenetic regulatory machinery such as H3K4me3.
The GO analysis of the differentially expressed genes between the two parities demonstrated that unique third-parity enrichment genes mainly involve in protein-DNA complex and some important immune processes. In the above analysis with cytokines and transcription factors cluster, several genes appeared as parity specific. We assumed that third-parity cows may be infected by subclinical mastitis by the following mechanisms: due to low level of H3K27me3 modification of IL10, more IL10 released by cytotoxic T-cells could inhibit the actions of NK cells during the immune response to bacterial infection [42]. These results suggested H3K27me3 levels of vital genes related to immune response could serve as a prognostic factor for udder health of dairy cattle [43].
We applied direct sequencing of ChIP DNA using Illumina Hiseq2000, which is an efficient method for mapping genomewide distributions of histone modifications and chromatin protein target sites. The technique is particularly well suited for analyzing protein locations and chromatin modification in large genomes, and allows researchers to survey more of the genome in less time. Indeed, we used these techniques to provide the most comprehensive genome-wide data for H3K27me3 marks in bovine lymphocytes. However, in the study, hundreds of peaks in bovine lymphocytes were found, which is far less than that in the human genome [44] and also in the maize genome [45] (see Figure S4). The reason may be associated with species, breeds or samples. As the H3K27me3 modification is highly complex, therefore, further studies are required to determine its function in different kinds of cells and tissues using larger sample sizes.
To the best of our knowledge, we have generated the first high resolution map of histone methylation for bovine lymphocytes. The results not only provides an initial blueprint of bovine H3K27me3 as an epigenetic mark that mediates gene silencing, but also offers a novel view of the complexity of regulating gene expression of lymphocytes. The down-expressed genes in the third parity cows regulated by H3K27me3 modification suggest that the regions of the up2Kb to up1Kb of the TSS are a potent target of epigenetic regulation. Further studies are warranted to better understand the regulation mechanism of the histone methyaltion on SCC increases or milk losses in the latter parity of dairy cattle.

Cows and Isolation of Lymphocytes
Four Holstein cows (C1, C2, C3 and C4) were randomly selected from a dairy herd in Beijing, China. They were fed on the same lactation diet according to energy recommendations for lactating Chinese Holstein cows and were handled in accordance with the guidelines of the Animal Care and User Committee. Blood samples were collected from the jugular vein for each animal. Performance testing data were provided by the official Dairy Data Center of China (Beijing, China), including daily milk yields, fat percentage, protein percentage and somatic cell counts (SCC) were obtained (Figure1B). In addition, peripheral blood lymphocytes were prepared by Lymphocytes Separation Medium (TBDsciences, Tianjin, China, PN.LTS1086) according to the manufacturer's instructions, and the purity was 90%-95%.

ChIP and ChIP-Seq
Lymphocytes from peripheral blood were digested with MNase (TAKARA, D2910) to generate mainly mononucleosomes with a minor fraction of dinucleosomes for histone modification mapping. Chromatin from 5610 7 cells were used for ChIP experiment and antibodies against histone H3K27me3 (millipore, Cat. #17-622) were used, that yielded approximately 30 ng of DNA. The enrichment efficiency of ChIP was detected using a standard Real-time Quantitative PCR approach with 1 ng ChIP DNA and non-enriched input DNA (control samples) as template. The ChIP DNA fragments were blunt-ended and ligated to the Solexa adaptors (Paired-End DNA Sample Prep kit, Illumina, Cat. #PE-102-1001), the ChIP DNA was then amplified using the adaptor primers for 15 cycles and the fragments around 100-300 bp (mononucleosome + adaptors) were isolated from agarose gel. The purified DNA was used directly for cluster generation and sequencing analysis using the Illumina Hiseq2000 following manufacturer protocols (BGI, China).

Solexa Sequencing Data Analysis of ChIP-seq
Sequence reads of 49 bp length of ChIP-seq were obtained using the Solexa Analysis Pipeline that were mapped to the bovine reference genome (Btau4.0) by SOAP 2.21 software and only uniquely matching reads were retained after filtering dirty reads. Unique reads numbers for each library were listed in Figure 2(A).
The output of the Solexa Analysis Pipeline was converted to browser-extensible data (BED) files for viewing the data in the UCSC genome browser. Data for H3K27me3 modifications were presented in both BED and graph format. BED file indicates the genomic coordinates of each tag. To make the graph file, we mapped tags into non-overlapping 50 bp windows of reference genome. The location of a tag on positive or negative strand was determined and represented with the beginning of the tag relative to reference genome. Based on these locations, the tags overlapping each 50 bp summary window were counted.
For ChIP-seq experiments, peak areas represent in vivo locations where proteins of interest (e.g. modified histones or transcription factors) were associated with DNA. We identified such peaks using the ''model-based analysis of ChIP-seq data'' MACS1.4.0 software, and its function is to isolate ChIP-enriched regions from non-enriched regions based on a dynamic Poisson distribution model. Detailed algorithms and models were described by Zhang et al [46]. We set up a bandwidth of 200 bp, mfold of 30, p-value of 1.00e-05 under a FDR cutoff of 1% to call peaks representing enriched H3K27me3 marks. The output result includes one EXCEL file containing the genome coordinates, summit, p-value, and fold_enrichment of each peak. Furthermore, for each gene, as long as there is 1 bp overlap between regions of a peak and a particular gene (includes regions of up20K, exon, intron and down20K) we consider that this peak is associated with this gene.

Digital Gene Expression and Analysis
Total RNA from primary lymphocytes was extracted with mirVana miRNA Kit (Ambion, PN.AM1560) according to the manufacturer's instructions. Approximately 6 mg of total RNA was transcribed into double-stranded cDNA through Reverse Transcription Kit (Applied Biosystems, PN.4368814). Sample preparation for digital gene expression profiles was achieved through the Gene Expression Sample Prep Kit (Illumina, Part #1004241) with details as follows. Restriction enzyme NlaIII was used to cut the CATG sites of cDNA strand, and Illumina adaptor A was ligated.
Mmel was used to digest at 17 bp downstream of the CATG site and Illumina adaptor B was ligated at 3' end. Primer GX1 and Primer GX2 were added for PCR for 15 cycles. Then, 95 bp fragments were regained through 6% TBE PAGE (see Figure S2). The DNA was purified followed by sequencing with Illumina Hiseq2000 Analyzer. Consequently, each tunnel generated millions of raw reads with sequencing length of 35 bp, but tags that include CATG site were of 21 bp length.
To determine the position of tags sequences, a more stringent mapping method (SOAP 2.21 software) were adopted [47,48] to map all tags to the bovine reference genome. Raw sequences were transformed into clean tags after certain steps of data-processing. All clean tags were mapped to the reference sequences and only tags with perfect matching or 1 bp mismatch were considered.
The expression level of one gene was represented by the total number of tags that uniquely aligned to this gene. In order to compare multiple samples differences for one gene, the number of unambiguous and clean tags for each gene was calculated and then normalized to TPM (number of transcript copies in per million clean tags), equaling to copy number of clean tags for this gene divided by total number of clean tags and multiplied by one million [49,50]. Sequencing quality evaluation is shown in Figure  S3.

Combination Analysis of ChIP-seq with DGE Data
To study the correlation of H3K27me3 modifications with gene expressions, transcriptional levels of genes in bovine lymphocytes were obtained by DGE analysis. These genes (,7000 genes) were then broken up into 10 sets of 700 genes by ranking these genes expression levels and assigning them to sets of 700 genes. Four out of the ten sets shown in Figures 4 and Figure 5 correspond to highly expressed, two degrees of intermediately expressed (medium and low) and silent genes. Tags detected were aligned in each gene set across transcription start sites (TSS) or gene bodies. For H3K27me3 modification near TSS, the tag density (number of tags per base pair) was calculated in 50 bp windows relative to the TSS. To calculate the H3K27me3 profiles across the gene bodies, the tag numbers detected in every 5% of the gene-body region and every 1 kb outside of the gene-body region were summed and normalized in the four expressed sets.

Quantitative Real-time PCR
ChIP-q-PCR reactions with SYBR green dye were carried out using Roche LightCycler 480 qPCR machine to confirm the enrichment of selective regions. PCR primer pairs were designed using Primer3 (http://fokker.wi.mit.edu/primer3/input.htm) and confirmed by Oligo 6.0. Primer sequences are given in Table S5. A standard curve was performed using genomic DNA prepared for each primer pair to determinate primer pair efficiency [51], and three serial tenfold dilutions were used ranging from 0.5 ng to 50 ng. Each 20 mL reaction contained 7 mL of nuclease-free water (Gibco), 10 mL of LightCycler 480 SYBR Green I Master mix (Roche, Cat. #04 707 516001), 1 mL of each primer (10 pmol/ mL) and 1 mL of DNA. The ChIP-qPCR reactions were done in duplicates for each site. To determine the relative fold enrichments, the 2 2ggCp method was used by comparing enrichment values for a given primer pair (totally five pairs in CD4 and IL10) to a negative primer pair (GAPDH_P1 or 18s rRNA_P1) between experimental (ChIP DNA) and reference (input DNA) samples. The P1 sites in promoters of GAPDH and 18s rRNA were selected as negative controls that has little H3K27me3 modification and no difference among the experimental samples and control samples ( Figure S8) [34].
For RT-qPCR of gene expression, 3 mg of total RNA were reverse transcribed into cDNA with High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems, PN.4368814) according to the recommended procedure. Triplicates were performed for RT-qPCR reactions. qPCR reaction was run using the program as follows: pre-incubation (95uC for 10 min), 45 cycles of amplification (95uC for 10 s, 60uC for 10 s, and 72uC for 10s), melting curves using a heat ramp and cool down. Crossing point values (Cp values) were obtained using the second derivative maximal analysis tool included in the Roche LightCycler480 software. The mRNA expression of CD4 and IL10 was normalized against three housekeeping genes (GAPDH, 18s rRNA and beta actin) cDNA in the corresponding samples. Figure S1 Scaling analysis of H3K27me3 peaks. The tag density on per peak of the sample was plotted. By increasing the fraction of tags selected for peak identification, a degree of saturation was estimated based on number of tags per base pair on the peak. (DOCX) Figure S2 Principle and procedure of Tag preparation for DGE pipeline [49]. Beads of Oligo(dT) are used to enrich mRNA from the total RNA, and then are synthesized to the first and second-strand cDNA by use of Oligo(dT) as primer. Subsequently restriction enzyme Nla? recognizes and cut off the CATG sites of cDNA strand, and the Illumina adaptor A is ligated to the sticky 5' end. In the same, Mme? which is a type of Endonuclease cut at 17 bp downstream of the CATG site, and the Illumina adaptor B is ligated to the 3' ends of tags, consequently, a tag library that tags with different adaptors of both ends is acquired. After 15 cycles of linear PCR amplification, 95 bp fragments are purified by 6% TBE PAGE Gel electrophoresis. The purified fragments were used directly for cluster generation and sequencing analysis using the Illumina Hiseq2000 following manufacturer protocols (BGI, China). (DOCX) Figure S3 Distribution of Tag Expression in DGE data. The upper panel is the total tag number of sample C1, C2, C3 and C4. For example, ''Tags Containing N (209886, 3.75%)'' means the number of tags containing N is 209886 and 3.75% of the total tags. The under panel is the distinct tags number of sample C1, C2, C3 and C4. For example, ''Tags Containing N (92427, 9.61%)'' means the number of the distinct tags containing N is 92427 and 9.61% of the total distinct tags; ''Only adaptors'' means the reads contain only the adaptors sequence; ''Copy Number ,2'' is the tags whose copy number is less than 2; ''Clean tags'' is the tags used to analysis after filtering the dirty tags. Raw sequences have 3' adaptor fragments as well as a few low-quality sequences and several types of impurities. Raw sequences are transformed into Clean Tags after certain steps of data-processing. (DOCX) Figure S4 H3K27me3 Peaks information. (A) to (C) Peaks Numbers, average lengths, and total lengths of epigenetically modified regions detected by MACS1.4.0 software, respectively. (DOCX) Figure S5 H3K27me3 modification in four different expressed sets. Profiles of the H3K27me3 covered the region of upstream 10 K to TSS for highly active (A), two kinds of intermediately active (medium (B) and low (C)) and silent gene (D) sets were shown. Each gene set included common genes shared by four individuals (C1, C2, C3 and C4), which were screened from 700 genes in Figure 4 or Figure 5. Here, the tag density (number of tags per base pair) was calculated in 50 bp windows in upstream 10 K regions to TSS (see Experimental Procedures). (DOCX) Figure S6 Annotation of DGE differential genes for two parities by WEGO. Gene Ontology Annotation Plotting. The BGI WEGO (Web Gene Ontology Annotation Plotting) was used to functionally categorize parity-specific differentially expressed genes. Of 53 differentially expressed genes, 43 genes with GO annotation that belong to two parities were grouped by cell component, molecular function and biological process based on the bovine GO annotation information (http://www. geneontology.org/GO.downloads.annotations.shtml). Gene numbers and percentages (on log scale) are listed for each category. (DOCX) Figure S7 The pathway analysis of down-regulated genes in the third parity. KEGG pathway analysis in DAVID for the 25 down-regulated genes was completed. Y axis indicated pathway terms of involving in these down-regulated genes.