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í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 Howea palms,
Andean Coeligena hummingbirds, Amphilophus sp. cichlids,
and 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 CaribbeanAnolis 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 ofGeospiza 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 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. 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 (C. brontotheroides )
and scale-eating (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.