A Novel Body Plan Alters Diversification of Body Shape and Genitalia in Live-Bearing Fish

Major evolutionary innovations can greatly influence subsequent evolution. While many major transitions occurred in the deep past, male live-bearing fishes (family Poeciliidae) more recently evolved a novel body plan. This group possesses a three-region axial skeleton, with one region—the ano-urogenital region—representing a unique body region accommodating male genitalic structures (gonopodial complex). Here we evaluate several hypotheses for the evolution of diversity in this region and examine its role in the evolution of male body shape. Examining Gambusia fishes, we tested a priori predictions for (1) joint influence of gonopodial-complex traits on mating performance, (2) correlated evolution of gonopodial-complex traits at macro- and microevolutionary scales, and (3) predator-driven evolution of gonopodial-complex traits in a post-Pleistocene radiation of Bahamas mosquitofish. We found the length of the sperm-transfer organ (gonopodium) and its placement along the body (gonopodial anterior transposition) jointly influenced mating success, with correlational selection favoring particular trait combinations. Despite these two traits functionally interacting during mating, we found no evidence for their correlated evolution at macro- or microevolutionary scales. In contrast, we did uncover correlated evolution of modified vertebral hemal spines (part of the novel body region) and gonopodial anterior transposition at both evolutionary scales, matching predictions of developmental connections between these components. Developmental linkages in the ano-urogenital region apparently play key roles in evolutionary trajectories, but multiple selective agents likely act on gonopodium length and cause less predictable evolution. Within Bahamas mosquitofish, evolution of hemal-spine morphology, and gonopodial anterior transposition across predation regimes was quite predictable, with populations evolving under high predation risk showing more modified hemal spines with greater modifications and a more anteriorly positioned gonopodium. These changes in the ano-urogenital vertebral region have facilitated adaptive divergence in swimming abilities and body shape between predation regimes. Gonopodium surface area, but not length, evolved as predicted in Bahamas mosquitofish, consistent with a previously suggested tradeoff between natural and sexual selection on gonopodium size. These results provide insight into how restructured body plans offer novel evolutionary solutions. Here, a novel body region—originally evolved to aid sperm transfer—was apparently co-opted to alter whole-organism performance, facilitating phenotypic diversification.

Major evolutionary innovations can greatly influence subsequent evolution. While many major transitions occurred in the deep past, male live-bearing fishes (family Poeciliidae) more recently evolved a novel body plan. This group possesses a three-region axial skeleton, with one region-the ano-urogenital region-representing a unique body region accommodating male genitalic structures (gonopodial complex). Here we evaluate several hypotheses for the evolution of diversity in this region and examine its role in the evolution of male body shape. Examining Gambusia fishes, we tested a priori predictions for (1) joint influence of gonopodial-complex traits on mating performance, (2) correlated evolution of gonopodial-complex traits at macro-and microevolutionary scales, and (3) predator-driven evolution of gonopodial-complex traits in a post-Pleistocene radiation of Bahamas mosquitofish. We found the length of the sperm-transfer organ (gonopodium) and its placement along the body (gonopodial anterior transposition) jointly influenced mating success, with correlational selection favoring particular trait combinations. Despite these two traits functionally interacting during mating, we found no evidence for their correlated evolution at macro-or microevolutionary scales. In contrast, we did uncover correlated evolution of modified vertebral hemal spines (part of the novel body region) and gonopodial anterior transposition at both evolutionary scales, matching predictions of developmental connections between these components. Developmental linkages in the ano-urogenital region apparently play key roles in evolutionary trajectories, but multiple selective agents likely act on gonopodium length and cause less predictable evolution. Within Bahamas mosquitofish, evolution of hemal-spine morphology, and gonopodial anterior transposition across predation regimes was quite predictable, with populations evolving under high predation risk showing more modified hemal spines with greater modifications and a more anteriorly positioned gonopodium. These changes in the ano-urogenital vertebral region have facilitated adaptive divergence in swimming abilities and body shape between predation regimes. Gonopodium surface area, but not length, evolved as predicted in Bahamas mosquitofish, consistent with a previously suggested tradeoff between natural and sexual selection on gonopodium size. These
Major transitions in body plan clearly set the stage for much subsequent morphological evolution. But since the Cambrian explosion, very few novel body plans have appeared (Gould, 1989;Levinton, 1992;Erwin et al., 1997;Davidson and Erwin, 2006;Maxwell and Wilson, 2013). For fish, primary literature and textbooks in ichthyology and comparative vertebrate anatomy describe the teleost axial skeleton body plan and vertebral formulae as a two-part system: anterior trunk and posterior caudal vertebral regions [see Rivera-Rivera et al. (2010) and references therein]. However, recent work revealed male livebearing fishes in the family Poeciliidae have altered this ancestral body plan with a novel three-part body plan composed of the two aforementioned regions plus a third region known as the ano-urogenital vertebral region (Rosa-Molinar et al., 1994, 1996Rivera-Rivera et al., 2010). This body plan appears to have evolved at least 44 million years ago (Hrbek et al., 2007). Because of the relatively little research conducted on this novel body plan to date, and the emphasis on Poeciliidae so far, its origination and extent is not fully understood-this body plan could be more widespread throughout much of Atherinomorpha.
While many organisms have a general body region encompassing the anus, urinary tract, and reproductive organs, the "ano-urogenital vertebral region" is unique to male poeciliid fishes. This region reflects a marked remodeling of the axial and appendicular skeleton to accommodate development of male genitalia and associated structures. The novel region includes a highly modified anal fin used for sperm transfer (the gonopodium) along with its musculoskeletal components. This unique region comprises vertebrae with ribs and hemal spines (Rosa-Molinar et al., 1994 and is thought to reflect an adaptation for rapid and effective sperm transfer and internal fertilization in the group Serrano-Velez et al., 2014). The ano-urogenital region is situated between the anterior trunk and posterior caudal regions, from vertebrae 11 to 16. Development of this region most notably involves the elongation, thickening, and anterior bending of hemal spines (typically of vertebrae 14-16; sometimes termed "gonapophyses") as well as the anterior transposition of the gonopodium and the gonopodial appendicular support ( Figure 1A). As a consequence, poeciliid fishes exhibit striking sexual dimorphism in adults, with the male anal fin (i.e., gonopodium) positioned more anteriorly than the female anal fin (Supplementary Figure 1).
Because the evolution of novel body plans can have important consequences for further morphological and behavioral evolution, it is perhaps not surprising that poeciliid fishes exhibit great variability in the various features of the gonopodial complex, such as the number of modified hemal spines (typically between 1 and 3), their size and extent of modification (e.g., thickening, bending, uncinate processes), the degree of gonopodial anterior transposition (i.e., relative placement of the gonopodium along the body), and the size of the gonopodium (Rosen and Gordon, 1953;Rosen and Bailey, 1963;Chambers, 1987;Ghedotti, 2000;Langerhans, 2011). Yet, to date, explanations for the origins of such diversity and the potential consequences for the evolution of other body regions are lacking (Rosa-Molinar, 2005;Martinez-Rivera et al., 2010;Langerhans, 2011). When novel innovations evolve, can we predict their subsequent evolution based on our functional understanding of the structures, and how might novel traits influence adaptive diversification after their initial emergence?
Here we address several outstanding questions regarding the evolution of this understudied, unique region and its role in facilitating adaptive diversification in poeciliid fishes. We directly tackled three broad questions: (Q1) How are various components of the gonopodial complex developmentally or functionally integrated, and how might they evolve in a correlated manner? (Q2) What forms of selection shape diversity of this region? and (Q3) Has the novel ano-urogenital region been co-opted for adaptive diversification in lateral body shape? To address Q1 (blue region in Figure 1), we used a laboratory mating experiment to test whether correlational selection favors particular combinations of two gonopodialcomplex traits ( Figure 1B, see below) and used comparative analyses to test predictions of correlated evolution of particular traits in the ano-urogenital region due to developmental linkages and selection on copulatory performance ( Figure 1C, see below). We addressed Q2 (yellow region in Figure 1) by testing among predatory environments specific predictions of how divergent selection on body shape and gonopodium size should drive evolutionary divergence of gonopodial-complex traits ( Figure 1D, see below). Finally, we addressed Q3 (green region in Figure 1) by directly testing the co-optation prediction of correlated evolution between gonopodial anterior transposition and lateral body shape during predator-driven diversification in an adaptive radiation ( Figure 1E, see below).
FIGURE 1 | Illustration of (A) the gonopodial complex in live-bearing fishes (Gambusia affinis depicted), and hypothetical depictions of our (B) hypothesis of correlational selection on two gonopodial-complex traits, (C) predictions of correlated evolution among traits, (D) predictions of trait divergence between predation (Continued) Frontiers in Ecology and Evolution | www.frontiersin.org 3 February 2021 | Volume 9 | Article 619232 FIGURE 1 | regimes, and (E) prediction of co-optation of the ano-urogenital region for the evolution of different body forms (see text for details). Labels in (A) indicate key components of the ano-urogenital region: vertebrae 11 and 16 (V11, V16; demarcating the ano-urogenital vertebral region), the gonopodium (G), and its appendicular skeletal support [hemal spines (HS), the uncinate process on the 15 th hemal spine (UNC), the columnar fusion of pterygiophores 2-5 (P2-5)]. Anterior/rostral is toward the left in (A). Fitness in (B) ranges from low (cool colors) to high (warm colors). The three colored panels reflect the three core questions addressed in the study (Q1-Q3, see text).

