Bugulina neritina adult colonies were collected from fishfarm rafts in the Daya Bay, Shenzhen, P. R. China (22°33′53.2″N, 114°32′09.0″E). The adult colonies were kept in dark, transported to the laboratory and maintained in seawater aquarium. To keep the activity of adult colonies, the container was constantly aerated and about one third of seawater was replaced daily. After applying bright light illumination for 30 min to the adult colonies, the released larvae were collected for the settlement bioassays. Butenolide powder was solubilized in 100% DMSO. Before the experiment, butenolide (BU) solutions were prepared in filtrated seawater (FSW) with the final concentrations of 2% DMSO and 2.5, 5, 10, 20, 40 μg mL^-1 of BU. Five hundred microliters of seawater containing about 10 to 15 larvae were added to each of the well of a 24-well polystyrene plates, followed by the addition of an equal volume of BU solution. In total, seven treatment groups, namely FSW control, 1% DMSO FSW control, 1.25, 2.5, 5, 10 and 20 μg mL^-1 BU were prepared in quadruplicate. The 24-well plates were kept in dark at room temperature (about 25°C) and the number of settled larvae was recorded after 3 hours using a dissecting microscope (Leica EZ4, Heerbrugg, Switzerland). Larval mortality was recorded at 3, 9, 24 hours. The full set of experiment was repeated three times using three different batches of B. neritina larvae released from different adult colonies.

To determine the tolerability and efficacy of the concentration of BU solution during treatment, recovery assays were performed. After adding 20-30 larvae (500 μL in volume) to each well of the 24-well plates, equal volume of BU stock solutions was applied to the final concentration of 1.25, 2.5, 5, 10, 20 μg mL^-1 with the seawater and DMSO (1%) as control. The plates were placed in dark for 3 hours and then the original solutions were replaced with same volume of fresh seawater twice. The numbers of settled larvae were recorded at 24 hours after incubation in dark.

To analyze the effect of BU treatment on larval behavior, three 10 seconds 12 frame per second (fps) video footages were taken using a digital camera (Olympus DP72, Tokyo, Japan) coupled with a stereo microscope (Leica MZ16, Heerbrugg, Switzerland) on each of the six replicates in the control (seawater) group and 5, 10 and 20 μg mL^-1 BU-treated groups at different time points. For the control group, video footages were taken at 1.5 hours after the initiation of the settlement bioassay; for the 5, 10 and 20 μg mL^-1 BU-treated groups, video footages were taken at 1.5, 6 hours and 1 hour after recovery from a 3 hours BU treatment. From each video footage, all frames were extracted into.jpeg format image using a python script (provided in Supplementary File 1 ). The extracted frames were then imported into Adobe Photoshop (version 2020). The initial frame image of a footage was imported firstly as background and all subsequent frame images were adjusted to dim-light mode and then overlaid on top of the initial frame image. Larval swimming behavior was then summarized based on the trajectory of each larva, as visible in the resulting overlaying image of each footage. Larvae with a circular swimming trajectory were concluded as exhibiting looping swimming behavior; larvae that exhibited spiral swimming trajectory were concluded as exhibiting forward looping behavior; larvae swimming trajectory pattern that did not match neither of the above swimming behavior were determined to be typical swimming behavior.

Five control larval biological replicate samples, four replicates of 5 μg mL^-1 and 10 μg mL^-1 butanolide treated samples were prepared for RNAseq experiment. To prepare RNAseq sample, approximately 300 larvae were treated with 5 or 10 μg mL^-1 butenolide and the FSW was used as negative control. After incubation in dark for 3 hours, the unsettled larvae were collected for RNAseq analysis, and all samples were frozen with liquid nitrogen immediately and then kept in -80°C until RNA extraction. All treatments were repeated four times.

Total RNA of the larval samples was extracted with TRIzol (Bioer, Hangzhou, China) and the quality of the total RNA was tested using Bioanalyzer (Aglient, Santa Clara, CA, USA). After enrichment and fragmentation of mRNA, cDNA library construction and normalization were conducted using the Illumina Truseq v3 kit according to the standard protocol. The cDNA libraries generated from each sample were sequenced using Illumina HiSeq 2500.

