\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[ngerman,english]{babel}



\usepackage{float}








\begin{document}

\title{The transcriptional state and chromatin landscape of cichlid jaw shape
variation across species and environments}



\author[1]{Emily Tetrault}%
\author[1]{John Swenson}%
\author[1]{Ben Aaronson}%
\author[1]{Chelsea  Marcho}%
\author[2]{Craig Albertson}%
\affil[1]{University of Massachusetts Amherst}%
\affil[2]{University of Massachusetts}%


\vspace{-1em}



  
  \date{May 24, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
Adaptive phenotypes are shaped by a combination of genetic and
environmental forces, and while the literature is rich with studies
focusing on either genetics or environment contributions, those that
consider both are rare. Here we utilize the cichlid oral jaw apparatus
to fill this knowledge-gap. First, we employed RNA-seq in bony and
ligamentous tissues important for jaw opening to identify differentially
expressed genes between species and across foraging environments. Our
foraging treatments were designed to force animals to employ either
suction or biting/scraping, which broadly mimic pelagic or benthic modes
of feeding. We found a large number of differentially expressed genes
between species, and while we identified relatively few differences
between environments, species differences were far more pronounced when
reared in the pelagic versus benthic environment. Further, these data
carried the signature of genetic assimilation, and implicated cell cycle
regulation in shaping the jaw across species and environments. Next, we
repeated the foraging experiment and performed ATAC-seq procedures on
nuclei harvested from the same tissues. Cross-referencing results from
both analyses revealed subsets of genes that were both differentially
expressed and differentially accessible in either the pelagic (n=15) or
the benthic environment (n=11), as well as loci where differences were
robust to foraging environment (n=13). All in all, these data provide
novel insights into the epigenetic, genetic, and cellular bases of
species- and environment-specific bone shapes, as well as the evolution
of phenotypic plasticity in this iconic model system.%
\end{abstract}%



\sloppy


\textbf{1 INTRODUCTION}

A major ongoing challenge in biological research is to understand the
origin and maintenance of biodiversity, with broad implications in
conservation, ecology and evolutionary biology. Traditionally, these
endeavors have involved characterizing the forces and mechanisms
operating above the organismal level (e.g., selection, environmental
change {[}Schluter 2000; Callaghan et al., 2004; Burns et al., 2009;
Siepielski et al., 2014{]}) or within the organism (e.g., genetic and
developmental mechanisms {[}Kawajiri et al., 2014; Hohenlohe, 2014;
Margres et al., 2015{]}). Understanding the intersection of extrinsic
and intrinsic forces (Mccairns \& Bernatchez, 2012; van Heerwaarden \&
Sgr\selectlanguage{ngerman}ò, 2017; Laitinen \& Nikoloski, 2019; Levis et al., 2020) holds
significant potential to advance the field.

African cichlids are a hyperdiverse group of fishes that have long been
used as an evolutionary model (Kocher, 2004; Seehausen, 2006), and have
been especially useful in revealing both the genetic and environmental
factors that contribute to biodiversity (McKaye et al., 1984; Sturmbauer
\& Meyer, 1992; Genner et al., 2004; Ding et al., 2015; Malinsky et al.,
2018). In Africa, approximately 2000 cichlid species have arisen over
the past \textasciitilde{}5 million years, which is unparalleled
compared to the speciation rates of other vertebrates (Seehausen, 2006).
Moreover, cichlid diversity is pronounced across several phenotypic
axes, including coloration (Seehausen et al., 1999; Maan et al., 2006;
Salzburger, 2009), activity levels (Lloyd et al., 2021), as well as
reproductive and foraging behaviors (Balshine-Earn \& Earn, 1998; Genner
et al., 1999; Lopez-Fernandez et al., 2014). Variation in feeding
architecture, which relates to the foraging niche exploited by each
species/population, is another critical axis of cichlid diversity (e.g.,
Cooper et al, 2010). Cichlid craniofacial variation is largely
continuous, but there are also examples of extreme or discontinuous
variants (reviewed by Powder and Albertson, 2016). In general terms,
cichlids partition their foraging niche along a benthic-pelagic
ecomorphological axis, with concomitant shifts in foraging anatomy
(Young et al 2009; Cooper et al. 2010; Conith and Albertson, 2021). For
instance, species inclined toward a benthic mode of feeding tend to have
steeply descending facial profiles, small eyes positioned toward the top
of their heads, and short, robust oral jaws with closely-spaced,
multicuspid teeth optimal for biting and scraping (e.g., Figure 1A,B).
On the opposite end of this spectrum, pelagic feeders tend to possess
longer, streamlined heads, large eyes, and long, up-turned oral jaws
with large, widely spaced teeth optimal for suction/ram feeding (e.g.,
Figure 1C,D) (Albertson et al., 2003; Cooper et al., 2010).

Significant efforts over the past 20 years have focused on
characterizing the genetic basis of cichlid craniofacial variation
(e.g., Albertson et al., 2005; Roberts et al., 2011; Powder et al.,
2014; Hu \& Albertson, 2017; Singh et al., 2017; DeLorenzo et al.,
2022). In addition, cichlids have long been a model of phenotypic
plasticity (Meyer 1987; Wimberger 1991; Huysseune 1995;
Machado-Schiaffino et al., 2014; Schneider et al., 2014; Meuthen et al.,
2018; Navon et al., 2020), which is defined as the ability of a single
genotype to produce a range of phenotypes in response to environmental
inputs. Plasticity is critical for organismal survival in an era of
rapid environmental change (Willis et al., 2008; Sih et al., 2011;
Gugger et al., 2015; Karasz et al., 2022; Morgan et al., 2022). It can
also influence the direction and/or speed of future evolutionary change
by exposing new phenotypic and genetic variants to natural selection
(Ledon-Rettig et al., 2010 ProcB; McGuigan et al., 2011 evol; Landy et
al., 2020 PNAS; Campbell et al., 2021). In spite of its importance
across a range of biological disciplines, there are many outstanding
questions about plasticity, including its genetic basis and evolutionary
potential (Gibert 2017). Plasticity is well documented in cichlids
across a range of morphological traits including full body,
craniofacial, oral jaw and pharyngeal jaw shapes (Huysseune, 1995;
Muschick et al., 2011; Gunter et al., 2013; Parsons et al., 2014; Navon
et al., 2020). A notable theme that has come from these data is that
closely related species can differ in either their magnitude or pattern
of plasticity in response to the same stimulus (Parsons et al., 2014;
Navon et al., 2020), suggesting that plasticity itself is an evolvable
trait. If true, then plasticity must also have an explicit genetic basis
(Kuttner et al., 2014; Lafuente et al., 2018; Diouf et al., 2020);
however, understanding plasticity at this level has proven challenging
(Gibert 2017).

Previous efforts in our lab have sought to describe the genetic basis of
plasticity, and have described roles for Wnt (Parsons et al., 2014) and
Hh (Hu \& Albertson 2017; Navon et al., 2020) signaling, respectively.
In addition, QTL analyses in cichlids have demonstrated the critical
importance of the environment in determining the genotype-phenotype
(G-P) relationship. Specifically, the genetic basis of variation in
multiple hard and soft tissue traits was shown to depend, almost
entirely, on the foraging environment in which the animals were reared
(Parsons et al., 2016; Zogbaum 2021). Such genetic mapping studies led
to the discovery of \emph{crocc2} as an environmentally-dependent
regulator of jaw shape (Gilbert \& Tetrault et al., 2021). \emph{Ciliary
rootlet coiled-coil 2} (\emph{crocc2} ) encodes a protein that is a
major structural component of the primary cilium's rootlet (Yang et al.,
2002). Primary cilia are important mechanosensors that help cells sense
and respond to environmental stimuli, but roles of the rootlet in
mechanosensing are less clear (Styczynska-Soczka \& Jarman, 2015).
Notably, this gene was only implicated in regulating cichlid jaw shape
in the mechanically demanding benthic/biting environment (Parsons et al
2016; Gilbert \& Tetrault et al., 2021), and functional analyses in
zebrafish showed that mutations in \emph{crocc2} led to degeneration of
cilia, decreased mechanosensing abilities, dysmorphic bone shapes, and
mis-regulation of gene networks in bone tissue (Gilbert \& Tetrault et
al., 2021). Together, this incipient literature has implicated a small
handful of genes that contribute to mechanosensitive signal transduction
pathways (e.g., Hh) and structural components of the cell (e.g.,
rootlets) in the evolution and plasticity of cichlid bone shape. Here we
seek to advance this research program by taking a genome-wide approach.

In particular, to address the question of genetic and epigenetic control
of plasticity in the cichlid feeding apparatus, we utilize two
complementary methods of assessing transcriptional output: RNA-seq to
analyze gene expression, and ATAC-seq to assess chromatin accessibility.
We focused on an important functional complex - i.e., the
interopercle-retroarticular (IOP-RA) complex - which (1) is part of the
opercle 4-bar linkage chain, (2) helps to drive lower jaw depression,
(3) is comprised of hard and soft tissues, (4) varies among Malawi
cichlids in a manner that predicts foraging mode/habitat, and (5) has
been shown to be plastic in previous research (Figure 1; Hu and
Albertson 2014; 2017; Navon et al 2021). Our goals are to identify genes
that are both differentially expressed (DE) and differentially
accessible (DA) between species and environments.

\textbf{2 MATERIALS AND METHODS}

\textbf{2.1 Fish Husbandry}

