The materials examined herein were collected during various surveys over different geographical areas and sampling programmes: the northeastern Arctic-Atlantic (BIOICE programme), the northeastern Atlantic (Scotia cruises, INDEMARES Chica), the Mediterranean Sea (INDEMARES Alborán, INDEMARES Cap de Creus), the southeastern Atlantic (BENGUELA VIII) and the Antarctic Peninsula and the eastern Weddell Sea (Polarstern cruises ANT XVII/3, ANT XIX/5, and ANT XXIII/8) (see Table 1).

During the different expeditions, the specimens were sorted and labelled on board. The colonies (or a tissue sample from each one) were directly fixed in 100% ethanol for further molecular studies. The remaining part of the colonies was fixed in hexamethylenetetramine-buffered 4% formalin-seawater or 70% ethanol. After the fixation period, all colonies were preserved in 70% ethanol. The sequenced voucher specimens are deposited in the Museu de Zoologia de Barcelona (MZB), in the Muséum National d’Histoire Naturelle (MNHN) in Paris and in the collection of the Biodiversidad y Ecología Acuática research group of the University of Seville (BECA).

The analysis of time divergence was undertaken within a Bayesian framework in BEAST 2.5.0 (Bouckaert et al. 2014), using the 45 taxa from which at least two of the three loci were available. BEAST allows topologies to be considered “fixed” or estimated to accommodate for phylogenetic uncertainty (Bouckaert el al. 2014, Drummond and Rambaut 2007). Here, we chose the second approach because the posterior node probabilities of Bayesian phylogeny were relatively low in few cases. Several works have emphasized the importance of the rigorous selection of an appropriate clock model (e.g. Duchêne et al. 2014), and to this end four models were compared: i) relaxed with an exponential local type distribution; ii) relaxed with a log-normal type distribution; iii) random local; and iv) strict. For each model, we estimated their marginal likelihoods using the nested sampling approach implemented in the NS package of BEAST 2 (Maturana et al. 2018) and used the Bayes factor to obtain the best-supported model. In each of these models, the three partitions were treated as linked, and we used a pure-birth (Yule) tree prior (Bouckaert el al. 2014, Drummond and Rambaut 2007, Drummond and Bouckaert 2015). The Yule model is a simple model of speciation that is generally more appropriate for considering sequences from different species (Drummond and Bouckaert 2015). The absolute estimates of divergence times were calculated after one-fossil calibrations (see Fossil calibration section below). For the best-supported model of molecular clock, two independent runs of 100 million generations, sampling every 10000 generations, were performed. Runs were considered complete with effective sample sizes greater than 200 for all parameters (Rambaut et al. 2018). LogCombiner (Rambaut and Drummond 2010) was used to combine the log files from the independent BEAST runs. TreeAnnotator (Rambaut and Drummond 2010) was used to summarize resulting tree samples into a single consensus tree using the maximum clade credibility, mean height options and discarding 10% of trees as burn-in. The maximum clade credibility summarized the 95% highest posterior density (HPD) limits of the node age.

Fossil calibration. The oldest undisputed pennatulacean fossil comes from several Late Cretaceous beds in Europe and North America (Reich and Kutscher 2011). These forms were assigned to the genera ‘Graphularia’ (4 spp.) and Glyptosceptron (1 sp.; see details in Reich and Kutscher 2011), which are from the lower to upper Maastrichtian (66-72.1 Ma) of the Netherlands, Germany, Belgium, USA and Ukraine. The other well-recognized fossil record of pennatulaceans is from the Early to Middle Eocene (41.2-47.8 Ma) of Trinidad (Pointe-a-Pierre Formation) based on the recognition of Virgularia presbytes (Bayer, 1955), a genus with extant representatives. The age of the most recent common ancestor of the genus Virgularia (ca. 41.2 Ma) was set as an offset, while the mean and standard deviation of the log-normal distribution was set as M=2 and S=1. We avoid using the age of the oldest undisputed pennatulacean fossil in our time calibration analyses (Reich and Kutscher 2011) because it would certainly restrict the early origin of this poorly fossiliferous clade. In the following paragraphs we use the standard abbreviation ‘Ma’ (= million years ago) for the age of a specific moment in the geological past (Aubry et al. 2009).

