RNA-Seq reveals divergent gene expression between larvae with contrasting trophic modes in the poecilogonous polychaete Boccardia wellingtonensis

The polychaete Boccardia wellingtonensis is a poecilogonous species that produces different larval types. Females may lay Type I capsules, in which only planktotrophic larvae are present, or Type III capsules that contain planktotrophic and adelphophagic larvae as well as nurse eggs. While planktotrophic larvae do not feed during encapsulation, adelphophagic larvae develop by feeding on nurse eggs and on other larvae inside the capsules and hatch at the juvenile stage. Previous works have not found differences in the morphology between the two larval types; thus, the factors explaining contrasting feeding abilities in larvae of this species are still unknown. In this paper, we use a transcriptomic approach to study the cellular and genetic mechanisms underlying the different larval trophic modes of B. wellingtonensis. By using approximately 624 million high-quality reads, we assemble the de novo transcriptome with 133,314 contigs, coding 32,390 putative proteins. We identify 5221 genes that are up-regulated in larval stages compared to their expression in adult individuals. The genetic expression profile differed between larval trophic modes, with genes involved in lipid metabolism and chaetogenesis over expressed in planktotrophic larvae. In contrast, up-regulated genes in adelphophagic larvae were associated with DNA replication and mRNA synthesis.


Álvaro Figueroa 1* , Antonio Brante 2,3 & Leyla Cárdenas 1,4
The polychaete Boccardia wellingtonensis is a poecilogonous species that produces different larval types. Females may lay Type I capsules, in which only planktotrophic larvae are present, or Type III capsules that contain planktotrophic and adelphophagic larvae as well as nurse eggs. While planktotrophic larvae do not feed during encapsulation, adelphophagic larvae develop by feeding on nurse eggs and on other larvae inside the capsules and hatch at the juvenile stage. Previous works have not found differences in the morphology between the two larval types; thus, the factors explaining contrasting feeding abilities in larvae of this species are still unknown. In this paper, we use a transcriptomic approach to study the cellular and genetic mechanisms underlying the different larval trophic modes of B. wellingtonensis. By using approximately 624 million high-quality reads, we assemble the de novo transcriptome with 133,314 contigs, coding 32,390 putative proteins. We identify 5221 genes that are up-regulated in larval stages compared to their expression in adult individuals. The genetic expression profile differed between larval trophic modes, with genes involved in lipid metabolism and chaetogenesis over expressed in planktotrophic larvae. In contrast, up-regulated genes in adelphophagic larvae were associated with DNA replication and mRNA synthesis.
Marine invertebrates exhibit contrasting developmental modes that may affect the speciation, extinction, and connectivity of species 1 . In species that encapsulate their offspring, the indirect developmental mode is characterized by embryos that develop partially inside capsules and hatch as planktotrophic larvae. This type of larvae shows a set of morphological traits that allows them to swim and feed on plankton. At the other end of the spectrum, species with direct development complete the whole larval development inside the capsule and hatch at the juvenile stage [1][2][3] . Within the direct development mode, larvae may show different intracapsular feeding strategies: lecitotrophy, in which larvae use an internal energetic source provided in the egg by way of yolk (lecitotrophic larvae); ovophagy, in which larvae feed on nurse eggs provided by the mother (ovophagic larvae); and adelphophagy, in which larvae feed on other embryos (adelphophagic larvae) 2 . From an ecological perspective, lecithotrophy and adelphophagy appear to be advantageous strategies due to the larger resulting hatching sizes, higher survival rates, and better competitive abilities of the hatchers [3][4][5] . On the other hand, planktotrophy allows a potentially greater larval dispersion range and the colonization of new habitats, thus reducing resource competition 6 .
The different nutritional strategies of the larval trophic modes have been associated with morphological evolutionary transitions, and in most marine invertebrate species they are considered taxonomic traits 7,8 . However, in some taxonomic groups with larval encapsulation, such as polychaetes and gastropods, the production of different larval types in the same species has been documented 9 . This reproductive strategy, called poecilogony, is a rare but interesting characteristic that depending on the species, different larval trophic modes can be observed www.nature.com/scientificreports/ within or between capsules, clutches, or populations; additionally, females may produce different larval types along the time. For example, in the poecilogonous polychaetes Boccardia wellingtonensis and Boccardia proboscidea, females lay both planktotrophic and adelphophagic larvae in the same clutch or capsule. These larvae are morphologically similar but differ in their feeding capacities; while both larval types can swim and feed on plankton if they are released into the water, only adelphophagic larvae may feed on nurse eggs and other larvae inside the capsules 10,11 . This suggests the existence of underlying mechanisms, other than morphological traits, that would explain the feeding behavior and capacity of marine invertebrate larvae. In this way, marine species that exhibit poecilogony provide an interesting biological model to study the evolutionary transitions between the different larval trophic strategies as well as the developmental modes in marine invertebrates because they enable intraspecific comparisons without underlying phylogenetic effects 12 .
The poecilogonous polychaete Boccardia wellingtonensis inhabits the intertidal zones of Chile, South Africa, and New Zealand. Females of this species reproduce throughout the year laying capsules that resemble the shape of a pearl necklace, which intensifies sibling competition for food at the larval stage. Oyarzun and Brante 10 describe two reproductive strategies in this species: Type I females produce only planktotrophic larvae while Type III females produce capsules containing both planktotrophic and adelphophagic larvae, which are indistinguishable from one another in the early stages of development, as well as nurse eggs. While planktotrophic larvae stop growth at mid developmental stages and do not appear to feed inside their capsules, adelphophagic larvae feed on nurse eggs and other embryos until the offspring hatch at an advanced larval stage or as juveniles 10 (Supplementary File S5). Previous works have found differences in the enzymatic activity of some digestive enzymes, such as hydrolases, hydrolases acids, and glycosidase, in the larvae of B. wellingtonensis that could be associated with the different larval trophic modes observed in this species 13 . Moreover, differences in growth rates of isolated embryos in vitro 10 and the differing expression of microRNAs between both types of larvae 14 suggest that the genetic mechanisms operating during early development could be driven trophic mode evolution in poecilogonous species.
Differences in gene expression profiles have been associated with contrasting early developmental patterns in marine invertebrates. For example, in echinoderms, transcriptomic analyses have expanded the knowledge of the cellular mechanisms responsible for developmental processes, for which a Gene Regulatory Network (GRN) has been described in detail 15,16 . The determination of GRNs has revealed differences in the expression of genes with known roles in skeletal, endomesoderm, and ectoderm control in different species with planktotrophic and lecithotrophic development 17 . Likewise, significant differences have been identified in gene expression in the early developmental stages between planktotrophic Clypeaster subdepressus larvae and facultative planktotrophic C. rosaceus larvae 18 . In the case of the poecilogonous polychaete Pygospio elegans, the gene-expression profiles of the two types of females, planktotrophic and benthic, have a reproductive bias 19 . Between the two types, more than 8000 differentially expressed unigenes have been identified. While benthic females were enriched in the transcription of genes associated with extracellular protein production including collagens, vitrin, and fibropellin-1, planktotrophic females showed enrichment in genes associated with metabolism and intracellular proteins 19 . From an evolutionary perspective, the plasticity in the aforementioned gene expression suggests its potential relationship with the evolutionary transition of different larval trophic modes. In this framework, exploration of these potential underlying mechanisms in B. wellingtonensis would improve the existing understanding of evolution and the diversity of developmental modes in marine invertebrates.
Here, we study the molecular mechanisms that could determine the different larval trophic modes in Boccardia wellingtonensis through the first transcriptome for this species developed using the RNA-seq approach. Transcriptomes were obtained from: (1) males, (2) females showing Type I and Type III reproductive strategies, (3) Type I planktotrophic larvae, (4) Type III planktotrophic larvae, (5) early adelphophagic larvae, and (6) late adelphophagic larvae. Knowing the transcriptome allows us to correctly interpret the functional elements of the genome and uncover molecular factors at the cellular, tissue, and individual levels 20 . From that point of view, in this work we examine whether transcriptional changes could explain the diversity of larval trophic modes in B. wellingtonensis.