All raw read data generated from the RNAseq of five control larval biological replicate samples, four replicates of 5 μg mL^-1 and 10 μg mL^-1 butanolide treated samples were deposited to the NCBI Short Read Archive (SRA) with the BioProject accession number PRJNA855486. Clean read data were firstly assembled using Trinity v2.6.2 (Haas et al., 2013) with the following specific parameters (–min_glue 3 –trimmomatic –max_memory 100G –group_pairs_distance 250 –path_reinforcement_distance 85 –min_kmer_cov 3). The obtained Trinity assembly was further clustered using Corset (Davidson and Oshlack, 2014) based on the result of read mapping and hierarchical clustering of the gene expression pattern of each Trinity transcript (unigene) across all the 12 samples. The clustered assembly was further collapsed at 95% sequence identity using CD-HIT-EST (Li and Godzik, 2006), generating the final B. neritina transcriptome assembly.

The final assembly was then annotated by BLASTx (NCBI-blast+ v.2.10.0) against NCBI non-redundant protein database (nr), NCBI nucleotide database (nt), EuKaryotic Orthologous Groups (KOG), Swiss-Prot, Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO). Protein domain assignment was also performed using HMMer against the Pfam HMM database (Finn et al., 2011). E-value cutoff was set as 1e^-5 for homology search by BLASTx and for protein domain assignment by HMMer. In parallel, Illumina read mapping was performed using Bowtie2 (Langmead and Salzberg, 2012) and subsequent transcript abundant estimation was performed using Salmon (Patro et al., 2017). Transcript expression level was normalized to transcripts per million reads (TPM). Differential expression analysis was performed using DESeq2 (Love et al., 2014). Two pairwise comparisons, 5 μg mL^-1 treatment vs control and 10 μg mL^-1 treatment vs control, were performed. The threshold for differential expression was set at abs(log₂(FoldChange))>1 (or at least 2 times differences) and adjusted p-value < 0.05. The differentially expressed unigenes were selected and further clustered by hierarchical clustering and k-mean clustering. Functional enrichment analysis of KEGG pathway was performed using KOBAS (2.0) (–fdr BH) (Xie et al., 2011), and the threshold for over-representation of a KEGG pathway was set at adjusted p-value < 0.05.

At 3 hours post-BU treatment, larval settlement rates at 10 μg mL^-1 or higher concentrations were significantly lower (ANOVA p<0.001) than the FSW or 1% DMSO FSW control, indicating that BU have a strong and significant inhibitory effect on B. neritina larval settlement ( Figure 1A ), which is in accord with our previous studies (Zhang et al., 2012b). While more than 80% larvae have settled in the control, 1.25 and 2.5 μg mL^-1 BU treatment groups in the first three hours of the bioassay experiment, persistent inhibitory effect on larval settlement was observed in 10 and 20 μg mL^-1 BU treatment ( Figure 1B ). No mortality was observed in all treatment or control groups in the first 3 hours of treatment. However, after 9 hours of BU treatment, more than 74% (± 19%) mortality was observed in 20 μg mL^-1 BU treatment group ( Figure 1C ). After 24 hours of BU treatment, all concentrations led to different degrees of larval mortality. Almost all larvae were dead in 10 and 20 μg mL^-1 BU treatment groups. But in lower BU concentration treatments, only the 5 μg mL^-1 BU treatment group resulted in significantly higher larval mortality rate than the seawater control group (ANOVA p<0.001) ( Figure 1C ). The mortality rates at 2.5 μg mL^-1 and lower BU concentration treatments were not significantly different to the SW or DMSO control group.

