Internal architecture of the mandibular condyle of rabbits is related to dietary resistance during growth

Plasticity allows an organism to ‘fine-tune’ its form to best fit its behavior and is thought to be important at the macroevolutionary scale (West-Eberhard, 1989, 2005; Scheiner, 1993; Agrawal, 2001; Dewitt and Scheiner, 2004; Pigliucci et al., 2006). One key component of phenotypic plasticity and morphological variability is bone functional adaptation (i.e. the modern generalization of Wolff's law; Ruff et al., 2006), succinctly defined as ‘over time, the mechanical load applied to living bone influences the structure of bone tissue’ (Cowin, 2001, pp. 30–31). This relationship between loading environment and bone tissue is well established and supported by an array of data from across vertebrates and in relation to a variety of factors, including age, health, genetics and hormone levels (Frost, 1987, 2003; Lanyon, 1996; Majumdar et al., 1997; Lieberman et al., 2001; Pearson and Lieberman, 2004; Pontzer et al., 2006; Glatt et al., 2007; Barak et al., 2011; Glass et al., 2016; Ravosa et al., 2016). Further, relationships between loading regimes and bone morphology have been examined across the skeleton and in multiple tissue types (e.g. trabecular versus cortical bone) in both intraspecific and comparative contexts (e.g. Fajardo and Müller, 2001; Griffin et al., 2010; Barak et al., 2013; Patel et al., 2013; Chirchir et al., 2015). Fewer studies have focused on bone functional adaptation, or plasticity more generally, in the masticatory apparatus.

Building on earlier work (Bouvier and Hylander, 1981; Bouvier, 1988), recent studies have highlighted the importance of understanding developmental plasticity in the masticatory apparatus. Diet-manipulation experiments conducted on growing rabbits and rodents demonstrate that masticatory morphology is sensitive to differences in food mechanical properties during development (e.g. Bresin et al., 1999; Taylor et al., 2006; Ravosa et al., 2007; Tanaka et al., 2007; Ravosa et al., 2008; Chen et al., 2009; Menegaz et al., 2009; Ravosa et al., 2010a,b; Scott et al., 2014a,b; Polur et al., 2015; Franks et al., 2016, 2017; Ravosa and Kane, 2017; Dutra et al., 2018), and diet-induced morphological differences within species can reach or exceed those observed between closely related species of living primates (Ravosa et al., 2016). Importantly, these studies consistently demonstrate that the bony dimensions and internal bone volumes of the masticatory apparatus decrease in relation to decreased loading (e.g. Chen et al., 2009; Rafferty et al., 2012; Scott et al., 2014a,b; Balanta-Melo et al., 2019a). Less clear is the importance of the timing of the introduction of mechanically challenging foods (i.e. before versus after weaning and/or skeletal maturity), the long-term (i.e. >16 weeks) effects of variation in loading patterns in the masticatory apparatus, and the role of more naturalistic diets (i.e. beyond disparities in food mechanical properties).

The temporomandibular joint (TMJ) is of particular interest for understanding plasticity and bone functional adaptation in the masticatory apparatus because of its role in accommodating joint reaction forces and in facilitating movements during biting and chewing. In the TMJ, joint reaction forces must be transferred via articular cartilage, subchondral bone and trabecular structure into the cortical shell of the condylar neck and ramus. Increasing condylar articular surface area, increasing biomineralization and/or altering bone structure to have more and/or thicker trabeculae will further help to withstand these forces. Work by Ravosa and colleagues referenced above has documented variation in condylar dimensions among rabbit treatment groups, though the magnitude and significance of differences varied throughout ontogeny and depended on experimental cohort (e.g. Ravosa et al., 2007; Scott et al., 2014a). A lack of observed differences in joint surface area in some experiments has been used to support the hypothesis (Lieberman et al., 2001; Reeves et al., 2016) that joint dimensions are more constrained and less subject to plasticity during development because joint congruence must be retained (but see Congdon et al., 2012). Even in the absence of changes in external joint morphology, internal structure may be altered in relation to loading, such that trabecular and/or cortical dimensions vary between groups because such responses likely have minimal impact on joint congruence.

Condylar trabecular structure has been examined by a number of researchers in humans (Giesen and van Eijden, 2000; van Ruijven et al., 2002; van Eijden et al., 2004, 2006; Kahn et al., 2019) and in a variety of other mammals (Teng and Herring, 1995, 1996; Tanaka et al., 1999; Mulder et al., 2005; Cornish et al., 2006; Willems et al., 2007; Chen et al., 2009; Rafferty et al., 2012; Matthys et al., 2015; Ben-Zvi et al., 2017). This prior research demonstrates that condylar trabecular structure varies across the articular surface, and that trabecular plates running anteroposteriorly in the condyle are most likely to be oriented perpendicular to the articular surface. Studies comparing dentate and edentate humans have documented clear differences in trabecular structure, with edentulous individuals having lower bone volume fractions (amount of bone in a given volume; bone volume/total volume, BV/TV) and lower bone density and elastic modulus than dentate individuals (Giesen et al., 2003, 2004). No differences in trabecular anisotropy (a measure of how highly oriented trabeculae are) or trabecular number were identified. Studies of pigs have also documented variation in trabecular structure among individuals, primarily in relation to development (Mulder et al., 2005; Willems et al., 2007). These analyses found that the amount of trabecular bone, trabecular thickness and bone mineralization increased during development, though directionality of trabeculae changed little. A variety of analyses demonstrate that reduced masticatory loading in rabbits and rodents, induced via a ‘soft’ diet (e.g. pulverized and watered pellets) or by decreased muscle output following botulinum toxin injections, results in lower cortical and trabecular bone volume, thinner trabeculae and fewer trabeculae in a relatively short time frame (e.g. Chen et al., 2009; Tsai et al., 2010; Rafferty et al., 2012; Kün-Darbois et al., 2015; Balanta-Melo et al., 2019b).

Here, we extend previous work on New Zealand white rabbits (Oryctolagus cuniculus) (Scott et al., 2014a,b; Franks et al., 2016, 2017; Ravosa et al., 2016; Ravosa and Kane, 2017) by investigating condylar trabecular structure in an experimental context. These earlier studies examined plasticity in multiple aspects of the masticatory apparatus in this same experimental sample, but none quantified trabecular structure of the mandibular condyle as we do here. Thus, the present study is the first to focus on adult condylar trabecular structure in an experimental sample of mammals raised from weaning into adulthood on diets that varied in their mechanical properties and differed in the timing of introduction of mechanically challenging foods, simulating seasonal variation in diet. We employed a novel method introduced by Sylvester and Terhune (2017) to explore variation in trabecular properties across the condylar articular surface. We hypothesized that trabecular properties will track diet such that rabbits raised on mechanically challenging foods will develop increased bone volume fraction, increased trabecular thickness and decreased trabecular spacing. We also determined how the timing of the introduction of mechanically challenging foods during ontogeny influences adult trabecular properties as well as differences in condylar size and shape among groups. This work addresses outstanding issues such as the role of adaptive plasticity in the masticatory apparatus and the implications of trabecular structure for understanding variation in diet among extant species and inferring dietary behaviors in fossil species.

All procedures used in this study were approved by the University of Notre Dame Institutional Animal Care and Use Committee (IACUC). Our sample originally consisted of male New Zealand white rabbits, Oryctolagus cuniculus (Linnaeus 1758) (N=40), obtained at weaning (age 5 weeks) from Harlan Laboratories (www.harlan.com) and housed at the University of Notre Dame's animal care facility, Friemann Life Science Center. Day-to-day care of the animals, including periodic health evaluations, was handled by trained veterinary staff. The animals were raised for 48 weeks, making them 53 weeks old at the conclusion of the experiment. In white rabbits, weaning typically occurs at 4–5 weeks of age, and skeletal and sexual maturity are attained at ∼26 weeks (Masoud et al., 1986; Isaksson et al., 2010). All rabbits were killed following approved IACUC protocols at the conclusion of the experiment (week 48). Oryctolagus cuniculus resembles other mammalian herbivores in ways that make it a good model organism for investigating diet-related plasticity in the mandibular condyle. The most salient features include: (1) a vertically deep facial skeleton, tall mandibular ramus and TMJ situated high above the occlusal plane; (2) a TMJ that is capable of rotational and translational movements, permitting transverse jaw movements during mastication; (3) intracortical bone remodeling; and (4) well-characterized patterns of covariation among dietary mechanical properties, jaw-muscle activity and jaw-loading regimes (Weijs and de Jongh, 1977; Weijs et al., 1989; Hirano et al., 2000; Langenbach and van Eijden, 2001).

We began the experiment at weaning for three reasons. First, it allowed us to examine how phenotypic sensitivity to environmental stimuli changes during development. Second, it mitigates the effects of other dietary influences that might confound comparisons among treatments. Finally, because mammals begin to adopt adult diets and chewing behaviors around weaning (Herring, 1985; Weijs et al., 1989; Iinuma et al., 1991; Westneat and Hall, 1992), initiation of diet manipulation just after weaning facilitates a more naturalistic experiment (Ravosa et al., 2007). The importance of the last strategy is particularly crucial when coupled with an investigation of dietary responses to protracted alterations in food mechanical properties, a likely condition in the wild.

Upon arrival at the University of Notre Dame, the rabbits were divided equally into four cohorts (N=10 each) (Table 1). Animals in the Control group (specimens ND1–ND10) were fed a diet consisting solely of Purina rabbit pellets throughout the experiment. Animals in the Early group (ND21–ND30) were each given three hay cubes (∼3.2×1.9×1.9 cm) per day, in addition to pellets, for the first 6 weeks of the experimental period and then switched to an all-pellet diet for the next 18 weeks (weeks 7–24); this simulated a seasonally variable diet with early developmental exposure to mechanically challenging foods (i.e. hay). Animals in the Late group (ND11–ND20) were put on the opposite feeding schedule: pellets-only for the first 18 weeks (weeks 1–18), then pellets with three hay cubes daily for the next 6 weeks (weeks 19–24); this simulated a seasonally variable diet with late developmental exposure to mechanically challenging foods. Thus, the Early rabbits consumed hay directly after weaning, whereas the Late rabbits were not exposed to hay until around the time of skeletal maturity. Animals in the Annual group (ND31–ND40) were given pellets and hay throughout the experimental period, with temporal variation in the number of hay cubes received per day (Table 1); this schedule simulated year-round consumption of mechanically challenging foods. The feeding schedules for the Early, Late and Annual groups were repeated in the second half of the experimental period (weeks 25–48) (Table 1). The amount of pellets received by each rabbit was determined by veterinary staff based on standard guidelines for caloric intake. Data on the mechanical properties of rabbit pellets and hay cubes demonstrate that hay cubes are stiffer (3335.6 MPa) and have higher fracture toughness (2759.8 J m⁻²) than pellets (29.2 MPa, 1030.6 J m⁻²) (Ravosa et al., 2007, 2015). Behavioral observations, combined with the lack of strain differences along the mandible between hay and pellet mastication (Weijs and de Jongh, 1977), suggest that greater repetitive loading is the most likely influence on post-weaning bone formation here via increased chewing investment (chews per gram) and chewing duration (Ravosa et al., 2015). Data for all four of these treatment groups were collected simultaneously as part of a single experiment, though prior research with this dataset (e.g. Scott et al., 2014a,b; Franks et al., 2016, 2017; Ravosa et al., 2016; Ravosa and Kane, 2017) has focused on subsets of these groups to answer separate but related research questions.

Two rabbits died prior to completion of the experimental period. Rabbit ND19, a member of the Late group, died during week 16 (21 weeks of age), before receiving its first bout of hay cubes and was thus excluded from all analyses. Rabbit ND7, a member of the control group, died during week 33 (at 38 weeks of age) after achieving adult size but 15 weeks before the end of the experiment. A third rabbit, ND32, assigned to the Annual group, stopped eating soon after the experiment began. This animal subsequently began eating pellets and survived until the end of the experiment, but never ate its ration of hay cubes. Though we retained ND7 here because it reached skeletal maturity by the time it died, we excluded ND32 because it cannot be considered an Annual rabbit. These exclusions gave us a final specimen total of 38 rabbits (Control N=10, Early N=10, Late N=9, Annual N=9).

Following removal of the mandibular condyles, micro-computed tomography (microCT) scans were collected for both sides of each specimen (total condyles for 38 specimens=76). MicroCT scans were conducted on a Scanco VivaCT 80 microCT housed at the University of Notre Dame Integrated Imaging Facility. Specimens were scanned at 70 kV with a slice thickness of 10.4 μm. Data were reconstructed to an in-plane voxel dimension of 10.4 μm and exported as 16-bit DICOM file format. Image stacks were examined for problems that may have occurred during scanning or with the specimen itself; in this process it was determined that the right condyle of ND1 was fractured postmortem. This condyle was excluded from further analysis (N=75).

Using these data, we generated a series of image stacks that allowed us to extract trabecular parameters (Fig. S1), including: (1) a binary stack, (2) a trabecular thickness map and (3) a trabecular spacing map. To generate the binary image stack, CT image stacks were segmented in ImageJ using the default automated threshold algorithm (Ridler and Calvard, 1978). The resulting binary image stack represents bone as white voxels and everything else (air, soft tissue) as black voxels. For three scans (ND7 left, ND8 right, ND34 right), the algorithm failed to converge on a threshold value; consequently, the average threshold value for all other scans was applied.

Next, to create a bone mask image stack from the segmented binary image stack, we employed a 3D morphological dilation of the image stack using a spherical kernel with a radius of three voxels (Blanchet and Charbit, 2006). Following dilation, the images were filled, and we performed a 3D morphological erosion using a spherical kernel of the same size (Blanchet and Charbit, 2006). This process resulted in an image stack in which all voxels internal to the bone surface were white. We then compared the bone mask images with the original binary images to generate a new stack where trabecular spaces inside the bone were white while all other voxels (bone and air outside the condyle) were black.

The BoneJ plugin (Doube et al., 2010) for ImageJ was used to calculate trabecular thickness values for the entire binary bone image stack (i.e. trabecular thickness map), as well as trabecular spacing using the trabecular space image stack (i.e. trabecular spacing map). Note that we chose not to separate the cortical and trabecular bone in these image stacks, as the cortical shell was extremely thin across the condylar articular surface (e.g. Fig. S1), making consistent identification of the boundary between these two types of bone too prone to error (Bryce et al., 2015). Further, both the subchondral bone and the trabecular structure are part of the load path for the mandibular condyle and therefore should be analyzed together. For ease of discussion we refer generally to variation in trabecular architecture throughout the paper but it should be noted that, at least in some locations on the condyle, this likely includes a small portion of cortical and subchondral bone as well. Finally, a surface model (PLY file format) of each mandibular section was generated in Avizo Lite 9.5 (Thermo Fisher Scientific) using the bone mask image stack. This surface model was used for the subsequent geometric morphometric (GM) analyses.

To quantify shape variation in the mandibular condyles, and for use in the multi-volume of interest (VOI) trabecular bone analysis, we carried out a 3D GM analysis of the mandibular condyle subchondral bone surface. Initially, we selected one individual as a template model and distributed 70 sliding semilandmarks across the articular surface and along the surface's edge. The landmarks were connected to make a triangulated mesh model of the joint surfaces, consisting of 28 landmarks along the edge of the articular surface and 42 landmarks distributed across the surface (Fig. 1). For the other specimens in the sample, we used the following procedure. First, we collected 75 landmarks representing the edge of the articular margin. These landmark positions were resampled to have the same number of edge landmarks as the template model. Because the template and the sample specimen had the same number of landmarks, and these landmarks were collected in the same anatomical order, correspondence between template and specimen landmarks could be established. This was used to generate a thin-plate spline (TPS) interpolation function that mapped the template edge coordinates to the exact position of the edge landmarks on the specimen. This TPS function was applied to all 70 template landmark positions, which then aligned the template roughly to the sample specimen. Next, surface normals for the template surface sliding semilandmarks were calculated and used to establish a correspondence with a surface vertex on the sample specimen. Once established, the correspondences were used to calculate a second TPS interpolation function that warped the template to the sample specimen. The warped templates provided starting positions for the sliding process.

For each specimen, landmarks were allowed to slide along tangent planes (surface landmarks) or tangent vectors (curve/edge landmarks) to minimize bending energy of the TPS function relative to a reference specimen (Gunz et al., 2005). For the initial round of sliding, one random specimen was selected to serve as the reference. Because landmarks slide along tangent planes and vectors, some landmarks ultimately slide off the curved surface of the anatomy during the sliding process; thus, following sliding, landmarks were projected back onto the closest anatomical position on the articular surface. This process was repeated until the landmarks attained stable positions (i.e. the bending energy was minimized and the sliding semilandmarks did not move), at which point we undertook a generalized Procrustes analysis (GPA). The resulting shape coordinates were used to estimate the average mandibular condyle shape, which was used as the reference specimen for the next round of sliding. This iterative process – sliding based on the TPS, followed by GPA to establish a new estimate of the average shape – proceeded until the average shape stabilized and the sliding semilandmarks no longer moved when subjected to the sliding algorithm.

We followed Sylvester and Terhune (2017) to locate the position of multiple VOIs based on the location of the sliding semilandmarks and the selected radius of the VOI (described below). Briefly, after establishing the final position of the sliding semilandmarks, a surface normal vector was calculated at each landmark location that was directed into the trabecular structure. One VOI was established for each of the 70 sliding landmarks by positioning the center of the VOI along its respective surface vector such that the most superficial surface of the VOI was coincident with the external surface of the bone. Three metrics of trabecular bone architecture – bone volume fraction, trabecular thickness and trabecular spacing – were measured from each VOI. Because each of the image stacks (e.g. segmented, trabecular thickness map) represents the exact same bone structure, the VOI positions could be used in each image stack. BV was extracted by counting the number of white voxels within each VOI; divided by the TV of the VOI that resides inside the condyle, this produced a measure of the fraction of bone in a given volume (BV/TV). Trabecular thickness (Tb.Th) and trabecular spacing (Tb.Sp) values were calculated based on thickness values determined using BoneJ and contained in the thickness and spacing map image stacks. These measures are reported in millimeters (mm). We did not quantify the degree of anisotropy in each of the VOIs, because this metric demonstrates variation upon repeatedly measuring the same VOI (Sylvester and Terhune, 2017). Additionally, we did not consider trabecular number (Tb.N) in our analyses, as previous work (e.g. Barak et al., 2013) suggests that this variable is primarily important when considering interspecific patterns in bony architecture across a wide range of taxa that vary in body size, which is not the case here.

We empirically determined the VOI position as described above for three different VOI diameters (1.0, 1.5 and 2.0 mm). We then randomly selected seven specimens and calculated all trabecular bone parameters for each VOI. We examined the relationship between VOI diameter and trabecular bone properties and determined that parameter values were similar using the 1.5 and 2.0 mm VOI, while the smaller 1.0 mm VOI produced different results. We therefore opted to use the smallest VOI that returned the most consistent results (in this case, 1.5 mm). Because all specimens are of nearly equal size, we did not apply any scaling factors for the VOI size across specimens (Lazenby et al., 2011; Kivell, 2016).

For each VOI, we extracted a measure of VOI quality, which quantifies the proportion of the VOI embedded within the bony structures. Because the condyles are small and irregularly shaped, portions of the VOI may protrude outside the bone after it has been positioned; the proportion of the VOI embedded in the bone structure was therefore calculated using the bone mask image stack. We examined the VOI quality values prior to subsequent analyses. Only 15 VOIs were more than 90% within the bone structure for all specimens (following Sylvester and Terhune, 2017), and all 15 VOIs were retained for further analyses. An additional 13 VOIs were retained for which only a few (≤2) specimens had VOI quality values of 80–90%. In order to maximize our sample of condyles and VOIs, right-side condyles for two individuals – ND31 and ND34 – were excluded during this assessment. These individuals are represented in subsequent analyses only by their left sides. As a result of this analysis of VOI quality, the sample used for analyses included 73 specimens/sides, each with 28 VOIs. All analyses described above were conducted in MATLAB R2018a (MATLAB and Statistics Toolbox 2018) and Avizo, with data exported to Microsoft Excel for analysis.

For the 28 VOIs selected for analysis, we first examined whether there was symmetry between the left and right sides of each individual by calculating correlation coefficients (r) between corresponding VOIs for each specimen. For nearly all individuals there was a significant [P≤0.00143; critical alpha (α) calculated as 0.05/35] relationship between left and right sides (BV/TV average r=0.95; Tb.Th average r=0.87; Tb.Sp average r=0.70) (Table S1). There was no significant relationship between left and right sides in Tb.Sp for four specimens, though the overwhelming pattern was still one of a significant relationship for this variable. Given these strong relationships, all subsequent analyses were performed on the average of left and right sides for each subject (N=38 specimens).

Using the side-averaged data for each variable, we ran a between-group principal components analysis (bgPCA; Yendle and MacFie, 1989; Mitteroecker and Bookstein, 2011) using treatment cohort as the group. This analysis allowed us to examine the aspects of trabecular bone variation that maximized group separation. To visualize variation in trabecular properties, we generated color maps showing the distribution of the relevant trabecular bone parameter for ±2 s.d. from the mean along each axis. We then examined the principal component (PC) scores from the bgPCA using one-way analysis of variance (ANOVA) and correlation coefficients to test for differences among the treatment groups, relationships among variables, and the relationship between trabecular bone parameters and bgPC scores. To test for differences between treatment groups, bgPC scores for BV/TV were subjected to a one-way ANOVA with Tukey honestly significant difference (HSD) tests for multiple post hoc comparisons. We predicted that, as the number of hay cubes consumed during the experimental period increased, bone volume fraction and trabecular thickness would increase (i.e. Control<Early=Late<Annual) and trabecular spacing would decrease (Control>Early=Late>Annual). We made no a priori predictions regarding differences between the Early and Late treatment groups. We also examined the correlations between the bgPC scores and average BV/TV, Tb.Th and Tb.Sp per specimen (i.e. mean of all VOIs for a specimen). Because we examined all three bgPC axes per variable, critical α for the ANOVAs and correlation analyses was set at 0.0167 (α=0.05/3). These analyses were conducted in the programs PAST and IBM SPSS (version 25; IBM Corp., 2017).

In addition, we conducted a GM analysis of condylar shape and size (measured as centroid size) to determine whether there were differences between treatment groups. Using the Procrustes shape coordinates described above (Fig. 1), we performed a bgPCA, again to examine how groups differed in shape space, but in this case with regard to condylar shape. We further calculated the mean Procrustes distance between groups to assess whether groups were significantly different in shape space; P-values of these distances were generated via permutation tests (10,000 iterations). We tested whether groups were significantly different in condylar size via one-way ANOVA with critical α set to 0.05. Finally, we performed a series of multivariate regression analyses where the Procrustes residuals were regressed on the natural log of a predictor variable (average BV/TV, Tb.Th, Tb.Sp or centroid size). These analyses were conducted in the program MorphoJ (Klingenberg, 2011).

Average values for all trabecular parameters are presented in Table 2 and variation across VOIs illustrated in Fig. S2 and provided in Table S2. In general, average BV/TV and Tb.Th followed the trend Control<Early<Annual<Late, while Tb.Sp was roughly opposite (Late<Annual=Early<Control). However, not all differences were statistically significant.

Color maps representing the average BV/TV distribution of all 38 rabbits, and averages for each treatment group, are in Fig. 2. In the overall and group averages, BV/TV was highest on the anteromedial portion of the condyle and followed a gradient in which the lowest BV/TV values were found towards the posterolateral portion of the condyle.

The bgPCA provides details about differences in BV/TV values among treatment groups, which can be appreciated in the group-specific color maps (Fig. 2). All groups overlapped to some extent in the bgPC plots (Fig. 3, top). Average BV/TV increased significantly along bgPC1 (r=0.994, P<0.0001; Fig. 3, top; Table 3). This axis was also positively correlated with average Tb.Th (r=0.728, P<0.0001) and negatively correlated with average Tb.Sp (r=−0.636, P<0.0001) (Table 3). ANOVA results (Table 4) indicated significant differences between Control versus Late and Annual groups on bgPC1 (BV/TV was lower in the Control group), and between Early and Late groups (BV/TV was lower in the Early group). bgPC2 (Fig. 3) describes differences in the pattern of BV/TV across the articular surface. Positive scores on bgPC2 indicate BV/TV values that are more similar across the articular surface (i.e. lower variance among BV/TV values for the VOI), while negatively loading specimens have higher BV/TV values along the anterior aspect of the condyle and lower values on the posterior aspect of the condyle (i.e. a larger range of values; Table 2). On bgPC2, there were significant differences between the Annual group and all other groups, indicating that Annual rabbits have greater variance in BV/TV values across the articular surface, and exhibit relatively high BV/TV values on the anterior condyle and relatively low values on the posterior condyle (Figs 2 and 3, Table 4). bgPC3 (not shown; 3.3% of sample variance) successfully differentiated between Early and Control groups (Table 4). However, this axis represents only minor differences in the distribution of BV/TV across the articular surface and is therefore challenging to interpret functionally.

In general, trabecular thickness followed a pattern similar to BV/TV in which thicker trabeculae were located on the anteromedial aspect of the condyle and trabeculae became thinner posteriorly and laterally (Fig. 2). This similarity is unsurprising given the significant positive correlation between average trabecular thickness and average BV/TV values (r=0.730, P<0.0001) (Table 3). Unlike BV/TV, however, the most posterolateral VOIs had slightly higher Tb.Th values than those in the central portion of the condyle (Fig. 2).

As is apparent in the bivariate plot of bgPC1 and bgPC2 (Fig. 3, middle) and reflected in the ANOVA results (Table 4), treatment groups did not separate along bgPC1. This axis was significantly positively correlated with average Tb.Th (r=0.996, P<0.0001) and to a lesser extent average BV/TV (r=0.758, P<0.0001) (Table 3). Thus, bgPC1 primarily reflects the magnitude of Tb.Th increase, which was roughly evenly distributed across the articular surface (Fig. 3) and varied considerably within groups. bgPC2 was more effective at differentiating between treatment groups and was significantly correlated with average Tb.Th (r=0.411, P=0.01) (Table 3). Only bgPC2 showed a significant ANOVA P-value after Bonferroni correction (Table 4), and indicated significant differences between Annual versus Control and Late groups, as well as between Early versus Control and Late groups (e.g. Annual and Early groups had more positive scores on this axis, which itself was positively correlated with average Tb.Th). Importantly, bgPC2 again appears to represent differences in the patterning of Tb.Th across the articular surface. Specimens loading more negatively on this axis had higher Tb.Th values focused on the VOI in the center and anterior aspects of the articular surface, whereas specimens loading more positively had higher trabecular thickness values on the edges of the articular surface (though these specimens also had overall lower Tb.Th values). There were no significant differences among groups in average Tb.Th values (Table 4).

Color maps showing the pattern of Tb.Sp across the condyle are provided in Fig. 2 for the total sample average and for group averages. In general, spacing was smaller on the anteromedial portion of the condyle and higher on the posterolateral portion of the condyle, as expected based on the patterns of BV/TV and Tb.Th.

The bgPCA of the Tb.Sp data (Fig. 3, bottom) indicated that groups were largely indistinguishable in Tb.Sp, though some outliers in the data are labeled in this plot. bgPC1 was positively correlated to average Tb.Sp (r=0.996, P<0.0001) (Table 3), indicating that spacing increased as values on this axis became more positive (Fig. 3). bgPC1 was also negatively correlated with average BV/TV (r=−0.670, P<0.0001) (Table 3), indicating an inverse relationship between spacing and BV/TV.

The GM analysis revealed no differences in condylar shape or size between groups. All Procrustes distances between groups returned non-significant P-values (P>0.05), and an ANOVA examining differences in centroid size between groups was not significant (P=0.304). Further, there were no significant relationships between centroid size and any measure of average BV/TV, Tb.Th or Tb.Sp, or between most bgPC scores for any of the trabecular parameters (Table 3). However, one interesting result was that bgPC2 from the Tb.Sp analysis was significantly negatively correlated with centroid size (r=−0.449, P=0.005). This appears to be a result, however, of one particular outlying specimen (ND31). When this specimen was removed, the correlation was lost (Table S3). Multivariate regressions of the Procrustes residuals of the condylar shape coordinates on centroid size, average BV/TV, Tb.Th and Tb.Sp found no significant relationship between condylar shape and centroid size or trabecular spacing, but there were significant relationships between condylar shape and average BV/TV (% predicted=7.69, P=0.02) and average Tb.Th (% predicted=8.2, P=0.016). Examination of shape variation related to these regressions indicated that for both BV/TV and Tb.Th, individuals with higher values tended to have slightly anteroposteriorly longer and mediolaterally narrower condyles. Notably, these patterns cut across groups.

One rabbit, ND31, was identified as a statistical outlier. Though no behavioral differences between this rabbit and the other Annual rabbits were observed during the experiment, this specimen fell on the far edges of the Annual distribution in the bgPC plots and had the second lowest average BV/TV of the entire sample (Table 2). This specimen also had particularly large mandibular condyles: its centroid size was 3.4 standard deviations above the mean. However, analyses where this specimen was excluded revealed a similar pattern of results to those above (Tables S3 and S4). One important difference was that BV/TV values (both average values and bgPC scores) were statistically different (Table S4) for Early versus Annual groups. Further, the bgPC1 scores for Tb.Th were statistically significantly different between Control versus Late and Annual groups.

The study presented here examined the link between loading environment and the bony architecture of the mandibular condyle. Using a novel experimental sample of rabbits raised from weaning well into mature adulthood on diets with different mechanical properties, we demonstrate that, while some architectural parameters differed significantly among experimental groups, treatment groups also exhibited a shared pattern of trabecular parameters across the articular surface of the mandibular condyle. We hypothesized that BV/TV and trabecular thickness would be highest in the Annual group compared with all other groups and that these values would be lowest in the Control group. Additionally, we anticipated this pattern would be reversed for trabecular spacing. We made no a priori predictions regarding differences between Early and Late groups. Observed differences among treatment groups were generally consistent with our predictions, though BV/TV showed the clearest signal among treatment groups. Indeed, depending on the variable examined, all groups differed significantly from one another, though the overarching pattern was one where Control and Early groups were most similar to one another, as were Annual and Late groups. These findings suggest that trabecular parameters (especially BV/TV) vary in relation to mechanical loading during development and that the timing of the introduction of mechanically challenging foods has an important impact on trabecular structure.

Our data indicate that some experimental groups of rabbits differed in the magnitude of BV/TV, while one group (the Annual rabbits) was distinguished by a subtle difference in the distribution of BV/TV across the articular surface. In general, rabbits that consumed a more mechanically challenging diet (i.e. Annual, Late) had elevated BV/TV values relative to Control rabbits. Average BV/TV per group demonstrated a clinal distribution (i.e. Control<Early<Annual<Late), but adjacent groups were not always significantly different in their average BV/TV values. Notably, the Late treatment group had the highest average BV/TV values, exceeding even those of the Annuals, which tended to be more variable. The results of the bgPC analysis further reinforced this clinal distribution, and across all three bgPC axes significant differences were present among groups in magnitude and/or pattern of BV/TV distribution.

The adaptive response of bone has been shown to be a function of several aspects of mechanical loading, including strain magnitude, strain rate, strain frequency, strain distribution, number of loading cycles and rest–recovery periods (Biewener, 1993; Lanyon, 1996; Turner, 1998; Hart et al., 2017). Hay is stiffer (3335.6 MPa) and has higher fracture toughness (2759.8 J m⁻²) than pellets (29.2 MPa, 1030.6 J m⁻²), and should therefore require higher magnitude bite forces and/or a greater number of chewing cycles to process. While muscle forces are challenging to estimate and quantify, Ravosa et al. (2015) documented an approximately threefold increase in chewing duration and investment (chews per gram) for hay relative to pellets in rabbits. Thus, even in the absence of higher bite forces (cf. Weijs and de Jongh, 1977), the increased number of loading cycles in Early, Late and Annual groups is expected to stimulate bone deposition (Bouvier and Hylander, 1981; Ozcivici et al., 2010; Ravosa et al., 2016). Indeed, the general pattern of increased BV/TV we observed in Annual and Late rabbits in comparison to Control, and to a lesser extent Early, individuals accords with the current understanding of bone functional adaptation.

Interestingly, the Late group was distinguished from the Early group in terms of BV/TV magnitude. These two experimental cohorts only varied in the timing of receiving hay cubes. The Early group was fed hay at the beginning and middle of the experiment – just after weaning and around the time of skeletal maturity, respectively. In contrast, Late rabbits were first provisioned with hay cubes during the middle and end of the experiment. The high BV/TV observed in the Late group indicates that bone adaptation is possible after attainment of skeletal maturity, which contrasts with prior evidence that plasticity in vertebrates decreases with age (Hinton and McNamara, 1984a; Meyer, 1987; Bouvier, 1988; Bouvier and Hylander, 1996a,b; Rubin et al., 1992). The low levels of BV/TV in the Early group – similar to those of Control rabbits – suggests that bone volume is lost when the loading environment is altered. Further, because the Late rabbits were provisioned with hay for the final 6 weeks of the experiment, their higher levels of BV/TV suggest that active addition of bone in the condyle was taking place at this time. Although we do not have ontogenetic data for the trabecular variables analyzed here, these interpretations of BV/TV in Late and Early rabbits are supported by previous analyses. Scott et al. (2014a,b) found that the Early rabbits diverged considerably from both Control and Late rabbits in corpus, symphyseal, palatal and condylar shape (bone cross-sectional area size-adjusted by cranial length) early in the experiment when the Early rabbits were receiving hay cubes. However, once the Early rabbits were switched to an all-pellet diet, they quickly developed morphologies similar to those observed in the Controls. We consider it likely that the Early rabbits would also have exhibited increased BV/TV relative to the Control and Late rabbits during this period, and that this relative increase was lost after the Early group transitioned to an all-pellet diet. These data therefore indicate that the consumption of mechanically challenging hay cubes results in increased BV/TV, but the osteogenic effect disappears after ceasing hay provisions and mechanical loads likely return to those experienced by Control rabbits.

One factor contributing to the ability of the condylar trabecular structure to distribute loads effectively and avoid failure during increased loading is bone mineralization or density. Although we did not examine bone density, prior work with this same experimental group by Franks et al. (2017) involving the Control and Annual (‘over-use’) cohorts documented a mosaic pattern of bone mineralization across the masticatory apparatus (i.e. corpus, symphysis and hard palate) and neurocranium. Specifically, Annual/over-use rabbits exhibited decreased mineralization in several masticatory sites versus Control rabbits except for the hard palate, where levels of mineralization were significantly higher than in Controls. This may indicate that structures experiencing higher loads (in this case, the mandibular corpus) exhibit greater bone turnover (i.e. increased rates of remodeling) and therefore are less mineralized because they are less mature (e.g. Meunier and Boivin, 1997; Cullen et al., 2001; Boivin et al., 2009). This could be particularly relevant for comparison of Annual and Late groups here; even though the differences in BV/TV between the Annual and Late groups did not reach statistical significance, the Late group had unexpectedly higher and less variable BV/TV values (Table 2) than the Annual rabbits. This could be a response to relatively short-term loading of the joint in the Late rabbits and perhaps also result from differences in bone density that could factor into how loads are distributed throughout the trabecular structure. One possibility is that, similar to processes identified for osteoarthritis (Burr and Gallant, 2012), the initial bone response to increased loading may be one of high bone volume but low bone mineralization. Ongoing remodeling in response to continued loading could then be focused on increasing mineralization with a reduction in bone volume (i.e. BV/TV) (e.g. Cullen et al., 2001).

In addition to differences among the groups in the magnitude of BV/TV, the Annual group was distinguished from the other groups in terms of a subtle pattern difference in how BV/TV is distributed across the articular surface. While the BV/TV values for the Annual rabbits on the anterior portion of the condyle were closer to values observed in the Late group, values on the central portion of the condyles were closer to those observed in the Early group. In other words, while the Annual group had high overall BV/TV values across the articular surface, bone seems to be disproportionately added along the anterior portion of the condyle, perhaps reflecting increased loading of this portion of the joint in the Annual rabbits. Interestingly, we did not observe a difference in the pattern of how BV/TV was distributed across the articular surface in the Late rabbits despite their overall high BV/TV values. This could potentially be due to habitual differences in how the Annual group loaded their condyle relative to the Late group, as well as the other groups. Further, the disparity between the Annual group and other rabbit cohorts might represent a threshold-based osteogenic response due to the protracted (versus ‘seasonal’) consumption of hay cubes. Moreover, because rabbit condyles are absolutely small, increased loading and addition of bone in one portion of the joint (i.e. the anterior aspect) may require the concomitant addition of excess bone in adjacent lower-strain areas to ensure the overall structural integrity of the condyle. Such a finding about the differentially greater role of bone adaptation is consistent with load-related decreases in the stiffness of TMJ articular cartilage in Annual rabbits (Ravosa and Kane, 2017). Viewed in a comparative and paleontological context, this pattern may contribute to the variable correspondence between diet and jaw form, particularly as some sites and parameters may exhibit pronounced and/or prolonged plasticity responses.

Differences between groups in trabecular thickness and spacing mostly did not reach statistical significance, though this result varied slightly when the outlier ND31 was removed from analysis. This general lack of statistically significant differences among groups is surprising given the differences in BV/TV among treatment groups, because BV/TV must, at some level, be a function of trabecular thickness and spacing. Regression analysis of the bgPC scores indicated that this was indeed the case here. Multiple regression of BV/TV bgPC1 scores on both Tb.Th bgPC1 and Tb.Sp bgPC1 scores returned an R² value of 0.89. Similar results (i.e. high R²) were obtained when we examined relationships between average BV/TV, Tb.Th and Tb.Sp, and for these variables when examined separately per VOI. Thus, variation in Tb.Th and Tb.Sp explained most of the observed variation in BV/TV across individuals in the sample, even in the absence of other measures of bony architecture (e.g. Tb.N). When assessed separately, Tb.Th bgPC scores explained ∼58% of the variation in BV/TV values (R²=0.58), while Tb.Sp bgPC1 scores explained roughly 43% of BV/TV variation (R²=0.43). The regression coefficients for Tb.Th were positive, while those for Tb.Sp were negative, demonstrating a direct relationship between BV/TV and Tb.Th, and an inverse relationship between BV/TV and Tb.Sp. Notably, there was no correlation between Tb.Th and Tb.Sp values.

One possible reason we did not observe consistent differences in trabecular thickness and spacing among treatment groups is that groups and/or individuals achieve variation in BV/TV values in different ways. In other words, individual rabbits may achieve high BV/TV values either by increasing trabecular thickness or by decreasing spacing, or some combination thereof. Thus, because BV/TV carries information about trabecular thickness and spacing simultaneously, this variable has more power to distinguish groups with different loading patterns. This multifarious potential for how individuals approach the same biomechanical solution even at the trabecular level warrants further examination and may represent another example of a many-to-one structure–function relationship (e.g. Wainwright et al., 2005).

Several rabbits serve as points of discussion. ND7 died during week 33 of the experiment but after reaching skeletal maturity. This individual had the lowest average BV/TV of the Control group but not the lowest BV/TV of all rabbits. In some of the bgPC plots, it fell on the margins of the convex hull for the Control rabbits but did not differ markedly from others of that group. These results suggest that, although this individual may have been skeletally mature (i.e. external skeletal dimensions were indistinguishable from adults), the low BV/TV values could reflect either a process leading to its early mortality or a slightly earlier developmental stage of trabecular bone development. As indicated above, one rabbit from the Annual group (ND31) was a strong outlier in our analyses, and had one of the lowest BV/TV and Tb.Th values for the entire sample. Though this individual appeared to be behaviorally consistent with the rest of the Annual rabbits, its mandibular condyles were particularly large in comparison to those of all other rabbits, which may explain its outlier status: condylar size may influence loading patterns in the TMJ, whether because of pathology, changes in stress and strain distributions in the joint, or some other factor.

The general patterns of trabecular bone parameters (BV/TV, Tb.Th, Tb.Sp) across the rabbit condyle demonstrated a non-uniform distribution of bone, but one that was consistent across experimental groups. This pattern was characterized by high BV/TV, thick trabeculae and narrow trabecular spacing on the anteromedial aspect of the mandibular condyles, with a general trend of decreasing BV/TV and thickness and increasing spacing laterally and posteriorly. These observations suggest that the anteromedial portion of the condyle is subjected to the highest loading during chewing cycles, while the posterior and lateral portions are loaded less heavily. In a review of the rabbit TMJ, King et al. (2010) referred to the mediolaterally wider anterior portion of the condyle as ‘articular surface’, and the narrow posterior portion of the condyle as ‘non-articular’. Because rabbits have reduced gape abilities (Weijs and Dantuma, 1981), it is likely that the more posterior joint surface may be used infrequently for increased gapes required for uncommonly large items. If true, examination of trabecular variation across the entire structure of the condyle, including the posterior aspect, may yield important clues to dietary adaptations, such as those involving the generation of muscle and bite forces at relatively large gapes.

As part of this analysis, we further examined variation in condylar size and shape across treatment groups. Like previous work by Scott et al. (2014a,b) that examined condylar cross-sectional area, there were no differences in condylar centroid size among groups. Further, we observed no significant differences in condylar shape among groups, though we did find some weak relationships between condylar shape and BV/TV and average Tb.Th. These findings suggest an association between mediolateral width and anteroposterior length of the condyle and trabecular parameters, although only a small amount of variation in shape is explained by these measures of trabecular structure. While these findings warrant further consideration, they are roughly consistent with previous work suggesting that because of requirements for maintaining function, joints may be highly canalized and/or developmentally stable, buffered from developmental perturbations that may accompany changes in the mechanical loading environment (Lieberman et al., 2001; Auerbach and Ruff, 2006; Reeves et al., 2016). Notably, however, other work has identified considerable plasticity in the dimensions of the mandibular condyle in relation to different dietary demands (e.g. Bouvier and Hylander, 1984; Bouvier, 1988; Ravosa et al., 2007). This is contrary to arguments that joints are buffered to environmental perturbations and requires further investigation, particularly with regard to the need for additional work on long-term plasticity in multiple skeletal joints across diverse mammals.

One major reason that our results differ from prior analyses in showing a clear link between trabecular parameters and loading environment is our use of multiple VOIs spread across the articular surface of the condyle. Multi-VOI approaches provide information about variation in trabecular bone (or, in this study, trabecular bone and a small portion of cortical bone along the articular surface) across the joint surface not accessible via single VOI approaches (e.g. Ryan et al., 2010), or analyses that treat all the trabecular bone within the mandibular condyle as a single structure (e.g. Coiner-Collier et al., 2018). Further, variation in BV/TV across the condyles highlights the methodological challenges of positioning a single VOI to investigate group differences. For example, in examining the BV/TV data for the total sample average, VOI23 and VOI24 (adjacent VOI on the anteromedial portion of the condyle, <1 mm apart; Fig. 1) differ in value by 14.9 percentage points (BV/TV: 64.5% versus 49.6%). By way of comparison, the maximum absolute difference between experimental group averages (across homologous VOI) is 16.8 percentage points (Control versus Late, VOI14), although the average of all such comparisons is 12.7 percentage points (range 9.3%–16.8%). Thus, variation within a small region of the mandibular condyle is of the same magnitude as, and often exceeds, variation between the experimental groups investigated here.

Given that the rabbits were fed diets that varied significantly in material properties, and yet variation in BV/TV within a condyle was similar to, or exceeded, variation between experimental groups, the position of the VOI is critically important. Without precise positioning, the signal that distinguished groups could have easily been lost in the variation of the trabecular structure. Locations of VOI become especially important when investigating potentially subtle differences among species that process food items that do not vary drastically in material properties. In addition, it is important to note that the center of the condyle (i.e. middle of a line connecting the most anterior and posterior projecting points on condyle) lies in the area of low BV/TV that King et al. (2010) characterized as non-articular. This further argues against constructing single VOI positions geometrically, as they may miss functionally relevant areas of the articular surface. This consideration could be especially important as overall joint size, and corresponding functional complexity of the joint, increases. Collectively, these results suggest that architectural properties of the mandibular condyle warrant reanalysis and may convey unappreciated clues to masticatory function.

Bone functional adaptation is often invoked as an interpretative paradigm for understanding skeletal variation in the context of interspecific behavioral diversity (e.g. Fajardo and Müller, 2001; Griffin et al., 2010; Morimoto et al., 2011; Ryan and Shaw, 2012, 2013; Barak et al., 2013; Chirchir et al., 2015). Despite its widespread use and potential application to the masticatory apparatus and dietary research, relatively little work has been conducted to establish the relationship between dietary adaptations, mechanical loading and the response of the masticatory apparatus to induced stresses. Prior work has demonstrated that variation in jaw musculature (e.g. Taylor et al., 2006) and mandibular cortical bone (e.g. Hylander, 1979; Bouvier, 1986; Daegling, 1989; Holmes and Ruff, 2011; Ravosa and Kane, 2017) are linked to variation in the mechanical challenges of food processing across taxa. However, only a handful of comparative analyses have examined variation in condylar trabecular structure between species, and the findings are largely equivocal. Ryan et al. (2010) examined trabecular structure in platyrrhine monkeys that gouge trees with their anterior teeth versus those that do not and found no significant differences among taxa for most trabecular variables. They suggested that their results indicate either that there are no substantive differences in TMJ load in gouging versus non-gouging behaviors or that trabecular architecture is not mechanically important and is perhaps ‘functionally uninformative’ (Ryan et al., 2010, p. 583). Coiner-Collier et al. (2018) identified a relationship between food toughness and trabecular anisotropy across a sample of 11 species of anthropoid primates. However, they did not find a link between food toughness and BV/TV or most other measures of trabecular structure. Thus, a clear relationship between condylar trabecular structure and loading environment across species has yet to be established, due in part to a limited number of studies. One possible reason for this disjunction between trabecular structure and masticatory function is that animals modulate loading of the TMJ behaviorally, obscuring a link with diet (sensu Ryan et al., 2010). Alternatively, it is possible that dietary stiffness plays a greater role in the development of certain masticatory parameters (Ravosa et al., 2015, 2016).

In addition to the methodological differences discussed above, one reason why our data reveal patterns of bone functional adaptation while others do not could be linked to the high-resolution nature of our dataset. Experimental studies such as the one presented here have more power to control and quantify factors that could otherwise influence bone morphology (e.g. diet, genetics, environment), which are not easily controlled or examined in analyses of wild species. In contrast, comparative studies of trabecular structure are frequently plagued by a host of potential confounding factors – some known, some unknown, some unknowable – that may influence analytical outcomes. For example, broad interspecific studies (e.g. Barak et al., 2013; Ryan and Shaw, 2013; Coiner-Collier et al., 2018) often focus on only limited specimens per species, and it is difficult to know whether those individuals experienced stereotypical loading regimes during their lifetime. Geographic and temporal variation within species may also play a critical role, given well-documented examples of how conspecific populations exploit foods with different material properties (e.g. Palombit, 1997; Chapman et al., 2004; Kamilar, 2006). Our data suggest that the season/timing of death is likely a critical determinant of trabecular structure, particularly for animals that experience large seasonal fluctuations in food material properties and therefore loading regimes. Finally, ontogenetic variation in the pattern and duration of loading events may further confound functional signals, especially as prior research suggests that trabecular structure, and, more generally, plastic responses to changes in loading environment, vary substantially in relation to age both before and after the attainment of skeletal maturity (Hinton and McNamara, 1984b; Parfitt et al., 2000; Lieberman et al., 2003; Pearson and Lieberman, 2004; Ravosa et al., 2008; Raichlen et al., 2015). Though it is unlikely that any single comparative study of wild species can control for all of these factors, we suggest that considerable scrutiny should be given to the size and composition of samples so as to minimize these sources of variation. Long-term experiments, such as those reported here, can advance our understanding of these factors, although further work is necessary to understand fully how dietary properties influence feeding behavior and jaw-loading patterns. Such information, including a hierarchical perspective on bone adaptation and a multi-tissue perspective for joints, is critical for the determination of form–function links throughout the skeleton (Ravosa et al., 2007, 2016; Ravosa and Kane, 2017).

Another complicating aspect of analyses examining adaptive plasticity and bone functional adaptation is that different skeletal regions exhibit variation in their sensitivity and responses to mechanical loading (Rawlinson et al., 1995; Pitsillides et al., 1999; Hsieh et al., 2001; Judex et al., 2004; Hamrick et al., 2006; Ravosa et al., 2010b, 2016; Franks et al., 2017). While a considerable body of work exists documenting variation in trabecular structure in relation to inferred patterns of loading of the postcranial skeleton (e.g. Ryan and Ketcham, 2002; Griffin et al., 2010; Ryan and Walker, 2010; Kivell et al., 2011; Chirchir et al., 2015), fewer studies have examined the masticatory system. As the mandible routinely experiences significant mechanical strain during masticatory and paramasticatory behaviors (Hylander et al., 1987), the feeding apparatus should be equally subject to bone functional adaptation as the postcranium (Ravosa et al., 2007, 2008, 2010a,b, 2016; Menegaz et al., 2009; Scott et al., 2014a,b; Ravosa and Kane, 2017). Our findings are consistent with patterns of cortical bone adaptation noted for the same experimental sample (Scott et al., 2014a,b), as well as other studies of adaptive plasticity in the jaws (Bouvier and Hylander, 1981; Beecher et al., 1983; Bouvier, 1988; Yamada and Kimmel, 1991; He and Kiliaridis, 2003; Lieberman et al., 2004; Ravosa et al., 2007, 2008, 2010a; Chen et al., 2009; Menegaz et al., 2009; Rafferty et al., 2012; Balanta-Melo et al., 2019b). The diet-induced responses in TMJ trabecular architecture observed among growing rabbits suggests that mammals are able to ‘fine-tune’ the internal structure of the TMJ to respond to changes in mechanical loading during their lifetime. This flexibility is perhaps best exemplified by our comparison of the Early and Late cohorts, where the timing of the onset of a mechanically challenging diet varied throughout the length of the experiment.

When coupled with prior work examining other aspects of the masticatory apparatus in this same experimental cohort, it is clear that different parts of the skull, and indeed different hierarchical levels of bone, respond to the same dietary treatment in different ways and to different extents (Scott et al., 2014a,b; Ravosa et al., 2016; Franks et al., 2016, 2017). For example, while increased cross-sectional area and cortical bone thickness was observed in the corpus, symphysis and palate in groups raised on more mechanically challenging diets, no significant differences were observed in zygomatic arch cortical thickness (Franks et al., 2016, 2017). Further, although biomineralization at the symphysis and corpus decreased in groups with more mechanically challenging diets, the opposite pattern was observed for the zygomatic arch and hard palate (Franks et al., 2016, 2017). Other experimental work has found similar site-specific variability in bone adaptation across the skull (reviewed by Ravosa et al., 2016). This mosaic plasticity response to loading suggests that analyses focusing on only one skeletal site and level of hard-tissue organization are more likely to mischaracterize the performance of a bony element, particularly if it is part of a multi-tissue joint (Ravosa et al., 2007; Ravosa and Kane, 2017). This especially complicates interpretations of fossil morphology, where behavioral inferences are typically based on a single structure (and hierarchical level) from relatively small samples.

In sum, this study demonstrates that changes in mechanical loading stimulate changes in the trabecular architecture of the TMJ, with some plasticity responses evident relatively late in ontogeny. Combined with previous work on this experimental sample, these data suggest that the masticatory apparatus responds to changes in mechanical loading at a variety of levels. Importantly, our data also suggest that the relationship between loading and bone morphology is complex and likely mediated by a variety of factors, including the developmental timing of the onset of particular loading regimes, the duration of the regimes and the region of the masticatory apparatus under investigation. Additional analyses are warranted to evaluate whether bony responses of trabecular bone to changes in loading track dietary adaptations across mammal species, and the extent to which diet-related plasticity contributes to ongoing adaptation among species.

We thank Christopher Carter, Kevin McAbee, Susan Coiner-Collier, Connor Shea and the staff of Friemann Life Science Center for help with animal handling, CT scan generation and processing.