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.