Results
De novo transcriptome assembly and quality assessment. We sequenced the RNA libraries of seven sample groups of B. wellingtonensis individuals including three groups of adults: Type I females (FTI), Type III females (FTIII), and males (M); and four groups of larvae: Type I planktotrophic larvae (PLTI), Type III planktotrophic larvae (PLTIII), early adelphophagic larvae (EAL), and late adelphophagic larvae (AL). Each group had three biological replicates. The total of 21 libraries produced a total of 1,157,058,616 PE reads with an average of 55.1 ± 6.2 (mean ± SD) million reads per sample. The raw reads were subjected to a quality trim (Phred = 0.001), and adapters were removed, resulting in a total of 624,918,240 PE post-trim reads with an average of 85.04 bp (Supplementary File S1). The trimmed reads were subjected to a de novo assembly using CLC Genomics, which produced 172,381 transcripts. To reduce the number of potentially spurious transcripts, we filtered based on the abundance and quality of the assembly. To do so, we used the Salmon selective alignment method (-validateMappings) to detect the levels of expression and redundancies. In the end, 133,314 highly expressed transcripts were kept with an N50 = 1532 bp, a % CG = 46.82 and size = 79 Mb. The longest transcript contains a sequence of 28,098 bp, and 18,805 transcripts measure more than 1 Kb. The length and distribution of all of the assembled transcripts are shown in Supplementary File S2. Meanwhile, a BUSCO (Universal Single-Copy Orthologs) analysis revealed a high rate of 862 (88.14%) orthologous genes from the Metazoa_odb9 database with few fragmented or missing BUSCO (Complete gene representation: 76.69% [Single copy: 73.21%, Duplicated: 3.48%], Fragmented: 11.45%, Missing: 11.86%), indicating a high-quality transcriptome. A summary of the assembly of the transcriptome and a quality assessment are shown in Table 1.

