Most material used in this study formed part of a previous taxonomical investigation in which Ptilella and Pennatula were compared (García-Cárdenas et al. 2019), and complete information about the surveys, the distribution area of the taxa considered and the data collected is contained therein. A preliminary analysis to ensure the existence of variability within and between colonies at polyp and sclerite measurements level was carried out using eight colonies of Ptilella grayi. These colonies were collected from several surveys in the northeast Atlantic during the period 2010-2014 (see Table 1), with a bathymetry range of 179 to 261 m depth. All Ptilella grayi colonies were collected using a demersal fish trawl with both the cod-end and the full body of the net being thoroughly examined for specimens after each deployment. Total length of preserved colonies ranged from 254 to 572 cm (see García-Cárdenas et al. 2019 for more details).

For the subindividual variability analyses of Ptilella and Pennatula, we selected three colonies of each of the following species: Ptilella grayi, Ptilella grandis, Ptilella inflata, Pennatula phosphorea, Pennatula rubra and Pennatula sp. (Table 1). Part of this material was collected over different geographical areas during various benthic surveys: Antarctica (ANT XVII/3, BIOROSS), the northeast Atlantic (Scotia cruises, INDEMARES Chica), the southeast Atlantic (BENGUELA VIII) and the Mediterranean (Bahía de Algeciras project).

During the expeditions, the specimens were sorted, labelled and fixed in buffered formalin (5% in seawater). After the fixation period, colonies were preserved in 70% ethanol. The terminology used here follows mainly Bayer et al. (1983). The total length of the colonies was considered from the base of peduncle to the distal top of the rachis. The rachis was divided into three zones of roughly similar length, namely the basal, medial and distal zones, following the methodology used in previous studies of octocorals (see for example Orejas et al. 2002, Soong 2005, Baillon et al. 2014b). The two quantitative morphological traits analysed were sclerite and polyp lengths. Twenty sclerites from the upper and lower area of the rachis-peduncle boundary (10 for each zone) were extracted, measured and compared. From each rachis zone, five autozooids were randomly selected, avoiding those located at the base of the leaf. Comparisons between polyps were carried out within each zone and between zones of different colonies. To estimate the measurement error (ME hereafter), three independent observations were made on different days. All measurements were obtained using the ImageJ 1.38x program (Wayne Rasband, National Institutes of Health, USA).

The statistical analyses were performed using the R computing environment v3.5.0. (R Core Team 2018). The libraries and functions used in each case are specified below. Subindividual variability in the two traits examined was estimated using two methods commonly used in plant studies: the intracolony coefficient of variation (CV hereafter) and variance partitioning (for more details see Herrera 2009). In the first method, the CV (calculated as the intracolony standard deviation/individual mean) is a relative dispersion measure that evaluates the proportion of standard deviation (sd) with respect to the mean (X) and enables a comparison of dispersion between different groups or variables (which may have different measurement units) (Herrera 2009, 2017). In order to test differences between the species and between the interaction of “species per zone of the colony”, an ANOVA test was applied to the model. The emmeans 1.4 package was used to calculate the marginal averages (Lenth 2019). The Levene test (Levene 1960) was used to compare the levels of variability between species, the type selected being the Levene median-log test (Schultz 1985). The function used was leveneTest included in the car package (Fox and Weisberg 2019), with the option centre = “median”. Variabilities between species were compared using “marginal averages”, i.e. the averages of the dependent variable (polyp size) for the different levels of one or more categorical predictor variables (see Herrera 2009).

