Chromatin dynamics enable transcriptional rhythms in the cnidarian Nematostella vectensis

ATAC-seq was applied to measure high-resolution chromatin accessibility in N. vectensis under two light regimes, resembling maximum and minimum temporal behavioural activity. Nuclear integrity [19] was achieved by (i) dissolving animal tissue into single cell suspensions and (ii) separating released nuclei from other organelles and cytoplasmic debris to reduce non-nuclear DNA contamination, as presented in Fig 1A. Each ATAC-seq library was prepared using a transposition reaction from ~400,000 cells that were sampled from one individual animal, from circadian time 13 (CT13) to circadian time 45 (CT45), in 8 hour intervals. Libraries were sequenced from two independent biological replicates from the LD treatment and compared to two independent biological replicates collected at the same time, under DD. After filtering the PCR duplicates, generated from the library amplification process and irregular reads between the biological replicates, our ATAC-seq libraries showed a median depth of 10 million unique, high-quality mapped reads per sample (Tables 1, 2 and Fig 1). As presented in S1 Table, the biological replicates were highly similar (average R² = 0.96 SD ±0.014 and p-value<0.01), demonstrating highly reproducible data from N.vectensis whole animal sampling. Furthermore, the significant peaks (8366–36582 peaks with an irreproducible discovery rate (IDR) cut-off of 0.05) from both LD and DD treatments were clustered around transcriptional start sites (TSSs, Fig 1B and 1C and S2 Fig).

Many ATAC-seq peaks, from both the LD and DD samples, were mapped within 1500 bp upstream to the TSS (S2 Fig), marking the accessible chromatin of extended promoters. ATAC-seq libraries (LD and DD) were enriched, on average, with 20.87% (SD ±1.76) promoter regions, 20.35% (SD ±1.9) intron regions, and 43.1% intergenic regions (Proximal– 16.6% SD ±1.05, Distal– 26.5% SD ±2.65) that may act as distant regulatory elements (Fig 2A and 2B). Comparing the DD and LD libraries revealed time based clustering, which shows that chromatin accessibility is maintained in N. vectensis under constant conditions (Fig 2C). Moreover, gene ontology (GO) enrichment analysis of genes with accessible promoters showed that nucleic acid binding and transcription regulation activities remain similar between the two light regimes. Interestingly, DD-specific accessible gene promoters were highly enriched with rhodopsin-like, G-protein-coupled receptors (GPCRs), which are related to external signal transduction [20] (Fig 2D).

Previously, we showed an association between circadian locomotor activity rhythm and transcriptional profile in N. vectensis. Using Fourier analysis, diel rhythmicity (i.e., 24-h periodicity) was identified in many genes. From these, 180 transcripts exhibiting significant oscillations (G-factor >0.5) were selected for further analysis. Through K-means clustering, these transcripts were divided into five clusters, each representing a different peak time of chronological expression (S1C Fig) [15]. To further study the relationship between promoter accessibility and gene expression, 139 genes were selected, with promoter regions (within 1500 bp upstream of TSS) showing higher accessibility at CT13 (the time point with the strongest change in gene expression between the five clusters) than the average promoter accessibility, genome-wide. To test the association between all the data, we conducted a Pearson correlation test, finding a significant correlation between expression and accessibility (R² = 0.343 and P-value<0.01 –see S2 Table). Moreover, comparing accessibility and expression within each time point showed that promoters’ accessible sites oscillated during a 24h cycle and correlated to the syn-expression pattern of circadian genes, as shown in Fig 3A and 3B and S2 Table (Note: at CT 37, the correlation is not significant, although the overall pattern is visible). However, we have to remain skeptical, as not all genes in the presented list (S4 Table) correspond to the rule of accessibility and expression correlation. For example, a subset of genes that share a common TF site (CEB/P) show no significant correlation between expression and accessibility. These genes show relatively high accessibility throughout the experiment, with an average–log₂(FPKM) rate of 2.5 and a SE of ±0.3. These genes are related to core clock mechanism and development, as we show later.