Juvenile-stage cichlids were purchased and housed in 40-gallon glass
aquaria at \textasciitilde{}28°C on a 14hr light/10hr dark cycle. Each
aquarium was part of a recirculating system, with automated daily water
changes and chemical dosing to ensure consistent water quality.
Juvenile-stages were chosen to mitigate the confounding effects of
dominance. It is also the period of development when all bony elements
have formed but animals are still actively growing (Fujimura and Okada,
2007; 2008a; 2008b). Cichlid husbandry follows a protocol approved by
the institutional animal care and use committee at the University of
Massachusetts.

\textbf{2.2 Experimental design}

To determine how alternative feeding regimes affected gene expression
and the \emph{cis} -regulatory system, we used a benthic diet or pelagic
diet to impose a power or speed demand, respectively, on the oral jaw.
This experimental design was used in two cichlid species,\emph{Tropheops
sp.} ``red cheek'' (TRC) and \emph{Maylandia zebra} (MZ) (Figure 1A-D)
to compare species differences in plasticity. Diet treatments for
RNA-seq ran for 17 days, and 28 days for ATAC-seq (Figure 1E). The
combination of approaches allowed us to examine both differential
expression and differential chromatin accessibility between species and
environments. That the experiments were performed at two different
time-points following the onset of foraging challenges, allowed us to
identify loci with effects over extended periods of time.

\textbf{2.3 Foraging challenges}

We separated each species into two diet treatments that would impose
different functional demands on the oral jaw. Each 40-gallon tank held
5-10 individuals of a single species that were provided either a benthic
diet or pelagic diet. For benthic feeders, we mixed ground high-quality
cichlid flake food and freeze dried brine shrimp into 1.5\% food-grade
agar to create a food ``paste'' that was spread onto store bought lava
rocks. We allowed the rocks to dry overnight and 2 rocks were placed
into each benthic tank each day. Pelagic feeders were given ground up
flake food and live brine shrimp daily. Animals were given 1 week to
train on each diet before the start of the experiment (Figure 1E).

\textbf{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).