The second method for assessing intracolony variation in quantitative traits consisted of partitioning the total variance of each trait (VAR_total) into its additive between- (VAR_between) and intracolony (VAR_within) components (Herrera 2009). This was carried out using a full random hierarchical mixed model (see below). One drawback of variance partitioning is that the within-individual component of variance may be inflated by ME unless individual reiterated structures are measured repeatedly, thus allowing proper estimation (Herrera 2009). In this approach a precise estimation of ME was performed using repeated measurements from each specimen in order to dissect true within-individual variance and measurement error. Thanks to the analysis of repeatability of measurements through the fully nested random design, the ME component was revealed to be different for both traits, being significantly low in Ptilella-Pennatula analyses (0.071%). This means that most “residual” variance in models is actually attributable to differences between polyp lengths, i.e. the within-individual component (see Discussion). As both methods (CV and VAR_within) focus on different aspects of intracolony variability which may or may not be related (Pearson 1901), they should ideally be used in combination (Herrera 2009), and this is the approach followed in the analyses carried out in this study.

The preliminary analysis performed in Pt. grayi to obtain initial information on the existence and degree of intracolony variability was based on a full random hierarchical mixed model. The required packages were readxl (Wickham and Bryan 2019), dplyr (Wickham et al. 2019), nlme (Pinheiro et al. 2018) and ape 5.0 (Paradis and Schliep 2018). The following more complete analysis using colonies of Ptilella and Pennatula was based on linear and full random hierarchical mixed models (Herrera 2009). The required packages were readxl 1.3.1 (Wickham and Bryan 2019), dplyr 0.8.3 (Wickham et al. 2019), nlme 3.1-137 (Pinheiro et al. 2018), ape 5.0 (Paradis and Schliep 2018), ggplot2 (Wickham 2016), car (Fox and Weisberg 2019) and emmeans 1.4 (Lenth 2019).

The distribution of mean values revealed that the upper zone of the rachis-peduncle boundary contained larger sclerites (200-300 µm) than the lower zone (100-150 µm). In most colonies the mean polyp length increased from the basal (4-6 mm) to the medial-distal zones of the rachis (approximately 6-8 mm). For both traits variances differed within and between colonies, as shown by differences in interquartile ranges (see Fig. S1). The CV of sclerite size ranged between 0 and 0.02, while the CV for polyp size ranged between 0 and 0.25 (Fig. 1). This indicates that the range of variability of polyp size was higher than that for sclerite size. The CV behaviour of polyp size by zone in the different colonies showed practically antagonistic patterns (e.g. between Pt.gy-1 or Pt.gy-2 and Pt.gy-3, Pt.gy-7 or Pt.gy-8). The colonies Pt.gy-4 and Pt.gy-6 were similar but differed from the rest. Regardless of the numerical results, which can be very marked by Pt.gy-1 and 2 with very high mid-zone CV, in general no clear patterns were observed using CV. However, the model confirmed significant differences for both traits within colonies (p-value <0.001; ANOVA test), and significant differences between colonies only in the case of polyp size (p-value =0.184 for sclerites, p-value <0.001 for polyps; ANOVA test) (see Table S1). Mean, standard deviation and CV values from each colony are summarized in Table S2.

The other approach to quantitatively assessing intracolony variation consisted of partitioning the total population-level variance of the trait (VAR_total) into its additive between-colony (VAR_between) and intracolony (VAR_within) components. For the two traits considered, the variance between colonies was negligible (VAR_between = 2.56e-04 for sclerite size; 6.20e-08 for polyp size). In contrast, variation of residual values was a substantial source of variation (VAR_within = 1.16e03 for sclerites; VAR_within = 0.92 for polyps; see Table S3). Expressing variances as percentages of the total (%VAR_total) confirmed the negligible variation between colonies for both traits (VAR_between <1%) and extensive intracolony variance (VAR_within 28% and 49% for sclerites and polyps, respectively) (Table S3).

