Morphological divergence rate tests for natural selection: Uncertainty of parameter estimation and robustness of results

Leandro R. Monteiro; José L. Gomes-Jr

Universidade Estadual do Norte Fluminense, Centro de Biociências e Biotecnologia, Laboratório de Ciências Ambientais, Campos dos Goytacazes, RJ, Brazil

Correspondence

ABSTRACT

In this study, we used a combination of geometric morphometric and evolutionary genetics methods for the inference of possible mechanisms of evolutionary divergence. A sensitivity analysis for the constant-heritability rate test results regarding variation in genetic and demographic parameters was performed, in order to assess the relative influence of uncertainty of parameter estimation on the robustness of test results. As an application, we present a study on body shape variation among populations of the poeciliine fish Poecilia vivipara inhabiting lagoons of the quaternary plains in northern Rio de Janeiro State, Brazil. The sensitivity analysis showed that, in general, the most important parameters are heritability, effective population size and number of generations since divergence. For this specific example, using a conservatively wide range of parameters, the neutral model of genetic drift could not be accepted as a sole cause for the observed magnitude of morphological divergence among populations. A mechanism of directional selection is suggested as the main cause of variation among populations in different habitats and lagoons. The implications of parameter estimation and biological assumptions and consequences are discussed.

Key words: divergence rate test, natural selection, Poecilia vivipara, statistical shape analysis, geographic variation.

Introduction

The inference of evolution or divergence rates has become an important tool in the assessment of causes and mechanisms of morphological diversification (Spicer, 1993; Hendry and Kinnison, 1999; Kinnison and Hendry, 2001). The development of neutral models of expected change (Lande, 1977; Bookstein, 1988; Turelli et al., 1988; Lynch, 1990) was also instrumental in the establishment of quantitative procedures for the investigation of proposed mechanisms of evolutionary change. The rate tests for neutral evolution hypotheses commonly used are the constant heritability of Lande (1977), the mutation-drift equilibrium of Turelli et al. (1988), and the delta of Lynch (1990). These tests compare the observed among-group (species or populations at the same or different times) variation with the variation expected by a neutral model of genetic drift, calculated from the within-population variation, genetic and demographic parameters. Whenever the rate of divergence or evolution is higher than expected by the neutral model, a mechanism of directional selection is suggested. On the other hand, if the rate is lower than expected, stabilizing selection is suggested as a cause for the lack of divergence (Spicer, 1993; Diniz-Filho, 2000). Studies on multidimensional shape usually lack estimates of evolutionary rates, probably because of the difficulty of calculation and comparison of different estimates among studies (Hendry and Kinnison, 1999), and because specific multivariate tests for selection and drift have been proposed (Lande, 1979; Lofsvold, 1988), which are supposedly independent of time and population sizes.

Because most studies published to date have investigated divergence or evolutionary rates among species or populations separated for a long time, stabilizing selection is commonly reported as a mechanism operating to reduce morphological diversification (Lynch, 1990; Spicer, 1993; Lemos et al., 2001). Genetic drift has not been rejected in some cases (Lynch, 1990; Cheetham et al., 1993), but studies claiming directional selection (such as Clegg et al., 2002) as a mechanism are rare, probably because of the reversals in directional selection patterns during longer periods of time.

One problem commonly encountered in divergence rate studies is the uncertainty in the estimation of parameters such as heritability, effective population size, number of generations since divergence, among- and within-group variances (Spicer, 1993). Because of this uncertainty, the rate tests have been considered as qualitative evidence of natural selection and a recommendation to use conservative estimates was set forth (Turelli et al., 1988; Spicer, 1993). In spite of the problems with parameter estimation, no study has attempted to quantify the influence of the different parameters on the test results. Sensitivity analysis of model results variation caused by parameter variation has been a widely used tool for the assessment of the influence of parameters in model results, particularly in ecological and evolutionary modeling of populations (Dowd, 1997; Vickery and Poulin, 1998). This analysis quantifies the magnitude of the influence of each parameter independently from the others on the model results and indicates which should be estimated more carefully.

Materials and Methods

Study system

