Assessing the Levels of Functional Adaptation: Finite Element Analysis Reveals Species, Hybrid, and Sexual Variation in the Biomechanics of African Cichlid Mandibles

To understand how adaptive divergence emerges it is essential to examine the function of phenotypic traits along a continuum. For vertebrates, the mandible provides a key link with foraging and other important activities which has made it highly relevant for investigations of biomechanical change. Variation in mandible shape is known to correspond with ecology but its function is often only investigated between distinct species. However, for such divergence to occur and be maintained selection likely draws from many sources of biomechanical variation. African cichlids represent an exemplar model for understanding how such processes unfold with mandible variation existing between species, sexes, and is likely generated in nature by the potential for hybridization. We explored such mandible variation through a finite element modelling approach and predicted that hybrids and females would have reduced functional capabilities, the former in line with disruptive selection and the latter due to potential trade-offs incurred by maternal mouthbrooding in Malawian haplochromines. We revealed evidence of structural adaptations between Tropheops ‘Red Cheek’ and Labeotrophues fuelleborni that impacted the dispersion of mechanical stress in ways that matched the foraging of these species. Also, hybrids showed higher stresses relative to both species across the mandible. Sexual dimorphism in stress handling was evident despite minor differences in shape with males showing enhanced load resistance. However, in hybrids it appeared that males were disadvantaged relative to females, and displayed asymmetry in load handling. Together, these results show evidence of species and sex based biomechanical variation, that could be targeted by divergent selection.


