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.