The quaternary plains of Northern Rio de Janeiro State (Figure 1) were formed by sediment deposition during the formation of the Paraíba do Sul River delta, and sea level fluctuations associated with paleoclimatic modifications during the Holocene (starting approximately 5100 years ago - Martin et al., 1997). The process of sediment deposition and the construction of the river delta continuously formed river arms, which were later abandoned and became long, narrow lagoons scattered all over the region (Martin et al., 1997; Sofiatti, 1998; Primo et al., 2002). Until the beginning of the XX century, there was considerable contact among most lagoons and the main river body. However, the construction of a large network of drainage and irrigation channels by the National Department of Work Against Drought (DNOS) after 1950, greatly decreased the amount of water in the region and isolated most lagoons from all others (Primo et al., 2002). The distance from the sea determined a major amount of variation in the environmental conditions of the different lagoons, and also large environmental gradients within lagoons (Suzuki et al., 1998; 2002). One of the oldest and largest lagoons in the region is the Campelo Lagoon, that is around 4000 years old and is distant 17 km from the sea (but close to the current Paraíba do Sul River channel). This is a freshwater lagoon with large banks of macrophytes (usually Typha) along the shallow margins. Although the Campelo Lagoon is large (12 km² water surface), it's environment is homogeneous throughout its extension (probably because of its distance from the sea). The fish species found in the Campelo Lagoon are characteristic of freshwater environments. The Açu Lagoon is a coastal lagoon and was formed from the ancient Iguaçu River, which once connected the Feia Lagoon (the largest freshwater lagoon in Brazil) and the sea (Primo et al., 2002). The Açu Lagoon is younger than the Campelo Lagoon (around 3000 years old) and is separated from the sea by a small sand bar. The closeness to the sea causes an environmental gradient along its extension, where the region closer to the sea (the sand bar region) presents higher amounts of salinity, no marginal vegetation or macrophytes, and fish species characteristic of marine or estuarine environments. The interior of the lagoon presents lower (but not freshwater) salinity, and mangrove vegetation or macrophytes such as Typha. Among the most recent (1500 years or less) are the coastal Grussaí and Iquipari Lagoons, which were formed by the closure (by a sand bar) of ancient fluvial channels connecting the Paraíba do Sul River and the sea (Sofiatti, 1998). The Grussaí and the Iquipari Lagoons present similar environmental gradients also. The sand bar region presents high salinity, particularly when the sand bar is artificially opened (Suzuki et al., 1998; 2002). The environmental gradient in these lagoons is similar to the one in the Açu Lagoon, but the salinity drops considerably towards the interior of the lagoons (going from saltwater to freshwater). The banks of macrophytes (mostly Typha) are absent in the sand bar regions, but abundant in the interior. The fish community is also different within the same lagoon (large freshwater predators such as Hoplias malbaricus are absent from sand bar regions; Bizerril and Primo, 2001).

The livebearing fish Poecilia vivipara Bloch and Schneider 1801 is widespread in the lagoons of northern Rio de Janeiro. Its high tolerance to environmental changes, particularly salinity and temperature (Trexler, 1989), makes it one of the few species present and abundant in all lagoon environments of the region (Bizerril and Primo, 2001). In fact, it can be found even in very small lentic environments with low levels of dissolved oxygen (Gomes-Jr and Monteiro, unpublished data). This system offers an interesting opportunity for the study of body shape differentiation processes among populations, and there is evidence of large components of shape variation, even among sites within the same lagoon (Neves and Monteiro, 2003).

Field and laboratory methods

From each site, a number of specimens was collected, fixed and stored in 10% formalin. Only adult (sexually mature) animals were included in the samples. All females were carrying developing embryos, and the males had a fully developed gonopodium and visible secondary sexual characters. The final sample size from each site was: Açu Sand bar = 60 males and 98 females; Açu Interior = 60 males and 97 females; Grussaí Sand bar = 60 males and 46 females; Grussaí Interior = 60 males and 96 females; Iquipari Sand bar = 60 males and 65 females; Iquipari Interior = 60 males and 97 females; Campelo = 60 males and 100 females. The specimens were photographed with a high-resolution digital camera (Pixera), and the coordinates of 12 landmarks (Figure 2A, B) were registered for each one of them, using the TpsDig program (Rohlf, 1998).

The landmark configurations were superimposed by generalized Procrustes analysis (a least squares superimposition - Monteiro and Reis, 1999), in order to remove effects of scale, position and orientation from the coordinate data. The aligned coordinates were projected into the shape space of partial warps and uniform components (Bookstein, 1991; Rohlf, 1996), which are shape variables that span a linear space tangent to the curved shape space generated by the Procrustes superimposition (Rohlf, 1999). The partial warps and uniform components are a choice of base for the statistical analysis of shape consistent with the commonly used tools of linear multivariate analysis (Rohlf, 1999). These shape variables describe partial differences between each shape in a sample and a reference configuration (usually a mean shape). The average shapes used as reference configurations are depicted in Figure 2C (females) and Figure 2D (males). In each case, the reference configuration (Rohlf, 1996) used for the calculation of partial warps was the grand mean of all populations. The partial warps and uniform components have been extensively described in the literature (Bookstein, 1991; Rohlf, 1996). To assess the structure of body shape variation among populations, the partial warps and uniform components were combined in a data matrix and submitted to canonical variate analyses (Johnson and Wichern, 1988). Previous analyses of a smaller subset of the same data have shown a significant interaction of shape differences between the lagoon and habitat factors (Neves and Monteiro, 2003). Therefore, the samples were grouped by the interaction term of lagoon and habitat (the collection sites described in the previous section), for an assessment of the direction of major variation among populations in shape space (while statistically controlling for within-group variation). The populations were represented by 95% confidence ellipses (with standard errors of centroids) generated by a parametric bootstrap resampling (Von Zuben et al., 1998), where 1000 replicates of the data matrices were constructed by random sampling, based on the original covariance matrix and the vector of means. Because of a marked sexual dimorphism, separate analyses were performed for males and females. The visualization of shape changes depicted by each canonical vector was obtained by multivariate regression of partial warps and uniform components on the canonical scores (Rohlf et al., 1996).

The aligned coordinates were also used to calculate geometric shape distances (Procrustes chord distances - Rohlf, 1999) between population mean shapes. The Procrustes metric is the base of the shape space used by geometric morphometric methods (Rohlf, 1999; Monteiro et al., 2000) and is defined as the summed squared distances between corresponding points in two configurations, after the least-squares superimposition. These distances serve as an overall measure of shape differentiation.

Evolutionary divergence tests and parameter estimation

