Genome-wide identification of histone methylation (H3K9me2) and acetylation (H4K12ac) marks in two ecotypes of switchgrass (Panicum virgatum L.)

Background Histone modifications play a significant role in the regulation of transcription and various biological processes, such as development and regeneration. Though a few genomic (including DNA methylation patterns) and transcriptomic studies are currently available in switchgrass, the genome-wide distribution of histone modifications has not yet been studied to help elucidate gene regulation and its application to switchgrass improvement. Results This study provides a comprehensive epigenomic analyses of two contrasting switchgrass ecotypes, lowland (AP13) and upland (VS16), by employing chromatin immunoprecipitation sequencing (ChIP-Seq) with two histone marks (suppressive- H3K9me2 and active- H4K12ac). In this study, most of the histone binding was in non-genic regions, and the highest enrichment was seen between 0 and 2 kb regions from the transcriptional start site (TSS). Considering the economic importance and potential of switchgrass as a bioenergy crop, we focused on genes, transcription factors (TFs), and pathways that were associated with C4-photosynthesis, biomass, biofuel production, biotic stresses, and abiotic stresses. Using quantitative real-time PCR (qPCR) the relative expression of five genes selected from the phenylpropanoid-monolignol pathway showed preferential binding of acetylation marks in AP13 rather than in VS16. Conclusions The genome-wide histone modifications reported here can be utilized in understanding the regulation of genes important in the phenylpropanoid–monolignol biosynthesis pathway, which in turn, may help understand the recalcitrance associated with conversion of biomass to biofuel, a major roadblock in utilizing lignocellulosic feedstocks. Electronic supplementary material The online version of this article (10.1186/s12864-019-6038-x) contains supplementary material, which is available to authorized users.


