Phylogenetic, ecological and biomechanical constraints on larval form: A comparative morphological analysis of barnacle nauplii

Barnacle naupliar larvae are differentiated from other zooplankton by their unique pair of frontal lateral horns, frontal filaments, and a pear-shaped cephalic shield. Their morphology impose constraints on their ecological functions and reflect their evolutionary history. To explore the potential functional basis underlying the similarities and differences in barnacle larval form, we conducted a meta-analysis on the shape of the barnacle nauplii’s cephalic shield and examined its relation to larval size, trophic mode, pelagic larval duration and habitat. Nauplii cephalic shield morphology of 102 species were quantified with normalized elliptic Fourier analysis. Most of the species were distributed around the center of the morphospace but a few extreme groups occupied the periphery: nauplii that were large and lecithotrophic. Subsequent principal component regression analyses showed that larval size was a good predictor of the first shape variations axis (aspect ratio). After allometry adjustment, nauplii from different trophic modes differentiated along the second axis of the major shape variations (relative frontal horn length). Habitat was a poor predictor of variations in naupliar body form, but it could be used to differentiate extreme morphology groups from other nauplii. Our result suggests that size-related biomechanical or developmental constraints and feeding requirements are important in shaping the evolution of the naupliar body form. Within the limitations of these functional constraints, habitat drives the divergence of extreme morphology groups from the majority of species. Our comparative morphometrics analysis demonstrated how variations in larval body form can be quantitatively linked to the functional needs that constrain or drive their diversity, and inform further empirical experiments on larval functional morphology.


Introduction
The study of larval morphology has important phylogenetic and ecological implications. For example, similarities in larval characters between taxonomic groups have been used to support PLOS ONE | https://doi.org/10.1371/journal.pone.0206973 November 8, 2018 1 / 16 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 common ancestry at higher taxonomic levels, e.g. in crabs [1]. Nauplius larva was historically used to assign barnacles (Cirripedia) as crustacea; and larval characters, such as the presence of frontal horns, was used to support the monophyly of the barnacles [2]. Despite the essential role of barnacle larvae in systematics and population connectivity, little is known about the functional implications of the presence or absence of larval morphological variations between different species. Functional constraints, the selective pressure exerted to perform various essential functions within the limits of development and biomechanics, are important forces that shape morphological evolution in marine invertebrate larvae [3]. Attempts have been made to connect overall larval morphology with ecological functions in barnacle larvae. For instance, the streamlined fusiform shape of cyprids was hypothesized to be an optimization for efficient swimming [4,5] and the lack of streamlining in the shape of naupliar cephalic shield was hypothesized to help increase feeding efficiency [6].
While common functional needs could reduce the morphological variations between species, adaptation to specific habitats could drive larval morphology to diverge. For example, the exaggerated morphology of Lepadidae nauplii with very long frontal horns and tail processes is hypothesized to be an adaptation to the oligotrophic pelagic habitat [7]. However, unlike adult barnacles for which adaptive significance of morphological variations has been studied [8,9,10], little is known about if the same constraints and driving forces apply to their planktonic nauplii. A comprehensive, comparative study is needed to first provide a summary of morphological variations. Such data would in turn inform future experiments on functional morphology of barnacle larvae, putting the hypothetical relationships between larval form and functions to a test.
The inter-specific morphological variation in barnacle nauplius larvae has been reviewed [11,12,13]. The goal of these studies were often classification, and therefore they, only provide very detailed descriptions of morphological features of species from a single (or a few) target taxonomic group(s) and are often limited to a single locale, probably due to limited sample availability and/ or difficulty in preparing numerous detailed descriptions. Not only are comparisons rarely made between barnacle species from different families and habitats, variations in the overall larval body form-i.e. the shape of cephalic shield-have not been addressed. Landmark-based geometric morphometrics tools have been applied to study of marine invertebrate larvae [14] to quantitatively summarize larval shape changes. Such landmark-based tools typically require well-defined morphological characters to serve as landmarks. In the case of the cephalic shields of larval barnacles, where landmarks are hard to define, outline analysis can be applied [15].
We present a quantitative analysis on the shape diversity of cephalic shield of barnacle nauplius II larvae of 102 different species. Using outline analysis, we identified the major directions of shape change for the cephalic shield of nauplii and tested for differences in habitats or biomechanical requirements to perform the ecological functions of swimming and feeding, i.e. if size and trophic modes account for the observed variations.