The choice of which model to use depends on the effective population size (N_e) and the number of generations (t) passed between ancestors and descendants, or the number of generations as separate populations for divergence rate tests. The constant heritability (CH) model (Lande, 1977) is designed for recently separated populations, which have not yet achieved mutation-drift equilibrium. It is appropriate if t < N_e/5. For populations in equilibrium, one should use the mutation-drift equilibrium model (MDE), which will be appropriate if t > 4N_e (Turelli et al., 1988). In general, one should use the CH model for recently separated populations and the MDE model for populations or species separated for a long time (for a discussion see Spicer, 1993).

The estimation of test parameters may be problematic, particularly of those related to quantitative genetics (heritabilities, mutational variances, effective population sizes, generation times and number of generations passed). A great number of parameter estimates for various characters and organisms can be found in the literature (Spicer, 1993). As a rule, reasonable and conservative estimates should be preferred (Diniz-Filho, 2000). In our case, the easiest-to-determine parameter was t, the number of generations passed since divergence. From geological dating results (Martin et al., 1997), we can infer that the oldest lagoon in our study (Campelo Lagoon) is around 4000 years old. The most recent lagoons (Grussaí and Iquipari) are closer to 1500 years old. Using a life table method, Reznick et al. (1997) estimated the number of generations/year as 1.74 in the guppy (Poecilia reticulata). Based on that estimate, we used a rounded estimate of 2 generations/year. Thus, we could set a maximum limit at 8000 generations passed since divergence for the population of the Campelo Lagoon in relation to the other lagoons. The minimum limit was more difficult to define, because colonization by P. vivipara needed not to have started immediately after the formation of the lagoons (although this is very likely). There is also a possibility that the populations of different lagoons were connected by gene flow until the construction of irrigation channels 50 years ago (100 generations). Fortunately, we could infer the maximum limit with greater confidence, and use it for conservativeness. The next step would be to calculate the effective population size (N_e), that can be estimated from the total population size (N), since the smallest N_e/N ratios range from around 0.5 (Spicer, 1993; Nunney, 1996) to 0.1 (Frankham, 1995). In our case, a standardized variance estimate of 0.42 fry/brood in female fecundity was available (Monteiro and Gomes-Jr, unpublished data). Using the relationship of N_e/N to standardized variance in fecundity published by Nunney (1996), we can calculate a conservative expectation of the N_e/N ratio as 0.3. From preliminary mean density estimates (based on our sampling) of number of fishes by marginal length (50 specimens per meter), we can conservatively estimate the population size in each one of the smaller lagoons (Grussaí and Iquipari) to be around 1,000,000 individuals. Of course, there should be spatial variation of density within lagoons, but the animals are really abundant and omnipresent in the region. Considering the N_e/N ratio to be 0.3, we can calculate the minimum N_e to be 300,000 in each lagoon (this number might be slightly larger in the Campelo and Açu Lagoons, which are larger). Considering that we were studying four lagoons, the estimated pooled N_e would be 1,200,000. According to these estimates of t and N_e, the relationship t/N_e would be 8,000/1,200,000 = 0.0006, smaller than 0.2, which would indicate the use of the constant heritability model for the directional selection test (Turelli et al., 1988).

The tests of evolutionary rates and evolutionary divergence are based on the assumption that, under neutral evolution (genetic drift) models, continuous phenotypic characters have an expected mean at generation t that is approximated by a normal distribution whose variance is proportional to s²h²t/N_e, where s² is the variance of the character, h² is the heritability of the character, and N_e is the effective population size (Lande, 1977; Turelli et al., 1988). In this model, the expected deviation from the mean at t = 0 should increase for larger variances, heritabilities and number of generations passed, but should decrease for large effective population sizes (as expected in genetic drift models). For morphological divergence tests, whenever the observed variation among populations is greater than expected by the neutral model, we can infer that a process of directional selection is taking place. When the observed variation is smaller than expected, we can infer that a process of stabilizing selection is maintaining the similarity among populations (Lande, 1977; Turelli et al., 1988). If the observed deviation is within the expected interval, we accept the null hypothesis of neutral evolutionary processes. The test statistics for evolutionary divergence follow an F-distribution with n - 1 degrees of freedom in the numerator (n = number of lineages or populations studied at once) and infinite degrees of freedom in the denominator. The statistics for the test of directional selection in the constant-heritability model is

where the difference between populations is measured by an among-populations mean square

where the term (t) corresponds to the grand average over all populations, and _i (t) corresponds to the mean of the ith population. The test for stabilizing selection is performed by the inversion of the equation for the F-statistics. Because of the inversion, the degrees of freedom in the numerator are infinite and in the denominator are equal to n - 1 (Lande, 1977; Turelli et al., 1988).

The amount of phenotypic divergence among populations () was measured by the scores of the first canonical axes for the analyses of males and females separately. The reduction of the multidimensional shape space using a variance-maximizing criterion allows for the univariate assessment of vectors of interest in shape space, in our case, the major axes of among-population variation (assuming, of course, that the percentage of among-group variation explained by the first canonical axis is significantly greater than one would expect at random). An advantage of this approach over the conventional multivariate statistical methods (Lande, 1979; Lofsvold, 1988) is that, as we shall see below, we do not need an estimate of the entire genetic covariance matrix (G), but only the inference of a significant additive genetic component of shape variation.

