Genomic insights into the sessile life and biofouling of barnacles (Crustacea: Cirripedia).

Members of the infraclass Cirripedia, commonly called barnacles, are unique among the subphylum Crustacea in that they exhibit a biphasic life cycle with a planktonic larval stage and a sessile adult stage. Understanding their unique sessile life and mechanisms of attachment are hampered by the lack of genomic resources. Here, we present a 746 Mb genome assembly of Lepas anserifera - the first sequenced stalked barnacle genome. We estimate that Cirripedia first arose ~495 million years ago (MYA) and further diversified since Mesozoic. A demographic analysis revealed remarkable population changes of the barnacle in relation to sea-level fluctuations in the last 2 MYA. Comparative genomic analyses revealed the expansion of a number of developmental related genes families in barnacle genomes, such as Br-C, PCP20 and Lola, which are potentially important for the evolution of metamorphosis, cuticle development and central nervous system. Phylogenetic analysis and tissue expression profiling showed the possible roles of gene duplication, functional diversification and co-option in shaping the genomic evolution of barnacles. Overall, our study provides not only a valuable draft genome for comparative genomic analysis of crustacean evolution, but also facilitates studies of biofouling control.

Introduction

Barnacles (Thecostraca: Cirripedia) are sessile crustaceans ubiquitously present in marine environments from the intertidal to the deep-sea, and from the tropics to the polar regions; some of them even live on other marine organisms such as whales and turtles (Chan et al., 2009; Buckeridge et al., 2018, 2019; Kim et al., 2020). The subclass Cirripedia is divided into three infraclasses that differ in morphology and way of life: Thoracica includes acorn and stalked barnacles with phosphatic or calcified shells (Chan and Høeg, 2015; Chan et al., 2021); Acrothoracica is characterized for not having shells but specifically burrowing into diverse calcareous substrates (Chan et al., 2014; Nielsen et al., 2016); Rhizocephala includes endoparasitic barnacles of other crustaceans (Anderson, 1994; Jung et al., 2019). Barnacles represent an excellent model system to study the evolution of body forms owing to their drastic changes from free-swimming planktonic naupliar and cyprid larval stages to sessile juveniles and adults with their forehead fixed to a solid substrate and their body placed upside down during the life cycle (Qiu et al., 1997; Qiu and Qian, 1999). Although many other crustaceans also possess the naupliar stage, other stages of the barnacle life cycle are unique. Unlike the body plan of many crustaceans with six thoracic segments and five abdominal segments, barnacles lack completely developed abdominal segments in the cyprid and sessile stages (Anderson, 1994). Among crustaceans, only barnacles have a cement gland that secrets a protein matrix for underwater adhesion, which is critical for permanent attachment (Liang et al., 2019).

The strong attachment of barnacles on underwater structures makes them an important group of biofouling organisms. Biofouling incurs economic loss by increase of fuel consumption on sailing vessels and damage of man-made submarine structures (Holm, 2012). Various antifouling coatings have been developed to control barnacle attachment, but many of them contain chemical ingredients that have serious toxic effects on non-target marine life (Qiu et al., 2005, 2008). For example, tributyltin, which had been widely used in antifouling paints before its ban by the International Convention on the Control of Harmful Anti-fouling Systems on Ships, is an endocrine disruptive chemical causing abnormal shell formation in oysters and imposex in neogastropods (Chan et al., 2008a; Qiu et al., 2011; Yi et al., 2012; Ip et al., 2016).

There has been a substantial amount of research on the synthesis and properties of the barnacle cement, especially the roles of barnacle cement proteins (BCPs) in forming the larval adhesive plaque on solid substrates as well as its expansion in the adult stage (Zhong et al., 2014; Liang et al., 2019). Previous studies of BCPs have mainly relied on conventional biochemical methods including identifying their corresponding genes and determining their tissue-specific expressions in different life stages. Assembly of barnacle genomes will not only shed light on the molecular mechanisms of adhesion and evolution of this group of highly derived crustaceans, but also facilitate the investigation of underwater adhesives materials, as well as the development of environmental-friendly antifouling agents (Liang et al., 2019).

In this study, we aimed to sequence and assemble the genome of Lepas anserifera (Cirripedia, Lepadidae; Figure S1), a stalked fouling barnacle widely distributed in temperate and tropical seas, which grows in clusters attached to driftwood and other floating objects by its peduncle (Isaac and Moyse, 1990). Despite their important roles in fouling communities, no species from the family Lepadidae has a sequenced genome. Within the whole infraclass Cirripedia, only the acorn barnacle Amphibalanus amphitrite (Balanidae) has a sequenced genome (Kim et al., 2019). We conducted comparative genomic analyses with other pancrustacean lineages to understand barnacle genomic structures that are likely associated with their unique sessile life with particular emphasis on gene family expansion and contraction. We also produced tissue-specific transcriptomes and analyzed them in order to gain insight into the functions of candidate genes that are likely critical for adaptation in the sessile lifestyle including underwater adhesion.

Materials and methods

Sample collection, library construction, and sequencing

The stalked barnacles were collected by hand from floating buoys in Keelung City, Taiwan (25°08′31.0″N 121°48′08.7″E) in November 2018. They were transported to and cultured with their attached substrate for one week in the Biodiversity Research Center, Academia Sinica. An individual was selected for genome sequencing. To reduce the contamination by bacteria or algae, the shells of the barnacle were cleaned using a brush and rinsed with Milli-Q water before dissection, and the deshelled somatic tissues were rinsed using Milli-Q water. A part of the peduncle (stalk) was cut and fixed in pure ethanol for DNA extraction and genome sequencing. The peduncle, cirri and other tissues were immediately fixed in RNAlater (Invitrogen, CA, USA) for transcriptome sequencing. All tissue samples were stored in a freezer at -80 °C until use for DNA and RNA extraction.

