Individual barnacles were collected in the Indian Ocean and the West Pacific from habitats such as mid-oceanic ridges, basins, trenches, hydrothermal vents, seamounts, shelf seas, and tidal zones. Of these, the deep-sea samples were collected during multiple expeditions led by the Institute of Oceanology Chinese Academy of Sciences, the Second Institute of Oceanography Ministry of Natural Resources, the Institute of Deep-sea Science and Engineering Chinese Academy of Sciences and Tongji University. The samples were obtained using a variety of remotely operated vehicles and manned submersibles. All samples were preserved in either 75% ethanol or at -80°C, and deposited in the Marine Biological Museum Chinese Academy of Sciences, Qingdao, China, and the Sample Repository of the Second Institute of Oceanography Ministry of Natural Resources, Hangzhou, China ( Table 1 ).

Total genomic DNA was extracted from pedunculate or abdominal muscle using a modified cetyltrimethylammonium bromide method (Allen et al., 2006), and a 500-bp paired-end library was constructed using the NEBNext Ultra DNA Library Prep Kit for Illumina (NEB, Ipswich, United States). Next-generation sequencing was performed on the Illumina NovaSeq 6000 platform (BIOZERON Co., Ltd, Shanghai, China).

The raw data, 150 bp paired-end reads, were filtered by removing the following: reads with adaptors, reads with a quality score of below 20 (Q<20), reads containing a percentage of uncalled bases (“N” characters) that was equal to or greater than 10%, and duplicated sequences. A combination of de novo and reference-guided assemblies was used to reconstruct the mitochondrial genome. The filtered reads were first assembled into contigs using GetOrganelle v1.6.4 (Jin et al., 2020). Next, the contigs were aligned with mitogenomes of closely related barnacle species using BLAST searching, and were also optimized using scaffolds from SPAdes 3.13.0 (Bankevich et al., 2012). Finally, the assembled sequence was reordered and oriented according to the reference mitochondrial genomes, thus generating the final assembled mitochondrion genomic sequence.