The estimated marginal likelihood for the time-tree analysis under different clock models indicated that the relaxed log-normal is the best supported one (Table S2). The result of this model is represented in Figure 2 (and Fig. S1). The divergence dates inferred by the Bayesian relaxed clock analyses indicated Lower Cretaceous age (Berriasian age; inferred mean age of 144.3 Ma; HPD 95%: 65.7–249.1 Ma) for the origin of Pennatulacea. In addition, the four main lineages within Pennatulacea were originated soon after, with Clade I and II diverged during the Hauterivian age (Lower Cretaceous; inferred mean age of 132.4 Ma; 95% HPD: 64.1–229.9 Ma). The groupings obtained in the BEAST analysis were similar to those from MrBayes, even considering the uncertain affinities of the genus Funiculina (Figs 2, S1). The last genus appears to begin its diversification in the middle Miocene (inferred mean age of 15.6 Ma; 95% HPD: 2.2–36.4 Ma). The most recent common ancestor of the genera Halipteris, Kophobelemnon and Gyrophyllum is from the Coniacian age (Upper Cretaceous; inferred mean age of 89.3 Ma; 95% HPD: 35.2–161 Ma). The beginning of the diversification of most genera with more than one sequence, excepting Virgularia, Gyrophyllum and Ptilella, occurred during Oligocene and Miocene times (ranging from 33.9 to 5.3 Ma). Finally, veretillids (including Cavernularia and Veretillum), appear to have diversified only after the late Eocene (inferred mean age of 36.9 Ma; 95% HPD: 10.3-70 Ma).

Hickson (1916:131) had already highlighted the difficulty of carrying out systematic studies on Pennatulacea based on features such as colonial symmetry because of the wide range of variation of those characters overlapping among genera and species. As has happened before with other octocoral groups, there is an evident lack of agreement between the current systematics of pennatulaceans based on morphological characters (see Kükenthal 1915, Hickson 1930, Williams 1997, among others) and phylogenetic hypotheses generated from molecular data (McFadden et al. 2006, Dolan et al. 2013, Kushida and Reimer 2018, among others). The presence of homoplastic characters (e.g. convergence and reversals) is a serious concern in the formulation of a reliable evolutionary hypothesis of pennatulaceans based only on morphology (Landing et al. 2015, Pérez et al. 2016).

However, several synapomorphies observed in sea pens have traditionally been used to differentiate pennatulaceans from other octocorals and support them as a monophyletic group (Hickson 1930, 1937, Williams 1994). In the last decades, despite the inclusion of unsuitable markers or doubtful sequences that questioned this common origin (see results of Berntson et al. 1999, 2001), most of the molecular studies carried out have corroborated the monophyletic origin of pennatulaceans with strong support (McFadden et al. 2006, Everett et al. 2016, Kushida and Reimer 2018).

The use of the two traditional suborders, Sessiliflorae and Subselliflorae (Kükenthal 1915), until the early 1990s (e.g. Hickson 1916, Tixier-Durivault 1965, Williams 1990) has been avoided in recent decades since Williams (1995: 136) questioned the validity of these groupings. In monographic works the different families were simply listed from the most structurally simple colonies (radial forms such as veretillids) to the most structurally complex ones (bilaterals bearing well-developed polyp leaves such as pennatulids), the latter often being considered from a morphological point of view as the most derived forms (Williams 1997, López-González and Williams 2002). In the present study, using a molecular barcode for octocorals and a wider taxa sampling range, this traditional classification into suborders has been clearly shown as a non-monophyletic one, in agreement with previous molecular studies (McFadden et al. 2006, Dolan et al. 2013, Kushida and Reimer 2018). Therefore, although Sessiliflorae and Subselliflorae are still mentioned from a practical point of view (e.g. Yesson et al. 2012, Hogan et al. 2019), it is evident that they are non-monophyletic groupings, the internal morphological evolution of each clade still being difficult to understand.

Morphological traits have been used to categorize the extant pennatulaceans into 14 families, 37 genera and approximately 200 considered valid species (Williams 2011, 2015, García-Cárdenas et al. 2019) according to the most important contributions concerning the order Pennatulacea (e.g. Kükenthal and Broch 1910, Kükenthal 1915, Hickson 1916, among many others). However, recent molecular studies have failed to recover the monophyly of some of these families, such as Umbellulidae, Pennatulidae, Virgulariidae, Protoptilidae, Scleroptilidae, Stachyptilidae and Kophobelemnidae (Dolan et al. 2013, Kushida and Reimer 2018). The present study, including seven genera not previously sequenced and additional mitochondrial sequences and molecular markers, confirmed the non-monophyletic nature of some of these families (such as Protoptilidae, Virgulariidae or Pennatulidae).