Differential expression analysis.
To compare the expression patterns among adults, larval trophic modes and developmental stages in B. wellingtonensis, a principal component analysis (PCA) was performed on the normalized transcripts ( Fig. 2A). The PCA showed that the adult and larval samples grouped separately within principal component 1 (PC1) with a 32.1% of the variance explained by this component (Fig. 2A). It should be noted that one of the AL replicates (AL1) showed a lower number of reads in comparison to the other samples resulting in an outlier sample, which could affect the statistical accuracy of the whole dataset (Supplementary File S1). Consequently, for further analyses, the AL samples were not included.
To identify up-regulated genes, we first compared the whole dataset; i.e., adults (FTI, FTIII, and M) versus larval stages (EAL, PLTI, and PLTIII). We identified 13,421 transcripts that exhibit differential expression (FDR < 0.01; Log2 Fold Change [> 1]) of which 5,221 are up-regulated in larvae while 8,200 are up-regulated in adults (Fig. 2B). The top ten differentially expressed genes (DEGs) in these comparisons are described in Table 2   in biological processes (BP) (FDR < 0.05 and collapsing by REVIGO) (Fig. 3). These terms include the oxoacid metabolic process (GO:0043436) along with enriched processes related to lipid transport (GO:0006869) and the regulation of the lipid metabolic process (GO:0019216). The organonitrogen compound catabolic process (GO:1901565) and the organonitrogen compound metabolic process (GO:1901564) are also included in this group (Fig. 3). Eleven GO terms were found to be related to molecular functions (MF), and four were related to cellular components (CC). Details of the GO terms enriched in BP, MF, and CC are listed in Supplementary File S4. In adelphophagic larvae, 29 GO terms were related to biological processes (BP) (FDR < 0.05 and collapsing by REVIGO) (Fig. 3). Other noteworthy RNA terms include RNA transport (GO:0050658), regulation of transcription from RNA polymerase II promoter (GO:0043436), mRNA processing (GO:0006397), tRNA metabolic process (GO:0006399), ncRNA metabolic process (GO:0034660), and post-transcriptional regulation of gene expression (GO:0010608). Other relevant GO terms were regulation of developmental process (GO:0050793), cell proliferation (GO:0008283), and mitotic cell cycle process (GO:1903047) (Fig. 3). Twenty-one GO terms were found to be related to MF, and eleven were related to CC. Details of the GO terms enriched in BP, MF, and CC are listed in Supplementary File S4.

Discussion
Here we assembled the first de novo transcriptome of the species B. wellingtonensis, a poecilogonous species, to characterize its genes and their expression profiles focused in the DEGs related to the different larval trophic modes. This study provides a new source of genetic information for research in polychaetes and reveals differences in the gene expression profiles among the different larval trophic modes. The results suggest that the regulation of gene expression could be a key mechanism underlying larval trophic strategies in the poecilogonous species. Early evidence suggests that the transcriptome profile could be associated with the developmental mode in poecilogonous species 19,21 . In our work, we use the recent development of sequencing technology to perform RNAseq experiments to test for the profile of differentially expressed genes of different individuals in a whole genome approach, which enables us to search for the genes particularly responsible for the evolution of their developmental modes.