The assembled mitogenomes were annotated using the online MITOS tool (http://mitos.bioinf.unileipzig.de). The default parameters were used to predict protein coding genes (PCGs), transfer RNA (tRNA) genes, and ribosome RNA (rRNA) genes. These genes were further verified using the Open Reading Frame Finder on NCBI (https://www.ncbi.nlm.nih.gov/orffinder), tRNAscan-SE 2.0 (Chan & Lowe, 2019) and rRNAmmer 1.2 (Lagesen et al., 2007), respectively. Manual corrections of genes for start/stop codons were performed in SnapGene Viewer (https://www.snapgene.com) by referencing the reference mitochondrial genome. The circular mitochondrial genome map was drawn using the CGview Server (Grant & Stothard, 2008).

The analyses were achieved using the mitogenomes of the 13 species obtained in the present study, as well as the open-access mitogenome datasets of other barnacle species, including an iblomorph, three asymmetrical barnacle species (including one chemosynthetic species), 17 stalked barnacle species (including three chemosynthetic species), and 31 acorn barnacle species (including two chemosynthetic species). A first-pass phylogenetic inference based on all the functional genes, (i.e., concatenated 13 PCGs as well as 22 tRNA and two rRNA sequences) was conducted using the Flowchart analyses of PhyloSuites v1.2.2 (Zhang et al., 2020) with the maximum-likelihood optimality criterion. However, this produced a phylogenetic tree that was poorly supported, possibly owing to noise from gene heterogeneity and mutational saturation. Next, we retained only the 13 PCGs. An inspection of the substitution status of these 13 PCGs using DAMBE v5.3 (Xia & Xie, 2001) indicated that all had mutational saturation at the third codon. Thus, we excluded the third codons in subsequent analyses on the concatenated nucleotide data. Bayesian inference (BI) and maximum likelihood (ML) methods were used to analyze the data from both the concatenated nucleotide and the putative amino acid sequences of the 13 PCGs.

Nucleotide and amino acid sequences of each 13 PCGs were aligned using MAFFT v7 (Katoh & Standley, 2013) with default parameters, and the nucleotide sequences were aligned in codon manner. MACSE v2 (Ranwez et al., 2018) and Gblocks (Talavera & Castresana, 2007) were used to optimize alignments and to eliminate segments that were divergent and ambiguously aligned. Finally, the trimmed nucleotide and amino acid sequences were concatenated into independent dataset respectively. The best partitioning schemes and best-fit substitution models of each partition in the concatenated dataset were determined using ModelFinder (Kalyaanamoorthy et al., 2017) and PartitionFinder 2 (Lanfear et al., 2017) according to the Bayesian information criterion. All the above analyses were performed using the integrated and scalable desktop platform PhyloSuites v1.2.2. The ML tree was reconstructed using RaxML-NG (Kozlov et al., 2019), and branch support was assessed by transfer bootstrap expectation (TBE) with 1,000 replicates (Lemoine et al., 2018). The Bayesian analysis was conducted using a parallel version of Mrbayes v3.2.6 (Ronquist et al., 2012). Two independent runs were executed with four Markov Chain Monte Carlo (MCMC) for 10 million generations (beginning with a random tree, chains were sampled every 10,000 generations). Convergence of the analyses was ensured by monitoring the likelihood values and the standard deviation of partition frequencies (<0.01). The resultant outputs were further validated by using Tracer 1.7.1 (Rambaut et al., 2018) to diagnose the effective sample size (ESS) values of all sampled parameters. The first 25% of trees generated prior to the achievement of stationarity of the log-likelihood values were discarded as burn-in. The remaining trees were used to construct a 50% majority rule consensus tree and to estimate the posterior probabilities (PP) simultaneously. FigTree (http://tree.bio.ed.ac.uk/software/figtree/) was used to visualize the phylogenetic trees. Linearized mitochondrial gene arrangement patterns were superimposed on the phylogenetic trees using the iTOL webserver (Letunic & Bork, 2021).

To test the hypothesis that sessile barnacle species (i.e., verrucomorphs and balanomorphs) form a monophyletic group, we performed a likelihood-based approximately unbiased (AU) test in IQ-TREE v2.0 (Nguyen et al., 2015) with 1,000 replicates to assess the P-value of the test topology.

The divergence times of barnacle lineages were estimated from the partitioned nucleotide dataset using BEAST v2.6.3 (Bouckaert et al., 2014). Variation in mutation rates amongst branches was allowed by assuming an uncorrelated relaxed lognormal molecular clock model. The birth-death model (Gernhard, 2008) was set as the tree’s priors. Unlinked site models were set for each partitioned dataset to allow independent substitution rates and base frequencies, and determined by bModelTest (Bouckaert & Drummond, 2017) implemented in the BEAST v2.6.3. Six fossil calibration points [C1: Praelepas jaworskii, 306.5–311.7 million years ago (MYA); C2: Calantica (Scillaelepas) ginginensis, 83.5–85.5 MYA; C4: Arcoscalpellum fossula, 70.6–89.3 MYA; C5: Pollicipes aboriginalis, 83.5–85.8 MYA; C9: Pachydiadema (Catophragmus) cretacea, 70.6–89.3 MYA; C11: Tetraclitella sp. cf. purpurascens, 20.4–23.0 MYA] were selected to calibrate the divergence time estimation; all points had been tested with high confidence in previous studies (Pérez-Losada et al., 2008; Linse et al., 2013; Herrera et al., 2015). Fossil ages were used as lower boundary constraints, and their prior exponential distribution were assumed with estimated mean values. The default prior distribution settings were assumed for all other parameters. The MCMC analysis was run for 10 million generations with a sampling frequency of 10,000 generations. Convergence diagnosis was examined by Tracer. After the first 25% of trees were discarded as burn-in, the maximum clade credibility tree with median nodal height was generated using TreeAnnotator 2.6.3 (Bouckaert et al., 2014). The time-scaled tree was mapped against geological time and visualized using the R packages phytools (Revell, 2012), PHYLOCH (http://www.christophheibl.de/Rpackages.html), strap (Bell & Lloyd, 2015), and CODA (Plummer et al., 2006).

Signatures of natural selection on the 13 PCGs were detected using the codon-based maximum likelihood (CodeML) implemented in the PAML 4 package (Yang, 2007), based on the tree that was reconstructed with the concatenated amino acid sequences. Both branch and branch-site models were used to identify the variation in selective pressures acting on the branch lineages and individual sites of PCGs. For the branch model, a one-ratio model and a two-ratio model were used to estimate the ω values across the whole lineage and foreground lineages (i.e., the chemosynthetic barnacles). The branch-site model assumed ω values to be alterable among sites in the proteins as well as across branches on the tree, and could detect any positive selection affecting a few sites along the foreground branches. Bayes Empirical Bayes (BEB) analysis was used to calculate the Bayesian posterior probability of the positively selected sites (ω>1). All the assumed models were compared against their corresponding null hypotheses using the Likelihood Ratio Test (LRT); the LRT P-values were then calculated using a Chi-square (χ2) analysis. Models with LRT P-values less than 0.05 and BEB larger than 95% were considered as the best-fit models for evolution and positive selection.

To achieve a more intuitive understanding of the significance of the sites under positive selection, we mapped the sites onto the three-dimensional (3D) structure of mitochondrial proteins. The SWISS-MODEL (Waterhouse et al., 2018) server and I-TASSER (Yang et al., 2015) server was used to model the 3D structure of the target proteins based on the sequence reconstructed using the PAML 4 package. Next, PyMOL (The PyMOL Molecular Graphics System, V2.5, Schrödinger, LLC.) was used for visualizing and depicting the positive selection sites in the proteins.

The final aligned datasets for nucleotide and amino acid sequences consisted of 7,178 bp and 3,563 characters, respectively. The best partitioning schemes and best-fit substitution models as determined by ModelFinder and PartitionFinder are listed in the Supporting Information, Supplementary Table 1 . Tree topologies resulting from the BI and ML analyses were identical for the same dataset, and only a few weakly supported internal nodes of Balanidae and Calanticidae were different between the two datasets. Furthermore, the topology that was revealed from the amino acid dataset ( Figure 2 ) received higher nodal support than that from the nucleotide dataset ( Supplementary Figure 1 ). But consistently, all orders received monophyletic support, and all families formed a monophyletic group except Balanidae in the phylogenetic analyses. Fortuitously, a conflict between the topologies of the nucleotide and amino acid datasets was found in the position of Smilium sinense. This species was recovered as a sister group of verrucids in the nucleotide topologies, which contrasted its coordinal species Calantica studeri in the amino acid topologies as well as in the topology generated by all the functional genes. Furthermore, the nodal support for S. sinense in nucleotide topologies was very poor (PP=0.55, BP=59) ( Supplementary Figure 1 ). Comparisons across different barnacle species revealed that they have a highly conserved gene order that corresponded to their families, except for Chthamalidae and Calanticidae, which displayed conspicuous translocations or inversions of some gene blocks ( Figure 3 ).

Most of the phylogenetic relationships among species recovered from the present analyses were consistent with the results of previous studies based on multilocus molecular markers (Pérez-Losada et al., 2008; Linse et al., 2013; Rees et al., 2014; Herrera et al., 2015; Lin et al., 2015). Barnacles from deep-sea chemosynthetic ecosystems can be divided into two strongly supported clades (PP=1, BP=100). Clade A, containing Ashinkailepas seepiophila, Leucolepas longa, Neoverruca intermedia and Vulcanolepas fijiensis, was well supported as the sister group to the scalpellids clade. Clade B, containing Eochionelasmus coreana and E. ohtai, was strongly supported as the basal clade of the balanomorphs ( Figure 2 , Supplementary Figure 1 ). This result was highly congruent with the findings of Herrera et al., (2015), who suggested that deep-sea chemosynthetic barnacles have colonized hydrothermal vents and/or cold seeps at least twice in the history of thoracican evolution. The verrucomorph clade occupied the basal position in the topologies of all the trees recovered ( Figure 2 ; Supplementary Figure 1 ), and the hypothesis that sessile barnacles (verrucomorph and balanomorph) form a monophyletic group was clearly rejected by the AU tests (P-value=3.73e-07 in the amino acid dataset, P-value=8.67e-06 in the nucleotide dataset).

The topology of the tree reconstructed from the divergence time estimation by BEAST2 was nearly identical to the topologies of trees recovered from the BI and ML analyses based on the nucleotide and amino acid data. Minor differences also appeared in a few shallow nodes of the Balanidae and Tetraclitidae. Divergence times of the barnacle lineages estimated under the prior settings are shown in Figure 4 . Nearly all the node times lay between the corresponding values estimated by Linse et al. (2013) and Herrera et al., (2015). The occurrence time of the most recent common ancestor (MRCA) of Clade A was estimated at 46.8 million years ago, with the 95% highest posterior density interval (HPD) occurring at 25.5–74.1 MYA. This clade diverged from its sister clade at 168.4 MYA (95% HPD: 126.3–220.4). Within Clade A, the origin of the Neolepadidae was estimated to have occurred 29.3 MYA (95% HPD: 13.7–49.2). The MRCA of Clade B was older than that of Clade A, had occurred 65.3 MYA (95% HPD: 34.6–101.7), where it split from other balanomorphs at 124.6 MYA (95% HPD: 91.2–158.9). In general, the extant deep-sea chemosynthetic barnacles most probably originated from the early Paleogene or upper Cretaceous, and speciated during the subsequent periods.

To test if any selective pressure occurred specifically in chemosynthetic barnacle lineages, we compared the alternative model (a two-ratios model, with different lineages assigned different ω values) and the null model (a one-ratio model, with all lineages assigned the same ω value) by labelling the foreground branches in two different ways. The first involved labelling the entire clades of Clade A and Clade B respectively; the second involved labelling the individual branches of Clade A and Clade B. Both approaches produced very significant results (P<10^-8) ( Table 2 ). The results also showed that chemosynthetic barnacle species had a higher ω value than other barnacle species ( Table 2 , Figure 5 ), suggesting they had accumulated more nonsynonymous substitutions throughout their evolution and their adaptation to extreme environments. Notably, the ω values of the common ancestor to Clade A and Clade B were distinctly high. This could indicate that the ancestors of chemosynthetic barnacles had undergone relatively high selective pressures when they colonized the chemosynthetic biotope during their emergence.

Studies suggest that positive selection always occurs episodically at a few sites or lineages in response to variation of the surrounding environment (Yang, 2007; Shen et al., 2010). Thus, we performed a series of branch-site model tests (ω>1 along the foreground branch vs. ω≤1 along the foreground and background branches) among individual branches of Clade A and Clade B. However, only the results of tests along branches leading to Clade A, Clade B, internal branches of Clade A and the Eochionelasmus ohtai branch received statistically significant support (P<0.05) ( Table 2 ). This pattern reflected the pattern observed in the results of the branch model analysis. BEB analysis showed that positive selection sites occurred in all mitochondrial PCGs of Clade A, except in cox2 and atp8. In contrast, for Clade B, positive selection sites occur in cox3, atp6, nad1, nad2, nad3, nad4 and nad5. No common positive selection site between the two lineages was found ( Table 2 ). The predicted 3D structure of the proteins revealed that these positive selection sites appeared as both conformation of α-helix structures and random coils ( Supplementary Figure 2 ).

Our findings pertaining to barnacle evolutionary history largely corroborated those of previous studies, particularly with regards to the independent origins of the two chemosynthetic barnacle lineages. However, our findings also suggested alternative evolutionary scenarios for other thoracicans. Most apparent among these was the phylogenetic position of verrucomorph barnacles. Our analysis indicated that verrucomorphs were basal to all barnacles except iblomorphs ( Figure 2 ; Supplementary Figure 1 ). This partly corresponded to the initial speculation that verrucomorphs and balanomorphs were of independent origins (Pilsbry, 1916; Newman & Ross, 1976). The sister relationship between verrucomorphs and balanomorphs was clearly rejected by the monophyletic test in our analysis. It had also not been confirmed in the mitochondrial genome analyses of Kim et al. (2018); Kim et al., (2019) and Lee et al. (2019), or the multiple gene analyses of Linse et al. (2013); Rees et al. (2014), Herrera et al., (2015), and Kim et al. (2022). Furthermore, the mitogenome arrangement of verrucomorph species, which displayed unique gene blocks ( Figure 3 ), was distinctly different from that of the balanomorphs and other barnacles. The ontogeny of verrucomorph species (Verruca stroemia) also differed from that of balanomorph species (Semibalanus balanoides) owing to a lack of lateral plates in individual development (Newman, 1989). Although limited evidence of morphology and fossil records to fully support current phylogenetic result, there were some clues to speculate that verrucomorphs might originate from archaeolepadomorphs. Archaeolepadomorphs (such as Archaeolepas and Bosquetlepas) and verrucomorphs only possess six plated capitula (i.e., the paired scuta and terga, a carina, and a scutum) and lack any lateral plates, the arrangement thus provides the most parsimonious evolutionary interpretation. Moreover, our analysis indicates that the origin of verrucomorphs is in the upper Triassic ( Figure 4 ), which corresponds with the divergence time of the archaeolepadomorphs indicated by Chan et al. (2021). However, considering that incongruent phylogenetic position of verrucomorphs was present in molecular phylogeny. The natural evolutionary history of verrucomorphs must be further verified with more evidence from added data analysis.

Another controversial issue in thoracican phylogeny concerns the calanticids and pollicipedids, which were traditionally placed in Scalpelliformes, until Chan et al. (2021) proposed the orders Calanticomorpha and Pollicipedomorpha for them, however, their phylogenetic status is uncertain owing to the incongruent results of separate molecular analyses. (Pérez-Losada et al., 2008; Linse et al., 2013; Rees et al., 2014; Lin et al., 2015; Herrera et al., 2015). Our mitogenome data indicate that the pollicipedids diverged immediately before the balanomorphs, and were followed closely by the calanticids ( Figure 2 , Supplementary Figure 1 ) with robust support. These findings coincide with the Bayesian results of Linse et al. (2013), as well as the relatively constant position occupied by the pollicipedids in all mitochondrial genome-based phylogenetic analyses (Song et al., 2017; Cai et al., 2018; Kim et al., 2018; Kim et al., 2019; Chen et al., 2019; Lee et al., 2019; Tian et al., 2020). The mitogenome arrangement of calanticids and pollicipedids more closely resembles the patterns of balanomorphs, except for Smilium sinense which shows remote inversion events of the gene blocks, C–Y–nad1–L1–rrnL–V–rrnS and I–M ( Figure 3 ). Anderson (1983) demonstrated that the evolution of balanomorphs from calanticids was functionally feasible. Gale (2014b) indicated that Etcheslepas was morphologically intermediate between Capitulum and Pycnolepas, which suggested that pollicipedomorphs and brachylepadomorphs shared a close relationship. He further proposed a phylogenetic hypothesis for balanomorphs based on their morphology and paleontology, with calanticomorphs and pollicipedomorphs located at the basal position (Gale & Sørensen, 2014). In the morphometric analyses on antennule of cyprid larvae, Al-Yahya et al. (2016) showed similar structures among pollicipedomorphs and balanomorphs. Collectively, the evidence indicates that balanomorphs, pollicipedomorphs and calanticomorphs might share a common ancestor before radiating into the respective lineages observed today.

Our result indicates the chemosynthetic barnacles diverged from the early Cenozoic ( Figure 4 ). A similar divergence time has been proposed by Linse et al. (2013) and Herrera et al. (2015), and the divergence time of neolepadoids is also in accordance with the fossil records (Carriol et al., 2016; Gale et al., 2020). During the same period, other marine chemosynthetic taxa also emerged and radiated. For instance, mytilid mussels 45 MYA (Lorion et al., 2013), alvinocarid shrimp 48.4–55.9 MYA, bythograeid crabs 51.5–69.7 MYA (Yang et al., 2013) and siboglinid tubeworms 60 MYA (Chevaldonné et al., 2002). These divergence times occurred just posterior to the Cretaceous-Paleogene mass extinction (K-Pg extinction), which caused a near 50% decline in marine generic diversity (Sepkoski, 1996; D’Hondt, 2005) and over 50% of cirripede genera extinction (Gale & Sørensen, 2014). Jacobs & Lindberg (1998) have proposed that the oceanic anoxic/dysoxic events occurring at the Cenomanian/Turonian boundary as well as the latest Paleocene led to the elimination of much of the contemporary deep-sea benthos, as well as vent faunas; this argument has also been supported by Little & Vrijenhoek (2003). Accordingly, an abundance of empty ecological niches would have been created following these mass extinctions and anoxic/dysoxic events. Hsü et al. (1982); Zachos et al. (1989), and D’Hondt (2005) showed that the flux of organic matter to the deep seafloor declined tremendously during the K-Pg extinction. All these ecological events may have driven food scarcity among the contemporary deep-sea benthos. Hence, it is possible that the colonization of deep-sea chemosynthetic habitats by chemosynthetic animals was promoted by ecological niche vacancy and food scarcity. Furthermore, hydrothermal vents or cold seeps are widespread, energetic and highly productive ecosystems supported by chemoautotrophic microbes rather than by the sinking of photosynthesis-based organics (Grassle, 1987; Little & Vrijenhoek, 2003; Baker et al., 2016). Nonetheless, these extreme habitats are also characterized by hypoxic, chemical release, and more extreme temperature gradients in comparison with the surrounding deep-sea environments (Hourdez & Lallier, 2007; Le Bris et al., 2017). These conditions may have promoted the development of specialized adaptations in the genomes of chemosynthetic species.

Hypoxic, sulfide, methane and heavy metals suffused in the chemosynthetic ecosystem are toxic for the mitochondrial electron transfer chain, and can inhibit aerobic metabolism (Grieshaber & Völkel, 1998; Belyaeva et al., 2008). Mitogenome encode key enzymes involved in the process of aerobic respiration. Changes in survival environments will create selective pressures on genome to confront alternative environmental conditions, and fix nonsynonymous mutations as adaptive substitutions for species to cope with the new demands. Our results show that the mitogenomes of chemosynthetic barnacles encountered relatively high selective pressures, especially in ancestral taxa, which were the first to colonize chemosynthetic habitats ( Figure 5 , Table 2 ). Remarkably, the ω value of the ancestor of Clade B (Eochionelasmus) is higher than that of Clade A (neolepadoids). It is possible that the MRCA of neolepadoids was already a deep-sea dweller before it colonized deep-sea chemosynthetic habitats ( Figure 2 ; Linse et al., 2013; Herrera et al., 2015). In contrast, Eochionelasmus is nested within the shallow-water clade ( Figure 2 ; Pérez-Losada et al., 2014), indicating that its ancestor invaded into deep-sea chemosynthetic habitats directly from shallow waters. This means that the MRCA of Eochionelasmus likely encountered more selective constrains than that of the neolepadoids. This may also have accounted for the apparent differences in species diversity between the two chemosynthetic barnacle lineages (i.e., three species vs. over 15 species). Positive selection in specific sites of a functional protein is typically taken as evidence of an adaptive substitution responding to environmental changes (Yang & Nielsen, 2002; Rausher & Huang, 2016). The branch-site model analysis obtained a similar pattern as that from the branch model analysis, most of the adaptive substitutions occurred in their MRCAs of chemosynthetic barnacles ( Figure 4 , Table 2 ). The 3D structure shows that these adaptive substitutions are scattered throughout both α-helix and random coil conformations, and that some sites are located precisely in the transformation point of these two conformations, such as 245Q of the Cytb in the neolepadoid lineage, and 210G of the Nad4 of the Eochionelasmus lineage ( Supplementary Figure 2 ). Meanwhile, some amino acids experienced radical substitutions in their physicochemical properties. For instance, the 104th site of the neolepadoids’ Nad2 arose as a substitution of a hydrophobic and nonpolar Valine with a hydrophilic and polar Threonine. These substitutions would have caused dramatic changes to protein conformation, in turn impacting their potency, just like the facts verified by physiological and biochemical experiments (Paulus et al., 2013; Laye et al., 2017). Interestingly, our findings also suggest that there is no identical positive selection site shared by neolepadoids and Eochionelasmus ( Supplementary Figure 2 ), which may imply that genomic response to environmental change not only depends on external factors but also relates to themselves (i.e., lineage-specific). Studies have shown that diversifying selection on physiological process can reduce intraspecific gene flow and promote population divergence, ultimately incurring speciation (Ellegren & Sheldon, 2008; Scott et al., 2011; Morales et al., 2015).

No studies have specifically discussed the morphological or physiological adaptations of chemosynthetic barnacles for surviving in extreme environments, but Yorisue et al. (2013) showed potential strategies of larval development to settle down in proper habitats developed by chemosynthetic barnacles. In comparison with their deep-sea or shallow-water relatives, chemosynthetic barnacles have a relatively larger number of setae on the intermediate segments of the posterior three cirri (i.e., generally more than six pairs vs. less than six pairs), and a higher ratio of setal length to the width of the segment where the setae are rooted (generally 2.8–14.1 vs. mostly 2.0–3.5) (primary data summarized from Yamaguchi and Newman, 1990; Buckeridge, 1998; Newman, 2000; Southward & Jones, 2003; Buckeridge et al., 2013; Watanabe et al., 2021). We speculate that these features represent morphological adaptations of chemosynthetic barnacles, as the dense and long setae would enhance the volume of water currents flowing into the mantle cavity when their cirri are beating for respiration and feeding. Buckeridge et al. (2013) have shown that the powerful morphological plasticity of Vulcanolepas scotiaensis was triggered by microenvironmental variation in chemosynthetic habitats. Several chemosynthetic barnacle species, such as V. osheai, V. buckeridgei and V. verenae, had developed ectosymbioses with sulfur-oxidizing bacteria for feeding (Southward & Newman, 1998; Watanabe et al., 2021). Furthermore, all barnacle somatotypes (e.g, sessile, stalked and asymmetric) were developed in chemosynthetic habitats. Together, these phenomena suggest that a subdivision of the ecological niche may have also been conducive to the speciation and diversification of chemosynthetic barnacles.