The same RNA samples used in previous studies including different developmental stages of A. digitifera (eggs, blastulae, gastrulae, planula larvae, early polyps, and adult branches) were used for cDNA synthesis (Shinzato et al., 2011; Shinzato et al., 2021a; Yoshioka et al., 2021). Samples were homogenized in TRIzol reagent (Invitrogen) on ice, and total RNA was extracted following the manufacturer's instruction and Chiu et al. (2020). Extracted total RNA was treated with recombinant DNase I (Roche) to remove contaminating genomic DNA. First-strand cDNA was synthesized from 5 μg of DNase-treated RNA using SuperScriptTM III reserve transcriptase (Thermo Fisher Scientific). Possible DL-L genes of A. digitifera were derived from predicted genes of A. digitifera (Shinzato et al., 2021a), and primers amplifying full length open reading frames (ORFs) were designed for each gene ( Supplementary Table 1 ). PCR products were cloned into the pGEM-T Easy cloning vector (Promega), transformed into competent cells (Escherichia coli, DH5α strain), and plasmids were sequenced at a commercial laboratory (Eurofin Genomics, Tokyo, Japan).

To investigate gene expression patterns of 18 DL-L genes predicted from the A. digitifera genome assembly (Shinzato et al., 2021a), we used previously reported RNA-Seq transcriptome data of A. digitifera from different developmental stages (egg, blastula, gastrula, planula, polyp, Shinzato et al., 2021a), adult stages (Yoshioka et al., 2021), and increased seawater temperature treatments (25°C to 30°C, Shinzato et al., 2021b). Briefly, Illumina adaptor sequences and low-quality reads (Quality score < 20, length < 25 bp) in the RNA-Seq data were trimmed with CUTADAPT v1.16 (Martin, 2011), and cleaned reads were mapped to A. digitifera gene models (Shinzato et al., 2021a) using KALLISTO v0.44.0 (Bray et al., 2016) or MiniMap v2.9 (Li, 2018) with default settings. Mapping counts were normalized using the trimmed mean of M values (TMM) method, and then converted to counts per million (CPM) using EdgeR v3.28.1 (Robinson et al., 2010; McCarthy et al., 2012) in R v3.6.3 (team 2015).

