2.4 RNA-seq
We terminated the diet treatment at 17 days and each individual was sacrificed. The interopercle, interopercle-mandibular ligament, and retroarticular (i.e., IOP-RA complex, Figure 1) was dissected from each fish for RNA-seq (left-side) or qPCR (right-side). It is to be noted that this complex is a mix of tissue types (e.g bone, ligament, epithelial tissue). For RNA-seq, 6 samples from each treatment and species were stored in Trizol (Invitrogen) at -80°C, homogenized using a Next Advance Bullet Blender and 5 UFO beads each, and processed using the phenol-chloroform method of RNA extraction, but did not undergo cDNA conversion. We standardized each sample to 500ng total RNA in 50uL and produced libraries using the TruSeq Stranded mRNA Library Prep Kit (Illumina). Any remaining RNA not used for RNA-seq was stored at -80°C. Libraries were sequenced at the University of Massachusetts Medical School Deep Sequencing Core with a HiSeq 4000 with 50 x 50 paired end reads.
Raw reads from RNA-Seq were assessed using FastQC (Andrews, 2010) and ends were trimmed accordingly using cutadapt (Martin, 2011). Cleaned reads were -mapped against the Maylandia zebra genome version UMD2a (Yates et al., 2020) with Bowtie2 (Langmead and Salzberg, 2012) and a matrix of read counts was generated from the alignments with HT-Seq-count (Anders et al., 2015).
We used edgeR (Chen et al., 2014; Robinson et al., 2010) to identify differentially expressed genes between treatments, as well as those that showed an additive effect between species and environment. Results were groundtruthed by visually comparing normalized counts (counts-per-million; cpm) among treatments. Gene ontology (GO) terms were assigned using biomaRt (Durinck et al., 2005, 2009) via the biomartr package (Drost and Paszkowski, 2017) and enrichment analysis was conducted via topGO (Alexa and Rahnenfuhrer) in the R environment (R Core Team, 2021) using the weight01 algorithm and Fisher’s exact test.
Tissues from individuals not used for RNA-seq, were stored in Trizol at -80°C, and homogenized as previously described. We followed the phenol-chloroform RNA extraction method and converted RNA to cDNA for gene expression validation purposes. These samples were standardized to an RNA concentration of 70ng/uL.
qPCR was used to measure gene expression of genes found to be differentially expressed and/or differentially accessible from the RNA-seq and ATAC-seq datasets, and those related to bone development/homeostasis using the comparative CT method. qPCR primer sequences are listed in the supplementary table (Table S1).