Introduction

The biphasic pelagobenthic life cycle is widespread in early-diverging animal lineages (sponges and cnidarians) and in molluscs, echinoderms, annelids and several other bilaterian phyla, suggesting it evolved before metazoan cladogenesis ^1–7. This ancient aquatic life cycle begins with endogenously-regulated embryogenesis that produces a ciliated planktonic larva, and culminates with settlement and metamorphosis into a benthic juvenile. In contrast to embryogenesis, metamorphosis is characterised by three defining features: environmental induction, coordinated transformation of one complex body plan into another, and a dramatic ecological transition ^1,3,6. Despite the scale of morphological and physiological change required, in diverse marine invertebrates this transition occurs strikingly quickly as a means to enhance post-settlement survival ^6,8–10.

Marine larvae are strongly influenced by planetary constants, including daily, tidal, lunar and seasonal rhythms, with gametes or larvae being released at a species-specific time. Most larvae also have developmentally-regulated behaviours that enable settlement to be induced by species-specific biochemical cues that are associated with a favourable settlement site for adult survival and reproduction ^11–16. The capacity to respond to these cues – widely called “larval competence” – develops while larvae swim in the plankton over timescales ranging from minutes to weeks, depending on the species ^{6,11–13,17,18}. Together, these traits indicate that larvae across disparate animal phyla have evolved to integrate constant planetary signals with local ecological cues to optimise dispersal, settlement success and subsequent survival.

Despite the prevalence and ecological importance of this life cycle, it remains unclear how the metazoan genome integrates planetary and ecological signals to regulate the acquisition of larval competence, the decision to settle, and the remarkably rapid remodelling of the body plan at the onset of metamorphosis. Here, we address this problem by exploiting the experimental tractability of the sponge Amphimedon queenslandica, in which larval release, competence, settlement and metamorphosis are well characterised and readily manipulated ^7,9,18–24. We experimentally induce competent larvae to settle on a naturally co-occurring coralline alga under natural and perturbed light conditions and profile chromatin accessibility and gene expression in the same individuals during the first hours of metamorphosis. This approach allows us to reveal that larvae anticipate impending settlement through diurnal light cues and the deployment of pioneer transcription factors that establish a permissive chromatin landscape prior to environmental induction. Given that sponges diverged from other animal lineages more than 700 million years ago, elucidating these mechanisms provides insight into how deeply conserved gene regulatory programs enabled the evolution of rapid, environmentally induced metamorphosis in the ancestral biphasic animal life cycle.

Results

Dramatic changes in gene expression in the first hour of metamorphosis

Amphimedon queenslandica larvae emerge from adult sponges in the early afternoon and become competent to settle on the coralline alga Amphiroa fragilissima 5-6 hours later, shortly after sunset ^7,18. We have shown previously that diminishing light at sunset is necessary for larvae to respond to the algal cue under normal conditions and to initiate settlement and metamorphosis ¹⁸. The first hour following settlement is characterised by extensive cellular reprogramming and the rapid dissolution of the larval swimming anteroposterior axis (Fig. 1A-C and Supplementary Video 1) ^22–26.

Fig. 1.

To investigate transcriptional and chromatin changes associated with the acquisition of competence and the onset of metamorphosis at sunset, we analysed transcriptomes from precompetent larvae, competent larvae, and postlarvae at 1, 3-4, 6-7 and 11-12 hours post-settlement (hps) (see Methods; Fig. 1D and Supplementary Table 1). Across this time course, more than 8,000 genes were differentially expressed (DESeq2; p-adj < 0.1). Swimming larvae, early postlarvae (1, 3-4 and 6-7 hps) and later postlarvae (11-12 hps) showed transcriptionally distinct states (Fig. 1E-G and Supplementary Tables 2 and 3).

The most pronounced transcriptional shift occurred within the first hour of metamorphosis, immediately following settlement of competent larvae. Of the 4,881 genes differentially expressed from competent larvae to 1 hps postlarvae, 63% were downregulated (Fig. 1G and Supplementary Table 3), and these repressed genes were enriched for larval metabolism (Fig. 1H and Supplementary Table 4). In contrast, genes involved in Wnt, TGF-β, mTOR and FoxO signalling, as well as endocytosis, were transiently upregulated during this interval, coincident with the onset of widespread reprogramming of larval cell types, and consistent with roles in patterning the postlarval and juvenile body plan. By 11-12 hps, gene expression shifted back towards metabolic processes, reflecting establishment of a new postlarval physiological state (Fig. 1H and Supplementary Table 4).

Weighted gene co-expression network analysis (WGCNA) assigned 69% of expressed genes to five major co-expression modules, typified by either up- or downregulation between competent larval and 1 hps postlarval stages (Fig. 1I, Supplementary Fig. 1 and Supplementary Table 5). These modules reinforce the conclusion that the first hour after settlement represents the most transcriptionally dynamic phase, marked by coordinated repression of larval metabolic programs and activation of developmental signalling pathways (Supplementary Fig. 1).

Conserved transcription factors are expressed at exceptionally high levels during metamorphosis