Data acquisition
Outlines of 102 species of barnacle nauplius larvae were analyzed. The majority of data were sourced from published descriptions (89 species from 63 publications) and the remaining 13 species were obtained from larvae reared in our laboratory (Table A in S1 Text). For data from published descriptions, digitization was performed on the outline drawings presented or, in a few cases, on the photos of larvae (n = 4). Publications without a scale bar for the drawings or photos, or without upright dorsal/ ventral view, were later excluded from analysis. For reared larvae, digitization was performed on photos of larvae preserved in 30% ethanol.
Stage II larvae were chosen as the focus of this study for two reasons: 1) they sufficiently captured the among-species variations (none of the stages differed in their shape disparity, see Fig A in S1 Text); and 2) data on stage II larvae are more readily available than those on other stages in both the literature and our collection due to the difficulty of rearing to later stages. For species with differentiation of nauplius sex, female nauplii were used.
Outlines of the cephalic shield of stage II nauplii excluding the appendages were traced manually on each drawing or photo. Tracing was performed in ImageJ [16] using polygon selection followed by spline fitting for smoothing. Outlines were manually adjusted to fit the original drawing as closely as possible, and were later saved as binary masks. In addition to nauplii outlines, linear measurements of larval length, larval width, and frontal horn length (recorded as mean of the pair) were also measured. These measurements were either directly retrieved from the literature or measured from the outline drawings in ImageJ if they were not available from the source literature. Mean values of individuals within each species were used for later analysis.

Outline processing and elliptic Fourier analysis
A standardized number of 200 co-ordinates from the binary mask of the larvae were selected for subsequent analysis. Starting from the tip of the dorsal thoracic spine (or the most posterior point), points were sampled along the outline at regular intervals ( Fig 1A).
Shape analyses were performed with normalized elliptic Fourier analyses (EFA; [17]). After normalization, the outlines were invariant to size, translation and rotation. EFA was computed according to Claude [18]. Coordinates of the outlines were described by series of harmonics functions in EFA and the coefficients of the harmonics were used in subsequent statistical analyses. The number of harmonics required to describe the outline sufficiently was usually much less than the number of points sampled, such that minute variations described by higher order harmonics were excluded to reduce data complexity. The number of harmonics used was determined by the number needed to explain > 99% of the shape data and visually inspected for a close fit between original and reconstructed outlines ( Fig 1B).

Variations in nauplii shapes
Principal component analysis (PCA) was used to explore shape variation in the nauplii. PCA was performed on the variance-covariance matrix of the shape data (i.e. coefficients of harmonics computed from the first 15 harmonics of EFA excluding the first three coefficients of the first harmonics which contain zero values after normalization). To identify the major shape changes pattern, shape changes along PC axes were visualized. To compute the hypothetical outlines on principal components, coefficients of harmonics were first calculated from the product of the PC score and its corresponding eigenvector. PC scores of mean ± 2SD on the first two principal axes were used for this calculation. Using inverse Fourier transformation, the hypothetical shapes were reconstructed from the calculated coefficients. Shape change was visualized with thin-plate spline deformation grids, shown as deformation from the overall mean outline.
Direct visualization suggested possible change in aspect ratio and relative frontal horn length (ratio of frontal horn to larval length) along the PC axes. The correlation between PC axes and changes in linear measurements of aspect ratio and relative frontal horn length were also tested. Vectors of coefficients from regression of shape data on aspect ratio and relative frontal horn length were first calculated. Then correlation (R) between the regression vectors and eigenvectors of PCA was calculated using the dot product of the vectors normalized by their magnitudes. We used the absolute value of R since the direction of PC axes are arbitrary, which give R value range from 0 to 1, with 0 and 1 representing orthogonal vectors and parallel vectors, respectively. Significance of R was determined by permutation test to test the null hypothesis of no correlation between the vectors. R was also expressed as an angle between the vectors with θ = arccos(R).

