Construction of P matrices (variance-covariance matrices) and statistical analyses
We used the CalculateMatrix function from the EvolQG 0.2-9 package in R (Melo et al., 2015) to construct a variance-covariance P matrix for the following groups: 1) SSI generalist-only, 2) SSI radiation, 3) Caribbean, 4) Crescent Pond F2 hybrids, and 5) Little Lake F2 hybrids. These matrices describe the variance both within a trait and the covariation between traits for a given group and are the primary unit of comparison in following analyses. We did not use any phylogenetic correction because within-lake populations of each species are sometimes more closely related than species across lakes (Martin and Feinstein 2014). To visualize similarities and differences between P matrices of 1) Caribbean, SSI generalist-only, and SSI radiating pupfish groups and 2) F2 hybrids from Crescent Pond and Little Lake we performed two separate principal component analyses. These PC analyses were covariance-based and included the variance and covariance estimates for the 18 craniofacial traits (prcomp function; R Core Team, 2021; Figures 4 & 8). For each of the PC analyses we also calculated the correlation between loadings and PC axes 1 and 2 to determine 1) which loadings mostly closely aligned with the variation along a given axis and 2) which loadings were most similar to one another. We also calculated the contributions of each group (i.e., Caribbean, generalist-only, radiation, Crescent Pond F2s, or Little Lake F2s) and trait towards the patterns observed in the analyses using the get_pca_var and get_pca_ind functions from the factoextra package (Kassambara and Mundt 2020)
We used the integration.Vrel and compare.ZVrel functions (geomorph 4.0.6; Baken et al. 2021; Adams D, Collyer M, Kaliontzopoulou A 2023) to quantify and compare the strength of integration between matrices (i.e., relative eigenvalue variance or Vrel ) and the MeanMatrixStatistics function (EvolQG 0.2-9; Melo et al., 2015) to estimate: autonomy, constraints, flexibility, mean squared correlation, and respondability for each matrix (detailed definitions of each measurement can be found in Table 2). We further estimated the similarity of matrices using the PCAsimilarity function (EvolQG 0.2-9; Melo et al., 2015) (Table 2). Most of the above functions provide point estimates to describe attributes of a single matrix or attributes of a comparison of matrices, however, point estimates do not allow for statistical inferences. We therefore used bootstrap resampling (iterations=100) to estimate 95% confidence intervals around point estimates and to make direct comparisons between the P matrices of each focal group. The exception to this is the results from comparing Vrel , which natively performs a two-sample z-test to determine significance.
Finally, to make inferences about the genetic architecture and relationship between craniofacial traits, we used the variance and covariance values from the P matrices of Crescent Pond and Little Lake F2 hybrids to compare: 1) estimates of variance for each population, 2) covariation across traits within and between populations, 3) regression coefficients across traits within and between populations and 4) squared correlation coefficients across traits within and between populations. We used a paired t-test to investigate differences in variance between populations and used linear models with either covariation, regression coefficients, or squared correlation coefficients as the response variable, and population, trait, and their interaction as fixed effects. We used AIC scores to compare linear models including and excluding the interaction term and moved forward with the best fitting model.