Genome-wide bovine H3K27me3 modifications and the regulatory effects on genes expressions in peripheral blood lymphocytes

Abstract

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.

Figure 1. Overview of the experimental design and samples characterization.

(A) Peripheral blood was collected from cow caudal vein, and primary lymphocytes were isolated. The genome-wide distribution of H3K27me3 was determined with ChIP-seq, and all genes expressions in lymphocytes were performed with Digital Gene Expression technique. (B) Performance testing information about the four cows. ^†Days_in_Milk is the days from calving to sampling. ^‡Daily milk yield, fat percentage and protein percentage were the average values of annual records. ^§Somatic cell counts.

Figure 2. Distribution of bovine H3K27me3 reads.

(A) 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 .

Figure 3. Distribution of bovine H3K27me3 peaks.

(A) H3K27me3 peaks and genes related to peaks. (B) Percentages of genes related to H3K27me3 peaks. The bovine genic region was divided into four kinds of regions: 20 kb upstream of TSS, exon, intron, and 20 kb downstream of TES. The bar chart described the percentage of genes associated with H3K27me3 peaks among kinds of genic regions.

Figure 4. Modifications of H3K27me3 near transcription start sites.

Profiles of the H3K27me3 indicated across the TSS for highly active (high), two kinds of intermediately active (medium and low) and silent gene sets were shown. Each gene set included 700 genes according to their expression levels in primary lymphocytes of cow peripheral blood. Here, up and down 20 kb regions of 700 genes in each group were aligned relative to their TSSs (x axis). The y axis shows the detected tag density.

Figure 5. Modifications of H3K27me3 across the gene bodies.

Profiles of H3K27me3 patterns of the four sets genes were shown across the gene bodies. The classification of four groups is the same as Figure 4 and density plots extend 20 kb of 5’ and 3’ of the gene bodies. 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.

Figure 6. Hierarchical clustering of all cytokines and transcription factors genes.

(A) Based on annotation for bovine DGE profiles, we identified 138 unique cytokine genes and their related genes. Normalized intensity values of genes (rows) and their H3K27me3 modifications were ordered using Centroid Spearman Rank Correlation and hierarchical clustering in Cluster3.0 software. The dendrogram showed the similarity (distance) of mRNA expression levels and H3K27me3 modifications of genes and was divided into sub-trees as distinguished from different colors. Arrays (columns) were grouped by four different individuals. Yellow and blue colors reflect the high and low intensities, respectively. (B) Based on annotation for bovine DGE profiles, 178 unique transcription factors and related genes were identified. The heatmap and the hierarchical cluster were generated as above described.

Figure 7. The genes expressions and H3K27me3 modification of DGE differential genes between the two parity cows.

A total of 53 differentially expressed genes between the two parity cows were tested (FDR ≤0.01 and | log2Ratio |≥1), in which 28 genes were up-regulated and other 25 genes were down-regulated in the third parity cows. (A) The heatmap of 25 down-expressed genes for the four individuals. Normalized intensity values of genes (rows) were ordered using Centroid Spearman Rank Correlation and hierarchical clustering in Cluster3.0 software. The dendrogram showed the similarity (distance) of mRNA expression levels and was divided into sub-trees as distinguished from different colors. Arrays (columns) were grouped by four different individuals. Yellow and blue colors reflect the high and low expression intensities, respectively. (B) Modifications of H3K27me3 in 25 down-regulated genes for the four individuals. Profiles of the H3K27me3 covered the region of up5K to TSS of the genes were shown. The tag density (number of tags per base pair) was calculated in 500 bp windows in upstream 5 K regions to TSS. (C) The heatmap of 28 up-expressed genes for the four individuals. (D) Modifications of H3K27me3 in 28 up-expressed genes for the four samples.

Figure 8. H3K27 trimethylations of CD4 and IL10 overlaying different regions and genes expressions in bovine lymphocytes.

(A–D) The H3K27me3 enrichment and expression for bovine CD4 gene were showed in four individuals (C1, C2, C3 and C4). (A) The histone modification profile was shown in custom track in IGV2.0 (Integrative Genomics Viewer). Regions enriched for H3K27me3 in bovine lymphocytes were shown (blue bar). Four sites of CD4 validated by real-time PCR were indicated (line boxes). The position of the gene was presented on the bottom of the panel, where blue boxes represent exons and the arrow means the transcriptional direction of the gene. (B–C) Real-time PCR results showing enrichment of indicated four sites of CD4 in H3K27me3 ChIP-seq results carried out in the four cows. The negative controls were P1 site in the promoter of GAPDH (Figure 8B) and P1 site in the promoter of 18s rRNA (Figure 8C), respectively. (D) Real-time RT-PCR were performed on the samples of the four cows for validation of CD4 expression (GAPDH, 18s rRNA and beta actin were housekeeping genes). The individual whose Cp value corrected was the minimum was selected as a control sample to calculate relative expression of CD4 for all four cows. Error bars indicate standard deviation of three technical replicates. DGE indicates the results of digital gene expression of CD4 (right Y axis for TPM). (E–G) The H3K27me3 enrichment and expression for bovine IL10 gene were showed in the four individuals. The details are similar to CD4 gene. A zero-height bar indicates this gene was not expressed in individual C3.