Selective activation of functional regulatory DNA elements defines where TFs may bind and act. Therefore, to predict the identity of active TFs in treatment-specific peaks, the enrichment of sequence motifs was computed using the HOMER motif analysis tool [21]. The identified motifs were divided into three groups: (i) common motifs—motifs that were enriched in both treatments, relative to their abundance within the genome (see Table 3), (ii) LD-enriched (see Table 4), and (iii) DD-enriched (see Table 5). Many of the identified motifs correspond to binding sites of TFs with previously-identified roles in regulating rhythmic processes [22,23]. For example, the MYB motif, enriched in LD-specific peaks, was shown to have an essential role in circadian rhythm maintenance in Arabidopsis [24]. Interestingly MYB motif was enriched around the promoters of differentially expressed (DE) genes under LD conditions, particularly within clusters 1, 2, and 5. Another example is the homeobox motifs, enriched in both the LD and DD treatments, as well as around genes from clusters 1, 3, 4 and 5. Homeobox factors, in particular members of the NK homeobox gene family, with motifs enriched in both treatments, have been shown to contribute to rhythmic regulatory processes in mammals [25]. Within cluster 5, 16 cyclic genes, including NvClock (a core component of the circadian clock machinery), contain the binding motif of C/EBP in their accessible promoter region (Fig 3C and S3 Table). C/EBP acts as an enhancer of promoter activation [26,27], and its association to NvClock promoter was predicted previously [17]. Functional analysis of these 16 genes reveals that they are related to the GO term’s developmental growth (GO:0048589) and cellular response to hormone stimulus (GO:0032870) (Fig 3D) [28]. These results should be treated with skeptical eyes, as it cannot be excluded that other, unidentified, TFs might be involved in circadian rhythm regulation in N. vectensis.

Promoter accessibility is essential for gene expression, so the proportion of promoters in the accessible chromatin loci is non-random. Many expression patterns found in our ATAC-seq data were visible in the RNA-seq data as well. For example, NvClock exhibits an expression peak from late day to early night (CT9-CT13), which overlaps our ATAC-seq results from CT13 that shows its promoter to be accessible (Fig 3C). In contrast, two cryptochromes (key components of the biological clock mechanism), NvCry1 and NvCry2, exhibit transcriptional rhythms, with peak expression during the day (CT4-CT11 for NvCry1 and CT0-CT4 for NvCry2) [15], prior to the ATAC-seq sampling at CT13. Concordantly, we did not find their promoters to be accessible in either treatment. Overall, the patterns observed for NvClock, NvCry1, and NvCry2 are aligned with the correlation between chromatin accessibility measured by ATAC-seq and expression profiles measured by RNA-seq (Fig 4 and S4 Table). Furthermore, the RNA data identified four minicollagen genes with strong rhythmicity, while three of these four gene promoters were identified as significantly accessible at CT13 (p-value < 0.05). Minicollagen is an important feature of the nematocyst structure and is expressed from the early stages of nematocyst morphogenesis until capsule maturation [29]. By identifying enriched TF motifs within the peak sequences, potential gene regulators can be revealed. Within the peaks, at these gene promoters, we have identified motifs for C/EBP, Sox1, Pax-4, and Pax-6, all of which have been shown to act as clock-controlled gene regulators in mammals [27].

Enhancers are distinct genomic regions containing binding site sequences for TFs that can regulate the transcription of a target gene. Along the linear genomic DNA sequence, active enhancers, marked with H3K27ac histon modification marker, can be located at a great distance from their target genes. About 70% of H3K27ac-marked enhancers in mammals are active and positively affect transcription in vivo [30,31]. Comparing our ATAC-seq profiles with previously published H3K27ac ChIP-seq [32] revealed that ~50% of ATAC-seq peaks overlap with H3K27ac sites in both treatments (S5 Table), indicating their potential enhancer activity. Out of ~5000 previously ChIP-seq-predicted enhancer elements in N. vectensis during different early development stages, our analysis identified 259 LD-treated and 333 DD-treated enhancers that overlap with the H3K27ac histone mark, and a total of 174 enhancers shared between the two treatments (S3 Fig).

Adult N. vectensis were kept in a plastic container filled with two liters of artificial seawater at a salinity of 12 PSU (Nematostella medium), under natural light and at a constant temperature of 18°C. Between 50 and 100 individuals were kept in each container in a recirculating water system. Animals were fed 5 times a week with freshly hatched brine shrimp (Artemia nauplii).

Female N. vectensis were incubated under two different light regimens: LD (12h light:12h dark) or DD (constant darkness), for 45 hours at a constant temperature of 18°C. Biological duplicates were sampled every 8 hours at CT13, CT21, CT29, CT37 and CT45 from both conditions and were processed as described in the “ATAC-seq nuclear isolation and library preparation” section. We did not sequence sample DD CT29 due to technical problems. These time points were chosen based on gene enrichment data from RNA-seq experiments previously published. The ATAC-seq sampling interval was due to the timing requirements of the ATAC-seq protocol (approximately 6–7 hours).