Relationship between habitats, PLD and nauplii shape variations
To explore whether adult and larval habitats are important factors in determining the shape of nauplii, information on adult and larval habitat and pelagic larval duration (PLD) were incorporated into the PCA plot. Adult barnacle habitats were categorized by the type of substrate or the locations in which these substrata are found. The adult habitats were categorized as borer (Acrothoracican barnacles, which bore into calcareous substratum and live in its own created burrows), parasitic (Rhizocephalan barnacles parasitic on decapod crustaceans), deep-sea (species from hydrothermal vents and deeper waters (>200m)), inter/sub-tidal including the coastal species, epibiotic (species that live on crab carapace and crab gills) and pelagic (the genus Lepas, which live on floating objects). Given the dual-phasic life history of barnacles, habitats of the larvae could differ significantly from that of the adult. Larval habitat (coastal water, coral reef, deep sea and open ocean) and length of planktonic dispersal (PLD: < 5 days, 5-30 days and > 30 days) were also tested for their impact on nauplii shape.

Effect of allometry and trophic modes on nauplii shape variations
Allometry (shape changes associated with size) and trophic mode (planktotrophy vs. lecithotrophy) were assessed as the factors associated with the major directions of nauplii shape change (PC1 and PC2). Allometry was assessed by regressing shape data on size (larval length). The one dimensional shape score proposed by Drake & Klingenberg [19] was calculated to check the fit of the regression between multivariate shape data and size because it is a more intuitive way to summarize the multivariate regression model. The shape score was calculated by projecting the shape data onto a line parallel to the regression vector with the equation s = yβ T � (ββ T ) −0.5 , with the regression model formulated as y = βx + ε (where y is the shape data; β is the vector of regression coefficients; β T is transposed β; x is the vector of size and ε is the vector of residuals). To visualize the predicted shape change along the size gradient, shapes predicted by the regression model were reconstructed from the predicted values of harmonic coefficients.
To check for the association between size and PC1, correlation between the allometric regression coefficients and PC1 was tested. Similarly, correlation between regression coefficients of shape against trophic mode and PC2 was tested. To tease apart the effect of allometry and trophic mode, 'allometry free' shapes, which represent shape variation not explained by size variation, were computed with the residuals from the regression model. Differences in shape between trophic modes were later compared with the allometry free shapes by 1) visualizing the mean difference between trophic groups, 2) projection onto the original PCA plot and 3) statistical testing on the significance of the shape difference using trophic mode as a predictor. We used bootstrapping for statistical estimation between trophic modes (mean shapes, relative frontal horn length). Permutation tests were performed to test the significance of regression of shape on size and shape difference between trophic modes using the geomorph R package [20]. For bootstrapping, equal numbers of species were sampled for both trophic modes in each bootstrap sample (using the number of the smaller group, i.e. lecithotrophy) in an effort to correct for the unequal number of group members. Case re-sampling was used in both bootstrapping and permutation tests for consistency.