The sharp transcriptional reorganisation that characterises the first hour of metamorphosis is accompanied by a striking and selective induction of transcription factors (TFs). Of the 127 TFs that are differentially expressed during larval development and early metamorphosis, 52 are expressed at levels higher than at least 95% of all expressed genes. These highly expressed TFs span multiple conserved families, including bZIP, bHLH, Tbox, Fox, Sp/KLF, homeobox and Sox (Fig. 1J, K and Supplementary Tables 3 and 6).

All but one of these highly expressed TFs (the exception being the constitutively expressed SRF) share three defining features: they are differentially expressed either during the acquisition of larval competency or within the first hour of metamorphosis; they belong to WGCNA modules that change sharply at the onset of metamorphosis; and they are significantly downregulated by 11-12 hps (Supplementary Table 5). This coordinated temporal pattern is consistent with a transient but potent regulatory phase associated with large-scale reprogramming of larval cell states.

As larvae acquire competence to settle at sunset, 11 TFs are upregulated and 14 are downregulated (Fig. 1J, K and Supplementary Table 6). Among the upregulated and highly expressed TFs are two bHLH factors (CLOCK and HIF-1), the bZIP ATF4/5, two Fox TFs (FoxN2/3 and sponge-specific Fox 2), TboxB and the zinc-finger (ZnF) TF Glass. CLOCK exhibits the most pronounced increase in expression at competence. In parallel, a sponge nuclear receptor (NR1) and multiple bZIP factors are downregulated, including the highly-expressed AP-1 components Fos and Jun, and several of their potential partners (CEBP, CREB, PAR, MAF and NFE2) ^20,27,28 . These TFs, together with most other highly-expressed TFs, are enriched in externally facing larval cell types – ciliated and non-ciliated epithelia, and flask cells – that are implicated in detecting the algal inductive cue and that undergo rapid morphological and functional transformations during the first hour of metamorphosis (Supplementary Fig. 1 and Supplementary Table 7) ^{19,22,23,25,26,29}.

In contrast to the widespread repression of gene expression that dominates the first hour of metamorphosis, the majority (63%) of differentially expressed TFs are upregulated during this period. Of these, 34 rank within the top 5% of all expressed genes (Fig. 1J, K and Supplementary Table 6). Notably, AP-1 and the bZIP factors that are downregulated at competence are now rapidly and strongly induced. In addition, TFs with deeply conserved roles in animal development (multiple Sox, bHLH, Fox, SMAD, Tbox and homeobox members, and NR1) and innate immunity (NF-κB, STAT and IRF) are transiently upregulated at the onset of metamorphosis. The rapid induction and exceptionally high expression of these TFs are consistent with a central role in coordinating environmentally-induced cell state transitions and establishing a robust regulatory environment capable of driving the rapid transformation of the larval body plan into that of the filter-feeding juvenile.

Similar to the dynamics observed for developmental signalling pathways, nearly all of these TFs are transiently upregulated during the first hour of metamorphosis and subsequently downregulated between 6–7 and 11–12 hps. Together, the scale, synchrony and tempo of TF activation and repression are consistent with the genome of competent larvae being poised to respond immediately to the algal inductive cue, enabling a rapid and coordinated transcriptional response that underpins extensive cell reprogramming to transform the larval body plan.

Open chromatin regions are enriched in binding motifs of highly expressed transcription factors

To determine how the A. queenslandica genome responds to the inductive algal cue, we profiled chromatin accessibility in larva and postlarvae using an assay for transposase-accessible chromatin with high-throughput sequencing (ATAC-seq) (Fig. 1D see Methods) ³⁰. We identified 60,716 non-overlapping open chromatin regions (OCRs; ATAC-seq peaks) that were highly reproducible across biological replicates and that overlapped with previously defined larval promoters and enhancers identified by ChIP-seq (Supplementary Fig. 2 and Supplementary Tables 8 and 9) ³¹. Consistent with patterns observed in A. queenslandica embryos and adults, as well as other animals ^32,33, approximately 80% of OCRs are located within protein-coding gene bodies or within 1 kb of transcription start or end sites.

Expressed genes are associated with substantially more OCRs than transcriptionally quiescent genes, with an average 3.5-fold enrichment. Between 67% and 86% of differentially expressed genes and genes belonging to the five WGCNA co-expression modules have proximal OCRs (Fig. 2A, B, Supplementary Fig. 3 and Supplementary Tables 10 and 11). Expressed genes without nearby OCRs may be regulated by more distal cis-regulatory elements or, given the high gene density of the A. queenslandica genome ³⁴, by OCRs associated with neighbouring genes.

Fig. 2.

