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,ccna2 is a cyclin that activates cdk2 and promotes cell cycle progression through both G1/S and G2/M phases (Pagano et al., 1992), and cks1b slows the progression of G1/S and can block entry to M phase (Westbrook et al 2007). In addition, cdc20regulates metaphase-anaphase transition during mitosis via activation of APC, which targets proteins for degradation (Visintin et al., 1997; Yu 2002). Finally, 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.,arrdc3b , Gordon et al., 2019; myoglobin, Beyer et al., 1984), and osteoblast proliferation via regulation of cyclins(e.g., 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 >0.05, might prove fruitful. For example, a transcript that was upregulated in benthic TRC (logFC = 5.44) corresponds to receptor transporting protein 3(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), poly(ADP-ribose) polymerase 14 (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.