Amino acid sequences were deduced from cloned cDNA sequences using the ExPASy translate tool (https://web.expasy.org/translate/; Gasteiger et al., 2003). To search for genes similar to DMSP lyase gene from transcriptomic databases, the deduced amino acid sequences of cloned A. digitifera DL-L genes, Emiliania huxleyi Alma1 (KR703620.1), and Symbiodinium A1 DMSP lyase (P0DN22) were used as query sequences in homology searches against the Transcriptome Shotgun Assembly (TSA) database at the National Center for Biotechnology Information (NCBI) using tBlastn (A. digitifera DL-L genes: e-value cutoff 1e^-50, query coverage >80%, Emiliania huxleyi Alma1 and Symbiodinium A1 DMSP lyase: e-value cutoff 1e^-5, query coverage >70%). Sequences satisfying these criteria were retrieved, and amino acid sequences were prepared using the ExPASy translate tool.

Sequences downloaded from NCBI TSA databases, including cnidarians, animals, phytoplankton and macroalgae were used for phylogenetic analyses, with bacterial DMSP lyase genes (WP028324949.1, WP012448288.1, WP034745099.1, WP020674609.1, and WP006965703.1) as an outgroup. First, all amino acid sequences were aligned using MAFFT v7.310 (Katoh et al., 2002; Katoh and Standley, 2013), and poorly aligned sequences were removed manually. Gaps in alignments were removed using trimAl v1.2 (Capella-Gutiérrez et al., 2009) with the “gappyout” option. After removing gaps, maximum likelihood analyses were performed using RAxML v8.2.10 (Stamatakis, 2014) with the “bootstrap 100” and “protgammaauto” options.

Among 18 DL-L genes predicted from the A. digitifera genome (Shinzato at el., 2021a), we successfully obtained 13 DL-L gene sequences with full-length mRNA coding regions from cDNA ( Figure 2 ; Supplementary Figure 1 ), demonstrating that multiple DL-L genes exist in the A. digitifera genome and that these are actually expressed as possible functional mRNAs without in-frame stop codons. The resultant 13 amino acid sequences had the same Asp/Glu/hydantoin racemase superfamily conserved domain as in Alma1 of E. huxleyi ( Figure 2 ; Supplementary Figure 1 ) (Alcolombri et al., 2015). These 18 genes showed different levels and patterns of gene expression in different developmental stages ( Supplementary Table 2 ) as shown by previous studies (Shinzato et al., 2021a; Shinzato et al., 2021b; Yoshioka et al., 2021). Gene expression of five genes was significantly changed (one upregulated and four downregulated) by increased sea water temperature (25°C to 30°C, q < 0.05, Supplementary Table 2 ).

Among the 6,258 eukaryotic TSA databases deposited in NCBI, we identified 123 gene sequences from 72 species that showed similarity to A. digitifera DL-L genes and satisfied the above criteria (tBlastn e-value cutoff 1e-⁵⁰, query coverage >80%). Most sequences were identified from the Cnidaria (stony corals, sea anemones, jellyfish and hydra), phytoplankton (Phyla Haptophyta, Myzozoa, Ochrophyta) and macroalgae (Phyla Chlorophyta and Rhodophyta) ( Table 1 ). Unexpectedly, four sequences were also identified from other animals, including a ctenophore (Beroe forskalii), a mollusk (Limacina retroversa) and arthropods (Dendroctonus frontalis and Pleuromamma xiphias) ( Table 1 ). We further investigated the amino acid sequence similarity of these sequences with the only known DMSP lyases from eukaryotes to date (Alcolombri et al., 2015), Emiliania huxleyi Alma1 and Symbiodinium A1 DMSP lyase ( Supplementary Table 3 ). Most eukaryote sequences similar to A. digitifera DL-L genes ( Table 1 ) also showed sequence similarities to both DMSP lyases, but percent identities were lower than those to A. digitifera DL-L genes ( Supplementary Table 3 ).

In order to infer evolutionary relationships of the identified sequences from eukaryotes and the four sequences from non-cnidarian animals, we performed molecular phylogenetic analyses with amino acid sequences of the 13 A. digitifera DL-L genes as representatives of Cnidaria and bacterial DMSP lyase genes from prokaryotes as an outgroup. Sequences from prokaryotes and eukaryotes separated into two distinct clusters, indicating that all of the sequences obtained from eukaryote TSA databases originated from eukaryotes, and were not contaminants from symbiotic bacteria ( Figure 3 ; Supplementary Figure 2 ). All sequences from A. digitifera clustered, demonstrating that they originated from the host coral ( Figure 3 ; Supplementary Figure 2 ). None of the sequences from non-cnidarian animals clustered with those from A. digitifera. Clustering of the D. frontalis sequence with dinoflagellate sequences displayed 100% bootstrap support, indicating that this sequence was derived from symbiotic or adhesive algae, not from the animal itself. However, the other three sequences from B. forskalii, P. xiphias and L. retroversa did not cluster with those of A. digitifera and other eukaryotes ( Figure 3 ; Supplementary Figure 2 ); thus, the evolutionary origin (animals or other eukaryotes) are unclear.

Among 137 Cnidarian TSA databases in NCBI, we found sequences similar to DL-L genes from 35 species ( Table 2 ). Most of these sequences were from Acropora millepora and A. tenuis and had been previously identified (Alcolombri et al., 2015; Shinzato et al., 2021a). Interestingly, possible DL-L genes were found not only in the subclass Hexacorallia, including stony corals, sea anemones, and zoanthids, but also in the Octocorallia, including soft corals and blue corals, and even from Hydroidolina, including jellyfishes and fire corals ( Table 2 ). We also investigated amino acid sequence similarity of these cnidarian sequences with known eukaryote DMSP lyases from E. huxleyi and Symbiodinium A1 (see above), and all sequences except one from Seriatopora caliendrum had significant similarities with these eukaryote DMSP lyases ( Supplementary Table 4 ).

Coral and/or cnidarian transcriptomic databases often contain sequences originating with symbiotic algae, as well as associated microorganisms. In order to identify which sequences were from animal hosts, we further performed molecular phylogenetic analyses of sequences from cnidarian TSA databases together with DL-L genes from the Symbiodiniaceae, coccolithophores and prokaryotes ( Figure 4A ). Six distinct clades were supported by strong bootstrap probabilities (93 ~ 100%). Most of the sequences belonging to the Hexacorallia clustered together in Clade 1 with A. digitifera DL-L genes, indicating that these originated from hexacorallian hosts. Interestingly, Clavularia sp. of the Octocorallia and Velella velella of the Hydroidolina were also relegated to Clade 1, indicating that these share a common ancestral gene with hexacorallian species. The family Symbiodiniaceae is genetically diverse, and includes multiple genera (LaJeunesse et al., 2018). As expected, some Cnidarian sequences clustered with sequences of symbiotic algae (Symbiodiniaceae) and were assigned to Clades 3-5 ( Figure 4A ; Supplementary Figure 3 ), indicating that these were contaminant from symbiotic algae. Among them, a sequence from Pocillopora acuta was clustered into Clade 3, sequences from Anthopleura elegantissima, Millepora squarrosa, Agaricia lamarcki, Heliopora coerulea, Zoanthus sp., Palythoa sp., M. alcicornis, Millepora sp. were clustered in Clade 4, and sequences from Fimbriaphyllia ancora, M. complanata, and Xenia sp. were clustered in Clade 5 ( Figure 4A ; Supplementary Figure 3 ; Table 2 ). Different cnidarian species harbor different types of symbiotic algae; thus, Clades 3-5, containing Symbiodiniaceae could reflect different types of symbiotic algae. All DL-L genes of the coccolithophore E. huxleyi within the Haptophyta clustered in Clade 2, and bacterial DMSP lyase genes clustered in Clade 6. Taken together, all DL-L genes of Cnidaria, including scleractinian corals, soft corals, and jellyfishes, evolved from an ancestral gene that already existed in the last common ancestor of Anthozoa and Hydrozoa, and have completely different evolutionary backgrounds from those of coccolithophores, the Family Symbiodiniaceae, or prokaryotes.

It has been reported that coral reefs are hotspots for DMSP, which have been attributed to symbiotic algae of corals (Broadbent et al., 2002; Broadbent and Jones, 2004; Broadbent and Jones, 2006; Swan et al., 2012). Recent studies have shown that adult and juvenile A. tenuis and A. millepora without symbiotic algae both produce DMSP (Raina et al., 2013). DL-L genes have also been discovered, not only in Acropora genomes (Alcolombri et al., 2015; Shinzato et al., 2021a), but also in genomes of Montipora, Astreopora, Goniastrea, and two corallimorpharians, Amplexidiscus and Discosoma (Shinzato et al., 2021a).

In this study, we successfully identified full-length ORFs of 13 DL-L genes, which are expressed in A. digitifera, and we identified DL-L genes from a variety of scleractinian lineages ( Figures 1 , 4A ; Supplementary Figure 3 ). The earliest coral fossil record of a Scleractinian dates to the middle Triassic (240 Ma) (Stanley, 2003). Most extant scleractinians are classified into two major clades, known as the Complexa (complex corals) and Robusta (robust corals), and are assumed to have diverged in the Late Carboniferous (300 Ma) (Romano and Palumbi, 1996; Romano and Cairns, 2000). In genomes of robust corals, a DL-L gene has been detected only in Goniastrea aspera to date (Shinzato et al., 2021a). However, in this study, we identified animal-type DL-L genes from four robust corals, including Cyphastrea serailia, Favites colemani, Favia lizardensis, and Seriatopora caliendrum ( Figure 4B ; Table 2 ), indicating that DL-L genes have been preserved in a variety of complex and robust corals, but were lost from some lineages.

Although a large proportion of DL-L genes were identified in corals and cnidarians, we also identified similar sequences from non-cnidarian animals, B. forskalii, P. xiphias and L. retroversa. However, none of these share a common ancestry with coral DL-L genes ( Figure 3 ). Although the origins of these sequences are not clear, we suggest that they may have come from symbiotic/adhesive algae, as in the case of a sequence from D. frontalis ( Figure 3 ). Based on our phylogenetic analyses using currently available transcriptomic databases, we conclude that, to date, only cnidarians possess DL-L genes that are clearly of animal origin. Further addition of more phytoplankton DL-L genes may reveal the origins of these from non-cnidarian species.

The Phylum Cnidaria contains three clades (subphyla): Anthozoa (comprising the Octocorallia, Hexacorallia, and Ceriantharia), Endocnidozoa (a clade of parasites) and Medusozoa (consisting of the Cubozoa, Hydrozoa, Scyphozoa, and Staurozoa) (Collins, 2009; Kayal et al., 2018). In this study, we identified eukaryotic DL-L genes, not only in the Scleractinia and Corallimorpharia but also in the Zoantharia (Hexacorallia), Alcyonacea (Octocorallia), and Hydrozoa (Anthoathecata, Hydroidolina) ( Figure 4B ; Table 2 ). This indicates that the last common ancestor of Anthozoa and Hydrozoa possessed a DL-L gene and that this gene has been inherited in a wide range of cnidarian species. The earliest cnidarian fossils occur in strata of the Ediacaran (560 Ma, Liu et al., 2014), and together with molecular and paleontological analyses, suggest that the phylum Cnidaria probably originated during the Precambrian Eon from the Cryogenian to the Ediacaran (700-595 Ma, Peterson et al., 2004; Erwin et al., 2011). Divergence of 2 major taxa (Anthozoaria and Medusozoa) may have occurred before the Cambrian (543 Ma) (Cartwright et al., 2007; Park et al., 2012). Consequently, we hypothesize that the DL-L gene may be have originated from an ancient gene that already existed in the last common ancestor of Cnidaria in the pre-Cambrian ( Figure 4B ).

How did cnidarians, including corals, obtain DL-L genes? Shinzato et al. (2021a) hypothesized that DL-L genes may have been acquired by the common ancestor of scleractinians and corallimorpharians, both of which are hexacorallians (Anthozoa), via horizontal gene transfer from symbiotic Symbiodiniaceae or Emiliania. However, as mentioned above, a variety of anthozoans and hydrozoans possess these genes ( Figure 4A ). McFadden et al. (2021) suggested that the Anthozoa probably arose in the Cryogenian to Tonian periods (648-894 Ma) and lacked photosymbionts. In the Devonian (383 Ma) from the Palaeozoic Era, anthozoans of the Scleractinia first formed associations with photosymbionts, followed by alcyonacean octocorals (318 Ma) and corallimorpharians (312 Ma), and photosymbioses have been gained and lost repeatedly in all orders through the Jurassic (199-151 Ma, McFadden et al., 2021). In addition, the fossil record shows that scleractinian corals had photosymbionts in most of the upper Triassic since the Mesozoic Era (Stolarski et al., 2011; Frankowiak et al., 2016). Recent molecular dating estimates suggest that the earliest diversification of the Symbiodiniaceae occurred in the Jurassic (∼160 Ma, LaJeunesse et al., 2018) and the first E. huxleyi appeared 270,000 years ago (Thierstein et al., 1977; Paasche, 2001), which were much later than the origin of the Anthozoa. Our phylogenetic analysis ( Figure 4 ) confirms distinct ancestries of DL-L genes in Symbiodiniaceae and cnidarians, indicating that they have completely different evolutionary backgrounds. Therefore, acquisition of DL-L genes via horizontal gene transfer from symbiotic algae (Shinzato et al., 2021a) may not have occurred. Based on the present findings, we propose an updated hypothesis: DL-L genes in cnidarians are “ancient genes” in the animal kingdom, dating back to the pre-Cambrian.

As of 10^th February 2022, 137 transcriptomic databases from 91 cnidarian species had been registered in NCBI TSA. Interestingly, we realized that all cnidarian species, except V. velella, possessing animal type DL-L genes (15 species) are limited to coral reefs ( Table 1 ). V. velella ( Figure 4 ; Table 2 ), known as “by-the-wind sailor”, is a cosmopolitan, free-floating, colonial hydrozoan that lives mainly at the water/air interface and floats in temperate and tropical seas (Araya and Aliaga, 2018). This implies that DL-L genes may be essential for animals to survive in coral reef or warm, shallow-water environments, although we acknowledge that publicly available transcriptomic databases of cnidarians to date may be limited to coral reef species.

Comparative genomic analysis using reported scleractinian coral genomes showed that DL-L genes are the most diversified gene family among gene families that significantly increased number of genes in the last common ancestor of Acropora (Shinzato et al., 2021a). Why did gene duplication events of DL-L gene specifically occur in the Acropora lineage? The earliest fossil records of Acropora are known from the late Paleocene (65-54 Ma) in Somalia and Austria (Carbone et al., 1993; Baron-Szabo, 2006), indicating that Acropora has existed for more than 50 million years. Acropora was the main reef builder in the Oligocene (28-23 Ma) of Greece and Early Miocene of Egypt (Schuster, 2002), and was widely distributed throughout the Miocene in different regions, including the Indo-Pacific and Caribbean (Budd, 2000). It is suggested that Acropora began to spread throughout the world in the Cenozoic and species diversification occurred in the Eocene and Oligocene (around 25-50 Ma). Molecular dating analysis using whole-genome data indicates that the Acropora ancestor survived warm periods without sea ice from the mid or late Cretaceous to the Early Eocene, when gene expansion of DL-L genes occurred specifically in the Acropora ancestor (Shinzato et al., 2021a). In addition, Acropora species also have high concentrations of DMSP compared to other corals (Broadbent et al., 2002; Guibert et al., 2020), suggesting that Acropora corals actively utilize DMSP and that high concentrations of DMSP may have triggered Acropora-specific duplication of DL-L genes. Several studies have shown that DMSP is involved in a wide range of coral stress responses, including responses to heat, sunlight, air exposure, and hyposalinity (Sunda et al., 2002; Raina et al., 2009; Deschaseaux et al., 2014; Aguilar et al., 2017). Not all A. digitifera DL-L genes (5 out of 18) responded to increased sea water temperature and response patterns of differentially expressed genes vary (1 upregulated and 4 downregulated, Supplementary Table 2 ), suggesting that the functions of the DL-L genes in A. digitifera have also diverged. Taken together, diversified DL-L genes in Acropora may have acquired new or different functions, not only mediating cleavage of DMSP into DMS, but that they may assist Acropora in adapting to environmental changes, for example from intense heat, light, and salinity (Shinzato et al., 2021a). Eventually Acropora may become the dominant coral genus in extant coral reefs.

The functions of DL-L genes in corals are completely undetermined at this stage. Interestingly, although overexpression of the Alma l genes from E. huxleyi and Symbiodinium A1 in E. coli cells had high DMSP lyase activities, an Alma1 homolog from A. millepora showed almost no DMSP lyase activity (Alcolombri et al., 2015), indicating that not all DL-L genes in eukaryotes retain their DMSP lyase function, and that acquisition of other functions or functional differentiation may have occurred during the process of gene duplication. In particular, duplicated DL-L genes in Acropora will need to be investigated to determine which genes have DMSP lyase activity. We identified expressed sequences of DL-L genes from A. digitifera RNA, and these could be used for molecular and functional characterization in DMSP lyase assays. We also found that expression levels of DL-L genes were indeed diverse in A. digitifera. Therefore, even if some DL-L genes in corals do not function as DMSP lyases, they may have other functions enabling adaptation to coral reef or warm, shallow-water environments. Identifying their actual biological functions in Acropora corals will be important to understand not only adaptation mechanisms to shallow and warmer environments, but also the impact of Acropora corals on the sulfur cycle in the oceans.