## Similar protocols

## Protocol publication

[…] In the first set of analyses, we aimed to compare the different suites of traits (i.e., gonopodial and nongonopodial) in terms of evolutionary rates and modes of trait evolution. All phylogenetic comparative analyses were performed using the best supported tree in Jones et al. () as a reference tree. This reference tree is based on a set of RAD markers and was estimated using maximum likelihood. This tree was transformed into an ultrametric tree using the chronopl function in the R package **ape** (Sanderson ), with smoothing parameter set to 1. Then, based on the multistate character datasets (gonopodial and nongonopodial), we computed a matrix of pairwise Gower's distances (Gower ) between species using the R package cluster (Maechler et al. ) and restricted the analyses to the traits scored in at least half of the species (we note that some species were not scored for all traits in the previously published data utilized here) (all traits listed in Table were included in these analyses). Next, we performed a principal coordinates analysis on each of these two matrices, retaining the score of each species on the first principal coordinate (accounting for 51.02% of total variation in the case of gonopodial‐related traits and 55.54% in the other set of traits) as a univariate measure of trait variation for the subsequent univariate analyses. We tested for the presence of **phylogenetic** signal in the multivariate datasets comprising the scores along all the principal coordinate axes for the two datasets (gonopodial and nongonopodial). This was accomplished using a method recently proposed by Adams (), which consists of a generalization of Blomberg's K statistic (Blomberg et al. ) to multivariate data and whose significance is tested through a permutational procedure (1000 permutations in our case; see Table for a summary of all analyses conducted in this study).Next, to determine the evolutionary dynamics of both the gonopodial and nongonopodial trait sets, eight models were fitted and the rates of trait evolution were compared between the two sets of traits (Adams ). We used Adams’ () method to compare the evolutionary rates between the first principal coordinate computed on the distance matrix based on gonopodial traits, and the first principal coordinate based on the other traits. We employed the R package **geiger** (Harmon et al. ) to fit different evolutionary models on each of the two principal coordinates. To identify the best‐fitting model, a model selection procedure was used. First, a likelihood ratio test was performed to compare a Brownian motion model (i.e., a random walk model with a constant rate of trait evolution; Felsenstein ) with a model of white noise to determine whether a phylogenetic model of trait variation represented a significant improvement over a model of random noise. Then, as the Brownian motion model was significantly better in both cases, the other models available in the function fit Continuous were fitted and compared to the Brownian motion model using a likelihood ratio test. These comprise the Ornstein–Uhlenbeck model (which is a random walk with an optimum in phenotypic space, toward which the evolution of the trait is “pulled”; Butler and King ), an early‐burst model (where evolutionary rates increase or decrease exponentially through time; Harmon et al. ), a trend model (where evolutionary rates increase or decrease linearly through time), and three models (lambda, kappa, and delta) based on tree transformations (Pagel ). The lambda model transforms the tree according to a parameter lambda, which ranges between zero (star‐like phylogeny, which implies that the evolution of the trait is not reflected by the phylogeny) and one (equivalent to a Brownian motion model). The kappa model differentially “stretches” longer and shorter branches; in its default implementation in geiger, it is a punctuational model of evolution, with values bounded to be comprised between zero (punctuational model, where the amount of evolution is independent of branch length) and one (no differential “stretching” of branches). In the delta model, based on a scaling of the path lengths, the rates of evolution can increase or decrease over time. When models fitted using default options in fitContinuous contained estimated parameters at their default bounds, the model was fit again increasing the range of the parameter used by the fitContinuous function. Among the models that fitted significantly better than the Brownian motion model (if any), the best was chosen using the version of the Akaike's Information Criterion (AIC; Akaike ) corrected for small sample sizes (AICc; Hurvich and Tsai ).With the aim of conducting a preliminary investigation of whether sexual selection is acting on gonopodial traits, we implemented the same analyses described above to compare the rates and modes of evolution in gonopodial traits and a subset of nongonopodial traits. We compared gonopodial traits with nongonopodial traits reasonably known to be under sexual selection (vertical bars and growth rate, e.g., Ryan and Causey ; Morris et al. ; Lampert et al. ) and for which data are available. We do not include the sword trait (known to be preferred by females) in this subset because the evolution of the sword involves a variety of factors. For example, in some species, this trait has been lost (Xiphophorus maculatus and Xiphophorus variatus); however, females of both species prefer males with a sword; therefore, it is difficult to accurately reflect this scenario in a presence/absence matrix, for example. We additionally compared this subset of nongonopodial traits putatively under sexual selection with a subset of nongonopodial traits where the selection mechanisms acting are unknown to date (head bump, multiple lateral stripes, solid mid‐lateral stripe at birth, body bicolored, dark subdermal dashes of pigment, two or more oblique lines behind pectoral base; Table ). This is a preliminary investigation as to date most morphological traits differentially exhibited among Xiphophorus species are yet to be identified as being under selection or evolving neutrally. [...] To determine whether the variation in specific gonopodial traits is correlated with habitat type, that is, sites with different water flow regimes such as ponds versus flowing rivers, we used habitat data descriptions from all existing studies where water flow has been characterized for Xiphophorus habitats (Rosen ; Rauchenberger et al. ; Meyer and Schartl ; Kallman et al. ; Kallman and Kazianis ; Jones et al. ), as well as from unpublished data collected and verified over 35 years of regular field studies (M. Schartl, unpublished data). We note that in some instances although different species have been recorded to inhabit the exact same rivers or streams, they have also been repeatedly observed to prefer different microhabitats of those waterways (M. Schartl, unpublished data). For example in the habitats where Xiphophorus kallmani and Xiphophorus milleri predominantly occur, the swordtails (X. kallmani) are always seen in the middle of the stream where the current is high, and they also court in this habitat (MS, pers. obs.). In contrast, the platyfish (X. milleri) are only found in the calm regions of the streams, generally close to the shore and under plants (MS, pers. obs.). The same holds true for X. variatus and the northern swordtails. In such cases, species repeatedly recorded in the faster‐flowing regions of rivers or streams were categorized as occurring in flowing habitat types, whereas species repeatedly recorded close to the banks and under plants in slower‐flowing regions of the waterways were categorized as occurring in still‐water habitats. We categorized all known habitat types as either flowing or still water and then used phylogenetic comparative methods to test for morphological differences between habitat types in traits deemed likely to be influenced by water flow (due to the fact that they are external structures on the gonopodium). Of the major clades, the claw character is present in 16 of 17 species from the two clades typified by flowing water environments, while the claw is present in only 1 of 9 species from the clades most commonly in still‐water environments (Fig. S1). We measured a further set of five morphometric traits on the putative holdfast traits, the claw and serrae (Fig. S2, these are linear measurements, different from the multistate gonopodial characters used as starting data above), computed species means, and adjusted for allometric variation using standard length (sample mean). We chose to utilize the claw and serrae for these analyses as these features are on the external part of the gonopodium and may have holdfast functions and contribute to copulatory compatibility. All the subsequent phylogenetic comparative analyses are based on the same ultrametric tree described above for the analyses using multistate characters as starting data.We tested for phylogenetic signal, that is, the tendency for evolutionary‐related organisms to resemble each other (Blomberg et al. ), in the morphometric traits on the putative holdfast traits using both a Mantel test and the adaptation of Bloomberg's K to multivariate data (Adams ). The Mantel test was used to test the significance of the correlation of allometry‐adjusted pairwise Euclidean morphometric distances with the matrix of **patristic** distances obtained from the phylogenetic tree: The same phylogeny was used for Adams’ method.We used phylogenetic generalized least‐squares method (Grafen ; Martins and Hansen ; Garland and Ives ; Rohlf ) to take into account phylogenetic nonindependence when comparing habitat types using the five morphometric measurements as dependent variables. For phylogenetic generalized least‐squares method, we used the expected covariance matrix under a Brownian motion model (with gamma parameter set to 1, obtained in ape) as the error covariance matrix. To ensure the consistency between the analyses here and those detailed below for tests of reproductive character displacement, we also obtained pairwise interspecific Euclidean morphometric distances based on the five morphometric traits (Fig. S2) after they had been subjected to a multivariate regression‐based allometric adjustment. We then used a partial Mantel test (Smouse et al. ; Oden and Sokal ) to test for the correlation between these distances and a binary matrix indicating whether two species live in the same environment or not. To account for phylogenetic nonindependence, we kept the matrix of pairwise patristic distances constant.Additionally, we asked whether genital evolution is influenced by the avoidance of interspecific hybridization. We addressed this question by comparing the differences in gonopodia of species pairs known to hybridize or not in nature and the laboratory. We asked whether or not those pairs that are sympatric in nature have more pronounced differences in gonopodial structure than pairs that are allopatric in nature. We utilize extensive interspecific hybridization records (both under laboratory conditions, Schartl et al. unpublished, and naturally hybridizing species, summarized in Kallman and Kazianis ()], as well as species geographical distribution information including sympatric and allopatric data (Tables S3, S4). We investigated sympatry and hybridization using, as outlined above, partial Mantel tests. These tests were implemented because sympatry and hybridization events can be expressed only as a property of species pairs and we could therefore not use the phylogenetic generalized least‐squares method to test for difference in the five morphometric traits. Specifically, we tested for the correlation between the matrix of pairwise morphometric distances (after allometric correction) and a binary matrix reflecting, respectively, if each pair of species lived in sympatry or not, if each pair of species hybridized under laboratory conditions, and if each pair of species hybridized under both laboratory and natural conditions (see Tables S3 and S4: data compiled from Rosen ; Meyer ; Kallman et al. ; Kallman and Kazianis ; M. Schartl pers obs.). As above, the matrix of patristic distances obtained from the phylogeny of Jones et al. () was used to account for phylogenetic nonindependence in all tests.We performed the above‐mentioned set of comparative analyses (phylogenetic generalized least‐squares test for comparing water flow regimes; partial Mantel tests for assessing the correlation of morphology with hybridization and sympatry), also on gonopodium length both accounting for allometric variation (using standard length as covariate) and using raw data.Phylogenetic comparative analyses were performed using the R (R Core Team ) packages ape (Paradis et al. ), nlme (Pinheiro et al. ), vegan (Oksanen et al. ), and **adephylo** (Jombart and Dray ). All analyses using partial Mantel tests are based on 1000 permutations. […]