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.