Role of phylogeny
The phylogenetic signal of the nauplii was evaluated to examine whether barnacle nauplii shape is a function of phylogenetic relatedness. A published phylogeny was taken directly from Pérez-Losada [21], which covered a subset (n = 36) of the species analyzed for their shapes. Tree was digitized from the paper with TreeSnatcher plus [22]. Phylogenetic signal was estimated with generalized version of Blomberg's K statistics [23] for multivariate shape data (K mult ) [24]. Calculation of K statistics and assessment of significance by permutation test were carried out with the geomorph R package. K statistics estimate phylogenetic signal of the trait relative to what is expected under neutral evolution with the given tree (Brownian motion model), where K = 1 indicates that the trait's evolution matches the model, K < 1 indicates that the trait are less similar among relatives than expected and K > 1 indicates that trait are more similar among relatives than expected. For visualization, 'phylomorphospace' was presented by overlaying phylogeny on the morphospace (PCA plot) of nauplii shape. If nauplii morphology carries a significant phylogenetic signal, a non-random branching pattern on the morphospace is expected. Positions of the phylogeny's internal nodes on the morphospace were estimated by ancestral state reconstruction using the maximum likelihood method with the phytools R package [25]. In addition, we also used the phylogeny to account for non-independence among species for analyses on role of size and trophic mode with the smaller dataset where phylogeny is available. The inclusion of phylogenetic covariance in the analyses enabled us to check whether the same conclusion could be reached after accounting for the effect of common ancestry. Phylogenetic ANOVA/regression was carried out with the geomorph R package.

Nauplii shape variations summarized by aspect ratio and relative frontal horn length
Principal component analysis summarizes the overall shape variations of stage II barnacle nauplii (Fig 2). Most of the shape variations (79% in total) can be explained by the first two principal components (PCs). Most species are scattered around the centroid of this "morphospace", with a small proportion of species with extreme morphologies occupying the periphery (Fig 2A). Direct visualization of changes along the PCs shows the major directions of shape changes (Fig 2B). For PC1, nauplius shape is more rounded in the positive direction while the shape is more elongated and slender in the negative direction. For PC2, the positive direction is characterized by the long frontal horns while the negative direction shows shape with short frontal horns and rounded body with very short tail processes. Some areas in the morphospace are not occupied, i.e. the positive extreme of PC2 at either extremes of PC1. These areas correspond to 1) slender shape, long tail processes with long frontal horns and 2) rounded shape, short tail processes and long frontal horns. Observations from direct visualization prompted us to associate PC1 with change in aspect ratio and PC2 with change in relative frontal horn length. Measurements of aspect ratio and relative frontal horn length overlaid on the plot of PC scores shows trends that agree with the observations from direct visualization (Fig 2A). Permutation tests confirm that PC1 significantly correlates with aspect ratio and PC2 significantly correlates with relative frontal horn length (Table 1).

Overall shape variation scales with size
Larval length is a good predictor for aspect ratio (linear regression, R 2 = 60.27%, Fig B in S1 Text). Multivariate regression of shape against size shows that larval length explains a large proportion of total variation (R 2 = 30.7%, Table 2). This relationship indicates that there is a strong allometric effect on nauplii shape variation. This conclusion holds when tested with a smaller dataset for which phylogenetic non-independence was considered ( Table B in S1 Text). Apermutation test on the correlation between the coefficients of this allometric regression model and PC1 is significant (Table 1), suggesting allometry as a good explaining factor for shape changes along PC1. The regression fits well linearly with no obvious problem except three possible outliers (Balanodytes habei, Scalpellum scalpellum and Tetraclita rufotincta, Fig  3B; regressions were also fitted with data with outliers removed, which improved the fit to a R 2 of 36%, see Table C and Fig C in S1 Text). Notably, the large nauplii from Poecilasmatidae and Lepadidae, which are both >90 th percentile in terms of larval length, are well separated from other nauplii along the allometric shape score axis, suggesting that allometry alone is enough to differentiate them from other nauplii. Visualization of shape change along the size gradient as predicted by the regression model shows a trend of increasing aspect ratio, agreeing with the conclusion that allometry is associated with the shape changes along PC1 (Fig 3C).