Genomic DNA was extracted using the CTAB method (Porebski et al., 1997). DNA quality was evaluated and its quantity measured using agarose gel electrophoresis and a Qubit fluorometer (Thermo Fisher Scientific, MA, USA), respectively. High-quality DNA (OD 260/280 = 1.8–2.0, OD 260/230 = 2.0–2.2) was sent to Novogene Bioinformatics Technology (Beijing, China) for library preparation and whole genome sequencing using both the Illumina and PacBio platforms (Table S1). In brief, 1 μg DNA was used to construct an Illumina library with a 350-bp insert size, and sequenced on an Illumina NovaSeq sequencer to generate ~120 Gb paired-end reads with a read length of 150 bp. In addition, 10 μg DNA was used to construct a 20-kb SMRTbell library, and sequenced on a PacBio Sequel sequencer to generate ~18 Gb data. For RNA sequencing, total RNA was extracted from the barnacle tissues using TRIzol reagent (Thermo Fisher Scientific, MA, USA). The quality of the samples was determined with 1% agarose gel electrophoresis and the quantity of RNA was determined using a NanoDrop 2000c Spectrophotometer (ThermoFisher, MA, USA). The RNA samples were sent to Novogene Bioinformatics Technology (Beijing, China) for transcriptome library preparation using an Illumina NEBNext Ultra RNA Library Prep Kit (New England Biolabs, MA, USA) and sequenced using an Illumina NovaSeq sequencer to produce 150-bp paired-end reads (Table S1).

Genome assembly

Illumine short reads were trimmed to remove adaptors and low-quality reads (quality score <20, length <40 bp) using Trimmomatic v0.38 (Bolger et al., 2014). A kmer histogram was then generated with the clean Illumina data using Platanus v1.2.4 (Kajitani et al., 2019; settings: “-k 17 -s 10 -u 0.2 -t 29”) to predict the genome size. GenomeScope (Vurture et al., 2017) was used to characterize the heterozygosity and contents of repetitive elements (Figure S2). FM-index Long Read Corrector (FMLRC) was used to correct the PacBio raw reads using Illumina reads (Wang et al., 2018; Wong et al., 2018). After trimming and error correction, a total of 788 million Illumina and 1.05 million PacBio reads (longest length = 120 kb, N50 = 18.2 kb) were retained (Table S1). Two bioinformatic approaches, including one based on the long-reads (PacBio data; wtdbg2 v2.5, Ruan and Li, 2020) and the other (hybrid one) based on both Illumina and PacBio data (MaSuRCA v3.4.1, Zimin et al., 2013) were applied to assemble the genome. A comparison of the assembly statistics showed that the wtdbg2 pipeline resulted in a better initial assembly, in term of genome size, length statistics and the percentage of complete Benchmarking Universal Single-Copy Orthologs (BUSCOs) (Table S2). To improve the initial assembly, the draft genome was further polished with PacBio and Illumina reads using the wtpoa-cns script implemented in wtdbg2. Due to the high heterozygosity rate of Lepas genome, Redundans (Pryszcz and Gabaldón, 2016) was performed to eliminate the haplotypic contigs. The curated contigs were scaffolded with PacBio data using SSPACE-LongRead v1.1 (Boetzer and Pirovano, 2014) and RNA-Seq data using P_RNA_scaffolder (Zhu et al., 2018), and then polished with wtpoa-cns. Possible contamination sequences (e.g., microbe, fungi, algae and virus) were removed by search using BLASTn against the NCBI database with an E-value threshold of 1e-20 and manual correction. The completeness of the final assembled sequences was assessed by analyzing the BUSCO scores (Simão et al., 2015) against the metazoan (metazoa_odb9) and arthropod (arthropoda_odb9) datasets.

Genome annotation and gene prediction