The family Veretillidae, which includes most of the structurally simplest and radial colonial forms observed in pennatulaceans, have been considered since the early 19th century as transitional forms between soft corals and pennatulaceans (Koch 1878). The idea that Veretillidae and Echinoptilidae were the earliest-diverging taxa in phylogenies based on morphology spread during the early (Kükenthal and Broch 1910, Niedermeyer 1913, Hickson 1916) and late (Williams 1994, 1997) 20th century and the transition to the 21st century (López-González and Williams 2002), before molecular approaches postulated a different origin and basal relationships of pennatulaceans from those of the rest of the octocoral taxa (McFadden et al. 2006, Dolan et al. 2013). In contrast to that traditional idea, and following later molecular postulates, our phylogeny does not support the hypothesis that Veretillidae and Echinoptilidae are located at a basal position, exemplifying the ancestral sea pen morphology.

Recent phylogenies based on the mitochondrial markers ND2 and mtMutS including sequences of the veretillid genera Veretillum and Cavernulina showed them to be closely related to sequences of the genus Sclerobelemnon (Kophobelemnidae) and placed all of them within Clade I with moderate support (Kushida and Reimer 2018: 9). Our research includes new additional veretillid sequences belonging to Cavernularia (C. pusilla) and Veretillum (including the type species V. cynomorium Pallas, 1766). In our phylogeny these veretillids (Cavernularia and Veretillum) are gathered in a clade with strong support (bst= 100%, PP= 1) within Clade I (Fig. 1). An additional exploratory maximum likelihood analysis based on mtMutS including all available mtMutS sequences from GenBank (Fig. S2) revealed that sequences attributable to the Pacific species Veretillum sp. 1 (MK133435) and Veretillum sp. 2 (MK133539, MK133545, MK133526) from Kushida and Reimer (2018) (although falling into the same Clade I) were not close to our sequence of the type species V. cynomorium, but rather close to Cavernulina and Sclerobelemnon sequences, in a different grouping. Thus, the inclusion of Pacific and Mediterranean (type species) sequences in a single genus would result in a paraphyletic taxon. This fact suggests the need for a detailed morphological study of Pacific colonies, as they could belong to a different genus, despite their possible similar morphological appearance. In short, with current information, the family Veretillidae seems to be a monophyletic grouping within Clade I (Fig. 1). The gathering of further molecular data on the remaining veretillid genera (such as Lituaria or Amphibelemnon) would be desirable in order to better delineate the monophyletic nature of the current list of genera included in the family Veretillidae.

The family Echinoptilidae (including Actinoptilum and Echinoptilum) was traditionally considered the other earliest-diverging taxon in the phylogeny of pennatulaceans, mainly because of its radially symmetrical rachis (or at least in its distal portion) (Kükenthal and Broch 1911, Niedermeyer 1913, Williams 1992). Kushida and Reimer (2018) considered radial symmetry only for those species included in Veretillidae (Veretillum, Cavernulina) within Clade I, the rest of taxa in Clade I and all taxa within the Clades II, III and IV being bilaterals (Kushida and Reimer 2018:5). Nevertheless, both echinoptilid genera exhibit radial symmetry in colonies to different degrees (Williams 1995). The phylogeny proposed by these authors showed the genus Echinoptilum as a derived taxon with strong support within Clade II, with Renilla-Pennatula as the sister group. In our study, the sequence attributable to the other echinoptilid genus, Actinoptilum molle, was also located as a derived taxon within Clade II as the sister group of a [(Renilla-Acanthoptilum)-(Ptilosarcus-Pennatula)] clade.

The family Protoptilidae, including the genera Protoptilum and Distichoptilum, is here recognized as a non-monophyletic taxon, in agreement with Kushida and Reimer (2018). The family Renillidae, constituted by the genus Renilla, was not close to Veretillidae (as was suggested from morphology, see Williams 1997, Pérez and Ocampo 2001), supporting the results of Dolan et al. (2013) and Kushida and Reimer (2018). In our phylogeny, Renilla is the sister group of Acanthoptilum (Virgulariidae), and both are a sister group to Ptilosarcus-Pennatula (Pennatulidae). The family Anthoptilidae, represented here by Anthoptilum sp. 1 MK919656 and A. grandiflorum MK919655, was placed within Clade I as the sister group of Umbellula, and not occupying an ancestral position as was suggested by Dolan et al. (2013, Fig. 1). As commented above, the two best-represented families of pennatulaceans Pennatulidae and Virgulariidae (Williams 1995) are here recognized as non-monophyletic groupings, in agreement with Dolan et al. (2013) and Kushida and Reimer (2018).