Comparison of DEGs in larvae compared to adults. Comparisons between the adult individuals and
the larval stages ( Fig. 2) identified DEGs in adults mainly related to physiological processes like reproduction, environmental responses, and the maintenance of the cellular matrix (Table 2). For example, we found an elevated expression of the vitellogenin 6 (VIT6), a the egg yolk precursor protein and for the transport of lipids; in females, it is a specific indicator of oocyte development 22 . Other up-regulated genes are CLIP1, HMCN1, www.nature.com/scientificreports/ VKT1, and ADGR1, which produce proteins for the extracellular matrix and are relevant to the integrity of the structural tissues 23 . In Caenorhabditis elegans, hemicentin (HMCN1) is an extracellular protein that forms regulatory networks in muscular and gonad tissues and facilitates the mechanosensory neurons' anchorage to the epidermis 23 . In the polychaete Schmidtea mediterránea, hemicentin contributes to the correct proliferation of stem cells 24 . Another interesting finding is the elevated expression of thyroglobulin (THYG), a glycoprotein that assists in the synthesis of thyroid hormones; in invertebrates, it can act as a sensor of iodinated tyrosine, an indicator of food availability 25 .
In the larval state, we observed the overexpression of JAG1A, a protein related to the Notch pathway (Table 2), and components of the Notch pathway such as NOTC1 and NOTC3 (Table 2). Recently, it was reported that the Notch pathway in the annelid Platynereis dumerilii are expressed in the structures that give rise to the nascent Table 2. The top ten up-and down-regulated genes of differential expression analysis. *FDR = 0.00 indicates values of 1E−16 or lower. www.nature.com/scientificreports/ chaetae and that deterioration causes malformations and abnormalities in chaetae-producing cells 26 . Therefore, the components of the Notch pathway in larvae of B. wellingtonensis could be related to chaetogenesis and/or to morphological development. Moreover, previous studies have found that the Notch pathway is involved in the regulation of settlement and metamorphosis of polychaetes 27 . In other annelids, such as Capitella teleta and Hydroides elegans, the segmentation process is controlled by Notch, the Wnt signaling pathway, and various Hes genes including Hes1 28  www.nature.com/scientificreports/ File S4). These genes likely regulate the segmented growth in B. wellingtonensis 28 . However, the cumulative evidence does not point to a common mechanism that regulates the segmentation in all annelids 26 . Therefore, it seems necessary to better understand the molecular control of larval development during their intracapsular development and to analyze the temporality and the special expression of genes such as Hes1 and Wnt1 and the components of the Notch signaling pathway.

Comparison of DEGs between adults. Low numbers of DEGs were detected in the comparison between
females (FTI and FTIII). Knowledge of these genes' functions is still poor, and it is therefore not easy to develop a clear picture of the relationships between the female developmental mode and egg production. In FTIII, high levels of expression in SLC3, GGT1, and ITB7 genes were detected. In the polychaete Osedax japonicus, a boneeating worm, the transporter SLC plays a crucial role in nutrient absorption 29 . In the invertebrate chordate Ciona intestinalis, the Gamma-glutamyltransferase 1 (GGT1) is part of the core pathway of genes contributing to glutathione biosynthesis, an important metabolic mechanism in the specie's environmental stress response 30 . In the sea cucumber Apostichopus japonicus, the integrin beta 7 (ITB7) is part of the family of cell adhesion molecules that plays an important role in a wide range of physiological processes including the immune response 31 .