Overview of Hypotheses and Predictions
All hypotheses evaluated in this study derive from previous research, with testable evolutionary predictions being derived from the hypotheses. We first tested a copulatory-performance hypothesis of functional integration between gonopodial anterior transposition and gonopodium length. We specifically hypothesized that when the gonopodium is located closer to the head (greater gonopodial anterior transposition), high performance during copulation should be achieved by a shorter gonopodium, and vice versa. This hypothesis describes positive correlational selection on the two traits, with a ridge-shaped selection surface where trait combinations along the top of the ridge have similarly high mating performance values ( Figure 1B). Copulation in most live-bearing fishes is very rapid, involving complex locomotor maneuvers, and selection favoring efficient sperm transfer may be strong in many species (Rosa-Molinar, 2005;Martinez-Rivera et al., 2010;Rivera-Rivera et al., 2010;Langerhans, 2011;Devigili et al., 2015;Head et al., 2015). As revealed in high-speed video (1000 Hz) of western mosquitofish (Gambusia affinis), the male circumducts the gonopodium, with the distal tip pointing anteriorly, and performs a rapid, complex torque-thrust behavior in an attempt to inseminate the female during copulation (Rosa-Molinar, 2005;Martinez-Rivera et al., 2010;Rivera-Rivera et al., 2010). When the gonopodium is positioned more posteriorly along the body (less transposition), a male may exhibit reduced maneuverability and acceleration [owing to a reduced caudal region (see below) and more posterior center of mass] and require a longer gonopodium to enhance insemination or fertilization success during copulation. Additionally, if males use visual cues for gonopodial placement during copulation (Langerhans, 2011), then an optimal position of the gonopodial tip relative to the eye would reinforce this functional integration. Thus, males' mating success should depend on the combination of the two traits, and we tested this prediction of correlational selection using a mating experiment. We then examined two predictions regarding correlated evolution of genital traits among species ( Figure 1C). First, we have an anatomical hypothesis that the magnitude of gonopodial anterior transposition is largely governed postdevelopmentally by variation in the extensive sexually dimorphic ligament system associated with the modified hemal spines (Rosa-Molinar et al., 1994). This anatomical hypothesis leads to the evolutionary prediction that species having a greater number of modified hemal spines will additionally exhibit greater gonopodial anterior transposition (gonopodium closer to the head, resulting in a smaller abdominal region). That is, an increased number of modified hemal spines should essentially push the gonopodium anteriorly to a greater extent through direct movement of the columnar fusion of pterygiophores 2-5 (see Figure 1A). Second, if our copulatory-performance hypothesis described above is accurate and captures the major aspects of selection on the traits (i.e., other aspects of natural and sexual selection are comparatively of little importance), this leads to an evolutionary prediction: correlational selection on gonopodial anterior transposition and gonopodium length will cause the two traits to exhibit correlated evolution among species/populations. Using phylogenetic comparative analysis of 10 closely related species in the genus Gambusia (mosquitofishes) known to exhibit variation in these genital traits, we tested the predictions that (1) a greater number of modified hemal spines is associated with greater gonopodial anterior transposition and that (2) a more posteriorly positioned gonopodium is associated with a longer gonopodium length ( Figure 1C).
Next, we tested three general predictions regarding evolutionary divergence of genital traits between predatory environments ( Figure 1D). Previous work has demonstrated the importance of two particular factors in generating diversity of gonopodium size and general body form in live-bearing fishes: (1) sexual selection via female mating preferences and (2) natural selection via predation from piscivorous fish [both reviewed in Langerhans (2010), Langerhans (2011)]. Specifically, females of some poeciliid species exhibit mating preferences for males with particular body shapes and larger gonopodia (Langerhans et al., 2005(Langerhans et al., , 2007Kahn et al., 2010;Langerhans and Makowicz, 2013), while natural selection in the presence of fish predators favors larger caudal regions and smaller gonopodia due to their positive effects on locomotor escape behaviors (Langerhans et al., 2005;Langerhans, 2009aLangerhans, , 2010Langerhans, , 2011. This has led to divergence in these traits between environments with and without piscivorous fish. For example, male poeciliids of a number of species exhibit smaller gonopodia and larger caudal regions in localities with higher intensities of fish predation [reviewed in Langerhans (2010), Langerhans (2011), Langerhans (2018]. To date, no study has investigated the role the gonopodial complex might play in these observed patterns of morphological divergence. Divergent selection between environments with and without predatory fish arising from this interplay of natural and sexual selection leads to three general predictions of phenotypic divergence in the gonopodial complex: males in high-predation environments should exhibit (1) smaller gonopodia, (2) more modified hemal spines that are larger and more strongly modified, and (3) a greater magnitude of gonopodial anterior transposition compared to males in low-predation localities ( Figure 1D). The first prediction derives from previous work demonstrating both female preference for larger gonopodia and detrimental effects of gonopodium size on escape performance and survivorship with a predatory fish. The latter two predictions derive from the functional and developmental hypotheses that increasing gonopodial anterior transposition should facilitate both shorter gonopodia and larger caudal regions. If selection on general body shape has driven evolution in the anourogenital region, we should see correlated divergence of gonopodial anterior transposition and lateral body shape across populations ( Figure 1E). Here we test these predictions using the model system of the post-Pleistocene radiation of Bahamas mosquitofish (Gambusia hubbsi) inhabiting inland blue holes.

Experimental Test of Functional Integration of Gonopodial-Complex Traits
Using laboratory-born Bahamas mosquitofish from a single population on Andros Island, The Bahamas (Cousteau's), we conducted 39 mating trials to test for correlational selection via copulatory performance on two male traits: gonopodial anterior transposition and gonopodium length. All fish were raised under common laboratory conditions in 10-L aquaria within a recirculating system (providing biological, mechanical and UV-filtration) at ∼25 • C in a temperature-controlled room and fed a varied diet of live brine shrimp, freeze-dried daphnia and bloodworms, and TetraMin Pro flakes. Each trial comprised one virgin female placed with one virgin male in a 2.5-L polycarbonate aquarium (22 × 9.5 × 12 cm, L × W × H) with no structure for 1 h. The first 30 min were video recorded from the side with a Sony DCR-SR68 camera (Sony, Tokyo, Japan). From the video, we counted the number of copulation attempts and the number of apparently successful copulations where appropriate genital contact seemed to occur. We watched videos frame-by-frame during copulation attempts, and recorded genital contact when unambiguous physical contact was made between the male gonopodium tip and the female body in the immediate vicinity of the urogenital opening. Fish were euthanized immediately following the trial, photographed for morphological measurements (standard length and the two aforementioned traits), and preserved in 95% ethanol. Using tpsDig2 (Rohlf, 2017), we measured gonopodial anterior transposition as the distance from the anterodorsal tip of anal fin ray 1 to the center of the eye orbit, and gonopodium length as the distance from the anterodorsal tip of anal fin ray 1 to the distal tip of the gonopodium. We estimated mating success in three ways, as it was not obvious which measure might be most appropriate: genital contact success (0 or 1), number of genital contacts, and genital contact efficiency (number of genital contacts divided by number of copulation attempts). Using sperm retrieval, we confirmed that these estimates provided adequate surrogates for insemination success (see Supplementary Material).
To test for positive correlational selection on gonopodium length and gonopodial anterior transposition, we employed three statistical approaches, one for each estimate of mating success. Prior to analysis, we regressed each log 10 -transformed trait on log 10 -transformed standard length and saved residuals to examine size-independent variables. Note that results are virtually identical if we measured relative gonopodial anterior transposition as its proportion of standard length, as described below (correlation among the two trait estimates: r = 0.999, P < 0.0001). In each case, we tested for significance of the interaction between the two traits on variation in mating success, and calculated the standardized correlational selection gradient (γ) using multiple regression of standardized trait values (mean = 0, standard deviation = 1) on relative fitness (mating success score divided by average mating success score) (Lande and Arnold, 1983). To accurately and thoroughly characterize selection on these traits, we additionally tested for directional selection (linear selection acting directly on each trait) and quadratic selection (i.e., stabilizing/disruptive selection on each trait). Note that the hypothesized ridge-shaped selection surface implies both positive correlational selection on the traits, as well as stabilizing selection acting on both traits. For significance testing, we conducted a logistic regression for genital contact success, a generalized linear model with a Poisson distribution and log link for the number of genital contacts, and a multiple regression for genital contact efficiency. We arcsin square-root transformed genital contact efficiency, as this transformation of the proportional values increased normality of model residuals. In each model, gonopodium length, gonopodial anterior transposition, their squared terms, and their interaction served as independent variables. The significance of the interaction term from these models provided the relevant tests of correlational selection. To interpret any observed correlational selection, we visualized the quadratic selection surface.

Correlated Evolution of Gonopodial-Complex Traits Among Species
To test predictions of correlated evolution among species, we chose to focus on the Gambusia genus, as this is the most species-rich poeciliid genus and has experienced the greatest scrutiny to date in understanding the development and anatomy of the gonopodial complex (Langerhans, 2011). We generated phylogenetic hypotheses for 10 Gambusia species using mitochondrial and nuclear gene sequences, measured the relevant traits using radiograph images, and calculated standardized phylogenetic independent contrasts. First, we used PCR to amplify a 975-bp fragment of the NADH subunit 2 mitochondrial gene (ND2) (following Langerhans et al., 2012), a 402-bp fragment of the cytochrome b mitochondrial gene (cyt b) (following Lydeard et al., 1995), and a 760-bp fragment of the first intron of the S7 ribosomal protein nuclear gene (S7) (following Langerhans et al., 2012) for each of 10 species of Gambusia, as well as two outgroup species, Heterophallus rachovii (Gambusia rachovii) and Belonesox belizanus. We sequenced a single specimen for most species (outgroups, G. clarkhubbsi, G. gaigei, G. geiseri, G. heterochir, G. hurtadoi, G. speciosa) but included multiple samples (each from geographically distinct populations) for four species (G. affinis: 4, G. aurata: 3, G. holbrooki: 3, G. quadruncus: 4). Sequences were aligned by eye and deposited in GenBank (Supplementary Table 1). We constructed phylogenetic trees using two separate Bayesian analyses with MrBayes 3.1.2 (Ronquist and Huelsenbeck, 2003), one for concatenated mitochondrial genes and one for the nuclear gene, as the different datasets revealed significantly different evolutionary histories (partition homogeneity test, P = 0.02). For additional details, see the Supplementary Material.
We obtained radiograph images of 4-50 adult males for each species to obtain estimates of the number of modified hemal spines on vertebrae 14-16, the magnitude of gonopodial anterior transposition, and the relative length of the gonopodium. The number of modified hemal spines, defined as hemal spines exhibiting obvious anterior bending, uncinate processes, or both, was counted on each radiograph. Although most species exhibit a fixed number of modified hemal spines (Rosen and Bailey, 1963;Rauchenberger, 1989), some species are polymorphic. Of the 10 species examined in the interspecific component of this study, two species (G. clarkhubbsi, G. holbrooki) were found to exhibit variation in modified hemal spine number (either 2 or 3) and were classified as 2.5 for the analyses performed here. We measured the magnitude of gonopodial anterior transposition as the relative placement of the columnar fusion of pteryogiophores 2-5 along the body: distance from the anterodorsal tip of the second pterygiophore to the center of the eye orbit divided by standard length (larger ratio values indicate relatively farther distances from the eye, and thus less transposition). We employed this ratio metric as our estimate of relative gonopodial anterior transposition because we were specifically interested in the relative placement of the gonopodium on the body (results are qualitatively similar using residuals of log-log regression on standard length instead, as the two estimates of relative gonopodial anterior transposition were highly correlated, r = 0.99, P < 0.0001). Gonopodium length was measured as the distance from the anterodorsal tip of anal fin ray 1 to the distal tip of the gonopodium. To control for the possible confounding effects of body size, we tested for associations between the three traits described above and standard length using linear regression and statistically removed effects of size where necessary. Relative gonopodial anterior transposition was not correlated with body size, and thus no additional size correction was performed for this variable. Correlation with body size was evident for the other two traits. We conducted phylogenetic size correction (following Revell, 2009) for both modified hemal-spine number and gonopodium length to eliminate size effects, i.e., calculated residuals from linear regression of hemal spine number and log 10 -transformed gonopodium length on log 10 -transformed standard length, controlling for non-independence due to phylogenetic history.
We calculated standardized phylogenetic independent contrasts using the ape package (Paradis et al., 2004) in the R statistical program (R Core Team, 2019) for both the mtDNA and nDNA phylogenetic hypotheses (after pruning the tree so that each species was represented by a single tip). Because normality of residuals was evident in all cases, we used linear regression through the origin to test for correlated evolution among (1) the number of modified hemal spines and gonopodial anterior transposition and (2) gonopodial anterior transposition and relative gonopodium length. For analyses using the nDNA phylogeny, we adjusted degrees of freedom based on the number of polytomies (Purvis and Garland, 1993;Garland and Diaz-Uriarte, 1999).

Predator-Driven Divergence in Gonopodial-Complex Traits in Bahamas Blue Holes
To evaluate the predictability of evolution in the gonopodial complex between divergent predation regimes, we investigated morphological divergence among male Bahamas mosquitofish inhabiting blue holes on Andros Island, The Bahamas. Bahamas mosquitofish colonized these relatively isolated blue holes (waterfilled vertical caves) during the past ∼15,000 years, with some blue holes being additionally colonized during this period by a larger predatory fish (bigmouth sleeper, Gobiomorus dormitor). Because no environmental factor is known to co-vary with the presence of the fish predator (e.g., salinity, water transparency, surface diameter, depth, resource availability; Langerhans et al., 2007;Riesch et al., 2013), this system represents an ideal "natural experiment" with which to test the effects of predation on evolutionary diversification (e.g., Langerhans et al., 2007;Langerhans, 2010Langerhans, , 2018. Based on hypothesized divergent selection between predation regimes, it is predicted that Bahamas mosquitofish should evolve larger caudal regions and smaller gonopodia in high-predation blue holes [see Langerhans (2010) for details], changes that might be facilitated by phenotypic modifications of the gonopodial complex (see Figure 1D). Prior work revealed that fish inhabiting blue holes with predatory fish have indeed evolved larger caudal regions, facilitating increased locomotor escape performance and survivorship during predator encounters (Langerhans, 2009a(Langerhans, , 2010(Langerhans, , 2018. Using the males from 18 populations examined in this study (see below), we confirmed prior studies' results regarding this pattern of morphological divergence using a partially overlapping dataset (see Supplementary Figure 2). Previous work also demonstrated a pattern of smaller gonopodia in high-predation localities of Bahamas mosquitofish across a range of habitat types (as measured by the gonopodium lateral surface area; Langerhans et al., 2005), but no direct test of differences in gonopodium size has yet been performed for inland blue hole populations. Here we perform such a test. It is important to note that although the hypothesis of divergent selection on gonopodium size focuses on effects of its surface area on fitness components (i.e., mating preferences for lateral projection area, frictional drag incurred by gonopodia during locomotion), the hypothesis of correlational selection on gonopodial anterior transposition and gonopodium size centers on its length. Thus, here we examine both gonopodium length (measured as described above) and its lateral surface area (measured as area inside the gonopodium's outer boundaries, comprising anal-fin rays 1-5) in our assessment of divergence in gonopodium size among predation regimes. The two aspects of gonopodium size are not very tightly correlated with one another (r = 0.34, P < 0.001). We test whether changes in the anourogenital region might be partially responsible for mediating responses to selection on lateral body shape and gonopodium size by comparing trait values among predation regimes (see Figure 1D).
We obtained radiograph images of 272 males from 18 blue holes: nine populations that evolved in the absence of fish predators and nine that evolved in the presence of predatory fish (Supplementary Figure 3). The following nine traits were measured on each image: number of modified hemal spines, length of the 14th hemal spine, length of the 15th hemal spine, length of the uncinate process on the 15th hemal spine, length of the 16th hemal spine, angle of the 16th hemal spine, magnitude of gonopodial anterior transposition, length of the gonopodium, and lateral surface area of the gonopodium. We selected these traits due to their hypothesized effects on gonopodial anterior transposition, and thus body shape, as well as the observation that these traits exhibited considerable variation in these fish. Length of hemal spines was measured using straight-line distances between 2 and 3 landmarks on each bone, depending on the presence of bone curvature. Angle of the 16th hemal spine (degrees) was measured from the angle formed by lines connecting each end of the bone with the vertex at the most posterior edge along the length of the bone (smaller angles depict greater curvature, with stronger anterior bending of the bone). Gonopodial anterior transposition was again measured as the distance from the second pterygiophore to the center of the eye orbit divided by standard length, to capture the relative position of the gonopodium along the body; however all results were again virtually identical if we used log-log residuals of the distance regressed on standard length rather than the ratio. Other traits were measured as previously described. We made all measurements using tpsDig2 software (Rohlf, 2017). Because prior work has shown that blue hole populations are relatively isolated and that fish from similar predation regimes are not more closely related to one another than to populations in different predation regimes (Schug et al., 1998;Langerhans et al., 2007;Heinen-Kay and Langerhans, 2013;Riesch et al., 2013), populations were treated as independent replicates in statistical analyses.
To test for differences among predation regimes, we first conducted a mixed-model nested multivariate analysis of covariance (MANCOVA) using all traits as dependent variables [e.g., see Riesch et al. (2013)], predation regime and population nested within predation regime (random factor) as independent variables, and log 10 -transformed standard length as a covariate to control for effects of body size. We found no evidence for heterogeneity of slopes among predation regimes (P = 0.55). This analysis provided the overall test for differences in gonopodialcomplex traits between predation regimes. To inspect traitspecific patterns of among-population variation, we additionally conducted mixed-model nested analysis of covariance (using restricted maximum likelihood estimation) separately for each trait, mirroring the model structure of the MANCOVA described above. Body size did not differ among predation regimes (mixedmodel nested ANOVA, P = 0.95), indicating that predation regime was not confounded with variation in body size. Prior to analysis, the length and area measurements were log 10transformed to meet assumptions of normality of residuals; the count data for hemal spines, the ratio metric for gonopodial anterior transposition, and the angle of the 16th hemal spine were left untransformed. For the univariate model examining variation in the number of modified hemal spines, we employed a generalized linear mixed model with a Poisson distribution and log link.
To more directly test the role of the ano-urogenital vertebral region in divergence of body shape and gonopodium length among populations, we calculated population means of our measurements of hemal-spine morphology (6 hemal-spine variables), gonopodial anterior transposition, gonopodium length, and the morphological divergence vector (d, describing multivariate lateral body shape variation, see Supplementary Material). Mean values represent least-squares means, statistically adjusted for effects of body size. We tested our anatomical prediction that hemal-spine morphology underlies divergence in gonopodial anterior transposition using canonical correlation analysis. We tested for an association between the six hemal-spine variables and gonopodial anterior transposition, and visualized the relationship using the canonical axis derived from the statistical model. Using linear regression between multivariate body shape (d) and gonopodial anterior transposition, we then tested our co-optation prediction that gonopodial anterior transposition partially underlies this primary axis of divergence in lateral body shape among populations, reflecting a novel mechanism for adaptive diversification of body form. Finally, we used linear regression to test the functional-integration prediction that populations with greater gonopodial anterior transposition will also exhibit shorter gonopodium length to enhance copulatory performance. We used the statistical program JMP (v. 14.2, 2018, SAS Institute Inc.) for most analyses, but used SAS (v. 7.15, 2017, SAS Institute Inc.) for the mixed-model MANCOVA, and the R package lme4 (Bates et al., 2011) for the generalized linear mixed model. In all cases, one-tailed P values are employed when we have a priori directional hypotheses (Rice and Gaines, 1994).
Because variation among populations in gonopodial-complex traits could reflect either genetically-based differentiation, phenotypic plasticity, or both, we tested for a genetic basis to trait divergence using a common-garden experiment. We raised lab-born offspring (n = 158) from eight populations (four without predators, four with predators) to adulthood under common laboratory conditions in a large, recirculating aquarium system comprising 72 115-L aquaria with biological, mechanical, and ultraviolet filtration. Water was maintained at ∼25 • C, a conductivity of 2,850 µS/cm, and a pH of 8.3, with a 14-h light/10-h dark photoperiod. Fish were fed daily with a mixture of TetraPro Tropical Crisps, Fluval Bug Bites for Tropical Fish, and Hikari freeze-dried Daphnia, blood worm, and brine shrimp. We collected the parental generation (F0) from the wild as newborns (in an attempt to minimize maternal and environmental effects), transported them to the laboratory at North Carolina State University, raised them to adulthood, and mated each female with one unique male from the same population (∼19 mated pairs per population). Offspring (F1) were then raised to adulthood, with sexes separated as they matured (males raised in 4-5 separate tanks per population; 5-15 males per tank).
For F1 fish, we measured all traits as described above for wild-caught fish. We then performed statistical analyses using a dataset that included for these eight populations both wild-caught fish that had been immediately euthanized and preserved after field collection and F1 lab-born fish that had been euthanized and preserved at ∼225 days of age. We conducted MANCOVA using all traits as dependent variables and ANCOVAs separately for each trait. In all models, log 10transformed standard length served as a covariate; population (testing for genetic differences), rearing environment (lab-born or wild-caught), and their interaction served as independent variables. Of particular interest in this study is whether fish born and raised in a common laboratory environment retained trait differences between populations consistent with those observed in wild-caught individuals (population term of model), indicating evolutionary divergence in gonopodial-complex traits among populations. In all cases, we did not adjust P values for multiple tests; however all significant results would remain significant when adjusted to maintain a false discovery rate of 5% (Benjamini and Hochberg, 1995).

Experimental Test of Functional Integration of Gonopodial-Complex Traits
We found significant positive correlational selection on gonopodium length and gonopodial anterior transposition when using either genital contact success or efficiency as our estimate of mating success, but not when using the number of genital contacts as the mating-success surrogate (Supplementary Table 2). We found a saddle-shaped selection surface rather than the predicted ridge-shaped selection surface (Figure 2). This shape arose from the presence of correlational selection and absence of stabilizing selection (Supplementary Table 2). Consistent with expectations, males with a long gonopodium and a long distance between their gonopodium and eye had relatively high mating success, as did males with a short gonopodium and a short distance between their gonopodium and eye. But males with intermediate values of both traits had slightly reduced mating success, while males with mismatched trait values (long gonopodia with short distance to eye or short gonopodia with long distance to eye) had especially low mating success (Figure 2). This latter effect was particularly pronounced for the combination of short gonopodia and a far distance from the eye.

Correlated Evolution of Gonopodial-Complex Traits Among Species
Bayesian analysis of mtDNA gene sequences produced a wellsupported phylogeny for the 10 Gambusia species examined (Figure 3A), while analysis of nDNA gene sequences yielded a less resolved phylogeny, differing in topology primarily with regard to the placement of G. aurata, G. speciosa, and G. holbrooki (suggesting possible historical mtDNA introgression; Figure 3B). FIGURE 2 | Quadratic selection surface for relative gonopodium length and gonopodial anterior transposition, depicting the predicted mating success of males across the range of trait values. Traits are mean-centered and in units of standard deviation. Genital contact success depicted on the Z axis (elevation), but results are similar for genital contact efficiency (see Supplementary Table 2). Surface constructed using all selection gradients (linear, quadratic, correlational).

Predator-Driven Divergence in Gonopodial-Complex Traits in Bahamas Blue Holes
In Bahamas mosquitofish, we found considerable amongpopulation variation in male genital traits, with males from divergent predation regimes differing in mean trait values for most measured traits ( Table 1). Consistent with predictions, male G. hubbsi inhabiting blue holes with predatory fish had a greater number of modified hemal spines (Figure 4A), a longer 14th, 15th, and 16th hemal spine  Figures 4A-C), a longer uncinate process on the 15th hemal spine (Supplementary Figure 4D), a smaller angle (greater anterior bending) of the 16th hemal spine (Supplementary Figure 4E), greater gonopodial anterior transposition (Figure 4B), and a smaller lateral surface area of the gonopodium (Figure 4C) than did males from blue holes without piscivorous fish. In contrast with our prediction, populations in blue holes with predatory fish did not exhibit relatively shorter gonopodia than those in predator-free blue holes (Supplementary Figure 4F). Using population means, we found support for our anatomical and co-optation predictions, but not for our functionalintegration prediction. Hemal-spine morphology was indeed significantly associated with gonopodial anterior transposition (F (6, 11) = 4.61, one-tailed P = 0.0070): populations with a greater number of modified hemal spines, longer hemal spines, a longer uncinate process on the 15th hemal spine, and greater curvature of the 16th hemal spine exhibited a gonopodium positioned more anteriorly (Figure 5A, Supplementary Table 3). Populations exhibiting greater gonopodial anterior transposition also had relatively larger caudal regions (one-tailed P = 0.0003, Figure 5B), matching predictions. However, populations with greater gonopodial anterior transposition did not have shorter gonopodia (one-tailed P = 0.93).

(Supplementary
Results from the common-garden experiment suggested that divergence in gonopodial-complex traits among populations, at least partially, reflects genetic differentiation. Our multivariate model revealed highly significant effects for all model terms, indicating that although genetically-based differences were evident, some traits for some populations exhibited differences between wild and laboratory conditions (Supplementary Table 4). Examining each trait separately revealed that some traits showed clear, genetically-based differences among populations with little-to-no effects of a laboratory environment on their expression, while others showed both genetic differences and influences of laboratory rearing (Supplementary Table 4). For instance, we found no evidence for consistent differences between wild and lab-raised fish for the length of the 15th hemal spine, length of the uncinate process of the 15th hemal spine, or gonopodial anterior transposition, but found strong effects of laboratory rearing on gonopodium size and number of hemal spines (smaller gonopodia and more hemal spines in lab-raised fish). In most cases, we found significant evidence of population-specific environmental effects (i.e., G × E), although its strength was always weaker than overall effects of population (see Supplemental Material). Thus, environmental influences on observed trait values in the wild likely exist for some traits, but genetically-based differentiation in gonopodial-complex traits was evident.

DISCUSSION
The evolution of novel traits and body plans can dramatically alter subsequent evolutionary trajectories (see Introduction).
Here we uncovered several important insights into the evolutionary causes and consequences of an understudied novel body region, the ano-urogenital vertebral region of male livebearing fish in the family Poeciliidae. Our findings revealed that developmental linkages among components of this gonopodial complex underlie correlated evolution among populations and species, but that correlational selection for functional integration of traits has not driven correlated evolution. We further discovered that predator-associated natural and sexual selection has driven evolutionary diversification of this region: novel traits intimately related to copulation have been co-opted to facilitate evolutionary changes in locomotor performance and body form. This work supports the notion that major evolutionary innovations fundamentally alter future evolutionary patterns, and one need not look into the deep past to find such examples.
In the mating experiment, we uncovered significant evidence for functional integration of two key gonopodial-complex traits. We predicted a ridge-shaped selection surface, where a combination of stabilizing selection and correlational selection would favor a positive association between relative gonopodium length and its distance from the eye-i.e., successful sperm insemination would benefit from longer gonopodia when positioned more posteriorly and from shorter gonopodia when positioned more anteriorly. The saddle-shaped correlational selection surface we observed is only partially consistent with this expectation. Our finding of positive correlational selection indeed suggests that insemination is highest for fish with the appropriate combination of gonopodium length and gonopodial anterior transposition, and particularly low for mismatched combinations of the traits; however males with intermediate values of both traits suffered slightly reduced fitness. This dip in fitness for intermediate trait values arose because of the lack of stabilizing selection on the traits-it is unclear why this fitness valley might occur. Moreover, when using the number of genital contacts as our estimate of mating success, we found no evidence  Table 1 for statistical details). Least-squares means ± one standard error depicted.
for correlational selection. This implies that this fitness estimate captures something different than the others (i.e., genital contact success and efficiency), perhaps reflecting male behavior, rather than mating success, more so than the others. For instance, males with more genital contacts might simply attempt more copulations or have larger sperm reserves, and not necessarily reflect sperm insemination the same way the other estimates do. Indeed, the number of genital contacts observed here was highly correlated with the number of coerced copulation attempts (P < 0.0001), while genital contact efficiency was not (P = 0.52). The functional integration of gonopodial anterior transposition and gonopodium length during mating is apparently more related to the ability and efficiency of achieving mating success than the overall number of times that a successful copulation might occur during a given mating encounter.
Although a number of studies have examined effects of gonopodium length and distal-tip elements on mating success in poeciliid fishes Gasparini et al., 2011;Kwan et al., 2013;Devigili et al., 2015;Head et al., 2015Head et al., , 2017Hernandez-Jimenez and Rios-Cardenas, 2017;Chung et al., 2020), no prior study has investigated the relationship examined here. Our findings add to the evidence that complex selection surfaces may characterize insemination and fertilization success in Gambusia fishes (Devigili et al., 2015;Head et al., 2015) and suggest we still have much to learn. That said, further work in additional populations and species, with larger sample sizes, will help us understand the generality and robustness of these findings. Overall, our results indicate that gonopodium length and its position along the body are functionally related to one another, and if this copulatory function captures the primary link between these traits and fitness, then this predicts correlated evolution of the two traits across populations and species.
Examining 10 closely related Gambusia species, and investigating 18 populations of the post-Pleistocene radiation of Bahamas mosquitofish, we did not find this predicted correlated evolution-we found no evidence that gonopodium length and gonopodial anterior transposition coevolve in the expected manner. The evolution of a more anteriorly positioned gonopodium apparently does not necessarily lead to the evolution of a shorter gonopodium, even though these trait combinations seem to enhance mating success. One possible explanation for this finding is that correlational selection might vary over time or among populations/species (e.g., Siepielski et al., 2009Siepielski et al., , 2013Weese et al., 2010). Further work is needed to evaluate this putative cause. However, we suggest a likely explanation for this finding is these traits experience complex forms of selection from multiple agents other than their co-dependency during copulation and their genetic basis probably involves different genes (i.e., little-to-no pleiotropy). One possible complicating factor is selection on testes size (e.g., via sperm competition), which could result in reduced gonopodial anterior transposition to accommodate larger testes. However, we have found no evidence supporting this notion (see Supplementary Material). Still, considering the disparate forms of selection, such as natural selection on locomotor performance and sexual selection via female preference, that might act on gonopodium length and gonopodial anterior transposition, it is perhaps not surprising that their evolutionary trajectories are not tightly correlated. On the other hand, at both micro-and macroevolutionary scales, we revealed correlated evolution of hemal-spine morphology and the magnitude of gonopodial anterior transposition, matching our a priori anatomical prediction. As Gambusia species evolve more modified hemal spines, they additionally evolve a positioning of the gonopodium closer to the head. Similarly, populations of Bahamas mosquitofish with more numerous, larger, and more curved modified hemal spines also show a more anteriorly positioned gonopodium. The anatomical/developmental linkage between growth of the modified hemal spines and the anterior transposition of the gonopodium has so far been described in only one species (G. affinis; Rosa-Molinar et al., 1994, but the findings here suggest this may represent an evolutionarily conserved mechanism of gonopodial transposition in poeciliid fishes. Thus, one of the most obvious external differences in anatomy between the sexes in the family Poeciliidae may be largely explained by the modified hemal spines (and their associated ligament system) present in males.
In the post-Pleistocene radiation of Bahamas mosquitofish inhabiting blue holes, we uncovered clear evidence for predator-driven change in the gonopodial complex, pointing to ecological agents partly responsible for the diversity of morphological variation associated with the novel ano-urogenital vertebral region in poeciliid fishes. Findings largely corresponded to our a priori predictions. Specifically, male Bahamas mosquitofish evolved predictable differences between predation regimes in hemal-spine structure, gonopodial anterior transposition, and lateral surface area of the gonopodium. We argue that these differences largely resulted from divergent natural and sexual selection on locomotor performance, body shape, and gonopodium size.
First, we found that diversification of male body shape in Bahamas mosquitofish-a well-known example of predatordriven parallel evolution-appears to be partially governed by changes in the ano-urogenital vertebral region. Natural selection favoring more energetically efficient steady-swimming abilities in low-predation environments where resource competition is particularly high has resulted in more streamlined bodies in these populations, while selection favoring greater acceleration during fast-start escape behaviors in the presence of predatory fish has resulted in bodies with larger mid-body caudal regions in high-predation populations (Langerhans et al., 2007;Langerhans, 2009a,b;Langerhans, 2018;Riesch et al., 2020; this study). Moreover, females show mating preferences within populations based on male body shape, with low-predation females preferring more streamlined males and high-predation females preferring males with larger caudal regions (Langerhans and Makowicz, 2013). We hypothesized that one way to partially develop these divergent morphologies in males involves modification of the hemal-spine morphology and transposition of the gonopodium: a more anteriorly positioned gonopodium creates a larger caudal region, while a more posteriorly positioned gonopodium creates a larger anterior region and may delay the separation of the boundary layer during steady swimming (Langerhans, 2010). Consistent with this notion, we found here that highpredation populations exhibited a suite of differences in hemalspine morphology that apparently resulted in greater anterior transposition of the gonopodium compared to low-predation populations. This highlights the role of predation regime in driving variation in this novel body region and suggests the region was co-opted to facilitate diversification of body shape after its initial evolution for genital support and mating functions.
The strongest evidence in support of this co-optation hypothesis was the observation that among populations there was a clear relationship between gonopodial anterior transposition and lateral body shape. The most parsimonious explanation for our results is that the ano-urogenital vertebral region, originally evolved for effective sperm transfer, served as an exaptation for body form and locomotor performance, and then evolved in response to selection on these attributes, reflecting secondary adaptation for a new function. That is, once the novel body region had evolved as an adaptation for rapid and effective sperm transfer, it inadvertently affected other attributes which selection later acted upon, leading to greater diversification in the region and novel evolutionary solutions. In this light, reexamination of past work on predator-driven morphological divergence in poeciliid fishes reveals that males often exhibit greater gonopodial anterior transposition in high-predation populations: e.g., G. affinis, G. caymanensis, Poecilia reticulata, B. rhabdophora Langerhans and Makowicz, 2009;Langerhans, 2010;Ingley et al., 2014). Meanwhile, female poeciliids also often exhibit morphological differences between predation regimes, but not a shift in the position of the anal fin (e.g., Langerhans, 2009b;Langerhans and Makowicz, 2009;Ingley et al., 2014), perhaps because they lack a readily available means of anteriorly transposing their anal fin. Combined with the widespread pattern of habitat-associated lateral morphology in poeciliid fishes (e.g., Hendry et al., 2006;Gomes and Monteiro, 2008;Tobler et al., 2008;Langerhans and Reznick, 2010;Torres-Mejia, 2011;Hassell et al., 2012;Ingley et al., 2014;Landy and Travis, 2015;Riesch et al., 2016Riesch et al., , 2018Langerhans, 2018;Moody and Lozano-Vilano, 2018), it seems that the co-optation of the novel ano-urogenital region could represent a common, but previously unrecognized, means of male body shape evolution in the family.
Prior work suggests females of several poeciliid fishes, including G. hubbsi, prefer males with larger gonopodia, but larger gonopodia increase drag and thus reduce fast-start escape performance (e.g., Langerhans et al., 2005;Kahn et al., 2010;Langerhans, 2011). Consistent with our prediction, and with empirical work in other habitat types and species, we found that low-predation populations exhibited a larger lateral surface area of the gonopodium. The evidence that a tradeoff between natural and sexual selection underlies this divergence was further strengthened by the fact that we did not find divergence in gonopodium length between predation regimes. This suggests that the results for the two estimates of gonopodium size were not confounded, and divergent selection acts on the surface area of the organ (via mating preference of lateral projection area and locomotor effects of frictional drag), not on its length. This further means that the changes observed between predation regimes in hemal-spine morphology and gonopodial anterior transposition cannot be explained as accommodations for changes in gonopodium length. Some prior studies have found longer gonopodia in high-predation environments [Kelly et al., 2000;Jennions and Kelly, 2002; but see Evans et al. (2011) for results similar to this study], perhaps explained by positive effects of length on insemination during coercive matings in some species Devigili et al., 2015;Head et al., 2015Head et al., , 2017. However, patterns of gonopodium length divergence between predation regimes appear inconsistent. The marked variation in gonopodium length across the family Poeciliidae has long attracted considerable research attention. Our findings add to this body of research, but suggest the primary causes of this variation are neither functional integration with gonopodial anterior transposition, nor the influence of either female mating preferences or natural selection via predation. Because most prior studies examining size of the gonopodium have only measured length, we caution researchers to consider whether in certain cases surface area might provide a more relevant measure.
We found that a combination of both genetically based differences and environmental influences appear to explain among-population variation in components of the gonopodial complex. Depending on the trait, we found moderate to strong population differences persisted in a common laboratory environment, suggesting evolutionary differentiation has occurred in all traits examined. Yet, the laboratory environment altered the expression of most traits, in at least some populations-the exception being gonopodial anterior transposition. Consistent with work in other poeciliid fishes (G. affinis: Langerhans, 2009b; Heterandria formosa: Landy and Travis, 2015), this particular trait showed strong genetically based population differences with little-to-no effects of laboratory rearing. On the other hand, the laboratory environment clearly had the greatest influence on gonopodium size, consistent with recent work showing that higher food levels result in shorter gonopodia relative to body size (Broder et al., 2020)-we observed relatively smaller gonopodia in the lab compared to the wild, and food was very likely in higher quality and quantity in the laboratory. But still, even gonopodium size showed genetically-based differences as well, matching prior work on heritability of gonopodium size within populations (Evans et al., 2013;Booksmythe et al., 2016) and genetically-based differences among populations in gonopodium morphology (Langerhans et al., 2005;Heinen-Kay and Langerhans, 2013;Broder et al., 2020). Future work can more carefully disentangle the roles of genetic divergence, plasticity, and population-specific plasticity in these traits, but it is clear all sources of variation may prove important for different traits and evolutionary changes in these traits have occurred among populations.

CONCLUSION
Among vertebrates, poeciliid fishes have evolved a unique body plan, and the ano-urogenital region appears to provide novel means of diversification. Subsequent to the evolution of the anourogenital region in poeciliids, the region has apparently been co-opted as a means of responding to natural and sexual selection on locomotor performance and lateral body shape, but the length of the gonopodium appears to evolve somewhat independently of the rest of the gonopodial complex. Changes in gonopodial anterior transposition may represent a common means of responding to selection on body shape in male live-bearing fishes, but peripheral consequences of such changes include a restricted trunk region which could lead to reshaping and repositioning of organs, decreased organ sizes, altered diets, and even changes in habitat use and foraging behaviors. These potential effects represent new hypotheses pointing to future research directions. Population divergence in body shape and genitalia can potentially influence reproductive compatibility, suggesting that evolution of diversity in the novel body region could play an important role in speciation. Evolution of the ano-urogenital region made new evolutionary solutions possible in poeciliid males. Poeciliid fishes exhibit marked variation in mating systems and habitat types, and future investigations at intraspecific and interspecific scales should yield important insights into the influence of the ano-urogenital region in diversification patterns in the family.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. GenBank accession numbers are provided in the article/Supplementary Material, and other datasets are available from the Dryad Digital Repository: https://doi.org/10.5061/ dryad.fbg79cnt8 (Langerhans and Rosa-Molinar, 2021).

ETHICS STATEMENT
The animal study was reviewed and approved by the Institutional Animal Care and Use Committee of North Carolina State University (protocols 13-101-O, 16-193-O).

AUTHOR CONTRIBUTIONS
RL and ER-M conceived the study and wrote the manuscript. RL collected and analyzed the data.