Previous studies have shown that sponge transcription factors can interact with bilaterian cis-regulatory elements and that sponge enhancers can drive developmental gene expression in vertebrates ^33,35,36. However, to avoid bilaterian-specific annotation biases – given that sponges and bilaterians diverged 700 million years ago – we restricted motif annotation to TF classes and families rather than individual TFs (see Methods). This analysis identified 173 TF binding motifs (TFBMs) corresponding to TF classes and families that are significantly enriched in OCRs associated with genes expressed during larval development and early metamorphosis (Homer p-value < 1e^-10; Fig. 2C and Supplementary Table 10). Of these motifs, 78% correspond to predicted targets of the highly expressed TFs identified above, with an additional 15.6% attributable to highly expressed non-Sp/Klf ZnF TFs (Fig. 2D, E). These include 35 variants of the bHLH E-box motif (CACGTG), 40 homeodomain motifs (including POU, LIM and Antp), 17 Fox, 14 Sp/Klf and 9 bZIP motifs. The enrichment of these 173 TFBMs varies dynamically across developmental stages and co-expression modules, consistent with stage-specific changes in chromatin accessibility during larval development and metamorphosis (Fig. 2D-F, Supplementary Fig. 3 and Supplementary Tables 10 and 11).

Across larval development and early metamorphosis, 3,482 chromatin regions exhibit significant changes in accessibility (differentially accessible chromatin regions; DACRs). More than 56% of differentially expressed genes have OCRs at the same developmental stage, while 4-8% have accessibility one stage prior to the corresponding change in gene expression (FDR < 0.05; Fig. 2B, F, Supplementary Fig. 3 and Supplementary Tables 3, 9 and 12). Notably, nearly one-third of DACRs (1,146) change accessibility between pre-competent and competent larvae, with 71.4% becoming less accessible at competence. This extensive chromatin remodelling contrasts sharply with the relatively modest number of transcriptional changes between these stages (Fig. 1G), and suggests that many chromatin changes precede settlement and may contribute to establishing a permissive genomic state for the rapid transcriptional activation required at the onset of metamorphosis.

Consistent with this interpretation, 62% of all genes that are differentially expressed in 1 hps postlarvae (3032) have proximal chromatin regions already accessible in competent larvae (Supplementary Tables 3 and 9). During the first hour of metamorphosis, 65% of the 1,555 DACRs become more accessible, a proportion closely matching the fraction of genes that are transcriptionally repressed at this stage (63%; cf. Figs. 1g and 2f). This correspondence is consistent with the involvement of rapidly induced, highly expressed TFs in transcriptional repression. It supports a model in which pioneer-like and other highly expressed TFs regulate both the acquisition of competence and the initiation of metamorphosis by interacting with their cognate binding motifs in OCRs proximal to differentially expressed genes.

CLOCK, Jun and Fos OCRs are enriched in binding motifs of highly expressed TFs

Given that more than 90% of TFBMs enriched in larval and postlarval OCRs correspond to highly expressed TFs, we next examined the regulatory architecture of three TFs that appear central to the acquisition of larval competence and the initiation of metamorphosis: the circadian regulator CLOCK ^20,37 and the AP-1 components Jun and Fos ^38–41. CLOCK is significantly and transiently upregulated in A. queenslandica competent larvae at sunset, consistent with a light-dependent role in establishing competence (Fig. 1K and Supplementary Table 6). bHLH-PAS and bHLH (E-box) family motifs consistent with CLOCK:BMAL/CLK-like activity are present in OCRs associated with 165 differentially expressed genes at competence and with 1,121 at 1 hps (Supplementary Table 13). This widespread distribution of bHLH-PAS/E-box family motifs is consistent with pioneer-like regulatory activity associated with CLOCK-related programs that may contribute to establishing a permissive genomic state for rapid transcriptional responses to the algal inductive cue.

In contrast, AP-1 components Jun and Fos, together with other potential bZIP partners (CEBP, CREB, PAR, MAF and NFE2) ^27,28, are significantly and transiently upregulated during the first hour of metamorphosis. At 1 and 6–7 hps, 470 and 340 differentially expressed genes, respectively, have AP-1 or related bZIP TFBMs within associated OCRs, including numerous other TFs and conserved developmental regulators (Supplementary Table 14). In addition, bHLH and bZIP TFBMs are abundant in OCRs associated with genes that remain consistently expressed across these stages (Supplementary Table 10), reinforcing the view that CLOCK and AP-1 occupy central positions in the regulatory networks governing both competence acquisition and environmentally induced cellular reprogramming at the onset of metamorphosis.

In a genomic context, CLOCK is flanked by protein-coding genes that do not show significant upregulation at competence (Fig. 3A, B). Multiple OCRs flanking and within the CLOCK locus are enriched for TFBMs of highly expressed TF families, including Sp/Klf and other zinc-finger, bHLH and Antp motifs (Supplementary Table 15). Notably, the two closest flanking OCRs become less accessible following settlement, coincident with CLOCK downregulation, suggesting that they function as activating cis-regulatory elements. CLOCK, together with many TFs predicted to bind these regulatory regions, is enriched in larval epithelial cell types that express light-detecting cryptochromes ⁴² (Supplementary Table 7), consistent with a role in responding to diminishing light at sunset and being associated with the acquisition of competence.

Fig. 3.