The remaining parameters of the model that have to be estimated are: phenotypic variance of shape (s²) and shape heritability (h²). The estimate of the phenotypic variance of shape was calculated (as the among-population shape variation) from the first canonical variate scores (as the within-population variance). Because of the standardization of within-group shape variation in canonical scores, the within-group variance is equal to 1. Heritability estimates are more difficult to obtain directly, because they represent the proportion of total variance that is caused by additive effects of genes (as opposed to dominance, epistasis and environmental variance) (Lynch and Walsh, 1998). As it depends on the environmental component of variance, the heritability of a character is also a variable parameter among populations in different environments and within populations at different times (Falconer and McKay, 1996). In the literature, the heritability estimates for morphological characters of fish species vary from 0.14 to 0.80, where even the lowest values are statistically significant (Smoker et al., 1994; Taniguchi et al., 1996; Jonasson and Gjedrem, 1997; Choe and Yomozaki, 1998; Nakajima and Taniguchi, 2002). Within Poeciliids, Mousseau and Roff (1987) reported a median heritability of 0.745 for morphological characters in the mosquitofish (Gambusia affinis). In our case, it would be important to define the heritability of the specific vector of among-group differences, for different directions in shape space can have different heritabilities (Monteiro et al., 2002; Klingenberg, 2003).

Because of the uncertainty associated with parameter estimation, the rate tests should not be considered rigorous quantitative tests, but qualitative indications of evolutionary processes (Turelli et al., 1988). In order to assess which parameters are more influential on the test results, we calculated a value of sensitivity (S) for each parameter when varying a standard (conservative) estimate within a range of values corresponding to 25%, 50%, 75%, 125%, 150% and 175% of the standard estimate for each parameter. The sensitivity for a given parameter S_p was calculated as in Dowd (1997):

where DF corresponds to the difference of F-statistics calculated with the variable parameter and the standard estimates of the parameters, F is the value of the F-statistics calculated with the standard estimates, Dp is the difference between the variable parameter and the standard estimate, and p is the standard estimate of the parameter. The parameters used for the sensitivity test were: (the among-group variance on canonical scores for the female data set, standard estimate = 3.32); s² (the within-group variance on canonical scores, standard estimate = 1); h² (the heritability of the shape variable, standard estimate = 0.6); t (the number of generations since divergence, standard estimate = 8,000); N_e (the effective population size, standard estimate = 300,000).

In order to assess the robustness of the results due to variation in the most uncertain parameters for this study, (N_e, h², t), we calculated a range of F-values from the F-statistics equation for a conservative series of parameter estimates (N_e = 1,000 to 100,000; h² = 0.1 to 0.9; t = 8,000 to 20,000) and checked which combinations of estimates would cause the acceptance of the null hypothesis of genetic drift. This second sensitivity analysis provided an in-depth examination of the effects and interactions among specific parameters on the test outcomes. The range for h² comprised a wide range of heritabilities possibly observed in nature; the range for t allowed for the possibility of error in the estimate of the number of generations/year (from our estimate of two up to the unrealistically conservative five generations/year); the range for N_e used was needed in order to find the borderline of non-significant results. First, a sensitivity analysis of two parameters (N_e and h²) was performed, fixing t at 8,000 generations (our best estimate). The results were visualized as a spline surface, where the parameters N_e and h² were designated to the X and Y-axes, respectively, and the F-statistics to the Z-axis. Second, in order to examine the effect of uncertainty on the three parameters jointly, we solved the F-statistics equation, to obtain N_e for the critical values of the F-statistics and a combination of values for h² and t. The results of the analysis using three variable parameters were visualized as a spline surface for critical N_e as a function of t and h².

Comparisons of observed divergence rates in males and females were also performed (considering each one of the seven sites), using the correlation of Squared Procrustes Distances with the number of generations since divergence for different lagoons. This latter matrix was built using geological dating information (Martin et al., 1997), assuming that the time of formation for each lagoon should be considered the maximum amount of time passed since its population diverged from the hypothesized main source of invaders, the Paraíba do Sul river (although for some lagoons the maintenance of gene flow until few years ago is possible). The significance of matrix correlations (shape distances and time) was assessed by a Mantel test (Diniz-Filho, 2000) using 9,999 permutations.

Results

The structure of among-population shape variation was assessed by canonical variate analyses of partial warps and uniform components. Because of the marked sexual dimorphism in body shape, a separate analysis was performed for each gender. The two first canonical variates for the female data set explain 85.44% of the total among-group variation. The ordination of population centroids (Figure 3) shows a considerable difference among the populations. The Açu and Campelo Lagoons appear in extreme positions of the first axis, whereas the Grussaí and Iquipari Lagoon populations (the lagoons geographically closer to each other) are overlapped. A conspicuous pattern related to the second canonical axis is the separation of sand bar populations from interior populations in the Grussaí and Iquipari Lagoons. The shape differences associated with the major axis of among-group variation (grids below the scatterplot in Figure 3) depict mostly differences in midbody height (relatively larger in positive scores, such as the Açu Lagoon populations), the head region (relatively smaller in positive scores), and the position of mouth opening (forwards in populations with positive scores and upwards in populations with negative scores). The two first canonical variates for the male data set explain together 81.78% of the total among-group variation. It is similar to the proportion of variance explained for the females, but for the males the variation is more concentrated on the first axis. The ordination of populations (Figure 4) shows a contrast between the Açu and Campelo Lagoons against the Grussaí and Iquipari Lagoons (similar to the pattern seen on the second axis in the female data set). The separation between the Açu and the Campelo Lagoons appears on the second axis. The morphological differences associated with the major axis of among-group variation are concentrated in the head region. It is mostly a difference in the relative position of the pectoral fin (more anteriorly positioned in populations with positive scores) and the shape of the opercle. The ordination patterns observed confirm the existence of an interaction between lagoon and habitat, and these two factors should not be considered separately. The components of among-population variance on the first canonical axis scores (used as estimates of ) were 3.25 for the female data set and 6.81 for the male data set.