The barnacle genome was annotated using MAKER v3.0 (Cantarel et al., 2008). Before gene prediction, repetitive contents in the genome were identified using RepeatMasker (http://www.repeatmasker.org/) with RepBase library v20150807 (Bao et al., 2015), and species-specific repeat libraries in RepeatModeler v1.0.8 (Smit and Hubley, 2008-2015). Divergences of transposable elements from the consensus sequences and insertion time were determined with the script parseRM.pl implemented in RepeatMasker with a mutation rate of 3.1 × 109 per site per year (Nunez et al., 2020). To provide transcriptomic evidence, de novo and genome-guided transcriptomes of the stalked barnacle were assembled using Trinity v2.5.1 (Haas et al., 2013) based on multiple tissues RNAseq data, and then merged with PASA pipeline v2.2.0 (Haas et al., 2003). Crustacean protein sequences from Swiss-Prot database and protein sequences for other crustacean genomes were used as protein evidence. Augustus 3.1 (Stanke and Morgenstern, 2005) and SNAP (Korf, 2004) were used for ab initio gene prediction. Results from different gene predictors were integrated into a consensus weighted annotation with EVidenceModeler v1.1.1 (EVM; Haas et al., 2008). The predicted genes were functionally annotated using Diamond v0.9.24 BLASTp (Buchfink et al., 2015) against NCBI non-redundant (nr) and UniProtKB databases under the “more-sensitive” option and an E-value threshold of 1e-5. Gene functional annotation was conducted using eggNOG-mapper v2 (Huerta-Cepas et al., 2017) for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, and InterProScan5 (Jones et al., 2014) for InterPro and Pfam.

Heterozygous SNP detection and estimation of population history

The pairwise sequentially Markovian coalescent (PSMC) analysis was used to estimate the history of effective population sizes (Li and Durbin, 2011). The Illumina sequences of Lepas anserifera and Amphibalanus amphitrite (SRR9595623) were aligned back to the respective reference genome using BWA-MEM v0.7.12-r1039 under default settings. Heterozygous variants were called with the pipeline of SAMtools mpileup v1.7 (-C 50) with bcftools call v1.2 and vcfutils.pl vcf2fq (-d 8, -D 100). The PSMC analysis was run under the following settings: -t 15 -r 5 -p 4+25∗2+4+6 and 100 bootstrap replicates assessed variance of the simulated results. Effective population size was calculated using a mutation rate of 3.1 × 10−9 substitutions per site per year (adopted from the acorn barnacle Semibalanus balanoides; Nunez et al., 2020) and a generation time of 1 year.

Phylogenomic and gene family analyses

Orthologous groups among nine arthropod genomes [i.e., Amphibalanus amphitrite (GCA_009805615.1), Daphnia pulex (GCA_900092285.2), Drosophila melanogaster (GCA_000001215.4), Lepas anserifera (this study), Parhyale hawaiensis (GCA_001587735.2), Penaeus vannamei (GCF_003789085.1), Portunus trituberculatus (GCA_008373055.1), Strigamia maritima (GCA_000239455.1), Tigriopus californicus (GCA_007210705.1)], and two barnacle transcriptomes [i.e., Semibalanus balanoides (GHBL00000000.1) and Pollicipes pollicipes (GGJN00000000.1)] and species tree were identified using OrthoFinder v2.2.7 (Emms and Kelly, 2015) with Diamond v0.9.24 BLASTp under “more-sensitive” mode, >50% identity and an E-value threshold of 1e-10 selected. All single-copy genes were aligned using MUSCLE v3.8.31 (Edgar, 2004) and trimmed using TrimAL v1.3 (Capella-Gutiérrez et al., 2009). The aligned sequences were concatenated and used for phylogenetic analyses. A maximum-likelihood method (ML) phylogenetic analysis was conducted using IQ-TREE v1.6.9 (Nguyen et al., 2014), with the best model of LG+F+G4 selected by ModelFinder (Kalyaanamoorthy et al., 2017) and 1000 bootstrapping replicates. A maximum parsimony (MP) analysis was conducted using PAUP4 (Swofford, 2001) with 1000 bootstrapping replicates. A Bayesian (BI) analysis was conducted using MrBayes v3.2.7 (Huelsenbeck and Ronquist, 2001) under the GTR+I+G model; four Markov chains were applied for 10 million generations, sampled every 1000 generations, and the first 25% discarded as the burn-in, respectively. MCMCtree implemented in PAML v4.8 (Yang, 2007) was used to estimate divergence times using the burn-in, sample frequency and number of samples as 10 million, 10,00 and 10,000, respectively. The node calibration among arthropods was based on the TIMETREE database (Hedges et al., 2006), i.e., 560–642 million year ago (MYA) for Arthropoda; 426–561 MYA for Pancrustacea; 467–494 MYA for Hexapoda. Based on the results of orthologous genes and estimated divergence times, gene family expansion and contraction were determined for each node using CAFÉ v4 (Han et al., 2013).

Transcriptomic analysis

Gene expression levels were normalized as Transcripts Per Kilobase Million (TPM) using Salmon v0.9.1 (Patro et al., 2017). One single RNA sample (n = 1) per tissue was prepared. Differentially expressed genes were determined using run_DE_analysis.pl script in Trinity with the edgeR method which utilizes statistical models of expected variation, i.e., Poisson and negative binomial distribution (Haas et al., 2013). Tissue-specific genes for the tissues were determined based on comparing expression levels across tissue types. Only genes with FDR <0.05 and TPM >5 were considered as overexpressed. Potential functions of these target genes were further assessed by GO enrichment analysis using GOseq v1.39.0 (Young et al., 2012).

Real-time PCR was employed to validate tissue expression patterns in twelve selected genes. Total RNA was extracted using TRIzol reagent (Thermo Fisher Scientific, MA, USA) and the first strand cDNA was synthesized by using High Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Massachusetts, USA). The primers were designed using the online NCBI Primer-BLAST tool (Table S11). Real-Time PCR was performed with KAPA SYBR® FAST qPCR Master Mix (2X) Kit (Kapa Biosystems, Massachusetts, USA) on StepOnePlus™ Real-Time PCR System (Applied Biosystems, Massachusetts, USA) through polymerase activation at 95 °C for 3 min, and annealing and extending at 60 °C for 20s with a total of 40 cycles. The specificity of primer pairs for the PCR amplification was checked by the melting curve method. Triplicates were applied for each gene, and the relative gene expression level was calculated based on the 2ΔΔCt method (Livak and Schmittgen, 2001). The cytochrome b gene was used as the internal standard gene as done in a previous study of barnacle gene expression (De Gregoris et al., 2009).

Results and discussion

General genome features

We generated approximately 118.2 Gb Illumina reads and 17.7 Gb PacBio reads (Table S1). Kmer analysis predicted the size of the L. anserifera genome to be 722 Mb (Figure S2A), which is relatively small among members of Thoracica (Cirripedia) whose genome sizes have been estimated (Figure S2B). The Illumina and PacBio sequences generated in this study corresponded to 162 and 24 times of the genome size, respectively. GenomeScope analysis of the Illumina reads predicted the heterozygosity and repetitive contents to be 4.21% and 33.8%, respectively (Figure S2A). Using a long-reads assembly approach coupled with correction using short reads, and elimination of haplotypic contigs with Redundance (Table S2), we produced a final haploid genome assembly of 746 Mb, which is very close to the estimated genome size. The final assembly consists of 2,220 contigs with an N50 of 0.67 Mb and GC content of 46.5% (Table 1). The contig size ranges from 10.1 kb to 3.16 Mb (Figure 1A). BUSCO analysis showed that the L. anserifera genome has 96.1% arthropod BUSCOs (94.6% complete +1.5% fragmented) and 95.0% metazoan BUSCOs (92.4% complete +2.6% fragmented). These L. anserifera BUSCO scores are comparable to those of the acorn barnacle A. amphitrite, and those of other pancrustaceans (Figure 1C; Table S3). The high quality of the L. anserifera assembly was also demonstrated by the high mapping rates of the Illumina DNA reads (95.2%), PacBio reads (98.3%) and de novo transcriptome (94.2%) (Table S4). A total of 28,040 protein-coding genes (PCGs) were predicted in the L. anserifera genome, with an average length of 1,517 bp. BUSCO analysis of the PCGs revealed 94.6% (84.8% complete +9.8% fragmented) arthropod BUSCOs and 94.5% (85.2% complete +9.3% fragmented) metazoan BUSCOs (Figure 1C; Table 1 and S3). Of the PCGs, 20,989 (74.9%) were found to be homologous to proteins in the NR/Uniport databases and 19,044 (67.9%) could be assigned functional domains using InterProScan5 (Table 1).