Jun and Fos similarly are flanked by genes with distinct expression dynamics (Fig. 3C-F). OCRs associated with both loci are enriched for binding motifs of TFs that are highly expressed during competence and early metamorphosis (Supplementary Tables 16 and 17). These include multiple bHLH and bZIP motifs, consistent with coordinated activity of bHLH and bZIP family transcription factors during the induction of Jun and Fos at the onset of metamorphosis, and with Jun, Fos and other bZIP factors contributing to their subsequent downregulation at competence and again at 6– 7 hps. With the exception of one flanking OCR at each locus that becomes inaccessible at 6–7 hps, chromatin accessibility in the vicinity of Jun and Fos remains largely stable during larval development and early postlarval stages, suggesting that transcriptional dynamics at these loci are primarily driven by changes in TF availability rather than wholesale chromatin remodelling.

CLOCK and other competence-associated transcription factors are repressed when larvae are prevented from experiencing sunset

To test the role of diminishing light at sunset in regulating gene expression during the acquisition of larval competence, we exposed larvae to constant light beyond sunset, a treatment that inhibits settlement and metamorphosis ¹⁸. Relative to larvae that experienced diminishing natural light, constant light-exposed larvae differentially up- and downregulated 779 and 879 protein-coding genes, respectively (DESeq2; p-adj < 0.1), and activated metabolic pathways that are not normally expressed in swimming larvae (Fig. 4A, B, and Supplementary Table 18). This transcriptional profile is distinct from both precompetent and competent larval states. Importantly, we have previously shown that larvae exposed to constant light beyond sunset swim normally and will rapidly settle and initiate metamorphosis as normal once they are transferred into darkness ¹⁸, suggesting that the constant light phenotype is not due to abnormal physiology or stress; no stress pathways are enriched in larvae exposed to constant light.

Fig. 4.

Strikingly, TF expression was broadly repressed under constant light: 45 TFs were significantly downregulated and only a single TF was upregulated (Fig. 4C). Consistent with this widespread decrease in TF expression, 88% of the OCRs present in larvae under natural light conditions became inaccessible in constant light-exposed larvae (Supplementary Table 18). Together, these results indicate that even brief exposure of larvae to constant light beyond sunset is associated with a novel transcriptional and chromatin state that is distinct from normal developmental trajectories, and extinguishes their capacity to respond to the algal cue, thereby preventing settlement and metamorphosis.

Several transcription factors that are normally upregulated and highly expressed in competent larvae – including CLOCK, HIF, TbxB, FoxN2/3, ATF4/5 and Glass (with CLOCK, TbxB, FoxN2/3, ATF4/5 and Glass ranking in the top 5% of expressed genes, and HIF in the top quartile) – were all significantly downregulated in constant-light–exposed larvae (Fig. 4D). In addition to HIF, the bHLH–PAS transcription factor SIMa, a potential partner of CLOCK, is also highly expressed in competent larvae and significantly repressed under constant light. Consistent with CLOCK and its partners having pioneer-like regulatory properties related to metamorphosis, a large proportion (24%) of the OCRs that close in constant light-exposed larvae have bHLH-PAS and E-box family motifs.

Examination of chromatin accessibility across the CLOCK locus and within 100 kb of its genomic neighbourhood revealed that proximal OCRs associated with CLOCK at competence remain largely unchanged under constant light, whereas six distal regulatory regions exhibit altered accessibility, with five becoming less accessible (Supplementary Table 18). Notably, the sole region that gains accessibility in constant-light larvae lies within an intron of an adjacent gene approximately 2.5 kb downstream of CLOCK (Fig. 4E). These rapid, light-dependent changes in chromatin state support the existence of a conditional gene regulatory program that is normally suppressed at sunset. When activated by aberrant light conditions, this program downregulates CLOCK and other competence-associated transcription factors, disrupts chromatin priming, and prevents the transcriptional and chromatin changes required for larval competence, settlement and metamorphosis (Fig. 4F).

Discussion

Metamorphosis represents one of the most dramatic and consequential transitions in animal life histories, transforming planktonic larvae specialized for dispersal and environmental sensing into benthic juveniles adapted for growth, survival, and reproduction. In Amphimedon queenslandica, this transition is both environmentally induced and remarkably rapid, unfolding within hours of settlement. Such speed implies that key regulatory processes must be initiated before the inductive cue is encountered. Here, by integrating time-resolved transcriptomic and chromatin accessibility profiling across larval competence, settlement, and early postlarval development, we show that larvae undergo extensive regulatory reorganisation while still swimming in the plankton, prior to metamorphic induction. In particular, diminishing light at sunset is associated with widespread changes in chromatin accessibility that precede large-scale transcriptional activation at settlement (Fig. 4F). These findings indicate that the larval genome enters an anticipatory state that is permissive for rapid, coordinated transcriptional reprogramming, rather than responding de novo to the inductive environmental cue.