Relative frontal horn length differed between trophic modes
Nauplii from different trophic mode are separated along PC2 for both shapes before and after allometry correction, suggesting that: 1) trophic mode is related to shape changes along PC2, i.e. relative frontal horn length and 2) shape differentiation between the trophic modes is independent of allometry (Fig 4A; Fig D in S1 Text for adjustments with alternative methods). Permutation tests show that shape distinction between trophic modes is correlated with PC2 (Table 1) and that there is significant shape difference between planktotrophic and lecithotrophic nauplii ( Table 2). This conclusion holds when tested with a smaller dataset for which phylogenetic non-independence is considered (Table B in S1 Text). However, Lepadidae are still separated from other nauplii along PC2. Visualization of mean shape difference between the trophic groups shows that planktotrophic larvae are characterized by longer frontal horns, a more slender body and a longer dorsal shield spine. Relative frontal horn length is longer in Table 1

Phylogeny and habitat account for extreme morphologies
The phylomorphospace plot (i.e., overlaying the phylogenetic tree on the PCA plot) shows a random branching pattern around the center of the morphospace with lots of branch Comparative morphometrics of barnacle nauplii crisscrossing, while the taxa at the periphery (Poecilasmatidae, Lepadidae and Rhizocephala) follow a non-random pattern (Fig 5). Direct comparison of the phylogenetic tree with the clustering tree of nauplii shape agrees with this observation that most species have morphologies that are not tightly coupled with phylogeny ( Fig E in S1 Text). In terms of the overall outline, nauplii shapes contain a significant but weak phylogenetic signal (K mult = 0.42 p = 0.001).
When the shape attributes of aspect ratio and relative frontal horn length were tested individually, aspect ratio was found to be phylogenetically conserved (K = 1.58, p = 0.001) while the relative frontal horn length has weak phylogenetic signal (K = 0.40, p = 0.02). Mapping of aspect ratio and relative frontal horn length on phylogeny is presented in Fig F in the S1 Text. Habitat signal analysis yielded results similar to that of phylogenetic signal analysis, i.e. neither habitat nor phylogeny alone is sufficient to explain shape variation of nauplii, and this works especially poorly for the majority of species at the center of the morphospace (Fig 5A). The center of the morphospace, where most the species are located, is dominated by species associated with intertidal, subtidal and coral reef habitats. Most species with extreme morphology from pelagic, epibiotic and parasitic habitats were readily identified as being from the family Lepadidae, the family Poecilasmatidae, and the superorder Rhizocephala, respectively.