In the recovery bioassay, BU solution was replaced by seawater after 3 hours BU treatment and the rate of larvae completed metamorphosis in the recovered wells was observed after 24 hours (3 hours BU+ 21 hours seawater) and compared to that in wells that experienced prolonged 24 hours BU treatment (no recovery). As expected, larvae in the control groups, 1.25 μg mL^-1 and 2.5 μg mL^-1 BU treatment groups completed metamorphosis within 24 hours in both the recovery and prolonged BU treatment experiments; larvae in 10 and 20 μg mL^-1 BU treatment groups could not complete metamorphosis in both the recovery and prolonged BU treatment experiments, suggesting that larvae were unable to recovered from a 3 hours 10 μg mL^-1 or higher BU concentration treatment; larvae in prolonged 5 μg mL^-1 BU treatment group could not complete metamorphosis, but removal of BU at 3 hours allowed 87% (± 18%) BU-treated larvae to recover and complete metamorphosis ( Figure 1B ).

To analyze the effect of BU treatment on B. neritina larval swimming behavior, we took video footages of B. neritina larvae in seawater (SW control), in BU treatments and in recovery experiments (replacement of BU solution with SW). By analyzing the larval movement trajectory, we observed circular looping (see Figure 2A ) and spiral swimming behavior (see Figure 2B ) among larvae treated with 5 μg mL^-1 or higher concentrations of BU, but seldom in 1.25 and 2.5 μg mL^-1 BU-treated or control larvae, which typically exhibited movement with large displacement distances (see Figure 2C ).

At 1.5 hours, the percentage of larvae that exhibited circular looping or spiral swimming behavior in 5, 10 and 20 μg mL^-1 BU treatment (55% ± 12%, 70% ± 20%, 68 ± 10%) groups were significantly higher than in 1.25 and 2.5 μg mL^-1 BU-treated [3 ± 8% and 4 ± 4%)] larvae and the control groups (SW:2 ± 2%; DMSO 3 ± 5%) ( Figure 2D ).

At 6 hours, more than 95% larvae in the control groups (SW and 1% DMSO), 1.25 and 2.5 μg mL^-1 BU treatment groups already settled and initiated metamorphosis and hence no video was taken for these experimental groups. Similar percentage of unsettled larvae in 5, 10 and 20 μg mL^-1 BU treatment groups (67% ± 8%, 59% ± 21%, 48 ± 18%) exhibited looping or spiral swimming behavior as at the 1.5 hours time point ( Figure 2E ).

In the recovery assay, after replacement of BU (after 3 hours treatments) by SW for 1 hour, the percentage of larvae exhibiting looping or spiral swimming behavior in 10 and 20 μg mL^-1 BU-treated groups (67% ± 8%, 63% ± 9%) remained at a similar level as at 1.5 and 6 hours. In contrast, percentage of larvae exhibited looping or spiral swimming behavior in the recovery assay was significantly lower in the 5 μg mL^-1 BU treatment group (6 ± 4%) ( Figure 2F ).

Based on the above bioassay results, we define 5 μg mL^-1 and 10 μg mL^-1 as “low” and “high” BU concentration (low [BU], high [BU]) treatments. To understand the influence of these BU treatments on larvae, four low and four high [BU] treated larval samples, as well as five untreated “control” larval samples were submitted for cDNA library preparation and HiSeq sequencing. After quality trimming, these 13 samples generated 635,547,702 paired-end reads, composing of 95.33 Gbp. De novo assembling of the RNAseq data of all the 13 samples generated a transcriptome assembly containing 138,636 transcripts, which were further clustered into 48,001 unigenes based on the result of read mapping and transcript abundance using the Corset algorithm (Davidson and Oshlack, 2014). The transcript length N50 of the assembly is 1,893 bp and the average transcript length is 1,214 bp ( Figure S1 ).

Among the 48,001 unigenes, 22,298 (or 46.45%) were annotated (by BLASTx, e-value cutoff at 1e-5 for homology search, OR, by HMMer for protein domain assignment) by at least one of the public databases, including NCBI nucleotide (nt) sequence database, the NCBI non-redundant (nr) database, SwissProt, Pfam, Gene Ontology (GO), EuKaryotic Orthologous Groups (KOG), KEGG Orthrology (KO) and the KEGG database ( Figure S2 ). Consistent to our previous RNAseq study (Xu et al., 2020), close to 22% of the nr-annotated unigenes were best-matched (share the highest sequence homology) to the brachiopod Lingula anatina ( Figure S3 ), once again indicating a potential close phylogenetic relationship between bryozoans and brachiopods, as suggested by an earlier study (Luo et al., 2015).

