\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[ngerman,greek,english]{babel}



\usepackage{float}








\begin{document}

\title{Phenotypic covariation predicts diversification in an adaptive radiation
of pupfishes}



\author[1]{Julia Dunker}%
\author[2]{Michelle St. John}%
\author[1]{Christopher Martin}%
\affil[1]{UC Berkeley}%
\affil[2]{University of Oklahoma}%


\vspace{-1em}



  
  \date{May 25, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
Phenotypic covariation among suites of traits may constrain or promote
diversification both within and between species, yet few studies have
empirically investigated this relationship. In this study we investigate
whether phenotypic covariation of craniofacial traits is associated with
diversification in an adaptive radiation of pupfishes found only on San
Salvador Island, Bahamas (SSI). The full radiation of pupfishes includes
generalist, durophagous, and lepidophagous pupfish species. We compare
phenotypic variation and covariation (i.e., P matrix) between 1)
allopatric populations of generalist pupfish from neighboring islands
and estuaries in the Caribbean, 2) SSI pupfish populations not
containing the full radiation of fishes, and 3) SSI pupfish populations
containing the full radiation in sympatry. Additionally, we interrogate
the patterns observed in the P-matrices of two independent F2 hybrid
crosses derived from the two most morphologically distinct members of
the radiation to make inferences about the underlying mechanisms
contributing to the variation in craniofacial traits in SSI pupfishes.
We found that the P matrix of SSI generalist populations not found in
sympatry with specialists exhibited higher levels of mean trait
correlation, constraints, and integration with simultaneously lower
levels of flexibility compared to generalist populations on other
Caribbean islands and sympatric populations of all three species found
on SSI. We also document that while many craniofacial traits appear to
be produced via additive genetic effects, variation in key traits such
as head depth, maxilla length, and lower jaw length may be produced via
non-additive genetic mechanisms. Ultimately, this study suggests that
differences in phenotypic covariation significantly contribute to
producing and maintaining organismal diversity.%
\end{abstract}%



\sloppy


\textbf{Introduction}

Understanding the factors that produce and maintain diversity is a
fundamental goal of evolutionary biology. Genetic and environmental
variation interact to produced observed phenotypic variation, and many
studies have successfully identified the specific genes, alleles, or
selective pressures which have resulted in divergence between closely
related groups (Occhialini et al., 2012; Zhan et al., 2010, Gao et al.,
2009; Colosimo et al., 2005; Steingr\selectlanguage{ngerman}ímsson et al., 2006; Torres-Dowdall
et al., 2012; Endler, 1995; Martin et al., 2019; Winchell et al., 2016;
Gomi et al., 2007). Yet, many organisms, such as \emph{Howea} palms,
Andean \emph{Coeligena} hummingbirds, \emph{Amphilophus} sp. cichlids,
and \emph{Vidua} indigobirds, display high degrees of phenotypic
variation and diversification but exhibit limited genetic variation or
live in similar environments where selective pressures are assumed to be
the same (Palacios et al., 2019; Papadapoulos et al., 2011; Barluenga et
al., 2006; Savolainen et al., 2006; Sorenson et al., 2003). The presence
of phenotypic variation and diversification across several taxa, in the
absence of obvious genetic or environmental divergence, suggests that
the relationship between diversity and additional axes of variation
should be explored.

The observed relationships between traits (hereafter referred to as
phenotypic variation and covariation) is an additional complex phenotype
which can provide valuable insight into diversification. Phenotypic
variance and covariance are typically quantified as a matrix (hereafter
referred to as the P matrix) which is an established method for
determining the multivariate phenotypes present within a population
(Cheverud et al., 1989; Lande \& Arnold, 1983). Estimates and
comparisons of P matrices have been used to make inferences about
various measures of organismal performance (e.g., feeding or locomotion)
and are particularly useful for quantifying levels of integration,
identifying modules, and understanding how populations will respond to
selective pressures (Pavlicev et al. 2009; Goswami and Polly 2010; Haber
2015; Kane and Higham 2015; Reichert and Höbel 2018). Although it was
previously assumed that P matrices were constant across time (Lande,
1979), more recent work--including empirical studies--now supports that
1) P matrices are subject to change, and 2) that evolutionary forces can
act on variation at the level of the P matrix to produce convergent or
divergent phenotypes (Steppan, Phillips, and Houle, 2002; Roff \&
Mousseau, 2005; Selz et al., 2014; Blankers et al., 2015; Evans et al.,
2017; Michaud et al., 2020). For instance, Závorka et al. (2017) found
weaker associations between traits in brown trout populations exposed to
a novel selective pressure (i.e., invasion of brook trout), which in
turn produced divergent phenotypes compared to unaffected populations.
Similarly, Kolbe et al. (2011) found that the P matrices of
Caribbean\emph{Anolis} lizard occupying the same habitat type (e.g.
trunk, trunk-crown, trunk-ground, and twig) were more similar to one
another than they were to more closely related species which occupied a
different habitat type. Together, these results suggest that
investigating the relationship between traits as an additional axis of
diversity may provide new avenues for uncovering the genetic and
environmental variation which produces it.

Investigating the relationship between variation in the P-matrix and
diversification can also provide insights into the underlying mechanisms
contributing to phenotypic variation. The P matrix has been used as a
proxy for the genetic variance-covariance matrix (i.e., G matrix) for
some time and is generally more easily attainable than direct estimates
of covariation between genes (Cheverud 1988; Roff 1995). Comparing P
matrices between groups can provide insight into if and how the genetic
architecture or shared developmental pathways of traits responds to
different selective pressures (Cheverud 1995; Schluter 1996; Marroig and
Cheverud 2001; Kolbe et al. 2011). Furthermore, investigating the P
matrix of hybrid offspring of divergent groups can provide more specific
information regarding the genetic basis of traits of interest because
the laws of independent assortment, segregation, and assumptions of
additivity provide a null expectation for how phenotypic variation
should shift across generations (Falconer 1996; Roff 1997). For
instance, departures from additive expectations in F1 hybrids and
backcrosses were used as evidence that non-additive genetic variation
was likely a large contributor to beak shape divergence in species
of\emph{Geospiza} finches. Deviations from null expectations observed in
hybrids has also been used to infer X-linkage of acoustic traits in
field crickets (Blankers, Lübke, and Hennig, 2015) and to infer the
contributions of gene-gene and gene-environment interactions in
diverging populations of flour beetle (Drury et al., 2013).

In summary, comparing P matrices, within and between species or groups,
can provide additional information on how we may expect populations to
respond to selective pressures. This information in turn can help us
make further inferences about the diversification process and can
provide us with clues as to the genetic architecture, developmental
pathways, or mechanisms which are responsible for variation in traits.
To gain these additional insights, we must empirically investigate 1) if
differences in variance and covariance between phenotypic traits is
associated with diversification and 2) if P matrix attributes of hybrids
deviate from the null expectations of the laws of independent
assortment, segregation, and assumptions of additivity.

The \emph{Cyprinodon} pupfish system is excellent for investigating
whether phenotypic covariation is associated with diversity for two
reasons: first, the pupfish system contains at least three species that
display extensive phenotypic divergence. \emph{Cyprinodon
variegatus}(hereafter referred to as the generalist pupfish) has an
extremely large range that stretches along the Atlantic coast from North
America, throughout the Caribbean, and into northern portions of South
America (Hildebrand, 1917). In the interior hypersaline lakes of San
Salvador Islands, Bahamas (hereafter referred to as SSI), however, there
is a radiation of pupfishes that not only includes the widespread
generalist pupfish species, but additional snail-eating (\emph{C.
brontotheroides} ) and scale-eating (\emph{C. desquamator} ) specialist
species. Previous work has documented behavioral, morphological, and
physiological diversity that characterizes each of these species, but
their most obvious axes of divergence is in their craniofacial
musculoskeletal elements (Martin \& Wainwright, 2013a; Martin, 2016; St.
John et al., 2020a, 2020b; Heras and Martin., 2022). Briefly, generalist
pupfish jaws are similar to other cyprinodontiformes, snail-eating
pupfish exhibit an expanded dorsal head of the maxillae, and
scale-eating pupfish have a significantly larger oral jaw apparatus
(Hernandez et al., 2018; Martin and Wainwright, 2011).

Second, despite the observed phenotypic diversity, obvious genetic and
environmental variation does not appear to completely explain the
observed pattern of diversity (Martin 2016). Although specialist pupfish
are endemic to the hypersaline lakes of SSI, generalist pupfish
populations are found across the Atlantic coasts of the Americas, in
lakes across other Caribbean islands, and even in allopatry on SSI,
suggesting that environmental differences alone are not solely
responsible for pupfish diversification (Martin, 2016; Martin \&
Wainwright, 2011). Furthermore, phylogenetic evidence suggests that
radiating and non-radiating populations of pupfish on SSI likely share a
common ancestor which resembles the Caribbean generalist pupfish (Martin
2016). While there is potentially adaptive genetic differentiation
between generalist, snail-eating, and scale-eating pupfish, a recent
study found similar levels of genetic diversity between radiating and
non-radiating lineages of pupfish and that nearly all of the adaptive
genetic variation associated with specialization on SSI is found in
other generalist-only populations across the Caribbean (Richards \&
Martin, 2017; Richards et al., 2021; Patton et al., 2022; Richards \&
Martin, 2022). The incredible diversification rates of the pupfish
radiation on SSI, paired with our current understanding of the available
genetic and environmental variation, indicate that this system is a good
candidate for investigating the relationship between phenotypic
covariation and diversification.

In this study, we use a comparative framework to 1) determine if
radiating lineages on SSI display unique multivariate phenotypes and
covariation between traits which may have promoted their
diversification, and 2) compare multivariate phenotypes and covariation
among F2 hybrid offspring to make inferences about the underlying
mechanisms of craniofacial traits in pupfishes. We calculated and
compared variance-covariance matrices for 18 craniofacial traits for 1)
allopatric populations of generalists from neighboring islands and
estuaries across the Caribbean, 2) SSI generalist populations that do
not occur in sympatry with specialist species, and 3) sympatric lake
populations of all three species found on SSI, to address our first
question. We further calculated variance-covariance matrices for F2
hybrid offspring of scale-eating and snail-eating crosses from two
radiating populations of pupfish on SSI to address our second question.
We predict that sympatric populations containing the full radiation of
pupfishes on SSI will differ in their multivariate phenotypes and
covariation structure between traits from other pupfish groups.
Additionally, we predict that the multivariate phenotypes and the
covariation between traits observed in F2 hybrids will not deviate from
the null expectations derived from the laws of independent assortment,
segregation, and assumptions of additivity.

\textbf{Methods}

\emph{Focal Populations}