The family Umbellulidae is shown here as a non-monophyletic grouping, as observed in previous works (Dolan et al. 2013, Kushida and Reimer 2018) (Figs 1, 2, S2). As the type species of Umbellula, U. encrinus (Linnaeus, 1758), falls into Clade I, that set of sequences has to be considered as the genus Umbellula. Although the genus Umbellula was involved in old hypotheses postulating the deep-sea origin of pennatulaceans (Kölliker 1880), other authors (Marshall 1883) recognized Umbellula as a highly modified form. In the 20th century, the Umbellulidae was considered as a highly derived family among the authors who supported the shallow-water origin of pennatulaceans (Kükenthal and Broch 1911, Williams 1997, among others). Our research supports a later divergence of Umbellula within Clade I (Fig. 2).

The inclusion in our study of the nuclear marker 28S was useful in order to reinforce the support of Clades I and II. However, it also highlighted the unstable location of the genera Funiculina and Halipteris (see Figs 1, S2-S5, and Table S1) depending on the phylogenetic approach used. Clade III and Clade IV identified by Kushida and Reimer (2018) were recovered in our study (with high or relatively low support, bst 96% and 60%, respectively) only when the maximum likelihood method was used, while Bayesian inference indicated a different phylogenetic hypothesis. Therefore, the inclusion of the 28S locus allows us, on the one hand, to reinforce the support of Clades I and II, but also provides evidence of the instability in clades proposed by previous studies, probably requiring greater taxonomic and molecular sampling.

In general, our results, based on a combination of mitochondrial and nuclear segments (mtMutS+Cox1+28S), indicate the monophyletic origin of most of sea pen genera with strong support (1 PP). This is one of the few common points in which morphology and molecular studies agree.

Time divergence estimation in pennatulaceans

Very little is known about the time frame of pennatulacean evolution, and our study is the first attempt to build a time-calibrated phylogeny based on Bayesian relaxed molecular clock analysis. The Lower Cretaceous (Berriasian, ~144 Ma) estimation for the divergence of the Pennatulacea from their sister group Ellisellidae is slightly older than those recently inferred (estimated age of its most recent common ancestor as 82-125 Ma; see Bilewitch 2014). However, it is highly consistent with previous studies regarding estimations of the divergence time of Calcaxonia (with which pennatulaceans have been related; 120-300 Ma; Park et al. 2012), as well as with the age of the oldest undisputed pennatulacean fossils so far recovered (Reich and Kutscher 2011).

The Cretaceous is well-known as a greenhouse period caused largely by increased CO₂ from elevated global volcanic activity (e.g. Takashima et al. 2006, Prokoph et al. 2008), but it also is the last stage in the Gondwana break-up, and displays high rates of seafloor spreading, high sea levels (as much as 260 m above the present) and high ocean temperatures (surface waters >35°C and deep-ocean water >20°C), as well as evidence for changes in global ocean circulation (Haq et al. 1987, Pearson et al. 2001, Friedrich et al. 2012). The global climatic conditions during the Berriasian, when the origin of pennatulaceans likely occurs, were generally arid and the sea level was rather low (Föllmi 2012). The Berriasian also coincides with an extinction phase (Tithonian/Berriasian), which, however, apparently affected continental life more strongly than marine life (Föllmi 2012). In addition, during the Lower Cretaceous several oceanic anoxic events occurred, representing time intervals of usually relatively short duration (<1 Ma) in which intermediate and bottom-water masses became depleted in oxygen (Föllmi 2012).

Our results therefore reinforce the notion that the Cretaceous was a pivotal time in octocoral evolution, because several clades appear to begin their diversification in this epoch (e.g. primnoids; Taylor et al. 2013). The main clades recovered here appear to have diverged quickly (Hauterivian-Berriasian) after their initial pennatulacean origins (Fig. 2), and this finding indicates that the early phase of their diversification took place in a greenhouse period and especially in periods of dramatic environmental changes (Prokoph et al. 2008, Föllmi 2012). Consequently, changes in sea temperature and oceanic anoxic events would have played a role in forcing pennatulacean diversification, and probably also vertical displacement (submersion and emersion of lineages) throughout the complete colonizable bathymetry, as it has been demonstrated that sea pens are one of the zoological groups present over the widest bathymetric range (0-6100 m depth) (Williams 2011). Our results also provide further support for a late origination of veretillids (late Eocene), contrasting with the early notion based on morphological features previously discussed.

Finally, it is evident that the single available fossil calibrations used on the present work have an impact on the quality of the resulting chronogram. We thus stress that this chronogram must be regarded as the first attempt to produce a time-calibrated sea pen tree. Only after achieving a wider sampling of molecular information will we be able to refine or even test the interpretations offered in the present work.