Outline analysis of larvae
We used the geometric morphometrics method to provide an updated summary and analysis on the morphological variations of the barnacle naupliar cephalic shield. Outline analysis is Enlarged symbols represent the group means. Lines were drawn to connect the data points before and after adjustment for allometry, and they are parallel to the vector of size gradient (Note: the angle of these lines to PC1 is smaller than the calculated angle in Table 1  sometimes referred to as the "landmark-free" morphometrics method, and is potentially useful for comparisons between taxa involving the gain or loss of character [15] (but also see [26]). For example, it is difficult to compare barnacle nauplii with other crustacean nauplii without frontal horns in a morphometrics study (see example of Ascothoracida in Fig 6). Outline analysis also allows the inclusion of more detailed shape variations, which may increase the effect sizes on shape differences [27]. However, the homology-free property of outline analysis also poses challenges to properly aligning outlines and to interpreting observed variations [28]. Another limitation of this analysis approach is that barnacle nauplii outlines are particularly susceptible to the problem of the "Pinocchio effect" [29] because of the large variations in the length of tail processes. Despite these limitations, when coupled with visualization of shape changes, outline analysis can provide an intuitive way to identify important patterns and summarize morphological changes for large numbers of species. Furthermore, output from these analyses can be used to inform computational fluid dynamics or dynamically scaled physical models. For instance, mean shapes can be used as morphology representative for all species of interest, and predicted shapes resulting from a regression model can be used for hypothesis testing concerning specific factors of interest.

Functional constraints imposed by swimming and feeding
Our results suggest that adaptation to adult habitats alone does not account for overall nauplii shape variation. Allometry and trophic mode, however, provide better explanations for the overall nauplii shape variations. In many animals, aspect ratio is a function of size and dominates the major direction of shape variations in various anatomical structures: e.g. fish body form [30], bird beaks [31] and mammal skulls [32]. Such a change in aspect ratio as a function of size is due to a combination of developmental and functional constraints [33]. A shared developmental program restricts the range of variations. However, within such developmental constraints, species can easily "explore" this specified range, possibly facilitating a rapid response to selection [34]. Birds' beak illustrated how functional constraints limit aspect ratio as it is correlated with the load the beak can take without risking fracture [35]. For barnacle nauplii, little is known about the developmental constraint on body form, but there is likely a functional constraint related to swimming. Hydrodynamic drag increases with larval size, and a higher aspect ratio can lower drag [36]. To swim, the propulsive force has to be higher than the drag, thus species with larger nauplii may be selected to have a higher aspect ratio.
Lecithotrophy repeatedly evolved within many phyla of marine invertebrates and is often accompanied by similar simplification of morphology, showing a high degree of convergence [36]. Regardless of the superorder to which the lecithotrophic nauplii belong to (Rhizocephala, Thoracica or Acrothoracica) or the habitat they originate from (deep sea, coral reef associated, parasitic or intertidal), they all occupy the same space morphospace, suggesting convergence (definition sensu Stayton [37]) of nauplii shape. Reduced body form in lecithotrophic marine  Fig 2) incorporated with information of adult habitats (Note: borer refers to Acrothoracican barnacles, which bore on calcareous substrate and live in its own created burrows; deep sea habitat refers to hydrothermal vents and open ocean >200m deep). Nauplii of most species have similar morphology and are concentrated around the center (blue shade,~77%) while species at the periphery are more correlated with certain habitats or taxonomic groups. Peripheral species not belonging to these taxonomic groups are annotated (Bh = Balanodytes habei; Cf = Chthamalus fragilis; Ll = Leucolepas longa, Ob = Octomeris brunnea; Ss = Scalpellum scalpellum; Tr = Tetraclita rufotincta). (B) Relationship between larval habitat and nauplii shapes. Larval habitat differs from adult habitat in that there is less emphasis on substrate type. (C) Relationship between pelagic larval duration (PLD) and nauplii shape. PLD was categorized into extreme and general groups, see Fig H in S1 Text. Data for species reared at low temperature were highlighted, as low rearing temperature may greatly prolong rearing duration. https://doi.org/10.1371/journal.pone.0206973.g005 Comparative morphometrics of barnacle nauplii invertebrate larvae has long been described [36], but not quantified. Our approach quantitatively shows that lecithotrophic nauplii have a reduced body form characterized by shorter frontal horns and tail processes, and a more rounded body (Fig 4).
In addition to the difference in relative frontal horn length between trophic modes, considerable variations were also observed between planktotrophic nauplii. In particular, Lepadidae nauplii have very long frontal horns. One plausible hypothesis is that frontal horns help increase hydrodynamic drags, which increase feeding efficiency of small planktons by acting Examples of naupliar cephalic shield outlines from Cirripedia (barnacles) categorized by their trophic mode. Larvae from Ascothoracida, the Thecostracan sister group to Cirripedia, were also included (outlines traced from photos in [46]). Frontal horns are only present in barnacle nauplii and are absent in other crustacean nauplii. Relative frontal horn length is likely related to feeding. https://doi.org/10.1371/journal.pone.0206973.g006 Comparative morphometrics of barnacle nauplii as partial tether [38]. However, the role of frontal horns in feeding has been contested as this character is present in both planktotrophic and lecithotrophic larvae [2]. If the transition to lecithotrophy occurred faster than the accompanied changes in morphology [39], rudimentary structures used for feeding in planktotrophic larvae could have been retained [40]. However, some taxa with lecithotrophic nauplii, such as the Rhizocephalan barnacles, diverged long time ago (625-596 mya [41]) and were likely to have ample time to lose the structure. Alternatively, frontal horns could serve other function(s) not unrelated to feeding that are retained in lecithotrophic nauplii. Frontal horns have gland cells, which are retained in lecithotrophic nauplii and also in the non-feeding cyprids [42]. The unknown function of the gland cells might have prevented the loss of the associated physical structure of frontal horns, perhaps due to linkages in developmental program. Interestingly, the presence of long frontal horns may also disrupt flow around larval bodies creating small eddies which could contribute towards chemosensing. The hydromechanical function of this unique morphological feature warrants further investigations.
Given both developmental and functional constraints, observed morphologies should reflect the outcome of functional trade-offs. Nauplii with very long frontal horns together with a very high or low aspect ratio are not observed, possibly because swimming ability could not be maintained with such morphology. Functional trade-offs in nauplii shape may have caused the clustering at the center of the morphospace with majority of species displaying a 'generalized', possibly 'compromised' shape [7]. Morphospace pattern with a species rich clade packing the center of the morphospace has been reported for other animals such as fishes and bivalves [43,44], and stabilizing selection due to functional constraint has been proposed as one of the factors driving such a pattern [45].