The sensitivity analysis performed clearly showed two classes of parameters with different influences on the results (Table 1). The parameters directly related with the F-statistics (, N_e) are less influential than the parameters inversely related to the F-statistics (h², t, s²). The latter ones also present a non-linear relationship of S and the percentage deviation from the standard estimate. The underestimation of parameters (25, 50 and 75%) had a greater influence on the test results than their overestimation (125, 150 and 175%). Each step of the analysis decreased the calculated value of S by approximately half the previous value.

The application of the constant heritability model clearly indicates a process of directional selection among populations. The critical F-value (one-tailed) for acceptance of the neutral model with an alpha of 0.05, (n - 1) = 6 degrees of freedom in the numerator and infinite degrees of freedom in the denominator is 2.1. This critical value was used to infer the parameter combinations that would cause acceptance of the null model. For the female data set, the analysis of contour plotted spline surfaces (Figure 5) showed that, for the amount of morphological divergence observed, the hypothesis of neutral evolution would be accepted for a heritability of around 0.6 and a maximum effective population size of around 5,000 individuals. To accept the neutral hypothesis, lower heritabilities would also lower the population sizes. For the male data set, the sensitivity analysis indicates even more conspicuous results (because the among-population variance in males is twice as large as in females). Hence, we do not show the graphical results for the males. Considering the possibility of error in the estimation of the number of generations/year, we have also varied t and calculated critical N_e (using the critical value of F) for a range of heritabilities and number of generations. The results for females and males were also very similar, and the resulting test is more dependent on heritability than on number of generations (Figure 6). The gradient observable in the isolines is more aligned with the heritability axis. This means that in the worst case scenario, where the real number of generations since population divergence was 20,000 (five generations/year) and the heritability of the shape differences between habitats was high (0.9), the minimum effective population size required for genetic drift to be accepted was around 30,000. Any smaller number of generations or smaller heritabilities would decrease the critical N_e. Since the minimum effective population size needed to accept the neutral evolution hypothesis was below the minimum estimated effective population size for each lagoon (300,000), we considered that the tests provided strong evidence that body shape diversification was faster than expected by the neutral model of genetic drift.

The comparison of morphometric distances (Procrustes Distances) and number of generations of (presumed) isolation (among the seven sites separately) shows a pattern of increasing shape distances as the time of separation increases (Figure 7). The pattern is not linear and there appears to be more variation for shape distances with long times of separation. There is, however, a significant correlation between shape distances and time (for females, R = 0.741; p = 0.008; for males, R = 0.716; p = 0.005). For all numbers of generations since divergence, the phenotypic divergence was larger for males than for females.

Natural selection and spatial structure of variation

The structure of shape variation among habitats and lagoons shows a clear pattern of spatial differentiation, both for males and females, although the ordinations were not the same for both sexes. A previous study focusing only the Grussaí and Iquipari Lagoons showed that, for females, there was more variation between habitats within the same lagoon than between lagoons, whereas for males the variation between lagoons was larger than between habitats (Neves and Monteiro, 2003). The difference in ordinations also indicates that the sexual dimorphism in this species goes beyond the morphological differences caused by internal fertilization (Constantz, 1989). The major axis of variation among lagoons was not the same for males and females (vector correlation -0.16), which may be considered evidence that different factors influence the shape divergence of the sexes among populations. This result might be associated with differences in behavior and migration patterns. Genetic evidence from other poeciliine species indicates that males have higher mobility among sites than females (Becher and Magurran, 2000; Johnson, 2001). Although there is a clear association of the shape differences among P. vivipara populations with a known environmental gradient on the scale of single lagoons (Neves and Monteiro, 2003), this larger-scale study does not show the same clear pattern. However, the shape differences observed (differences in body height and the relative placement of pectoral, pelvic and dorsal fins) might have consequences in the functional performance during foraging and escape behaviors (Walker, 1997).

The nonlinear pattern in the relationship of shape distances and number of generations apart indicates that the rate of evolutionary divergence among populations of P. vivipara in the region was not constant. A noticeable aspect was that populations within the same lagoon were as different (or more different, in the case of females) as populations separated for at least 3,000 generations. This result is probably associated with a large amount of phenotypic differentiation among habitats within the same lagoon for the females, as observed by Neves and Monteiro (2003). A second interesting observation was that males diverged more than females (their rate of shape divergence was twice as high). This might suggest that the body shape has been more influenced by selection in males than in females. In fact, the males are more conspicuous than the females and might be under stronger predation. Mortality in males is certainly greater, since most sex ratios observed in the samples were 3 females to 1 male.