Table 1

Comparison of assembly and annotation statistics between the stalked barnacle genome assembly and the published acorn barnacle assembly.

	Lepas anserifera	Amphibalanus amphitrite
A. Genome assembly
Estimated genome size (Mb)	722.3	481
Assembly total length (Mb)	756.4	613.4
Repeat content (Mb)	342.1 (45.2%)	274.1 (44.7%)
GC content (%)	46.5	49.8
Unknown Sequences (%)	0	0.46
No. of contig (scaffold)∗	2,220	3,308 (2,240)
N50 (Mb)	0.67	0.31 (0.46)
Average length (Mb)	0.34	0.18 (0.27)
Max. length (Mb)	3.16	1.45 (2.22)
Min. length (kb)	10.1	1.7 (1.9)
No. of contig >10 Kb	2,220	3,168 (2,184)
BUSCO arthropoda#	C:94.6%, F:1.5%, M:3.9%	C:94.1%, F:2.3%, M:3.6%
BUSCO metazoa#	C:92.4%, F:2.6%, M:5.0%	C:93.4%, F:1.1%, M:5.5%
B. Gene prediction
No. of genes	28,040	28,182
Average gene length (bp)	1,517	1,628
Average no. of exons per gene	4.84	5.79
Average exon length (bp)	313.6	280
BUSCO arthropoda#	C:84.8%, F:9.8%, M:5.4%	C:86.6%, F:6.6%, M:6.8%
BUSCO metazoa#	C:85.2%, F:9.3%, M:5.5%	C:85.6%, F:5.2%, M:9.2%
C. Annotation
With hits	22,132	22,867
NCBI Nr	20,705	--
UniProtKB	20,706	--
eggNOG	18,141	--
InterPro	18,990	--
Pfam	16,449	--
GO	10,538	--
KEGG	11,073	--
COG	18,146	--

Lepas anserifera genome only contains contigs, and Amphibalanus amphitrite genome contains both contigs and scaffolds.

C, number of complete BUSCOs; F, number of fragmented BUSCOs; M, number of missing BUSCOs.

Figure 1

Genome characteristics of Lepas anserifera. (A) Length distributions of the assembled contigs. (B) Composition of repetitive sequences (doughnut plots) and transposable element (TE) landscape (histogram) of L. anserifera and Amphibalanus amphitrite. Percentage of divergence from consensus (by bins of 1%) is used as a proxy for age: the aged TEs have accumulated mutations with higher percentage of divergence; young TEs show little divergence from consensus. (C) Arthropoda BUSCO scores (top bar: genome; bottom bar: protein) of eight pancrustaceans.

Repetitive sequences and genome size evolution

Repeat content analysis showed that 341.3 Mb (45.2%) of the L. anserifera genome assembly consists of transposable elements (TEs), whereas the A. amphitrite assembly comprises 272.7 Mb (44.7%) repetitive sequences (Table S5). These sequences, which amounted to 68 Mb more than A. amphitrite, accounted for most of the difference in size (48%) between the two assembled genomes. The two barnacles also differed in their TE composition, with L. anserifera containing more DNA transposons (7.1% vs 5.6% of total genome size) and long interspersed nuclear elements (LINEs, 5.7 vs 3.6%), and A. amphitrite having higher proportions of long terminal repeats (LTRs, 3.6% vs 2.0%) and unclassified TEs (28.0% vs 24.5%) (Figure 1B; Table S5). Since TEs have been shown to be a major contributor to the evolution of genomic size (Kidwell, 2002), we compared the TE divergence and insertion times. Higher levels of divergence were determined in most of the major TE classes (i.e., LINEs and DNA transposons) of L. anserifera when compared with A. amphitrite (Figure 1B). TEs insertions were active in L. anserifera since ~100 million year ago (MYA) (~30% divergence), whereas the burst of TEs in A. amphitrite occurred in ~50 MYA (~15% divergence). These results indicate that the expansion of TEs in acorn and stalked barnacles occurred after the divergence of these two major barnacle lineages (Figure 3C), and such expansion may have contributed to the differences in their genome sizes.

Figure 3

Lepas anserifera gene contents, and phylogenetic relationship, divergence times and Hox genes of pancrustaceans. (A) Venn diagram of the orthologous gene families among L. anserifera, barnacles and other crustaceans. (B) Phylogeny of arthropods inferred from 171 one-to-one orthologous genes (79,437 amino acid positions with 1.0–24.6% gaps among species) using ML, MP and BI methods with support values on nodes, respectively. Black circles on nodes indicate 100% bootstrap support from all three analyses. (C) Estimated divergence time in million year ago (MYA) among 11 pancrustaceans. Numbers on the branch show the lineage-specific expanded/contracted families. (D) Hox clusters among nine assembled genomes. The first letter of the genus name and the first two letters of the species name are used in abbreviations, e.g., Lan for Lepas anserifera.

Genomic variation and population history