Extreme morphology and the role of habitat
Habitats of the sessile adults have a less direct influence on barnacle planktonic larvae; therefore, it is not surprising that barnacle naupliar morphology generally does not change with adult habitats. Our results partially support the idea that larval morphology is not tightly related to adult habitat as is often assumed [47] and that the inter-specific morphological diversity is achieved in post larval stages [48].
Larvae with extreme morphology-i.e. those occupying the periphery of the morphospace-are from the family Lepadidae, the family Poecilasmatidae, and the superorder Rhizocephala. The observed phylogenetic signal is likely a reflection of habitat-associated selective pressure exerted over the species at some part of their evolutionary history.
For Lepadidae species, their nauplii possess complex setation and large labrum to capture sparse food particles in the open ocean [6,7]. This morphological change is distinct from other groups for which feeding was lost and larval morphology became "simplified" (see discussion above): e.g. Leucolepas longa and Neoverucca sp. found in deep-sea hydrothermal vents and coastal Tetraclita rufotincta which reproduce during the food-limited winter.
For nauplii of Poecilasmatidae and Rhizocephala, the substrate type of their adults (epibiotic and parasitic) might play a more important role in shaping their overall form. Scarcity of settlement substrate could select for extreme pelagic larval duration: long larval duration (Poecilasmatidae) may enhance chances for opportunistic encounter with substrate, while short larval duration (Rhizocephala) may enhance chances for settling in favorable substrate near the parent. We observed a positive relationship between extremely long PLD and extremely large larval size (Fig 5C; Fig G in S1 Text), which could suggest that increased larval size favors long distance dispersal and small size favors rapid development. This association could also explain why other species from epibiotic habitat, such as Chelonibia testudinaria and Coronula diadema, do not display extreme shape as they do not possess extreme larval duration.
As evolution progresses, yesterday's adaptation may become the present day's constraint [49]. Given that morphology is confounded by phylogenetic relatedness, we observed high phylogenetic signal in species with extreme morphology. The scenario discussed here only represents one of the many possibilities, and future studies on the evolution of barnacle larval morphology will benefit from modeling studies that explicitly evaluate all possible evolutionary models and clarify the relations among the inter-related factors of habitat, phylogeny, development and functional constraints.

Conclusions
Of the 102 species surveyed, the overall morphology of barnacle naupliar larvae shared a high degree of similarity suggesting size-related biomechanical or developmental constraints, such that feeding and other functional requirements are important in shaping the evolution of larval form. Larvae from family Lepadidae, family Poecilasmatidae, and superorder Rhizocephala, however, have extreme morphologies with increased aspect ratio and/or frontal horn length. These extreme characters are likely associated with their larval ecology: scarcity of food or settlement habitats. Quantitatively summarizing variations in larval morphology and testing their link to phylogeny and ecology help better inform evolutionary constraints and understand functional morphology of an important life stage.