Beyond Pairwise Interactions: Multispecies Character Displacement in Mexican Freshwater Fish Communities

Competition has long been recognized as a central force in shaping evolution, particularly through character displacement. Yet research on character displacement is biased, as it has focused almost exclusively on pairs of interacting species while ignoring multispecies interactions. Communities are seldom so simple that only pairs of species interact, and it is not clear whether inferences from pairwise interactions are sufficient to explain patterns of phenotypes in nature. Here, we test for character displacement in a natural system of freshwater fishes in western Mexico that contains up to four congeneric species of the genus Poeciliopsis. We analyzed body shape differences between populations with different numbers of competitors while accounting for confounding environmental variables. Surprisingly, we found evidence for convergent character displacement in populations of P. prolifica, P. viriosa, and P. latidens. We also found that the convergence in body shape was not consistently in the same direction, meaning that when three or more competitors co-occurred, we did not find more extreme body shapes compared with when there were only two competitors. Instead, when three or more competitors co-occurred, body shape was intermediate between the shape found with a pair of species and the shape found with no competitor present. This intermediate shape suggests that evolution in multispecies communities likely occurs in response to several competitors rather than to simple pairwise interactions. Overall, our results suggest that competition among multiple species is more complex than simple pairwise competitive interactions.


Introduction
Species interactions are known agents of selection that can drive evolution. Several studies have demonstrated that trait evolution can be caused by species interactions, such as predation (Reznick and Endler 1982;Johnson and Belk 2001;Reznick et al. 2001;Ingley et al. 2014), parasitism (Buckling and Rainey 2002;Spottiswoode and Stevens 2010), mutualism (Gracia-Lázaro et al. 2018;Keller and Lau 2018), and competition (Ellis et al. 2015;Roth-Monzón et al. 2017). Competition is also known to drive speciation, in some cases leading to species radiations, and can be an important factor in the assembly of ecological communities (Pfennig and Pfennig 2012b). It has thus received much attention from both community and evolutionary ecologists (Pfennig and Pfennig 2010;Stuart and Losos 2013). In ecological communities natural selection is thought to act in a way that reduces competition, thereby promoting species coexistence (Pfennig and Pfennig 2012b).
Trait evolution that reduces competitive interactions is called "character displacement" (Brown and Wilson 1956;Grant 1972;Schluter and McPhail 1993;Schluter 2000). Three forms of character displacement have been described, each defined by a different type of competitive interaction: ecological character displacement (Brown and Wilson 1956;Pfennig and Pfennig 2010), reproductive character displacement (Butlin 1995;Gröning and Hochkirch 2008), and agonistic character displacement (Grether et al. 2009(Grether et al. , 2017. Here, we focus on ecological character displacement, which is caused by indirect competition for resources (food, habitat, etc.) and usually causes divergence between species to avoid niche overlap (Pfennig and Pfennig 2010). However, ecological character displacement can also cause convergence, leading to increased niche overlap between species (Grant 1972;Pfennig and Pfennig 2012c). Most ecological character displacement studies focus on pairwise interactions between species (but see Lemmon and Lemmon 2010;Miller et al. 2014;Grant 2017). Studies of complex competitive interactions (defined here as competition among three or more species) have received less attention, and it is still unclear how multiple species interactions drive trait evolution (terHorst et al. 2018).
Understanding complex species interactions is important because most ecological communities are composed of multiple interacting species. Moreover, recent evidence suggests that multispecies competition can be fundamentally different from pairwise species competition (terHorst et al. 2018). For example, one species can alter the effect that a second species has on a third species, creating a nonadditive interaction such that trait evolution cannot be inferred by simply knowing the effect of each pairwise interaction. Interestingly, the effects of multispecies competition can remain even when pairwise interactions are the dominant selective force because of the indirect effects of other co-occurring species (terHorst et al. 2018).
How then should character displacement occur where more than two competitors interact? Interestingly, two opposing predictions emerge, each dependent on whether the presence of additional competitors causes niche saturation (i.e., when the niche space is filled by species; MacArthur and Levins 1967;Scheffer and van Nes 2006). When the addition of competitors causes niche saturation, traits should converge or not shift at all. This is because when the niche is saturated by co-occurring competitors, it will constrain trait divergence, as the niche space will be unavailable for the focal species to diverge (Scheffer and van Nes 2006;Fox and Vasseur 2008;terHorst et al. 2010). If convergence occurs, the magnitude of phenotypic variation between species should be smaller when they co-occur with a competitor than when they are found alone. Moreover, this degree of convergence should decrease as the number of competitors increases (Hubbell 2006;Scheffer and van Nes 2006). Alternatively, if there is available niche space, then theory predicts that traits will diverge and that divergence will increase with the addition of more competitors; thus, the magnitude of phenotypic variation will increase between species when they co-occur with competitors compared with when they are found alone (Slatkin 1980;Abrams 1983;terHorst et al. 2018). Additionally, it has recently been proposed that in populations with high intraspecific variation, extinction or divergence may be more common if the populations are in multispecies communities (Barabás and D'Andrea 2016). In this final case, once again, the magnitude of phenotypic variation between species will be greater when species co-occur with competitors, either because of divergence or because of the extinction of some of the species (Hardin 1960;Vadas 1990; Barabás and D'Andrea 2016).
Here, we examine a natural system of four live-bearing fish species from the genus Poeciliopsis to determine whether and how character displacement occurs in a multispecies context. We evaluate body shape traits that have long been used in studies of character displacement and are known to be affected by competitive interactions (Schulter and Mc-Phail 1992;Adams and Rohlf 2000;Adams 2004;Husemann et al. 2014). Changes in body shape are expected in response to resource use and habitat use (Rüber and Adams 2001;Aguirre and Bell 2012;Meyers and Belk 2014). We focus on three fundamental questions in this system: (1) do these co-occurring species of Poeciliopsis show character displacement, (2) what is the pattern of character displacement among these Poeciliopsis species, and (3) how does the number of co-occurring species affect the magnitude of character displacement?

Study System
Poeciliopsis prolifica, P. viriosa, P. latidens, and P. presidionis co-occur throughout western Mexico on the Pacific slope from the Rio Yaqui, Sonora, south to near Las Varas, Nayarit (Miller et al. 2005; fig. 1). Although these four species belong to the same genus, they are not sister taxa. Both P. prolifica and P. viriosa belong to a strictly northern clade, whereas P. presidionis and P. latidens belong to predominantly southern clades. Northward dispersal of the ancestor of P. presidionis and P. latidens is hypothesized to be relatively recent (2.8-6.4 million years ago; Mateos et al. 2002). General accounts indicate that these four species are ecologically similar-they all inhabit the midwater column in streams and small rivers, they are similar in body size and shape, and they are omnivorous, consuming plant and animal matter (Miller et al. 2005). In fact, basic ecological descriptions of these four species coincide in their frequent association with green algae, which may be an important component of the diet of all these species (Miller et al. 2005). Furthermore, we collected all four species from the same microhabitats while conducting fieldwork. To our knowledge there are no published accounts of the degree of similarity or of the competitive overlap among these four species. However, because of their ecological similarity, we conclude that there is potential for competitive interactions. Hence, we use the number of co-occurring species as an indicator of the level of potential competition in our study.

Study Sites
We collected females of all species of Poeciliopsis with a handheld seine net (1.3 m # 5 m; 8-mm mesh size) during the dry season: October 2007 and November 2015. Our intent was to obtain as many independent replicates and combinations of co-occurrence of all four species as possible. This said, the number of replicates and combinations of species that were available in the field differed among species ( fig. 1; table 1). For each location we also gathered data on two potentially influential environmental variables: canopy cover and stream slope (table S1; tables S1, S2 are available online). Canopy cover and stream slope serve as proxies for resource availability and stream velocity, which can affect body shape (Langerhans and Reznick 2010;Scharnweber et al. 2013). We estimated canopy cover with the use of a handheld densitometer at each collection location. We calculated stream slope for each location in Arcmap (ESRI 2014). We calculated stream slope as the difference between upper elevation and lower elevation of a 2-km segment of stream (table S1). The locations evaluated here were similar in terms of other abiotic properties. For example, we found no differences among locations with different numbers of competitors for pH (F 2 p 0:719, P p :51), temperature (F 2 p 3:105, P p :08), and conductivity (F 2 p 3:278, P p :09). We also found no potential pisciv-orous predators in the locations sampled. However, in all locations (with the exception of site 3) we found another species of live-bearer (Poecilia butleri). We also found one  location (site 2) with a very low density (only 16 individuals collected) of an introduced live-bearer (Gambusia affinis).

Quantifying Body Shape
To quantify body shape we used a geometric morphometric approach (Rohlf and Marcus 1993). We digitized 14 biologically homologous landmarks and semilandmarks on a lateral image for each fish included in the analysis using the computer software tpsDig2 (Rohlf 2016a). Landmarks were defined as follows: (1) anterior tip of the snout, (2) anterior extent of the eye, (3) posterior extent of the eye, (4) semilandmark midway between landmarks 1 and 5, (5) semilandmark midway between landmarks 4 and 6, (6) anterior insertion of the dorsal fin, (7) dorsal origin of the caudal fin, (8) ventral origin of the caudal fin, (9) semilandmark midway between landmarks 8 and 10, (10) posterior insertion of anal fin, (11) semilandmark midway between landmarks 10 and 12, (12) semilandmark midway between landmarks 11 and 13, (13) anterior extent of the eye orbit, and (14) semilandmark midway between landmarks 12 and 4. Landmarks 1, 6, and 12 were used to allow digital unbending of the specimens (tpsUtil; Rohlf 2016b). We used generalized Procrustes analysis to remove all nonshape variation for each fish and to generate affine (i.e., shape variation that affects all landmarks in the same way) and nonaffine shape variables (W matrix; Mitteroecker and Gunz 2009). We summarized shape variation from the W matrix using a principal components analysis to generate relative warps in the package geomorph (Adams et al. 2017) in R programming software (R Core Development Team 2010; Adams et al. 2017). The 14 original landmarks yielded 24 reative warps. To account for the reduced dimensionality from the use of sliding semilandmarks and to avoid including shape variables that explain only minute amounts of shape variation, we used the first 10 relative warps for subsequent analysis, which combined accounted for 195% of observed shape variation. To characterize body shape variation among these four Poeciliopsis species, given that there are no previous data on body shape that can inform our study, we plotted the means and confidence intervals (CIs) for each species at each location on the first two relative warps (these two relative warps accounted for 66% of total shape variation). These plots showed two locations that appear to be outliers (location 11 for P. presidionis and location 14 for P. viriosa; fig. 2). These locations may have inordinately large effects on patterns of body shape. Thus, we conducted all analyses, including the generalized Procrustes analysis, with and without these two locations. Both sets of analyses yielded the same interpretation, so we present here the results from the analyses that include all locations.

Statistical Analyses
We used a two-step approach to evaluate variation in body shape that is due to competitor co-occurrence for all three species. First, for all populations of a single species we used a mixed model approach to test for variation in body shape due to the number of competitors. If results of the mixed model showed significant variation in body shape in response to the number of competitors, we proceeded with the second test-a phenotypic change vector analysis (PCVA)to determine whether the significant result found in our For the first step in testing for variation in the shape of focal species in response to the number of co-occurring competitor species, we used a mixed model multivariate analysis of variance (analyzed in SAS ver. 9.2; SAS Institute 2008). Incorporation of random effects in a multivariate analysis of variance requires vectorization of the response variable matrix such that each response variable occupies a single row of the data matrix (Anderson 2003). Thus, in our data set each individual is represented by 10 rows with a single response variable for each row (i.e., the first row includes relative warp 1 for the first individual, the second row includes relative warp 2 for the first individual, and so forth until all 10 warps are represented for the first individual and all remaining individuals). This vectorization creates a univariate model structure that can then incorporate random effects in a way similar to that of a standard univariate mixed model. The vectorization process requires including an index variable as a predictor in the model (e.g., Ingley et al. 2014;Meyers and Belk 2014;Heins and Baker 2017). We created the index variable by using the identifying order number of the relative warps (i.e., 1-10). Given that relative warps are orthogonal and ordered according to the amount of variation they explain, they can be treated as sequential descriptors of overall shape. To retain the order inherent in the relative warps, we use the index variable as a predictor variable (i.e., analogous to using time as an ordering variable in a traditional repeated-measures analysis).
We ran a separate model for each focal species. Shape was the response variable and was represented by the first 10 relative warps. The predictor variables were number of co-occurring competitor species, index variable, and the interaction between the number of co-occurring species and the index variable. Because relative warps are essentially principal components of shape variables, the direction of divergence of a species for one relative warp is independent of the direction of divergence for other relative warps. Thus, the index variable and the number of co-occurring competitor species by themselves have no useful interpretation.
Instead, it is the interaction term between the number of co-occurring competitor species and the index variable that tests the hypothesis of interest-that is, that the shape of the focal species varies for at least some of the relative warps in response to the number of co-occurring competitor species.
In addition, we included canopy cover and stream slope as covariates and their interaction with the index variable. Canopy cover has been used as a proxy for resource availability (Johnson 2002), and higher canopy cover has been shown to be correlated with lower resources in small tropical streams similar to ours (Grether et al. 2001). We used stream slope as a proxy for stream velocity (Zúñiga-Vega et al. 2007;Meyers and Belk 2014). Stream velocity is known to affect body shape in poeciliids-typically, higher-flowing streams favor a more streamlined body shape (Langerhans et al. 2004;Zúñiga-Vega et al. 2007;Langerhans 2008Langerhans , 2009Langerhans and Reznick 2010;Ingley et al. 2014). We view both canopy cover and stream slope as general descriptors of variation in the environment that can influence body shape. Our goal in including these two variables was to remove variation not related to the number of co-occurring competitor species to provide a more precise test of our hypothesis. Unfortunately, the sampling locations available did not allow a more thorough characterization of the influence of these two environmental variables on shape beyond what we present here.
Finally, we incorporated two random effects into the model. The first was a random effect to account for the structure among locations where fish were collected. Each collection location represented a specific level of competi-tion: no co-occurring competitors (allopatry), co-occurring with one species (sympatry), or co-occurring with two or more species of competitors (sympatry). The second random effect was included to account for the structure among the 10 relative warps (measures of shape) within individuals. Although among all individuals relative warps are orthogonal and independent, we still must account for a potential covariance structure within individuals among all relative warps. It is important to note that when using a vectorization approach to a multivariate analysis of variance, each shape measure provides a degree of freedom for each individual. Thus, the denominator degrees of freedom for the hypothesis tests are a product of the number of individuals and the number of shape variables (i.e., first 10 relative warps), with some adjustment depending on the specific method used for calculating degrees of freedom (Kenward and Roger 1997).
Because of differences in the number of locations for each of the species in this study, we coded the number of competitors in two different ways to allow replicates for each competition level. First, we designated the number of competitors as 0, 1, or 21 (meaning two or more competing species), and second, we designated the number of competitors as 0 or 11 (meaning one or more competing species). For P. prolifica and P. latidens we used two separate models utilizing each of the categorization schemes for competition. For P. viriosa we were able to analyze only the model with the simpler scheme (0 and 11) as the main effect because of the number of sampled locations (table 1). For all three focal species in both competition schemes there were no differences in the inclusion of competitors in the competitive levels; all three species were part of each of the competitive levels. Unfortunately, because of the lack of allopatric localities (0 competitors) for P. presidionis, we were unable to conduct a formal test of the number of competitors for this species (table 1). However, we did include P. presidionis samples for general shape comparison in another analysis to show its position in shape space. When we found a significant interaction between the number of competitors and the index variable in the models-which indicated that there was significant variation in body shape due to number of competitors-we proceeded with the second step: a PCVA to determine the specific nature of the interaction (Adams andCollyer 2007, 2009;Collyer and Adams 2007). We calculated the magnitude of shape difference between species in each competitive level using all 10 relative warps. This allowed us to test whether the significant interaction resulted in convergence, divergence, or no change in all three species. When the magnitude of shape differences among species when species were alone (allopatry) was larger than the magnitude of shape differences among species when they co-occurred with competitors (sympatry), we concluded that convergence in body shape had occurred. However, when the magnitude of shape difference in sympatry was larger than the magnitude of shape difference in allopatry, we concluded that divergence in body shape had occurred. Finally, no difference in the magnitude of shape change among species between sympatric and allopatric populations was indicative of no character displacement. Additionally, we conducted a PCVA between the species to understand whether all three focal species had a similar response in the magnitude of shape change to the number of co-occurring competitor species (Adams andCollyer 2007, 2009;Collyer and Adams 2007). Thus, if the magnitude of shape change was the same for all species between sympatric and allopatric populations, then we consider all three focal species to have responded similarly to the presence of a competitor. We conducted all of the PCVAs using ASReml-R version 3.0 (Butler et al. 2009). To visualize the differences in body shape for all analyses we plotted least squares means (95% CI) on the first two relative warp axes, which accounted for more than half of the shape variation (65.7%). We generated thin-plate spline visualizations (Zelditch et al. 2012) to represent the observed changes in body shape and to aid with the interpretation of the differences found for the number of competitors.
Before conducting all of these analyses, we tested for a potentially confounding effect of year, given that we collected specimens in 2007 and 2015. To do so we implemented linear mixed models using each relative warp as the response variable, year as a fixed effect, and populations nested within years as a random factor. We used a nested design because each population was sampled only once; therefore, we could not test for interannual differences within each population. Instead, these mixed models allowed us to test for overall changes in shape between collection years for each species. That said, only two of our four focal species were collected in two different years (P. prolifica and P. viriosa) and were amenable to interannual comparisons. In all but one case we found no effect of collection year (P 1 :07 in all of these cases). The only exception was the seventh relative warp for P. viriosa, which explained only 3% of the total variability in body shape. Therefore, we concluded that shape did not differ between collection years. This is consistent with prior work in which short-term temporal changes (within less than 10 years in our case) in body shape were uncommon compared with spatial (interpopulation) variation (Zúñiga-Vega et al. 2007;Langerhans and Makowicz 2009).
Although all three focal species were present in both competition schemes such that there were no differences in the species identity of competitors in all of the competitive levels, we conducted a test to see whether there is variation in body shape due to the particular identity of the cooccurring conspecific. This test was possible only for P. prolifica. We used only the localities in which P. prolifica occurred alone or with a single competitor. We followed the same mixed model approach but considered the identity (ID) of the conspecific (none, P. latidens, P. viriosa, and P. presidionis) and its interaction with the index variable as the predictor variable of interest. This allowed us to test whether there was significant shape variation due to the identity of the co-occurring conspecific in at least some of the 10 shape variables (i.e., relative warps). All 10 relative warps served as the response variable. In addition to conspecific ID and the index variable, we included both canopy cover and stream slope as covariates, and location was included as a random effect as in previous models. When a significant effect of conspecific ID (i.e., the interaction between conspecific ID and index variable) was found, we conducted pairwise post hoc tests to understand which species caused the difference in the shape of P. prolifica. We adjusted the P value for the post hoc tests for multiple comparisons (Tukey-Kramer adjustment). However, our test was limited by the sampling because there is only one location where P. prolifica co-occurs only with P. viriosa and there are some missing data for the covariates (table S1). Thus, these results should be interpreted with caution.

Results
There was high overlap in body shape among all four species across all locations ( fig. 2). When we removed the two outlier locations, three of the four species showed nearly complete overlap with the other species in shape space, as seen for the first two relative warps ( fig. 2). Yet despite this overlap among species, variation in body shape within species and locations was quite broad.
Each species showed a shift in body shape when it occurred with competitors. For both P. prolifica and P. latidens we found a significant effect of the interaction between the index variable and the number of competitors (coded both as 0, 1, or 21 and as 0 or 11). Both canopy cover and stream gradient also showed significant effects on body shape (tables 2, 3). For both P. latidens and P. prolifica body shape converged in the presence of competitors ( fig. 3). Both species are less similar to each other in shape in allopatry (MD 0 p 0:0422, P ! :0001) than when they co-occur with competitors (MD 1 p 0:0192, P p :7989; MD 21 p 0:0223, P p :2193). Overall, the body shape of both P. prolifica and P. latidens when they cooccurred with other species was less streamlined and more robust ( fig. 4) than when they occurred alone. However, some differences in body shape depended on whether they co-occurred with one competitor or with two or more competitors ( fig. 4). When P. prolifica and P. latidens cooccurred with one competitor, they had a more robust shape with a deep belly, a more dorsally oriented mouth, and a smaller head and eye compared with their allopatric shape. However, when they co-occurred with two or more competitors, the belly shape was intermediate between the shape with zero competitors and the shape with one competitor, but the eye was larger than that of the allopatric shape ( fig. 4).
Poeciliopsis viriosa also showed a shift in body shape when it occurred with competitors-there was a significant effect of the interaction between the index variable and the number of competitors for P. viriosa after accounting for environmental differences. Additionally and similar to the previous analysis, both P. latidens and P. prolifica differed significantly in shape when compared between zero competitors and one or more competitors (table 3). Similar to previous results, all three species converged in body shape when they co-occurred with at least one competitor (fig. 5). The magnitude of shape differences among all three species in allopatry (MD 0 p 0:0429, P p :0001) was greater than the magnitude of shape differences among species when they co-occurred with competitors (MD 11 p 0:0236, P p :5421). The body shape for all species when they co-occurred with one (or more) competitor was more robust with a deeper belly, a larger eye, and a less pronounced snout and a more dorsally oriented mouth compared with the shape when no competitor was present ( fig. 6).
Not all species responded to competition with the same magnitude of shape change (figs. 3, 5). If we consider the body shape without competition as a baseline compared with the body shape in the presence of at least one competitor, P. prolifica had a greater magnitude of change in body shape (MD 0-11 p 0:0335, P p :0012) compared with  Note: See "Methods" for an explanation of the 0 and 11 categories. Boldface type indicates statistical significance (P ! :05). P. latidens (MD 0-11 p 0:0295, P p :0740), and P. latidens had a greater magnitude of change in body shape compared with P. viriosa (MD 0-11 p 0:0237, P p :6213; figs. 3, 5). This order of the amount of change in body shape was also found when the outlier locations were removed ( fig. S1, available online).
The identity of the co-occurring species had a significant effect on the body shape of P. prolifica (conspecific ID# index: F 27, 2,253 p 36:30, P ! :0001). Canopy cover also had a significant effect on the body shape of P. prolifica independent of the effect of conspecific ID (canopy cover: F 1, 1,354 p 10:84, P p :001; canopy#index: F 9, 1,431 p 9:04, P ! :0001). No other main effects or covariates were significant predictors of body shape in P. prolifica (conspecific ID: F 3, 1,354 p 2:34, P p :0716; stream slope: F 1, 1,354 p 0:83, P p :3612; stream slope#index: F 9, 1,431 p 1:21, P p :2849). For at least some of the relative warps all of the pairwise contrasts testing for differences in the shape of P. prolifica in allopatry and when co-occurring with a given competitor species were significant (table S2). The body shape of P. prolifica alone was most different from the body shape of P. prolifica when it co-occurred only with P. latidens. However, the body shape of P. prolifica in the presence of either P. presidionis or P. viriosa differed little from the body shape of P. prolifica alone (table S2).
Discussion Surprisingly, we found that all species converge in body shape when they occurred with potential competitors. We found a smaller magnitude of shape differences among species when they co-occurred with other competitors compared with the variation in shape when they were found alone (allopatry), indicating convergence for all three focal species. Convergence is uncommon in studies of character displacement (Grant 1972;Pfennig and Pfennig 2012a) and is usually explained by shared environmental effects (Verde Arregoitia et al. 2018). However, the convergence we observed in these Poeciliopsis fishes remained even after accounting for potential environmental drivers of body shape (resource availability and stream velocity). It is possible that other environmental variation that we did not measure-such as nonpscivorous predators, microhabitat structure, and fine-scale water velocity-may have caused some of the variation that we observed in body shape. However, in light of our current evidence the pattern appears to be predominantly driven by the co-occurrence of competitors, a phenomenon known as convergent character displacement. Furthermore, the shift in body shape that we observed is not consistent with the expected effects of either resources or stream velocity. For example, we found that fish from populations with two or more competitors had a less distended abdomen even though these sites tended to have a more open canopy cover (more resources), which should lead to a more distended abdomen (Olsson et al. 2007).
Theory suggests that the pattern of convergence we found is consistent with a niche that is saturated with insufficient space to allow species to diverge (Hubbell 2006;Scheffer and van Nes 2006;terHorst et al. 2010terHorst et al. , 2018. Interestingly, the convergence in shape we found was not consistently in the same direction as the number of competitors increased, meaning that when a species occurred with two or more competitors, we did not find a more extreme shape than when that species occurred with only one competitor. Instead, when two or more competitors co-occurred, we found a body shape that was intermediate to the shape found with one or no competitors. This intermediate shape matched the result of a previous study with multiple competing species wherein species with multiple competitors tended to have intermediate competitive phenotypes (Miller et al. 2014). It is possible that in multispecies communities pairwise interactions independent of other species are rare, simply because species interact with several competitors and respond with phenotypes suited to deal with all of these competitors (Connell 1980;Strauss et al. 2005). Thus, the multiplicity of interactions could limit adaptation to any particular co-occurring species, instead resulting in a compromised phenotype suited to dealing with multiple competitors. Body shape in fishes is known to be affected by several biotic and abiotic factors (Langerhans and Reznick 2010). Each of those factors can create specific morphological responses (Langerhans and Reznick 2010). What does the body shape found in our study tell us about these three Poeciliopsis species? In both P. latidens and P. prolifica we found that when co-occurring with one competitor, individuals had a deeper body and a shorter and wider caudal peduncle. We also found that the eye was smaller and the mouth appeared to be more dorsally located. This shape, in general, does not match predictions of benthic versus limnetic phenotypes (Robinson and Wilson 1995;Ruehl and DeWitt 2005;Palkovacs et al. 2011), lower resources (Olsson et al. 2007, or higher stream velocities (Langerhans et al. 2004;Langerhans 2008Langerhans , 2009Ingley et al. 2014). However, all shape differences that were found (except for the size of the eye) show great similarity to the body shape that has been associated with the ability to have a higher resource acquisition rate (Palkovacs et al. 2011). Furthermore, when both P. latidens and P. prolifica co-occur with two or more competitors, the only changes are a larger eye and a shallower body, but the body still remained deeper than the body shape when no competitors co-occur-this result is also consistent with this ability for a higher rate of resource acquisition. These differences in body shape were also found in P. viriosa when the comparison was only between no competitors and one or more competitors. Thus, these results all suggest that competition could be driving body shape to allow individuals to improve their ability to acquire more resources. These results lend further support to the theory that convergence is due to increasing competitive abilities (Pfennig and Pfennig 2012b), so that all species are converging in a body shape that is associated with allowing a higher acquisition of resources.
Our results also support contemporary coexistence theory in that sympatry can result in reduced fitness differences rather than increased niche differences (Chesson 2000;Germain et al. 2018

Zero competitors
One plus competitors P. viriosa Figure 6: Body shape for Poeciliopsis latidens, P. prolifica, and P. viriosa for zero versus one or more (0 and 11) competitors. Thin-plate spline grids are shown at a scale of 3#, with lines connecting landmarks added to assist interpretation.
focused on head morphology, particularly on mouth characteristics that can be affected by resource consumption or by the particular type of food items that these fishes consume, should yield further insights to test these hypotheses (Tobler 2008;Palkovacs et al. 2011). Finally, the asymmetric degree of displacement in body shape found here suggests that each species could be affected differently by competition, which could explain the observed differences in the degree of shape response. Asymmetries in character displacement are usually posited to be due to differences in competitive abilities among the species considered (Pfennig and Pfennig 2012b). It is possible that the greater shift of P. prolifica when co-occurring with competitors, as opposed to when occurring alone, is due to the effects of the presence of other competitors. Other cooccurring competitors may lower the competitive effect of each species on P. prolifica (perhaps through a nonhierarchical interaction) but not with the other competitors. Another possibility is that each species has a species-specific effect on its competitors. Our analysis of conspecific ID supports this hypothesis. However, it is unclear with such a limited sample and the potential confounding variation (table S1; single locality of P. prolifica co-occurring with P. viriosa) how strong the effect is and whether all three species cause such differences. Hence, future studies focusing on competitive abilities and each pairwise interaction, as well as on different combinations of multispecies interactions, will be needed to more fully explore the causes of these asymmetries.
Future studies could also consider laboratory tests of performance in the food consumption of each species from the different locations with different numbers of co-occurring species. This may shed some light on the intermediate shape of P. prolifica and P. latidens when co-occurring with two or more competitors. Furthermore, such experimental results could strengthen the argument that when more competitors coexist, responses will be intermediate to account for multispecies interactions instead of pairwise interactions.
One caveat is that we cannot exclude the possible effects of correlated selection on other traits, most notably life history (Langerhans and Reznick 2010). Studies have shown that coexistence can be maintained in competitive environments through trade-offs with life-history traits (Hutchinson 1957;Levine and Rees 2002;Leibold et al. 2004;Calcagno et al. 2006;Chapuis et al. 2017). Life history can contribute to differences in body shape in poeciliids (Plaut 2002). In particular, a deeper body with a larger belly could be the by-product of adaptation to a large reproductive investment (Plaut 2002;Ghalambor et al. 2004;Langerhans and Reznick 2010). Nevertheless, all species included in this study have superfetation (i.e., the ability to carry multiple broods simultaneously), which has been proposed as an adaptation to maximize reproductive investment without a trade-off in morphology (Zúñiga-Vega et al. 2007, 2010. Furthermore, reproductive allocation itself can be molded by trade-offs between selection for reproduction and the mobility required for foraging (Ghalambor et al. 2003(Ghalambor et al. , 2004. However, if there is an effect of life history on body depth, it would not discount the effects found here that point to competition as a driver of body shape. We assume, in accordance with several studies of evolutionary change of body shape, that differences in body shape are underpinned, at least to some degree, by genetic changes among populations. Although we did not directly test this and cannot exclude the possibility of plasticity, shape is typically driven to some extent by heritable variation in poeciliid fishes (Reznick et al. 2008). Moreover, studies of female guppies (closely related to our species) have failed to induce a plastic response on body shape due to different levels of resource availability, while such responses did occur in male guppies (Robinson and Wilson 1995).

Conclusion
We found surprising evidence for convergent character displacement in body shape in three species of Poeciliopsis. Overall, our results suggest that competition in multispecies interactions is more complex than that in simple pairwise comparisons. The important conclusion from our work is that when multiple competing species co-occur, simply considering pairwise competitive interactions to explore trait evolution is likely insufficient. In our work we found a nonintuitive pattern of trait convergence that makes best sense only in the context of multiple competing species.