The combination of geometric methods with evolutionary models was highly informative regarding possible processes influencing the variation of body shape in the populations of P. vivipara. Given that the estimates of uncertain parameters were very conservative, we can indicate directional selection as a probable process responsible for the observed shape differences in populations inhabiting different lagoons and even in different habitats within lagoons. The magnitude of the observed shape differences and the conservative estimates of effective population sizes were too large for genetic drift to explain the observed pattern as a sole cause. Important questions have been raised by our study, and further research should be conducted in order to answer them. Information from genetic studies is needed for a better understanding of the genealogical relationship of the populations, the actual magnitude of genetic variation, the presence and magnitude of gene flow among lagoons and populations among lagoons. This information would provide a framework where the relationship of morphological and neutral genetic divergence could be compared (Diniz-Filho, 2000). Quantitative genetic experiments are also needed, in order to partition the effects of genetic and environmental differences within and among populations (Lynch and Walsh, 1998), and to determine the genetic, rather than the phenotypic, morphological divergence rates (Hendry and Kinnison, 1999). Further research in this area will also require the extension of the data collection to other lagoons, including wider environmental gradients and geographical area. A better documentation of the fish community structure among lagoons would also increase our understanding of larger-scale patterns (Magurran and Phillip, 2001) influencing the evolution of P. vivipara in the region.

Sensitivity analysis and robustness of results

Sensitivity analysis indicates that the most influential parameters of the test (and those which should be estimated with more accuracy) were the ones inversely related to the F-statistics (h², t, s²). Unfortunately, these are also the more difficult to estimate accurately. For these parameters, the choice of conservative values should preferably overestimate the actual parametric values, for high estimates of the parameters will have less influence than low ones. Although the sensitivity (S) calculated was the same for h², t, and s², the influence of each one of these parameters will not be the same for this specific example, because the magnitude of their influence also depends on the range of conservative estimates assumed for each parameter. The parameters with the largest relative scale will ultimately have the greatest influence on the test results. This is one of the reasons why we should consider the outcomes of a range of parameter values, within sensible intervals (Turelli et al., 1988; Spicer, 1993). In our study, the parameters estimated with less certainty were effective population size (N_e) and heritability. The estimates of N_e presented lack a specific method (such as demographic and molecular variation information - Clegg et al., 2002) for their calculation. However, even if we consider the lower expectation of the N_e/N ratio in the literature (Frankham, 1995) as 0.1 (100,000 individuals breeding in each lagoon), we would still be above the minimum effective size for genetic drift to be significant (5,000 individuals). There is also a possibility that the populations have experienced bottlenecks after invasion of the lagoons. However, there is no geological evidence of climatic or geological phenomena that could have caused such events (Martin et al., 1997). It is difficult to set narrow confidence intervals for heritability estimates, because they require larger sample sizes than most studies can provide (Falconer and MacKay, 1996; Lynch and Walsh, 1998). However, given that any value of heritability would cause rejection of the neutral hypothesis (in fact, smaller heritabilities would yield more significant results), the important question would be: Is there a significant genetic component of body shape variation? Given the significance of all heritability estimates for morphological variables in the literature, it is not likely that complex variables such as shape or size will have zero parametric genetic variation (no matter how large the reaction norms might be). The estimate of time of formation of the lagoons was based on geological dating (Martin et al., 1997), and a large error component seems unlikely. The number of generations/year, however, is not known for this species. Because a similar well-known congeneric species (Poecilia reticulata) had reliable, time-table-based estimates of generation times, we expect the maximum number of generations/year to be not distant from two. The Sensitivity analysis showed that, even if the actual number of generations of divergence was much larger, the conclusions would be the same. It is more likely that the actual number of generations is smaller than our conservative estimate. Because of the uncertainty in parameter estimation, the rate tests have been considered a qualitative evidence of natural selection in nature (Spicer, 1993). This is the first work that assesses the effect of parameter estimation error on the results by means of a sensitivity analysis using conservative parameter intervals. An interesting result was the observation of which parameters are more important and how they interact with each other to influence the outcome of the test, for example, smaller heritabilities require smaller effective population sizes for selection to overcome genetic drift. As the number of generations varied on a smaller relative scale (the smaller value in the range differed from the larger value by a factor of 2.5, compared to a factor of 9 for heritability and 10 for Ne), it did not have a considerable influence on the results.

The approach presented here has an advantage over traditional multivariate tests for directional selection (Lofsvold, 1988), because it allows to perform a rate test even if the G (genetic covariance) matrix is not known. The use of phenotypic covariances as surrogates for genetic covariances in multivariate selection tests has been proposed (Cheverud, 1988; Roff, 1995). However, words of caution have been set forth by a number of researchers against that procedure (Phillips et al., 2001; Steppan et al., 2002; Willis et al., 1991), because the genetic and phenotypic covariances can differ by a number of reasons, and more research on that subject is certainly needed before genetic covariances can be replaced by phenotypic covariances with confidence (Lynch and Walsh, 1998). The multivariate selection tests used are based on the fact that, under genetic drift, the G matrix will be modified by a constant proportional to t/N_e (Lande, 1979), and the G matrices of diverging populations or species should remain proportional (Cheetham et al., 1993; Lofsvold, 1988). Recent evidence indicates that genetic drift might change the G matrices by more than a proportionality constant (Phillips et al., 2001), and a theory is still needed to describe the expected evolution of genetic covariances under genetic drift. When facing uncertainty of parameter estimation, it is simpler to perform a sensitivity analysis for a smaller number of parameters, as in the divergence rate tests, than for all genetic variances and covariances (as required by the traditional multivariate tests). It was also shown that, if at least some of the parameters are estimated with certainty, it is possible to conduct a test with sensitivity analysis to consider genetic drift as a possible explanation for observed morphological divergence patterns.

Acknowledgments

References

Becher SA and Magurran AE (2000) Gene flow in Trinidadian guppies. J Fish Biol 56:241-249.        [ Links ]

