Geometric morphometric investigation of craniofacial morphological change in domesticated silver foxes

To test the effects of domestication on craniofacial skeletal morphology, we used three-dimensional geometric morphometrics (GM) along with linear and endocranial measurements to compare selected (domesticated) and unselected foxes from the Russian Farm-Fox Experiment to wild foxes from the progenitor population from which the farmed foxes are derived. Contrary to previous findings, we find that domesticated and unselected foxes show minimal differences in craniofacial shape and size compared to the more substantial differences between the wild foxes and both populations of farmed foxes. GM analyses and linear measurements demonstrate that wild foxes differ from farmed foxes largely in terms of less cranial base flexion, relatively expanded cranial vaults, and increased endocranial volumes. These results challenge the assumption that the unselected population of foxes kept as part of the Russian Farm-Fox experiment are an appropriate proxy for ‘wild’ foxes in terms of craniofacial morphology and highlight the need to include wild populations in further studies of domestication syndrome to disentangle the phenotypic effects of multiple selection pressures. These findings also suggest that marked increases in docility cannot be reliably diagnosed from shape differences in craniofacial skeletal morphology.


Estimation strategy
We analyzed linear and volumetric measurement data using linear regression with a generalized least squares (gls) estimator. The gls estimator of the coefficients of a linear regression is a generalization of the ordinary least squares (ols) estimator used in the general linear model (glm). It is used to deal with situations in which the ols estimator is not the best linear unbiased estimator (blue) because one of the main assumptions of the Gauss-Markov theorem, namely that of homoskedasticity and absence of serial correlation, is violated. In such situations, provided that the other assumptions of the Gauss-Markov theorem are satisfied, the gls estimator is the blue. The general linear model is: where #" y is an n × 1 vector of responses (n is the sample size), X is a n × k design matrix of regressors (k is the number of regressors), #" β is a k × 1 vector of weights (regression coefficients) to be estimated, and #" is an n × 1 vector of errors. We assume that: 1. X has full rank, 2. E[ #" y | X] = X #" β , 3. E[ #" | X] = 0, 4. Var[ #" | X] = σ 2 I, where σ is a constant and I is the n × n identity matrix.
These are the assumptions made in the Gauss-Markov theorem in order to prove that ols is the blue. 1 Assumption 4 means that the errors of the regression are homoskedastic (they all have the same variance) and uncorrelated (their covariances are all equal to zero). When we use the gls estimator, we generalize assumption 4 so that: where V is an n × n symmetric positive definite matrix. 1 This allows for heteroskedasticity (the errors can have different variances) and correlation (the covariances between errors can be different from zero): So, the gls estimator makes fewer assumptions than the ols estimator. If the off-diagonal elements of V are restricted to be zero (i.e., covariances between errors are not allowed), then we have the weighted least squares estimator (wls).
In this study, for each specimen, 6 different variables were recorded using linear measurements and 1 variable using volumetric measurement. Our primary interest was to regress each of these variables on population identity and sex. However, analyzing these data by estimating a series of separate glms (e.g., anova, ancova, linear regression etc.) would result in biologically unrealistic inferences. This is because the 7 variables are correlated (see Fig. S3) for two related reasons: 1) because they are measured on the same specimens, and 2) because they represent non-independent aspects of shape variation. Since a series of 7 glms would involve fitting 7 separate, independent equations, there would be no mechanism to account for the correlations among the responses. Using the glm approach, the response correlations would be fixed at zero, rather than estimated within the model. The gls estimator, on the other hand, extends the functionality of the glm so that the correlations between responses can be estimated and accounted for within a single model. Failing to account for these correlations between responses would result in inferences that are overly optimistic (i.e., interval estimates that are too narrow) and therefore biologically unrealistic.

