Bipartite life cycle of coral reef fishes promotes increasing shape disparity of the head skeleton during ontogeny: an example from damselfishes (Pomacentridae)

Background Quantitative studies of the variation of disparity during ontogeny exhibited by the radiation of coral reef fishes are lacking. Such studies dealing with the variation of disparity, i.e. the diversity of organic form, over ontogeny could be a first step in detecting evolutionary mechanisms in these fishes. The damselfishes (Pomacentridae) have a bipartite life-cycle, as do the majority of demersal coral reef fishes. During their pelagic dispersion phase, all larvae feed on planktonic prey. On the other hand, juveniles and adults associated with the coral reef environment show a higher diversity of diets. Using geometric morphometrics, we study the ontogenetic dynamic of shape disparity of different head skeletal units (neurocranium, suspensorium and opercle, mandible and premaxilla) in this fish family. We expected that larvae of different species might be relatively similar in shapes. Alternatively, specialization may become notable even in the juvenile and adult phase. Results The disparity levels increase significantly throughout ontogeny for each skeletal unit. At settlement, all larval shapes are already species-specific. Damselfishes show high levels of ontogenetic allometry during their post-settlement growth. The divergence of allometric patterns largely explains the changes in patterns and levels of shape disparity over ontogeny. The rate of shape change and the length of ontogenetic trajectories seem to be less variable among species. We also show that the high levels of shape disparity at the adult stage are correlated to a higher level of ecological and functional diversity in this stage. Conclusion Diversification throughout ontogeny of damselfishes results from the interaction among several developmental novelties enhancing disparity. The bipartite life-cycle of damselfishes exemplifies a case where the variation of environmental factors, i.e. the transition from the more homogeneous oceanic environment to the coral reef offering a wide range of feeding habits, promotes increasing shape disparity of the head skeleton over the ontogeny of fishes.


Background
A primary aim of evolutionary biology is to explain the origin, structure and temporal patterns of phenotypic diversity. Many studies have proposed adaptive explanation for phenotypic variation and have focused on the role of selection in shaping patterns of diversification. Divergent selection is expected to enhance adaptive differences through time (evolutionary and ontogenetic scales). Selection may occur at any life stage and be specific to one life stage, leading to the potential for stage-specific adaptation [1]. Despite its obvious importance, very little is known about the role of ontogeny in adaptive divergence and rigorous studies focusing on the variation of phenotypic diversity over ontogeny in various zoological groups are needed.
Coral reef fishes represent one of the most diverse assemblages of vertebrates, and moreover studies dealing with their morphological diversity are numerous [2][3][4][5][6][7][8][9][10]. These studies have mainly focused on the adult stages and fewer have addressed ecomorphological variation through the ontogeny of reef fishes [e.g. [11,12]]. The majority of coral reef fishes have a complex life-cycle with two distinct phases: (1) a dispersive pelagic larval phase and (2) a sedentary demersal adult phase associated with the coral reef environment. The larval phase ends at reef settlement [13]. In contrast to reefs, the open water environment is thought be more homogenous, especially with regard to the diversity of habitats. Conversely, intrinsic factors of the coral reefs such as high productivity, high spatial and ecological complexity, and high trophic diversity may be involved in promoting the high standing levels of fish diversity [14]. It is thus surprising that quantitative studies on the ontogeny of the radiation of coral reef fishes are lacking. Such studies dealing with the variation of disparity, i.e. the diversity of organic form [15], over the ontogeny could be a first step in detection of evolutionary mechanisms of diversification in these fishes. By studying their ontogeny, via changes in shape with size (i.e. allometry), it is possible to gain a clearer understanding of the timing of selective pressures in coral reef fishes [1].
To our knowledge, very little quantitative data exists on adaptive allometry and on the ontogeny of shape diversification within zoological groups in general and fish clades especially [but see [1,16,17]]. Zelditch et al. [18] have focused on the dynamic of body shape disparity over ontogeny in piranhas and highlighted that the disparity decreases significantly and substantially over ontogeny within this particular clade. Adams and Nistri [19] investigated ontogenetic trajectories of foot morphology in eight species of European plethodontid cave salamander and showed the disparity of adult foot morphology was significantly lower than in juveniles. Two broad categories of factors have been suggested for explaining these ontogenetic patterns of disparity [20]: the external or "ecological" constraints such as the availability of ecological space [21,22], and the internal constraints such as developmental, genetic or functional constraints [18,23,24]. New studies in various taxa selected according to their life-cycle and their ecological diversity are needed. They should give further insights into the roles of external and internal constraints shaping the levels and patterns of disparity over ontogeny.
Including more than 300 species living in coral reef environments, the damselfishes (Pomacentridae) represent one of the most successful radiations of coral reef fishes [3,14,25,26]. Similar to the great majority of coral reef fishes, they have a bipartite life-cycle. In the pelagic environment, all damselfishes larvae feed on planktonic copepods [27]. At the adult stage, three trophic guilds are commonly recognized [25,28]: the pelagic feeders sucking planktonic copepods; the benthic feeders grazing filamentous algae or biting coral polyps; and an intermediate group feeding on planktonic prey, small benthic invertebrates and algae in variable proportions. Consequently, we can consider that damselfishes have a higher trophic diversity at their adult stage, and thus we can expect that functional demands lead to an increase in morphological disparity over ontogeny. Indeed, while all species might be expected to be adapted to catch small planktonic prey during the pelagic larval phase, adult damselfishes would be expected to have specific morphological adaptations allowing optimal prey catching and processing. For example, planktivorous and algivorous species have muscles and skeletal shapes optimizing suction feeding and grazing, respectively [3,4,29].
The evolution of damselfishes exemplifies a case of reticulate adaptive radiation [3] in which morphological divergence at speciation has been associated with the repeated convergence on a limited number of ecotypes: algivory, omnivory and planktivory. Among the Pomacentridae, algivory and omnivory have both arisen seven times, planktivory four times, and feeding on scleractinian coral polyps twice [3,25,26].
The study of ontogenetic allometry can reveal differences in developmental patterns among species underlying their morphological differentiation. How differences in developmental patterns generate species divergence has been successfully addressed in many groups, including trilobites [e.g. [30]], fishes [e.g. [17]], newts [e.g. [31]], reptiles [e.g. [32]], rodents [e.g. [33]], primates [e.g. [34]] and humans [e.g. [35]]. The pattern of ontogenetic shape changes can be described by the allometric trajectory of an organism plus the rate at which it proceeds along the trajectory [36]. Morphological divergence among taxa sharing a common allometric trajectory could be due to changes in the rate or the duration of development (i.e. rate and event heterochrony [37]). Divergence could result from the directional change in the allometric trajectories (i.e. allometric repatterning [37]), coupled or not to alterations in the rate and/or the timing of development. In the present study, we explore whether and how differences in allometric trajectories can shape the ontogenetic dynamic of morphological disparity in damselfishes.
We used geometric morphometrics to examine the patterns of morphological diversification among damselfishes throughout their ontogeny. Four functional units of the head skeleton were separately studied: the neurocranium, the unit «suspensorium and opercle», the mandible and the premaxilla. All these skeletal units are movable elements involved in prey catching [38]. The aims of the present study were to: (1) test the hypothesis that larval morphologies are more similar (less disparate) than the morphologies of juveniles and adults, (2) compare the patterns of shape disparity, i.e. the distribution of shapes and the dimensions along which shapes are most disparate, at different ontogenetic stages, (3) compare ontogenetic trajectories to identify the evolutionary changes in developmental parameters (i.e. allometric patterns, rate of shape changes, amount of shape changes undergone over the course of ontogeny) shaping the level and the pattern of morphological disparity, and (4) explore the relationships between the developmental parameters, phylogenetic data and ecological data.