Bizerril CRSF and Primo PBS (2001) Peixes de Águas Interiores do Estado do Rio de Janeiro. FEMAR-SEMADS, Rio de Janeiro, 417 pp.        [ Links ]

Bookstein FL (1988) Random walk and the biometrics of morphological characters. Evolutionary Biology 23:369-398.        [ Links ]

Bookstein FL (1991) Morphometric Tools for Landmark Data. Geometry and Biology. Cambridge University Press, New York, 999 pp.        [ Links ]

Cheetham AH, Jackson, JBC and Hayek LC (1993) Quantitative genetics of bryozoan phenotypic evolution. I. Rate tests for random change vs. selection in differentiation of living species. Evolution 47:1526-1538.        [ Links ]

Cheverud, JM (1988) A comparison of genetic and phenotypic correlations. Evolution 42:958-969.        [ Links ]

Choe MK and Yomozaki F (1998) Estimation of heritabilities of growth traits, and phenotypic and genetic correlations in juvenile masu salmon Oncorhynchus masou. Fish Sci 64:903-908.        [ Links ]

Clegg SM, Degnan SM, Moritz C, Estoup A, Kikkawa J and Owens IPF (2002) Microevolution in island forms: The roles of drift and directional selection in morphological divergence of a passerine bird. Evolution 56:2090-2099.        [ Links ]

Constantz GD (1989) Reproductive biology of poeciliid fishes. In: Meffe A and Snelson FF (eds) The Ecology and Evolution of Poeciliid Fishes (Poeciliidae). Prentice Hall, New Jersey, pp 33-50.        [ Links ]

Diniz-Filho JAF (2000) Métodos Filogenéticos Comparativos. Holos Editora, Ribeirão Preto, 162 pp.        [ Links ]

Dowd M (1997) On predicting the growth of cultured bivalves. Ecol Modell 104:113-131.        [ Links ]

Falconer DS and Mackay TFC (1996) Introduction to Quantitative Genetics. Longman Group, Harlow Essex, 464 pp.        [ Links ]

Frankham R (1995) Effective population size / adult population size ratios in wildlife - A review. Genet Res 66:95-107.        [ Links ]

Hendry AP and Kinnison MT (1999) The pace of modern life: Measuring rates of contemporary microevolution. Evolution 53:1637-1653.        [ Links ]

Johnson JB (2001) Hierarchical organization of genetic variation in the Costa Rican livebearing fish Brachyrhaphis rhabdophora (Poeciliidae). Biol J Linn Soc Lond 72:519-527.        [ Links ]

Johnson RA and Wichern DW (1988) Applied Multivariate Statistical Analysis. Prentice-Hall, New Jersey, 607 pp.        [ Links ]

Jonasson J and Gjedrem T (1997) Genetic correlation for body weight of Atlantic salmon grilse between fish in sea ranching and land-based farming. Aquaculture 157:205-214.        [ Links ]

Kinnison MT and Hendry AP (2001) The pace of modern life II: From rates of contemporary microevolution to patterns and processes. Genetica 112-113:145-164.        [ Links ]

Klingenberg CP (2003) Quantitative genetics of geometric shape: Heritability and the pitfalls of the univariate approach. Evolution 57:191-195        [ Links ]

Lande R (1977) Statistical tests for natural selection on quantitative characters. Evolution 31:442-444.        [ Links ]

Lande R (1979) Quantitative genetic analysis of multivariate evolution, applied to brain: Body size allometry. Evolution 33:402-416.        [ Links ]

Lemos B, Marroig G and Cerqueira R (2001) Evolutionary rates and stabilizing selection in large bodied opossum skulls (Didelphimorphia, Didelphidae). J Zool Lond 255:181-189.        [ Links ]

Lofsvold D (1988) Quantitative genetics of morphological differentiation in Peromyscus. II. Analysis of selection and drift. Evolution 42:54-67.        [ Links ]

Lynch M (1990) The rate of morphological evolution in mammals from the standpoint of the neutral expectation. Am Nat 136:727-741.        [ Links ]

Lynch M and Walsh B (1998) Genetics and Analysis of Quantitative Traits. Sinauer, Sunderland, 980 pp.        [ Links ]

Magurran AE and Phillip DAT (2001) Evolutionary implications of large scale patterns in the ecology of Trinidadian guppies, Poecilia reticulata. Biol J Linn Soc Lond 73:1-9.        [ Links ]

Martin L, Suguio K, Dominguez JML and Flexor JM (1997) Geologia do Quaternário Costeiro do Litoral Norte do Rio de Janeiro e do Espírito Santo. CPRM, Belo Horizonte, 104 pp.        [ Links ]

Monteiro LR and Reis SF (1999) Princípios de Morfometria Geométrica. Holos Editora, Ribeirão Preto, 188 pp.        [ Links ]

Monteiro LR, Bordin B and Reis SF (2000) Shape distances, shape spaces and the comparison of morphometric methods. Trends Ecol Evol 15:217-220.        [ Links ]

Monteiro LR, Diniz-Filho JAF, Reis SF and Araújo ED (2002) Geometric estimates of heritability in biological shape. Evolution 56:563-572.        [ Links ]

Mousseau TA and Roff DA (1987) Natural selection and the heritability of fitness components. Heredity 59:181-197.        [ Links ]

Nakajima M and Taniguchi N (2002) Genetic control of growth in the guppy (Poecilia reticulata). Aquaculture 204:393-405.        [ Links ]