The functional annotation of the unigenes exhibited different tendency in different functional annotation methods. GO classification of the unigene highlighted the dominance of genes relating to “Binding” and/or “Catalytic activity”, in terms of molecular function; “Macromolecule complex”, in terms of Cellular Component; “Signaling” and “Regulation of biological process” in terms of Biological Process ( Figure S4 ). On the other hand, KEGG classification indicated the dominance of genes relating to signal transduction and translation ( Figure S5 ).

The transcriptome profiles of high [BU] and low [BU] treated larvae exhibited major differences comparing to the control untreated larvae, as revealed by the Pearson correlation coefficient analysis of the gene expression profile of each sample replicates. However, high [BU] and low [BU] treatments have led to a very similar larval transcriptome profile, as indicated by the outstanding Pearson correlation coefficient (r² = 0.845 to 0.911) between high [BU] and low [BU] treated larval samples ( Figure S6 ).

To further characterize the effect of BU on larval gene expression, the larval gene expression levels at high [BU] and low [BU] treatments were compared to that of the control. 8,040 genes were differentially expressed in both high [BU]-to-control and low [BU]-to-control comparisons while only 2,995 were specific to high [BU]-to-control ( Figure S7 ) and 1,757 DEGs to low [BU]-to-control comparisons ( Figure S8 ).

To further understand the larval response to BU treatment, the functional enrichment analysis was performed based on the over-representation of enriched KEGG and GO terms among the list of DEG in high [BU]-to-control and low [BU]-to-control comparisons. In the GO enrichment analysis, the GOBP term “DNA integration” and the GOMF term “ATP activity” were both over-represented in both high [BU]- and low [BU]-to-control DEGs ( Figures 3A, B ), and the fold of enrichment of both GO terms were more significant in the up-regulated DEGs than in the down-regulated DEGs. The GOBP term “xenobiotic metabolic process”, which is mainly related to the metabolism of foreign compounds, or in other words, detoxification, was over-represented in the list of up-regulated DEGs in the low [BU]-to-control comparison ( Figure 3A ).

The list of enriched KEGG terms among the DEG in the low and high [BU]-treated larvae were highly similar. However, the enriched KEGG terms in the up- and down-regulated DEGs were remarkably different. For the up-regulated DEGs, the nervous system related term “Serotonergic synapse”, the detoxification related term “ABC transporter”, the lipid biosynthesis related term “fatty acid biosynthesis”, and the development pathway “Notch signaling pathway” were over-represented ( Figures 3C, D ); for the down-regulated DEGs, some of the essential functional terms such as “Cell cycle”, “Proteasome” and “Ribosome biogenesis in eukaryotes” were over-represented ( Figures 3E, F ). The results of GO and KEGG enrichment analyses indicated both low and high [BU] treatments suppressed de novo synthesis, cell proliferation and developmental processes, and at the same time triggered detoxification and increased energy production.

It is worth to note that B. neritina is a non-model organism and some of the differentially expressed genes that have important functions may not be included in the GO or KEGG functional categories. Therefore, we decided to investigate unigenes involved in nitric oxide signaling and other important functional categories that were implicated in larval settlement.

The enrichment of the KEGG term “serotonergic synapse” among up-regulated DEGs in both low and high [BU]-treated larvae suggested the effect of BU treatment on larval nervous system. The list of up-regulated DEGs relating to “serotonergic synapse” included the synaptic enzyme aromatic-L-amino-acid decarboxylase (AADC, or DOPA decarboxylase) which catalyzes the formation of dopamine from L-DOPA, neurotransmitter receptors Gamma-aminobutyric acid receptor (subunit beta) and 5-hydroxytryptamine receptor. The selected list also included a number of intracellular components directly relating to intracellular calcium signaling, including inositol 1,4,5-triphosphate receptor, protein kinase A, phosphatidylinositol phospholipase C and various types of voltage-dependent calcium channels ( Tables 1 , S1 ). The above unigenes were up-regulated by 1.2 to 4.3 folds in low and/or high [BU] treated larvae ( Table S1 ).

