3.4 Species-specific differences in chromatin accessibility are influenced by foraging conditions
Areas of the genome that contain open chromatin are more accessible to transcriptional machinery and therefore able to increase gene expression, whereas closed chromatin sites are less accessible for transcription. The differences in accessibility (e.g., between species or environments) are considered differentially accessible (DA), and may be due to either genetic (e.g., deletion of a TF binding site) or epigenetic (e.g., methylation changes) processes. When comparing species across both environments, we identified 10770 areas of accessible chromatin, and of these 297 were DA. Note that many genes have contain more than one accessible chromatin peak (e.g., Figure 6B,D). The number of DA loci was therefore considerably less than the number DE loci, which may reflect differences in the timing and/or nature of each experiment. In addition, we did not observe a marked bias in one environment versus the other in terms of the number of differentially accessible genes (DAGs), with 114 DAGs identified in pelagic fishes and 157 in benthic fishes (Table 1). These data suggest that species differences in DA are not being driven by one environment at the time when tissues were collected.
The overlap between RNA-seq and ATAC-seq datasets implicates loci acting over extended periods (i.e., at both timepoints) following the onset of foraging trials (Figure 6A, Table S3). In all, we identified 15 genes that were both DE and DA in fishes reared in the pelagic foraging environment, and 11 from the benthic environment. We also identified 13 genes that were DE and DA in both environmental conditions, which suggest that their expression is robust to differences in the environment The direction of DE and DA across all genes was generally consistent (Table S3). In particular, out of 39 loci, the polarity of DE and DA was similar in 34. Differences in the other 5 could be due to DA being associated with the binding of a repressive transcription factor. Alternatively, given that RNA-seq and ATAC-seq experiments were performed at different time points, it is also possible that expression of these factors may oscillate over time.
A few of the genes in these overlapping datasets have been implicated in craniofacial development, including KIAA0586 (also known astalpid3 ), which is essential for primary cilia formation and Hedgehog (Hh) signaling (Schock et al., 2016). The identification ofKIAA0586 is especially notable given previous work from our lab, which has demonstrated important roles for genes associated with the primary cilium-Hedgehog molecular mechanism in species-specific shaping and plasticity of craniofacial bones in cichlids and zebrafish (Hu & Albertson, 2014; Navon et al., 2020; Gilbert et al., 2021; Zogbaum et al., 2021). While KIAA0586 is well-studied in the context of early craniofacial patterning (reviewer by, Schock et al., 2016), roles at latter stages, including growth and plasticity, have not been explored. Taken together, multiple independent experiments in the cichlid system support the thesis that the primary cilium-Hedgehog “signal transduction machinery” is an important and evolvable mechanism for shaping the craniofacial skeleton. As opposed toKIAA0586 , other genes in this dataset are largely new to the field of craniofacial biology, but implicate biological processes important in bone patterning, formation, growth and homeostasis, including chromatin remodeling (e.g actr6 [Yoshida et al., 2010]), cell signaling (e.g asb5 [Yoshioka et al., 2006],etv4 [Mao et al., 2009]), and cell growth (e.g impdh1[Chang et al., 2015]).
Determining whether or not these factors are DA due to genetic or epigenetic factors would be a fruitful line of future study. Resequencing a panel of cichlids around DA peaks would allow us to assess whether any indels or SNPs might underlie differences noted here. In addition, the co-occurrence of DA peaks and CpG islands would suggest that differential DNA methylation might be driving differences in expression. For example, the DA peak associated with KIAA0586expression is in intron 7-8 (Table S3), which is large and contains several predicted CpG islands, although none overlap with the DA peak (Figure 6C).