Nuclei were isolated from adult Nematostella that were incubated in different lighting treatments. From each sample, tissue was suspended in 500 μL PBS-NAC 2% (N-acetyl-cysteine, sigma) by pipetting in a 1.5 mL tube [41]. The suspension was centrifuged at 1500 xg for 5 minutes, at 4°C. The pellet was re-suspended in 500 μL PBS and cells were counted. 400,000 cells were then re-suspended in 500 μL PBS and centrifuged at 1500 xg for 5 minutes, at 4°C. The pellet was suspended in 50 μL of ATAC-seq lysis buffer (10mM TRIS-Cl pH 7.4, 10mM NaCl, 3mM MgCl2, 0.1% IGEPAL CA630) and centrifuged at 300 xg for 10 minutes, at 4°C. The supernatant was collected and kept in a 1.5 mL tube on ice. The pellet was re-suspended in 50 μL and centrifuged at 300 xg for 10 minutes, at 4°C. The supernatant was combined with the supernatant from the previous step. Then 9 μL of isolated nuclei were stained with DAPI to verify the isolation of intact nuclei. The isolated nuclei were then centrifuged at 1500 xg for 10 minutes, at 4°C. Immediately following this centrifuge step, the pellet was re-suspended in the transposase reaction mix (25 μL 2× TD buffer, 2.5 μL transposase (Illumina REF: 15028212) and 22.5 μL nuclease-free water). The transposition reaction was carried out for 30 minutes, at 37°C. Directly following transposition, the sample was purified using an Invitrogen PureLink PCR purification kit (REF: K310001). Following purification, library fragments were amplified using 1× NEBnext PCR master mix (#M0541S) and 1.25 μM of custom Nextera PCR primers, forward and reverse, using the following PCR conditions: 72°C for 5 minutes, 98°C for 30 seconds and a variable number of cycles as needed (we added 4–9 cycles) at 98°C for 10 seconds, 63°C for 30 seconds and 72°C for 1 minute. To reduce GC and size bias in our PCR, we monitored the PCR reactions using qPCR to stop amplification before saturation. To do this, we amplified the full libraries for 5 cycles, after which we took a 4-μl aliquot of the PCR reaction and added 6 μl of the PCR cocktail with Sybr Green (Promega, REF: A6001), at a final concentration of 0.6×. We ran this reaction for 20 cycles to determine the additional number of cycles needed for the remaining 46-μl reaction. The libraries were purified using Agencourt AMPure XP beads (cat. No. 63881) and analyzed on a TapeStation. Primers used to amplify ATAC-seq libraries see S6 Table (Note: Ad1_noMX is a global primer used in all libraries).

Samples of whole N. vectensis were prepared using single-end 50bp reads from a single Illumina HiSeq run. Treatments were run on one lane of Illumina HiSeq2000. On average, ~50 million single-end reads were obtained for each sample.

Sequenced reads were aligned to the Nemve1 Nematostella vectensis genome using bowtie [42]. Only unique mapped reads were used. Peaks were called by applying MACS2 [43] with the following parameters: -g 450000000—nomodel—extsize 75—shift -30. TF-binding motifs enrichment were identified within the peaks using scripts within HOMER [21]: findMotifsGenome.pl and annotatePeaks.pl were used with default parameters and the Nemve1 genome was used as background [44].

To compare chromatin accessibility with circadian patterns in gene expression, we evaluated a set of 180 transcripts that were previously shown to exhibit a diel expression pattern in N.vectensis [15]. The genes had been sorted into 5 clusters with similar temporal expression patterns, using a K-means clustering, implemented in MatLab as described by Oren et al., [45]. Pearson correlation tests were apllied to assess correlation between gene expression and DNA accessibility. We have compared all genes’ expressions to promoter accessibility within a single time point (for example, gene expression at CT13 was compared to promoter accessibility at CT13). In addition, we compared the total RNA-seq data from our 139 gene list to the total ATAC-seq data at all time points (see S2 Table). To identify possible distal enhancer sites in the ATAC-seq data, ChIP-seq data was used from a previous study of histone markers in N.vectensis [32]. Reads were downloaded from NCBI GEO (accession number: GSE46488↗) and aligned to the genome using bowtie2. Only uniquely mapped reads were used. Peaks were called by applying MACS2 [43] with default parameters. Further analysis was performed using the BamTools and BEDTools suites [46,47]. Gene promoters found within treatment specific peaks were defined and subsequently, analyzed for enriched GO terms using the metascape suit [28].