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