The functional nature of sclerites and polyps is directly related to the structure, feeding and fitness of the colony (Hickson 1916, Williams et al. 2012). However, as shown here, spatial distribution of variability in some of their features (e.g. length) was not homogeneous along the vertical axis of colony, suggesting a certain “sectoriality”, as has also been reported in plants (Herrera 2009). For example, sclerites of Pt. grayi exhibited a higher variation above the rachis-peduncle boundary (mud-water gradient) in the basal portion of the rachis. Menezes et al. (2013) suggested that the rachis is a higher environment pressure zone than the peduncle (buried) due to competitive contact with surrounding organisms. In a similar way, Herrera (2009) observed differences in variability distribution of certain traits in aquatic plants related to the water-air gradient. In most colonies analysed here, the polyp length showed higher variation levels at the medial-distal than the basal zone of the rachis. Among the factors that may modulate intracolonial variation, some authors have suggested the feeding source, fecundity investment, continuous burial by sediment or competitive contact with surrounding organisms (Soong and Lang 1992, Goffredo et al. 2011, Menezes et al. 2013). This has been suggested for Anthoptilum grandiflorum, Balticina finmarchica and Pennatula aculeata (Baillon et al. 2016). Thus, polyps located in the medial-distal zone of the rachis, far from the substrate, might be subject to higher environmental stress conditions than polyps located in the basal zone.

The study of the continuous variation in quantitative features of homologous organs within the same plant remains relatively unexplored from an ecological perspective despite being a quintessential plant feature. However, the few studies that have been carried out have provided an enormous amount of information and suggested manifold biological and ecological implications (Herrera 2009, 2017). The study of marine modular organisms, including sea pens (such as number of features, individuals and species analysed), is even more limited, but some noteworthy parallelisms can be inferred when the distribution of variability is compared.

When the levels of subindividual variability of reiterated structures obtained here for Ptilella and Pennatula colonies are compared with the subindividual variability reported for certain plant organs (Figure 4), we can see that polyps (more intracolonially variable than sclerites) show lower levels of variability than leaves (fruits or seeds) but comparable levels to flowers (Fig. 4). Like polyps, flowers are the organ responsible for reproduction (perhaps the burden of reproduction may be exerting more influence on the distribution of variability levels of polyps than the capture of resources). In any case, in both reiterated structures (flowers and polyps), there is an important source of variability that must be considered. What might this parallelism involve for pennatulaceans from an ecological point of view? The similar organization of phenotypic variation in pennatulaceans to that observed in plants supports the hypothesis that subindividual variability is an emergent property of individual organisms caused by their modular construction and the reiteration of homologous structures with the same function (leaves, fruits, flowers or seeds in plants, sclerites or polyps in sea pens), regardless of their evolutionary origin. Parallelisms like this between plants and animals allow us to determine the consequences of certain natural processes that are known to occur in one of the two kingdoms and infer them in the other (Herrera et al. 2015, Alonso et al. 2018). Following this idea, the influence exerted by the associated fauna (e.g. animal consumers) in the organization of within-individual variance in plants (Herrera 2009), could have similar effects between predators and sea pens (Clippele et al. 2015), opening a new field of study. For example, it has been suggested that subindividual variability in some functional traits of leaves may be advantageous to individuals by enhancing whole-plant photosynthetic performance and optimizing the exploitation of environmental variation (Osada et al. 2014, Herrera 2017). It is known that some animal consumers have the opportunity to discriminate between not only individual plants but also the multiplicity of non-identical organs borne by each of them, affected by the distribution of variability of reward offered by plants (Herrera 2009, 2017). In the same way, subindividual variability of polyp traits could be related to the optimization of resource exploitation via “division of labour” within the individual (Sides et al. 2014, Siefert et al. 2015). In this hypothetical case, sea pen predators (such as Tritonia sp. on Pteroeides sp. in Duncan 1998, García-Matucheski and Munian 2011; or Armina sp. on Veretillum cynomorium or Ptilosarcus guerneyi in Birkeland 1974, Jones et al. 2000, Buhl-Mortensen et al. 2010, among others) would exert a certain discrimination between colonies and could likewise be influenced by levels of subindividual variation, an interesting hypothesis that could be tested by future studies on pennatulaceans. Perhaps in this hypothetical scenario, following the parallelism with plants, the term heterozooidy (analogous to the term heterophylly, see Herrera 2009) may be considered appropriate to refer to this subindividual variation within the same zooid type.