Background
In eukaryotes, gene expression is a complex, coordinated process orchestrated by transcription factors (TFs) and chromatin proteins. As a result, chromosomal DNA is either compacted or relaxed to repress or activate the transcription process, respectively. At the whole-genome level, genetic and epigenetic factors together regulate gene expression. The accessibility of DNA sequences and chromatin sites associated with nucleosome packaging is known to determine transcription efficiency. Epigenetic factors such as DNA methylation, histone modifications, and chromatin remodelors play a significant role in the regulation of gene expression and various biological processes, such as development [1], regeneration [2], and oncogenesis [3]. In plants, with recent advances, genome-wide epigenomic analysis methods are increasingly supporting the role of chromatin modifications in gene expression [4].
Transcription factors are important in gene regulation; they bind to cis-regulatory sequences in the promoter region of their target genes to fine-tune gene expression [5]. In plants, transcription factors and histone modifications demonstrably play a role in development, vernalization, and response to biotic and abiotic stressors. [6][7][8][9][10]. The application of high-throughput technologies such as chromatin immunoprecipitation followed by genome-wide sequencing (ChIP-Seq) provides an opportunity to identify binding patterns for a protein(s) of interest in the entire genome. Genome-wide association and correlation of modifications, such as histone methylation and acetylation with biotic and abiotic stress responses in plants have been reported [8,[11][12][13].
Switchgrass (Panicum virgatum L.) is a native perennial and warm-season grass, widely used for biofuel, livestock feed, erosion control, and wildlife habitat. Carbon sequestration by switchgrass, a C4 plant with high cellulosic conversion rates and drought tolerance, facilitates its use as a biofuel directly or by ethanolic conversion, making it a valuable, natural renewable energy resource [14]. To date, studies in switchgrass have pertained to molecular breeding [15,16], genomic [17,18] and transcriptomic [19,20] analyses, and response to biotic and abiotic [21,22] stresses. Using Bisulfite sequencing (BS-Seq) and Methylated DNA immunoprecipitation sequencing (MeDIP-Seq), similar DNA methylation patterns between AP13 and VS6 have been revealed, and the methylation levels were higher in Transposable elements (TEs) when compared to genic regions [23]. The recent bisulfite sequencing study suggested small interfering RNAs (siRNAs) positively regulate DNA methylation and thereby interferes with gene and long non-coding RNAs (lncRNAs) expression in switchgrass [24]. However, the genome-wide distribution of histone modifications in switchgrass has not been reported. The histone marks (a repressive mark, H3K9 me2 , and an activation mark, H4K12 ac ) utilized in this study have been studied and reported previously in several plant species [7,8,25,26]. Methylation of histone H3K9 is preferentially localized to heterochromatin in Arabidopsis [25], and its profiling has been reported in other plants, for instance, rice [7], maize [26], and recently in the common bean [8]. On the other hand, H4K12 ac is primarily localized in active coding regions of the genome and therefore creates binding sites for regulatory factors to promote transcription [8].
Switchgrass has two main ecotypes, lowland and upland, with distinct morphological and genetic characteristics and specific geographical niches [27]. Two important switchgrass genotypes (AP13 and VS16) were selected to represent each of these two contrasting ecotypes. Among these, AP13, a lowland ecotype, was used for the developing the switchgrass reference genome [28]. A female parent AP13 and male parent VS16 (an upland ecotype) were utilized for the construction of a genetic linkage map [18] and the identification of quantitative trait loci (QTL) associated with biomass yield and traits related to recalcitrance [22,29]. Due to their adaptability and ability to withstand selection pressure from diverse environmental stressors, the inherent genetic potential of these genotypes at the level of gene expression has been recently explored [30]. However, little is known about the epigenetic regulation of gene expression associated with these genotypes. The H3K9 me2 and H4K12 ac marks selected in this study regulate the closed and open chromatin states, respectively, during transcription. In this study, we report the first genomewide profiling of histone methylation H3K9 me2 and histone acetylation H4K12 ac marks by using ChIP-Seq in upland (VS16) and lowland (AP13) switchgrass ecotypes.

Plant materials and sample collection
The two contrasting, heterozygous, tetraploid genotypes (2n = 4x = 36) utilized in this study (AP13 and VS16) were derived from two ecotypes of switchgrass: the lowland cultivar, Alamo and the upland cultivar, Summer, respectively. The leaf materials were obtained from our collaborator and co-author on this manuscript Malay Saha at the Noble Research Institute (Ardmore, OK). Both AP13 and VS16 actively grow in summer, and hence the conditions (29/22°C day/night temperatures and a 16 h photoperiod) considered in the greenhouses of the Noble Research Institute were uniform. Fully expanded leaves were randomly collected and pooled from five plants of one-month-old clonal ramets, cryopreserved in liquid nitrogen, and stored at − 80°C. This experiment included two genotypes (AP13 and VS16), two marks (H3K9 me2 and H4K12 ac ), two conditions (input control and actual sample), and three replicates resulting in (2 X 2 X 2 X 3) a total of 24 samples. Of significance is that this current study on global ChIP-seq analyses utilized the same leaf samples that were collected, stored (− 80°C), and used for RNA-seq analyses in our previous studies [30]. Since our goal is to understand the genome-wide histone modifications in regulating global gene expression rather than a specific focus on pathways related to lignin formation, we have used leaves instead of stems.

Isolation and immunoprecipitation of chromatin
All the lyophilized leaf samples were ground to a fine powder in liquid nitrogen, collected in frozen 50 mL falcon tubes (Corning Inc., Corning, NY). Next, 25 ml (mL) of cold nuclear isolation buffer containing 1% formaldehyde, and 20 μl (μl) of proteinase inhibitor (Thermo Scientific, IL, USA, Product #87786) were added. The samples were incubated for 5 min at room temperature. Glycine (2 M) was used to terminate the crosslinking reaction. The lysate was filtered before nuclei pelletization at 3000 g for 10 min at 4°C. The pellet was resuspended in cold nuclear isolation buffer (no formaldehyde included). An equal volume of 15% Percoll solution was added to the lysate, and the nuclei were pelleted again at 3000 g for 5 min at 4°C. The pellet was then mixed in nuclear lysis buffer. This lysate was then sonicated using a Soniprep 150 (MSE UK Ltd., London, UK) for 6 pulses with 15 s between each pulse. The sheared DNA fragments collected after sonication ranged between 100 and 350 bp when visualized by agarose gel electrophoresis. ChIP dilution buffer was then added to the nuclear lysate to attain 10X dilution, and chromatin samples were then divided equally into two tubes. One of the samples was incubated with rabbit polyclonal antibody in a 20:1 mixture that specifically adhered to the N-terminus of histone subunit H4 acetylated at lysine 12 (H4K12 ac ) [31]. Also, the chromatin sample was incubated with a mouse antibody in a 20:1 mixture, which adhered precisely to the epitope region of histone H3 di-methylated at lysine 9 (H3K9 me2 ) [32]. Twenty microliters of Pierce Protein A/G magnetic beads (ThermoFisher Scientific, Grand Island, NY) were added to the nuclear lysate. The samples were then incubated overnight at 4°C. The chromatin samples that had no antibodies served as positive controls. Later, the centrifuged and pelleted magnetic beads were washed thrice with 100 μl ChIP dilution buffer. The magnetic beads bound to chromatin and antibody complex were then eluted with freshly prepared elution buffer. To terminate the crosslinking, 20 μl of 5 M NaCl was added and subsequently incubated at 65°C overnight. After this, 5 μl of 2 M Proteinase-K (Thermo Scientific) was added and incubated for 2 h (h) at 45°C. DNA was then collected by phenol-chloroform extraction followed by ethanol precipitation. The ChIP-DNA was resuspended in 20 μl of TE buffer. After verification by ChIP-PCR, each of these purified DNA samples were then used to generate a ChIP-Seq library using Illumina HiSeqTM 2500 (Illumina Inc., San Diego, CA). Libraries were sequenced at the Delaware Biotechnology Institute (DBI, Newark, DE).

Library construction and sequencing
The quality, purity, and size of the immunoprecipitated DNA samples were determined by an AATI Fragment Analyzer (Advanced Analytical, Ames, IA). The ChIP-Seq libraries (101 bp paired-end) were prepared by using Illumina TruSeq ChIP Sample Preparation kit (IP-202-1012; Illumina Inc., San Diego, CA). The NGS data generated in this study was submitted to the SRA section of NCBI with the project number PRJNA373877 (https:// www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA373877).

ChIP-Seq data analysis workflow
Sequence quality was assessed using FastQC (v0.11.2) [33] and quality trimming was carried out using the FASTX toolkit (v0.0.13) (http://hannonlab.cshl.edu/fastx_toolkit) to filter bases with less than the Phred33 score of 30 and to remove sequences with less than 50 bases. The quality trimmed reads were mapped to the AP13 v1.1 reference genome, Panicum virgatum, (https://phytozome. jgi.doe.gov/pz/portal.html#!info?alias=Org_Pvirgatum), using bowtie2 (v2.2.6) [34] with default parameters, except the number of mismatches was set to 1. Peaks marked by H3K9 me2 and H4K12 ac were identified using Spatial Clustering for Identification of ChIP-Enriched Regions (SICER; v1.1) [35] for both genotypes, using the following parameters: redundancy threshold of 10; window size of 200; gap size of 600; effective genome size as a fraction of reference genome of 0.7; and false discovery rate (FDR) threshold controlling significance at 0.01, for each replicate/sample with respect to their input control. The set of peaks present in at least two replicates was used for further analysis. A gene was considered methylated (H3K9me2-marked) or acetylated (H4K12ac-marked) only if it was overlapping (based on known annotated genes) with peak coordinates by at least one base. To study the pattern of histone binding in different regions of the genome, we also calculated distribution of peaks identified in different genomic regions (5′ UTR, exon, intron, 3′ UTR, upstream, downstream, and non-genic regions using PAVIS (v1.8) [36] based on the annotation of known genes. For each peak, the distance from the peak to the TSS was determined and plotted using annotations from Phytozome for the AP13 v1.1reference genome. The percentage of genes for each modification was calculated as the total number of genes with peaks divided by the number of genes for each modification and multiplied by 100. As routinely used in expression profiling of genes associated with histone peaks in other biological systems, once we obtained the nearest annotated genes, we performed Gene Ontology (GO) analysis to assign gene functions such as biological processes, molecular functions, and cellular components, before conducting Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. Among thousands of genes that were differentially marked or enriched using ChIP-Seq analysis, five genes (1-5) were selected for functional validation based on: i) their role in phenylpropanoid-monolignol pathway; ii) fold change (Log2FC) > 2; and iii) FDR < 0.01.

ChIP-quantitative real time (RT)-PCR (qPCR)
To quantitatively measure the amplification of ChIP-DNA, quantitative real-time (RT)-PCR (qPCR) validation was performed with an ABI 7500 real-time PCR (Applied Biosystems, Foster City, CA). The five genes selected here belong to the phenylpropanoid pathway that represents the active (H4K12 ac ) and inactive (H3K9 me2 ) modifications from ChIP-Seq analysis. The primers were designed using the TaqMan qPCR primer design tool (GenScript USA Inc., Piscataway, NJ), and the primer pairs used to amplify the ChIP-enriched regions are presented in Additional file 1: Table S1. qPCR was carried out in a 25 μL reaction that contained 10 ng of immunoprecipitated DNA, 10 μM each of forward and reverse primers, and 12.5 μl of SYBR Green PCR Master Mix (Germantown, MD). The qPCR reaction settings were: 95°C for 10 min, followed by 40 cycles of 95°C for 15 s and 65°C for 1 min.
We used three replicates to test the repeatability of qPCR validations. The universal cons7 primer was used as a constitutive control to determine normalized expression for all samples. Primer efficiency was determined by using 2-ΔΔCT method [37], where ΔΔCT (CT of genes -CT of cons7) = tissue to be observed -leaf tissue (CT of genex -CT of cons7). Minitab 17 [38] was used to analyze the normalized CT values (ΔΔCT), and the expression results were presented as mean ± SE. To test the multiple comparisons between the sample means, one-way ANOVA was conducted on normalized CT values collected from qPCR analysis.

Distribution of histone modifications
Overall distributions of H3K9 me2 and H4K12 ac peaks were obtained using reads that mapped to all chromosomes and contigs greater than 5 kb, and contigs between 2 kb and 5 kb that have predicted genes while calling peaks (Additional file 2: Figure S1). A representative graph that shows the distribution of histone marks on phenylalanine ammonia lyase (PAL) gene in chromosome 7b was visualized in AP13 and VS16 for H3K9 me2 and H4K12 ac modifications, using input (AP13 v1.1) as a positive control through the Integrative Genomics Viewer (IGV) [39] (Fig. 1). For H3K9 me2 and H4K12 ac in both AP13 and VS16, numerous peaks were identified in non-genic regions (over 40%), and less than 40% of peaks were in genic regions (Fig. 2). The range of nongenic region for AP13 for H3K9 me2 (399 bp -27,799 bp) and H4K12 ac (399 bp -32,399 bp) modifications and VS16 for H3K9 me2 199) and H4K12 ac (199 bp -29,399 bp) modifications has been observed in this study.

Genes related to peaks
By analyzing genes within 10 kb from the peaks, 24,163 and 21,918 genes were identified in AP13 and VS16, respectively. Of these, 5852 (AP13) and 3918 (VS16) geneassociated peaks were marked with H3K9 me2 , while 18,311 (AP13) and 18,000 (VS16) gene-associated peaks were marked with H4K12 ac ( Table 2). In AP13 and VS16, the percentage of methylated peaks located on genic regions was 24.22 and 17.88%, respectively. Also, the percentage of acetylated peaks located on genic regions in AP13 and VS16 was 75.78 and 82.12%, respectively. Furthermore, 71.86 and 63.31% of these genes carried both of the modifications (H3K9 me2 and H4K12 ac ) in AP13 and VS16, respectively ( Table 2).

Distances from peaks to TSSs
Based on nearest gene distribution using BEDTools (v2.21.0) [40], over 18,000 methylated and acetylated peaks spanned between 0 to + 2000 bases from TSSs in AP13. In addition, over 2000 acetylated peaks spanned between − 1 to − 2000 bases from TSSs in AP13. Over 12, 000 methylated and 18,000 acetylated peaks spanned between 0 to + 2000 bases from TSSs in VS16. (Fig. 3). Most enrichment was seen in the region 2 kb upstream or downstream of a TSS. The enrichment levels from both upstream and downstream of 2 Kb from TSS was seen to gradually decrease in both modifications for both genotypes (Additional file 3: Figure S2A-D). To understand the relationship between histone modifications (H3K9 me2 and H4K12 ac ) and TF binding, the peak distances to nearest genes were binned into 5 kb intervals (Additional file 4: Figure S3A-D). Of 30 bins with a nearest gene distance of 15 kb, 26 contain at least 100 peaks for H3K9 me2 modification in AP13 and VS16 (Additional file 4: Figure S3A and B). All 30 bins with a nearest gene distance of 15 kb for the H4K12 ac modification contain at least 100 peaks for both genotypes (Additional file 4: Figure S3C and D).

Epigenomic factors
The five major categories of epigenetic factors identified, along with the number of genes marked with H3K9 me2 in both genotypes (AP13/VS16), were methyl groups (186/159), methyl CpG-Binding domains (73/69), acetyltransferase groups ( Figure S4). In this study, epigenomic factors, methyl CpG binding domains, and acetyltransferase groups were most dominant.

Pathways and genes associated with peaks
Although several pathways and pathway-related genes were identified in the study, we only selected 16 pathways that had greater than or equal to10 genes associated with both H3K9 me2 and H4K12 ac in each genotype (Table 3). Of these, secondary metabolite biosynthesis pathways were predominant, with over 100 genes methylated and acetylated in both AP13 and VS16 genotypes. Therefore, greater emphasis was placed on genes, pathways, and TFs associated with biomass/biofuel production as well as biotic and abiotic stresses. Functions for the nearest genes to these peaks were assigned using KEGG pathway analysis, and most of these genes were functionally associated with cellular and metabolic processes, responses to stimuli, other biological regulation, etc. (Additional file 1: Tables S2, S3, S4 and S5).

C4 photosynthesis-related enzymes
An important intermediate in C4 carbon-fixation pathway, phosphoenol pyruvate (PEP) and two key pathway enzymes; i.e., phosphoenolpyruvate carboxylase (PEPC), that catalyzes PEP and carbonic anhydrase (CA), and two transporter genes viz., tonoplast dicarboxylate transporter  (DIT) and bile acid sodium symporter/pyruvate transporter (BASS) were identified as marked in both AP13 and VS16 for both modifications, whereas thiophosphate isomerase was only identified in VS16 for both modifications ( Table 4). The other genes, nicotinamide adenine dinucleotide phosphate (NADP)-dependent oxidoreductase, phosphoenolpyruvate carboxykinase (PEP-CK), and thiophosphate isomerase showed varied patterns in histone methylation and acetylation. For instance, in AP13, NADP-dependent oxidoreductase and PEP-CK were only methylated, while thiophosphate isomerase was neither methylated nor acetylated.

Photorespiration pathway-related genes
Marked by H3K9 me2 and H4K12 ac in both the genotypes were eight key photorespiratory enzymes: 2-oxoglutarate  Table 4). Of these eight genes, 2-OG, GLU, and GLN showed binding affinity to H3K9 me2 and H4K12 ac in both genotypes.

Phenylpropanoid pathway-related genes
Utilizing shikimate pathway intermediates, phenylpropanoid metabolism synthesizes a diverse array of secondary metabolites. Many genes and metabolites associated with phenylpropanoid pathways play a key role in feedstock quality and they were identified. Among these, the abundantly marked genes/metabolites in AP13 and VS16 for both H3K9 me2 and H4K12 ac include hydroxycinnamoyl-CoA:shikimate hydroxycinnamoyl transferase (HCT), cinnamyl alcohol dehydrogenase (CAD), phenylalanine    Table 4).

Peaks associated with transcription factors
In total, 653 and 650 significantly enriched peaks were annotated as TFs and were methylated and acetylated, respectively, in AP13, while 431 and 579 significantly methylated and acetylated peaks, respectively, were annotated as TFs in VS16. Among these, 17 TFs were differentially marked (H3K9 me2 and H4K12 ac ) between the two genotypes ( Table 6; Additional file 1: Tables  S34, S35,

Quantitative PCR validation
Genome-wide ChIP-Seq analysis revealed 17,364 and 13, 878 genes associated with peaks between AP13 and VS16, respectively, for both modifications (H3K9 me2 and H4K12 ac ). Among these, five genes (PAL, 4CL, HCT, CAD, and CCoAOMT) from the phenylpropanoidmonolignol pathway, differentially marked with histone modifications, were selected to confirm the expression level by qPCR. The significant (P < 0.05) variations, in the relative expressions, of all five genes indicated that they were more marked in AP13 than in VS16. Further, for both H3K9 me2 and H4K12 ac in these two genotypes the levels of expression were higher for CAD and HCT, while lower for CCoAOMT (Fig. 4).

Discussion
Histone modifications, such as H3K9 me2 , H4K12 ac , H3K27 me2 / me3 , H3K27 me1 , H3K4 me3 , and H3K9 me3 , bind either to heterochromatic or euchromatic regions and orchestrate epigenetic regulation of gene expression in a tissue-specific or genome-wide manner. These modifications have been suggested to play an important role in growth, development, biotic stresses, and abiotic stresses in plants [8,[11][12][13]. Similar to our previous study in common bean [8] and other studies in rice [7], Arabidopsis [41], cotton [42], and a diatom (Phaeodactylum tricornutum) [43], in this study, a genome-wide screening to understand the binding of two histone marks associated with gene inactivation (H3K9 me2 ) and activation (H4K12 ac ) was performed using ChIP-Seq analyses in two diverse switchgrass genotypes (AP13 and VS16)  belonging to two ecotypes. The percentage of mapped reads to the reference genome in rice ranged between 91.70-98.90% [7], while in switchgrass we obtained mapping rates of 82.32 and 73.42% for H3K9 me2 as wells as 81.20 and 68.43% for H4K12 ac in AP13 and VS16, respectively, which are attributed to the genotype-specific adaptations, variations in GC-content, repeat regions, genome duplications and ploidies [6].

Histone modification patterns in switchgrass
With only slight variations between the two genotypes, genome-wide histone modification patterns were similar (Additional file 2: Figure S3A-D) when compared to gene expression data from our recent effort [30]. More genes were associated with peaks marked with H4K12 ac than H3K9 me2 in both genotypes, which is attributed to the activation of genotype-specific genes (Additional file 2: Figure S1). Even though the overall distribution pattern of histone marks is similar, slight variations existed between the two genotypes, which is reflected in the distribution of peaks in the respective genomes. Also, more genes were marked with H4K12 ac than H3K9 me2 in both the genotypes. Moreover, bivalent histone modifications were commonly marked on specific genes for both genotypes, suggesting a significant degree of overlap between these two epigenetic marks. Our findings are in agreement with genome-wide histone modification studies in rice [7,44] and barley [45]. For example, in rice both activating (H3K4 me3 ) and repressive (H3K27 me3 ) histone modifications formed strong peaks within 1 kb downstream of the TSSs [7,46]. ChIP-Seq analysis in Populus revealed that~90% of ARK1 and ARK2 binding regions were within 1 kb of the TSSs [47]. In rice the enrichment of histone acetylation and methylation marks showed at 500 bp downstream from the TSSs [48]. The H3K4 me3 distribution showed a clear peak around 300 bp downstream of the TSS in Arabidopsis [49]. Additionally, in Arabidopsis H3K9 me2 and H3K4 me enrichment spanned around the entire gene body [50,51]. Previous reports in rice [46] and humans [52] demonstrated that active histone marks are more enriched in bidirectional promoters compared to unidirectional promoters. A detailed analysis of the genomic locations associated with the peaks for H3K9 me2 and H4K12 ac indicated higher affinity levels of H4K12 ac to switchgrass genes rather than H3K9 me2 , suggesting the role of active histone marks in transcriptional activation of switchgrass genes. However, a greater number (> 10%) of peaks associated with H4K12 ac and H3K9 me2 marks spanned 5 kb upstream and downstream from TSSs in AP13 when compared to VS16.
Recently, the preferential binding of TFs (~86%) between − 1000 to + 200 bp from a TSS has been identified in A. thaliana [53]. In Populus, the majority of binding regions for all TFs analyzed have shown a general preference for binding near genes, and they were specifically enriched within 5 kb of the TSSs [47]. Histone-mediated transcription initiation/inactivation occurs by binding to promoter and enhancer regions and also by modulating transcription elongation [48]. Analysis of histone modification profiles with bins of increasing gene length in Arabidopsis revealed that most peaks have been localized around 480 bp downstream of the TSSs. The shape and width of the H3K36 ac peaks were independent of gene length [12]. Similarly, in this study most H3K9 me2 and H4K12 ac peaks were localized between 0 and 499 bp downstream, supporting the idea that peak position, shape, and length are independent of gene length. Transcription factors and DNA binding proteins were highly enriched in the cluster of genes associated with H3K27 me3 in tissue-specific repression as reported in Arabidopsis (6). Genome-wide comparative analyses using two histone marks, H3K9 ac and H3K4 me3 , in A. thaliana and A. arenosa revealed coordination between histone modifications and gene expression within and between these species, including allopolyploids [6].
Pathways and genes related to biomass/biofuel production Further, we narrowed our search to identify histone markings on genes related to C4-photosynthesis andphotorespiration pathways, and phenylpropanoid-monolignol pathway, as they play a key role in biomass and biofuel production [54]. Using three histone acetylation marks (H3K9 ac , H4K5 ac , and H3K18 ac ) and two histone methylation marks (H3K4 me3 and H3K4 me2 ) and targeting five C4 genes (CA, PPDK, ME, PEPCK, and RBCS2), comparative ChIP-Seq analyses reported in maize, sorghum (Sorghum bicolor), and foxtail millet (Setaria italica) revealed a role for histone modifications in light-specific and cell type- Fig. 4 Differential expression of genes between AP13 and VS16 in H3K9 me2 and H4K12 ac modifications using Quantitative-Real-Time PCR (qPCR) analysis. The normalized CT values (ΔΔCT) from qPCR analysis were collected and analyzed by using Minitab 17, and the expression results were presented as mean ± SE. One-way ANOVA was performed on qPCR experiments for multiple comparisons between the mean of samples specific responses [55]. Supporting this idea, our analyses here showed preferential binding of histone marks (H3K9 me2 and H4K12 ac ) between AP13 and VS16 in five prominent C4-pathway genes (DIT, PEP, PEPC, BASS and CA) five photorespiration-related genes (2-OG, Glu, 3-PGA, RuBisCo, and GLN); and, six monolignol pathwayrelated genes (C4H, PAL, 4CL, HCT, CAD, and CCoAOMT). Genes that mediate phenylpropanoid pathways play a key role in feedstock quality. A recent histone binding study in Eucalyptus grandis reported that about 8% of phenylpropanoid pathway genes were marked by H3K4 me3 [56]. In this study, the histone modification patterns for 430 monolignol pathway-related genes between two ecotypes were mostly correlated (Fig. 5), supporting the idea of a common histone modification code between genotypes and possibly across species, as has been suggested.
Interestingly, this study identified several abiotic and biotic stress-related genes differentially marked by H3K9 me2 and H4K12 ac in VS16 and AP13, which may contribute to differing genetic potential of these two contrasting ecotypes (upland and lowland). Here, we identified H3K9 me2 and H4K12 ac modifications on different stress-related genes abundantly marked in both genotypes such as disease-resistance, heat tolerance, and flood tolerance genes. Our recent transcriptome analyses in switchgrass genotypes AP13 and VS16 identified over 200 significantly enriched stress-related genes [8]. Similarly in this study, prominently associated with H3K9 me2 and H4K12 ac were genes that encode the leucine rich repeats (LRR)-family of proteins, NBS-LRR, NB-ARC, LRR-NB-ARC, and CC-NBS-LRR (disease resistance); the heat shock proteins (HSPs), DNAJ and DNAK (heat tolerance); and the subtilase family protein, LRR-transmembrane protein kinase and adenine nucleotide alpha hydrolase-like superfamily protein (flooding tolerance) ( Table 5). Even though the overall binding patterns were similar, slight variations in markings were observed between the two ecotypes. Therefore, changes noticed in histone bindings represent inherent differences between the upland (VS16) and lowland (AP13) switchgrass ecotypes.

Quantitative PCR validation
It is important to validate differentially marked genes in switchgrass, a prominent biofuel crop, especially in regard to the phenylpropanoid-monolignol biosynthesis pathway. Previous reports suggested that over ten different genes play an important role in the phenylpropanoidmonolignol pathway in switchgrass [57][58][59], including the five genes selected in this study (PAL, 4CL, HCT, CAD, and CCoAOMT). Though transcriptional networks associated with the phenylpropanoid pathway have been reported [60], studies regarding the regulation of phenylpropanoid genes by histones in plants are limited. In maize, acetylation of H3K9/K14 in the promoter region resulted in activation of the anthocyanin biosynthetic gene, A1, by interacting with EMSY-related factor, thus linking TFs with chromatin modifications [61]. Our Fig. 5 Differentially marked histone modifications (H3K9 me2 and H4K12 ac ) on 430 selected phenylpropanoid-monolignol pathway-related genes between AP13 and VS16 using heatmap analysis results suggested that all five genes selected for ChIP-qPCR from the phenylpropanoid pathway were significantly marked by both histone modifications in both genotypes. Recently, comparative RNA-Seq analyses in switchgrass suggested that phenylpropanoid pathway genes were abundantly expressed in crowns and rhizomes of lowland cultivar, Kanlow, when compared to upland cultivar, Summer [62]. Similarly, comparative ChIP-Seq and ChIP-qPCR analyses as reported here in switchgrass suggest that differentially methylated and acetylated histone modification regions associated with the phenylpropanoid pathway genes were relatively abundant in lowland cultivar, AP13, when compared to upland cultivar, VS16. This suggests the role of active and inactive histone modification marks in regulating genes associated with biomass yield and biofuel production. Further, it is reasonable to assume that phenylpropanoid pathway genes may undergo a combinatorial transcriptional regulation involving TFs and histone modifications.

Conclusions
Our genome-wide screening for histone methylation (H3K9 me2 ) and acetylation (H4K12 ac ) marks in two contrasting ecotypes of switchgrass identified similarities and disparities between AP13 and VS16 at the epigenome level. The distribution of significantly enriched peaks spanned 0-2000 bp around TSSs, and most of these peaks were observed in non-genic regions. Our comprehensive analysis of annotated peaks facilitated the identification of key regulatory marks, genes, and TFs involved in biosynthetic pathways, including C4 photosynthesis, photorespiration, phenylpropanoid-monolignol, as well as biotic and abiotic stresses.
In summary, we integrated ChIP-seq and RNA-seq datasets to explore the interactions between histone modifications and TFs/transcripts related with gene expression and their patterns across two contrasting genotypes of switchgrass, AP13 and VS16. Furthermore, the integrated methodology and results highlighted here will guide future studies in exploring similar interactions in gene expression and regulation.