A defining feature of A. queenslandica metamorphosis is the immediate and substantial activation of a large suite of deeply conserved developmental transcription factors (TFs). To our knowledge, such a diverse suite of TFs reaching expression levels within the top 5% of all genes has never before been observed in animal developmental transitions. Their exceptionally high expression at competence and at the onset of metamorphosis likely confers robustness, ensuring sustained occupancy of cis-regulatory elements and enabling rapid, coordinated shifts in gene expression. Among these, the AP-1 components Fos and Jun, together with their potential bZIP heterodimeric partners ^27,28, are rapidly and strongly induced following settlement. AP-1 is a conserved effector of environmentally responsive signalling pathways, frequently activated downstream of MAPK cascades to regulate cell differentiation, proliferation and death ^38,39. The strong enrichment of bZIP motifs within open chromatin regions associated with genes that change expression immediately after settlement supports a central role for AP-1 in coordinating the transcriptional reprogramming underlying metamorphosis. Consistent with this interpretation, inhibition of MEK/ERK signalling blocks settlement and metamorphosis in A. queenslandica and other marine invertebrate larvae ^8,19,43–45. Thus, although the timing of competence acquisition and the identity of inductive settlement cues are largely species specific, the initiation of metamorphosis itself may rely on a conserved regulatory logic involving MAPK signalling and AP-1–mediated transcriptional activation.

Crucially, the induction of metamorphosis in A. queenslandica is conditional on the prior acquisition of competence, which develops after a period of planktonic swimming that includes exposure to darkness at sunset ^7,18. We show here that competence acquisition is itself accompanied by extensive remodelling of gene expression and chromatin accessibility. These changes not only appear to permit settlement in response to the algal cue, but also establish a permissive chromatin landscape that enables the rapid and large-scale transcriptional response observed during the first hour of metamorphosis. Larvae prevented from experiencing sunset fail to settle and metamorphose, exhibit markedly reduced expression of CLOCK and multiple other highly expressed developmental TFs, and display widespread loss of chromatin accessibility. Together, these findings support a model in which CLOCK and associated pioneer TFs mediate anticipatory chromatin priming that is required for the subsequent induction of AP-1 and other developmental regulators at settlement.

Diminishing light intensity has been implicated in regulating larval settlement in a range of marine species ⁴⁶, raising the possibility that circadian entrainment of competence via CLOCK represents a widespread and ancestral feature of the pelagobenthic life cycle. Notably, A. queenslandica larvae maintained in constant light do not simply remain in a pre-competent state but instead transition into an alternative transcriptional and chromatin configuration that is distinct from both pre-competent and competent larvae. This observation suggests the existence of sensitive, rapidly deployable and conditional gene regulatory programs that are normally suppressed under predictable sunset conditions. Activation of such alternative programs may allow larvae to delay settlement and extend dispersal when environmental conditions are unfavourable, thereby modulating connectivity and colonisation potential.

Together, our results support a model in which environmentally predictable planetary cues, such as the daily light–dark cycle, are integrated during the planktonic larval phase to establish a permissive genomic state that enables rapid induction of metamorphosis upon encounter with a suitable settlement cue. In A. queenslandica, this state is characterised by extensive chromatin remodeling prior to settlement, transient but exceptionally high expression of conserved developmental transcription factors at competence and induction, and enrichment of their binding motif families in regulatory regions associated with genes activated during early metamorphosis. Preventing larvae from experiencing sunset does not simply delay development, but instead diverts them into an alternative transcriptional and chromatin configuration that is incompatible with settlement and metamorphosis, highlighting the conditional nature of this regulatory program. Although functional dissection of individual regulators is currently unfeasible in this system, the temporal ordering of chromatin accessibility changes, transcription factor expression, and environmental induction is consistent with a model in which light-entrained transcriptional activity—potentially involving CLOCK-associated and AP-1/bZIP family factors— contributes to chromatin priming and rapid transcriptional reprogramming. We propose that such anticipatory chromatin regulation may represent a general solution to the problem of speed and robustness in metamorphosis, enabling biphasic animal life cycles to integrate predictable planetary signals with species-specific ecological cues across evolutionary timescales.

The existence of an alternative, inducible regulatory program that is conditional on the environmental experience of the sponge larva provides a potential mechanistic explanation for how larvae modulate dispersal range and timing of settlement. The magnitude and tempo of these regulatory changes are comparable to the rapid metamorphic transitions observed across many marine invertebrates ^6,16, suggesting that the genomic control mechanisms identified here may be broadly relevant. However, in the absence of equivalent chromatin-resolved datasets from parahoxozoan lineages (e.g. echinoderms, hemichordates, annelids, molluscs, bryozoans or cnidarians), the extent to which these mechanisms are conserved remains an open question. We propose that, despite the extraordinary diversity of larval forms and settlement cues across the animal kingdom, metamorphosis is governed by a shared set of genomic control principles that integrate planetary and species-specific environmental signals through deeply conserved transcriptional regulators.

Methods

Larval and postlarval collection and cell isolation

Adult Amphimedon queenslandica were collected from the intertidal reef flat of Shark Bay, Heron Island Reef, Southern Great Barrier Reef (23° 27’S, 151° 55’E) and were immediately transferred to aquaria at Heron Island Research Station, where they were maintained in flow-through, unfiltered ambient seawater drawn from the adjacent reef flat ⁴⁷. Naturally emerging larvae were collected in the early to mid-afternoon, and maintained in ambient conditions as previously described ¹⁸. Randomly-selected larvae were transferred at sunset (5-6 hours post-emergence; hpe) to a 6-well dish with 10 ml of 0.22 µm filtered seawater (FSW) (10 per well), and then exposed 30 min later to three small branches of the coralline algae Amphiroa fragilissima as previously described ¹⁸. These dishes were maintained under ambient water temperature and light conditions ¹⁸. For collection of postlarval stages, dishes were checked every 30 mins to determine the time of settlement, defined by the larva remaining attached to the alga when the dish was given a gentle shake.