Comparison of DEGs between larvae.
The comparative analyses between planktotrophic (PLTI, PLTIII) and adelphophagic (EAL) larvae suggest that the enrichment of DEGs observed in planktotrophic larvae is related to the regulation of the lipid metabolic process (GO:0019216) and lipid transport (GO:0006869) (Fig. 3). Lipids are the main reserve of metabolic energy in animals, and they play a central role in the larval development of marine invertebrates 32 . Of the total upregulated genes found in PLTI, four genes contain desaturase domains in PLTI: Protocadherin FAT2, Protocadherin FAT4, FAT atypical cadherin 4, and Fat acyl-CoA hydrolase medium chain (Fragment) ( Table 2, Supplementary File S4). Desaturase enzymes are directly involved in the biosynthesis of polyunsaturated fatty acids (PUFA) in invertebrate marine animals 33 . The desaturases are widely distributed across animal phyla, and in fact, multiple invertebrates can actively produce PUFA de novo 34 . For example, the annelid Riftia pachyptila has been shown to have a desaturase that enables it to biosynthesize α-Linolenic acid 35 . In the specific case of, the FAT2 desaturase catalyzes the first step for the biosynthesis of PUFA 36 . Our data suggests that in planktotrophic larvae the desaturase enzymes can selectively control the PUFA metabolism and thus regulate lipid composition. Fourteen other strongly over-regulated genes are lipases and lipoxygenase domain-containing proteins. Interestingly, some genes that PLTI and PLTIII have in common, such as PLB1, Lipasa-GDSL, LIPP, and LOXH1, possess lipase domains that may, in combination with various transport lipoproteins, control the composition of structural lipids such as phospholipids and cholesterol. Based on our results, we hypothesize that in B. wellingtonensis, planktotrophic larvae should maximize the use of lipids as nutrition sources provided by the mother within the egg. It has been observed that in echinoderms with planktotrophic development, the utilization of lipids as a source of energy is dynamic and highly regulated in order to sustain the development of the egg, the growth, and the survival of the larvae to ensure successful settlement and subsequent metamorphosis 37 . Specifically, larvae use lipovitellin, a component found in the internal reserves of embryos and larvae that contains high levels of phospholipids in polychaetes 38 . We observed more DEGs between PLTI and EAL than we did between PLTIII and EAL (Supplementary File S3A,B). In contrast to the congeneric species B. proboscidea, the Type III reproductive strategy in B. wellingtonensis is characterized by a wider range of larval sizes (PLTIII), which is probably influenced by intracapsular food availability 10 . This causes a high variance in the larval size, making it likely to find larvae in different developmental stages in the same capsule. We also hypothesized that maternal imprinting could be another important factor. To arrive at a final conclusion, more experiments should be done to determine which external factors influence the transcript expression during the species' development. It is important to note that over 80% of the up-regulated DEGs in PLTIII are shared with PLTI (Supplementary File S3C). Further investigation of this overlapped expression pattern of 378 genes up-regulated in planktotrophic larval in other Polydora-complex species could shed light on the molecular and genetic mechanisms explaining the diversity of larval traits and the evolution of developmental modes in marine invertebrates 39 .
In PLTI larvae, five up-regulated transcripts can be classified into the Cytochrome P450 family, with CYP 450 3A24 being one of the top DEGs ( Table 2). The CYP 450 family plays a pivotal role in the physiology of marine invertebrates, particularly in annelids, by catalyzing the biosynthesis of signal molecules including steroids 40 . Further, in the embryos of mussels, these enzymes are important for the metabolism of xenobiotics such as carcinogens, pesticides, and drugs 41 . In particular, the expression of CYP 450 3A24 is associated with the palmitoleic acid-induced stress response in the Argopecten irradians 42 . Moreover, other multifunctional enzymes that detoxify the xenobiotics, such as the glutathione S-transferases (GST), are also expressed in PLTI's DEGs, (Supplementary File S4). In this context, the up regulation of CYP450 and GST genes could be fundamental in PLTI's response to environmental stress.
Another remarkable find is the presence of TRPM8 and several acid-sensing ion channels (ASICs) in PLTI ( Table 2 and Supplementary File S4). TRPM8 is an important sensor in cold environments that requires the presence of phospholipids in the membrane 43 . In larvae of marine invertebrates, temperature regulates metabolism, growth rate, and the duration of the larval development 44 . ASICs belong to the family of degenerin/epithelial sodium channels, which are presumably involved in mechanosensation 45 . This mechanosensory response is key during feeding, tactic responses or predator-prey interactions in marine animals 46 . In general, our results suggest that genes with the highest expressions in planktotrophic larvae are related firstly to energy production from endogenous sources, primarily phospholipids, and secondly to the formation of structures and the sensory system that enable the capture of exogenous food to complete the plankton's development. www.nature.com/scientificreports/ Our analyses of EAL compared to planktotrophic larvae (PLTI, PLTIII) found significant increases in the number of GO terms related to: regulation of the developmental process (GO:0050793), cell proliferation (GO:0008283), cellular component organization or biogenesis (GO:0,071,840), cellular process (GO:0009987), and the mitotic cell cycle process (GO:0007267). Furthermore, RNA polymerase II promoter (GO:0006357), mRNA processing (GO:0006397), and mRNA metabolic process (GO:0016071) were also detected (Fig. 3). These results align with an adelphophagic development that, in comparison with planktotrophic development, is characterized by a rapid increase in size until hatching at the juvenile stage. Oyarzun and Brante 10 explained that shortly after the first segment appears in adelphophagic larvae of B. wellingtonensis, larvae begin to feed on nurse eggs reaching an average length of 13.3 setigers after 15 days at 20 °C 10 . Studies in the nemertine Lineus ruber provided evidence that cells in a proliferative stage are widespread throughout fast adelphophagic development, from gastrulation until metamorphosis 47 . Cell proliferation is highly regulated by transcription factors that contain bHLH (basic helix-loop-helix) domains 48 . We observed a high expression of 37 transcription factors in the EAL of B. wellingtonensis; of them, five DEGs contained the bHLH domain (Supplementary File S4).
Meanwhile, another important group of 39 DEGs were related to DNA replication (Supplementary File S4). Nevertheless, despite the fast growth of EAL, Gibson et al. 49 specifically described that early adelphophagic larvae present a delay in the differentiation of the endoderm-and mesoderm-derived tissues 50 . Therefore, at this point in their development, the heterochrony observed in adelphophagic larvae may be related to the high expression of the thyroid hormone receptor β (THB) gene ( Table 2). In echinoderm species, it has been demonstrated that the thyroid hormone (TH) accelerates development to metamorphosis 51 . Typically, larvae obtain TH through consumption of microalgae; however, the complete intracapsular development of EAL indicates that these larvae may produce TH by themselves or could obtain it from nurse eggs or even via cannibalism 10 .
Lastly, another group of DEGs that could be particularly important in EAL of B. wellingtonensis are HLH6, t-bx2, and t-bx20 (Supplementary File S4). HLH6 is a transcription factor in C. elegans that acts as a development regulator in the cells of the pharyngeal gland and is necessary for feeding 52 . The deletion of HLH6 in this species produces a phenotype similar to that observed in starvation conditions (slow growth rate and small size) 52 . The t-bx2 gene regulates the development of the pharyngeal muscles 53 , and t-bx20 plays a role in the development of the hindgut 54 . Chowdhuri et al. observed arrested development in tbx-2 mutant larvae of C. elegans associated with an affected pharynx that made ingestion impossible 53 . In annelids, t-bx2 is expressed in the dorsal side of the embryo 55 , and in more advanced stages of development, the expression of t-bx2 is clearly associated with the intestine and the pharynx 56 . This suggests a conserved expression of the HLH6 and tbx genes during the development of the intestine and the pharynx in invertebrates. Additional research on these genes will be required to understand their evolution and function in polychaetes' larval feeding patterns.
In conclusion, our results indicate that the different larval trophic modes of B. wellingtonensis present transcriptomic plasticity related to changes in the expression of specific genes that could be related to the evolution of larval developmental modes in marine invertebrates. These genes could be tested in future studies by modifying their nutrition sources and lab-controlled environmental factors, and results could be refined as more functional data are accumulated. Comparative RNAseq analyses to other species displaying similar reproductive modes could be a productive avenue in distinguishing whether these differences are a product of adaptive evolution specific of B. wellingtonensis or are due to the evolutionary history of poecilogonous species.