Neves FM and Monteiro LR (2003) Body shape and size divergence among populations of Poecilia vivipara in coastal lagoons of Southeastern Brazil. J Fish Biol 63:928-941.        [ Links ]

Nunney L (1996) The influence of variation in female fecundity on effective population size. Biol J Linn Soc Lond 59:411-425.        [ Links ]

Phillips PC, Whitlock MC and Fowler K (2001) Inbreeding changes the shape of the genetic covariance matrix in Drosophila melanogaster. Genetics 158:1137-1145.        [ Links ]

Primo PBS, Bizerril CRSF and Sofiatti A (2002) Lagoas do Norte Fluminense - Perfil Ambiental. SEMADS, Rio de Janeiro, 148 pp.        [ Links ]

Reznick DN, Shaw FH, Rodd FH and Shaw RG (1997) Evaluation of the rate of evolution in natural populations of guppies (Poecilia reticulata). Science 275:1934-1937.        [ Links ]

Roff DA (1995) The estimation of genetic correlations from phenotypic correlations: A test of Cheverud's conjecture. Heredity 74:481-490.        [ Links ]

Rohlf FJ (1996) Morphometric spaces, shape components and the effect of linear transformations. In: Marcus LF, Corti M, Loy A, Naylor G and Slice DE (eds) Advances in Morphometrics. Plenum Publishing Corp, New York, pp 117-130.        [ Links ]

Rohlf FJ (1998) TPSDig. Version 1.17. Department of Ecology and Evolution, State University of New York at Stony Brook, Stony Brook.        [ Links ]

Rohlf FJ (1999) Shape statistics: Procrustes superimpositions and tangent spaces. J Class 16:197-223.        [ Links ]

Rohlf FJ and Slice DE (1990) Extensions of the Procrustes method for the optimal superimposition of landmarks. Syst Zool 39:40-59.        [ Links ]

Rohlf FJ, Loy A and Corti M (1996) Morphometric analysis of old world Talpidae (Mammalia, Insetivora) using partial warp scores. Syst Biol 45:344-362.        [ Links ]

Smoker WW, Gharret AJ, Stekoll MS and Joyce JE (1994) Genetic-analysis of size in an anadromous population of pink salmon. Can J Fish Aquat Sci 51:9-15.        [ Links ]

Sofiatti A (1998) Aspectos históricos das lagoas do Norte do Estado do Rio de Janeiro. In: Esteves FA (ed) Ecologia das Lagoas Costeiras do Parque Nacional da Restinga de Jurubatiba e do Município de Macaé, RJ. NUPEM-UFRJ, Rio de Janeiro, pp 1-35.        [ Links ]

Spicer GS (1993) Morphological evolution of the Drosophila virilis species group as assessed by rate tests for natural selection on quantitative caracters. Evolution 47:1240-1254.        [ Links ]

Steppan SJ, Phillips PC and Houle D (2002) Comparative quantitative genetics: Evolution of the G matrix. Trends Ecol Evol 17:320-327.        [ Links ]

Suzuki MS, Ovalle ARC and Pereira EA (1998) Effects of sand bar openings on some limnological variables in a hypertrophic tropical coastal lagoon of Brazil. Hydrobiologia 368:111-122.        [ Links ]

Suzuki MS, Figueiredo RO, Castro SC, Silva CF, Pereira, EA, Silva JA and Aragon GT (2002) Sand bar opening in a coastal lagoon (Iquipari) in the northern region of Rio de Janeiro State: Hydrological and hydrochemical changes. Braz J Biol 62:51-62.        [ Links ]

Taniguchi N, Yamasaki M, Takagi M and Tsujimura A (1996) Genetic and environmental variances of body size and morphological traits in communally reared clonal lines from gynogenetic diploid ayu, Plecoglossus altivelis. Aquaculture 140:333-341.        [ Links ]

Trexler JC (1989) Phenotypic plasticity in poeciliid life histories. In: Meffe A and Snelson FF (eds) The Ecology and Evolution of Poeciliid Fishes (Poeciliidae). Prentice Hall, New Jersey, pp 201-213.        [ Links ]

Turelli M, Gillespie JH and Lande R (1988) Rate tests for selection on quantitative characters during macroevolution and microevolution. Evolution 42:1085-1089.        [ Links ]

Vickery WL and Poulin R (1998) Parasite extinction and colonisation and the evolution of parasite communities: A simulation study. Int J Parasitol 28:727-737.        [ Links ]

Von Zuben FJ, Duarte LC, Stangenhaus G, Pessoa LM and Reis SF (1998) Bootstrap confidence regions for canonical variates: Application to studies of evolutionary differentiation. Biom J 40:327-339.        [ Links ]

Walker JA (1997) Ecological morphology of lacustrine threespined stickleback Gasterosteus aculeatus L. (Gasterosteidae) body shape. Biol J Linn Soc Lond 61:3-50.        [ Links ]

Willis JH, Coyne JA and Kirkpatrick MA (1991) Can one predict the evolution of quantitative characters without genetics? Evolution 45:441-444.        [ Links ]

Correspondence to
Leandro R. Monteiro
Universidade Estadual do Norte Fluminense, Centro de Biociências e Biotecnologia, Laboratório de Ciências Ambientais
Av. Alberto Lamego 2000
28013-600 Campos dos Goytacazes, RJ, Brazil
E-mail: lrmont@uenf.br

Received: May 5, 2004; Accepted: February 15, 2004.