The goal of this study was to compare the P matrices of allopatric
generalist pupfish populations to sympatric populations of pupfish
containing both generalist and specialist species, and to estimate the
P-matrices of F2 hybrid pupfish to make inferences about the underlying
mechanisms of craniofacial traits on SSI. To that end, we measured and
compared craniofacial traits of fish from: five Caribbean populations of
pupfish that contained only generalist species (Lake Cunningham (New
Providence Island, Bahamas), Flamingo Pond (Acklins, Bahamas), Lake
George (Rum Cay Island, Bahamas), Etang Saumatre (Dominican Republic),
and Laguna Oviedo (Dominican Republic)); six SSI populations that do not
contain the full radiation of pupfishes, and instead contain only or
mostly generalist species (Wild Dilly Pond, Reckley Field Pond, Pain
Pond, Moon Rock Pond, Six Pack Pond, and Mermaid Pond); and three SSI
populations of pupfish that contained the full radiation of pupfishes in
sympatry (Crescent Pond, Little Lake, and Oyster Pond), although only
generalist and scale-eater specimens were available for measurement from
Oyster Pond (Table 1; Figure 1). A portion of the specimens from the
above ponds were initially collected, measured, and analyzed in a
previously published article (Martin, 2016), and complete collection
details can be found there. Additional justification and information
about the use and categorization of these populations can be found in
appendix A. For simplicity, we collectively refer to the five Caribbean
populations of pupfish as the `Caribbean' group/population (N = 61), the
six populations that do not contain the full pupfish radiation as the
`generalist-only' group/population (N=85), and the three populations
containing all three pupfish species in sympatry as the `radiating'
group/population (N=42).

To make inferences about the potential mechanisms underlying
craniofacial traits on SSI, we measured traits of F2 scale-eater X
snail-eater hybrid offspring and estimated phenotypic variation and
covariation of these traits. We produced F2 hybrids by first crossing a
single male snail-eater with a single female scale-eater to produce F1
offspring. We repeated this process independently for the Crescent Pond
and Little Lake populations. At least four F1 offspring from each
population were then allowed to interbreed to produce F2 hybrids. In
total we produced and measured 301 F2 hybrids from Crescent Pond and 194
hybrids from Little Lake (Table 1). These measurements were used for
separate QTL mapping study of these crosses (St. John et al. 2021).

When comparing F2 hybrid phenotypes and P matrices we assume a simple
additive model of inheritance (Falconer 1996; Roff 1997). The
assumptions of this model include: 1) that there are two alleles per
locus, 2) that mendelian laws of segregation are adhered to, 3) that
loci across the genome are in linkage equilibrium, and 4) that there are
only additive genetic effects (i.e., no dominance or epistatic effects).
With these assumptions we expect that the distribution of phenotypic
traits in the F2 generation should follow a 1:2:1 ratio, where
intermediate phenotypes are most common and phenotypes aligning with
either parental phenotype are less common and equally distributed, and
we expect covariation between traits to correspond to recombination
events, which for the sake of simplicity we assume to uniformly
distributed across the genome. To investigate deviations from these
assumptions we estimate expectations of additivity by calculating the
average of parental traits from the F0 generation and calculate
expectations of variation within traits using the parental trait values
and hybrid population sample sizes (Appendix E: Tables E1, E2, \& E3).

\emph{Clearing and staining}

Clearing and staining specimens promotes observation of specimens'
skeletal structures and specifically allows us to measure internal
skeletal traits. We therefore cleared and stained 683 fish specimens, in
preparation for future measurements. We fixed pupfish specimens in 95\%
ethanol, skinned them, and then immersed the specimens in 5\% buffered
formalin and double-stained the specimens with alizarin red and alcian
blue as outlined in (Dingerkus \& Uhler, 1977). After, we suspended the
specimens in glycerin and took lateral and dorsal photos with a Canon
EOS 60D digital SLR for downstream digital measurements.

\emph{Collection of Morphological Data}

We measured morphological traits to compare the phenotypes and P
matrices of pupfish across species and populations. We measured 18
craniofacial traits detailed in Martin et al. (2016) (Figure 2). We
specifically focused on craniofacial traits for this study because jaw
diversity is the primary axis of diversification in this system and is
hypothesized to be adaptive for each specialist's unique ecological
niches (Martin et al., 2017; Hernandez et al., 2018). We measured traits
using the program DLTdv8a (Hedrick, 2008), which outputs X and Y
coordinates for each landmark and subsequently calculated linear
distances from these coordinates. We standardized each cranial trait
measurement by performing a linear model with each trait measurement as
the response variable and log-transformed standard length as the
predictor variable using the lme4 package in R 4.1.1 (Bates et al. 2014,
R Core Team, 2021). We extracted the residuals for each individual and
used these values in all further analyses.

\emph{The effect of sample size on estimates of covariation}

While the sample sizes presented in this study are modest, they are
within the ranges suggested as sufficient in the literature. For
example, Cheverud's analysis (1988) suggests that 40 individuals is a
large enough sample size for estimating genetic variance from phenotypic
variance. More recent simulation-based studies suggest that adequate
sample sizes likely vary depending on the statistics being used and the
properties of one's dataset, but show that samples sizes of 16-32 are
large enough to accurately estimate most common statistics associated
with variance-covariance datasets (Grabowski and Porto 2017; Watanabe
2022). There are also several previously published empirical studies
examining the properties of variation-covariation with samples sizes
similar to our own (Polly 2005; Goswami and Polly 2010; Blankers et al.
2015). To further ensure that our sample size was adequate for detecting
differences in means between groups we performed power analyses to
calculate the minimum sample size required across a range of
correlations and effect sizes with 80\% power. We also performed Monte
Carlo simulations across a range of hypothetical trait values and
standard deviations to determine at which sample sizes estimates of
variation stabilized. Full methodologies and results can be found in the
appendix B, but overall, we found that our samples sizes were sufficient
for detecting differences between groups even at small effect sizes
(\selectlanguage{greek}η\selectlanguage{english}\textsuperscript{2}~\textgreater{}0.02) and that estimates of trait
variation stabilized at 20 individuals.

\emph{Investigating differences in multivariate phenotypes between
groups and hybrid populations}

We used a Bayesian mixed-effects model to investigate if trait means
varied between Caribbean, SSI generalist-only, and SSI radiating groups
while accounting for the effects of species and population. We used the
brms function from the BRMS package to fit a model with a student's T
distribution, and uninformative priors (B\selectlanguage{ngerman}ürkner 2017). The model
included the 18 trait measurements as the multivariate response
variable, group type and species designation as fixed effects, and
population ID as a random effect.

We fit a Bayesian fixed-effects model with a student's T distribution,
and uninformative priors to investigate differences in trait values
between Crescent Pond and Little Lake F2 hybrids. This model included
the 18 trait measurements as the multivariate response variable and
population as a fixed effect. While testing for differences in trait
values between hybrid groups is not a primary goal of this study,
detected differences could indicate population level differences in the
genetic architecture of these traits. Furthermore, the output of this
model provides estimates of uncertainty around average trait values
which were used to investigate departures from null expectations.
Specifically, we used these estimates to determine which phenotypes
deviated from additive trait value expectations (Appendix Table E4)

For both models we extracted posterior samples and estimated the median
values and 95\% highest density intervals (HDI) to make inferences about
differences in trait values between the levels of our fixed effects
using functions from the brms and emmeans packages (Lenth 2023; Figures
3 \& 7).

\emph{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 V\emph{rel} ) 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
V\emph{rel} , 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.

\textbf{Results}

\emph{SSI radiating populations have larger craniofacial trait values
than Caribbean or SSI generalist-only populations.}

The results of Bayesian mixed-effects model suggest that on average the
SSI radiating group has larger trait values (mean trait value: 0.06,
95\% HDI: 0.03, 0.09) than either the Caribbean (mean trait value:
-0.09, 95\% HDI: -0.13, -0.05) or SSI generalist-only groups (mean trait
value: -0.07, 95\% HDI: -0.11, -0.03; Figure 3). Model estimates for
each trait indicate that this pattern was reflected in 13 of the 18
traits, while 5 of the 18 traits showed no significant differences
between groups (Appendix Table D1).

We also included species as a fixed effect in this model to account for
the expected differences in trait values between species. We found that
on average generalists (mean trait value: 0.039, 95\% HDI: 0.020, 0.058)
had larger trait values than both snail-eaters (mean trait value: -0.11,
95\% HDI: -0.17, -0.046) and scale-eaters (mean trait value: -0.031,
95\% HDI: -0.075, 0.013; Appendix Figure D1.A), however, this pattern is
variable depending on the trait. For instance, scale-eaters have larger
estimates of lower jaw length and lower estimates of maxillary head
protrusion compared to both generalists and snail-eaters (Appendix Table
D2).

Finally, we investigated the effects of population in our model and
found that populations did not have significantly different trait value
estimates (Appendix Figure D1.B). We also examined the variability of
estimates due to population ID for each craniofacial trait. The standard
deviation reported for 16 of the 18 traits was below 0.1 indicating that
populations generally have similar trait values. The exceptions being
maxillary head height and nasal tissue protrusion, which had higher
standard deviations of 0.25 and 0.16 respectively (Appendix Figure
D1.C).

\emph{SSI generalist-only populations exhibit unique P matrix properties
that are not observed in Caribbean or SSI radiating populations.}

The visualization of the PCA suggests that the P matrix of SSI
generalist-only populations is distinct from both the Caribbean and SSI
radiating groups (Figure 4). In direct comparisons, we found that the
SSI radiating, and Caribbean groups had the highest level of PCA
similarity indicating that the P matrices of these two populations are
extremely similar (Figure 5; Appendix Table C1). Direct comparisons of
each of these groups with SSI generalist-only populations also indicates
that their P matrices are significantly different. Specifically, the P
matrix of the SSI generalist-only group has significantly higher levels
of integration, constraints and mean squared correlation between traits,
with simultaneously lower estimates of flexibility than the P matrices
of SSI radiating and Caribbean populations (Figure 6; Appendix Table
C2). In fact, there are only two P matrix features for which the SSI
generalist-only group is similar to other pupfish populations:
respondability and autonomy. For respondability, both SSI radiating and
SSI generalist-only populations exhibited similar values
(respondability= \textasciitilde{}0.09), while the Caribbean population
value was approximately half of this estimate (respondability= 0.048).
On the other hand, autonomy estimates for the Caribbean and SSI
generalist-only group were extremely similar
(autonomy=\textasciitilde{}0.07), while the P matrix of SSI radiating
populations had significantly lower values (autonomy=0.034) (Figure 6;
Appendix Table C2).

\emph{Differences in patterns of variation and covariation for specific
traits}

Our PCA investigating patterns across the P-matrices of Caribbean, SSI
generalist-only, and SSI radiating groups indicates that jaw closing
in-lever, maxilla length, and palatine height have strong negative
correlations with PC1(\textasciitilde{} -0.96 for all three traits),
while nasal tissue protrusion has the only positive correlation value
for PC1 (0.014; Figure 4; Appendix Table D3, Appendix Figure D2). For PC
axis 2 we observed a strong positive correlation associated with
maxillary head protrusion (0.91), moderate positive correlations for jaw
opening in-lever (0.56) and head depth (0.51), and moderate negative
correlations for lower jaw length (-0.59) and dentigerous arm width
(-0.43; Figure 4; Appendix Table D3, Appendix Figure D2). Similar
correlation values between traits along respective PC axes suggest that
these traits share similar patterns of variation and covariation, while
similar values with opposite signs indicate antagonist patterns. Further
investigation shows that most of the variation along PC1 is contributed
by the SSI generalist-only P matrix (52.5\%), while most variation along
PC2 is contributed by the SSI radiating P matrix (84.8\%; Appendix Table
D4).

Along PC1, we see that the major loadings of jaw closing in-lever,
maxilla length, and palatine height have a strong relationship in the
negative direction with SSI generalist-only traits, a variable
relationship with SSI radiating traits, and strong relationship in the
positive direction with Caribbean traits (Figure 4; Appendix Table D4).
This suggests that as the covariation of jaw closing in-lever, maxilla
length, and palatine height increases the covariation of other traits in
SSI generalist-only groups also increase, while the covariation in
Caribbean groups decreases. The relationship is more complicated in SSI
radiating groups as the covariation of some traits, such as lower jaw
length and dentigerous arm width, have a positive relationship with the
covariation of jaw closing in-lever, maxilla length, and palatine
height, while others such as nasal tissue protrusion and maxillary head
protrusion have a negative relationship.

The major loadings for PC2 include maxillary head protrusion, lower jaw
length, and jaw opening in-lever (Appendix Table D3). In general, the
variation along PC2 corresponds to the SSI radiation P matrix (Figure 4;
Appendix Table D4). For instance, for the SSI radiating group the
covariation between maxillary head protrusion and jaw opening in-lever
have positive relationship with one another but a negative relationship
with lower jaw length, dentigerous arm width, and dentigerous arm base.
These patterns are either significantly weaker, or opposite in the
Caribbean and SSI generalist only matrices.

\emph{F2 hybrid crosses differ in their deviations from null
expectations but exhibit similar average trait values.}

The results of Bayesian model suggest that there is no difference in
average trait values between Crescent Pond (mean trait value: -0.00096,
95\% HDI: -0.0044, 0.0023) and Little Lake (mean trait value: 0.00032,
95\% HDI: -0.0040, 0.0050) F2 hybrids. Estimates of 95\% HDI for each of
the 18 traits confirmed that there were no significant differences
between F2 populations (Appendix Table E4). Despite the lack of
differentiation in trait values between these two populations, we found
significantly different patterns of deviation from our null expectations
between Crescent Pond and Little Lake. For example, the estimates for 16
of the 18 craniofacial traits in Little Lake met the assumptions of
additivity, while only 7 of the 18 traits in Crescent Pond met the same
assumptions (Figure 7, Appendix Table E2). Specifically, we found that
estimates of cranial height and head depth were greater than expected in
Crescent Pond along with estimates of orbit diameter for both F2
populations. Estimates of maxillary head protrusion were smaller than
expected for Little Lake, and estimates of dentigerous arm width,
dentigerous arm base, lower jaw length, palatine height, jaw closing
in-lever, suspensorium length, nasal tissue protrusion and maxilla
length were all smaller than expected in Crescent Pond (Figure 7).

\emph{Similarities and Differences between for Crescent Pond and Little
Lake F2 Hybrid P-matrices}

The PCA visualization of F2 hybrid P matrices shows a large amount of
overlap between populations suggesting that their P matrices are quite
similar (Figure 8). This interpretation is supported by a large estimate
of PCA similarity between the two populations (\textasciitilde{}0.83),
and the relative similarity in mean squared correlation, autonomy, and
integration between populations (Appendix Tables C1 \& C2). Despite this
overall appearance of similarity there are still some distinct
differences between populations. For instance, variation in the P-matrix
of Crescent Pond aligns very closely with PC axis 2 while variation in
the Little Lake P-matrix is uniformly distributed across PC1 and PC2
(Figure 8). Furthermore, we found smaller estimates of respondability
(Crescent Pond: 0.022, Little Lake: 0.045) and flexibility (Crescent
Pond: 0.58, Little Lake: 0.62) for the Crescent Pond P matrix compared
to the Little Lake P matrix (Appendix Table C2). Together this suggests
that there are unique relationships between traits in Crescent Pond and
Little Lake.

\emph{Patterns of variance and covariance for specific traits}

Suspensorium length, palatine height, and ascending process length have
strong positive correlation with PC axis 1 (0.75,0.75, 0.71
respectively) while pelvic girdle length, dentigerous arm depth and
lower jaw length have strong negative correlations (-0.79, -0.74, -0.72
respectively; Appendix Table E5). PC2 is defined by positive
correlations with ectopterygoid (0.45) and jaw opening in-lever (0.22)
and strong negative correlations with head depth (-0.85), dentigerous
arm width (-0.77), and cranial height (-0.76; Appendix Table E5). Most
variation across both axes can be attributed to patterns in the Little
Lake P matrix (93\% of variance across PC1, and 63\% variance across
PC2), however, we can still identify clear differences in the
relationship between traits for each population. For instance, there is
a strong positive relationship between covariation patterns of lower jaw
length and jaw opening in-lever in Little Lake, but this relationship is
non-existent in Crescent Pond (Figure 8; Appendix Table E6). In fact,
many relationships between traits captured by PC1 are strong in Little
Lake but weak in Crescent Pond. Maxillary head height is the only trait
in Crescent Pond to show a strong positive relationship with head depth,
dentigerous arm width, and cranial height, while several traits in
Little Lake share this relationship (Figure 8).

Finally, we found that overall variation within traits was similar
between Crescent Pond and Little Lake (T test: t = -1.97, df = 27.45,
p-value = 0.059) but variances in Little Lake were greater than
variances in Crescent Pond for all traits. Although similar between
populations, only the variation in Crescent Pond met null expectations
associated with the laws of independent assortment and segregation
(expected variation: 0.14, observed variation: 0.013; T test: t = -0.14,
df = 33.72, p-value = 0.89). Variation in Little Lake was significantly
higher than expected (expected variation: 0.007, observed variation:
0.029; T test: t = 3.05, df = 21.37, p-value = 0.006). We found that
covariation estimates between traits were not significantly different
from one another (LM: F=1.39, df=17, p-value=0.13), but that overall
Little Lake had higher estimates of covariation than Crescent Pond (LM:
F=30.1, df=17, p-value\textless{}0.01).

The interaction between population and trait significantly affected
regression coefficient values (LM: F= 1.93, df=17, p-value=0.013).
Within Crescent Pond, head depth, maxilla length, and palatine height
had significantly higher regression coefficients than 50\% of measured
traits. Within Little Lake, head depth, lower jaw length, and
dentigerous arm width had higher regression coefficients than 50\% of
measured traits. Comparing traits between populations regression
coefficients for orbit diameter (Tukey's HSD: t ratio=2.67,
p-value=0.0078), palatine height (Tukey's HSD: t ratio=3.45,
p-value=0.0078), jaw closing in-lever (Tukey's HSD: t ratio=2.22,
p-value=0.027), and maxilla length (Tukey's HSD: t ratio=2.18,
p-value=0.003) were all higher in Crescent Pond than in Little Lake.
Lastly, we found that squared correlation coefficients significantly
varied between traits (LM: F=10.19, df=17, p-value\textless{}0.01) where
dentigerous arm width, maxilla length, lower jaw length, and palatine
height, had higher correlation coefficients than 50\% of traits in the
data set.

\textbf{Discussion}

\emph{Phenotypic covariation among traits significantly differs between
SSI radiating and non-radiating groups.}

The first aim of this study was to investigate if unique multivariate
phenotypes and covariation between traits was associated with pupfish
diversification on San Salvador Island, Bahamas. While pupfish have
radiated in several lakes on SSI, specialist species are notably absent
in other lakes on the same and nearby islands, despite their close
proximity and environmental similarity. We initially predicted that the
pattern of diversification would be explained by shifts in the mean
trait values and patterns of covariation in the SSI radiating group. Our
comparisons between mean trait values support this prediction, as SSI
radiating groups had larger trait values than SSI generalist-only and
Caribbean groups. This result was not surprising as the hallmark traits
of the specialist species within the radiation are the expanded dorsal
head of the maxillae and the larger jaw apparatus (Hernandez et al.,
2018; Martin and Wainwright, 2011).

Comparisons of the relationship between traits, however, suggest that
shifts in patterns of covariation primarily occur in the SSI
generalist-only group, whereas the Caribbean and the SSI radiation have
high levels of similarity between their P matrices (Figures 5 \& 6).
This could indicate that the P matrices of Caribbean populations of
pupfish are ``primed'' for diversification but have not been exposed to
the specific selective pressures that may induce such diversification.
On the other hand, SSI generalist-only populations are exposed to
similar environmental pressures, yet have not fully radiated, suggesting
that the phenotypic matrix structure for these populations lacks an
important quality for diversification. Indeed, we found that the P
matrices of these non-radiating populations differed from those of SSI
radiating and Caribbean populations (Figure 4), and our specific results
indicate that diversification may be limited by stronger associations
between traits, higher constraints, and lower flexibility (Figures 5 \&
6).

\emph{Estimates of integration, mean trait correlation, constraints, and
flexibility may limit pupfish diversification on San Salvador Island}

Our results suggest that greater independence in trait variation is
associated with diversification in this system. The P matrix of the SSI
generalist-only group had significantly greater estimates of integration
and mean squared correlation between craniofacial traits, higher
estimates of constraints, and lower estimates of flexibility-- all
measurements that can be interpreted as proxies for independence between
traits (Table 2; Figure 6). For example, we observed the highest
correlation between cranial traits in the SSI generalist-only group,
which indicates that the mechanisms producing shifts in trait values for
one trait are likely to affect the values for many others, which
ultimately may limit the phenotypic space that these populations can
access (Goswami \& Polly, 2010). Furthermore, the constraint and
flexibility estimates for the SSI generalist-only group suggest that
they are limited, specifically in comparison to Caribbean and SSI
radiating group, in their ability to respond to and align with a wide
range of selection gradients. This limited ability to respond to a range
of selection gradients, may explain why these populations are exposed to
similar environmental conditions, yet have not evolved the snail- and/or
scale-eater specialists found in other lakes.

Other empirical studies also suggest that independence across many
biological levels may be positively associated with diversification.
Greater amounts of gene duplication, larger ratios of
nonsynonymous:synonymous mutations, and lower levels of covariation
between morphological traits, have all been suggested as integral to
diversification in some African cichlid radiations (Brawand et al.,
2014; Machado et al., 2014; Selz et al., 2014), and represent
independence at both the genetic and morphological levels. Similarly,
Ravinet et al. (2014) suggest that populations of three-spined
stickleback from the Japan Sea have not colonized freshwater habitats
and subsequently diversified, as many of the closely related Pacific
freshwater populations have done, specifically due to the increased
correlation between their morphology and dietary niche. Lastly, the
flexible stem hypothesis suggests that the divergence within adaptive
radiations emerges from plastic phenotypic variation in an ancestral
population (West-Eberhard, 2003). Through this lens, the properties of
the Caribbean P-matrix could be viewed as being plastic (i.e., low
levels of constraints, high levels of flexibility, and low
respondability)---estimates which remain largely the same in radiating
groups and which have shifted in non-radiating groups. Other empirical
studies, however, have found mixed results regarding the relationship
between plasticity and diversification. Navolón et al. (2020) observed
that higher integration, which may imply lower levels of plasticity,
likely helped produce the diversity seen in the adaptive radiations of
Darwin's finches and Hawaiian honeycreepers. On the other hand, lower
levels of integration in brown trout allowed them to diversify into a
new ecological niche when faced with a new selective pressure (Závorka
et al., 2017). These results support that independence and flexibility,
either within or between, molecular, morphological, and/or behavioral
traits contributes to a group's ability to diversify. The findings
presented in this study further support this conclusion.

\emph{Phylogenetic relationships between populations of pupfish}

A caveat for the presented methods and results is that we did not
correct for phylogenetic relationships when estimating or comparing P
matrices across populations and species. Comparative methods are
important tools for making comparisons between samples that violate
assumptions of independence within comparison groups, such as comparing
closely related species that vary in their evolutionary relatedness to
each other. While the pupfish species and populations used in this study
are indeed closely related, we would argue that the patterns and
interpretations presented here are still informative for several
reasons. SSI generalists from both radiating and generalist-only ponds
have generally been shown to be more closely related to one another than
to specialist pupfish species-- especially the scale-eaters-- or to the
Caribbean pupfish populations used in this study (Richards et al. 2021,
Martin and Feinstein 2014, Martin 2016), a pattern that is reflected in
appendix figure A1. This information suggests that we would expect a
high level of similarity in the P matrices of SSI generalist-only and
SSI radiating groups due to their phylogenetic relationship, yet our
methods detected significant differences between these two groups. A
similar argument can be made to justify the interpretation that the P
matrix of the SSI generalist-only group is different from the P matrix
of the Caribbean group. Finally, the comparisons between the Caribbean
and SSI radiating groups indicates relative similarity between the two P
matrices. Additionally, the expectations of similarity in P matrices due
to evolutionary history are the same for generalist-only and radiating
groups. The fact that we detected differences when comparing SSI
generalist-only and Caribbean matrices suggests that the observed
similarity between the SSI radiation and Caribbean groups is not a false
positive.

\emph{Which trait relationships may promote or inhibit diversification
on San Salvador Island}

Within the above phylogenetic context, parsimony predicts that trait
relationships in the SSI radiating group which diverge from both
Caribbean and SSI generalist-only groups may be associated with
diversification. Indeed, we find that dentigerous arm width, lower jaw
length, palatine height, and maxilla length---traits that have high
correlations with one another in our PCA--- display either increasing or
decreasing strength in their relationships across pupfish groups which
supports this inference. To easily illustrate this pattern, we visualize
the R\textsuperscript{2} values from linear models with the covariation
values of dentigerous arm width, lower jaw length, palatine height, or
maxilla length as the response variable and the covariation values of
all other traits as the predictor variable and display the relationships
which may promote or inhibit divergence on SSI (Appendix Figure D3). For
example, dentigerous arm width, lower jaw length, palatine height, and
maxilla length all have higher associations with suspensorium length and
ectopterygoid in the Caribbean, weakened but intermediate relationships
in the SSI generalist-only group, and the weakest relationships in the
SSI radiating group. This pattern could suggest that increased
independence between these traits occurred on initial colonization of
SSI and further shifted during the radiating process. On the other hand,
all four traits have a weak relationship with maxillary head protrusion
in the Caribbean and stronger relationships on SSI.

The SSI generalist group also exhibits unique relationships between
dentigerous arm width, lower jaw length, palatine height, and maxilla
length and other traits that are not present in either the Caribbean or
SSI radiating group. We see a stronger relationship with other jaw
apparatus traits such as maxillary head height and the jaw closing and
opening in-levers in SSI generalist groups, but also see stronger
relationship with traits associated with general body shape and size
such as pelvic girdle length, orbit diameter, head depth, and cranial
height (Appendix Figure D3). We also see a stronger association between
maxilla length, palatine height, and nasal tissue protrusion in the
Caribbean and SSI radiation that is weakened in the SSI generalist-only
group.

Dentigerous arm width, lower jaw length, palatine height, and maxilla
length are traits with clear connections to feeding kinematics across
many fish species (Hulsey and García De León 2005; Westneat 2005; Muñoz
et al. 2018; DeLaurier 2019). Interpreting the patterns of associations
with these traits can provide us with potential modules which may be
important for feeding performance in SSI radiating groups and which
explain the lack of diversification in SSI generalist-only populations.
Four-bar linkages have been previously used to describe the sequence of
steps involved in fish feeding kinematics (Grubich \& Westneat, 2006;
Westneat 1990, 2005; Wainwright et al., 2004; Hulsey \& Wainwright,
2002), and our measurements of dentigerous arm width, lower jaw length,
maxilla length, and palatine height correspond to portions of the oral
four bar system. The adductor mandibulae complex also connects to the
premaxilla, maxilla and the dentary via several buccopalatal ligaments
aiding in the control of the jaw important for suction feeding
(Hernandez et al. 2009; Datovo and Vari 2013). Similarly, the
interoperculomandibular ligament is involved in opening the jaw during
the first stage of feeding and rotation around the joint of the lower
jaw and quadrate contributes to increasing gape opening during feeding
(Martin et al., 2017; Westneat, 2005). The similar patterns of
covariance between these traits along with their strong relationship to
one another suggests they may form a module in pupfishes. While all
these traits seem to act together across groups, the SSI generalist-only
group is the only one that displays a strong relationship between this
module and other traits such as head depth, maxillary head height, and
jaw opening in-lever. This could mean that these additional
relationships are constraining the jaw apparatus in SSI generalist-only
groups and preventing them from diversifying into additional ecological
niches.

Overall, we observed that the P matrix of generalist-only group on SSI
was distinct, and generally more constrained, than those of the
Caribbean and SSI radiating groups, a pattern that could be explained by
several mechanisms: First this pattern may be due to founder effects
(Barton \& Charlesworth, 1984). For example, if the founding population
of pupfish which first populated Wild Dilly Pond, Reckley Field Pond,
Pain Pond, Moon Rock Pond, Six Pack Pond, and Mermaid Pond on SSI had
more constrained P matrices than those that populated Crescent Pond and
Little Lake, then available variation may have been insufficient to
produce the scale-eating and/or snail-eating phenotypes. Alternatively,
the generalist-only ponds may have once contained the full adaptive
radiation, but the specialist species went extinct. We therefore may be
observing differences in the P matrix attributes due to the Baldwin
effect. The Baldwin effect predicts an initial burst of plasticity when
an environment changes, as in when pupfish first colonized SSI, and
subsequently predicts that plasticity will decrease as selection acts on
the population to bring them to a stable fitness optimum (Baldwin, 1896,
1902; Crispo, 2007). For instance, if selection in Wild Dilly Pond,
Reckley Field Pond, Pain Pond, Moon Rock Pond, Six Pack Pond, and
Mermaid Pond populations led to divergence into generalist, snail-eater
and scale-eater species and plasticity was then lost, an extinction
event of the specialists would leave the generalist species alone in the
pond with reduced plasticity. It should be noted that the Baldwin effect
can occur without species divergence, and the loss of plasticity may
occur following selection in the new environment. Further work needs to
be done to understand which scenario is most likely for this system.

\emph{Inferences about the underlying mechanisms of craniofacial traits
from F2 hybrid P-matrices}

In addition to making inferences about how variation in P-matrices may
contribute to diversification in the wild, this study also sought to
investigate the underlying mechanisms associated with these traits on
SSI by investigating the similarities and differences in P-matrices and
phenotype values of F2 hybrids from two separate crosses. We found that
6 of 18 traits met the null expectations of additivity in both
populations (Figure 7), suggesting that the mechanisms underlying these
traits generally fit the simple models often used in quantitative
genetics which assume Mendelian segregation, independent assortment, and
the presences of additive genetic effects (i.e., no dominance or
epistatic effects; Falconer 1996; Roff 1997). Of the remaining 12
traits, we found that only orbit diameter failed to meet additive
expectations in both crosses. Crescent Pond and Little Lake F2 hybrids
both exhibited larger mean estimates of this trait than additive
expectations, a pattern which could be produced through dominance or
epistatic effects. Ten of the remaining traits met null expectations of
additivity in Little Lake but not in Crescent Pond, and one trait,
maxillary head protrusion, met expectations in Crescent Pond but not in
Little Lake. This pattern could again be driven by dominance or
epistatic effects which are likely different between ponds, however, it
may also be driven by other factors such as new mutations in F2 hybrids,
recessive alleles, lethality of specific allele combinations (i.e.,
Dobzhansky-Muller or constitutive incompatibilities), or gene drive
(Rick and Smith 1953; Burkart-Waco et al. 2012; Chevin et al. 2014;
Andersen and Rockman 2022). There is indeed some evidence for this in
our current dataset. For example, many traits, especially in Little
Lake, have much higher levels of variation in phenotype trait values and
exhibit phenotypes that fall well outside those observed in the F0
parental types. Little Lake also shows higher levels of covariation
between traits than we observed in Crescent Pond. The increased
variation in trait values in both ponds may imply that new combinations
of alleles or epistatic effects are responsible for this additional
variation, while the high levels of covariation in Little Lake may
suggest that these traits share a common mechanism that is not present
in the Crescent Pond cross.

Increased variation and covariation within and between hybrid traits is
commonly observed in nature, and there is even evidence that
hybridization may drive speciation events in many systems including
pupfish (Seehausen 2004; Bell and Travis 2005; Selz et al. 2014;
Richards et al. 2019, 2021). On the other hand, F2 hybrids in Crescent
Pond do not display phenotype value that meet or exceed the scale-eater
parental values for traits such as palatine height, lower jaw length and
dentigerous arm width (Figure 7), implying that dominance, segregation
distortion, or lethality may be at play. While lethal combinations of
alleles is possible it is more commonly expected to occur when species
undergone divergence for longer periods of time than what is observed in
the pupfish system (Coyne and Orr 1997). Therefore, future work needs to
1) identify the loci associated with these traits and 2) investigate
whether these loci suffer from segregation distortion as a first step to
untangling the underlying mechanisms of these traits.

\emph{What can the F2 hybrid P-matrix tell us about craniofacial traits
on SSI?}

In this study we include calculations of the regression coefficient and
the squared correlation coefficient between traits. The regression
coefficients calculated in this study describe how much change is
observed in Trait B when Trait A increases by one unit. Values close 0
suggest that Trait A has little effect on Trait B, values close to
positive or negative one suggests that the magnitude of change A and B
are similar, and values outside of negative to positive one range
suggest that change in Trait B exceeds that of Trait A (although the
sign dictates if the change is in the same direction as Trait A; Kelly
2009). The formula for squared correlation coefficients is often used
between genetic loci to calculate the strength of association between
alleles (i.e. r\textsuperscript{2} as a measure of linkage
disequilibrium; VanLiere and Rosenberg 2008). In the context of this
study, high values suggests that individuals with large trait values for
Trait A also have large values for Trait B and vice versa. Interpreting
the patterns of both estimates together may give us insight into the
traits which are integral to producing the unique craniofacial pupfish
traits found on SSI. For instances, head depth in both Little Lake and
Crescent Pond produced the largest amount of change in other traits
which could mean that the underlying mechanisms responsible for
variation in this trait are reused across most other craniofacial traits
in pupfishes. Maxilla length and palatine height in Crescent Pond, and
lower jaw length and dentigerous arm width in Little Lake also produced
large changes in other traits further supporting this idea. The average
squared correlation coefficients were also significantly higher for
these traits compared to most of the remaining traits in the dataset.
From a mechanistic view point this pattern may reflect high levels of
linkage disequilibrium or it may suggest that these traits utilize a
shared underlying mechanism either through pleiotropy or shared genetic
or developmental pathways. Previous work in pupfishes has found
differential gene expression and fixed single nucleotide polymorphism
differences between species for regions of the genome associated with
the wnt signaling pathway which has ties to craniofacial diversity
across a wide array of organisms (Lencer et al. 2017; Lencer and Mccune
2020; Richards and Martin 2022). While much more work is needed to
determine if this specific pathway produces the observed craniofacial
diversity on SSI, these previous results support at least the
possibility of the reuse of a conserved evolutionary pathway as a
mechanism for producing pupfish diversity.

\textbf{Conclusion}

In this study, we sought to more fully understand which factors
contribute to diversity in the \emph{Cyprinodon} pupfish system. Our
data support the idea that phenotypic variation and covariation, or the
relationship between traits, is an additional and important axis of
variation that contributes to diversification in the pupfish system.
When investigating the factors driving diversification between pupfish
groups, we find that the lack of a flexible P matrix may impede the
ability to diversify. We also used F2 hybrid P-matrices to make
inferences about underlying mechanisms which produce variation in
craniofacial traits in the pupfish system. We found that many traits
adhere to the expectations of a simple additive genetic model. Yet, key
traits such as head depth, maxilla length, palatine height, lower jaw
length, and dentigerous arm width appear to be influences by nonadditive
genetic effects and produce large phenotypic changes in other traits.
Future work should examine whether this relationship between flexibility
and diversification is generalizable across many species or is specific
to \emph{Cyprinodon} pupfish.

\textbf{Acknowledgements}

We thank the Gerace Research Centre and Troy Day for logistical support
and the government of the Bahamas for permission to collect and export
samples. This research was funded by the National Science Foundation DEB
CAREER grant \#1749764, National Institutes of Health grant
5R01DE027052\selectlanguage{english}-02, Binational Science Foundation, the Miller Institute for
Basic Research in Science, the University of North Carolina at Chapel
Hill, and the University of California, Berkeley to CHM and an NSF
Postdoctoral Research Fellowship \#2209219 to MSJ. We also thank Dr.
Emilie Richards, Dr. Mallory Ballinger, the Martin Lab, and the Stein
Lab for giving feedback on the manuscript.

\textbf{Data Accessibility}

Upon acceptance, data will be deposited in dryad and scripts will be
posted on GitHub. For submission and review data and scripts can be
viewed here: https://github.com/stjohn3/Ecology\_Evolution\_Submision

\textbf{References}

Adams D, Collyer M, Kaliontzopoulou A, B. E. 2023. Geomorph: Software
for geometric morphometric analyses.

Andersen, E. C., and M. V. Rockman. 2022. Natural genetic variation as a
tool for discovery in Caenorhabditis nematodes. Genetics 220.

Baken, E. K., M. L. Collyer, A. Kaliontzopoulou, and D. C. Adams. 2021.
geomorph v4.0 and gmShiny: Enhanced analytics and a new graphical
interface for a comprehensive morphometric experience. Methods Ecol.
Evol. 12:2355--2363.

Bates, D., M. Machler, B. Bolker, and S. Walker. 2014. Fitting Linear
Mixed-Effects Models using lme4. J. Stat. Softw. 67:1--48.

Bell, M., and M. Travis. 2005. Hybridization, transgressive segregation,
genetic covariation, and adaptive radiation. Trends Ecol. Evol.
20:358--361.

Blankers, T., A. K. Lubke, and R. M. Hennig. 2015. Phenotypic variation
and covariation indicate high evolvability of acoustic communication in
crickets. J. Evol. Biol. 28:1656--1669. Blackwell Publishing Ltd.

Burkart-Waco, D., C. Josefsson, B. Dilkes, N. Kozloff, O. Torjek, R.
Meyer, T. Altmann, and L. Comai. 2012. Hybrid incompatibility in
Arabidopsis is determined by a multiple-locus genetic network. Plant
Physiol. 158:801--812.

Burkner, P.-C. 2017. brms\selectlanguage{english} : An R Package for Bayesian Multilevel Models
Using Stan. J. Stat. Softw. 80.

Cheverud, J. M. 1988. A comparison of genetic and phenotypic
correlations. Evolution (N. Y). 42:958--968.

Cheverud, J. M. 1995. Morphological Integration in the Saddle-Back
Tamarin (Saguinus fuscicollis) Cranium. Am. Nat. 145:63--89.

Chevin, L. M., G. Decorzent, and T. Lenormand. 2014. Niche
dimensionality and the genetics of ecological speciation. Evolution (N.
Y). 68:1244--1256.

Coyne, J. A., and H. A. Orr. 1997. ``Patterns of Speciation in
Drosophila'' Revisited. Evolution (N. Y). 51:295.

Datovo, A., and R. P. Vari. 2013. The Jaw Adductor Muscle Complex in
Teleostean Fishes: Evolution, Homologies and Revised Nomenclature
(Osteichthyes: Actinopterygii). PLoS One 8.

DeLaurier, A. 2019. Evolution and development of the fish jaw skeleton.

Falconer, D. S. 1996. Introduction to quantitative genetics. Pearson
Education India.

Goswami, A., and P. D. Polly. 2010. Methods for Studying Morphological
Integration and Modularity. Paleontol. Soc. Pap. 16:213--243.

Grabowski, M., and A. Porto. 2017. How many more? Sample size
determination in studies of morphological integration and evolvability.
Methods Ecol. Evol. 8:592--603.

Haber, A. 2015. The Evolution of Morphological Integration in the
Ruminant Skull. Evol. Biol. 42:99--114.

Hernandez, P. L., A. C. Gibb, and L. Ferry-Graham. 2009. Trophic
apparatus in cyprinodontiform fishes: Functional specializations for
picking and scraping behaviors. J. Morphol. 270:645--661.

Hulsey, C. D., and F. J. Garcia De Leon. 2005. Cichlid jaw mechanics:
Linking morphology to feeding specialization. Funct. Ecol. 19:487--494.

Kane, E. A., and T. E. Higham. 2015. Complex Systems Are More than the
Sum of Their Parts: Using Integration to Understand Performance,
Biomechanics, and Diversity. Integr. Comp. Biol. 55:146--165.

Kassambara, A., and F. Mundt. 2020. factoextra: Extract and Visualize
the Results of Multivariate Data Analyses.

Kelly, J. K. 2009. Connecting qtls to the G-matrix of evolutionary
quantitative genetics. Evolution (N. Y). 63:813--825.

Kolbe, J. J., L. J. Revell, B. Szekely, E. D. Brodie, and J. B. Losos.
2011. Convergent evolution of phenotypic integration and its alignment
with morphological diversification in caribbean Anolis ecomorphs.
Evolution (N. Y). 65:3608--3624.

Lencer, E., and A. R. Mccune. 2020. Differences in Cell Proliferation
and Craniofacial Phenotype of Closely Related Species in the Pupfish
Genus Cyprinodon. J. Hered. 111:237--247.

Lencer, E. S., W. C. Warren, R. Harrison, and A. R. McCune. 2017. The
Cyprinodon variegatus genome reveals gene expression changes underlying
differences in skull morphology among closely related species. BMC
Genomics 18:424.

Lenth, R. V. 2023. emmeans: Estimated Marginal Means, aka Least-Squares
Means.

Marroig, G., and J. M. Cheverud. 2001. A comparison of phenotypic
variation and covariation patterns and the role of phylogeny, ecology,
and ontogeny during cranial evolution of New World Monkeys. Evolution
(N. Y). 55:2576--2600.

Martin, C. H. 2016. The cryptic origins of evolutionary novelty:
1000-fold faster trophic diversification rates without increased
ecological opportunity or hybrid swarm. Evolution (N. Y). 70:2504--2519.

Munoz, M. M., Y. Hu, P. S. L. Anderson, and S. N. Patek. 2018. Strong
biomechanical relationships bias the tempo and mode of morphological
evolution. Elife 7.

Pavlicev, M., J. M. Cheverud, and G. P. Wagner. 2009. Measuring
morphological integration using eigenvalue variance. Evol. Biol.
36:157--170.

Polly, P. D. 2005. Development and phenotypic correlations: The
evolution of tooth shape in Sorex araneus. Evol. Dev. 7:29--41.

Reichert, M. S., and G. Hobel. 2018. Phenotypic integration and the
evolution of signal repertoires: A case study of treefrog acoustic
communication. Ecol. Evol. 8:3410--3429. John Wiley \& Sons, Ltd.

Richards, E. J., and C. H. Martin. 2022. We get by with a little help
from our friends: shared adaptive variation provides a bridge to novel
ecological specialists during adaptive radiation. Proc. R. Soc. B Biol.
Sci. 289.

Richards, E. J., J. A. McGirr, J. R. Wang, M. E. St. John, J. W.
Poelstra, M. J. Solano, D. C. O'Connell, B. J. Turner, and C. H. Martin.
2021. A vertebrate adaptive radiation is assembled from an ancient and
disjunct spatiotemporal landscape. Proc. Natl. Acad. Sci.
118:e2011811118.

Richards, E. J., M. R. Servedio, and C. H. Martin. 2019. Searching for
Sympatric Speciation in the Genomic Era. BioEssays 41:1900047.

Rick, C. M., and P. G. Smith. 1953. Novel Variation in Tomato Species
Hybrids. Am. Nat. 87:359--373.

Roff, D. A. 1997. Evolutionary Quantitative Genetics. Springer US,
Boston, MA.

Roff, D. A. 1995. The estimation of genetic correlations from phenotypic
correlations: a test of Cheverud's conjecture. Heredity (Edinb).
74:481--490.

Schluter, D. 1996. Adaptive radiation along genetic lines of least
resistance. Evolution (N. Y). 50:1766--1774.

Seehausen, O. 2004. Hybridization and adaptive radiation. Trends Ecol.
Evol. 19:198--207.

Selz, O. M., K. Lucek, K. A. Young, and O. Seehausen. 2014. Relaxed
trait covariance in interspecific cichlid hybrids predicts morphological
diversity in adaptive radiations. J. Evol. Biol. 27:11--24. John Wiley
\& Sons, Ltd.

VanLiere, J. M., and N. A. Rosenberg. 2008. Mathematical properties of
the r2 measure of linkage disequilibrium. Theor. Popul. Biol.
74:130--137.

Watanabe, J. 2022. Statistics of eigenvalue dispersion indices:
Quantifying the magnitude of phenotypic integration. Evolution (N. Y).
76:4--28.

Westneat, M. W. 2005. Skull Biomechanics and Suction Feeding in Fishes.

Tables

\textbf{Table 1:} \textbf{Sample sizes of measured individuals across
Caribbean populations, San Salvador Island generalist-only populations,
and San Salvador Island radiating populations.} Overall, we measured
craniofacial traits for 188 wild-caught individuals across the Caribbean
(including San Salvador Island and 496 lab-reared F2 hybrids. We
measured 61 individuals from five populations (i.e., `Caribbean'), 85
individuals from six populations on San Salvador Island that do not
contain the full radiation of pupfish species (i.e., `SSI
generalist-only'), and 42 individuals from three of the radiating
populations on San Salvador Island (i.e., `SSI radiation`).\selectlanguage{english}
\begin{longtable}[]{@{}llllll@{}}
\toprule
Population & Location & Generalist & Snail-eater & Scale-eater & F2
Hybrids\tabularnewline
\midrule
\endhead
Caribbean & Caribbean & Caribbean & Caribbean & Caribbean &
Caribbean\tabularnewline
Etang Saumatre & Dominican Republic & 16 & & &\tabularnewline
Flamingo Pond & Acklins, Bahamas & 7 & & &\tabularnewline
Lake Cunningham & New Providence Island, Bahamas & 10 & &
&\tabularnewline
Lake George & Rum Cay, Bahamas & 15 & & &\tabularnewline
Laguna Oviedo & Dominican Republic & 13 & & &\tabularnewline
Total: 61 & Total: 61 & Total: 61 & Total: 61 & Total: 61 & Total:
61\tabularnewline
SSI Generalist-only & SSI Generalist-only & SSI Generalist-only & SSI
Generalist-only & SSI Generalist-only & SSI
Generalist-only\tabularnewline
Mermaid Pond & San Salvador Island, Bahamas & 16 & & &\tabularnewline
Moon Rock pond & San Salvador Island, Bahamas & 21 & & &\tabularnewline
Pain Pond & San Salvador Island, Bahamas & 14 & & &\tabularnewline
Reckley Field Pond & San Salvador Island, Bahamas & 11 & &
&\tabularnewline
Six Pack Pond & San Salvador Island, Bahamas & 15 & & &\tabularnewline
Wild Dilly Pond & San Salvador Island, Bahamas & 8 & & &\tabularnewline
Total: 85 & Total: 85 & Total: 85 & Total: 85 & Total: 85
&\tabularnewline
SSI Radiation & SSI Radiation & SSI Radiation & SSI Radiation & SSI
Radiation & SSI Radiation\tabularnewline
Crescent Pond & San Salvador Island, Bahamas & 4 & 3 & 5 &
301\tabularnewline
Little Lake & San Salvador Island, Bahamas & 4 & 8 & 7 &
194\tabularnewline
Oyster Pond & San Salvador Island, Bahamas & 3 & & 8 &\tabularnewline
Total: 42 & Total: 42 & Total: 42 & Total: 42 & Total: 42
&\tabularnewline
\bottomrule
\end{longtable}

\textbf{Table 2:} \textbf{Descriptions of the different statistics used
to compare the P matrices of focal groups.} These statistics come from
the EvolQG package (Melo et al. 2015) and the geomorph package (Baken et
al. 2021; Adams D, Collyer M, Kaliontzopoulou A 2023), and complete
details can be found there. Briefly, the autonomy, constraint,
flexibility, mean squared correlation, respondability, and integration
statistics provide point estimates that serve as descriptions of a
single group's P matrix. On the other hand, PCA similarity statistics
provide estimates of similarity between two P matrices.\selectlanguage{english}
\begin{longtable}[]{@{}ll@{}}
\toprule
Statistic & Description\tabularnewline
\midrule
\endhead
Autonomy & Proportion of variance that is aligned with a selection
gradient and is separate from variation in other
directions\tabularnewline
Constraints & Mean correlation between the response vector to random
selection variants and the matrix's first principal component, which
largely constrains the phenotype due to it accounting for much of a
species' phenotypic variation\tabularnewline
Flexibility & Measure of how close the response of a species aligns with
the different random selection gradients\tabularnewline
Mean squared correlation & Average of the correlation coefficients,
corresponding to the overall association between the traits
measured\tabularnewline
Respondability & Pace of the change in population mean when a population
is under directional selection\tabularnewline
PCA similarity & Measure of similarity of the matrices' variation by
incormporating the relative similarity of the matrices' corresponding
principal components\tabularnewline
Integration & Measure of the strength of morphological integration
quantified using the relative eigen index\tabularnewline
\bottomrule
\end{longtable}

Figures\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/image1}

\end{center}
\end{figure}

\textbf{Figure 1:} \textbf{Map highlighting the location of study
populations and species.} Inset designated by the blue triangle provides
a large view of San Salvador Island, Bahamas. Circles filled with color
indicate species sampled from each location (green: generalist pupfish,
orange: scale-eating pupfish, blue: snail-eating pupfish). Grey circles
outlined in color represent species that may be present in each
population but were not sampled in this study. Map images from Google
Earth (2023).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image2/image2}

\end{center}
\end{figure}\selectlanguage{english}
\begin{longtable}[]{@{}llllll@{}}
\toprule
Landmark 1 & Landmark 2 & Trait & Landmark 1 & Landmark 2 &
Trait\tabularnewline
1 & 2 & Lower Jaw Length & 9 & 10 & Ascending Process
Length\tabularnewline
2 & 3 & Jaw Closing In-Lever & 11 & 13 & Maxillary Head
Height\tabularnewline
2 & 4 & Jaw Opening In-Lever & 11 & 14 & Ectopterygoid\tabularnewline
2 & 11 & Palatine Height & 12 & 13 & Maxillary Head
Protrusion\tabularnewline
2 & 18 & Suspensorim Length & 12 & 19 & Nasal Tissue
Protrusion\tabularnewline
5 & 8 & Dentigerous Arm Width & 14 & 15 & Orbit Diameter\tabularnewline
\bottomrule
\end{longtable}

\textbf{Figure 2:} \textbf{Landmark location and craniofacial trait
names displayed on the lateral view of a cleared and stained pupfish
specimen.} Landmarks from Martin (2016).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image3/image3}

\end{center}
\end{figure}

\textbf{Figure 3. The Multivariate trait values of SSI radiating pupfish
are on average larger than the multivariate trait values of SSI
generalist-only and Caribbean pupfish groups.} Violin plots display the
posterior draws from a Bayesian mixed effects model with the 18
craniofacial trait values as the response variable, pupfish group and
species designations as fixed effects and population ID as a random
effect. Black points and lines represent the median and 95\% HDI for
each group. Dots show average trait values for individuals.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image4/image4}

\end{center}
\end{figure}

\textbf{Figure 4: Results of a principal components analysis calculated
with the variance and covariance estimates for 18 craniofacial traits
for the Caribbean, SSI generalist-only, and SSI radiating
populations.}PC1 and PC2 explain 59\% and 16\% of the variance,
respectively. Shapes and colors indicate the three focal groups (dark
grey triangle: Caribbean group, \textbf{light grey} square: SSI
generalist-only group, gold circle: SSI radiating group), and ellipses
represent 95\% confidence intervals per group. Arrows display four
craniofacial traits that are most strongly correlated with PC1 (jaw
closing in-lever and lower jaw length) and PC2 (maxillary head
protrusion and nasal tissue protrusion). San Salvador Island
generalist-only group clusters separately from the Caribbean and San
Salvador Island radiating groups along PC1.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image5/image5}

\end{center}
\end{figure}

\textbf{Figure 5: Caribbean pupfish populations and SSI radiating
populations have the highest estimates of similarity, while SSI
generalist-only populations have the lowest estimates of
similarity.}Results of matrix comparisons using the PCA similarity
method. The squares along the longest diagonal show the mean values of P
matrix metrics for the corresponding comparisons. The remaining squares
show the differences between two comparisons. Fill colors represent low
(light purple) and high (dark purple) values within each metric.
Asterisks indicate that metrics are significantly different between
groups (Appendix Table C1).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image6/image6}

\end{center}
\end{figure}

\textbf{Figure 6:} \textbf{SSI generalist-only populations have lowest
estimate of flexibility, and the highest estimates of constraints, mean
squared correlation, and integration.} The squares along the longest
diagonal show the mean values of P matrix metrics for the corresponding
groups. The other squares show the differences between means of groups.
Fill colors represent low (light purple) and high (dark purple) values
within each metric. Asterisks indicate significant differences between
groups (Appendix Table C1).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image7/image7}

\end{center}
\end{figure}

\textbf{Figure 7.} \textbf{The estimated trait values of F2 hybrid
pupfish (Snail-eater X Scale-eater ) compared to F0 parental species
from Crescent Pond and Little Lake.} Grey dots forming a histogram
beneath each trait show the distribution of F2 phenotypes; The colored
density graph above each histograms depict the posterior draws from a
Bayesian model with the 18 craniofacial trait values as the response
variable and population as the fixed effect. Blue circles represent the
F0 snail-eater parental phenotype for each trait, while orange triangles
represent the F0 scale-eater parental phenotype for each trait. Finally,
grey diamonds depict the mid-parent trait value which represents the
expected phenotype trait value under the assumptions of additivity.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image8/image8}

\end{center}
\end{figure}

\textbf{Figure 8:} \textbf{Results of principal components analysis
including the variance and covariance estimates for 18 craniofacial
traits for F2 hybrids from Crescent Pond (blue triangles) and Little
Lake (green squares).} Colors represent species designations and
ellipses represent 95\% confidence intervals per group. Arrows display
four craniofacial traits that are most strongly correlated with PC1
(pelvic girdle length and suspensorium length) and PC2 (ectopterygoid
and head depth). Variation in Crescent Pond F2 hybrids aligns with PC2,
while variation in Little Lake is spread across both axes.

\textbf{Appendix}

\emph{Appendix A: Population groupings}

For comparisons, we formed three focal groups: SSI generalist-only, SSI
radiating, and the Caribbean. The PCA (Figure A1) indicates that
generalists, snail-eaters, and scale-eaters on SSI form distinct
clusters based on their genetic variation. The PCA visualization also
suggests that: 1) Oyster Pond scale-eaters and generalists fall well
within the 95\% confidence intervals for their species designations,
justifying their inclusion in the SSI radiating population analyses and
2) Mermaid Pond ``scale-eaters'' (designated as such based on their
morphological traits) do not fall within the 95\% confidence interval
for genetically clustering scale-eaters, and are instead closer to the
centroid of generalists. Although previous studies suggest that Mermaid
Pond contains scale-eating individuals (Martin and Wainwright. 2013),
more recent work suggests that these individuals are distinct from
scale-eaters (Richards et al. 2021. Richards and Martin. 2022).
Furthermore, there are no published sightings of snail-eating pupfish in
Mermaid Pond. For these reasons, we included Mermaid Pond in the
generalist-only population calculations as opposed to the SSI radiating
populations. Finally, snail-eaters have previously been documented in
Moon Rock and Wild Dilly Pond, but scale-eating specialists have never
been found in these locations (Martin 2014). Considering that we did not
have access to any of the snail-eater specimens from these locations,
and that these populations have never contained the full radiation of
pupfish,

we included Moon Rock and Wild Dilly Pond in the generalist-only
population calculations as opposed to the SSI radiating populations.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image9/image9}

\end{center}
\end{figure}

Figure A1. Results of a principal component analysis investigating the
genetic relationships between generalist, snail-eater, and scale-eater
pupfish species in lakes across San Salvador Island, Bahamas. Points
represent sequenced individuals (n=110), shapes represent species
designations based on morphological traits of each specimen
(generalists: circles, snail-eaters: triangles, scale-eaters: squares),
colors and labels correspond to 14 lakes on SSI that were sampled for
this dataset, and ellipses represent 95\% confidence intervals per
group. Figure modified from Richards et al. (2021).

\emph{Appendix B: Power analyses and Monte Carlo simulations}

To support our claim that samples sizes are sufficient, for detecting
differences between the groups of interest (i.e., SSI radiation, SSI
generalist-only, and Caribbean), we include the results of power
analyses, which calculate the required minimum sample size across the
range of effect sizes (eta2: 0.0002 - 0.12), number of individuals per
population (n=7-21), and correlations of trait estimates for individuals
observed in our dataset (corr: 0.00046- 0.50). From these iterative
power analyses, we were able to determine for which combination of our
parameters we would have sufficient sample sizes to detect differences
between groups with 80\% power. We compared these cut offs to the
observed values for each trait and population (Figure B1). Ultimately,
we found that we have enough samples to detect differences between
groups for all traits except for ascending process length, orbit
diameter, and cranial height.

We also performed Monte Carlo simulations to investigate at which sample
sizes estimates of variation stabilize based on observed means and
standard deviations of traits from our dataset. To perform these
simulations, we first calculated the average trait values (range: -0.71
- 0.42) and standard deviations (range: 0.01-0.8) for each trait across
the potential groups of our data set (19 groups representing unique
species per population, plus the additional three parameter estimates
associated with Caribbean, SSI generalist-only, and SSI radiating
groups). We used these observed means and standard deviations to
generate normal distributions to represent hypothetical populations to
sample from. We resampled from these distributions using sample sizes
from 1-100 for 1000 iterations, recalculating the mean each time. For
each of our 22 groups we ended up with 1800000 independent estimates of
means across the sample size range of 1-100. At each sample size we then
calculated the standard deviation as a metric representing the amount of
variation introduced due to sample size. As expected, as sample size
increased standard deviation decreased, however, the amount of variation
at low sample sizes was dependent on the original parameters of the
dataset (Figure B2 A\&B).

To determine at which sample sizes we observed a stabilization of
variation (i.e., when adding more samples did not significantly affect
the estimates of SD) we grouped sample sizes into sets of 5 (e.g., 1-5,
6-10, 11-15, etc.) and calculated the derivatives for each of these
groups across the sampling range. We used a one-sample t-test to
determine at which sampling range the derivative no longer significantly
differed from zero using a conservative cut off of non-significance
with\emph{p} =\textgreater{}0.1. Here, a derivative of zero represents
when variation in SD is no longer significantly affected by the addition
of more samples (i.e., when variation is stable). We determined at which
sample size we first observed non-deviance from zero for each of the
observed parameters from our dataset and visualized these results using
a histogram (Figure B2.C) Overall we found that the median sample size
where stabilization occurred was at 20 individuals (range: 5-35
individuals; Figure B2.C).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image10/image10}

\end{center}
\end{figure}

Figure B1: Visualization of the results of power analyses for a one-way
repeated measures analysis of variance. We used a range of different
correlations, effect sizes, and number of individuals within each
population to estimate the total sample sizes necessary to detect
differences between groups with 80\% power. The black line indicates the
cut off where the sample size of this dataset (n=188) is sufficient for
detecting differences between groups given the observed correlation
between traits and estimated effect size. Colored triangles represent
the mean correlation for each trait, while grey dots show the
correlations for each trait combination.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image11/image11}

\end{center}
\end{figure}

\textbf{Figure B2: Visualization of Monte Carlo simulations
investigating how sample size affects estimates of variation.} Panels A
and B display samples of simulations using observed means and standard
deviation observed for maxillary head height and head depth in our
dataset. The red line in each panel represents the mean trait value at
each sample size. The black ribbons show the minimum and maximum trait
value at each sample size, and the blue ribbons show the standard
deviation around the mean at each sample size. The vertical line in each
graph indicates where increasing sample size no longer significantly
impacts estimates of standard deviation (i.e., where estimates of
variation stabilize). Panel C shows histograms of the sample sizes
associated with stabilization across traits.

\emph{Appendix C: P-Matrix Statistical Analyses}

Table C1: PCA similarity estimates, including upper and lower 95\%
confidence intervals (bootstrapping; iterations=100). The PCA similarity
metrics estimate the similarity between two matrices ranging from 0 to
1.\selectlanguage{english}
\begin{longtable}[]{@{}ll@{}}
\toprule
Population & PCA Similarity\tabularnewline
SSI Radiation/Caribbean & 0.71 (0.69, 0.72)\tabularnewline
SSI Generalist-only/Caribbean & 0.61 (0.59, 0.62)\tabularnewline
SSI Generalist-only/SSI Radiation & 0.6 (0.59, 0.62)\tabularnewline
CP Hybrid/LL Hybrid & 0.83 (0.83, 0.84)\tabularnewline
\bottomrule
\end{longtable}

Table C2: Mean matrix statistics estimates, including upper and lower
95\% confidence intervals (bootstrapping; iterations=100), for focal
groups. Each of the six statistics represents a descriptive
characteristic of the P matrix for the indicated group. A full
description of each statistic can be found in Table 2 of the main text.\selectlanguage{english}
\begin{longtable}[]{@{}lllllll@{}}
\toprule
Population & Respondability & Autonomy & Mean Squared Correlation &
Flexibility & Constraints & Integration\tabularnewline
Caribbean & 0.048 (0.046, 0.049) & 0.07 (0.068, 0.072) & 0.15 (0.15,
0.16) & 0.51 (0.51, 0.52) & 0.61 (0.6, 0.62) & 0.19 (0.19,
0.2)\tabularnewline
SSI Generalist-only & 0.095 (0.09, 0.1) & 0.075 (0.073, 0.078) & 0.46
(0.44, 0.47) & 0.42 (0.42, 0.43) & 0.78 (0.77, 0.79) & 0.33 (0.32,
0.34)\tabularnewline
SSI Radiation & 0.099 (0.095, 0.1) & 0.034 (0.033, 0.036) & 0.14 (0.14,
0.15) & 0.53 (0.52, 0.53) & 0.6 (0.6, 0.62) & 0.17 (0.16,
0.17)\tabularnewline
CRP Hybrid & 0.022 (0.022, 0.022) & 0.15 (0.15, 0.15) & 0.11 (0.11,
0.11) & 0.58 (0.58, 0.59) & 0.56 (0.56, 0.57) & 0.13 (0.13,
0.13)\tabularnewline
LL Hybrid & 0.045 (0.045, 0.046) & 0.14 (0.14, 0.14) & 0.1 (0.098, 0.1)
& 0.62 (0.62, 0.62) & 0.54 (0.54, 0.55) & 0.1 (0.1, 0.1)\tabularnewline
\bottomrule
\end{longtable}

\emph{Appendix D: Estimates for Caribbean, SSI generalist-only, and SSI
radiating group analyses}

Table D1: Estimated marginal mean trait values and 95\% HDI for
Caribbean, SSI generalist only, and SSI radiating groups from a Bayesian
mixed effects model with the 18 craniofacial trait values as the
response variable, pupfish group and species designations as fixed
effects and population ID as a random effect.\selectlanguage{english}
\begin{longtable}[]{@{}llll@{}}
\toprule
Trait & Caribbean & SSI Generalist-Only & SSI Radiation\tabularnewline
Dentigerous Arm Width & -0.075 (-0.16, 0.007) & -0.076 (-0.16, 0.00044)
& 0.076 (-0.014, 0.17)\tabularnewline
Dentigerous Arm Base & -0.081 (-0.13, -0.025) & -0.067 (-0.13, -0.015) &
0.078 (0.035, 0.12)\tabularnewline
Lower Jaw Length & -0.14 (-0.18, -0.084) & -0.095 (-0.14, -0.049) &
0.062 (0.023, 0.11)\tabularnewline
Dentigerous Arm Depth & -0.16 (-0.23, -0.091) & -0.13 (-0.19, -0.062) &
0.061 (-0.00067, 0.12)\tabularnewline
Palatine Height & -0.067 (-0.12, -0.018) & -0.054 (-0.1, -0.0024) &
0.079 (0.03, 0.12)\tabularnewline
Jaw Closing In Lever & -0.04 (-0.11, 0.037) & -0.063 (-0.13, 0.0086) &
0.077 (0.0056, 0.14)\tabularnewline
Cranial Height & -0.12 (-0.17, -0.065) & -0.14 (-0.19, -0.083) & 0.015
(-0.021, 0.055)\tabularnewline
Orbit Diameter & -0.1 (-0.15, -0.057) & -0.11 (-0.15, -0.064) & 0.0025
(-0.043, 0.052)\tabularnewline
Suspensorium Length & 0.0014 (-0.083, 0.087) & -0.021 (-0.094, 0.063) &
0.086 (0.0039, 0.18)\tabularnewline
Ascending Process Length & -0.049 (-0.17, 0.074) & -0.073 (-0.18, 0.051)
& 0.011 (-0.1, 0.12)\tabularnewline
Jaw Opening In Lever & -0.13 (-0.24, -0.012) & -0.072 (-0.18, 0.041) &
0.15 (0.059, 0.25)\tabularnewline
Nasal Tissue Protrusion & -0.027 (-0.29, 0.22) & 0.077 (-0.16, 0.32) &
-0.06 (-0.29, 0.19)\tabularnewline
Ectopterygoid & -0.11 (-0.22, 0.0047) & -0.048 (-0.15, 0.066) & 0.076
(-0.025, 0.18)\tabularnewline
Head Depth & -0.11 (-0.15, -0.052) & -0.094 (-0.14, -0.045) & 0.0015
(-0.053, 0.057)\tabularnewline
Maxillary Head Protrusion & -0.089 (-0.19, 0.021) & -0.026 (-0.13,
0.074) & 0.14 (0.043, 0.24)\tabularnewline
Maxilla Length & -0.064 (-0.11, -0.012) & -0.061 (-0.11, -0.013) & 0.076
(0.032, 0.12)\tabularnewline
Maxillary Head Height & -0.18 (-0.51, 0.1) & -0.2 (-0.49, 0.091) & 0.19
(-0.11, 0.51)\tabularnewline
Pelvic Girdle Length & -0.057 (-0.13, 0.0099) & -0.055 (-0.12, 0.0095) &
-0.027 (-0.085, 0.03)\tabularnewline
\bottomrule
\end{longtable}

Table D2: Estimated marginal mean trait values and 95\% HDI for
generalist, snail-eater, and scale-eater pupfish species from a Bayesian
mixed effects model with the 18 craniofacial trait values as the
response variable, pupfish group and species designations as fixed
effects and population ID as a random effect.\selectlanguage{english}
\begin{longtable}[]{@{}llll@{}}
\toprule
Trait & Generalists & Snail-Eaters & Scale-Eaters\tabularnewline
Dentigerous Arm Width & 0.041 (-0.0044, 0.087) & -0.22 (-0.33, -0.13) &
0.11 (0.032, 0.19)\tabularnewline
Dentigerous Arm Base & 0.038 (0.014, 0.063) & -0.21 (-0.3, -0.13) & 0.11
(0.041, 0.17)\tabularnewline
Lower Jaw Length & 0.037 (0.014, 0.06) & -0.34 (-0.41, -0.27) & 0.13
(0.08, 0.19)\tabularnewline
Dentigerous Arm Depth & 0.056 (0.023, 0.089) & -0.28 (-0.36, -0.18) &
-0.0061 (-0.081, 0.06)\tabularnewline
Palatine Height & 0.034 (0.0092, 0.058) & -0.17 (-0.25, -0.1) & 0.094
(0.039, 0.15)\tabularnewline
Jaw Closing In Lever & 0.041 (0.0056, 0.078) & -0.11 (-0.21, -0.0032) &
0.042 (-0.036, 0.12)\tabularnewline
Cranial Height & 0.053 (0.033, 0.077) & -0.15 (-0.23, -0.064) & -0.14
(-0.2, -0.083)\tabularnewline
Orbit Diameter & 0.036 (0.011, 0.06) & -0.14 (-0.2, -0.074) & -0.11
(-0.16, -0.062)\tabularnewline
Suspensorium Length & 0.0089 (-0.034, 0.054) & -0.024 (-0.12, 0.082) &
0.082 (0.0034, 0.18)\tabularnewline
Ascending Process Length & 0.039 (-0.022, 0.097) & -0.17 (-0.34, 0.0085)
& 0.016 (-0.12, 0.15)\tabularnewline
Jaw Opening In Lever & 0.038 (-0.015, 0.089) & -0.024 (-0.19, 0.15) &
-0.064 (-0.19, 0.06)\tabularnewline
Nasal Tissue Protrusion & -0.036 (-0.17, 0.083) & 0.18 (-0.18, 0.52) &
-0.15 (-0.42, 0.094)\tabularnewline
Ectopterygoid & 0.051 (-0.004, 0.1) & -0.098 (-0.25, 0.064) & -0.032
(-0.16, 0.077)\tabularnewline
Head Depth & 0.038 (0.011, 0.068) & -0.1 (-0.17, -0.044) & -0.13 (-0.18,
-0.085)\tabularnewline
Maxillary Head Protrusion & 0.056 (0.0058, 0.11) & 0.24 (0.097, 0.39) &
-0.27 (-0.38, -0.16)\tabularnewline
Maxilla Length & 0.04 (0.016, 0.064) & -0.17 (-0.24, -0.095) & 0.078
(0.022, 0.13)\tabularnewline
Maxillary Head Height & 0.093 (-0.072, 0.25) & -0.072 (-0.45, 0.27) &
-0.21 (-0.51, 0.06)\tabularnewline
Pelvic Girdle Length & 0.031 (0.0022, 0.064) & -0.086 (-0.19, 0.019) &
-0.084 (-0.16, -0.013)\tabularnewline
\bottomrule
\end{longtable}

Table D3: Correlation of Loadings with PC axes one and two from a
principal components analysis calculated with the variance and
covariance estimates for 18 craniofacial traits for the Caribbean, SSI
generalist-only, and SSI radiating populations.\selectlanguage{english}
\begin{longtable}[]{@{}lll@{}}
\toprule
Trait & PC1 & PC2\tabularnewline
Dentigerous Arm Width & -0.86 & -0.43\tabularnewline
Dentigerous Arm Base & -0.90 & -0.33\tabularnewline
Lower Jaw Length & -0.80 & -0.59\tabularnewline
Dentigerous Arm Depth & -0.89 & -0.36\tabularnewline
Palatine Height & -0.96 & -0.24\tabularnewline
Jaw Closing In Lever & -0.97 & -0.091\tabularnewline
Cranial Height & -0.70 & 0.47\tabularnewline
Orbit Diameter & -0.91 & 0.29\tabularnewline
Suspensorium Length & -0.86 & 0.15\tabularnewline
Ascending Process Length & -0.69 & 0.034\tabularnewline
Jaw Opening In Lever & -0.58 & 0.56\tabularnewline
Nasal Tissue Protrusion & 0.14 & 0.28\tabularnewline
Ectopterygoid & -0.75 & 0.40\tabularnewline
Head Depth & -0.83 & 0.51\tabularnewline
Maxillary Head Protrusion & -0.074 & 0.91\tabularnewline
Maxilla Length & -0.97 & -0.21\tabularnewline
Maxillary Head Height & -0.25 & 0.13\tabularnewline
Pelvic Girdle Length & -0.80 & 0.17\tabularnewline
\bottomrule
\end{longtable}

Table D4: Contribution of the variation present in the P-matrices of the
Caribbean, SSI generalist-only, and SSI radiating groups to PC axes one
and two.\selectlanguage{english}
\begin{longtable}[]{@{}lllllllll@{}}
\toprule
Caribbean & Caribbean & Caribbean & SSI Generalist-Only & SSI
Generalist-Only & SSI Generalist-Only & SSI Radiation & SSI Radiation &
SSI Radiation\tabularnewline
Trait & PC1 & PC2 & Trait & PC1 & PC2 & Trait & PC1 & PC2\tabularnewline
Dentigerous Arm Depth & 1.16 & 0.216 & Jaw Opening In Lever & 7.63 &
1.78 & Maxillary Head Protrusion & 7.57 & 28.5\tabularnewline
Jaw Closing In lever & 0.773 & 0.186 & Ectopterygoid & 4.17 & 1.19 &
Lower Jaw Length & 0.759 & 15.2\tabularnewline
Lower Jaw Length & 0.839 & 0.177 & Jaw Closing In lever & 6.23 & 1.03 &
Dentigerous Arm Width & 0.219 & 10.0\tabularnewline
Palatine Height & 0.770 & 0.145 & Suspensorium Length & 2.56 & 0.951 &
Dentigerous Arm Base & 0.172 & 6.47\tabularnewline
Jaw Opening In Lever & 1.42 & 0.145 & Palatine Height & 3.35 & 0.706 &
Dentigerous Arm Depth & 0.203 & 5.17\tabularnewline
Dentigerous Arm Base & 0.620 & 0.138 & Maxilla Length & 3.22 & 0.671 &
Palatine Height & 0.0211 & 4.72\tabularnewline
Maxilla Length & 0.835 & 0.131 & Maxillary Head Protrusion & 1.47 &
0.657 & Maxilla Length & 0.00304 & 4.17\tabularnewline
Pelvic Girdle Length & 1.78 & 0.112 & Dentigerous Arm Base & 3.04 &
0.608 & Jaw Opening In Lever & 1.69 & 3.40\tabularnewline
Suspensorium Length & 0.681 & 0.106 & Head Depth & 1.62 & 0.595 & Jaw
Closing In lever & 0.0232 & 2.75\tabularnewline
Orbit Diameter & 0.931 & 0.0824 & Maxillary Head Height & 7.01 & 0.574 &
Nasal Tissue Protrusion & 9.46 & 1.21\tabularnewline
Dentigerous Arm Width & 0.306 & 0.0706 & Ascending Process Length & 2.40
& 0.568 & Ascending Process Length & 0.0592 & 0.915\tabularnewline
Nasal Tissue Protrusion & 2.85 & 0.0664 & Dentigerous Arm Width & 3.49 &
0.535 & Suspensorium Length & 0.569 & 0.882\tabularnewline
Head Depth & 1.34 & 0.0587 & Dentigerous Arm Depth & 2.60 & 0.447 &
Cranial Height & 0.423 & 0.510\tabularnewline
Maxillary Head Height & 4.18 & 0.0459 & Lower Jaw Length & 1.88 & 0.377
& Ectopterygoid & 0.149 & 0.442\tabularnewline
Ectopterygoid & 0.111 & 0.0203 & Cranial Height & 0.192 & 0.330 & Head
Depth & 0.691 & 0.368\tabularnewline
Cranial Height & 1.49 & 0.0104 & Pelvic Girdle Length & 1.02 & 0.315 &
Pelvic Girdle Length & 0.375 & 0.0599\tabularnewline
Ascending Process Length & 0.375 & 0.00276 & Orbit Diameter & 0.435 &
0.198 & Orbit Diameter & 0.513 & 0.00618\tabularnewline
Maxillary Head Protrusion & 1.11 & 0.0000133 & Nasal Tissue Protrusion &
0.147 & 0.0893 & Maxillary Head Height & 1.22 & 0.00372\tabularnewline
Total Variation & 21.6\% & 1.71\% & & 52.5\% & 11.6\% & & 24.1 &
84.8\%\tabularnewline
\bottomrule
\end{longtable}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image12/image12}

\end{center}
\end{figure}

Figure D1: Results of Bayesian mixed model displaying A) differences in
average trait value for generalists, snail-eaters, and scale-eaters, b)
average trait value and standard deviations for each population, and c)
the estimated standard deviation due to the random effect of population
on each of the traits. In panel A, \textbf{v} iolin plots display the
posterior draws from a Bayesian mixed effects model with the 18
craniofacial trait values as the response variable, pupfish group and
species designations as fixed effects and population ID as a random
effect. Black points and lines represent the median and 95\% HDI for
each species. For panel B, points, represent mean value for each
population while the colored lines show the 95\% confidence interval
around these estimates. For panel C, points show the average standard
deviation estimate for each trait while the black lines show the 95\%
confidence intervals around these estimates.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image13/image13}

\end{center}
\end{figure}

Figure D2. Visualization of the loadings for PC axes one and two from a
principal components analysis calculated with the variance and
covariance estimates for 18 craniofacial traits for the Caribbean, SSI
generalist-only, and SSI radiating populations. The colors of each
loading represent their overall contribution (in percentages) to both
axes.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image14/image14}

\end{center}
\end{figure}

Figure D3. Visualization of the R\textsuperscript{2} values from linear
models investigating the relationship between the covariances of
dentigerous arm width, lower jaw length, maxilla length, and palatine
height and the traits labeled in each graph for Caribbean, SSI
generalist-only, and SSI radiating groups. The first column of graphs
shows a relationship between traits that shift across groups in a manner
that suggests that they promote diversification. The second column shows
relationships between traits which may inhibit diversification
specifically in SSI generalist-only groups, as they deviate from those
seen in the Caribbean and SSI radiation.

\emph{Appendix E: Estimates for F2 Hybrid Analyses}

Table E1: Trait values for F0 snail-eater and scale-eater parental types
for the Crescent Pond and Little Lake crosses.\selectlanguage{english}
\begin{longtable}[]{@{}lllll@{}}
\toprule
\begin{minipage}[t]{0.19\columnwidth}\raggedright\strut
Trait\strut
\end{minipage} & \begin{minipage}[t]{0.19\columnwidth}\raggedright\strut
Crescent Pond Scale-eater\strut
\end{minipage} & \begin{minipage}[t]{0.19\columnwidth}\raggedright\strut
Crescent Pond Snail-eater\strut
\end{minipage} & \begin{minipage}[t]{0.19\columnwidth}\raggedright\strut
Little Lake Scale-eater\strut
\end{minipage} & \begin{minipage}[t]{0.19\columnwidth}\raggedright\strut
Little Lake Snail-eater\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
Dentigerous Arm Width\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.36\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.0024\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.18\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.19\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
Dentigerous Arm Base\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.33\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.19\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.098\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.081\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
Lower Jaw Length\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.41\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.19\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.19\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.21\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
Dentigerous Arm Depth\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.11\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.11\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.13\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.19\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
Palatine Height\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.26\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.058\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.099\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.097\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
Jaw Closing In Lever\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.14\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.0081\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.10\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.099\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
Cranial Height\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.047\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.12\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.013\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.012\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
Orbit Diameter\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.055\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.051\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.017\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.055\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
Suspensorium Length\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.12\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.047\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.013\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.0054\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
Ascending Process Length\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.29\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.20\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.046\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.012\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
Jaw Opening In Lever\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.076\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.19\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.042\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.13\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
Nasal Tissue Protrusion\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.052\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.36\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.16\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.0080\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
Ectopterygoid\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.015\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.016\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.0017\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.023\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
Head Depth\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.11\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.066\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.0083\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.019\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
Maxillary Head Protrusion\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.16\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.42\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.21\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.36\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
Maxilla Length\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.27\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.043\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.11\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.095\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
Maxillary Head Height\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.098\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.36\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.013\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.021\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
Pelvic Girdle Length\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.043\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.0032\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
0.052\strut
\end{minipage} & \begin{minipage}[t]{0.17\columnwidth}\raggedright\strut
-0.11\strut
\end{minipage}\tabularnewline
\bottomrule
\end{longtable}

Table E2: Average trait value estimates for Crescent Pond and Little
Lake F2 hybrids under the assumptions of additivity. Estimates are the
mid parent values calculated as (Trait Value \textsubscript{Parent 1} +
Trait Value \textsubscript{Parent 2}) / 2.\selectlanguage{english}
\begin{longtable}[]{@{}lll@{}}
\toprule
Trait & Crescent Pond F2 Additive Expectations & Little Lake F2 Additive
Expectations\tabularnewline
Dentigerous Arm Width & 0.18 & -0.0038\tabularnewline
Dentigerous Arm Base & 0.073 & 0.0081\tabularnewline
Lower Jaw Length & 0.11 & -0.011\tabularnewline
Dentigerous Arm Depth & 0.00054 & -0.034\tabularnewline
Palatine Height & 0.10 & 0.0012\tabularnewline
Jaw Closing In Lever & 0.065 & 0.00044\tabularnewline
Cranial Height & -0.083 & -0.012\tabularnewline
Orbit Diameter & -0.053 & -0.019\tabularnewline
Suspensorium Length & 0.083 & 0.0091\tabularnewline
Ascending Process Length & 0.045 & 0.017\tabularnewline
Jaw Opening In Lever & 0.055 & 0.042\tabularnewline
Nasal Tissue Protrusion & 0.20 & -0.075\tabularnewline
Ectopterygoid & 0.016 & 0.013\tabularnewline
Head Depth & -0.090 & -0.014\tabularnewline
Maxillary Head Protrusion & 0.13 & 0.076\tabularnewline
Maxilla Length & 0.11 & 0.0057\tabularnewline
Maxillary Head Height & 0.13 & 0.0037\tabularnewline
Pelvic Girdle Length & -0.023 & -0.029\tabularnewline
\bottomrule
\end{longtable}

Table E3: Estimates of expected variation for Crescent Pond and Little
Lake F2 hybrids if crosses follow the laws of independent assortment and
segregation. We calculated expected variation using the parental trait
values, expected average F2 trait values, and sample sizes.\selectlanguage{english}
\begin{longtable}[]{@{}lll@{}}
\toprule
\begin{minipage}[t]{0.32\columnwidth}\raggedright\strut
Trait\strut
\end{minipage} & \begin{minipage}[t]{0.32\columnwidth}\raggedright\strut
Crescent Pond F2 Trait Variation Expectations\strut
\end{minipage} & \begin{minipage}[t]{0.32\columnwidth}\raggedright\strut
Little Lake F2 Trait Variation Expectations\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
Dentigerous Arm Width\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.016\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.017\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
Dentigerous Arm Base\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.034\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.0040\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
Lower Jaw Length\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.045\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.019\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
Dentigerous Arm Depth\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.0057\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.013\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
Palatine Height\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.013\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.0048\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
Jaw Closing In Lever\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.0027\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.0049\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
Cranial Height\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.00066\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.00000021\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
Orbit Diameter\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.0000017\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.00063\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
Suspensorium Length\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.00063\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.0000068\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
Ascending Process Length\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.031\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.00043\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
Jaw Opening In Lever\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.0086\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.0035\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
Nasal Tissue Protrusion\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.012\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.0034\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
Ectopterygoid\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.000000056\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.000058\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
Head Depth\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.00029\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.000014\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
Maxillary Head Protrusion\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.043\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.041\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
Maxilla Length\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.012\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.0050\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
Maxillary Head Height\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.026\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.00014\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
Pelvic Girdle Length\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.00020\strut
\end{minipage} & \begin{minipage}[t]{0.30\columnwidth}\raggedright\strut
0.0032\strut
\end{minipage}\tabularnewline
\bottomrule
\end{longtable}

Table E4: Estimated marginal mean trait values and 95\% HDI for F2
hybrids from a Bayesian model with the 18 craniofacial trait values as
the response variable and population ID as the fixed effect.\selectlanguage{english}
\begin{longtable}[]{@{}lll@{}}
\toprule
Trait & Crescent Pond & Little Lake\tabularnewline
Dentigerous Arm Width & -8e-04 (-0.01, 0.0082) & -0.00075 (-0.013,
0.011)\tabularnewline
Dentigerous Arm Base & -0.00063 (-0.01, 0.0089) & -0.0069 (-0.021,
0.0065)\tabularnewline
Lower Jaw Length & -0.0015 (-0.0092, 0.0069) & -0.0029 (-0.013,
0.0074)\tabularnewline
Dentigerous Arm Depth & 0.0013 (-0.0081, 0.01) & -0.002 (-0.014,
0.01)\tabularnewline
Palatine Height & 0.00066 (-0.008, 0.0091) & 8e-04 (-0.011,
0.014)\tabularnewline
Jaw Closing In Lever & 0.0018 (-0.0085, 0.012) & 0.00069 (-0.015,
0.016)\tabularnewline
Cranial Height & -0.003 (-0.013, 0.0067) & -0.0046 (-0.018,
0.0095)\tabularnewline
Orbit Diameter & -0.0028 (-0.01, 0.0046) & 0.011 (-5e-04,
0.024)\tabularnewline
Suspensorium Length & -0.0021 (-0.0098, 0.006) & 0.0034 (-0.0071,
0.015)\tabularnewline
Ascending Process Length & -0.0026 (-0.022, 0.016) & 0.019 (-0.0075,
0.043)\tabularnewline
Jaw Opening In Lever & -0.0095 (-0.028, 0.0069) & -0.011 (-0.037,
0.014)\tabularnewline
Nasal Tissue Protrusion & 0.0075 (-0.023, 0.036) & -0.0038 (-0.04,
0.034)\tabularnewline
Ectopterygoid & -0.0035 (-0.019, 0.012) & 0.0024 (-0.021,
0.023)\tabularnewline
Head Depth & -0.0025 (-0.009, 0.0037) & -0.0046 (-0.013,
0.0037)\tabularnewline
Maxillary Head Protrusion & -0.0025 (-0.019, 0.015) & -0.0023 (-0.026,
0.021)\tabularnewline
Maxilla Length & -4.1e-05 (-0.0083, 0.0088) & -0.00045 (-0.012,
0.012)\tabularnewline
Maxillary Head Height & 0.0021 (-0.028, 0.031) & 0.01 (-0.033,
0.047)\tabularnewline
Pelvic Girdle Length & 0.00067 (-0.0093, 0.011) & -0.00048 (-0.014,
0.012)\tabularnewline
\bottomrule
\end{longtable}

Table E5: Correlation of loadings with PC axes one and two from a
principal components analysis calculated with the variance and
covariance estimates for 18 craniofacial traits for F2 hybrids from
Crescent Pond and Little Lake.\selectlanguage{english}
\begin{longtable}[]{@{}lll@{}}
\toprule
Trait & PC1 & PC2\tabularnewline
Suspensorium Length & 0.75 & -0.25\tabularnewline
Palatine Height & 0.75 & -0.26\tabularnewline
Ascending Process Length & 0.71 & -0.21\tabularnewline
Maxilla Length & 0.68 & -0.48\tabularnewline
Jaw Closing In Lever & 0.53 & -0.30\tabularnewline
Maxillary Head Protrusion & 0.35 & -0.56\tabularnewline
Cranial Height & 0.35 & -0.76\tabularnewline
Ectopterygoid & 0.10 & 0.45\tabularnewline
Orbit Diameter & -0.058 & -0.55\tabularnewline
Head Depth & -0.077 & -0.85\tabularnewline
Nasal Tissue Protrusion & -0.49 & 0.048\tabularnewline
Jaw Opening In Lever & -0.57 & 0.22\tabularnewline
Dentigerous Arm Width & -0.59 & -0.77\tabularnewline
Maxillary Head Height & -0.61 & -0.44\tabularnewline
Dentigerous Arm Base & -0.69 & -0.46\tabularnewline
Lower Jaw Length & -0.72 & -0.41\tabularnewline
Dentigerous Arm Depth & -0.74 & -0.29\tabularnewline
Pelvic Girdle Length & -0.79 & 0.089\tabularnewline
\bottomrule
\end{longtable}

Table E6: Contribution of the variation present in the P-matrices of F2
hybrids from Crescent Pond and Little Lake to PC axes one and two.\selectlanguage{english}
\begin{longtable}[]{@{}llllll@{}}
\toprule
Crescent Pond & & & Little Lake & Little Lake & Little
Lake\tabularnewline
Trait & PC1 & PC2 & Trait & PC1 & PC2\tabularnewline
Ascending Process Length & 1.47 & 2.26 & Maxillary Head Height & 23.6 &
12.6\tabularnewline
Suspensorium Length & 0.527 & 3.55 & Ascending Process Length & 15.0 &
1.93\tabularnewline
Cranial Height & 0.453 & 0.559 & Nasal Tissue Protrusion & 13.4 &
0.481\tabularnewline
Maxillary Head Protrusion & 0.375 & 0.524 & Jaw Opening In Lever & 9.58
& 3.23\tabularnewline
Maxillary Head Height & 0.308 & 3.47 & Suspensorium Length & 4.15 &
0.100\tabularnewline
Jaw Closing In lever & 0.290 & 1.94 & Pelvic Girdle Length & 4.09 &
1.60\tabularnewline
Head Depth & 0.260 & 1.38 & Jaw Closing In lever & 3.99 &
0.957\tabularnewline
Palatine Height & 0.194 & 0.586 & Palatine Height & 3.85 &
4.49\tabularnewline
Maxilla Length & 0.142 & 0.405 & Maxillary Head Protrusion & 3.21 &
9.47\tabularnewline
Dentigerous Arm Depth & 0.0765 & 1.45 & Dentigerous Arm Base & 2.54 &
2.23\tabularnewline
Pelvic Girdle Length & 0.0650 & 2.77 & Cranial Height & 2.52 &
5.11\tabularnewline
Lower Jaw Length & 0.0438 & 0.799 & Maxilla Length & 2.17 &
5.44\tabularnewline
Orbit Diameter & 0.0328 & 3.25 & Dentigerous Arm Depth & 2.03 &
0.871\tabularnewline
Dentigerous Arm Width & 0.0104 & 0.0856 & Dentigerous Arm Width & 1.27 &
6.00\tabularnewline
Ectopterygoid & 0.00823 & 2.98 & Lower Jaw Length & 1.24 &
0.352\tabularnewline
Jaw Opening In lever & 0.00341 & 3.76 & Ectopterygoid & 0.140 &
4.35\tabularnewline
Nasal Tissue Protrusion & 0.00122 & 3.98 & Head Depth & 0.130 &
1.77\tabularnewline
Dentigerous Arm Base & 0.00116 & 0.471 & Orbit Diameter & 0.120 &
2.04\tabularnewline
Total Variation & 7\% & 37\% & & 93\% & 63\%\tabularnewline
\bottomrule
\end{longtable}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image15/image15}

\end{center}
\end{figure}

Figure E1: Visualization of the loadings for PC axes one and two from a
principal components analysis calculated with the variance and
covariance estimates for 18 craniofacial traits for F2 hybrids from
Crescent Pond and Little Lake. The colors of each loading represent
their overall contribution (in percentages) to both axes.

\selectlanguage{english}
\FloatBarrier
\end{document}