We identified 2,976,254 and 2,106,553 heterozygous single-nucleotide polymorphisms (SNPs) in the L. anserifera and A. amphitrite genome, respectively. With these heterozygous SNPs, the demographic history of the two barnacle species from 2 MYA to 10 thousand year ago (KYA) was reconstructed using the Pairwise Sequentially Markovian Coalescent (PSMC) model (Figure 2). The population size of L. anserifera was estimated to be larger than that of A. amphitrite. The population histories of the two barnacle species showed a remarkable relation to the eustatic sea-level fluctuation. A drastic decline in effective population size was apparent ~120 KYA, coinciding well with the inception of the Last Glacial Period (115 KYA – 12 KYA). For coastal marine species (including barnacles), the rapid cooling and expansion of glaciers caused a drastic reduction in the sea level and led to a decrease in available coastal habitat, resulting in fragmented populations, genetic bottleneck and population restructuring (Ludt and Rocha, 2015; Tsang et al., 2012). In addition, since adult L. anserifera usually attaches to floating objects and are transported by ocean currents which may enhance their genetic connectivity (Buckeridge, 2012; Schiffer and Herbig, 2016), the genetic variation and population size of L. anserifera was estimated to be larger than those of the A. amphitrite, which usually inhabits rocky shores.

Figure 2

Demographic history inferred from barnacle genomes. Pairwise sequentially Markovian coalescent (PSMC) analysis was conducted using variant data from the sequenced individual. The blue and pink line represents the population size changes in Lepas anserifera and Amphibalanus amphitrite, respectively. The green and light purple lines, around the blue and pink lines, are the PSMC estimates on 100 rounds of bootstrapped sequences. The orange line denotes the fluctuation of the global sea level relative to the present day (de Boer et al., 2014). The light blue shading represents the Last Glacial Period in the Late Pleistocene. g, generation time; μ, mutation rate.

Phylogenetic relationships and divergence times

Phylogenetic relationships of pancrustaceans were constructed by including the two available barnacle genomes, six other pancrustacean genomes and two published barnacle transcriptomes with high completeness based on BUSCO analysis (Table S3). Strigamia maritima (Myriapoda) was used as the outgroup. Orthologous analysis identified that the L. anserifera assembly contained 10,358 gene families, in which 7,662 are shared with other crustaceans, and 2,669 are unique to L. anserifera (Figure 3A). The phylogenomic trees, reconstructed using Maximum Likelihood (ML), Maximum Parsimony (MP) and Bayesian Inference (BI) methods based on 171 single-copy orthologs with 79,437 amino acid positions, consistently recovered with strong support L. anserifera as an early diverging Cirripedia. It is sister to a clade including A. amphitrite, Pollicipes pollicipes and Semibalanus balanoides. Together, these four species form the Cirripedia clade, which is sister to the Malacostraca clade. Cirripedia and Malacostraca form a clade that is sister to a clade comprising Copepoda and Hexapoda, although this conclusion is not fully supported in ML (bootstrap value 86%) or MP (bootstrap value 85%) analysis (Figure 3B). Furthermore, our analysis recovered Branchiopoda as sister to a clade comprising Cirripedia, Malacostrea, Copepoda and Hexapoda, but again this conclusion is not fully supported by ML or MP analysis.

Despite the use of genome-scale data, previous studies did not fully resolve the phylogenetic relationships of the major groups of pancrustaceans, but suggested the two clades: i) Malacostraca, Copepoda, and Cirripedia; and ii) Branchiopoda, Remipedia, and Hexapoda (von Reumont et al., 2012; Lozano-Fernandez et al., 2019). To evaluate our tree topology, we performed the Approximately Unbiased (AU) test for all 105 possibilities among Branchiopoda, Copepoda, Hexapoda, Malacostraca and Cirripedia with the outgroup of Myriapoda. The AU topology test rejected the sister relationship between Copepoda and (Malacostraca, Cirripedia) at a significance level of 0.05. The statistical results supported the sister group relationship between Malacostraca and Cirripedia, but the phylogenetic positions of Branchiopoda, Copepoda and Hexapoda remain unsolved (Table S6). The divergence time estimation based on all single-copy orthologs showed that the Pancrustacea emerged between the Ediacaran and the Cambrian, and Cirripedia split with Malacostraca at the end of the Cambrian (~495 MYA). Cirripedia became diversified in Mesozoic, and Lepas separated from other groups of barnacles 154 MYA (95% confidence interval of 110–212 MYA; Figure 3C). Our divergence time estimation of Cirripedia (Sessilia—Pedunculata) is different from previous estimation of >500 MYA based on 18S genes (Pérez-Losada et al., 2004), but provides strong statistical support for the result of ~220 MYA in a multiple-gene analysis (Linse et al., 2013). Within the Cirripedia clade, our result showed that Lepas is sister to the clade containing Pollicipes, Amphibalanus and Semibalanus.

Dynamic evolution of Hox genes in crustaceans