Methods
Sampling locations and sample processing. Samples for transcriptome sequencing were collected from two localities in Chile: Coliumo (36° 33 11 S; 72° 57 24 W) and Los Molinos (39° 50 46 S; 73° 23 54 W). Samples were found in the high intertidal zone perforating sedimentary rock. They were transported alive to the laboratory and processed within 24 h. Individuals and larvae were classified under a binocular stereo microscope. For subsequent analyses, the samples were classified into seven groups: Type I planktotrophic larvae (PLTI), Type III planktotrophic larvae (PLTIII), early adelphophagic larvae (EAL), late (juvenile) adelphophagic larvae (AL), Type I Female (FTI), Type III Female (FTIII), and male (M). Each group contains three biological replicates. Determination of the developmental stage of the larvae was based on the number of setigers (segment bearing chaetae), a method commonly used in polychaetes 13 . Larvae of the groups EAL, PLTI, and PLTIII used for the experiment had 5 setigers while the AL had at least 10 setigers. Adelphophagic larvae with less than 5 setigers cannot be differentiated from PLTIII. Supplementary File S5 depicts a summary of the larval development and reproductive strategies of B. wellingtonensis. In order to have sufficient tissue for the analyses, larvae were collected from capsules of different mothers and pooled depending on their trophic type and developmental stage.
RNA extraction. Samples were homogenized in QIAzol Lysis Reagent (Qiagen). Total RNA extraction was performed using the Universal RNeasy Mini Kit with a gDNA Eliminator Solution (Qiagen) according to the manufacturer's instructions. We used a NanoDrop 2000 Spectrophotometer (NanoDrop Technologies) and Fragment Analyzer to assess quantity, purity, and integrity of the RNA samples. All samples used for subsequent analyses had to have RIN > 8, as previously described 14,57 . Illumina sequencing and data quality control. Starting  Assembly and annotation. The sequenced raw data was processed using the CLC Genomics Workbench software v.12.0.3 (Qiagen Bioinformatic) in a four-step pipeline: (i) removal of adapter sequences, (ii) removal of ambiguous nucleotides ('N'), (iii) trimming base pairs with a Phred quality score 0.001 from the 3′-end of each sequence, and (iv) removal of reads shorter than 50 bp after trimming. The trimmed reads were subjected to a quality assessment using FastQC and MultiQC. High quality reads were assembled using the CLC de novo assembly option, and the following parameters were applied: mapping mode = Map read back to contig (slow), bubble size = 93, minimum contig length = 200 nt, word size = 25, perform scaffolding = yes, auto-detect paired distances = yes. The quality and completeness of the new draft assemblies were analyzed using the software QUAST for assembly statistics 58 . Moreover, comparisons were made to the Benchmarking Universal Single-Copy Orthologs (BUSCO) database using the Metazoan (metazoa_odb9), which consists of 978 single-copy genes that are present in at least 90% of Metazoan 59 . TransDecoder v5.3.0 was used to identify putative protein coding regions, and homology options were included as retention criteria for the candidate ORFs. Predicted ORFs were used for BLAST searches and annotation against the UniProt (Swiss-Prot) and Uniref90 using the BLASTX algorithm with an e-value cutoff of 1e−10. Annotated unigenes were further searched for Gene Ontology (GO) terms using the Blast2GO software in the main GO categories of molecular functions, biological processes, and cellular components 60 . Complementary annotations were done with the InterProScan v.5 software, which provides functional analysis of proteins by classifying them into families and predicting domains and important sites. The annotation results were further fine-tuned with the EggNOG and GO slim functions of the Blast2GO software in order to improve and summarize the functional information of the transcriptome dataset. Supplementary File S6 depicts the bioinformatics workflow used in this study.
Differential gene expression analysis. Trimmed reads of each sample were mapped back separately on the high-quality set transcripts by applying the following parameters (mismatch cost = 2, length fraction = 0.8, similarity fraction = 0.8 and maximum number of hits for a read = 10) and using the RNA-Seq tool of the CLC Genomics Workbench v12.0.3 (Qiagen Bioinformatic). Thus, the expression of each transcript in each sample was calculated as Transcript Per Million (TPM). Subsequently, a DEA between adults, larval samples, and intro larvae was carried out in order to identify DEGs using the Differential Expression for RNA-Seq tool of the CLC Genomics Workbench v12.0.3 (Qiagen). The normalization factors were calculated using the trimmed mean of M-values (TMM) method 61 . The threshold FDR < 0.01 was adjusted to identify the differentially expressed genes by log2 fold change [≥ 1].
Gene ontology and enrichment analysis. A GO enrichment analysis was conducted to reveal GO terms significantly enriched (p value < 0.01) in DEGs using the CLC Genomics Workbench v12.0.3 (Qiagen). The redundant enriched GO terms were removed using the web server REVIGO (http:// revigo. irb. hr) 62 The SimRel semantic similarity was used with small-allowed similarity (0.5), and all other parameters were kept as default. Resulting .cvs tables were exported and used to graphically plot the condensed GO terms using ggplot2 in Rstudio 63 .