Five and six nitric oxide synthase (NOS) unigenes were up-regulated in low and high [BU]-treated larvae, with a magnitude of up-regulation ranged from 1.5- to 2.9-fold. In addition to the up-regulation of NOS unigenes, the unigenes encoding downstream components cGMP-dependent protein kinase (PKG), soluble guanylate cyclase (sGC) and guanylate cyclase (GC) were also up-regulated in BU-treated larvae. The soluble guanylate cyclase unigene (Cluster-12452.419) was found to be up-regulated by 1-fold only in low [BU]-treated while the guanylate cyclase unigene (Cluster-13761.0) was up-regulated by 3.5- and 3.1-fold in low and high [BU]-treated larvae; the unigene of PKG was up-regulated by 1.4 and 1.5 fold in low and high [BU]-treated larvae. The expression of the unigene of arginine kinase, which uses the same substrate (L-arginine) as NOS and hence affects one of the upstream components the NO signaling pathway, was not found to exhibit significant variation (data not shown).

The developmental KEGG term “Notch signaling process” was over-represented among up-regulated DEGs, implying that developmental processes in the larvae were likely to be affected. We therefore investigated in detail the expression pattern of developmental pathway genes. Surprisingly, except Notch, which was up-regulated by almost two folds in both low and high [BU] treated larvae, other major developmental morphogens were all significantly down-regulated by BU treatment. Wnt ligands wnt10 (Cluster-12452.10908), wnt16-like (Cluster-9233.0), an WNT antagonist sfrp-5 (Cluster-12452.10592), and the Wnt pathway key player beta-catenin (beta-catenin-like protein 1, Cluster-12452.22208) were either not detected in low [BU] or significantly down-regulated in both high [BU] and low [BU] treated larvae. In addition, the BMP ligand BMP5 (Cluster-12452.1991), the Shh ligand ihh (Cluster-12452.18656), were also not detected in low [BU] or significantly down-regulated in the high [BU] treated larvae ( Table 1 ).

Proteasome is one of the over-represented KEGG term in the list of down-regulated DEGs. In fact, there were 29 proteasome DEGs down-regulated in both low and high [BU]-treated larvae and one proteasome DEG down-regulated exclusively in high [BU] treated larvae. These proteasome DEGs were annotated as different subunits, regulatory subunits or non-ATPase regulatory subunit of 20S and 26S proteasomes, proteasomes inhibitor subunits, or proteasome assembly chaperone. The magnitude of down-regulation ranged from 1- to 2.18-fold.

Although the KEGG term “fatty acid biosynthesis” was enriched among up-regulated DEGs in both low and high [BU]-treatment groups, it was rather surprising to find that only two lipid metabolic proteins acetyl-CoA carboxylase and long-chain-fatty-acid–CoA ligase were up-regulated. While acetyl-CoA carboxylase catalyzes the production malonyl CoA, the building block for fatty acid biosynthesis, the unigene of fatty acid synthase (Cluster-12452.14482) was down-regulated by 1- and 2-fold in low and high [BU]-treated larvae; the unigene of elongation of very long chain fatty acids protein (Cluster-12452.16696, Cluster-12452.7219) was down-regulated by 1.2-fold in high [BU]-treated larvae; the unigene of fatty acid 2-hydroxylase-like, which is an important lipid metabolic protein in neural tissue, was drastically down-regulated by 5.5- and 5.6-fold in low and high [BU]-treated larvae.

More surprisingly, four vitellogenin-like unigenes were drastically down-regulated in both low and high [BU] treated larvae. The magnitude of down-regulation ranges from 5.6- to 7-fold. These four B. neritina vitellogenin-like unigenes were two short sequences (Cluster-12452.14875, Cluster-12452.15808) that were shorter than 500 bp and two long sequences (Cluster-12452.20127, Cluster-12452.15767) that were longer than 4,500 bp. Nevertheless, the deduced ORFs of all four unigenes contained the lipoprotein domain (PF01347) (or Vitellogenin N-terminal domain).