Six randomly-sampled larvae and postlarvae were collected from each of the following stages for cell dissociation, and CEL-seq2 and ATAC-seq analyses: precompetent larvae (0-1 hpe); competent larvae (5-6 hpe); 1 hps postlarvae; 3-4 hps postlarvae; 6-7 hps postlarvae; and 11-12 hps postlarvae (Fig. 1D). All individuals were treated individually (i.e. not pooled). All postlarvae that had settled on A. fragilissima were removed rapidly from the alga with a sharp tungsten needle and transferred into FSW on ice. Precompetent and competent larvae were treated the same but did not need to be removed from the alga. All individual larvae and postlarvae were washed three times for 2 min in ice cold FSW and then centrifuged through a 25 µm metal mesh fitted into a 1.5 ml Eppendorf at 500 rpm for 10 min at 4°C. Pelleted fragments and cells were resuspended and washed with ice cold 0.22 µm filtered calcium-magnesium-free artificial seawater (CMFSW; 449 mM NaCl, 33 mM Na₂SO₄, 2.15 mM NaHCO₃, 9 mM KCl, 10 mM Tris-HCl pH 10, and 2.5 mM EGTA) for 5 min and recentrifuged at 500 rpm for 10 min at 4°C. The cell pellet was resuspended in 20 µl ice cold CMFSW. This procedure took less than 25 minutes at 4°C, thus minimising changes in gene expression or chromatin state associated with these perturbations.

For four developmental stages – precompetent and competent larvae, and 1 and 6-7 hps postlarvae – we separately assayed both gene expression and chromatin accessibility in the same six individuals by splitting their resuspended cells into two pools. For each individual of these four stages, 5 and 15 µl of the suspended cells were used for CEL-seq2 and ATAC-seq library construction, respectively. For CEL-seq2, 5 µl of the cell suspension was immediately mixed with 50 μl ice-cold TRIZol and stored in -80°C. For ATAC-seq, 15 µl of the cell suspension was immediately mixed with 50 μl ice cold ATAC-seq cell lysis buffer and processed as described below. The other two stages – 3-4 and 11-12 hps postlarvae – were assayed only for gene expression, so that total 20 µl of the cell suspension was immediately mixed with 200 μl ice-cold TRIZol, resulting in the same dilution as above, and stored in -80°C.

Imaging early metamorphosis

