2.5 ATAC-seq
We designed an ATAC-seq protocol optimized for bony/ligamentous fish
tissues (adapted from Buenrostro et al., 2013; Corces et al., 2017). The
ATAC-seq diet treatment was terminated 28 days after foraging treatment
began (Figure 1E). Similar to the RNA-seq experiment, each animal had
the IOP-RA complex from the left side of the face removed after being
euthanized. Briefly, each sample was placed in 3% collagenase II in 5%
FBS/DMEM for cell collection for 2 hours. To ensure we collected cells
of the appropriate size, we filtered the cells through a 70um strainer.
Cell quality and count was confirmed by inverted light microscopy. We
collected up to 500,000 cells for each sample. Cells were then lysed to
isolate nuclei and underwent a transposition reaction to cut chromatin,
and the resulting DNA fragments were purified using a Qiagen MinElute
Cleanup Kit. We constructed libraries from the transposed DNA and
performed double-sided bead purification to remove large
(>1000bp) and small (e.g primer dimer) DNA fragments. The
detailed protocol is attached in the supplementary information.
Libraries were sequenced in the same manner as RNA-seq libraries.
Raw reads were again aligned against the Maylandia zebra genome using
bowtie2, and data were converted to appropriate formats for downstream
analyses using samtools (Li et al., 2009) and bedops (Neph et al.,
2012), with parallelization enabled by gnu parallel (Tange, 2018). Peaks
were called twice using Genrich (https://github.com/jsh58/Genrich): once
where we combined replicates from each treatment to provide more power
for peak calling for occupancy analysis and once where peaks were called
for each replicate individually to improve statistical power for
downstream differential accessibility analysis. In both instances, flags
were set to filter mitochondrial reads and PCR duplicates before peak
calling.
Occupancy and differential accessibility were assessed using DiffBind
(Stark and Brown), with the latter analysis calling the edgeR algorithm.
Peaks that were enriched for occupancy in one treatment or
differentially accessible between any two treatments were annotated by
intersecting the genomic coordinates of the peaks with the Maylandia
zebra gtf file using bedtools (Quinlan and Hall, 2010), with the
requirement that 30% of the peak must overlap with the annotated
feature. Peaks identified from the occupancy analysis were then filtered
down to those that were enriched in a single treatment only, and
examined for enrichment of GO terms using biomaRt (Drost and Paszkowski,
2017) and topGO (Alexa and Rahnenfuhrer) with the latter invoking the
classic algorithm and Fisher’s exact test.