Sample and data collection
Our study species belong to 6 of the most specious genera of Pomacentridae [25], the phylogenetic relationships of which are illustrated on Figure 1. They represent every trophic group known at the adult stage [25,28]: zooplanktivorous species (Abudefduf sexfasciatus, Chromis sp and Dascyllus trimaculatus), algivorous species (Stegastes nigricans, Chrysiptera brownriggii and Chrysiptera glauca) and species belonging to the intermediate group (Dascyllus aruanus and Pomacentrus pavo). Within the intermediate group, D. aruanus is a mainly carnivorous species feeding on planktonic and benthic copepods [39] when P. pavo feeds on zooplankton and filamentous algae [25]. A total number of 406 specimens of damselfishes were analyzed ( Table 1). The sample of each species represents a complete ontogenetic series from larval (i.e. larvae settling reef) to adult specimens. The samples of two Chromis species, Chromis viridis (n = 15) and Chromis atripectoralis (n = 21), were pooled to build a single ontogeny referred as Chromis sp in this study. These two very close species are named the blue green damselfishes, differing only by the coloration of the pectoral fin base [40]. All juvenile and adult specimens were collected by the authors in the lagoon or on the outer reef slope at Toliara (Mozambique Channel,  [41]. At Moorea Island, crest nets were placed on the reef crest. At Rangiroa atoll, the crest net was positioned in a hoa or channel through which water entered the lagoon. Fish larvae captured in crest nets during the night were collected at dawn. Fishes were preserved in neutralised and buffered 10% formalin for ten days and then transferred to 70% alcohol. All specimens were cleared and stained with alizarin red S [42] in order to display the osseous skeleton. One of us (B.F.) dissected fish heads and collected two-dimensional anatomical landmarks (LMs) on lateral views of the neurocranium, the unit «suspensorium and opercle», the mandible and the premaxilla using a Leica M10 binocular microscope coupled to a camera lucida. Lucida images drawn on sheets were then scanned and the x, y coordinates of LMs were digitized using TpsDig (version 1.40). The configuration of LMs of each skeletal unit used for the analysis is shown on Figure 2. Sixteen homologous LMs were defined on the neurocranium, 12 on the unit «suspensorium and opercle», 12 on the mandible and 6 on the premaxilla, forming four separate data sets. Landmarks are described in [4] despite the fact that the number of LMs used for defining the unit «suspensorium and opercle» is lower than in our previous study of adult damselfishes [4] because some LMs could not be easily and precisely discriminated in larvae. The size of samples differed slightly between structures because some units were damaged before or during the dissection (Table 1).

Geometric morphometrics
We used Procrustes-based geometric morphometrics to study shape variation and dynamics of shape disparity throughout ontogeny [43][44][45][46]. For each skeletal unit, the digitized landmark configurations were subjected to a Generalized Procrustes Analysis (GPA) in order to remove non-shape variation (location, orientation and scale) [47,48]. The "grand mean", i.e. the consensus shape of all specimens and shape variables were then generated as partial warp scores (PWs) including both uniform and non-uniform components [44,49]. The centroid size (CS) of the structure was also computed as the square root of the sum of the squares of the distances from all LMs to their centroid [44].  Figure 1 Phylogenetic relationships among the most specious genera of Pomacentridae. The tree topology drawn according to molecular phylogenetic studies [26] is simplified because all genera are not shown. size and shape [36]. Age information is not available for our specimens, so we use size as ontogenetic scale (i.e. ontogenetic allometry).

Comparing models of ontogenetic allometries
The allometric patterns of shape variation were analyzed using linear multivariate regression of PWs on log-transformed size (ln-CS) [48,[50][51][52]. The null hypothesis that shape develops isometrically was tested in all species using TpsRegr (Version 1.34). The fit of the regression models was evaluated by the explained variance of the model and by a permutation test based on a Generalized Goodall's F-Test with 10,000 permutations.
Differences in allometric models among species were tested by a MANCOVA, testing the null hypothesis of homogeneity of linear allometric models. In these tests, shape variables (PWs) are considered as dependent variables, size (ln-CS) as covariate and species is grouping factor. As suggested by Zelditch et al. [53] and Webster and Zelditch [37], two factors can explain differences in allometric models: (a) the divergence of allometric trajectories and (b) the rate of shape changes. Consequently, we estimated and compared both factors plus two other parameters: shape at settlement (i.e. the starting point of the ontogenetic trajectories) and the length of the ontogenetic trajectory, which is a function of the rate and duration of development. As these parameters are estimated by linear multivariate regressions, they assume a linear relationship between shape and ln-CS. PLD, pelagic larval duration [73]; SL, standard length; SL settlement , standard length at settlement [74]; SL max , maximum standard length [25] in millimeter; n, number of specimens for each structure. Plots of Procrustes distance from the mean larval shape (see below), and the variance explained by the models are used to assess the validity of log-linear models (Goodall's F-test). This is reasonable for our data because the regression of shape on ln-CS explains a large proportion (up to 88%) of shape variance for each structure in every species (see Additional File 1).
1.a. Comparing larval shapes at settlement (starting point of ontogenetic trajectories) The samples of settling larvae are limited (≤ 10 specimens). Consequently, we used a standardized regression residual analysis to estimate and compare larval shapes at settlement [18,54,55]. From the multivariate regression of shape on ln-CS, the non-allometric residual fraction is standardized by Standard6 (IMP-software). «Standardized» data sets of larvae with their respective SL (see SL settlement in Table 1), which are the predicted shapes of the entire population at these sizes, are generated. Multivariate analyses of variance (MANOVA) were performed using partial warp scores, followed by a posteriori tests of pairwise comparisons to test shape differences among species. The statistical significance of the pairwise differences was tested by a resampling-based Ftest and the results of CVA assignment tests were also examined. MANOVAs were performed using Statistica 8.1 (Statsoft 2007). Pairwise F-tests were done in TwoGroup6 (IMP-software) and misclassification rates were given by CVAGen (IMP-software).

1.b. Comparing ontogenetic trajectories in shape space
The differences in trajectories of shape changes were analyzed by comparing the angle between the speciesspecific multivariate regression vectors using VecCom-pare6 (IMP-software). This test is described in detail in [48,53] and was already exemplified in a previous study of allometry in damselfishes [50]. Here, a within-species vector is composed of all regression coefficients of the shape variables (PWs) on the log-transformed CS. The range of angles between such vectors within each species is calculated using a bootstrapping procedure (N = 400). This range was then compared with the angle between the vectors of both species. If the between-species angle exceeds the 95% range of the bootstrapped within-species angles, the between-species angle is considered significantly different, and thus the allometric trajectories are different. Angles were computed pairwise between allometric vectors, and the resulting interspecific dissimilarity relationships (angles) between the allometric trajectories were summarized with scatter plots calculated using nonmetric multidimensional scaling (NMDS). Differences between a reference shape (i.e. larval shape at settlement) and a target shape (i.e. adult shape) can be illustrated with deformation grids (interpolating function "thin plate spline" or TPS [44]). Multivariate regression models of shape on size for each unit were used to show graphical illustrations of ontogenetic allometries in each species. Deformation grids illustrating these ontogenetic allometries were obtained from TpsRegr (Version 1.34).

1.c. Comparing rates of shape changes during ontogeny
The dynamics of shape change (developmental rate), defined as the rate of shape changes per unit of size in this study, was estimated for each species using the Procrustes distance (PD), the metric defining shape dissimilarity in the Kendall shape space [45]. PD between each specimen and the average larvae were regressed on ln-CS in all species separately. The rate of divergence away from the average larval shape was compared among species using the slope of the regressions with Regress6 (IMP-software) ([for detailed explanations on this methodology, see [48,56]). Because the relationships between PD and ln-CS are close to linear (see Additional File 2), these can be statistically compared by ANCOVA.

1.d. Comparing lengths of ontogenetic trajectories
The length of the ontogenetic trajectories is used to compare the net amount of shape change undergone over post-settlement ontogeny in the eight species. These lengths were calculated as the PD between the average larval shape at settlement and at the maximum adult body size. Confidence intervals are placed on these lengths by a bootstrapping procedure detailed in [18]. Because our analyses are based on size-standardized data, so the bootstrapping procedure also takes the uncertainties of the regression into account. The calculations were done using DisparityBox (IMP-software). We focused on the relationships between the studied developmental parameters, phylogenetic data and the diet of each species. Moreover, correlation analyses (Pearson r) were used to test the relationships between some developmental parameters, i.e. the dynamic of shape changes or the length of ontogenetic trajectories, and other species-specific characters, i.e. the pelagic larval duration (PLD) and the size variation undergone during post-settlement growth.

Measuring and analyzing shape disparity
We used the same methodology as [18] to measure morphological disparity (MD) at three ontogenetic stages: (1) at the mean size observed at settlement (Table 1), (2) at a common size of 60 mm SL for every species and (3) at the species-specific maximum adult body size (MAX SL). The intermediate size (60 mm SL; at such a size, some species are already adults and others are always at the juvenile stage) was arbitrarily chosen to estimate the MD when fishes are already settled on coral reef and when size difference among species are eliminated. For these analyses, we used "standardized" data sets, which are the predicted shapes of the entire population (see Table 1) at these stages.
MD was calculated with the following equation: where d j = PD between the mean shape of species j and the grand mean shape (i.e. consensus shape). N = total number of species. The level of disparity among the eight studied species was calculated for each skeletal unit (i.e. disparity at the family level). Then, we calculated the disparity level at the three studied ontogenetic stages within each trophic group. As only two species of the intermediate group were studied (i.e. D. aruanus and P. pavo), we could not calculate the disparity level of this trophic group. Consequently, both species were related to the two other trophic groups. As P. pavo commonly grazes algae, it was grouped with algivorous species. Conversely, D. aruanus was groupd with the zooplanktivrous species. All the calculations were done by DisparityBox (IMP-software), which also uses a bootstrapping procedure to place confidence intervals on this measure. We refer to [18,48] for a detailed explanation of this methodology. In addition to measuring the level of shape disparity, we also examined its pattern of variation, i.e. the distribution of shapes and the dimensions along which shapes are most disparate. This exploration is informative about the dynamic nature of disparity. For example, larvae and adults can have the same level of disparity but they can occupy two different sub-spaces (hyperplanes, i.e. "flat" surfaces of more than two dimensions embedded in higher dimensional space) of the shape space. Biologically, this would mean that the patterns of variation have been re-organized or re-structured to lie along very different pathways. The structures of disparity in different ontogenetic stages can be compared by comparing variance-covariance matrices. First, visual exploration of shape variations in the sub-space defined by the first two relative warps (Relative Warps analysis is equivalent to a Principal Components Analysis [PCA] of shape variables when the scaling factor α = 0, [49]) allow to check if the three ontogenetic stages (i.e. settling larvae, 60 mm SL, MAX SL) share the same hyperplane. Then, we used the program SpaceAngle (IMP-software) to compute the angle between hyperplanes and determine if the between-ontogenetic stages angle is no larger than the within-ontogenetic stages range. This program allows to specifically test if every ontogenetic stage occupied the same subspaces of the morphospace. SpaceAngle (IMP-software) used a method developed by [57] based on PCA. See [48] for detailed explanations on this approach but the angle between two subspaces embedded in a common higher dimensional space can be defined as the angle through which one subspace must be rotated to match the other. If the angle between two subspaces (i.e. between-ontogenetic stages angle) does not significantly differ from zero, we may not reject the null hypothesis (H 0 ) in which specimens of two ontogenetic stages occupy the same subspace. The significance of the angle was determined by a bootstrapping procedure (N = 400) similar to that applied for comparisons of allometric vectors (see above). Relative warps analyses of standardized shape data at the three ontogenetic stages were performed using PCAGen (IMP-software) and the angles between subspaces defined by the first two PCs were calculated by Space-Angle (IMP-software).
Geometric morphometric analyses were performed using computer programs from the TPS series (TpsDig and TpsRegr), written by F.J. Rohlf

Comparing ontogenetic allometries
The post-settlement ontogeny of the cephalic region is highly allometric in damselfishes. Shape variation in each skeletal unit is significantly correlated with logtransformed size (ln-CS) (all p levels of the Generalized Goodall's F test < 0.05, see Additional File 1) and a large proportion (up to 88%) of variation in total shape change during damselfish post-settlement development is explained by ontogenetic allometries.
The eight damselfishes show highly different allometric models. Indeed, interspecific differences in allometric models were highly significant for each structure (MANCOVA, Table 2). MANOVA revealed significant differences among species at the settlement stage (for each skeletal unit, p < 0.0001). All pairwise comparisons revealed statistically significant differences in larval shapes (F-tests, p < 0.0015); thus all species differ significantly from each even after Bonferroni adjustments for 26 comparisons. Moreover, the CVA assignment rates were high for each structure (from 91 to 98% depending to the skeletal unit). The differences in the allometric models of damselfish species could also be attributed to their different rates in shape changes, or to their allometric patterns (i.e. the direction in which the ontogenetic trajectories point in the shape space), or to both. The analysis of the angles between multivariate regression vectors of ontogenetic allometries shows that the great majority of species have a species-specific allometric pattern for each skeletal unit (Table 3). Although belonging to the same genus, each angle between C. glauca and C. brownriggii was low but significantly different from zero (Table 3). On the other hand, the angle between D. aruanus and D. trimaculatus was not significantly different from zero for the neurocranium and the mandibule. The angle between D. aruanus and P. pavo was not significantly different from zero for the unit «suspensorium and opercle». Pairwise comparisons of angles between allometric vectors of each skeletal unit were summarized in the NMDS plots of Figure 3. The plot topology varied according to each skeletal unit although some constants were present. Chromis sp. is highly divergent from the other studied species for the majority of skeletal units, except for the neurocranium. According to the plot topology, a parallelism between the allometric patterns and the diet is present for the neurocranium, the mandible and the premaxilla ( Figure 3A, C, D). Indeed, the first axis allowed good discrimination between mainly planktivorous and mainly algivorous species for the neurocranium and the mandible ( Figure 3A, C). Concerning the premaxilla ( Figure 3D), mainly algivorous species shared similar trajectories showing the highest values along axis 2 although Chromis sp seemed to be an exception. Finally, all large species (i.e. A. sexfasciatus, D. trimaculatus, S. nigricans; Table 1) showed similar allometric patterns for the unit «suspensorium and opercle». Using deformation grids, differences in allometric patterns among species are often visually striking (see Additional File 3). However, at large-scale, some ontogenetic shape changes are common to all species. Adults had shorter and higher neurocranium than larvae while these differences are limited in the two Chrysiptera species and P. pavo. In all other species, the heightening is mainly explained by the growth of the supraoccipital crest (LMs 5-7). Adults showed higher opercle and suspensorium than larvae. The suspensorium of adults is always shorter in its central part. Indeed, the distance between the articulation of the palatin (LM 2) and the two articulations of the hyomandibular (LMs 4, 5) is always shorter in adults. All larvae had a maxillary process of the palatin (LMs 1, 2) rostro-dorsally directed and an articulation quadrate-mandible (LM 8) behind a hypothetical vertical bar passing through the articulation of the palatin. The mandible of all larvae is less high showing a shorter symphisis mandibulae than the adult one. Except in Chromis sp., adults showed a shorter dentigerous process of the premaxilla (LMs 1, 6) than larvae. The ascending process of the premaxilla (LMs 1, 2) is always longer in adults.
For each skeletal unit, the dynamic of shape changes (Procrustes distance Vs ln-CS) significantly differs among species (ANCOVA, Table 4). Figure 4 displays all rates of shape changes for each skeletal unit. Pairwise    Table 1 for abbreviations of the species. comparisons revealed that the developmental rate of the neurocranium was statistically significantly lower in D. aruanus and C. glauca than in the other species. Some significant differences existed among the other species but they cannot be subdivided into mutually exclusive subsets. Consequently we regard them as having the same developmental rate. The dynamic of shape changes varied little for the unit «suspensorium and opercle» although Chromis sp. and S. nigricans had the slowest rates, and D. trimaculatus and P. pavo the highest ones. The highest interspecific differences of developmental rates were observed for the mandible and the premaxilla in which Chromis sp had the significant lowest rate of shape changes for both units. The two Chrysiptera species had the highest developmental rates for the premaxilla. The rate of development of these two species was three times higher than the rate of Chromis sp. Dascyllus aruanus and P. pavo shared a lower developmental rate than both Chrysiptera but a higher one than A. sexfasciatus, S. nigricans and D. trimaculatus. For the mandible, A. sexfasciatus had a rate of shape change higher than Chromis sp but lower than the developmental rates of the other species. In conclusion, the dynamics of shape change were not very different among species and there was little evidence of a   Table 1 for abbreviations of the species. hypothetical relationship between the rate of shape changes and the diet of the species (Figure 4). For the neurocranium, the unit «suspensorium and opercle», the mandible and the premaxilla, all studied species can be subdivided in 2, 3, 3 and 4 groups, respectively. No correlation existed between the rate of development of a skeletal unit and the pelagic larval duration (PLD) of the species (0.52 < p < 0.95 according to the skeletal unit), nor between the developmental rate and the size variation undergone during the post-settlement growth of the species (0.12 < p < 0.83 according to the skeletal unit). The length of ontogenetic trajectories differs among damselfishes for each skeletal unit ( Figure 5). These lengths, expressed in Procrustes distance, are directly proportional to the amount of shape changes undergone by every species during growth. For the neurocranium, A. sexfasciatus showed the longest ontogenetic trajectory and the two Chrysiptera species had the shortest ones. According to the confidence limits (95% CI) obtained by the bootstrapping procedure, Chromis sp, D. aruanus, P. pavo and S. nigricans could be differentiated statistically from each other in pairwise comparisons, but they cannot be subdivided into mutually exclusive subsets by length. So we regard them as having the same length. Dascyllus trimaculatus had an intermediate length of ontogenetic trajectory between these four species and A. sexfasciatus. For the unit «suspensorium and opercle», the lengths of the ontogenetic trajectory of A. sexfasciatus and D. trimaculatus were significantly longer than the other species. On the other hand, Chromis sp, the two Chrysiptera species and S. nigricans underwent the lowest amount of shape changes during growth. The length of the ontogenetic trajectory of Chromis sp was almost twice as short as the other seven studied species. Dascyllus aruanus and P. pavo had an intermediate  Table 1 for abbreviations of the species.
length. Chromis sp had the shortest length of ontogenetic trajectory for the mandible. For both mouthparts, all other species did not show significant differences in the length of ontogenetic trajectories. In conclusion, these comparisons suggest that the length of ontogenetic trajectory is a moderately variable developmental trait. Indeed, all studied species can be subdivided into 4, 3, 2 and 1 group(s) for the neurocranium, the unit «suspensorium and opercle», the mandible and the premaxilla, respectively. The amount of shape changes undergone by each species is not correlated to the size variation observed during the post-settlement growth, except for the neurocranium (neurocranium r = 0.84, p < 0.01; suspensorium and opercle r = 0.66, p = 0.07; mandible r = 0.25, p = 0.55; premaxilla r = -0.16, p = 0.70). No correlation existed between the length of ontogenetic trajectory and the mean PLD of species (0.26 < p < 0.83 according to the skeletal units).
Comparing levels and patterns of shape disparity at three ontogenetic stages The shape disparity among species significantly increased over ontogeny for each skeletal unit ( Figure 6A). The neurocranium and the mandible showed the most spectacular variation. Indeed, the neurocranium and the mandible shapes were respectively five and four times more disparate at the adult stage (MAX SL) than at the larval stage (settlement). The increase of shape disparity is the lowest for the unit «suspensorium and opercle». Calculation of the levels of shape disparity at an intermediate size, common for each studied species, revealed differences in the ontogenetic dynamic of disparity among structures. The level of shape disparity for the neurocranium in fishes sizing 60 mm SL accounts for 64% of the level of disparity at the MAX SL while, at the same intermediate size, the disparity of the premaxilla is 88% of the maximal disparity. The two other skeletal units show intermediate values: the disparity level of the unit «suspensorium and opercle» and the mandible at 60 mm SL is 76% and 78% of the maximal disparity, respectively. The shape disparity significantly increased over ontogeny within both main trophic groups (i.e. mainly zooplanktivorous and mainly algivorous species, Figure 6B-C). Within each trophic group, the increasing of disparity is most important for the neurocranium, the mandible and the premaxilla. Conversely to zooplanktivorous species, the group of mainly algivorous species did not show significant variation of shape disparity of the unit «suspensorium and opercle» over ontogeny. The patterns of shape disparity also varied throughout ontogeny. The distribution of shapes in the plane defined by the first two relative warps shows that larvae and adults do not occupy the same sub-space of shape spaces (Figures 7,8,9,10), indicating a redistribution of variance over ontogeny.
At settlement, the first two relative warps RW1 and RW2 explain 54.7% of the total shape variance of the neurocranium (Figure 7). Having the lowest scores on the RW1 axis, the two Dascyllus species showed a higher neurocranium than the others. The second component (RW2), explaining 11% of the shape variance, highlights small differences in the shape of the supraoccipital crest and in the relative length of the orbital and the post-orbital regions. The variance explained by RW1 and RW2 is higher at the size of 60 mm SL and at the maximum adult body size (~72% of total shape variance) (Figure 7). Three groups can be easily distinguished along the RW1 axis at 60 mm, a first including the two Chrysiptera species and P. pavo, a second grouping A. sexfasciatus, Chromis sp and S. nigricans, and a third including the two Dascyllus species. Pomacentrus pavo and both Chrysiptera have a proportionally longer neurocranium. Abudefduf sexfasciatus and Chromis sp show a larger supraoccipital crest than the other species of the second group (RW2). At the maximum size, S. nigricans shares a more similar shape to the two Dascyllus species allowing the clear distinction of three groups in the shape space defined by RW1 and RW2. The two Chrysiptera species, P. pavo and S. nigricans, have a shorter and ventrally directed vomer (low RW1 and RW2 scores, see LMs 1-2).
The shape variance of the unit «suspensorium and opercle» is relatively low at the settlement stage ( Figure  8). The two Dascyllus species have the lowest RW1 values and P. pavo the highest. Both Dascyllus species show the narrowest opercle and the highest unit «suspensorium and opercle». The second component (RW2) opposes S. nigricans, which shows a horizontal maxillary process of the palatin, to A. sexfasciatus and Chromis sp., in which this process is rostro-dorsally oriented. At 60 mm, the first component (RW1) directly opposes the two Dascyllus species to the two Chrysiptera species and P. pavo; A. sexfasciatus and S. nigricans occupy an intermediate position (Figure 8). The second component (RW2) allows the discrimination of Chromis sp. having a lower suspensorium and a maxillary process rostrodorsally oriented. At the maximum size, S. nigricans has a ventrally bent maxillary process of the palatin (-RW2) (Figure 8).
At the settlement stage, the first two components explain nearly equal amounts of variance of mandible shape (RW1 = 39% and RW2 = 30%) (Figure 9). Only Chromis sp. and S. nigricans show distinct shapes from the other species. Chromis sp. has the highest RW1 scores, revealing a lower mandible with a short retroarticular. Stegastes nigricans differs from the others along RW2, showing a shorter angular in its ventral part. At 60 mm and maximum body size, the first two relative warps seem to express almost the same percentage of total shape variance and the same pattern of shape variation (Figure 9). For both stages, RW1 (~50% of the total shape variance) mainly differentiates species according to the height of the mandible and the length of the symphisis mandibulae, distinguishing A. sexfasciatus and Chromis sp. from the other species. The second relative warp (24% of the shape variance) distinguishes the three herbivorous species (the two Chrysisptera and S. nigricans) from P. pavo and the two Dascyllus by their massive mandible (Figure 9).
The dynamic of the pattern of premaxilla shape disparity seems similar to that of the mandible (Figure 10). At settlement stage, Chromis sp. shows a highly divergent shape along RW1 revealing a thin dentigerous process and an obtuse angle between the ascending and the dentigerous processes. The second relative warp distinguishes A. sexfasciatus and S. nigricans from the two Chrysiptera, the two Dascyllus species and P. pavo. Indeed, the first ones have a shorter dentigerous process and a longer ascending process than the others. The pattern of shape disparity is the same for the intermediate and the adult stages ( Figure 10). Chromis sp. is highly divergent from all other species along RW1; and A. sexfasciatus and the two Chrysiptera species are distinguished from the two Dascyllus species, S. nigricans and P. pavo along RW2.
The calculation of the angles among sub-spaces defined by the first two relative warps corroborates the previous visual explorations (Table 5): larvae and adults occupy different sub-spaces of morphospaces. For the neurocranium, the angle between the sub-spaces of the settlement stage and the maximum adult body size is 86°and thus largely higher than the ranges of the within-hyperplane angles (64°for the settlement stage and 14°for the maximum size). Consequently, these hyperplanes are significantly distinct from each other in the total shape space, revealing that the pattern of shape disparity of larvae at settlement differs from that of adults. The same conclusions can be drawn for (1) all comparisons between the hyperplanes of settlement stage and maximum adult size and for (2) all comparisons between the hyperplanes of settlement stage and intermediate stage (60 mm) for each skeletal unit (Table  5). On the other hand, the results differ for the comparisons between the hyperplanes of 60 mm stage and maximum size stage according to the skeletal unit, highlighting differences in the dynamics of the patterns of shape disparity. For the neurocranium and the unit     «suspensorium and opercle», the analyses show that the shapes at 60 mm and at the maximum adult body size occupy different sub-spaces. Conversely, the angle between these hyperplanes is not significant for the mandible and the premaxilla (Table 5); thus the patterns of shape disparity are similar at these two ontogenetic stages for both skeletal units.

Discussion
The level and the pattern of head skeleton shape disparity significantly vary throughout the post-settlement ontogeny of damselfishes. The disparity level increases over ontogeny both at the family level and within each trophic group ( Figure 6). Adults are more disparate in shapes and occupy different sub-spaces of the total shape spaces than settling larvae, implying reorganization of variance with growth. Although larval shapes are more similar, they are already species-specific at the time of the coral reef settlement probably due to differences in the larval growth and/or the pelagic larval duration (Table 1). Indeed, all species do not have the same age when they settle on the coral reef [58].
According to the study of developmental parameters, the increasing disparity seems to be mainly due to the divergence of allometric patterns. This parameter is clearly not constrained: nearly every species has their own allometric pattern for each skeletal unit (Table 3).
Eight specific allometric patterns exist for the premaxilla and 7 for the other units (Table 6). In the shape space, every species has its own ontogenetic trajectory. Only the allometric vectors of the neurocranium and the mandible in D. aruanus and D. trimaculatus point in the same direction (Table 3). Both species are dissimilar at each stage (i.e. larvae and adult) but have parallel ontogenetic trajectories in the size-shape space, illustrating a case of lateral transposition where dissociation had occurred in an earlier period of reef settlement. The conservation of the same allometric pattern in these two closely related species belonging to the same genus may be explained by genetic similarities. More surprisingly, the allometric patterns of the unit «suspensorium and opercle» in D. aruanus and P. pavo reveal a third case of parallel trajectories. All others show a total divergence of allometric trajectories (i.e. allometric repatterning). Generally speaking, the range of the divergence of allometric patterns (9.9°-53.6°, Table 3) is similar to the results obtained from other geometric morphometric studies of different structures of the head skeleton in vertebrates [31,33,35,59,60], except for the premaxilla which is highly divergent (18.7°-101.3°). The angles among the ontogenetic trajectories seem to be fairly related to the diet (Figure 3). On the other hand, the incongruence between allometric trajectories and phylogeny (Figures 1 & 3) may be linked to the reticulate  The figures represent the diversity of studied allometric parameters (i.e. larval shapes, allometric trajectories, rates of shape changes and lengths of allometric trajectory) among damselfishes. As eight damselfish species were studied, 8 is a maximum stating that each species has a specific developmental parameter.
adaptive radiation (i.e. evolutionary patterns characterized by rapid and repeated shifts between a limited number of ecomorphological states) of damselfishes suggested by Cooper and Westneat [3]. Indeed, our data probably illustrate successive allometric repatternings allowing ecomorphological shifts during the evolution of Pomacentridae. This result implied that the ontogenetic trajectories of the damselfish species are evolutionarily labile. The evolutionary changes seen in ontogenetic trajectories across the phylogeny may reflect the complex interplay between developmental processes and selective pressure that are more intense at various developmental stages (e.g. reef settlement stage and/or adult stage) [19]. The variability of the angles within the same genus (e.g. Chrysiptera and Dascyllus) remains limited in spite of some differences according to the species and the skeletal units. In the two Chrysiptera species, the divergence of allometric trajectories is low for every studied structure. On the other hand, the two Dascyllus species show two identical allometric patterns. Likewise the phylogenetic proximity of the genera Pomacentrus and Chrysiptera suggested by molecular data (Figure 1) is confirmed by the similarity of allometric trajectories. Finally, a functional demand related to the head size could explain the similarity of allometric patterns for the unit «suspensorium and opercle» among the large species: D. trimaculatus, A. sexfasciatus and S. nigricans ( Figure 3). The length of ontogenetic trajectory and the rate of shape change appear to be conservative developmental parameters. Depending on the skeletal unit studied, 1 to 4 values have been highlighted for 8 species (Table 6). Thus some species share the same length of ontogenetic trajectory or the same dynamics of shape change. However the variation of both parameters could be underestimated by our methodology. Indeed, having no age information on the studied specimens, compensatory changes in growth rate (defined as the relationships between age and size) could not be explored. These two developmental parameters are not correlated to the duration of the larval phase or the size variation of post-settlement ontogeny. The reasons for such a low variability of these parameters are not obvious. During settlement, the great majority of damselfishes directly recruit to the adult populations and thus directly use the same habitat of their congeners [61,62] unlike other taxa which undergone one or some habitat changes before their recruitment (e.g. Labridae, Acanthuridae, Serranidae). Among the studied species, only A. sexfasciatus deviates from this general trend: the juveniles settle in micro-atolls of the fringing reef and adults live on the barrier reef [61]. This consistent behavior pattern observed in pomacentrids could partially explain the low inter-specific variability of the rate of development. Both parameters are not related to diet and phylogeny. Intuitively, we hypothesized that planktivorous species at the adult stage, thus having no ontogenetic change of diet, undergo lower amount of shape changes than herbivorous species but that is not verified. The ontogenetic shape changes that are related to a shift of feeding strategies [38,50] and other functional demands might also explain shape changes related to respiration [63] or sound production [64].
Allometry is an important factor in modeling of the head shape during the post-settlement growth of damselfishes. Our results reveal a higher degree of similarity at reef settlement, followed by a significant spreading of species within morphospace over ontogeny. Different factors could explain a closer similarity among settling larvae than among adults, including the internal (e.g. functional and developmental) and external (e.g. ecological) constraints. The latter might seem predominant. Indeed, as all larvae live in the pelagic zone of the ocean [13], the shape similarity seems to be linked to a common environment and very similar diets, exclusively composed of copepod nauplii and adults [27]. On the other hand, the lagoon and the barrier reef offer a multitude of resources and feeding habits. Three major modes of feeding have been identified in Teleost fishes: sucking, biting and ram feeding (for a review see [38]). These three modes may be predicted from the functional design of the damselfish head [3,4] although they possess skulls of intermediate design that will also allow some opportunism in their diet [28,65]. The allometric patterns detailed in the present study strengthen previous data revealing an improvement of suction feeding system during damselfish ontogeny [50]. The main allometric shape changes include: a heightening of the suspensorium and the opercle, an elevation of the supraoccipital crest, a forward displacement of the mandibular articulation with respect to the neurocranium, a mandible becoming proportionally shorter and a lengthening of the ascending process of the premaxilla. Figure 8 summarizes the variation of trophic diversity and morphological disparity throughout ontogeny of damselfishes, highlighted during the present study and others [4,28,29,50,66,67]. As demonstrated in other Teleosts such as Labridae [68,69] and Centrarchidae [70], morphological (shape) diversity may not always be an accurate predictor of biomechanical diversity. An overview of our previous studies dealing with feeding habits [28] and trophic morphology [4,66,67] in adult damselfishes shows that the higher level of shape disparity corresponds to a higher level of ecological and functional diversity in comparison to the larval stage. Nevertheless, future studies will have to specifically test the relationship between the morphological (shape) disparity and the mechanical disparity in damselfishes, as done by Hulsey and Wainwright [71] in labrids. Although larval shapes are already species-specific, some similarities suggest that settling larvae can be considered as ram-suction feeders (see discussion in [50]). Adult damselfishes show either a short neurocranium type (S. nigricans, A. sexfasciatus, Chromis sp, both Dascyllus species) or a more lengthened type (both Chrysiptera species, P. pavo). The main character allowing the differentiation of grazing species is the presence of a ventrally bent vomerian region ( Figure 11). The other skeletal shapes literally vary between the two main functional constraints: suction feeding and grazing/biting ( Figure  11). The pelagic feeders (A. sexfasciatus and D. trimaculatus) and the benthic feeders (both Chrysiptera species and S. nigricans) show skeletal shapes optimizing suction feeding (e.g. high suspensoria and opercles, short mandible) and grazing (e.g. robust mandible, short and massive dentigerous process of the premaxilla, ventrally bent maxillary process of the palatin, broad hyomandibular allowing a large place for the muscle adductor mandibulae) respectively [3,4]. Species belonging to the intermediate group such as P. pavo have intermediate skeletal shapes or their head is a mix of skeletal units typical of suction feeding or grazing species. Among the pelagic feeders, Chromis sp (i.e. Chromis viridis and C. atripectoralis, see Materials & Methods) showed some striking divergences (Figures 3, 7,8,9,10). The shapes and the short lengths of ontogenetic trajectories (Figure 4) are suggestive of functional similarity with larvae. Furthermore our morphological and developmental data corroborate kinematic studies [72] showing that Chromis viridis use ram-suction feeding to capture highly evasive prey such as copepods. In summary, one main morphofunctional group can be recognized at the pelagic larval stage and three at the demersal adult stage (Figure 11).

Conclusion
The present study describes qualitatively and quantitatively the diversification of the head skeleton during the post-settlement growth of damselfishes. The morphological disparity of the head skeleton is higher at the adult stage associated with the coral reef environment in comparison with the settling larval stage. At reef settlement, larval shapes are already species-specific probably due to differences in age. The process of diversification throughout the ontogeny of damselfishes is mainly explained by Morpho-functional groups Figure 11 Schematic representation of trophic diversity and morpho-functional disparity throughout ontogeny of damselfishes. See the text for details. Arrows express gradient of shapes between the two extremes of morpho-functional group.
the divergence of allometric trajectories in shape space. The length of ontogenetic trajectory and the rate of development seem to be more constrained developmental parameters. The bipartite life-cycle of damselfishes exemplifies a case where the environment (i.e. the coral reef) promotes the increasing shape disparity of the head skeleton over ontogeny. Linked to different developmental parameters, we show adult head morphology diverged on various phenotypes across species, suggesting that functional demands for varied feeding strategy are more intense in adults than in larvae. Our study demonstrates that both selection and developmental processes have influenced phenotypic evolution in this fish group.

Additional material
Additional File 1: Fit of regressions of shape versus log-values of centroid size (CS) for each skeletal unit. Tabular data (.xls). Results of regression analyses using TPSRegr (Version 1.34).
Additional File 2: Plots of Procrustes Distance (PD) between each specimen and the average shape of larvae on log-transformed centroid size (ln-CS). Example of Abudefduf sexfasciatus, Dascyllus aruanus and Pomacentrus pavo for each studied skeletal unit: the neurocranium, the unit «suspensorium and opercle», the mandible and the premaxilla.
Additional File 3: Illustrations of ontogenetic shape changes. TPS deformation grids illustrating the ontogenetic shape changes for each skeletal unit in every studied species