The homeobox (Hox) genes are conserved transcription factors that play crucial roles in the development and segmentation of the anterior–posterior axis in metazoans (Sun et al., 2019; Ip et al., 2021). The ancestral arthropod has been hypothesized to possess ten Hox genes (eight with homeotic function and two Hox-like genes) (Grenier et al., 1997). To reveal the Hox cluster in Cirripedia, a BLASTp search of arthropod Hox genes was performed for the two barnacle genomes. The L. anserifera genome assembly encodes nine Hox genes in five contigs including lab, pb, Hox3, Dfd, Scr, Diva, Antp, Ubx and Abd-B (Figure 3D). Among them, Diva, the homologue of arthropod ftz, is located between Scr and Antp (contig669). Neither the L. anserifera assembly nor the A. amphitrite assembly has Abd-A. The missing of this gene was further confirmed by a phylogenetic analysis of the pancrustacean bithorax complex (Figure S3). In Cirripedia, the abdominal segments are reduced (Anderson, 1994). Previous studies of Hox genes in three cirripedes using molecular cloning and in situ hybridization revealed the expression of Abd-B, but not Abd-A – a gene required for the development of anterior abdomen in other non-barnacle crustaceans (Mouchel-Vielh et al., 1998; Blin et al., 2003). Previous studies of Rhizocephalan Sacculina barnacles, which has no distinct abdominal structures in adults but has an extensive root system, revealed no expression of Abd-A related Hox and para-Hox genes (Antp, Ubx and caudal) in the vestigial abdomen during development (Rabet et al., 2001; Mouchel-Vielh et al., 2002; Blin et al., 2003). However, Ascothoracida, the sister-group of Cirripedia possessing a well-developed abdomen, does have an Abd-A. Therefore, Abd-A might have been lost during the evolution of cirripedes, and this may be related to the conspicuous change in the body plan (reduction of the abdomen) as a result of their adoption of sessile life (Deutsch and Mouchel-Vielh, 2003).

In arthropods, ftz has been implicated in the early development (Hughes and Kaufman, 2002), but its function and expression profiles are highly divergent in different crustacean lineages (Gibson, 2000). For instance, in the water flea Daphnia pulex ftz is expressed in the posterior mandibular segment (Papillon and Telford, 2007), while the ftz homolog is missing in the amphipod Parhyale hawaiensis (Serano et al., 2016). In Cirripedia, Diva (the homologue of ftz) has lost its Hox-like expression pattern, and it is expressed solely in the central nervous system (CNS), suggesting its function in neuronal differentiation (Mouchel-Vielh et al., 2002). Unlike other crustaceans, the barnacle cyprid has a well-developed CNS and sense organs. This larva detects settlement cues and coordinates the substrate exploratory behaviors using a pair of antennules equipped with terminal and subterminal sensory setae (Harrison and Sandeman, 1999; Chan et al., 2018). During metamorphosis, the CNS system is largely reconstructed in adult barnacles with a high diversity of sensory setae on the cirri for detection of prey and water currents (Chan et al., 2008b). Therefore, the evolution of Diva and its CNS specific transcriptional profile might have facilitated the development of CNS – a key adaptation of barnacles to the sessile adult life.

Gene family expansion provides insight into environmental adaptation

Reciprocal BLAST analysis of L. anserifera and A. amphitrite gene sets identified >85% genes with homologs and 10,212 one-to-one orthologs in these two species (Table S7). As changes in gene family size are one of the primary drivers of evolution and adaptation (Demuth and Hahn, 2009), we compared gene family sizes between cirripedes and other pancrustaceans. To ensure the reliability of the analysis, the genome of Portunus trituberculatus with relatively low BUSCO values and two de novo transcriptomes were excluded. Results showed that the Cirripedia clade had much more expanded (119) than contracted (4) gene families (Figure 3C). The L. anserifera genome had 111 expanded and 23 contracted gene families. The A. amphitrite had 77 expanded and 47 contracted gene families (Figure 3C).

The expanded gene families unique to each of the two barnacles differ in their enriched functional categories (Table S8). In the stalked barnacle, they mainly include neurogenesis, such as cerebellar granule cell differentiation, synapse assembly and neuromuscular process; whereas in the acorn barnacle the primary functions are protein metabolism, such as proteolysis and protein deubiquitination. These differences likely reflect adaptations of the two barnacles to their way of life: the enhanced proteolytic capability in A. amphitrite may be advantageous for cementation and attachment under strong waves in the intertidal and shallow subtidal areas (Dickinson et al., 2009); whereas in L. anserifera, which are attached to floating objects, enhanced neurogenesis may help them better sense and respond to the dynamic environmental conditions in the water column.

To understand the evolution of Cirripedia, we examined the 119 gene families that are expanded in both the acorn and stalk barnacle genomes. Among them are gene families associated with nervous system development, settlement and biofouling – key traits required for their successful sessile lifestyle (Table S9). One of the most expanded families is broad-complex core protein isoforms 1/2/3/4/5-like (Br–C), which is a transcription factor playing a critical role in ecdysone and juvenile hormone (JH) responsive processes such as molting and metamorphosis in insects (Bayer et al., 1996; Reza et al., 2004). In holometabolans, Br–C is involved in JH signaling, and knockdown of Br–C resulted in failure of metamorphosis from larva to adult (Konopova and Jindra, 2008). Little information of Br–C function is available in crustaceans, except that its regulation of Br–C resulted in the action of methyl farnesoate that functionally resembles JH in crustaceans (Cheong et al., 2015), and it was involved in larval metamorphosis (Sin et al., 2015). A previous phylogenetic study of arthropods discovered that Br–C was specifically gained in the ancestral pancrustaceans and it possibly coevolved with the JH pathway for regulation of metamorphosis (Qu et al., 2015). A total of 80 and 63 Br–C genes were found in L. anserifera and A. amphitrite, respectively. In contrast, only zero to four copies of Br–C were found in each of the other pancrustaceans. Phylogenetic analysis revealed only one origin of Br–C in pancrustaceans (Figure 4A), and in the Cirripedia lineage it has undergone massive duplication, with remarkable tandem duplications of Br–C observed on contig3. There are 19 copies of Br–C in contig3 within a 191-kb length in the L. anserifera genome (Figure 4A). To further examine the massive expansion of Br–C in cirripeds, we preformed the ortholog search in two barnacle transcriptomes and identified more than 20 paralogs of Br–C in each species (Figure S4), suggesting its expansion in all cirripeds. During the life cycle of barnacles, there are two major metamorphosis events – one from nauplius to cypris and the other cypris to juvenile. Larval settlement involves the attachment of the cypris antennules onto the substratum by a cement, followed by ecdysis of the cyprid exoskeleton and metamorphosis (Høeg et al., 2012; Wong et al., 2018). Consistent with its key function in metamorphosis, the expansion of Br–C genes may facilitate the metamorphosis and settlement, enabling the drastic change from the larval to juvenile body plan in barnacles. The roles of Br–C genes should be investigated in the future by comparative gene expression analyses at different developmental stages.