\textbf{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
(\textgreater{}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.

\textbf{3 RESULTS \& DISCUSSION}

\textbf{3.1 Extensive species differences in gene expression}

The cichlid genome has over 25000 protein-coding genes (Conte et al.,
2019). Of these, 17525 were expressed in focal tissues, and 5318 were
differentially expressed (DE) between our two species, across both
foraging conditions, with 2667 upregulated in MZ (relative to TRC - red
in Figure 2A) and 2651 upregulated in TRC (relative to MZ - blue in
Figure 2A) (Table 1). Cluster analysis using the top 500 most
differentially expressed genes (DEGs) illustrates a clear separation of
species (Figure 2B), but differences between foraging environments are
less apparent. While MZ shows some separation by environment, TRC does
not (more on this below).

Gene ontology (GO) analyses were performed on DEGs, using the annotated
MZ genome (ensembl.org), to determine what biological processes were
enriched (Figure 2C). Whereas processes related to cell cycle regulation
were among the most enriched in MZ, a diversity of other processes were
enriched in TRC. These data suggest that there may be divergent modes of
bone growth operating in these species at the time when tissues were
collected.

\textbf{3.2 Differential gene expression is rare between foraging
environments within species}

Craniofacial plasticity is well documented in cichlids (Meyer 1987;
Wimberger 1991; van Snick Gray \& Stauffer, 2004; Schartau et al., 2009;
Navon et al., 2020), and, to an extent, skeletal plasticity has been
associated with changes in gene expression (Parsons et al., 2014;
Parsons et al., 2016; Navon et al., 2020). However, most work on the
genetics of craniofacial plasticity in cichlids has focused on the lower
pharyngeal jaw (Gunter et al., 2013; Schneider et al., 2014). Our work
complements this body of literature by focusing on skeletal and
soft-tissue elements critical for lower oral jaw depression (i.e.,
IOP-RA complex, Figure 1). While 38 DEGs were detected between foraging
environments in MZ, none were detected for TRC (Table 1, Figure 3A,B),
suggesting that MZ is more plastic than TRC, and/or plasticity arises
earlier in MZ compared to TRC. Of the nearly forty DEGs within MZ, 25
were upregulated in the pelagic environment, whereas 13 were upregulated
in the benthic environment (Table 1, Figure 3B), which is consistent
with previous data showing greater rates of bone deposition in MZ
exposed to a pelagic environment (Navon et al 2020).

Within MZ, genes upregulated in the pelagic environment (relative to the
benthic environment) seemed to be largely associated with proliferation,
which was reflected by the GO analysis (Figure 3C). For
example,\emph{ccna2} is a cyclin that activates \emph{cdk2} and promotes
cell cycle progression through both G1/S and G2/M phases (Pagano et al.,
1992), and \emph{cks1b} slows the progression of G1/S and can block
entry to M phase (Westbrook et al 2007). In addition,
\emph{cdc20}regulates metaphase-anaphase transition during mitosis via
activation of APC, which targets proteins for degradation (Visintin et
al., 1997; Yu 2002). Finally, \emph{cdca5} , also known as sororin,
plays an important role in cell proliferation (Fu et al., 2020), and
more specifically in the binding of Cohesin to chromatin during cell
division (Rankin et al 2005; Schmitz et al. 2007). Data presented here
strongly suggests that cell proliferation is critical in mediating
skeletal plasticity in this species.

For MZ in the benthic environment, the only GO term enriched for
upregulated genes (relative to pelagic) was stress response. In
addition, a diversity of other processes are implicated by the other
genes, including skeletal muscle changes in response to stimuli
(e.g.,\emph{arrdc3b} , Gordon et al., 2019; \emph{myoglobin,} Beyer et
al., 1984), and osteoblast proliferation via regulation of
\emph{cyclins}(e.g., \emph{per2} , Fu et al., 2005). Our previous work
suggested that benthic foraging might be a non-preferred environment for
MZ, at least in terms of environmentally-stimulated bone growth (Navon
et al., 2020). Here, this assertion is supported by the observation that
genes upregulated in the pelagic environment all seemed to contribute to
the same biological process -- e.g., cell proliferation -- with known
roles in bone formation/growth (Capecchi et al., 2018; Shekhar et al.,
2019; Gao et al., 2020; Du et al., 2021). Alternatively, genes
upregulated in the benthic environment were involved in a diversity of
processes, including stress response.

Craniofacial plasticity has been noted in both MZ and TRC (Parsons et al
2014; Navon et al., 2020), and so the lack of more extensive DEGs
between foraging treatments, especially in TRC, was somewhat surprising.
Previous studies have used time points measured in months (Schneider et
al., 2014) or years (Gunter et al., 2013) to examine DEGs between
foraging environments in the cichlid lower pharyngeal jaw, and so it is
possible that we did not allow enough time for a plastic response to
manifest. However, we have previously demonstrated a measurable plastic
response in bone matrix deposition after 5 weeks (Navon et al., 2020),
which should be underlain by an earlier transcriptional response.
Indeed, mechanical load has been shown to induce gene expression changes
in bone cells in a matter of hours (Raab-Cullen et al., 1994; Mantila
Roosa et al,. 2011; Govey et al., 2015; Kelly et al., 2016). It is
therefore also possible that a more robust plastic response in gene
expression might occur earlier than the time point sampled here. In
addition, FDR is a stringent metric, and it is possible that
biologically relevant changes in expression have occurred but are not
detected by standard pipelines. In this regard, examining genes with
high fold changes, but FDR-values \textgreater{}0.05, might prove
fruitful. For example, a transcript that was upregulated in benthic TRC
(logFC = 5.44) corresponds to \emph{receptor transporting protein
3}(\emph{rtp3} ), which has been linked to human femoral cortical
thickness and buckling ratio, as well as hip fractures (Zhao et al.,
2010), and another upregulated benthic gene (logFC = 4.83),
\emph{poly(ADP-ribose) polymerase 14} (\emph{parp14} ), has been shown
to regulate cell-cycle progression via Cyclin D1 (O'Connor et al.,
2021). Thus, viable candidates for craniofacial bone plasticity may be
found just under the threshold set by RNA-seq protocols.

Finally, we stress that a relatively low number of DEGs within species
does not preclude more general roles for the environment in determining
species-specific bone shapes. As a next step we therefore assessed the
effect of foraging environments on DE between species.

\textbf{3.3 The pelagic environment drives species-specific differences
in gene expression and reveals signatures of genetic assimilation.}

We have shown previously that foraging conditions can have a marked
impact on the genotype-phenotype map. Specifically, quantitative trait
loci (QTL) for the same trait map to largely distinct regions of the
genome when animals are reared under alternate benthic/pelagic foraging
conditions (Parson et al., 2016; Zogbaum et al., 2021). We therefore
examined expression differences between species within each environment,
and documented a marked imbalance in DEGs. Specifically, when only
considering animals exposed to pelagic conditions, we found over 3500
DEGs between species, whereas fewer than 1000 DEGs were detected between
species when only comparing animals reared under benthic conditions
(Table 1, Figure 4A-C). Additionally, when comparing
environment-specific DEGs to the total dataset (i.e., combining both
environments), we found that over 2500 genes from pelagic animals were
represented in the global comparison, whereas only 155 genes from
benthic animals overlapped between datasets (Table 1, Figure 4C). These
data underscore the importance of environmental context in determining
the genetic basis of species-specific bone shapes (Parsons et al 2016;
Zogbaum et al, 2021). More specifically, they suggest that the pelagic
foraging environment is driving species differences in gene expression
within the IOP-RA functional complex.

This trend is drawn out when comparing genes from an additive model
(S+E), whereby DEGs were detected at the level of both species and
foraging environment (Table S2; Figure 5A). When illustrated in a
heatmap, these data support the assertion that species differences in
gene expression are driven by the pelagic environment, and reveal
patterns consistent with either genetic accommodation or assimilation.
Genetic assimilation is a mechanism by which plasticity is lost over
evolutionary time as genetic variation that facilitated plasticity in an
ancestral population becomes fixed as descendent populations adapt to a
specific environment (reviewed by Pigliucci et al., 2006). If we assume
that plasticity is ancestral, evidence for genetic assimilation is
apparent in several gene clusters (denoted by pink dots, Figure 5A),
whereby TRC expression levels are indistinguishable between foraging
environments and match those of benthic MZ. Consistent with previous
data many of the DEGs identified by this model contribute to cell cycle
regulation -- e.g., \emph{cdc20} , \emph{cdca5} , \emph{cdca8}
,\emph{ccne2} , \emph{ccnb1} , \emph{ccnb2} , \emph{ccnf} . A list of
all the DEGs in this model can be found in Table S2.

Alternative to genetic assimilation is genetic accommodation, or an
increase in genetic plasticity over evolutionary time. We cannot rule
out that this is the case, as it is possible that plasticity has been
enhanced beyond the ancestral condition in MZ. Regardless, the main
conclusion to be drawn from these data is that the evolution of
plasticity in this system may be traced to divergent patterns of gene
expression associated with cell cycle regulation.

We next performed GO analyses for DEGs between species in each foraging
environment. When considering animals reared in the pelagic foraging
environment, GO analysis revealed a diversity of biological processes;
however, those associated with cell division were among the most
enriched in MZ, whereas translation and cell differentiation were among
the most enriched processes in TRC (Figure 4D). For animals reared in
the benthic environment, comparatively fewer biological processes were
enriched in general, consistent with fewer DEGs being identified.
Similar to pelagic fishes, this analysis found enrichment of cell cycle
genes in MZ, and cell differentiation in TRC (Figure 4E). While a
greater number of DEGs, contributing to a larger number of biological
processes, underlie species-specific differences in the pelagic
environment, there are notable consistencies between environments.
Specifically, an increase in cell number seems to be important for
skeletal growth in MZ, whereas cell differentiation may shape growth in
TRC at the time points when tissues were collected.

Unsurprisingly, enriched GO terms for the additive model are similar to
those for MZ in the pelagic environment, and include cell cycle, cell
division, and chromosome segregation (Figure 5B). In addition, this
analysis found enrichment of cytoskeleton organization, which is
critical to many cellular functions relevant to bone formation and
plasticity, including mechanotransduction (Gunst \& Zhang, 2008), and
primary cilia formation (Mirvis et al., 2018), which we and others have
found to be necessary for load-induced bone formation (Chen et al.,
2016; Moore et al., 2019; Gilbert \& Tetrault et al., 2021).

\textbf{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 \emph{KIAA0586} (also known
as\emph{talpid3} ), which is essential for primary cilia formation and
Hedgehog (Hh) signaling (Schock et al., 2016). The identification
of\emph{KIAA0586} 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 \emph{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
to\emph{KIAA0586} , 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 \emph{actr6} {[}Yoshida et al.,
2010{]}), cell signaling (e.g \emph{asb5} {[}Yoshioka et al.,
2006{]},\emph{etv4} {[}Mao et al., 2009{]}), and cell growth (e.g
\emph{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
\emph{KIAA0586}expression 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).

\textbf{3.5 Validation of select candidate genes by qPCR}

We validated a subset of genes from this analysis with qPCR, focusing on
those that overlapped in the RNA-seq/ATAC-seq datasets (Figure S1A-E).
Trend across transcript counts (i.e., counts per million, cpm) and qPCR
expression were generally equivalent, but we note that qPCR picked up
significant differences in gene expression that did not meet the
threshold for RNA-seq significance. For example, \emph{actr6,} which
plays important roles in heterochromatin formation in
yeast,\emph{Drosophila} and vertebrates (Ohfuchi et al., 2006), was DE
expressed between benthic and pelagic TRC according to qPCR, but not by
cpm (Figure S1A). In addition, \emph{gnmt} , a methyltransferase
involved in the methionine pathway, which is linked to bone density
(Wang et al., 2011; Wang et al., 2011; Ables et al., 2012; Vijayan et
al., 2014), exhibits an expression pattern that is similar to cpm, but
it is also plasticity expressed among Maylandia reared in different
environments (Figure S1D). We also validated the expression of
\emph{cdc20} , a gene that is involved in cell-cycle regulation
(Visintin et al., 1997; Yu 2002), and was significantly DE between
benthic and pelagic MZ. Expression of this gene via qPCR did not quite
reach significance at the 0.05 alpha level, but exhibited a lot of
variation across samples, and showed a similar pattern that was trending
toward significance (p=0.17, Figure S1E).

\textbf{4 CONCLUSIONS}

A number of general themes emerge from these data. First, they provide
clear support for the hypothesis that foraging environment influences
the genotype-phenotype map for craniofacial skeletal traits (Parsons et
al., 2016; Navon et al., 2021; Zogbaum et al., 2021). More specifically,
our data suggest that pelagic foraging ``drives'' species- and
environment-specific DE. This may seem counterintuitive as diets that
involve large/hard prey items are generally considered to be the more
mechanically demanding compared to small/soft food (Muschick et al.,
2011; Gunter et al., 2013; Hulsey et al., 2020). However, Navon et al
(2020) showed that in MZ bone matrix was deposited at a fast rate under
pelagic foraging conditions, and speculated that suction feeding imposes
mechanical load on the feeding apparatus as animals repeatedly open and
protrude their jaws. Our data support this assertion, and thus we
consider the foraging treatments utilized here to challenge the feeding
apparatus in two distinct ways (compared to a ``standard'' flaked food
diet); our benthic treatment was designed to impose high amplitude but
low frequency loading onto the feeding apparatus as animals scrapped
food from rocks, whereas our pelagic treatment translated to higher
frequency but lower amplitude loading as animals repeatedly protruded
their jaws to gather small food items.

Our data also detected evidence for genetic assimilation. In particular,
when considering loci that were DE between species + environments,
patterns in MZ benthic fish resembled those across TRC.
\emph{Tropheops}species, including \emph{sp.} ``red cheek'', are
generally found in a benthic environment (Ribbink et al., 1983), and may
have lost a degree of plasticity as they evolved to specialize on
benthic food items. MZ on the other hand are true generalists in the
sense that they routinely foraging from both the benthic and pelagic
zones (Ribbink et al., 1983). While plasticity has been noted in TRC
(Parsons et al., 2014; Navon et al., 2020), our data suggest that MZ may
be more plastic than TRC in that they mount a more pronounced
transcriptional response, at least at the time point analyzed in this
study.

Cell cycle regulation consistently appeared in GO analyses, describing
species differences, as well as plasticity within MZ. This implicates
cell proliferation as an important biological mechanism of species- and
environment-specific bone growth in cichlids. This observation is
notable as our previous work has implicated Hedgehog signaling in the
evolution and plasticity of the cichlid jaw, including the IOP-RA
complex (Hu \& Albertson, 2014; Parsons et al., 2016; Navon et al.,
2020). While canonical members of the Hedgehog signaling pathway were
not significantly DE or DA in this dataset (although
\emph{KIAA0586}regulates the signal, Schock et al., 2016), cell
proliferation is well-known to be regulated by this pathway (St. Jacques
et al., 1999; Tiet et al., 2006; Sun \& Deng, 2007; Zaman et al., 2019),
providing a potential cellular mechanism through which variation in
Hedgehog signaling leads to differences in bone shape among and within
cichlid species.

Finally, with these large overlapping genome-wide datasets, we were able
to narrow down thousands of DEGs to roughly two dozen that were both DE
and DA. Given that each experiment was conducted at a different time
point, this reduced dataset points to loci whose expression is important
for species divergence over extended periods of time. Among these were
genes that were both sensitive and robust to the environment. Notably,
nearly all of these genes are new to the field of bone biology, and
while some encode known effectors of well studied signaling pathways
(e.g., interleukin/Wnt, Talpid/Hh) and cell behaviors (e.g.,
Casp6/apoptosis, Impdh1b/cell-cycle), others implicate largely novel
mechanisms (e.g., Gnmt/methionine cycle). Thus, this work establishes a
robust foundation for future studies into how genotype and the
environment combine to influence bone formation, remodeling, and
evolution.

\textbf{ACKNOWLEDGEMENTS}

We would like to thank current and former members of the Albertson Lab
for helpful comments and suggestions on various aspects of this project.
This project was funded by an NIH grant awarded to R.C.A (DE026446).

\textbf{DATA ACCESSIBILITY}

All data are available on GitHub at
https://github.com/tetra22e/Genomics2022

\textbf{BENEFITS-SHARING}

Benefits from this research accrue from the sharing of our data and
results via open access publishing, and on public databases as described
above.

\textbf{AUTHOR CONTRIBUTIONS}

R.C.A. and E.M.T. designed the experiments and wrote the manuscript.
E.R.T. and C.M. designed the ATAC-seq protocol and executed all
sequencing experiments. E.R.T. and B.A. performed the qPCR validation.
E.R.T. and J.S. analyzed the RNA- and ATAC-seq data. R.C.A., E.R.T.,
J.S., and C.M. interpreted the results. J.S. and C.M. contributed to
editing the manuscript.

\textbf{REFERENCES}

Albertson, R. C., Streelman, J. T., \& Kocher, T. D. (2003). Genetic
basis of adaptive shape differences in the cichlid head. \emph{Journal
of Heredity} , \emph{94} (4), 291--301.
https://doi.org/10.1093/jhered/esg071

Albertson, R. C., Streelman, J. T., Kocher, T. D., \& Yelick, P. C.
(2005). Integration and evolution of the cichlid mandible: The molecular
basis of alternate feeding strategies. \emph{Proceedings of the National
Academy of Sciences of the United States of America} , \emph{102} (45),
16287--16292. https://doi.org/10.1073/pnas.0506649102

Alexa, A., and Rahnenfuhrer, J. Gene set enrichment analysis with topGO.
26. .

Anders, S., Pyl, P.T., and Huber, W. (2015). HTSeq---a Python framework
to work with high-throughput sequencing data. Bioinformatics \emph{31} ,
166--169. https://doi.org/10.1093/bioinformatics/btu638.

Andrews, S. (2010). FastQC: a quality control tool for high throughput
sequence data (Babraham Bioinformatics, Babraham Institute, Cambridge,
United Kingdom).

Balshine-Earn, S., \& Earn, D. J. D. (1998). On the evolutionary pathway
of parental care in mouth-brooding cichlid fish. \emph{Proceedings of
the Royal Society B: Biological Sciences} , \emph{265} (1411),
2217--2222.\emph{https://doi.org/10.1098/rspb.1998.0562}

Beyer, R.E., \& Fattore, J.E. (1984). The influence of age and endurance
exercise on the myoglobin concentration of skeletal muscle of the
rat.\emph{Journal of Gerontology} , \emph{39} (5),
525-530.\emph{https://doi.org/10.1093/geronj/39.5.525}

Buenrostro, J. D., Giresi, P. G., Zaba, L. C., Chang, H. Y., \&
Greenleaf, W. J. (2013). Transposition of native chromatin for fast and
sensitive epigenomic profiling of open chromatin, DNA-binding proteins
and nucleosome position. \emph{Nature Methods} , \emph{10} (12),
1213--1218.\emph{https://doi.org/10.1038/nmeth.2688}

Burns, J. G., Di Nardo, P., \& Rodd, F. H. (2009). The role of predation
in variation in body shape in guppies Poecilia reticulata: A comparison
of field and common garden phenotypes. \emph{Journal of Fish Biology}
,\emph{75} (6),
1144--1157.\emph{https://doi.org/10.1111/j.1095-8649.2009.02314.x}

Callaghan, T. V., Björn, L. O., Chernov, Y., Chapin, T., Christensen, T.
R., Huntley, B., \ldots{} Zöckler, C. (2004). Biodiversity,
distributions and adaptations of arctic species in the context of
environmental change. \emph{Ambio} , \emph{33} (7),
404--417.\emph{https://doi.org/10.1579/0044-7447-33.7.404}

Campbell, C. S., Adams, C. E., Bean, C. W., Pilakouta, N., \& Parsons,
K. J. (2021). Evolvability under climate change: Bone development and
shape plasticity are heritable and correspond with performance in Arctic
charr (Salvelinus alpinus). \emph{Evolution and Development} ,\emph{23}
(4), 333--350. https://doi.org/10.1111/ede.12379

Capecchi, G., Baldassarri, M., Ferranti, S., Guidoni, E., Cioni, M.,
Nürnberg, P., \ldots{} Grosso, S. (2018). CKAP2L mutation confirms the
diagnosis of Filippi syndrome. \emph{Clinical Genetics} , \emph{93} (5),
1109--1110.\emph{https://doi.org/10.1111/cge.13188}

Chang, C. C., Lin, W. C., Pai, L. M., Lee, H. S., Wu, S. C., Ding, S.
T., \ldots{} Sung, L. Y. (2015). Cytoophidium assembly reflects
upregulation of IMPDH activity. \emph{Journal of Cell Science}
,\emph{128} (19), 3550--3555.\emph{https://doi.org/10.1242/jcs.175265}

Chen, J. C., Hoey, D. A., Chua, M., Bellon, R., \& Jacobs, C. R. (2016).
Mechanical signals promote osteogenic fate through a primary
cilia-mediated mechanism. \emph{FASEB Journal} , \emph{30} (4),
1504--1511. https://doi.org/10.1096/fj.15-276402

Chen, Y., Lun, A.T.L., and Smyth, G.K. (2014). Differential Expression
Analysis of Complex RNA-seq Experiments Using edgeR. In Statistical
Analysis of Next Generation Sequencing Data, S. Datta, and D. Nettleton,
eds. (Cham: Springer International Publishing), pp. 51--74.

Conith, A. J., \& Albertson, R. C. (2021). The cichlid oral and
pharyngeal jaws are evolutionarily and genetically coupled. \emph{Nature
Communications} , \emph{12} (1), 1--11.
https://doi.org/10.1038/s41467-021-25755-5

Conte, M. A., Joshi, R., Moore, E. C., Nandamuri, S. P., Gammerdinger,
W. J., Roberts, R. B., \ldots{} Kocher, T. D. (2019). Chromosome-scale
assemblies reveal the structural evolution of African cichlid
genomes.\emph{GigaScience} , \emph{8} (4),
1--20.\emph{https://doi.org/10.1093/gigascience/giz030}

Cooper WJ, Parsons K, McIntyre A, Kern B, McGee-Moore A, Albertson RC
(2010) Bentho-Pelagic Divergence of Cichlid Feeding Architecture Was
Prodigious and Consistent during Multiple Adaptive Radiations within
African Rift-Lakes. \emph{PLoS ONE} \emph{5} (3): e9551.
https://doi.org/10.1371/journal.pone.0009551

Corces, M. R., Trevino, A. E., Hamilton, E. G., Greenside, P. G.,
Sinnott-Armstrong, N. A., Vesuna, S., \ldots{} Chang, H. Y. (2017). An
improved ATAC-seq protocol reduces background and enables interrogation
of frozen tissues. \emph{Nature Methods} , \emph{14} (10),
959--962.\emph{https://doi.org/10.1038/nmeth.4396}

Ding, B., Curole, J., Husemann, M., \& Danley, P. D. (2015). Habitat
complexity predicts the community diversity of rock-dwelling cichlid
fish in Lake Malawi, East Africa. \emph{Hydrobiologia} , \emph{748} (1),
133--143.\emph{https://doi.org/10.1007/s10750-014-1932-3}

Diouf, I., Derivot, L., Koussevitzky, S., Carretero, Y., Bitton, F.,
Moreau, L., \& Causse, M. (2020). Genetic basis of phenotypic plasticity
and genotype × environment interactions in a multi-parental tomato
population. \emph{Journal of Experimental Botany} , \emph{71} (18),
5365--5376. https://doi.org/10.1093/jxb/eraa265

Drost, H.-G., and Paszkowski, J. (2017). Biomartr: genomic data
retrieval with R. Bioinformatics \emph{33} ,
1216--1217.\emph{https://doi.org/10.1093/bioinformatics/btw821}.

Du, Y., Zhang, M., Liu, X., Li, Z., Hu, M., Tian, Y., \ldots{} Zhou, Y.
(2021). CDC20 promotes bone formation via APC/C dependent ubiquitination
and degradation of p65. \emph{EMBO Reports} , \emph{22} (9),
1--21.\emph{https://doi.org/10.15252/embr.202152576}

Durinck, S., Moreau, Y., Kasprzyk, A., Davis, S., De Moor, B., Brazma,
A., and Huber, W. (2005). BioMart and Bioconductor: a powerful link
between biological databases and microarray data analysis.
Bioinformatics \emph{21} , 3439--3440.
https://doi.org/10.1093/bioinformatics/bti525.

Durinck, S., Spellman, P.T., Birney, E., and Huber, W. (2009). Mapping
identifiers for the integration of genomic datasets with the
R/Bioconductor package biomaRt. Nat Protoc \emph{4} ,
1184--1191.\emph{https://doi.org/10.1038/nprot.2009.97}.

Fu, G., Xu, Z., Chen, X., Pan, H., Wang, Y., \& Jin, B. (2020). CDCA5
functions as a tumor promoter in bladder cancer by dysregulating
mitochondria-mediated apoptosis, cell cycle regulation and PI3k/AKT/mTOR
pathway activation. \emph{Journal of Cancer} , \emph{11} (9),
2408--2420.\emph{https://doi.org/10.7150/jca.35372}

Fu, L., Patel, M. S., Bradley, A., Wagner, E. F., \& Karsenty, G.
(2005). The molecular clock mediates leptin-regulated bone
formation.\emph{Cell} , \emph{122} (5),
803--815.\emph{https://doi.org/10.1016/j.cell.2005.06.028}

Fujimura, K., \& Okada, N. (2007). Development of the embryo, larva and
early juvenile of Nile tilapia Oreochromis niloticus (Pisces:
Cichlidae). Developmental staging system. \emph{Development Growth and
Differentiation} , \emph{49} (4), 301--324.
https://doi.org/10.1111/j.1440-169X.2007.00926.x

Fujimura, K., \& Okada, N. (2008a). Bone development in the jaw of Nile
tilapia Oreochromis niloticus (Pisces: Cichlidae). \emph{Development
Growth and Differentiation} , \emph{50} (5),
339--355.\emph{https://doi.org/10.1111/j.1440-169X.2008.01032.x}

Fujimura, K., \& Okada, N. (2008b). Shaping of the lower jaw bone during
growth of Nile tilapia Oreochromis niloticus and a Lake Victoria cichlid
Haplochromis chilotes: A geometric morphometric
approach.\emph{Development Growth and Differentiation} , \emph{50} (8),
653--663. https://doi.org/10.1111/j.1440-169X.2008.01063.x

Gao, Y., Guo, C., Fu, S., Cheng, Y., Song, C. (2021). Downregulation of
CDC20 suppressed cell proliferation, induced apoptosis, triggered cell
cycle arrest in osteosarcoma cells, and enhances chemosensitivity to
cisplatin. \emph{Neoplasma} , \emph{68} (2),
382--390.\emph{https://doi.org/10.4149/neo}

Genner, M. J., Taylor, M. I., Cleary, D. F. R., Hawkins, S. J., Knight,
M. E., \& Turner, G. F. (2004). Beta diversity of rock-restricted
cichlid fishes in Lake Malawi: Importance of environmental and spatial
factors. \emph{Ecography} , \emph{27} (5), 601--610.
https://doi.org/10.1111/j.0906-7590.2004.03824.x

Genner, M. J., Turner, G. F., \& Hawkins, S. J. (1999). Foraging of
rocky habitat cackled fishes in Lake Malawi: Coexistence through niche
partitioning? \emph{Oecologia} , \emph{121} (2),
283--292.\emph{https://doi.org/10.1007/s004420050930}

Gibert, J. M. (2017). The flexible stem hypothesis: evidence from
genetic data. \emph{Development Genes and Evolution} , \emph{227} (5),
297--307. https://doi.org/10.1007/s00427-017-0589-0

Gilbert, M. C., Tetrault, E., Packard, M., Navon, D., \& Albertson, R.
C. (2021). Ciliary Rootlet Coiled-Coil 2 (crocc2) Is Associated with
Evolutionary Divergence and Plasticity of Cichlid Jaw
Shape.\emph{Molecular Biology and Evolution} , \emph{38} (8),
3078--3092. https://doi.org/10.1093/molbev/msab071

Gordon, B. S., Rossetti, M. L., \& Eroshkin, A. M. (2019). Arrdc2 and
Arrdc3 elicit divergent changes in gene expression in skeletal muscle
following anabolic and catabolic stimuli. \emph{Physiological Genomics}
,\emph{51} (6),
208--217.\emph{https://doi.org/10.1152/physiolgenomics.00007.2019}

Govey, P. M., Kawasawa, Y. I., \& Donahue, H. J. (2015). Mapping the
osteocytic cell response to fluid flow using RNA-Seq. \emph{Journal of
Biomechanics} , \emph{48} (16),
4327--4332.\emph{https://doi.org/10.1016/j.jbiomech.2015.10.045}

Gugger, S., Kesselring, H., Stöcklin, J., \& Hamann, E. (2015). Lower
plasticity exhibited by high- versus mid-elevation species in their
phenological responses to manipulated temperature and
drought.\emph{Annals of Botany} , \emph{116} (6),
953--962.\emph{https://doi.org/10.1093/aob/mcv155}

Gunst, S. J., \& Zhang, W. (2008). Actin cytoskeletal dynamics in smooth
muscle: A new paradigm for the regulation of smooth muscle
contraction.\emph{American Journal of Physiology - Cell Physiology} ,
\emph{295} (3),
576--587.\emph{https://doi.org/10.1152/ajpcell.00253.2008}

Gunter, H. M., Fan, S., Xiong, F., Franchini, P., Fruciano, C., \&
Meyer, A. (2013). Shaping development through mechanical strain: The
transcriptional basis of diet-induced phenotypic plasticity in a cichlid
fish. \emph{Molecular Ecology} , \emph{22} (17),
4516--4531.\emph{https://doi.org/10.1111/mec.12417}

Hohenlohe, P. A. (2014). Ecological genomics in full
colour.\emph{Molecular Ecology} , \emph{23} (21),
5129--5131.\emph{https://doi.org/10.1111/mec.12945}

Hu, Y., \& Albertson, R. C. (2014). Hedgehog signaling mediates adaptive
variation in a dynamic functional system in the cichlid feeding
apparatus. \emph{Proceedings of the National Academy of Sciences of the
United States of America} , \emph{111} (23),
8530--8534.\emph{https://doi.org/10.1073/pnas.1323154111}

Hu, Y., \& Albertson, R. C. (2017). Baby fish working out: An epigenetic
source of adaptive variation in the cichlid jaw. \emph{Proceedings of
the Royal Society B: Biological Sciences} , \emph{284}
(1860).\emph{https://doi.org/10.1098/rspb.2017.1018}

Hulsey, C. D., Meyer, A., \& Streelman, J. T. (2020). Convergent
Evolution of Cichlid Fish Pharyngeal Jaw Dentitions in Mollusk-Crushing
Predators: Comparative X-Ray Computed Tomography of Tooth Sizes,
Numbers, and Replacement. \emph{Integrative and Comparative Biology}
,\emph{60} (3), 656--664. https://doi.org/10.1093/icb/icaa089

Huysseune, A. (1995). Phenotypic plasticity in the lower pharyngeal jaw
dentition of astatoreochromis alluaudi (teleostei:
cichlidae).\emph{Archives of Oral Biology} , \emph{40} (11),
1005--1014.\emph{https://doi.org/10.1016/0003-9969(95)00074-Y}

Karasz, D. C., Weaver, A. I., Buckley, D. H., \& Wilhelm, R. C. (2022).
Conditional filamentation as an adaptive trait of bacteria and its
ecological significance in soils. \emph{Environmental Microbiology}
,\emph{24} (1), 1--17.\emph{https://doi.org/10.1111/1462-2920.15871}

Kawajiri, M., Yoshida, K., Fujimoto, S., Mokodongan, D. F., Ravinet, M.,
Kirkpatrick, M., \ldots{} Kitano, J. (2014). Ontogenetic stage-specific
quantitative trait loci contribute to divergence in developmental
trajectories of sexually dimorphic fins between medaka
populations.\emph{Molecular Ecology} , \emph{23} (21), 5258--5275.
https://doi.org/10.1111/mec.12933

Kelly, N. H., Schimenti, J. C., Ross, F. P., \& van der Meulen, M. C. H.
(2016). Transcriptional profiling of cortical versus cancellous bone
from mechanically-loaded murine tibiae reveals differential gene
expression. \emph{Bone} , \emph{86} ,
22--29.\emph{https://doi.org/10.1016/j.bone.2016.02.007}

Kocher, T. D. (2004). Adaptive evolution and explosive speciation: The
cichlid fish model. \emph{Nature Reviews Genetics} , \emph{5} (4),
288--298. https://doi.org/10.1038/nrg1316

Küttner, E., Parsons, K. J., Easton, A. A., Skúlason, S., Danzmann, R.
G., \& Ferguson, M. M. (2014). Hidden genetic variation evolves with
ecological specialization: The genetic basis of phenotypic plasticity in
Arctic charr ecomorphs. \emph{Evolution and Development} , \emph{16}
(4), 247--257.\emph{https://doi.org/10.1111/ede.12087}

Lafuente, E., Duneau, D., \& Beldade, P. (2018). Genetic basis of
thermal plasticity variation in Drosophila melanogaster body
size.\emph{PLoS Genetics} , \emph{14} (9), 1--24.
https://doi.org/10.1371/journal.pgen.1007686

Laitinen, R. A. E., \& Nikoloski, Z. (2019). Genetic basis of plasticity
in plants. \emph{Journal of Experimental Botany} , \emph{70} (3),
795--804.\emph{https://doi.org/10.1093/jxb/ery404}

Landy, J. A., Oschmann, A., Munch, S. B., \& Walsh, M. R. (2020).
Ancestral genetic variation in phenotypic plasticity underlies rapid
evolutionary changes in resurrected populations of
waterfleas.\emph{Proceedings of the National Academy of Sciences of the
United States of America} , \emph{117} (51),
32535--32544.\emph{https://doi.org/10.1073/pnas.2006581117}

Langmead, B., and Salzberg, S.L. (2012). Fast gapped-read alignment with
Bowtie 2. Nat Methods \emph{9} ,
357--359.\emph{https://doi.org/10.1038/nmeth.1923}.

Ledón-Rettig, C. C., Pfennig, D. W., \& Crespi, E. J. (2010). Diet and
hormonal manipulation reveal cryptic genetic variation: Implications for
the evolution of novel feeding strategies. \emph{Proceedings of the
Royal Society B: Biological Sciences} , \emph{277} (1700),
3569--3578.\emph{https://doi.org/10.1098/rspb.2010.0877}

Levis, N. A., Reed, E. M. X., Pfennig, D. W., \& Burford Reiskind, M. O.
(2020). Identification of candidate loci for adaptive phenotypic
plasticity in natural populations of spadefoot toads. \emph{Ecology and
Evolution} , \emph{10} (16), 8976--8988.
https://doi.org/10.1002/ece3.6602

Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N.,
Marth, G., Abecasis, G., and Durbin, R. (2009). The Sequence
Alignment/Map format and SAMtools. Bioinformatics \emph{25} ,
2078--2079.\emph{https://doi.org/10.1093/bioinformatics/btp352}.

Lloyd, E., Chhouk, B., Conith, A. J., Keene, A. C., \& Albertson, R. C.
(2021). Diversity in rest-activity patterns among Lake Malawi cichlid
fishes suggests a novel axis of habitat partitioning. \emph{Journal of
Experimental Biology} , \emph{224}
(7).\emph{https://doi.org/10.1242/jeb.242186}

López-Fernández, H., Arbour, J., Willis, S., Watkins, C., Honeycutt,
R.L., Winemiller, K.O. (2014) Morphology and Efficiency of a Specialized
Foraging Behavior, Sediment Sifting, in Neotropical Cichlid
Fishes.\emph{PLoS ONE} , \emph{9} (3): e89832.
https://doi.org/10.1371/journal.pone.0089832

Maan, M. E., Hofker, K. D., van Alphen, J. J. M., \& Seehausen, O.
(2006). Sensory drive in cichlid speciation. \emph{The American
Naturalist} , \emph{167} (6),
947--954.\emph{https://doi.org/10.1086/503532}

Machado-Schiaffino, G., Henning, F., \& Meyer, A. (2014).
Species-specific differences in adaptive phenotypic plasticity in an
ecologically relevant trophic trait: Hypertrophic lips in midas cichlid
fishes. \emph{Evolution} , \emph{68} (7),
2086--2091.\emph{https://doi.org/10.1111/evo.12367}

Malinsky, M., Svardal, H., Tyers, A. M., Miska, E. A., Genner, M. J.,
Turner, G. F., \& Durbin, R. (2018). Whole-genome sequences of Malawi
cichlids reveal multiple radiations interconnected by gene
flow.\emph{Nature Ecology and Evolution} , \emph{2} (12),
1940--1955.\emph{https://doi.org/10.1038/s41559-018-0717-x}

Mantila Roosa, S. M., Liu, Y., \& Turner, C. H. (2011). Gene expression
patterns in bone following mechanical loading. \emph{Journal of Bone and
Mineral Research} , \emph{26} (1), 100--112.
https://doi.org/10.1002/jbmr.193

Mao, J., McGlinn, E., Huang, P., Tabin, C. J., \& McMahon, A. P. (2009).
Fgf-Dependent Etv4/5 Activity Is Required for Posterior Restriction of
Sonic hedgehog and Promoting Outgrowth of the Vertebrate
Limb.\emph{Developmental Cell} , \emph{16} (4),
600--606.\emph{https://doi.org/10.1016/j.devcel.2009.02.005}

Margres, M.J., Wray, K.P., Seavy, M., McGivern, J.J., Sanader, D.,
Rokyta, D.R. (2015). Phenotypic integration in the feeding system of the
eastern diamondback rattlesnake (\emph{Crotalus adamanteus}
).\emph{Molecular Ecology} , \emph{24} (13), 3405-3420.

Martin, M. (2011). Cutadapt removes adapter sequences from
high-throughput sequencing reads.

Mccairns, R. J. S., \& Bernatchez, L. (2012). Plasticity and
heritability of morphological variation within and between parapatric
stickleback demes. \emph{Journal of Evolutionary Biology} , \emph{25}
(6), 1097--1112.\emph{https://doi.org/10.1111/j.1420-9101.2012.02496.x}

Mcguigan, K., Nishimura, N., Currey, M., Hurwit, D., \& Cresko, W. A.
(2011). Cryptic genetic variation and body size evolution in threespine
stickleback. \emph{Evolution} , \emph{65} (4),
1203--1211.\emph{https://doi.org/10.1111/j.1558-5646.2010.01195.x}

McKaye, K.R.; Kocher, T.; Harrison, R.; Kornfield, I. (1984). Genetic
evidence for allopatric and sympatric differentiation among color morphs
of a Lake Malawi cichlid fish. \emph{Evolution} , \emph{38} (1),
215--219.

Meuthen, D., Baldauf, S. A., Bakker, T. C. M., \& Thünken, T. (2018).
Neglected patterns of variation in phenotypic plasticity: Age- and
sex-specific antipredator plasticity in a cichlid fish. \emph{American
Naturalist} , \emph{191} (4),
475--490.\emph{https://doi.org/10.1086/696264}

Meyer, A. (1987). Phenotypic Plasticity and Heterochrony in Cichlasoma
managuense (Pisces, Chichlidae) and their Implications for Speciation in
Cichlid Fishes. \emph{Evolution} , \emph{41} (6), 1357.
https://doi.org/10.2307/2409100

Mirvis, M., Stearns, T., \& Nelson, W. J. (2018). Cilium structure,
assembly, and disassembly regulated by the
cytoskeleton.\emph{Biochemical Journal} , \emph{475} (14),
2329--2353.\emph{https://doi.org/10.1042/BCJ20170453}

Moore, E. R., Chen, J. C., \& Jacobs, C. R. (2019). Prx1-Expressing
Progenitor Primary Cilia Mediate Bone Formation in response to
Mechanical Loading in Mice. \emph{Stem Cells International} ,\emph{2019}
.\emph{https://doi.org/10.1155/2019/3094154}

Morgan, R., Andreassen, A. H., Åsheim, E. R., Finnøen, M. H., \&
Dresler, G. (2022). \emph{Reduced physiological plasticity in a fish
adapted to stable conditions} . 1--37.
https://doi.org/10.1073/pnas.2201919119/-/DCSupplemental.Published

Muschick, M., Barluenga, M., Salzburger, W., \& Meyer, A. (2011).
Adaptive phenotypic plasticity in the Midas cichlid fish pharyngeal jaw
and its relevance in adaptive radiation. \emph{BMC Evolutionary Biology}
, \emph{11} (116).

Navon, D., Hatini, P., Zogbaum, L., \& Albertson, R. C. (2021). The
genetic basis of coordinated plasticity across functional units in a
Lake Malawi cichlid mapping population. \emph{Evolution} , \emph{75}
(3), 672--687.\emph{https://doi.org/10.1111/evo.14157}

Navon, D., Male, I., Tetrault, E. R., Aaronson, B., Karlstrom, R. O., \&
Craig Albertson, R. (2020). Hedgehog signaling is necessary and
sufficient to mediate craniofacial plasticity in
teleosts.\emph{Proceedings of the National Academy of Sciences of the
United States of America} , \emph{117} (32),
19321--19327.\emph{https://doi.org/10.1073/pnas.1921856117}

Neph, S., Kuehn, M.S., Reynolds, A.P., Haugen, E., Thurman, R.E.,
Johnson, A.K., Rynes, E., Maurano, M.T., Vierstra, J., Thomas, S., et
al. (2012). BEDOPS: high-performance genomic feature operations.
Bioinformatics \emph{28} ,
1919--1920.\emph{https://doi.org/10.1093/bioinformatics/bts277}.

O'Connor, M. J., Thakar, T., Nicolae, C. M., \& Moldovan, G. L. (2021).
PARP14 regulates cyclin D1 expression to promote cell-cycle
progression.\emph{Oncogene} , \emph{40} (30),
4872--4883.\emph{https://doi.org/10.1038/s41388-021-01881-8}

Pagano, M., Pepperkok, R., Verde, F., Ansorge, W., \& Draetta, G.
(1992). Cyclin A is required at two points in the human cell
cycle.\emph{EMBO Journal} , \emph{11} (3),
961--971.\emph{https://doi.org/10.1002/j.1460-2075.1992.tb05135.x}

Parsons, K. J., Concannon, M., Navon, D., Wang, J., Ea, I., Groveas, K.,
\ldots{} Albertson, R. C. (2016). Foraging environment determines the
genetic architecture and evolutionary potential of trophic morphology in
cichlid fishes. \emph{Molecular Ecology} , \emph{25} (24), 6012--6023.
https://doi.org/10.1111/mec.13801

Parsons, K. J., Trent Taylor, A., Powder, K. E., \& Albertson, R. C.
(2014). Wnt signalling underlies the evolution of new phenotypes and
craniofacial variability in Lake Malawi cichlids. \emph{Nature
Communications} , \emph{5} ,
1--11.\emph{https://doi.org/10.1038/ncomms4629}

Pigliucci, M., Murren, C. J., \& Schlichting, C. D. (2006). Phenotypic
plasticity and evolution by genetic assimilation. \emph{Journal of
Experimental Biology} , \emph{209} (12),
2362--2367.\emph{https://doi.org/10.1242/jeb.02070}

Powder, K. E., \& Albertson, R. C. (2016). Cichlid fishes as a model to
understand normal and clinical craniofacial
variation.\emph{Developmental Biology} , \emph{415} (2),
338--346.\emph{https://doi.org/10.1016/j.ydbio.2015.12.018}

Powder, K. E., Cousin, H., McLinden, G. P., \& Craig Albertson, R.
(2014). A nonsynonymous mutation in the transcriptional regulator lbh is
associated with cichlid craniofacial adaptation and neural crest cell
development. \emph{Molecular Biology and Evolution} , \emph{31} (12),
3113--3124. https://doi.org/10.1093/molbev/msu267

Quinlan, A.R., and Hall, I.M. (2010). BEDTools: a flexible suite of
utilities for comparing genomic features. Bioinformatics \emph{26} ,
841--842. https://doi.org/10.1093/bioinformatics/btq033.

R Core Team (2021). R: A language and environment for statistical
computing.

Raab-Cullen, D.M., Thiede, M.A., Peterson, D.N., Kimmel, D.B., \&
Recker, R.R. (1994). Mechanical loading stimulates rapid changes in
periosteal gene expression. \emph{Calcified Tissue International}
,\emph{55} , 473-478.

Rankin, S., Ayad, N. G., \& Kirschner, M. W. (2005). Sororin, a
substrate of the anaphase- promoting complex, is required for sister
chromatid cohesion in vertebrates. \emph{Molecular Cell} , \emph{18}
(2), 185--200.\emph{https://doi.org/10.1016/j.molcel.2005.03.017}

Ribbink, B.J., Marsh, A.J., Marsh, B.A., Ribbink, A.C., \& Sharp, A.C.
(1983). A preliminary survey of the cichlid fishes of rocky habitats in
Lake Malawi: results-The Mbuna-Pseudotropheus, \emph{African Zoology}
,\emph{13} (3).

Roberts, R. B., Hu, Y., Albertson, R. C., \& Kocher, T. D. (2011).
Craniofacial divergence and ongoing adaptation via the hedgehog
pathway.\emph{Proceedings of the National Academy of Sciences of the
United States of America} , \emph{108} (32),
13194--13199.\emph{https://doi.org/10.1073/pnas.1018456108}

Robinson, M.D., McCarthy, D.J., and Smyth, G.K. (2010). edgeR: a
Bioconductor package for differential expression analysis of digital
gene expression data. Bioinformatics \emph{26} ,
139--140.\emph{https://doi.org/10.1093/bioinformatics/btp616}.

Salzburger, W. (2009). The interaction of sexually and naturally
selected traits in the adaptive radiations of cichlid
fishes.\emph{Molecular Ecology} , \emph{18} (2),
169--185.\emph{https://doi.org/10.1111/j.1365-294X.2008.03981.x}

Schartau, J. M., Sjögreen, B., Gagnon, Y. L., \& Kröger, R. H. H.
(2009). Optical Plasticity in the Crystalline Lenses of the Cichlid Fish
Aequidens pulcher. \emph{Current Biology} , \emph{19} (2),
122--126.\emph{https://doi.org/10.1016/j.cub.2008.11.062}

Schluter, D. (2000). \emph{The ecology of adaptive radiation} . Oxford
University Press.

Schmitz, J., Watrin, E., Lénárt, P., Mechtler, K., \& Peters, J. M.
(2007). Sororin Is Required for Stable Binding of Cohesin to Chromatin
and for Sister Chromatid Cohesion in Interphase. \emph{Current Biology}
,\emph{17} (7),
630--636.\emph{https://doi.org/10.1016/j.cub.2007.02.029}

Schneider, R. F., Li, Y., Meyer, A., \& Gunter, H. M. (2014). Regulatory
gene networks that shape the development of adaptive phenotypic
plasticity in a cichlid fish. \emph{Molecular Ecology} , \emph{23} (18),
4511--4526.\emph{https://doi.org/10.1111/mec.12851}

Schock, E. N., Chang, C. F., Youngworth, I. A., Davey, M. G., Delany, M.
E., \& Brugmann, S. A. (2016) Utilizing the chicken as an animal model
for human craniofacial ciliopathies. \emph{Developmental Biology} .
415(2):326-337. doi: 10.1016/j.ydbio.2015.10.024

Seehausen, O. (2006). African cichlid fish: A model system in adaptive
radiation research. \emph{Proceedings of the Royal Society B: Biological
Sciences} , \emph{273} (1597),
1987--1998.\emph{https://doi.org/10.1098/rspb.2006.3539}

Seehausen, O., Mayhew, P. J., \& Van Alphen, J. J. M. (1999). Evolution
of colour patterns in East African cichlid fish. \emph{Journal of
Evolutionary Biology} , \emph{12} (3), 514--534.
https://doi.org/10.1046/j.1420-9101.1999.00055.x

Siepielski, A. M., Morrissey, M. B., Buoro, M., Carlson, S. M., Caruso,
C. M., Clegg, S. M., \ldots{} Maccoll, A. D. C. (2017). \emph{Science}
.\emph{355} , 959--962.

Sih, A., Ferrari, M. C. O., \& Harris, D. J. (2011). Evolution and
behavioural responses to human-induced rapid environmental
change.\emph{Evolutionary Applications} , \emph{4} (2),
367--387.\emph{https://doi.org/10.1111/j.1752-4571.2010.00166.x}

Shekhar, R., Priyanka, P., Kumar, P., Ghosh, T., Khan, M., Nagarajan,
P., \& Saxena, S. (2019). The microRNAs miR-449a and miR-424 suppress
osteosarcoma by targeting cyclin A2 expression. \emph{Journal of
Biological Chemistry} , \emph{294} (12),
4381--4400.\emph{https://doi.org/10.1074/jbc.RA118.005778}

St-Jacques, B., Hammerschmidt, M., \& Mcmahon, A. P. (1999).\emph{Indian
hedgehog signaling regulates proliferation and differentiation of
chondrocytes and is essential for bone formation} . Retrieved from
www.genesdev.org

Stark, R., and Brown, G. DiffBind: Differential binding analysis of
ChIP-Seq peak data - April 2021. 71. .

Sturmbauer, C., \& Meyer, A. (1992). Genetic divergence, speciation and
morphological stasis in a lineage of African cichlid
fishes.\emph{Nature} , \emph{358} (6387),
578--581.\emph{https://doi.org/10.1038/358578a0}

Styczynska-Soczka, K., \& Jarman, A. P. (2015). The Drosophila homologue
of Rootletin is required for mechanosensory function and ciliary rootlet
formation in chordotonal sensory neurons. \emph{Cilia} , (SUPPLEMENT 1),
1--11.\emph{https://doi.org/10.1186/s13630-015-0018-9}

Sun, J., \& Deng, W. M. (2007). Hindsight Mediates the Role of Notch in
Suppressing Hedgehog Signaling and Cell
Proliferation.\emph{Developmental Cell} , \emph{12} (3),
431--442.\emph{https://doi.org/10.1016/j.devcel.2007.02.003}

Tange, O. (2018). Gnu Parallel 2018 (Zenodo).

Tiet, T. D., Hopyan, S., Nadesan, P., Gokgoz, N., Poon, R., Lin, A. C.,
\ldots{} Wunder, J. S. (2006). Constitutive Hedgehog signaling in
chondrosarcoma up-regulates tumor cell proliferation. \emph{American
Journal of Pathology} , \emph{168} (1),
321--330.\emph{https://doi.org/10.2353/ajpath.2006.050001}

van Heerwaarden, B., \& Sgrò, C. M. (2017). The quantitative genetic
basis of clinal divergence in phenotypic plasticity. \emph{Evolution}
,\emph{71} (11), 2618--2633.\emph{https://doi.org/10.1111/evo.13342}

van Snick Gray, E., \& Stauffer, J. (2004). Phenotypic plasticity: its
role in trophic radiation and explosive speciation in cichlids
(Teleostei: Cichlidae), \emph{Animal Biology} , \emph{54} (2), 137-158.
doi:\emph{https://doi.org/10.1163/1570756041445191}

Visintin, R., Prinz, S., \& Amon, A. (1997). CDC20 and CDH1: A family of
substrate-specific activators of APC- dependent
proteolysis.\emph{Science} , \emph{278} (5337),
460--463.\emph{https://doi.org/10.1126/science.278.5337.460}

Westbrook, L., Manuvakhova, M., Kern, F. G., Estes, N. R., Ramanathan,
H. N., \& Thottassery, J. V. (2007). Cks1 regulates cdk1 expression: A
novel role during mitotic entry in breast cancer cells. \emph{Cancer
Research} , \emph{67} (23),
11393--11401.\emph{https://doi.org/10.1158/0008-5472.CAN-06-4173}

Westneat, M. W. (1991). Linkage Biomechanics and Evolution of the Unique
Feeding Mechanism of Epibulus Insidiator (Labridae:
Teleostei).\emph{Journal of Experimental Biology} , \emph{159} (1),
165--184.\emph{https://doi.org/10.1242/jeb.159.1.165}

Willis, C. G., Ruhfel, B., Primack, R. B., Miller-Rushing, A. J., \&
Davis, C. C. (2008). Phylogenetic patterns of species loss in Thoreau's
woods are driven by climate change. \emph{Proceedings of the National
Academy of Sciences of the United States of America} , \emph{105} (44),
17029--17033.\emph{https://doi.org/10.1073/pnas.0806446105}

Wimberger, P. H. (1991). Plasticity of jaw and skull morphology in the
neotropical cichlids Geophagus Brasiliensis and G.
Steindachneri.\emph{Evolution. International Journal of Organic
Evolution} ,\emph{45} (7), 1545--1563.

Yang, J., Liu, X., Yue, G., Adamian, M., Bulgakov, O., \& Li, T. (2002).
Rootletin, a novel coiled-coil protein, is a structural component of the
ciliary rootlet. \emph{Journal of Cell Biology} , \emph{159} (3),
431--440. https://doi.org/10.1083/jcb.200207153

Yates, A.D., Achuthan, P., Akanni, W., Allen, J., Allen, J.,
Alvarez-Jarreta, J., Amode, M.R., Armean, I.M., Azov, A.G., Bennett, R.,
et al. (2020). Ensembl 2020. Nucleic Acids Research \emph{48} ,
D682--D688.\emph{https://doi.org/10.1093/nar/gkz966}.

Yoshida, T., Shimada, K., Oma, Y., Kalck, V., Akimura, K., Taddei, A.,
\ldots{} Harata, M. (2010). Actin-related protein Arp6 influences
H2A.Z-dependent and -independent gene expression and links ribosomal
protein genes to nuclear pores. \emph{PLoS Genetics} , \emph{6} (4),
10--17.\emph{https://doi.org/10.1371/journal.pgen.1000910}

Yoshioka, M., Boivin, A., Ye, P., Labrie, F., \& St-Amand, J. (2006).
Effects of dihydrotestosterone on skeletal muscle transcriptome in mice
measured by serial analysis of gene expression. \emph{Journal of
Molecular Endocrinology} , \emph{36} (2), 247--259.
https://doi.org/10.1677/jme.1.01964

Young KA, Snoeks J, Seehausen O (2009) Morphological Diversity and the
Roles of Contingency, Chance and Determinism in African Cichlid
Radiations. PLoS ONE 4(3):
e4740.\emph{https://doi.org/10.1371/journal.pone.0004740}

Yu, H. (2002). Regulation of APC -- Cdc20 by the spindle checkpoint
Hongtao Yu. \emph{Current Opinion in Cell Biology} , 706--714.

Zaman, F., Zhao, Y., Celvin, B., Mehta, H. H., Wan, J., Chrysis, D.,
\ldots{} Sävendahl, L. (2019). Humanin is a novel regulator of Hedgehog
signaling and prevents glucocorticoid-induced bone growth
impairment.\emph{FASEB Journal} , \emph{33} (4),
4962--4974.\emph{https://doi.org/10.1096/fj.201801741R}

Zhao, L. J., Liu, X. G., Liu, Y. Z., Liu, Y. J., Papasian, C. J., Sha,
B. Y., \ldots{} Deng, H. W. (2010). Genome-wide association study for
femoral neck bone geometry. \emph{Journal of Bone and Mineral Research}
,\emph{25} (2), 320--329.\emph{https://doi.org/10.1359/jbmr.090726}

Zogbaum, L., Friend, P. G., \& Albertson, R. C. (2021). Plasticity and
genetic basis of cichlid gill arch anatomy reveal novel roles for
Hedgehog signaling. \emph{Molecular Ecology} , \emph{30} (3),
761--774.\emph{https://doi.org/10.1111/mec.15766}

\textbf{TABLE 1:} Number of genes that are DE or DA from RNA-seq and
ATAC-seq datasets, respectively. In each comparison/model, the total
numbers are provided, followed by the numbers that are upregulated in a
particular species or environment.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/image1}

\end{center}
\end{figure}

\textbf{FIGURES}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image2/image2}

\end{center}
\end{figure}

\textbf{FIGURE 1: African cichlids from Lake Malawi display differences
in foraging anatomy.} (A) \emph{Tropheops sp} . ``red cheek'' (TRC) has
a more downturned oral jaw apparatus, and a steeply sloping craniofacial
profile, adapted for benthic foraging via scraping/biting/nipping. (B) A
schematic of the opercle 4-bar linkage chain, which is critical for
lower jaw depression, is shown for TRC. (C) \emph{Maylandia zebra} (MZ)
is characterized by a more upturned oral jaw, better suited for pelagic
feeding via fast jaw rotation. (D) A schematic of the opercle 4-bar
linkage chain is shown for MZ. Relative to TRC, MZ possesses a longer
coupler link (CL) and shorter output link (OL). (A, C) Images courtesy
of Ad Koning at Cichlid Press. (B, D) Red bars indicate movable linkages
- input link (IL), CL, and OL - while the black bar represents the fixed
link. Blue depicts the interopercle (IOP) bone, and green shows the
retroarticular process (RA) of the lower jaw, while orange is the
interopercle-mandibular (ioml) ligament that connects the two bones.
Throughout the text we refer to this as the IOP-RA functional complex.
(E) Experimental schematic in which we fed cichlids either flake or rock
food for 17 days (RNA-seq) or 28 days (ATAC-seq), with a 1 week training
period. Inset shows an external view of the IOP-RA complex dissected for
all experiments.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image3/image3}

\end{center}
\end{figure}

\textbf{FIGURE 2: Differential expression between species is more robust
than between foraging environments.} (A) Volcano plot of the pairwise
comparison between MZ and TRC across environments, showing a large
number of DEGs (n=5318) with roughly equal numbers of significantly
upregulated genes between species (n=2667 MZ; n=2651 TRC). Given the
nature of the comparison, genes considered upregulated in MZ are
downregulated in TRC, and vice versa. Red indicates upregulated genes
for MZ, blue depicts genes upregulated in TRC, while black represents
genes that do not meet the significance threshold of \textless{}0.05
FDR. (B) Heatmap of the top 500 most variable genes in the RNA-seq
dataset is shown. Individuals from the pelagic foraging treatment are
labeled with gray triangles, while benthic individuals are labeled with
black triangles. Species cluster together, but there is less obvious
structuring by foraging environment, although MZ segregate by
environment more so than TRC. Photographs courtesy of Ad Koning at
Cichlid Press. (C) Enriched GO terms associated with genes upregulated
in the MZ:TRC comparison across environments. Colors are representative
of the -log10(p-value). Gray bars indicate no significance. Gene counts
are given as a total for both species along with the corresponding GO
term.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image4/image4}

\end{center}
\end{figure}

\textbf{FIGURE 3: Little differential expression is detected between
foraging environments.} (A) Within TRC, the volcano plot shows no DEGs
between environments. (B) Within MZ, a small number of DEGs (n=38) were
detected between environments, with more upregulated in the pelagic
(n=25; blue) versus the benthic (n=13; red) environment. Black dots
represent genes that are not significantly DE at FDR \textless{}0.05.
(C) Within MZ, more GO terms were returned for animals exposed to
pelagic versus benthic environments. Cell-cycle regulation features
prominently in the pelagic environment.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image5/image5}

\end{center}
\end{figure}

\textbf{FIGURE 4: The pelagic environment drives differences in
expression between species and implicates cell cycle regulation as a
mode of increased bone deposition.} (A) Volcano plot depicting DE
between species in the benthic environment (n=984 total; n=549 MZ in
red; n=435 TRC in blue). (B) Nearly four times the number of DEGs are
detected in the pelagic environment (n=3761 total; n=1927 MZ in red;
n=1834 TRC in blue). (C) The venn diagram shows that not only are there
more DEGs in the pelagic versus benthic environment, but that most
(813/984 = 83\%) of the DEGs detected in the benthic environment are
also DE in the pelagic environment. Alternatively, only 22\% (813/3761)
DEGs detected in the pelagic environment are also DE in benthic fish.
(D, E) Enriched GO terms are shown for MZ versus TRC in the pelagic (D)
and benthic (E) environments. Colors are representative of the
-log10(p-value). Gray bars indicate no significance. Gene counts are
given as a total for both species along with the corresponding GO term.
Many terms enriched for in MZ in both environments are associated with
cell cycle (e.g., cell cycle, cell division, etc).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image6/image6}

\end{center}
\end{figure}

\textbf{FIGURE 5: Differential expression across species plus
environment reveals signatures of genetic assimilation and further
supports a role for the cell cycle.} (A) We constructed an additive
model to identify genes that were differentially expressed between
species and environment (S+E), which identified 128 DEGs. With the
exception of a single MZ individual, species clustered together.
Clustering by environment was observed in MZ but not TRC. Assuming that
plasticity represents the ancestral condition, this analysis also
provides evidence for genetic assimilation, whereby pelagic MZ exhibited
relatively high gene expression compared to benthic MZ, which generally
resemble TRC in terms of cpm counts. Pink dots on the cladogram to the
left of the heatmap denote gene clusters that exemplify this pattern. We
note that 3 MZ exhibited especially robust expression levels. (B)
Enriched GO terms associated with genes in the S+E model. Colors are
representative of the -log10(p-value). Gray bars indicate no
significance. The nature of the additive model precludes us from having
an up- vs down-regulated analysis of genes, because it takes into
account both species and environment at the same time. These GO terms
are similar to previous analyses in returning processes involved in cell
cycle regulation.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image7/image7}

\end{center}
\end{figure}

\textbf{FIGURE 6: The overlap of RNA-seq and ATAC-seq datasets narrows
the list of candidate genes.} (A) In total, 89 genes were identified
that were both differentially expressed (DE) and differentially
accessible (DA) between species. Of these, 15 overlapped with genes
identified in the pelagic dataset (light blue), 11 overlapped with genes
from the benthic dataset (dark blue), and 13 were identified in both
foraging environments. These overlapping datasets identified a
relatively small subset of genes where expression differences may be due
to differences in the \emph{cis-} regulatory region. Boxes to the right
of the venn diagram display the specific genes from each area of
overlap. (B-D) Diagrams of genes are shown (each panel has its own scale
bar), as well as the location of ATAC-seq peaks based on DA analyses,
and CpG islands (from the UCSC genome browser CpG island track based on
at least 50\% GC content, \textgreater{}200 bp, and \textgreater{}0.6
ratio of observed number of CG dinucleotides). Shown are representative
data from the benthic dataset (B, \emph{actr6} on LG7), the pelagic
dataset (C, \emph{kiaa0586} on an unlinked contig), and a gene that was
significantly DE and DA in both environments (D, \emph{capn1-like} on
LG10). Beneath each gene model, ATAC-seq peaks are colored based on
concentrations, with lighter colors indicating lower concentrations and
darker colors indicating higher concentrations.

\selectlanguage{english}
\FloatBarrier
\end{document}