Images of A. queenslandica postlarvae were taken with a Nikon SMZ25 microscope using the program Nikon Elements and 6X zoom every 3 min, starting from ∼10 min after the initiation of settlement on A. fragilissima through to ∼ 2 hps. Each time point consisted of a z-stack with 14 images at an interval of 25 µm, which were then stacked into a single image using the focus stacking program Zerene Stacker with the P Max algorithm ⁴⁸. These were cropped, brightened and denoised, then unsharp masked using Adobe Photoshop. To generate the timelapse movie, the frames were compiled into an avi using ImageJ and converted to mp4 (or mkv) using AnyVideoConverter (https://www.any-video-converter.com/en8/for_video_free/).

CEL-seq2 library construction, sequencing and analysis

RNA isolation and quality-assessment, and CEL-seq2 library construction all were performed as previously described ^49–51. Libraries were sequenced on two lanes of Illumina HiSeq-PE 150 platform at Novogene Co., China and raw sequencing data were processed using a publicly available pipeline (https://github.com/yanailab/celseq2) ⁵⁰ with a built-in a suite of python-based scripts for demultiplexing, mapping to the genome, and counting reads. Raw reads were mapped to A. queenslandica Aqu3.1 genome (NCBI Bioproject PRJNA668660) ⁵² using Bowtie2 ⁵³, and assessed for coverage and quality as previously described ^18,51 (Supplementary Table 1). Count tables for the two sequencing lanes were merged using a custom R script (https://github.com/hfyuanuq/thesis.scripts/blob/main/CEL-Seq2-combine%202%20lanes) before the differential expression analysis. CEL-seq2 raw counts from the two larval ages in two different light regimes were obtained from https://github.com/tahshasay/TS1015_DESeq2/tree/master ¹⁸ and larval single-cell RNA-seq raw counts were obtained from https://github.com/tanaylab/ ²⁹.

Analysis of differentially expressed genes

Only genes with at least 4 read counts in at least 4 samples were considered expressed; those not meeting these criteria were filtered out from all subsequent analysis. PCA of all expressed genes was performed on variance-stabilizing transformation (VST) transformed counts and visualized using the function ‘plotPCA’ to examine sample-to-sample distances and identify potential outliers with 95% confidence ellipses. Differential gene expression was conducted in R ⁵⁴ using the Bioconductor package DESeq2 ⁵⁵. Pairwise comparisons were conducted between each of the six developmental stages to produce a differentially expressed gene (DEG) list for each pair of stages, and an extra comparison was conducted between 1 hps and 6-7 hps, with a False Discovery Rate threshold (FDR) < 0.1 (Supplementary Table 3). This FDR allowed us to detect substantial changes in transcript abundance within one hour (that is, between competent larvae and newly settled 1 hps postlarvae), taking into account known eukaryotic rates of RNA transcription and degradation ⁵⁶. Venn diagrams were generated using the R package ggVennDiagram ⁵⁷. Heatmaps were produced in R using packages pheatmap (https://CRAN.R-project.org/package=pheatmap) and RColorBrewer (https://CRAN.R-project.org/package=RColorBrewer) to visualize expression profiles with rlog normalised counts for CEL-seq2 datasets. Single-cell RNA-seq reads were normalised with the function CPM (cpm(data, log=TRUE)) using the Bioconductor package edgeR (version 3.14.0) implemented in R ⁵⁸. Alluvial (Sankey) diagrams were generated using SankeyMATIC (vBeta https://sankeymatic.com/build/).

Co-expression analysis

Weighted Gene Co-expression Network Analysis (WGCNA; version 1.61 ⁵⁹ was performed on the VST-transformed developmental dataset using R (version 4.0) ⁵⁷. The soft-thresholding power was determined using the powerEstimate function. Modules were identified via the dynamic tree cut method with a minimum module size of 30 genes. To retain the distinct expression patterns of each module, a merging distance threshold of 0.1 was applied.

ATAC-seq library preparation, sequencing and analysis

ATAC-seq libraries were prepared as previously described ³². Briefly, we added ATAC-seq cell lysis buffer (10 mM Tris-HCl, pH 7.4, 10 mM NaCl, 3 mM MgCl₂, 0.1% v/v IGEPAL CA-630) to the 15 µl of dissociated cells (15-20,000 cells) and centrifuged at 500 g for 10 min at 4°C. The pellet was resuspended in 50 μl transposition reaction and incubated at 37°C for 30 min. The resultant DNAs were amplified, purified and quality assessed as previously described ³². Libraries were sequenced on a 2 x 101 bp NovaSeq SP flow cell platform. Raw reads were trimmed and low-quality reads removed using FastQC and Trim-Galore (parameters: -q 30 --phred33 --stringency 4 --length 20 -e 0.1) (https://github.com/FelixKrueger/TrimGalore). Reads were mapped to the A. queenslandica Aqu3.1 genome using Bowtie2 (parameters: -X 1000) ⁵³. Paired reads were identified using SAMtools (parameter: -f 2) ⁶⁰ and duplicated reads were identified using Sambamba and removed ⁶¹. MACS2 was used to call peaks (parameters: -g 1.2e8 --nomodel --shift -100 --extsize 200) ⁶². This yielded over 46 million reads that corresponded to over 20,000 peaks from each library (Supplementary Table 8). ATAC-seq library quality was determined by calculating fraction of reads in peaks (FRiP) scores using the ENCODE ATAC-seq pipeline (https://www.encodeproject.org/atac-seq/); a FRiP score over 0.2 was deemed acceptable (Supplementary Table 8). Peaks from the three replicate ATAC-seq libraries for each stage were pooled and overlapping peaks were intersected using bedtools ⁶³, also in the ENCODE ATAC-seq pipeline (https://www.encodeproject.org/atac-seq/).

Analysis of open chromatin regions

ATAC-seq peaks were annotated using the Bioconductor package ChIPseeker ⁶⁴, with the promoter region set as ± 500 bp (Supplementary Table 9). The Bioconductor package DiffBind (V2.16.0) was used for differential binding analysis of peaks ⁶⁵. Significantly different peaks (FDR < 0.05) between the four developmental stages were identified by pairwise comparisons on successive stages using DESeq2 ⁵⁵. Consistent peaks and overlapping peaks were generated in DiffBind using the function dba.peakset (parameters: consensus=DBA_FACTOR, minOverlap=0.33) ⁶⁵. Library correlation analysis was performed using ChromVAR ⁶⁶ to determine peaks across all 12 libraries. Genes associated with these peaks were compared to the list of genes associated with ChIP-Seq peaks ⁶⁷.

Motif identification, annotation, clustering, and accessibility analysis

Homer ⁶⁸ was used to identify TFBMs with parameters of -size 500 -len 6, 8, 10, 12 -h. All identified motifs were combined into a single motif list, and we excluded all the potential A. queenslandica specific, sequence bias motifs, and non-significant identified motifs (p-value > 1e^-10). Motifs were matched to the database (https://jaspar2022.genereg.net/downloads/) that we customised by removing fungi- and plant-specific motifs (available at https://github.com/hfyuanuq/atac-pub) using MEME-TomTom (version 5.5.8) with default threshold 0.5 ⁶⁹. After motif matching, all putative TFBMs were assigned with JASPAR motif identifiers. Then motifs were clustered into transcription factor classes and families based on JASPAR annotations using the Bioconductor package TFBSTools ⁷⁰.

TF binding sites analysis was performed using the GimmeMotifs scan tool ^71,72. OCR sequences of expressed genes and differentially expressed genes were scanned against the JASPAR non-redundant database. Motif scan significance was determined using log-odds scores calculated relative to the A. queenslandica genomic background model. Putative motif occurrences were filtered using a FDR threshold of 0.01. Following the scan, all putative TFBMs were clustered into transcription factor class and family as previously described. Chromatin accessibility and peak profiles were visualized using the R package Gviz ⁷³ and the Integrative Genomics Viewer (IGV) software ⁷⁴.

KEGG analysis

All the KEGG enrichment analyses were performed using the online tool STRING (https://string-db.org/) with default settings.

Light perturbation experiment

From adults collected maintained in flow-through seawater aquaria at Heron Island Research Station, as described above, > 50 naturally-released larvae were collected within 30 min of emerging, and split equally among two 500 ml glass beakers of sea water. One beaker was exposed to natural light and the other to constant light, both at ambient temperature, as described previously ¹⁸. At sunset (5–6 hpe), eight larvae from each beaker were randomly collected. Their cells were dissociated and CEL-seq2 and ATAC-seq libraries were constructed, sequenced and analysed as described above ⁷⁵.

CEL-seq2 and ATAC-seq sequencing were performed using the Illumina NovaSeq X Series sequencer (PE150, NovaSeq Control Software v1.2.2.48004 and Real Time Analysis v4.6.7). Precompetent and competent larval, and constant and natural light CEL-seq2 datasets were combined and integrated, with competent larval and natural light samples being combined into a single group. Batch effects were corrected using ComBat-Seq ⁷⁶. Differential OCRs and associated TFBMs were identified, characterised and visualised as described above. For each motif, the average deviation difference between constant and natural conditions was computed to determine the direction and magnitude of activity change.

Data availability

All data generated or analyzed during this study have been deposited in NCBI BioProjects database. CEL-seq2 and ATAC-seq data for normal light condition are available at accessions PRJNA1162246 and PRJNA1161795, respectively, and for constant light condition at accessions PRJNA1394668 and PRJNA1394227, respectively. Source data files have been provided for all figures as Supporting Files.

Acknowledgements

This study was supported by funds from the Australian Research Council to S.M.D and B.M.D (DP210100703 and DP230102109). Biological samples were collected with the support of the Heron Island Research Station and Bin Yang, and ATAC-seq was performed under guidance from Miloš Tanurdžić. Computational analyses were enabled by support from the Queensland Cyber Infrastructure Foundation, with a special thanks to Nick Rhodes.

Additional information

Author contributions

S.M.D, B.M.D. and H.Y. conceived and designed the project. H.Y. isolated cells and prepared all libraries. H.Y. performed all gene expression and chromatin state analyses. O.B. and Z.P. visually documented postlarval morphogenesis. H.Y, S.M.D. and B.M.D wrote the manuscript with minor contributions from the other authors.

Materials & Correspondence

Materials & Correspondence should be addressed to S.M.D or B.M.D.

Funding

Australian Research Council (DP210100703)

• Bernard M Degnan

• Sandie M Degnan

Australian Research Council (DP230102109)

• Bernard M Degnan

• Sandie M Degnan

Additional files

Supplementary figures

Supplementary Table 1. Information on CEL-seq2 libraries.

Supplementary Table 2. Normalised and average CEL-seq2 counts for all expressed genes using DESeq2.

Supplementary Table 3. Genes that are differentially expressed during larval development and early metamorphosis.

Supplementary Table 4. KEGG enrichments for genes that are upregulated at each larval and early postlarval stage.

Supplementary Table 5. KEGG enrichments for genes in each WGCNA modules.

Supplementary Table 6. TF genes that are expressed and differentially expressed during larval development and early metamorphosis.

Supplementary Table 7. Raw and normalised MARS-seq counts 29 of differentially expressed TF genes in larval cell types.

Supplementary Table 8. Information on ATAC-seq libraries: sequencing reads, mapping ratios, peaks, and fractions of reads in called peak regions.

Supplementary Table 9. Genomic location and TFBMs enrichment in OCRs present during larval development and early metamorphosis.

Supplementary Table 10. OCRs and their TFBMs associated with expressed genes at all four stages.

Supplementary Table 11. OCRs and their TFBMs associated with WGCNA modules.

Supplementary Table 12. OCRs and their TFBMs associated with DEGs between precompetent and competent, and competent and 1 hps.

Supplementary Table 13. OCRs and their PAS related (bHLH) TFBMs associated with differentially expressed genes during larval development and early metamorphosis.

Supplementary Table 14. OCRs and their Fos/Jun-related (bZIP) TFBMs associated with differentially expressed genes during larval development and early metamorphosis.

Supplementary Table 15. OCRs and TFBMs found in and around CLOCK during larval development and early metamorphosis.

Supplementary Table 16. OCRs and TFBMs found in and around Jun during larval development and early metamorphosis.

Supplementary Table 17. OCRs and TFBMs found in and around Fos during larval development and early metamorphosis.

Supplementary Table 18. Comparison of gene expression, chromatin state and accessible TFBMs in larvae exposed to natural or constant light.

Supplementary Video 1. Time lapse video (mp4) of the first two hours of A. queenslandica metamorphosis after settling on an algal fragment of A. fragilissima. The posterior