Bioluminescence is a widespread phenomenon that has evolved multiple times across the tree of life, converging among diverse fauna and habitat types. The ubiquity of bioluminescence, particularly in marine environments where it is commonly used for communication and defense, highlights the adaptive value of this trait, though the evolutionary origins and timing of emergence remain elusive for a majority of luminous organisms. Anthozoan cnidarians are a diverse group of animals with numerous bioluminescent species found throughout the world’s oceans, from shallow waters to the light-limited deep sea where bioluminescence is particularly prominent. This study documents the presence of bioluminescent Anthozoa across depth and explores the diversity and evolutionary origins of bioluminescence among the Octocorallia – a major anthozoan group of marine luminous organisms. Using a phylogenomic approach and ancestral state reconstruction, we provide evidence for a single origin of bioluminescence in Octocorallia and infer the age of occurrence to around the Cambrian era, approximately 540 million years ago- setting a new record for the earliest timing of emergence of bioluminescence in the marine environment. Our results further suggest this trait was largely maintained in descendants of a deep-water ancestor and bioluminescent capabilities may have facilitated anthozoan diversification in the deep sea.

The recently published fossil-calibrated phylogeny (Quattrini et al. 2020) (1,729 loci) was used to map bioluminescent traits and depth categories across Anthozoa. 

For the phylogenomic analyses, new fully-resolved phylogeny for Octocorallia reconstructed using target-capture sequence data (from 739 loci from 185 octocoral taxa) (McFadden et al. 2022), was used to investigate the evolution of bioluminescence in octocorals. The tree was constructed using maximum likelihood in IQTree v2.1 using the best model (GTR+F+R10). Divergence-dating of this recently revised octocoral tree was conducted using BEAST v 2.6 as outlined in Quattrini et al. (2020). First, a time-calibrated starting tree of the phylogeny in McFadden et al. 2022 was generated using penalized likelihood method (chronopl, based on 26) in the R package ape. The calibration points included the minimum ages of two fossils and a minimum age at the crown Octocorallia and at the root of the tree, obtained from Quattrini et al. (2020). We then selected a 25-locus alignment to input into the BEAST v.2.6 analysis. In the BEAST analysis, loci were partitioned so that a GTRGAMMA model (initial 1.0, 0 to infinity bounds) was applied to each locus. Exponential priors were used for fossil calibration points and a normal prior was applied to the crown Octocorallia and root of the tree, based on Quattrini et al. (2020). Two separate runs of 200M generations were conducted. Log and tree files from each run were combined in LogCombiner, with a 10% burnin. The combined log file was assessed for convergence of parameter values and age estimates by inspecting traces and effective sample sizes in Tracer v.1.7. Trees were combined from both runs, and resampled by selecting one out of every 10K trees, resulting in 75K trees. TreeAnnotator was then used to produce a maximum clade credibility tree based on mean ages. 

We expanded the highly-supported tree generated using ultraconserved elements (UCE tree) by McFadden et al. 2022, by combining mitochondrial data, mutS-like DNA repair gene (mtMutS), for an additional 107 octocoral taxa (McFadden et al. 2022) to produce a hybrid (UCE/MutS) phylogeny. Phylogenomic analyses were conducted using maximum likelihood and a partitioned analysis for multi-gene alignments with IQTree v2.1 using both the concatenated UCE alignment (739 loci) and the alignment for mtMutS (McFadden et al. 2022). ModelFinder was used to find the best substitution model for each partition (UCE/mtMutS). The partition model was given discrete substitution models for each gene/character set and each partition was allowed to have its own evolution rate. Ultrafast bootstrapping (-bb 2000) and the Sh-like approximate likelihood ratio test (-alrt 1000) were conducted. The reconstructed tree was pruned using the phytools (v1.2-4) package in R (v3.5.0) to remove four aberrantly placed taxa with poorly supported nodes, based on the recently published phylogeny (McFadden et al. 2022). The final hybrid (UCE/MutS) octocoral tree consisted of 270 operational taxonomic units (or OTUs). The Scleralcyonacea clade was rooted to the Malacalcyonacea clade for downstream analysis uin R based on the findings of McFadden et al. (2022).

For Ancestral State Reconstruction (ASR), trait tables were generated for octocoral taxa represented in both the time-calibrated UCE tree (185 OTUs) and the UCE/MutS hybrid tree (270 OTUs) based on a scale of 0-1 (0: Non-bioluminescent or unlikely; 1: Likely bioluminescent or confirmed; NA or 0.5: unknown). Ancestral states of bioluminescence were calculated using stochastic character mapping, sampling ancestral states from posterior probability distributions generated from 100 stochastic character maps for the bioluminescence trait using the make.simmap function (nsim=100) in the R package phytools. The best fit model was determined by fitting and comparing the extended Mk models (fitMk) for discrete character evolution. Akaike information criterion (AIC) values for the fitted models were then compared to determine the model with the lowest AIC score. This was determined to be the ARD model assuming different rates of trait gain/loss, which was utilized for downstream analysis. Phytools was then used to plot one stochastic character map for the bioluminescence trait on both the time-calibrated phylogeny for Octocorallia (UCE tree) and the more inclusive hybrid tree, along with the posterior probabilities (pie charts) at each node.

To estimate the ancestral depth ranges, a Bayesian dispersal-extirpation-cladogenesis (DEC) model implemented in RevBayes was conducted on the dated, maximum clade credibility tree. Ancestral ranges of shallow (< 200 m) and deep-sea (> 200 m) depths were estimated following the guidelines in the online tutorial on simple analysis of historical biogeography (https://revbayes.github.io/tutorials/biogeo/biogeo_simple.html). Current depth ranges were coded as binary values (0,1) whether a species was recorded from deep, shallow, or both depth ranges. Depths were obtained from obis.org and supplemented with recent museum (NMNH) and geome (geome-db.org) records and the literature. For species identified to a rank above genus, the depth of collection was used. In RevBayes, the number of generations (MCMC) was set to 5000. The R package RevGadgets was then used to plot the probability of ancestral states (as pie charts) at the nodes.