Introduction
The form of an organism often reflects its function and ecology (Schluter, 2000;Skúlason et al., 2019;Wainwright, 1996). Adaptive radiations, which result in a multitude of species diverging from a common ancestor to fit different ecological niches, provide an ideal context to investigate form/function relationships (Schluter, 2000;Skúlason et al., 2019). These relationships can be driven by competition for food, habitat, and other resources that enacts selection to change form/function relationships (Cooper et al., 2011a;Schluter & McPhail, 1993;Skúlason & Smith, 1995). However, phenotypes are often multipurpose causing interactions with other factors such as mate competition, sexual conflict, and parental care (Cooper et al., 2011b;Parsons et al., 2015). The observation that divergence is the consequence of cascading effects that build from a range of characters to precipitate speciation is vital but often overlooked (Foster et al., 1998). Consequently, while divergence in trophic morphology has been suggested as key to adaptive radiation (Highham et al., 2016;Hulsey et al., 2007;Streelman & Danley, 2003) the ultimate outcomes of form and function can be 1 3 driven by multiple processes and the potential trade-offs they incur (Cooper et al., 2011b;Skúlason et al., 2019). In other words, multiple functional needs may not be simultaneously optimized by selection and there may be a sequential order of change during the formation of adaptive phenotypes (Cooper et al., 2011b;Danley & Kocher, 2001;Foster et al., 1998).
Adaptive radiations are often investigated through phenotypic and functional differences between distinct species (e.g., Cooper et al., 2017;Holzman & Hulsey, 2017;Mitchell et al., 2018;Smith et al., 2015). While such studies are useful, it is important to note that adaptive divergence should normally occur through a process that essentially causes continuous phenotypic variation to become discontinuous through the removal of intermediates (i.e. disruptive selection). Following this, selective maintenance is required for divergence to persist over generations, especially when reproductive isolation between species is incomplete, and hybridization can occur regularly. While hybridization can cause adaptive divergence to collapse (Taylor et al., 2005), it can create a number of functionally novel transgressive phenotypes useful for filling unoccupied niches (Holzman & Hulsey, 2017;Parsons et al., 2011), as well as phenotypes intermediate between species that should be functionally inferior (Highham et al., 2016). Support for this idea comes from field experiments which generally show reduced ecological performance within phenotypic intermediates (Martin & Pfenning, 2012;Robinson et al., 1996;Schluter, 2000;Schluter & McPhail, 1993;Svanback & Bolnick, 2007) but a functional perspective on this phenomenon is less clear. Previous research has suggested that relationships between form and function may not be as strong as is commonly believed with different phenotypes providing equivalent performance through many-to-one mapping (Wainwright et al., 2005). Therefore, obtaining intermediate phenotypes for performance testing could provide key comparative opportunities, but this is challenging for many cases of adaptive radiation because reproductive isolation is essentially complete with intermediates unavailable.
Additionally, other finer scales of inquiry are needed to provide a more complete understanding for how function evolves at different levels. As mentioned, functional trophic variation at the species level may not be the only driver of adaptive change. Sexual differences in trophic morphology are known to occur both within, and as a form of adaptive divergence unto itself (Shine, 1989). Through a phenomenon known as ecological sexual dimorphism (ESD) adaptive variation can arise between sexes (McGee et al., 2013;Selander, 1972;Shine, 1989;Temeles & Roberts, 1993;Temeles et al., 2010). While ESD and conventional adaptive divergence are predicted to be exclusive outcomes because they can provide equal fitness benefits (Bolnick & Doebeli, 2003) evidence suggests that patterns of ESD and adaptive radiation can actually be aligned . Indeed, evidence from African cichlids shows that the genetic basis of variation in trophic traits can differ between sexes (Parnell et al., 2012) suggesting that there is potential for ESD to contribute toward overall patterns of adaptive radiation. Given that sex-based biomechanical variation is rarely investigated (but see McGee et al., 2013) there could be significant gaps in our understanding of how morphological function evolves within adaptive radiations.
For vertebrates the craniofacial apparatus often forms a direct physical link between the organism and its prey intake. Indeed, the shape of the craniofacial apparatus is a useful predictor of function, feeding ecology, and performance (Albertson et al., 2005;Grant, 1986;Holzman & Hulsey, 2017;Mitchell et al., 2018;Wainwright et al., 2004;Westneat, 1995). The mandible, as part of the feeding apparatus, is among the first structures to make contact with prey and its form can directly limit prey shapes and sizes as well as the movements required for capture (Westneat, 1995). For example, in fish short, wide mandibles which also typically have heavy bone deposition facilitates biting, whereas suction feeding, where water and food are drawn into the buccal cavity, is associated with a long mandible (out-lever) relative to the ascending arm (opening in-lever) (Albertson et al., 2005;Holzman & Hulsey, 2017). These shapes should relate to the mechanical advantage of the mandible, measured as the ratios between the out-lever and opening in-lever (e.g. Albertson et al., 2005;Hulsey et al., 2010;Wainwright & Richard, 1995;Wainwright et al., 2004;Westneat, 1995). A higher mechanical advantage usually occurs with a short out lever relative to the opening in lever. This in part reflects the degree of muscle attachment to the ascending arm of the mandible, and thus the amount of force that can be applied to it by the adductor mandibulae muscle. A short out lever, but relatively long opening in lever, should result in strong but slow movements for biting (often used for benthic bottom feeding), while conversely, suction feeding requires a low mechanical advantage typically associated with a long outlever relative to opening in lever for rapid movements (often used for limnetic feeding) (Albertson et al., 2005;Hulsey et al., 2010;Westneat, 1995). This understanding of mandible function comes largely from investigation of outward anatomical variation (i.e. lever mechanics) between species with very different ecologies. Given that adaptive radiations should occur through selection acting upon a range of scales (between species, against intermediate phenotypes, between sexes) this provides a need for finer scale investigation of this complex structure with deeper consideration of internal biomechanics.
To this end, the adaptive radiation of African cichlids from Lake Malawi provide an exceptional model to examine functional variation across a range of scales. Across species they show extensive craniofacial variation related to foraging ecology including variation in their mandibles (Hulsey et al., 2010;Parsons et al., 2012;Ribbink et al., 1983). While the major axis of divergence in these cichlids is characterized by morphologies related to suction and biting modes of feeding (Cooper et al., 2010), it is very likely that finerscale divergence has occurred to allow species with similar ecologies to coexist (Ribbink et al., 1983). The functional significance of the mandible for prey capture may also be particularly enhanced in cichlids as its evolution is often posited to be decoupled from food processing which primarily occurs in their pharyngeal jaws (Hulsey et al., 2006;Liem, 1973). Their recent divergence and shared genetic background allow cichlids to hybridize to create transgressive but also phenotypes that are intermediate between species (Albertson et al., 2005;Parsons et al., 2011;Seehausen, 2004). Such intermediate phenotypes comprise the majority of individuals in hybrid crosses and allow for the rare ability to perform functional investigations where reduced performance in intermediates relative to parental species could be expected (Higham et al., 2016).
Building on previous research our understanding of morphological function now has the potential to be greatly enhanced through the adoption of newer imaging techniques (Rayfield, 2007). This could be especially helpful for investigating finer-scale levels of variation where more sensitive tools could be insightful. For example, micro-computed tomography (µCT) is now allowing biologists to examine the internal structure of anatomical elements with a high degree of three-dimensional detail (Polly et al., 2016). Coupled with finite element analysis (FEA), an engineering method for modelling the transmission of mechanical loads within complex biological structures (Polly et al., 2016;Ross, 2005). This combination of techniques has been applied to different ecological scenarios such as the ability to resist bites from predators (e.g. Rivera & Stayton, 2011), biting variation (e.g. Dumont et al., 2005;Pierce et al., 2009) and structural adaptations for feeding (e.g. Hulsey et al., 2008). Validation studies have shown that carefully designed FEA models that incorporate knowledge of specific biology can be tentatively trusted to predict the mechanical behaviour of organismal structures, even with simplified model conditions (Bright & Gröning, 2011;Metzger et al., 2005;Stayton, 2009Stayton, , 2018. This holds promise for FEA but to our knowledge these detailed approaches have yet to be specifically applied toward intermediate phenotypes, nor used to understand the functional consequences of sexual dimorphism, within the context of an adaptive radiation.
Here we examine variation in mechanical loading of the mandible between African cichlids at a relatively fine scale. Previous research between cichlid species has already revealed biomechanical variation that corresponds with foraging tactics (Cooper et al., 2011a(Cooper et al., , 2011bHolzman & Hulsey, 2017;Parsons et al., 2012) but we aim to add new insights through investigation of additional levels of variation. While we will deal at the level of species, the two we focus on have quite similar ecological niches (invoking many-to-one mapping from an ecological perspective). We also include investigation of hybrid intermediates between these species, and variation between sexes. We predict (1) that intermediate phenotypes should show reduced function (in terms of load handling) relative to parental species. This would provide a proximate mechanism for divergent selection that complements and reinforces previous ecological performance studies (Grant & Grant, 1992;Hatfield & Schluter, 1999;Mayr, 1963). At a finer scale we also predict, (2) that males will show functionally superior biomechanical abilities relative to females within species. This is suggested by previous research showing that male cichlids possess steeper craniofacial profiles relative to females, in line with a morphology associated with the generation of power for feeding and defence of territories . Finally, we predict (3) that the biomechanical properties of intermediate hybrid phenotypes will also differ between sexes. This is based on previous findings suggesting that the genetic basis of cichlids mandibles shows evidence of sex linkage (Parsons et al., 2012. If supported these factors would suggest a means for different levels of variation to interact during adaptive radiation and provide an inroad for further functional investigations that focus on variation below the species level.

Study Species
To test how cichlid mandibles cope with the external compressive loads we focused on two species of Malawi cichlids: Tropheops "Red Cheek" (TRC) and Labeotropheus fuelleborni (LF), their hybrids, and different sexes for each. LF is an algae-scraping specialist which performs a stereotypic biting mode of feeding to remove algae from rock surfaces (Ribbink et al., 1983). Functional investigation of this species based on lever mechanics suggests that is has a low degree of kinematic transmission relative to other species in the radiation, which should correspond to a low ability to produce suction force (Holzman & Hulsey, 2017). Indeed, among a large sample comprising 80% of cichlid genera from the Malawi, Tanganyikan, and Victoria radiations LF is the species possessing the shortest jaws, forming the extreme boundary of morphospace (Cooper et al., 2010).
TRC is also an algal grazer and is the immediate neighbour of LF in morphospace in terms of short jaws (Cooper et al., 2010), but it also has a much narrower mandible relative to LF (Parsons et al., 2011). However, this species feeds by plucking targeted threads of algae rather than the scraping mode employed by LF (Albertson, 2008;Albertson & Pauers, 2018). Lever mechanic measures of the mandible suggest that Tropheops species differ functionally from LF by possessing a higher KT coefficient. This suggests that while these species would be considered 'biters' suction forces are used to draw food into its buccal cavity, which could aid in this type of targeted foraging (Holzman & Hulsey, 2017). LF has repeatedly shown an ability to hybridize with other Malawi species in lab studies (Albertson et al., 2005;Parsons et al., 2015). Hybridization in Malawi cichlids is thought to be pervasive, has been observed in the lake (Streelman & Danley, 2003), and is considered part of the evolutionary history for Malawi cichlids (Malinsky et al., 2018). Given that LF is a cosmopolitan species within the rocky areas of the lake it should be in regular contact with TRC raising the likelihood of hybridization events. Indeed, both AFLP and mtDNA data show support for TRC and LF belonging to the same major phylogenetic branch of 'mbuna' rock-dwelling species within Malawi suggesting a close relationship (Joyce et al., 2011). While whole genome sequencing data confirms this relationship suggesting introgression has occurred recently even among distantly related genera (Malinsky et al., 2018).
In line with ESD, for each of LF and TRC there are suggestions of sex differences in foraging locations as males are highly territorial and remain higher on rocky shoals (Ribbink et al., 1983). In addition, males use their mouths for biting and aggressive interactions, including locking their jaws together (Danley, 2011;Ribbink et al., 1983). This is in contrast to the functional constraints that could be imposed on females through mouthbrooding.

Specimen Collection and Preparation for 3D Morphometrics
To identify representative individuals for FEA we first performed a 3D morphometric approach. To quantify variation in mandible shape between LF, TRC, their hybrids, and sexes a high-resolution 3D morphometric approach was conducted. Ten adults for each of LF and TRC were selected from aquarium populations held at the University of Glasgow and sacrificed with an overdose of benzocaine solution following UK Home Office guidelines. Fish were then frozen at − 20 °C until dissected. Samples were comprised of an equal number of each sex confirmed by colouration, internal dissection of gonads and vent size (Moore & Roberts, 2017). The mandible was disarticulated from the maxilla and the quadrate of the skull and cleaned of soft tissue before storage in a 1% phosphate buffered saline solution (PBS). Similarly, to enable the identification of intermediate phenotypes, 120 F2 hybrid mandibles were obtained from a previous cross between LF and TRC (Albertson et al., 2014;Parsons et al., 2015).
To assess mandible shape, µ-CT scanning was conducted using a Bruker Skyscanner 1172 (Bruker, Billerica MA). The µ-CT scanning parameters were 70 kV and 10 µm resolution using the largest camera with a width of 4000 pixels for all mandibles. Post-processing of the scans and registration was carried out using NRecon (v. 1.6.9.18). One TRC model was removed from the sample due to poor quality, leaving 9 specimens (female = 4 and male = 5), while for LF there were 10 specimens (female = 5 and male = 5). Also, 120 F2 hybrids were scanned for 3D morphometrics to inform the choice of intermediate specimens for FEA (details for this approach below). While some F2s would potentially show segregation our previous research on F2 head morphology demonstrated that (even with transgressive individuals) the majority of individuals displayed phenotypes intermediate between the parental species (Parsons et al., 2011).

3D Model Generation for Morphometrics
To create 3D models for surface morphometrics, ScanIP (Version 7.0) was used on the processed µ-CT images. To reduce noise and smooth the model surface, the "Recursive Gaussian" filter was used. Each model was then segmented using the "Interactive Threshold" function to highlight regions of bone based on greyscale values from the scanning process, followed by the "Flood Fill" function which created a mask that included only connected areas of the model to remove scanning artefacts. The final steps of model creation involved smoothing the surface topology and reducing the number of triangles before exporting as an STL (stereolithography) file which then converted to a PLY file to enable 3D morphometrics using MeshLab (meshlab.sourceforge. net).

Morphometrics and Analysis for Species Differences and Representative Specimen Selection
Mandible shape for TRC, LF, and hybrids were quantified using 3D morphometrics; all analyses were conducted using R version 3.4.1 (R Core Team, 2017) using the geomorph package (Adams & Otárola-Castillo, 2013;Adams et al., 2019) unless otherwise stated. Firstly, using Landmark Editor (graphics.idav.ucdavis.edu), 20 homologous landmarks, including those that would be used to measure traditional lever mechanic measures, were manually placed on each model surface (e.g. Albertson & Kocher, 2001;Parsons et al., 2012) (Fig. 1). Prior to analysis in R, this landmark data (in pts format), was converted to a TPS file format using Simple3D (Zelditch et al., 2012a(Zelditch et al., , 2012b. Shape analysis began with a Procrustes superimposition which rotated, translated and scaled landmark data to a common size using the gpagen function. To ensure that downstream analysis of shape only encompassed the bilaterally symmetric component of shape variation, asymmetry was also removed using a Procrustes ANOVA from the bilat.symmetry function (Klingenberg, 2015). Also, because the two different species displayed different allometric trajectories for mandible shape (tested by comparing regression slopes using proc.allometry), no allometric correction was performed (Klingenberg, 2016). This meant that some of our results could be attributable to species-specific and common allometric variation, but the potential for this was reduced by similarity in the size of selected specimens.
As the number of hybrid samples greatly outnumbered those of pure species a statistical model including both species and hybrids would be unbalanced (making effect size calculations unreliable). Such a model could therefore be dominated by hybrid variation, thus, as a first principle we aimed to simply determine whether species differences in mandible shape existed irrespective of hybrids. Such evidence of species differences in mandible shape would in turn justify the identification of hybrid intermediates. Similarly, testing for sex differences across species, and determining whether such differences interacted with species effects was viewed as a step to establish an understanding of sexual dimorphism, and whether it differed between species, prior to broader testing that included hybrids.
Therefore, direct comparisons (excluding hybrids) were done to identify the effect of species and sex on mandible variation. To start, a grouping variable designating TRC and LF was modelled with sex using a Procrustes ANOVA on landmark coordinates. Following this, we further specified the effects of species and sex on shape through the use of two separate discriminant function analyses (DFA) on our LF and TRC samples. To begin, a principal component analysis (PCA) using plotTangentSpace was conducted on Procrustes coordinates. However, because of the small sample size relative to the number of variables we reduced the data set to the first 5 PCs (comprising 90% of the variation) in line with guidelines from Zelditch and Swiderski (2018). Each of sex and species-based DFAs on these PCs was conducted using lda from the MASS package (Venables & Ripley, 2002). Similarly, we performed an additional DFA on F2 hybrids using sex as a grouping variable, but using the full set of PCs obtained. Finally, differences in mandible morphology relative to DFA axes were predicted using a linear least squares regression and visually depicted as landmark points using the shape.predictor, and plotRefToTarget functions.

Representative Specimen Selection for Finite Element Analysis
Due to computational demands our FEA analysis is often limited to representative individuals. In our case, given the multiple levels of variation under investigation, we identified representative individuals from our shape analysis. With knowledge of mandible shapes representative of the main shape differences between parental species from the analyses above, we now examined F2 hybrid mandibles within a new DFA model that now simultaneously included TRC, LF, and F2 hybrids. This DFA provided two sets of canonical root scores which best separated the two species and hybrids. Using the mean values from these scores we then selected a total of six specimens (nearest to these means) including a single male and female mandible from each of TRC, LF, and intermediate F2 hybrids. We specifically adopted this DFA based approach because it maximized differences between our groups. In other words, for our choice of representative specimens this was intended to minimize interindividual aberrations and identify individuals which best exemplified the differences among groups. However, while most F2 hybrids were expected to intermediate between species we confirmed this through pairwise comparisons of Euclidean distances between average TRC, LF, and F2 hybrids using dist. Variation in mandible morphology relative to the two DFA axes were again predicted the shape.predictor, and plotRefToTarget functions as above.

Mesh Creation and Finite Element Analysis
To enable finite element modelling the 3D models for selected specimens were adjusted for assessment through additional steps (Panagiotopoulou, 2009;Peterson & Müller, 2018;Richmond et al., 2005;Ross, 2005). For each specimen meshes were created from the original geometry of the 3D model following a voxel-based approach using the software ScanIP (Version 7). Four node Linear Tetrahedral elements (C3D4) were used for each model. Elemental shape was checked visually and by using the mean in-out and edge length aspect ratios of elements that were calculated to ensure mesh quality (Stayton, 2009). Prior investigation in cichlids has shown that bone density does not differ between LF and a similar species . Therefore, we assigned each model element to homogenous material properties (Young's modulus and Poisson's ratio). To reduce computational demands, this also included assigning the teeth the same material properties which is common practice for vertebrate jaws (Dumont et al., 2005), but it is noteworthy that some variation in material properties has been observed within cichlid teeth, and between their bone and teeth (Hulsey et al., 2008). Nonetheless, variation due to tooth shape and number were incorporated into our models. The Young's modulus was set for both teeth and bone at 6 GPa (Cohen et al., 2012;Horton & Summers, 2009) and the Poisson's ratio at 0.3 (Dumont et al., 2005;Hulsey et al., 2008). Although asymmetry was removed for morphometric analysis, it was retained for FEA models as it is a natural component of mandible variation that likely influences mechanical loading (Stewart & Albertson, 2010).
To examine variation in the ability of the mandible to cope with loading, each FEA model was tested with compressive loading on the teeth and the transfer of forces to the jaw at levels likely present during the moment of initial contact through scraping and plucking of algae. Therefore, bite forces originating directly from the muscles attached to the mandible were excluded from our FEA model because they should be initiated following the first contact with prey. Further, our models were not intended to mimic the rotation of the jaw on the articular joint with the quadrate. Thus, the boundary conditions for FEA models were assigned along each ascending arm of the mandible using nodal forces where the adductor mandibulae attach, and immediately dorsal to the articular joint (Fig. 2). In other words, our models were intended to mimic a stable stationary jaw the moment is comes into contact with a surface, with the adductor mandibulae bracing, but not yet closing the jaw.
Evidence suggests that mbuna cichlids bias their biting to one side (Stewart & Albertson, 2010). Therefore, to provide realism the left side of the mandible was constrained in all three directions (x, y and z) whereas the right side was constrained along y and z so as to allow lateral movement for flexibility and deformation during loading (Fig. 2). To reflect species variation due to the different foraging tactics that LF and TRC employ, a range of loading scenarios were applied to all meshes. Specifically, during algal scraping LFs probably use the full width (75-100%) of the mandible during scraping due to their foraging tactics. Unlike most other algae eating mbuna species which use a bite-and-jerk/ twist mode (Konings, 2007;Ribbink et al., 1983), Labeotropheus species employ a bite-and-roll mode of feeding whereby they protrude their upper jaw ventrally toward the rocky substratum, take bites of tough filamentous algae, and then roll forward onto their snouts (serving as a fulcrum) while retracting their upper jaw to allowing algae to be leveraged off the rocks by their lower jaw (Concannon & Albertson, 2015;Konings, 2007; KM personal observation). This foraging mode likely places more overall force on the lower jaw, but during initial contact it is likely spread over the very wide and square-shaped mandibles of Labeotropheus which provide more contact (~ 75-100% in contact). Conversely, for TRC the more typical bite-and-jerk/twist foraging mode it uses should place more focused forces on the anterior midline of the mandible during initial contact. This is likely why the TRC mandible is overall much rounder and narrower than in LF with only a small proportion of the width (25-50%) likely to be in surface contact when plucking targeted algal threads. Therefore, to compare species our FEA models included force being applied equally across four different widths of the mandible (25%, 50%, 75% and 100%) (Fig. 2) in an effort to reciprocally mimic their foraging tactics. The FEA was carried out for each selected individual using Abaqus (v. 6.13, Simulia).
After preliminary tests on our models a 1 N compressive force was chosen to simulate external loading on the mandible, then scaled depending on the surface area of the mesh to accurately compare performance differences attributable to shape rather than size. This followed the recommendations of Dumont et al. (2009), ensuring that the applied force to surface area ratio was constant. Therefore, the applied load from the model with the lowest surface area was used as a reference, and the force for each model was scaled using Eq. 1: where SA is the surface area of the reference (SA A ) and target (SA B ) models, F A is the force for the reference model (1 N), and F B is the force applied to the target model. Surface area, rather than volume, was used as the jaws differ in shape as well as size (Dumont et al., 2009). Because a set of comparisons was the intention, the actual loading values are less important than the stress patterns (Rayfield, 2007) with visualization of von Mises stresses being considered a reliable predictor of fracture (Dumont et al., 2009).

3D Morphometrics
Mandible shape differed between species, as well as between sexes. However, there was no interaction between these factors which suggested that sexual dimorphism was similar in both species (Table 1). Species explained a much greater degree variation in mandible shape (57%) compared to sex (7%), and the basis of R-squared, although this was less distinct on the bases of Z-scores (Table 1). A priori groupings of species were also supported as the DFA showed a correct classification rate of 100%. Evidence for sexual dimorphism was also strong across species with the combined correct classification rate being 100% and 89% for males  and females, respectively across groups (Fig. 3). The major shape difference between all groups appeared to be width, with LF mandibles appearing wider than TRC with the angle between the left and right sides of the jaw being larger (Fig. 3a). Similarly, males appeared to have wider mandibles than females with more widely angled left and right sides (Fig. 3b). Finally, the DFA including F2 hybrids and both species showed a correct classification rate for all groups, with F2 hybrids appearing to be intermediate between species on axis 2 which appeared to be driven by differences in jaw width (Fig. 4). The intermediate morphology our hybrids was confirmed with Euclidian distance being larger between the consensus forms of species (10.3), than between TRC and F2s (8.9), and LF and F2 (8.1) consensus forms.

Finite Element Analysis
Model meshes for all specimens were of sufficient quality according to previously published values (Gislason et al., 2010;Ito et al., 2006). Mesh statistics and quality parameters (Mean In-Out Aspect ratio and the Mean Edge Length Aspect Ratio) can be seen in Table 2.
Considering the foraging modes of each parental species the models focused on mimicking the moment of first contact with the surface, and showed that TRCs had a greater area of high loads than the LF when loading was applied across 100% and 75% of the jaw width (Figs. 5,6,Suppl. Figs. S1,S2). This suggested that LF foraging tactics applied to TRC mandible models caused detrimental loading outcomes. Additionally, all four species models showed stress laterally on the mandible to differing degrees, particularly at the edge of the articular web where the values of stress were 10 MPa and higher. The LF models showed relatively low stress across the dentary region of the mandible whereas the TRC models showed slightly higher stress values across this region when loaded as an LF (i.e. 75-100% loading width). Notably, the lateral projections or "wings" of the tooth-bearing dentary region unique to the LF mandible showed minimal stresses when force was loaded across 100% and 75% of the mandible width (Figs. 5, 6). While these wings could simply be used to increase the surface area of the mandible to increase prey intake, this finding suggested that shape of LF mandibles could also provide the ability to efficiently dissipate loading across the full width. Placing the loads across 50% and 25% of the width of the mandible simulated how a TRC would apply loading during feeding (Figs. 2, 5, Suppl. Fig. S2). Under these conditions, the TRCs showed higher stress values across the dentary region and laterally on the jaw but appeared to have minimal stress across the lateral midline. This midline "ridge" area of low stress at the suture between the left and right sides of the mandible was not present to the same extent in LFs. Comparatively, LFs exhibited higher stress laterally and in the dentary region than for the simulated LF scenarios, but the stress patterns in this area were not as intense as for TRC (Fig. 6, Suppl. Fig. S1). Relative to LFs and TRC males, F2 hybrids showed higher stresses across the articular web on both sides of the mandible for all loading scenarios and the stress patterns for each scenario for both sexes of hybrids were similar to the TRC female (Figs. 5, 6).
Females generally experienced higher loading stresses than males, this pattern was particularly evident in the TRC female, which showed > 7 MPa across the mandible for all loading scenarios. Both sexes were similar in how they dissipate stress with both LF males and females having low levels of stress across the dentary region and the "wings" during loading scenarios using 75% and 100% of the jaw width. Also, both TRC males and females had a reinforced Fig. 4 Scatterplot depicting the scores from a discriminant function analysis (DFA) on mandible shape used to aid specimen selection for finite element analysis (FEA). Tropheops "Red Cheek" (TRC) is depicted in green, Labeotropheus fuelleborni (LF) in pink, and F2 hybrids are in yellow. The shape changes that occur across each dis-criminant function is depicted for each axis using landmark points with × 5 magnification. The inset box indicates a magnified view of the F2 hybrid specimens; the two circles in the magnified box indicate the models which were selected for the FEA (Color figure online) "ridge" which likely dissipated stress along the midline. However, the intermediate hybrid female showed very low stress across the dentary region when loading was across the full width of the jaw but displayed a similar stress pattern to the TRC male for the loading scenarios where the force was applied at 50% and 25% of the width. Also, within intermediate hybrids the female showed lower stress values across the midline whereas the hybrid male showed notable asymmetry in stress patterns across the mandible.

Discussion
The ability of organisms to evolve to accommodate functional demands should explain much about adaptive radiations. However, functional data from multiple levels of organization (species, hybrids, sexes) are rare within a single study, but as we have suggested, could lead to insights about the fine scale processes that occur along a continuum of variation, and perhaps even sequentially, during adaptive evolution (Foster et al., 1998;Parsons et al., 2015). Our data provides evidence that species, hybrids, and sexes, differ in biomechanical properties. The specific nature of these differences provides support for how selection could target functional variation to favour and maintain species differences (Higham et al., 2016). Further, our data suggest functional differences could occur through sexual dimorphism, making it an avenue for adaptive divergence. Below we discuss how each level we investigated could provide implications for adaptive divergence.

Species Differences
Our FEA models provided evidence the morphology of cichlid mandibles corresponds with ecological function. Specifically, LF possessed a much wider mandible than TRC with FEA models indicating a greater ability to dissipate mechanical loading as evidenced by relatively low levels of maximum stress (Fig. 4; Albertson & Pauers, 2018). This was evident at the anterior of the mandible near the toothbearing region where LF models did not display stresses much above 2.5 MPa, whereas TRC showed stresses above 5 MPa extending well into this region (Fig. 4). Such differences matched our expectations that the mode of feeding used by LF (bite-and-roll) would require an enhanced ability to handle mechanical loads at the point of contact with food (Concannon & Albertson, 2015). Overall, species differences accounted for the majority of mandible shape variation, while FEA models also showed relatively greater differences between species than other levels ( Table 1, Figs. 4,5). Therefore, this suggested that shape evolution could be most strongly driven by species-specific foraging tactics and their inherent biomechanical demands.
Species-specific attributes of mandible morphology based on lever mechanics suggest that function and stress handling abilities play a strong role across a wide range of cichlid species (Holzman & Hulsey, 2017;Martinez & Sparks, 2017), but our data suggests that unique finer scale structures, with biomechanical variation captured by our FEA approaches, could also play a role. In particular, the lateral expansions of the dentary ("wings") of LF, would seem to be a novel trait that aids in the dissipation of stress (Fig. 5). In homologous regions of the TRC mandible stress levels were higher suggesting that these lateral expansions of the dentary (which are relatively flat across their width) achieve reduced stress by simply increasing the amount of surface area in contact during biting. Conversely, TRC with a rounded anterior tip of the dentary, should make less surface contact which is likely suited toward targeted algae plucking by focusing more bite force on the midline. Indeed, the apparent reinforcement and overall increase in bone volume in the midline of TRC relative to LF seems to dissipate stress before being transferred laterally (Fig. 6). This trait would seem to result from developmental changes in how the suture between the left and right sides of the dentary form. Notably, several FGF receptors have been implicated for the induction of bone fusion and subsequent enhanced but localized bone growth at suture points in coordination with BMPs (Kim et al., 1998;Rice et al., 2003) and should provide tractable targets for future developmental investigation.
Biomechanical differences in the mandible likely occurs commonly for cases of adaptive divergence. For example, in postglacial fishes morphological variation in the jaw corresponds with dietary differences between benthic and limnetic phenotypes, which respectively represent habitats that should favour biting and suction modes of feeding (Robinson et al., 1996;Schluter & McPhail, 1993). This habitat axis appears to drive the most major patterns for both postglacial fishes and African cichlids (Cooper et al., Skúlason et al., 2019) with biting being employed by both LF and TRC, although it would appear that TRC also employs aspects of suction feeding as inferred from our data and previous study (Holzman & Hulsey, 2017). Indeed, adaptive variation is also common off this major axis (Skúlason et al., 2019), especially within cichlid radiations, but the drivers for this are far less understood. FEA could provide such insights into specific mechanical implications of shape. In particular, the LF mandible, relative to TRC, showed a greater propensity to dissipate stress under the wider loads mimicking their foraging mode (see Fig. 5 whereby Von Mises values in the LF dentary and articular web are markedly lower than in TRC). Conversely, under the narrower loads mimicking TRC foraging tactics, we observed a marked ability to dissipate stress through the aforementioned midline of the TRC mandible relative to LF (Fig. 6). However, some relative difficulties in dealing with stress seemed to be imposed by reciprocally applying FEA species models. While TRC exposed the toothbearing region to more stress under all model scenarios, LF was particularly susceptible to midline stress under

Sexual Dimorphism in the Biomechanics of the Mandible
Sexual dimorphism was evident through reduced female performance relative to males (Figs. 5, 6). Notably, sexual dimorphism in handling stress seemed to be enhanced by matching the loading scenario toward a given species. For example, for the 100% (LF mimic) loading scenario LF males appeared more able to dissipate stress than LF females, while for the 50% (TRC mimic) sexual differences were less pronounced in terms of the regions facing stress (although stress levels were still lower in males) (Fig. 5).
Similarly, the 50% loading tests showed increased sexual differences in TRC, but reduced differences in LF (Fig. 6). This suggests that sexual dimorphism in function should Fig. 6 Results for each of the models using the 50% loading scenario from multiple angles (ventral, left lateral, dorsal, and posterior views). Colours on each mandible model represent gradients of stress from finite element analysis going from 0 (blue) to 15 MPa (red); grey represents the maximum stress value (in MPa). The colour scale for the Von Mises stress is consistent across all models, but the maximum stress value (grey) varies between models. Black arrow indicate the lateral "wings" in Labeotropheus fuelleborni as well as the midline "ridge" present in Tropheops "red cheek" (Color figure online) be enhanced within species, perhaps as a consequence of species-specific foraging tactics. Following this, ESD might therefore be enhanced, or driven by intraspecific factors. Curiously, previous research on cichlids showed that patterns of sexual dimorphism in craniofacial shape from a hybrid cross between TRC and LF were aligned with the main trajectory of divergence across the Malawi radiation . However, while aligned, the degree of differentiation for sex was several orders of magnitude smaller than that found across the radiation. Our FEA results complement these patterns by suggesting that such sexual dimorphism could hold functional ecological relevance. This property may facilitate rapid evolution of adaptive radiations that can coexist with ESD and needs further testing, perhaps through foraging performance assays and FEA model validations that take account of sex (see Stayton, 2018).
However, several reasons could explain sexual dimorphism in mandible functional abilities. As mentioned, males show high levels of aggression and fighting involving the locking of jaws between rival males should increase mechanical loads (Ribbink et al., 1983). For females mouthbrooding likely places different functional constraints on the craniofacial apparatus not faced by males. Indeed, mouthbrooding cichlids from Lake Victoria show that the volume needed for mouth-brooding in females has resulted in a "biting" phenotype (Tkint et al., 2012). Therefore, both reproductive and ecological differences could have additive and interactive effects on adaptation (Herler et al., 2010;Oliveira & Almada, 1995).

Hybrid Performance
Through asymmetry of mechanical stress dispersion, and areas of extreme stress (especially in the articular web) our FEA models suggested reduced hybrid performance (Figs. 5, 6). Thus, functional variation within intermediate morphologies likely provides a target for disruptive selection. While intermediate hybrid inferiority is a long-standing and often supported concept (Abrams et al., 2006;Arnegard et al. 2014;Burke & Arnold, 2001), it has often been the focus of genetic studies rather than for its role in biomechanical variation. Through intermediate hybrids the mixing of two developmental systems could mismatch form and function (Higham et al., 2016). Here, the higher stresses modelled in the hybrid dentary and articular web relative to parentals represent key areas of contact with prey, and musclegenerated force by the adductor mandibulae suggesting a decoupling of form and function (Parsons et al., 2012).
The male hybrid incurred higher stresses than the female (Figs. 5, 6) which was particularly prevalent on the left articular web. While this individual was representative of an intermediate form this pattern of stress may simply reflect particular aspects of this specimen. However, this could also be reflective of how the inheritance of speciesspecific functional features occur. For example, the hybrid female showed evidence for the presence of both a midline "ridge", which is prominent in TRC and protrusions like the lateral dentary extensions or "wings" present and specific to LF. These features were not as pronounced in the male. Thus, it may be that some variation occurs through an epistatic interaction with sex determining loci (Parnell et al., 2012). Indeed, findings from several QTL studies in cichlids show a strong tendency for sex linkage in a range of traits (Albertson et al., 2014;Parsons et al., 2015Parsons et al., , 2018) that could be impacting functional variation. Further research is needed but we predict that some trait combinations will not be possible for both sexes.

Developmental Implications for Mandible Biomechanics
As mentioned, the extended width of the LF mandible is facilitated by lateral bony outgrowths from the dentary ("wings") that harbour teeth and increase the amount of surface area contact. Our models suggest that this could enhance load handling in LF relative to the narrower TRC (Figs. 5, 6). The unique foraging tactics of LF have likely promoted the formation of other novel features, including an exaggerated snout (serving as a fulcrum during the biteand-roll motion) comprised of a hypertrophied intermaxillary ligament, which connects the right and left sides of the upper jaw, and the overlying loose connective tissue (Conith et al., 2018). We suggest that both of these novel traits are likely to be functionally and developmentally integrated as both mandible shape and snout variation have been linked to the transforming growth factor β pathway (Albertson, 2008;Conith et al., 2018).
Load handling abilities may contribute to novel anatomical features by 'shielding' development from environmental influence. Indeed, mechanical loading on bone induces its remodelling through Wolff's Law, or more generally phenotypic plasticity (Wolff, 1892). Given that a biting tactic should create more external stress on the mandible, it could be expected that LF would experience more bone remodelling . However, with their efficient ability to dissipate external loading from feeding, remodelling might be negated in LF and could explain why they possess lower levels of plasticity in craniofacial traits relative to TRC (Navon et al., 2020;Parsons et al., 2015). Therefore, we suggest that plasticity, which is often predicted to be reduced through genetic assimilation as adaptive divergence progresses (Skúlason et al., 2019), could also be linked to the optimization of anatomical function.

Conclusions
Using a series of FEA model conditions mimicking the different foraging tactics of two species we have gained insights into how form translates into function at multiple biological levels. While validation of these models through a range of approaches could aid in further supporting our interpretations our findings indicate that finer-scale phenotypic variation could be an important target of selection that contributes to broader patterns of divergence. Indeed, even when outcomes appear to be clearcut (i.e. species) there are likely multiple avenues for function to change and even facilitate further adaptive radiation.