The symbiotic relationship between cnidarians and dinoflagellates is one of the most widespread endosymbiosis in our oceans and provides the ecological basis of coral reef ecosystems. Although many studies have been undertaken to unravel the molecular mechanisms underlying these symbioses, we still know little about the epigenetic mechanisms that control the transcriptional responses to symbiosis.
Here, we used the model organism Exaiptasia diaphana to study the genome-wide patterns and putative functions of the histone modifications H3K27ac, H3K4me3, H3K9ac, H3K36me3, and H3K27me3 in symbiosis. While we find that their functions are generally conserved, we observed that colocalization of more than one modification and or DNA methylation correlated with significantly higher gene expression, suggesting a cooperative action of histone modifications and DNA methylation in promoting gene expression. Analysis of symbiosis genes revealed that activating histone modifications predominantly associated with symbiosis-induced genes involved in glucose metabolism, nitrogen transport, amino acid biosynthesis, and organism growth while symbiosis-suppressed genes were involved in catabolic processes.
Our results provide new insights into the mechanisms of prominent histone modifications and their interaction with DNA methylation in regulating symbiosis in cnidarians.
Coral reefs are often considered the rainforests of the sea, as they form marine-biodiversity hotspots. Reef ecosystem health directly depends on symbiotic cnidarians, such as corals and anemones, that provide essential habitats for a myriad of marine organisms. To thrive in the oligotrophic environment of tropical oceans, corals, and other symbiotic cnidarians, depend on an intimate endosymbiosis with photosynthetic dinoflagellates of the family Symbiodiniaceae, also known as zooxanthellae [1,2,3]. Living within the host’s gastrodermal cells, the symbionts provide their hosts with over 90% of their total energy demands , making these symbiotic relationships vital for the functioning of the coral reef ecosystem. The disruption of this host-symbiont relationship, also known as bleaching, can result in extensive mortality and subsequent degradation and loss of entire coral reefs [2,3,4]. Significant efforts have been made to understand the molecular mechanism underlying this relationship . However, there are still substantial knowledge gaps in our understanding of the molecular underpinnings of these relationships, especially pertaining to the role of epigenetic mechanisms in regulating the interactions between the host and the symbionts, which remain elusive in cnidarian symbiosis research . The uptake and maintenance of symbionts require specific host responses, such as the suppression of the immune system [7,8,9,10] and the regulation of nutrient fluxes to control symbiont proliferation , to maintain a stable symbiotic relationship. Such responses are mediated through transcriptional changes that are known to be regulated via epigenetic mechanisms in other organisms [9, 12,13,14,15], and some endosymbionts have even been shown to evoke such responses by directly modifying the epigenome of their hosts [9, 10]. However, while many recent studies have highlighted the importance of epigenetic mechanisms in maintaining symbiotic relationships in plants and animals [16,17,18], only one study looking at the role of DNA methylation in symbiosis has been conducted in zooxanthellate cnidarians [6, 19].
Eukaryotic genomes are packaged in the form of a DNA-protein complex termed chromatin. The structural subunit of chromatin is known as the nucleosome, which consists of a core protein octamer and a stretch of ~147 bp of DNA that is wound around it. The protein octamer comprises two of each of the core histones H2A, H2B, H3, and H4. Linker histones such as H1 sit externally at the base of nucleosomes, providing further stability to it and contributing to chromatin compaction [20, 21]. This essential organization of histones aids in the folding of the DNA into a higher-order chromatin structures [20,21,22]. The N-terminal tail of histone proteins can include reversible covalent changes termed post-translational modifications (PTMs), which participate in the formation, folding, and structural/functional regulation of chromatin structure and have thus a profound role in epigenetic regulation of gene expression and genome stability. These PTMs (along with several histone variants) are part of the epigenetic mechanism known as “histone code” . Along with these modifications, other epigenetic mechanisms such as DNA methylation and small RNAs collectively influence the chromatin structure, and consequently, the accessibility of the genetic information [20, 21].
Histone modifications can affect gene regulation differently depending on their type and location in the genome (Table 1). Among the different PTMs, acetylation and methylation of specific histone tail residues have been most extensively studied [12, 13], and they have been found to promote repressing and activating roles in the regulation of gene expression . In general, activator complexes methylate or acetylate-specific amino acid residues in tails of histones bound to gene promoter regions, thereby destabilizing the nucleosome-DNA interaction and facilitating the assembly of the transcriptional machinery at the promoter. However, repressor complexes demethylate/deacetylate histone tails and strengthen the DNA-histone interaction, resulting in hindered accessibility of the respective genomic regions for the transcriptional machinery .
Methylation of histone tails occurs mainly at lysine (K) and arginine (R) residues, most commonly observed as mono-, di-, or trimethylation of the lysine residues on H3 and H4 histone tails [23,24,25]. Methylation of H3K36 and H3K4, for instance, act as activating histone modifications, while H3K27 methylation has a role in repressing the gene expression [26,27,28,29]. Similarly, H3K27ac and H3K9ac modifications are associated with active transcription and are predominantly associated with promoter and enhancer regions [26,27,28,29]. Histone methyltransferase (HMT) are histone-modifying enzymes that catalyze the transfer of methyl groups to the targeted residues (lysine and arginine) through a domain known as SET domain [30, 31]. Acetylation and deacetylation of histone tails, on the other hand, are catalyzed by histone acetyltransferases (HATs) and histone deacetylases (HDACs), respectively. Different chromatin-modifying enzymes, including histone deacetylases and histone lysine methyltransferases, function through multi-protein complexes that can also interact with methyl-CpG-binding proteins, thereby linking mechanisms of histone modifications to the biochemical mechanism that maintains and modifies DNA methylation [32, 33]. This implies that the enzymatic control of different epigenetic mechanisms is linked via crosstalk, and hence mutually interactive in regulating gene expression . In general, enzymes involved in depositing these chemical modifications (acetyl, methyl, etc.) onto the chromatin at the specific location are known as writers. In contrast, those which remove such modifications are called erasers [35, 36].
Despite the importance of the cnidarian Symbiodiniaceae relationship for ecosystem functioning, we still know very little about the role of epigenetic mechanisms, and specifically histone modifications, in the regulation of host-symbiont interactions. Here, we profiled the genome-wide association of the histone modifications H3K27me3, H3K36me3 and H3K4me, H3K27ac, and H3K9ac, in the cnidarian symbiosis model Exaiptasia diaphana (Aiptasia). We describe their genetic context, their correlation with CpG methylation (mCpG), and gene expression, as well as their association with and putative regulation of symbiosis genes.
Genome-wide distributions of histone modifications in E. diaphana and their correlation with CpG methylation
To understand the regulatory function of prominent histone modifications and their role in symbiosis, we performed ChIP-seq experiments in symbiotic E. diaphana and analyzed the genomic distribution of five major modifications; H3K27me3, H3K36me3, H3K27ac, H3K4me3, and H3K9ac (Additional file 1: Table ST1-ST5), as well as their correlations with respect to DNA methylation  and gene expression . To do this, we first called “peaks” from all ChIP-seq data, which are regions in the genome that are actively bound by a respective histone modification (see “Methods” for more details). A circular visualization plot of all five histone modifications with gene and repeat density over the E. diaphana genome is shown in Additional file 2: Fig. S1. We used the term “peak” throughout the manuscript to refer to these bound regions as well as their signal intensity relative to the input control.
For initial validation purposes, we compared the distribution of histone modifications between active and inactive regions of the E. diaphana genome . Repeat regions in the genome are mostly silenced  and are known to differ in bound histone modifications in comparison to non-repeat, i.e., genic, regions . Genome-wide analysis showed that the modifications H3K27me3 and H3K36me3 had significantly higher peaks (T-test; p < 0.0001) in repeat regions compared to genic regions. In contrast, H3K4me3, H3K27ac, and H3K9ac had significantly higher peaks (T-test; p < 0.0001) in genic regions (Fig. 1A). This suggested that the transcriptional suppression of repeat elements in the E. diaphana genome aligns with higher H3K27me3 and H3K36me3 signals, and simultaneously lower peak signals for H3K4me3, H3K27ac, and H3K9ac.
Next, we determined the peaks for each histone modification around all protein-coding genes in the E. diaphana genome. Peaks of H3K27me3 and H3K36me3 were prevalent in the gene body and promoter regions, but not the transcriptional start site (TSS). In contrast, peaks of H3K4me3, H3K27ac, and H3K9ac exhibited a bimodal peak pattern in the TSS region, with a smaller peak around the promoter region and a prominent peak around the first exon, while gene bodies featured comparatively lower peaks (Fig. 1B, Additional file 2: Fig. S2A and S2B). To investigate the relationship between the different histone modifications and mCpG, we classified all genes from E. diaphana as either methylated (n=8018) or non-methylated (n=21,322) based on their methylation density and methylation level (see “Methods” for more details). We found that H3K27me3 was predominantly present in the gene body of non-methylated genes, while H3K36me3, H3K4me3, H3K27ac, and H3K9ac were associated with methylated genes (adjusted p < 0.01) (Additional file 2: Fig. S2C – S2G, Additional file 1: Table ST6).
Interestingly, however, analysis of the core nucleosome regions of all five histone modifications showed either very low or no CpG methylation (Fig. 1C), suggesting that the DNA wound around the histone octamer core containing these modifications is mCpG depleted.
To determine how the different histone modifications are correlated with CpG methylation and gene expression, we analyzed their associations with mCpG, GC content, and gene expression  using a principal component analysis (Fig. 1D). All parameters, i.e., histone peak score, mCpG ratio, and GC content, were averaged for each gene. We observed that all histone modifications aligned on the same plane of the second principal component, along with gene expression. This suggests a tighter relationship between histone modifications and gene expressions compared to mCpG or GC content. Conversely, only histone modifications exhibiting gene body prevalence aligned with mCpG, GC content, and gene expression along the first principal component. This is likely because mCpG and GC content also display strong gene body prevalence. Further, gene body prevalent histone modifications (H3K27me3 and H3K36me3) and TSS-prevalent modifications (H3K4me3, H3K27ac, and H3K9ac) clustered together, respectively (Fig. 1D). This suggests that the distribution patterns of mCpG and GC content are more similar to those of gene body prevalent histone modifications. To further confirm this, we performed linear regression analyses to identify potential interactions between histone modifications at each gene (Additional file 2: Fig. S2H – S2K). We found a positive correlation between H3K4me3 and H3K9ac (R2= 0.27, p = 0.0023), H3K27ac and H3K9ac (R2= 0.36, p = 0.0014), and H3K27ac and H3K4me3 (R2= 0.43, p = 0.004), which suggests that a substantial number of genes could potentially be bound by more than one of the three TSS-dominated histone modifications (Additional file 2: Fig. S2H, S2I). In contrast, the correlation between gene body-dominated histone modifications (H3K27me3 and H3K36me3) was very weak (R2= 0.052, p = 0.045) (Additional file 2: Fig. S2J, S2K). Actual examples of the distribution of all five histone modifications over a methylated (AIPGENE18521) and non-methylated (AIPGENE18575) gene are shown in Fig. 1E.
Histone modification positioning and number as a function of transcription
Nucleosome positioning and spacing has previously been shown to correlate with gene expression levels . To study this in E. diaphana, we analyzed how the positions of the different histone modifications in the promoter and gene body regions vary and correlate with gene expression. To do this, we used the peak width as a parameter for position, with more precise nucleosome positioning reflecting tighter localization, and hence, narrower peak widths at a given position. For each modification, we compared the average peak width in the promoter and gene body separately (Fig. 2A) across lowly, intermediately, and highly expressed genes. We found that the average peak width of activating histone modifications in the promoter region and the gene body significantly increased with gene expression. Only the repressive modification H3K27me3 did not show a similar trend, and peak width remained constant across all expression levels.
To investigate the effect of multiple peaks on gene expression, we selected genes with no peaks and compared their average expression level against genes with a single peak from one modification. We found that genes that only have a peak for the repressive modification (H3K27me3) exhibited a significantly lower average expression than genes without any peaks. At the same time, all genes with active modifications showed higher gene expression values than genes without any peak or with H3K27me3 (Fig. 2B). Furthermore, we found a general increment in gene expression with an increasing number of peaks for active histone modifications. Meanwhile repressive histone modification peak counts, H3K27me3, showed a negative and weak relation with expression (Fig. 2C). We also examined the cumulative effect of multiple histone modifications on gene expression. Interestingly, we observed that the average gene expression level was higher if genes were associated with more than one histone modification, with every additional modification resulting in significantly higher expression levels of associated genes, as long as the repressive modification H3K27me3 was not included (Fig. 2D). Inclusion of H3K27me3 consistently correlated with significantly lower gene expression levels, further confirming its repressive effect.
Histone modifications in E. diaphana correlate with gene expression
To further analyze the correlation of the different histone modifications with mCpG and gene expression, we divided all the protein-coding genes from E. diaphana into six categories based on their transcription level and DNA methylation status. We observed a positive correlation between histone peak height and gene expression for H3K36me3 and for all TSS-prevalent histone modification (H3K4me3, H3K27ac, and H3K9ac), and only H3K27me3 displayed a negative correlation (Fig. 3A). This finding affirmed our previous analysis showing that H3K27me3 peak counts decreased with increasing gene expression (Fig. 2C).
To further investigate the potential interactions of histone modifications and DNA methylation, we plotted the average peak heights for every histone modification for low, intermediate, and highly expressed genes, each with and without DNA methylation respectively (Fig. 3B). We found that H3K27me3 showed a significant negative correlation (p < 2.2 × 10−16) with gene expression both when present in the promoter or the gene body, and this effect was even more pronounced in methylated genes. In contrast, H3K36me3 showed a positive correlation with DNA methylation and gene expression (Fig. 3B), with H3K36me3 peak height positively correlating with increasing expression in methylated genes. Similarly, TSS-prevalent histone modifications, i.e., H3K4me3, H3K27ac, and H3K9ac, also showed a positive correlation with gene expression and methylation (p < 2.2 × 10−16), and this effect was more pronounced for peaks in the promoter region (Fig. 3A, B).
Histone modifications regulate the transcriptional response to symbiosis
To investigate the role of histone modifications in the regulation of the symbiotic relationship between E. diaphana and its dinoflagellate symbionts, we analyzed the correlation between peak occupancy for each histone modification and gene expression across 731 previously identified symbiosis-associated genes . We first categorized these 731 symbiosis-associated genes into symbiosis-repressed (n=365) and symbiosis-induced (n=366) genes. We found that most symbiosis genes (544; 74.4% of the 731 genes) were associated with at least one of the five histone modifications we analyzed (Fig. 4A, Table 2, Additional file 1: Table ST7 – ST11). Active histone modifications (i.e., H3K36me3, H3K4me3, H3K27ac, and H3K9ac) showed a significantly higher association with symbiosis-induced genes (p < 0.01) while the repressive histone modification (H3K27me3) had an almost equal number of peaks in both categories.
To confirm the relationship between peak height and gene expression, we compared the average histone peak height of every modification across symbiosis-repressed and symbiosis-induced genes. We found that the repressive modification H3K27me3 had significantly higher peaks in symbiosis-repressed genes (p < 2.2 × 10−16, Fig. 4B), while all active modifications (H3K36me3, H3K4me3, H3K27ac, and H3K9ac) had significantly higher peaks (p < 2.2 × 10−16) in symbiosis-induced genes (Fig. 4C–F). This finding confirmed that histone modifications play an active role in the regulation of the symbiotic relationship between E. diaphana and its dinoflagellate symbiont.
As a final step of validation, we analyzed the correlation between histone modifications and changes in gene expression in response to symbiosis, we divided both symbiosis-induced (n=366) and symbiosis-repressed genes (n=365) into two groups based on their median gene expression fold change. We compared the profiles of the upper and lower 50th percentile for each of the histone modifications separately. Similar to our previous observation, we found that the repressive modification H3K27me3 showed a higher prevalence in symbiosis-repressed genes irrespective of the fold change (p < 2.2 × 10−16, Additional file 2: Fig. S3A). Conversely, our analysis on the active histone modifications; H3K36me3, H3K4me3, H3K27ac, and H3K9ac, also confirmed our previous findings of significantly higher peaks in symbiosis-induced genes (p < 2.2 × 10−16) (Additional file 2: Fig. S3B – S3E).
Histone modifications are involved in symbiosis-induced nutrient metabolism
To understand the role of histone modifications in the regulation of the symbiotic relationship, we performed gene ontology (GO) enrichment analyses for both symbiosis-repressed and symbiosis-induced genes for each histone modification individually (Additional file 1: Table ST12 – ST21, Additional file 2: Fig. S4). While each modification had several unique enriched biological functions, we found a considerable number of categories that were enriched across two or more histone modifications (Additional file 1: Table ST22 – ST23). Many of these shared categories within the symbiosis-induced genes were involved in amino acid metabolic processes, such as the regulation of cellular amino acid and protein metabolic process (Table 3). Furthermore, we found processes involved in the response to glucose and ammonium transport to be associated with multiple modifications, suggesting that the central function of this symbiotic relationship in driving host amino acid biosynthesis is regulated through histone modifications . In line with this, we also found shared enrichment of amino acid biosynthesis-related processes, including L-serine, glutamine, and methionine, as well as categories involved in amino acid transport. Similarly, we found several shared biological processes associated with organism growth to be enriched across multiple histone modifications, most notably GO terms related to central growth pathways such as the insulin, hippo, and the TORC1 pathways.
While many of the enriched categories in the symbiosis-induced genes were associated with anabolic processes, we found symbiosis-repressed genes associated with histone modifications to be predominantly involved in catabolic processes. These included the catabolism of molecules like xanthine and glycerol-3-phosphate but also amino acids such as sulfur amino acids, L-phenylalanine, betaine, arginine, and L-threonine, among others.
Interestingly, we also found enrichment of GO terms involved in various metabolic pathways that generate metabolites important for epigenetic modifications, such as acetyl-CoA/fatty-acyl-CoA, S-adenosylmethionine, methylation process, and lactate. Previous studies have shown that these metabolites serve as cofactors for the enzymes responsible for depositing the chemical modifications (acetyl and methyl) onto chromatin; chromatin writers [35, 36]. In addition, we found GO terms related to metabolites such as α-ketoglutarate and NAD+, which are essential cofactors for certain enzymes that remove chemical modifications; chromatin erasers [35, 36]. This suggests that the chromatin changes induced to regulate gene expression in response to symbiosis might be supplied by these processes and ultimately established through the respective writers and erasers.
The process of symbiosis establishment and maintenance requires changes in the cnidarian host’s cell function and specialization. Epigenetic mechanisms have been shown to play critical roles in symbiotic relationships of eukaryotic and bacterial cells . The general importance of histone modifications in host-microbe interactions has been acknowledged in plants, humans, and other invertebrates [16,17,18]. Through chemical signals and metabolites, endosymbionts can influence epigenomes of host cells and directly enable communication between the two partners [18, 41]. Interestingly, histone acetylase and deacetylase activity have been shown to be influenced by microbes and dietary factors [10, 42, 43]. Although epigenetic studies in cnidarians remain scarce, there is evidence that histone modifications may play a critical role in host-algae symbiosis mechanisms [6, 44,45,46]. Here, we report the first genomic landscape of five histone modifications, H3K27me3, H3K36me3, H3K4me3, H3K27ac, and H3K9ac, in a symbiotic cnidarian.
We find that their genomic distribution and putative primary functions align with observations made in other organisms , suggesting functional conservation of these histone modifications in E. diaphana. Further, our results revealed strong correlations between the histone modifications analyzed and transcriptional changes observed in response symbiosis. These findings collectively suggest a direct role for histone modifications in regulating the host’s response to symbiosis.
Conserved roles of histone modifications in regulating gene expression
The general explanation for the ability of histone modifications to enhance or repress transcription is that they affect the DNA-histone association and, thus, promote or suppress access for transcription factors and the transcriptional machinery to the DNA. As such, these modifications represent an essential mechanism for the epigenetic control of transcriptional responses in eukaryotes [13, 48,49,50]. In line with this, our analyses revealed highly significant correlations between histone modifications and gene expression. Analysis of the genomic distribution of H3K27me3 and H3K36me3 showed enrichment in repeat regions [48, 49] while the activating modifications H3K4me3, H3K27ac, and H3K9ac showed enrichment in the genic regions (Fig. 1A), as expected based on observations in other organisms [12, 51].
It is interesting to note, however, that we found a substantial number of genes associated with more than one histone modification, suggesting that several histone modifications might act on the same gene simultaneously (Additional file 2: Fig. S2H – S2K). Such a cooperative interaction in regulating gene expression was further supported by the finding that the number of active histone modifications present on genes positively correlated with gene expression levels, suggesting an additive effect. However, it needs to be pointed out that ChIP-seq data cannot inform if the modifications were present on the same DNA molecule or if they were just associated with the same gene but in different cells of the organism. This limitation is evident when looking at the lower average expression observed for genes that were bound by an activating histone modification as well as the repressive modification H3K27me3. Since an actively expressed gene is unlikely to be simultaneously associated with activating and repressive modification, it is more likely that the H3K27me3 association stems from cells where this gene was silenced. Since these cells would not contribute any transcripts for this gene to the whole organism RNA pool, this would reduce the observed overall expression level of the gene in the organism.
Crosstalk between histone code and DNA methylation
DNA has a determined nucleotide sequence that cannot be changed. However, it has been postulated that the transcription of the genetic information is partly regulated by epigenetic mechanisms such as the underlying histone modifications and DNA methylation. Our analyses of potential interactions of histone modifications and DNA methylation in regulating gene expression revealed strong correlations for all activating modifications that suggest crosstalk between these epigenetic mechanisms in E. diaphana. We observed that the average expression of genes associated with activating histone modifications was generally higher if they were also methylated (Fig. 3B), suggesting a cooperative interaction between activating histone modifications and DNA methylation. In contrast, we found that genes associated with the repressive modification H3K27me3 showed the opposite trend for methylated genes if the histone modification was found in the gene body. However, it should be noted that H3K27me3 was predominantly present in the gene body of non-methylated genes, while H3K36me3, H3K4me3, H3K27ac, and H3K9ac were associated with methylated genes (Additional file 2: Fig. S2C – S2G). While these results suggested a cooperative interaction, analysis of the core nucleosome regions of all five histone modifications showed either very low or no CpG methylation (Fig. 1C), indicating that they are present on the same genes but that their precise locations within the gene are mutually exclusive. In summary, our results are indicative of crosstalk between active histone modifications and DNA methylation in modulating gene expression, while repressive modifications associate predominantly with non-methylated genes to suppress their expression.
A model for the regulation of gene expression via histone modifications in E. diaphana
Based on our results, we propose a model for how the histone modifications analyzed here could regulate gene expression in E. diaphana. In addition, the model demonstrates how histone modification and DNA methylation crosstalk may be functioning in symbiotic cnidarians.
When a gene is silenced, it is bound by H3K27me3 in the promoter and gene body (Fig. 5A), which promotes repression through the polycomb complex. H3K27me3 has been shown to recruit PRC1 (polycomb repressive complex), which contributes to the compaction of the chromatin, leading to the formation of heterochromatin and its inaccessibility for transcription factors and the transcriptional machinery .
However, when a gene needs to be activated, it requires the preinitiation complex (PIC) to assemble at the promoter region of a gene and to recruit RNA Pol II to the promoter to build the transcription initiation complex [53,54,55]. Based on our results, and in line with previous findings [12, 28, 56], we propose that the presence of the activating histone modifications H3K27ac and H3K9ac around the promoter and the TSS promote access for transcription factors of the PIC to the gene promoter (Fig. 5B). Once the PIC is assembled, RNA Pol II can be recruited to form the transcription initiation complex and H3K27ac, H3K9ac, and H3K4me3 (TSS-dominated modifications) can act as a pause-release signal for Pol II to initiate transcription. The transcriptional elongation process is then supported by H3K4me3 and H3K36me3 . Meanwhile, low CpG methylation at the promoter further favors the attachment of the assembly of the PIC and the transcriptional complex , while high CpG methylation in the gene body prevents the assembly of the transcriptional machinery at cryptic promoter sequences within the gene body, which would lead to spurious transcripts and the production of truncated proteins . The crosstalk between histone modifications and DNA methylation is brought about through the interaction of histone-modifying enzymes. For instance, the histone methyltransferase Set2D is recruited by the active transcriptional complex and trimethylates H3K36 along the gene body. H3K36me3, in turn, is then actively bound by DNA methyltransferase 3b which methylated CpG within the gene body of actively transcribed genes . Together, histone modifications and DNA methylation create a chromatin landscape conducive of high gene expression and the production of full-length transcripts, while at the same time reducing transcriptional noise and spurious transcripts (Fig. 5B) .
The role of histone modifications in symbiosis
Our analyses revealed that genes associated with activating histone modifications were significantly enriched in the fraction of symbiosis-induced genes. This suggests that their increased expression in response to symbiosis is promoted via their association with activating histone modifications. Interestingly, we did not see the opposite trend for the repressive modification H3K27me3 which is associated with the same number of symbioses induced and repressed genes. However, analysis of H3K27me3 peak heights did show significantly higher peaks in symbiosis-repressed genes compared to symbiosis-induced ones. The fact that H3K27me3 peaks were significantly higher in symbiosis-repressed genes suggests that the association of these genes with H3K27me3 was evident in more host cells, which increased the number of ChIP-seq reads obtained, and thus the relative peak heights.
Analyses of the biological functions enriched in 2 or more histone modifications highlighted that the histone modifications studied associated with anabolic functions in symbiosis-induced genes and catabolic functions in symbiosis-repressed genes. In particular, processes associated with amino acid biosynthesis and growth suggested that these histone modifications are involved in regulating the metabolic response and growth in symbiotic anemones. However, these processes are also involved in the regulation of the symbiotic relationship itself. Endosymbiotic relationships are usually driven by synergies arising from the complementation of the host’s metabolic capabilities that enable the resulting metaorganism to thrive in nutrient-poor environments or to use previously inaccessible diets [62, 63], as is also the case for symbiotic anemones and corals. However, this intimate form of symbiosis requires maintaining a delicate balance of nutrient fluxes to provide nutrients to the symbionts, to keep benefiting from them, but at the same time ensure they do not over proliferate at the expense of the host. The maintenance of this balance in symbiotic cnidarians is achieved through the regulation of genes involved in ammonium assimilation and amino acid biosynthesis . We find that genes involved in the assimilation of waste ammonium and amino acid biosynthesis are predominantly associated with activating histone modifications. Further, the observed crosstalk of activating histone modifications and DNA methylation in driving higher expression of symbiosis-induced genes suggest a multi-layer epigenetic regulatory mechanism that may be critical for cnidarian symbiosis.
In summary, our analysis of the genome-wide distribution of five histone modifications in the symbiotic sea anemone Aiptasia identified generally conserved patterns for all the modifications. Correlations of histone modifications, DNA methylation patterns, and gene expression further suggest a cooperative function of histone modifications and DNA methylation in regulating gene expression. Analysis of previously identified symbiosis genes revealed a functionally consistent association with active and repressive histone modifications and DNA methylation, indicating that the maintenance of symbiosis-associated gene expression is likely mediated by the synchronous action of both epigenetic mechanisms. However, the biological interpretations of the results presented here are only first insights into research that clearly requires further expansion. We acknowledge that to further decipher the role of histone modifications in symbiosis, more ChIP-seq studies including aposymbiotic individuals (E. diaphana without its dinoflagellate symbionts) will be required. Nonetheless, our results on the biological functions align with recent observations of regulated gene accessibility and chromatin states in symbiotic anemones .
E. diaphana culture and maintenance
E. diaphana of the clonal strain CC7 , originating from North Carolina, was used in this study. Anemones were maintained in polycarbonate tubs with autoclaved seawater at 25°C. Animals were exposed for 12-h light/dark cycle at 20–40 μmol photons m2s1 light intensity. The anemones were fed twice weekly with freshly hatched Artemia nauplii (brine shrimp). For the experiment, three independent biological replicates were taken, each consisting of two individual anemones to provide enough material for all immunoprecipitations. Each of the three biological replicates was then processed independently for chromatin extraction, resulting in 3 independent chromatin samples.
The process of establishing a reproducible ChIP-seq protocol in E. diaphana, which so far has primarily been optimized for human, mice, and plant cell studies, included many quality control and optimization steps that require attention. In hopes of streamlining future attempts at ChIP-seq in other cnidarians, especially corals, we opted to optimize pre-immunoprecipitation steps to the point that kits could be confidently used thereafter. We used Zymo-Spin ChIP Kit (Zymo Research) to extract histone-bound DNA fragments; however, we applied minor adjustments to the pre-IP steps. In a recent publication , we published a summarized version of the protocol. In another publication , a summarized version of the protocol was presented. Here, we provide a detailed description of the protocol (see also Additional file 2); steps of validation and optimization are described in greater detail to, hopefully, allow future research to progress and further advance the field of epigenetic research in cnidarians.
Corresponding input controls for each of the three replicates were generated as suggested by the manufacturer. After validation of various histone antibodies (Additional file 2: Fig. S4), immunoprecipitation was conducted using a target-specific antibody to histone 3 acetylation at lysine 27 – H3K27ac (ab4729, Abcam), histone 3 trimethylation at lysine 4 – H3K4me3 (ab8580, Abcam), histone 3 acetylation at lysine 9 – H3K9ac (ab10812, Abcam), histone 3 trimethylation at lysine 36 – H3K36me3 (ab9050, Abcam), and histone 3 trimethylation at lysine 27 – H3K27me3 (ab6002, Abcam). Each of the three independent chromatin extractions (biological replicates) was divided into 6 aliquots to perform immunoprecipitations for each of the 5 histone modifications as well as retaining an input control for each biological replicate. This resulted in 15 immunoprecipitations and 3 input controls, one for each biological replicate. Upon validation of immunoprecipitations, using High Sensitivity DNA Reagents (Agilent Technologies, California, USA) on a Bioanalyzer, ChIP libraries were constructed using TruSeq Nano HT DNA kit (Illumina, California, USA).
Paired and single-end sequencing was performed for all 24 libraries at the Bioscience Core Lab (BLC) at the King Abdullah University of Science and Technology, Thuwal, KSA, with NextSeq 500. The ChIP-seq mapped files are deposited in NCBI SRA under accession number PRJNA826667.
ChIP-seq library sequencing resulted in 10 to 40 million single- and paired-end reads per replicate. Reads from 24 libraries (3 biological replicates for each of the 5 histone modifications; upon which 3 (H3K27ac, H3K27me3, and H3K4me3) modifications have common input controls library while 2 (H3K9ac and H3K36me3) have separate input controls) were checked for their quality with FASTQC toolkit . Adapters were removed and reads were cleaned using Trimmomatic . Subsequently, clean reads were uniquely mapped on the E. diaphana genome (http://exaiptasia diaphana.reefgenomics.org/)  using bowtie 1.1.2 with default parameters .
Identification and annotation of histone modification peaks
Bam files resulting from the mapping of the 24 libraries were used to call histone peaks separately for each modification using the Model-based Analysis of ChIP-Seq (MACS3: https://github.com/macs3-project/MACS/tree/master/MACS3) . Genomic regions for the five different modifications (3 replicates each) were identified through “macs3 callpeak -t treatment.bam -i input.bam -f BAM -g 2.7e+8-B --nomodel --d-min 10 --call-summits” parameters. Subsequently, the peak output of the replicates from each sample were combined for a replicate aware analysis using the Multiple Sample Peak Calling tool (MSPC) (https://github.com/Genometric/MSPC) , with “./mspc -I rep*.bed -r bio -w 1e-4 -s 1e-8” parameters. Two filters were applied for final peak calling; (1) in MACS3, we filtered each called peak by adjusted p < 0.01, and (2), in MSPC, consensus regions of peaks from MACS3 were called across the three replicates of each histone modification with adjusted p < 0.01. MSPC generated a single file from the three biological replicates with genome-wide peaks for each histone modification, which was used for all downstream analyses.
Each MSPC generated histone enrichment file with peak locations was annotated based on the E. diaphana genome (GFF3 file)  using ChIPseeker: An R/Bioconductor package for ChIP peak annotation, comparison, and visualization . This package annotates peaks based on genome features, e.g., genic or intergenic region, as well as distances to the 5′ and 3′ ends of each genomic feature (promoter/exon/intron). The annotated tables of all five histone enrichment peaks are provided in Additional file 1: Tables ST1–ST5.
Analysis of previously published RNA-seq and DNA methylation data
Gene expression data were obtained from a previously published study in our lab . This data was generated through a meta-analysis with random effects across four independent differential gene expression data sets. We identified a robust set of 731 genes involved in symbiosis (Additional file 1: Table ST24). Based on the median gene expression fold change, we classified these 731 genes into symbiosis-induced (n=366) and symbiosis-repressed (n=365) genes.
For DNA methylation, we used BS-seq raw reads from symbiotic E. diaphana previously published by Li et al. . We used the six replicates from the symbiotic condition to determine the methylation level following the steps described in Li et al. . After filtering the reads for their quality, we first mapped the reads to E. diaphana genome with bowtie 1.1.2  and performed methylation calling using Bismark-0.22.3 . We applied three filters to reduce false positives as mentioned in Li et al.  and annotated methylated locations using ChIPseeker.
GO enrichment of histone peaks
To analyze the functional enrichment of the histone peak-bound genes, we obtained the GO annotation from the genome annotation and analyzed it using topGO  using default settings. In order to test for the potential role of all five histone modifications in symbiosis, we used the list of 731 identified symbiosis genes as identified in Cui et al.  and matched them with the binding sites of each histone modification.
Average enrichment scores, plots, and heatmaps at genomic features of interest were generated using deepTools . Histone peak width and read counts for each associated gene and modification were determined and plotted in boxplots, and statistical tests to compare the means of two or more groups were performed in R (version 3.5.1 and R CRAN package: dyplr). In boxplots, the bottom and top of the box indicate the 25th and 75th percentile, respectively. The horizontal thick bars in the boxplot denote the medians. Whiskers indicate the 1.5X interquartile range (IQR).
Integrated Genome Browser - 9.1.8 (IGB)  has been used to explore and visually analyze different histone modification peaks and DNA methylation on the E. diaphana genome. Circular visualization plots were generated in R by circlize package .
Zhou J, Wang X, He K, Charron JB, Elling AA, Deng XW. Genome-wide profiling of histone H3 lysine 9 acetylation and dimethylation in Arabidopsis reveals correlation between multiple histone marks and gene expression. Plant Mol Biol. 2010;72:585–95.
Nagymihály M, Veluchamy A, Györgypál Z, Ariel F, Jégu T, Benhamed M, et al. Ploidy-dependent changes in the epigenome of symbiotic cells correlate with specific patterns of gene expression. Proc Natl Acad Sci USA. 2017;114:4543–8.
Bednar J, Horowitz RA, Grigoryev SA, Carruthers LM, Hansen JC, Koster AJ, et al. Nucleosomes, linker DNA, and linker histone form a unique structural motif that directs the higher-order folding and compaction of chromatin. Proc Natl Acad Sci USA. 1998;95:14173–8.
Nishida H, Suzuki T, Kondo S, Miura H, Fujimura Y, Hayashizaki Y. Histone H3 acetylated at lysine 9 in promoter is associated with low nucleosome density in the vicinity of transcription start site in human cell. Chromosom Res. 2006;14:203–11.
Young MD, Willson TA, Wakefield MJ, Trounson E, Hilton DJ, Blewitt ME, et al. ChIP-seq analysis reveals distinct H3K27me3 profiles that correlate with transcriptional activity. Nucleic Acids Res. 2011;39:7415–27.
Charlet J, Duymich CE, Lay FD, Mundbjerg K, Sørensen KD, Liang G, et al. Bivalent regions of cytosine methylation and H3K27 acetylation suggest an active role for DNA methylation at enhancers. Mol Cell. 2016;62:422–31.
Sunagawa S, Wilson EC, Thaler M, Smith ML, Caruso C, Pringle JR, et al. Generation and analysis of transcriptomic resources for a model system on the rise: the sea anemone Aiptasia pallida and its dinoflagellate endosymbiont. BMC Genomics. 2009;10:258.
M.A. and M.J.C. conceived the project. M.J.C. performed the ChIP-seq experiments; K.N., M.J.C., K.G.M., and G.C. analyzed the data. K.N., M.J.C., and M.A. wrote the first manuscript draft. K.N. and M.A. revised the manuscript. All authors read and approved the final manuscript.
Description of all the peaks estimated from three replicates for H3K27me3, using Multiple Sample Peak Calling (MSPC) program. Supplementary Table ST2. Description of all the peaks estimated from three replicates for H3K36me3, using Multiple Sample Peak Calling (MSPC) program. Supplementary Table ST3. Description of all the peaks estimated from three replicates for H3K4me3, using Multiple Sample Peak Calling (MSPC) program. Supplementary Table ST4. Description of all the peaks estimated from three replicates for H3K27ac, using Multiple Sample Peak Calling (MSPC) program. Supplementary Table ST5. Description of all the peaks estimated from three replicates for H3K9ac, using Multiple Sample Peak Calling (MSPC) program. Supplementary Table ST6. Description of all the peaks estimated from three replicates for each histone modification, using Multiple Sample Peak Calling (MSPC) program. Supplementary Table ST7. List of Symbiosis-repressed and -induced genes with H3K27me3 peaks. Supplementary Table ST8. List of Symbiosis-repressed and -induced genes with H3K36me3 peaks. Supplementary Table ST9. List of Symbiosis-repressed and -induced genes with H3K4me3 peaks. Supplementary Table ST10. List of Symbiosis-repressed and -induced genes with H3K27ac peaks. Supplementary Table ST11. List of Symbiosis-repressed and -induced genes with H3K9ac peaks. Supplementary Table ST12. GO IDs and terms of symbiosis induced genes associated with H3K27me3. Supplementary Table ST13. GO IDs and terms of symbiosis repressed genes associated with H3K27me3. Supplementary Table ST14. GO IDs and terms of symbiosis induced genes associated with H3K36me3. Supplementary Table ST15. GO IDs and terms of symbiosis repressed genes associated with H3K36me3. Supplementary Table ST16. GO IDs and terms of symbiosis induced genes associated with H3K4me3. Supplementary Table ST17. GO IDs and terms of symbiosis repressed genes associated with H3K4me3. Supplementary Table ST18. GO IDs and terms of symbiosis induced genes associated with H3K27ac. Supplementary Table ST19. GO IDs and terms of symbiosis repressed genes associated with H3K27ac. Supplementary Table ST20. GO IDs and terms of symbiosis induced genes associated with H3K9ac. Supplementary Table ST21. GO IDs and terms of symbiosis repressed genes associated with H3K9ac. Supplementary Table ST22. GO IDs and terms of symbiosis induced genes associated with one or more histone marks. Supplementary Table ST23. GO IDs and terms of symbiosis repressed genes associated with one or more histone marks. Supplementary Table ST24. List of all symbiosis induced and repressed genes list from Cui et al., 2019.
Circular visualization of histone modifications, gene and repeat density over the 7 largest scaffolds of E. diaphana genome. Figure S2. Genome-wide distribution of histone modifications in E. diaphana and their correlations. Figure S3. Upper and lower 50% percentile: histone modifications change within symbiosis induced and repressed genes. Figure S4. Gene Ontology bubble plots. Figure S5. Sequences alignment of histone 3 (H3) across species. Figure S6. Western blot of histone specific antibodies on total protein content of Aiptasia. Table S1. Table showing fixation buffer chemicals and its concentrations. Table S2: Table showing nucleic preparation buffer chemicals and its concentrations. Figure S7. Schematic representation of ChIP-seq protocol using Aiptasia.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
Nawaz, K., Cziesielski, M.J., Mariappan, K.G. et al. Histone modifications and DNA methylation act cooperatively in regulating symbiosis genes in the sea anemone Aiptasia.
BMC Biol20, 265 (2022). https://0-doi-org.brum.beds.ac.uk/10.1186/s12915-022-01469-y