Figure 4

Maximum likelihood trees of expanded (A) broad-complex core protein isoforms 1/2/3/4/5-like (Br–C) and (B) longitudinals lacking protein-like (Lola). (A) Br–C with 80 in Lepas anserifera and 63 in Amphibalanus amphitrite compared to 0–4 copies in other pancrustaceans (substitution model: Blosum62+G4), in which 19 copies (purple color) located in L. anserifera scaffold3 with three non-Br-C genes (gray color). (B) Lola with 29 in L. anserifera and 30 in A. amphitrite compared to 2–17 copies in other pancrustaceans (substitution model: LG+G4). Expression data of L. anserifera is shown in Figure S5. The first letter of the genus name and the first two letters of the species name are used in abbreviations, e.g., Lan for Lepas anserifera. Nodes with ≥70 bootstrap support are labelled with a black dot. Gene ID labeled of barnacles are in blue, other pancrustaceans in green and outgroup in black color.

In addition, cypris are highly mobile and equipped with organs for environmental sensing. We found an expansion of the longitudinals lacking protein-like (Lola) gene families in the L. anserifera genome (39 copies) and the A. Amphitrite genome (29 copies), compared to only two to seven copies in each of the other five compared pancrustacean genomes. Phylogenetic analysis of these Lola genes revealed several duplication events in the Cirripedia clade (Figure 4B). Lola was first identified in Drosophila as a broad-complex, tramtrack and bric-à-brac - zinc finger (BTB-ZF) family transcription factor regulating nervous system development. In Drosophila embryos, Lola regulated the axon tract development in the CNS, and axon guidance and growth in the peripheral nervous system (Giniger et al., 1994; Crowner et al., 2002). Although little is known about the function of this protein in crustaceans, previous phylogenetic analysis revealed that Lola is a new arthropod specific BTB protein that is only present in pancrustaceans (Quijano et al., 2016). A previous study of A. amphitrite cypris revealed a complex connection between peripheral sense organs and the CNS by a large number (2,000) of neurons, and the antennular nerve is the major peripheral nervous system (Harrison and Sandeman, 1999). Among the four antennular segments, the third segment contains an attachment disc with dense villi, and the fourth segment contains setae for detection of the physical and chemical properties of the substrate surface (Chan et al., 2018). Given the common origin of insect and crustacean Lola and the conserved nervous system features between insects and crustaceans (Osorio et al., 1997), the expansion of Lola genes and their expression in barnacles during different life stages are consistent with their purported roles in response to settlement cues and environmental signaling.

Unlike other pancrustaceans which grow by replacing the whole chitinous exoskeleton through periodic molting, cirripedes retain the outer and armored shell and grow by shedding only the exoskeleton between the soft tissues and the attachment base (Burden et al., 2014). The pupal cuticle protein 20-like (PCP20) homolog may have contributed to the unique growth pattern in cirripedes (Skinner et al., 1992). PCPs are proteins known to play an important role in the molting of various insects (Chihara et al., 1982; Nakato et al., 1990). For instance, in the silkworm Bombyx mori, it is synthesized in the epidemics during the metamorphosis from larvae to pupae, and then secreted into the pupal cuticle, with a peak of abundance before molting (Nakato et al., 1990). Under the search criteria of sequence identity >50% and E-value of 1e-10, we found 3 and 20 copies of PCP20 in L. anserifera and A. amphitrite, respectively, but no homolog of this gene in other pancrustaceans. Using more relaxed criteria (i.e., query coverage >30%, identity >30%, E-value of 1e-5), we found 90 PCP20 genes in L. anserifera, 231 in A. amphitrite, and 16 to 128 genes in each of the other pancrustaceans. Phylogenetic analysis revealed lineage-specific expansion of PCP20 in Cirripedia (Figure S6). The duplication of PCP20 and their active transcription (Figure S5) are consistent with the roles of these proteins as novel molting proteins uniquely found in barnacles.

Tissue-specific expression profiles and functions

To gain insight into tissue-specific functions, we analyzed the 16,356 genes expressed (TPM >1) in the peduncle (stalk), 17,861 genes in the cirri and 16,922 genes in the somatic body part (body). Among them, 13,804 genes were shared among the three tissues, which included many house-keeping genes such as elongation factor 1, ribosomal proteins and actin (Table S10). There were 879, 364 and 714 tissue-specific genes (TSGs) in the peduncle, cirri and body, respectively (Figure 5A and 5B). As our RNA-Seq analysis doses not have biological replication, the gene expression pattern of selected DEGs was further confirmed with real-time PCR (Figure 5C; Table S11). GO enrichment analyses were performed on the TSGs to gain insight into their biological functions (Figure 5D; Table S12). Scatterplots were also used to show tissue enriched biological processes. The flexible peduncle is armed with powerful longitudinal muscles, cement glands and covered by an unarmored cuticle (Power et al., 2010). Attachment of the peduncle to the substrate is maintained through continuous secretion from the cement gland (Anderson, 1994; Power et al., 2010). The enrichment of muscle development and biosynthetic related GO terms (e.g., muscle structure development and regulation of cellular macromolecule biosynthetic process) is consistent with the anatomy of the peduncle with strong muscle structure, and its role in secretion of adhesives for attachment to the substrate. The cirri were substantially enriched in cuticle development and metabolic pathways, which are related to its chitin-based exoskeleton and its role in filter-feeding. The somatic body parts which contained multiple organs, including the digestive, reproductive and nervous systems were enriched in an array of biological functions, such as innate immunity, metabolic process, response to stress and reproduction.