Model specification
To include multiple responses in a gls model framework, the values for each response variable must be 'stacked' on top of one another to form a single response vector #" y . The basic explanatory variables are a set of dummy variables (i.e., with values of 0 or 1) that indicate which response variable is present. Further explanatory variables can be included by multiplying these dummy variables (i.e., creating interactions), so that separate effects can be estimated for each response.
In the design matrix X, categorical variables with g levels were converted into g − 1 dummy variables. The reference level for the j = 7 − 1 'response' dummy variables (cranial vault width, upper jaw width, total skull length, snout length, cranial vault height, endocranial volume) was set to 'zygomatic width'. The reference level for the k = 3 − 1 'population' dummy variables (wild, unselected) was set to 'selected', while for the single 'sex' dummy variable (male) it was set to 'female'. The choice of the reference level for all dummy variables is arbitrary and does not change the fitted model, only its parameterization.
Our gls model specification therefore included 73 specimens and 7 response variables for a total of n = 73 × 7 = 511 observations, with 28 fixed effects parameters, 21 covariance parameters, and 7 variance (error) parameters (the variance-covariance matrix V was stratified by specimen, so that only observations for a given specimen could covary). For economy of notation, the j = 7 individual response variables are referred to here as response with a j subscript.
where #" y ij is a vector of i = 73 specimens for each of j = 7 responses, β 1 is the expected difference between wild and selected foxes for response j = 1 (zygomatic width), β 2 is the expected difference between unselected and selected foxes for response j = 1, β 3 is the expected difference between male and female foxes for response j = 1, #" β 4−9 is a vector of effects for selected foxes, #" β 10−15 is a vector of effects for selected versus wild foxes, #" β 16−21 is a vector of effects for selected versus unselected foxes, #" β 22−27 is a vector of effects for female versus male foxes, and #" is a vector of errors for each combination of specimen and response. We can visualize this model specification using a path diagram (see Fig. S4). In addition, when testing hypotheses related to sexual dimorphism we extended the above specification to include the 3-way interaction (and all subordinate terms) between response, male, wild, and unselected.
For models using a gls or wls estimator, we used the Bayesian information criteron (bic) to select models with an appropriate variance-covariance structure for V. The bic is a measure of goodness-of-fit that penalizes model complexity and is typically used for explanatory models to provide a balance between parsimony and over-fitting. 2 Model assumptions were checked using diagnostic plots of standardized residuals versus fitted values and sample quantiles versus Gaussian quantiles (see Fig. S5).

Quantities of interest
We extracted quantities of interest, such as predicted marginal means and marginal effects, from the fitted model using linear combinations of coefficients in the R package emmeans v. 1.5.1. 3 The estimands of particular interest to us were the predicted marginal means of population identity and sex for each response variable and the pairwise marginal effects between these means. Point and interval (95% confidence) estimates of each size normalized effect size are reported.  Maine ?

MCZ52832
Maine 1948 Table S2: List of landmark definitions used in 3D geometric morphometric analyses. * denotes that this landmark was measured at the midline of the skull. Type 1 landmarks are supported by histology, while type 2 are supported by geometry. Type 3 landmarks are those with at least one deficient coordinate. The most posterior point of the palatine foramen 2 Table S3: Measurement repeatability analyses on individual craniofacial variables. Intra-specimen standard deviation (sd) represents the average sd of 15 separate measurements across three different domesticated female fox skulls. Inter-specimen sd represents the sd across the domesticated females in the entire study sample. The sensitivity ratio is calculated as follows: Inter-specimen SD − Intra-specimen SD /Intra-specimen SD. The ratio gives a rough estimate of the magnitude of difference between the intra-specimen sd and inter-specimen sd, with larger numbers indicating a more precise measurement.

Variable
Intra-specimen SD Inter-specimen SD Sensitivity Ratio   Table S5: Pairwise comparisons of shape variance (dispersion of residuals) around mean population shape. D: domesticated, U: unselected, W: wild. d: Effect size expressed as differences between population variances (first population minus second population). ucl (95%): 95% upper bound confidence limit for the null hypothesis of no difference in shape variance between populations. Z: Z-score of the variance of the second population with respect to the variance of the first population. P-values have been adjusted for family-wise error using the sequential Bonferroni method.