Figure 5

Tissue expression profile in Lepas anserifera. (A) Venn diagram of tissue-specific genes (TSGs). (B) Heatmaps showing expression of all TSGs with the hierarchical clustering. (C) Real-time PCR analysis. Relative fold changes of selected DEGs in stalk, cirri and other body part of L. anserifera were examined (data in Table S11). The heatmap shows the Z-scores (n = 3 for qPCR; n = 1 for RNA-Seq, in each tissue). (D) Biological processes (GO enrichment analysis, Table S12) enriched in TSGs. In the scatterplots, each bubble indicates a significantly enriched term (FDR <0.05) in a two-dimensional space derived by applying multidimensional scaling to a matrix of the GO terms' semantic similarities (Supek et al., 2011). Bubble color represents the Log10 FDR value. The circle size indicates the frequency of the GO term in the UniPort database (more general terms are represented by larger size bubbles).

Adaptation to sessile life by settlement inducing protein complex and cement proteins

Apart from gene family expansion, formation of new genes also provides genetic novelty. Previous studies have reported two settlement related proteins in barnacles, i.e., the settlement inducing protein complex (SIPC) and cement protein (CP). SIPC, a large glycoprotein acting as a pheromone, is critical for larval settlement and substance/barnacle species recognition (Dreanno et al., 2006a, 2007). Barnacle CP is used as a glue for settlement of the cyprids and permanent underwater attachment of the adults (Liang et al., 2019). Multiple adult CPs have been identified and named based on protein masses (i.e., CP19k, CP43k, CP52k and CP100k) (Power et al., 2010). Our phylogenetic analysis of pancrustacean SIPCs revealed a specific barnacle SIPC clade, including one L. anserifera gene and six A. amphitrite genes; these barnacle SIPCs are clustered with reported barnacle SIPCs (Figure 6A; also see Yorisue et al., 2012). These SIPC genes exhibited tissue-specific expression profiles (Figure 6C). Apart from its role in larvae settlement, previous study of SIPC in Balanus barnacle found its transcription in different developmental stages and tissues, and showed that it functioned in cuticle secretion (Dreanno et al., 2006b). Therefore, the phylogeny and tissue expression of SIPCs indicated that this gene has duplicated many times during the evolution of barnacles, and different paralogs of this gene may have co-opted into different organs with functions that are related to settlement and cuticle formation.

Figure 6

Maximum likelihood trees of (A) settlement inducing protein complex (SIPC) and (B) cement protein (CP) in pancrustacean genomes and barnacle sequences downloaded from UniProt. (A) SIPC with 18 in Lepas anserifera and 29 in Amphibalanus amphitrite compared to 3–29 copies in other pancrustaceans (substitution model: WAG+F+G4). Barnacle-specific clade highlighted in orange color. (B) CP with 9 in L. anserifera and 18 in A. amphitrite, but only one homolog in Daphnia (substitution model: WAG+F+G4). (C) Tissue expression data of SIPC and CP in L. anserifera. The first letter of the genus name and the first two letters of the species name are used in abbreviations, e.g., Lan for Lepas anserifera. Nodes with ≥70 bootstrap support are labelled with a black dot. Gene ID labeled by barnacles in blue, other pancrustaceans in green and outgroup in black color.

Phylogenetic analysis showed that the CPs are clustered by their protein masses (Figure 6B), and there are long branch lengths (substitution) between these CP clades, which indicated episodes of historical duplication and diversification of this gene in Cirripedia. Among pancrustaceans, only Daphnia possesses one homolog of CP - a protein called Filaggrin involved in filament-aggregation, epidermal barrier formation against antimicrobial attack and damage by UV radiation (Harding et al., 2013). Our transcriptome analysis showed that the CPs were specifically expressed in the peduncle (Figure 6C) where the cement gland is located, indicating its function in permanent attachment. Thus, our results indicated that barnacle CPs may have been co-opted from a crustacean gene with multiple functions including cuticle secretion and immunity. This hypothesis can be tested using phylogenetic methods when more genomes from barnacles and their relatives become available.

Conclusions

In summary, we have provided a comprehensive genomic resource and established a phylogenetic framework for understanding the molecular adaptations of cirripedes to the sessile lifestyle. Gene family analysis revealed genetic mechanisms for environmental sensing, CNS development and cuticle formation. The formation of SIPCs and CPs as new genes and their duplications may have enabled barnacle underwater attachment and survival in habitats with strong wave action. Our results will be valuable for the development of antifouling strategy against barnacle settlement, and underwater adhesives that mimic the functions of barnacle cement proteins. Moreover, given the general lack of genomic resources in pancrustaceans, the L. anserifera genome will facilitate comparative studies aiming to understand the evolution of this large group of arthropods.

Declarations

Author contribution statement

Jack Chi-Ho Ip: Performed the experiments; Analyzed and interpreted the data; Wrote the paper.

Jian-Wen Qiu & Benny K. K. Chan: Conceived and designed the experiments; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Funding statement

Jian-Wen Qiu was supported by the National Key R&D Program, Ministry of Science and Technology, China (2018YFC0310005), Southern Marine Science and Engineering Guangdong Laboratory (Guangzhou) (GML2019ZD0409), Hong Kong Branch of Southern Marine Science and Engineering Guangdong Laboratory (L20190005), General Research Fund from the (12302917), and State Key Laboratory of Environmental and Biological Analysis, Hong Kong Baptist University (SKLP_1920_P04).

Benny K. K. Chan is supported by a grant in Senior Investigator Award, , Taiwan (AS-IA-105-L03).

Data availability statement

The genomic DNA and RNA sequences have been deposited in NCBI Sequence Read Archive under Bioproject No. PRJNA678024. The genome assemblies and annotations are available from the Figshare DOI:10.6084/m9.figshare.12831926.

Declaration of interests statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.