Samples of the deep-sea barnacles from a hydrothermal vent Longqi hydrothermal field on the Southwest Indian Ridge (49°64’E, 37°17’S) at a depth of 2,759 m in the Southwest Indian Ocean collected on board of “Explorer 1” scientific research ship of Institute of Deep-sea Science and Engineering (Chinese Academy of Sciences) by using “Deep-sea Warrior” manned submersible during the winter of 2018, we thank the crew for their help during the cruise. Samples of the shallow-water barnacles were obtained from the Guangxi Beibu Gulf (109° 12 ‘ E, 21° 04 ‘ N), depth of about 300 m. In contrast to deep-sea species, the depth of 300m can be considered relatively shallow. Therefore, we provisionally identified the two barnacle species as the deep water species Vulcanolepas fijiensis (Lee et al., 2019) and the shallow water Scalpellum stearnsi (Chan et al., 2010). As a result of the Chen et al. Commentary Paper (Chan et al., 2024) we became aware that our initial identification of V. fijiensis was incorrect. We have confirmed that the barnacle is actually Neolepas marisindica and have revised this paper to include the correct species indentation. After collection, they were immediately stored at -80°C.

According to the manufacturer’s directions, the TIANamp Marine Animal DNA Kit (TIANGEN Biotech, China) was used to extract DNA from the abdominal muscle. Then, 1% gel electrophoresis and ultraviolet spectroscopy (Eppendorf BioPhotometer D30, Eppendorf, Germany) were used to estimate and quantify the genomic DNA. Then, the DNA was randomly cut into 350 bp fragments, which were repaired at the end and amplified by a polymerase chain reaction to construct the sequencing library. Next, Sanger dideoxy sequencing (N. marisindica) was conducted on an ABI-3730XL (Thermo Fisher Scientific, USA), and high-throughput sequencing (S. stearnsi) was conducted on Illumina NovaSeq 6000 platform (TSINGKE Biotechnology Co, Ltd, China).

For each species of barnacle (N. marisindica and S. stearnsi), a random selection of three individuals was made, and all RNA was extracted using TRIzol reagent (Invitrogen, USA). The quality of the RNA was assessed through 1% agarose gel electrophoresis. The purity and concentration of the RNA were determined using a NanoPhotometer spectrophotometer (IMPLEN, Germany) and Qubit 2.0 Fluorimeter (Life Technologies, USA), respectively. The integrity of the RNA was evaluated using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, USA). The sequencing library was established according to the manufacturer’s instructions of the NEBNext Ultra RNA Library Prep Kit (NEB, USA). The mRNA was enriched using poly-T oligo-attached magnetic beads, and an Agilent Bioanalyzer 2100 system (Agilent Technologies) was used to evaluate the synthesis of the complementary DNA (cDNA), the construction of a pair-end library, and the quality of the library. Finally, the two different kinds of species’ cDNA libraries were clustered and sequenced using the Illumina HiSeq 2500 platform by NovoGene (Beijing, China).

The initial dataset of the mitochondrial genome consisted of 150 base pairs at the paired ends. Reads containing an adapter, reads with a mass score below 20 (Q< 20), and duplicated sequences were removed. A reference genome (Lepas anserifera) was used to retrieve the reads belonging to the mitochondrial genome for assembly. Then, de novo assembly was performed using SPAdes 3.13.0 (Bankevich et al., 2012), and the assembly results were corrected for mutual alignment to confirm the accuracy of the assembled genome sequence. The assembled mitochondrial genomes were subjected to analysis using the online MITOS software (http://mitos.bioinf.unileipzig.de) with default parameters. This software was utilized to predict the PCGs, transfer RNA (tRNA) genes, and ribosomal RNA (rRNA) genes. Manual corrections to the genes for the start/stop codons were performed in Sequin 15.10 (http://www.ncbi.nlm.nih.gov/Sequin/index.html) by referencing the mitochondrial genome. On this basis, BLAST (https://blast.ncbi.nlm.nih.gov/Blast.cgi) technology was used to compare the existing research data including the PCGs and rRNA to verify and correct the accuracy of the analysis results. The circular mitochondrial genome map was drawn using OGDRAW 1.3.1 (Greiner et al., 2019), and the complete mitogenomes of N. marisindica and S. stearnsi have been registered as MZ772032 and OP345466 on GenBank.

The biases for the codon usage of the mitogenome can regulate genes and reflect the evolutionary relationship between different species (Wei et al., 2014). Consequently, the usage of codons was analyzed for the PCGs and rRNA. Using the following formula, the inclination related to the nucleotide bases was calculated:

The low-quality reads and reads containing adaptors were removed from the raw reads using the internal software (ng_qc) with custom script of NovoGene, and the transcriptome clean reads were obtained. The de novo transcriptome assembly was performed using Trinity v2.5.1, with the minimum k-mer parameter coverage of 2, while default parameters were employed for other sections (Grabherr et al., 2011). The longest transcript of each transcription group was regarded as a unigene for the subsequent analysis. Next, BUSCO v1.22 was used to verify the quality of the assembly that was performed using the arthropoda_odb9 database (Simão et al., 2015). Additionally, all the unigenes were annotated by blasting the common databases, such as the National Center for Biotechnology Information non-redundant protein (Nr), Swiss-Prot, and euKaryotic Ortholog Group using DIAMOND v0.8.22 (Buchfink et al., 2015), and the NCBI nucleotide database using BLAST v2.2.28+ (Altschul et al., 1990). The Kyoto Encyclopedia of Genes and Genomes (KEGG) classification was performed using the KEGG Automatic Annotation Server (Moriya et al., 2007). The protein family (Pfam) alignments were put into practice using HMMER v3.0 (http://hmmer.org/) and were based on the results of Nr and Pfam, and Blast2GO v2.5 was used to classify the Gene Ontology (GO) (Götz et al., 2008). The Nr database was searched for unigenes with BLASTx, and then SOLAR was used to compare the link fragments (Rota-Stabelli et al., 2010), so all the open reading frames (ORF) of the unigenes were obtained and converted into amino acid sequences.

The target sequences for this analysis consisted of the 13 PCGs from N. marisindica and 10 other barnacle species, aligning using clustalW (codons) in MEGA v7.0. Manual removal of gap and stop codons was performed to ensure accurate alignment (Kumar et al., 2016). The resulting sequence alignment file was then utilized for the analysis of positive selection. A molecular phylogenetic tree was constructed using the maximum likelihood method in MEGA v7.0 to assess selective pressure. Subsequently, the reconstructed tree was subjected to EasyCodeML in the PAML software (Gao et al., 2019) to identify natural selection characteristics within the 13 PCGs. Among them, the 13 PCGs from N. marisindica were treated as the foreground branch, while PCGs from other barnacle species served as the background branch. Site models were employed to identify the individual sites of PCGs. Changes in the different sites of the 13 PCGs were explored by calculating the ω values (ω = K_a/K_s or dN/dS). For the positive selection sites (ω > 1), the Bayes Empirical Bayes (BEB) method was used to estimate the Bayesian posterior probabilities (Biswas and Akey, 2006). The likelihood Ratio Test (LRT) was applied to compare the hypothetical models with the corresponding original hypotheses, and the LRT P-value was determined using the Chi-square method. The sites of the PCGs with an LRT P-value< 0.05 or BEB > 95% were considered to be under evolution and positive selection pressure. The ProtParam program (http://www.expasy.ch/tools/protparam.html) software was utilized to compute the relative molecular weight and isoelectric points of the sequence, and SignalIP 4.1 (http://www.cbs.dtu.dk/services/SignalP/) was employed to predict the signal peptides. Additionally, the SMART program (http://smart.embl-heidelberg.de/) was utilized to predict the conserved functional domains of proteins (Letunic and Bork, 2021). On this basis, we matched the sites to the three-dimensional (3D) conformation of the PCGs to observe their role under positive selection more clearly. The 3D structure of the target proteins was simulated using the SWISS-MODEL (Waterhouse et al., 2018) server, which visualized the positive selection sites in the proteins.

Transcriptome sequence set from the two species (N. marisindica and S. stearnsi), internal reference (Amphibalanus Amphitrite) and external reference species (Eurytemora affinis) were utilized for the analysis of the orthologous genes and nonsynonymous (dN)/synonymous (dS) substitution rates. Subsequently, adaptive evolution was assessed by comparing the ω value (dN/dS). BLAST was employed to pair the orthologs with the best hit probabilities. Following this, sequence alignment was conducted using MUSCLE, and all gaps were eliminated. The ω numerical calculation for each ortholog was performed using PAML v4.48a. Genes with ω values exceeding 1 were considered as positively selected genes (PSGs). The PSGs were further subjected to GO and KEGG enrichment analysis using Omicshare CloudTools (http://www.omicshare.com/tools/), with a significance threshold set at a q-value less than 0.05.

The complete mitogenomes of N. marisindica and S. stearnsi are 16,960 bp and 15,462 bp in length, respectively. They each contain the general 13 PCGs (COI, COII, COIII, ATP8, ATP6, ND1, ND2, ND3, ND4, ND4L, ND5, ND6 and COB), 2 rRNA genes (16S and 12S rRNA), 22 tRNA genes, and a D-loop. With the exception of the nicotinamide adenine dinucleotide hydride (NADH) genes (ND5, ND4, ND4L, ND1) and eight tRNA genes (trnF, trnQ, trnH, trnP, trnC, trnY, trnL, trnV), the majority of genes are encoded by the heavy strand (H-strand). Furthermore, apart from a few gene overlaps, the majority of nucleotide gaps are observed between the tRNA. In addition, 16S and 12S are located between trnL₁ (TAG), trnV (TAC) and trnW (TCA) respectively ( Figure 1 ; Table 1 ).

The nucleotide compositions (PCGs and rRNA) of N. marisindica and S. stearnsi are 20.9~39.8% A, 5.3~24.1% C, 32.9~48.3% T, and 3.1~24.2% G, with an AT content ranging from 62.4% to 76.8% ( Supplementary Tables 1, 2 ). Their majority of the PCGs have a negative AT skew, except for ND1, ND4, ND4L, and ND5, which display a positive GC skew ( Figure 2 ).

Except for COI, which is encoded by the ACG codon, the start codons of the PCGs predominantly follow canonical sequences (ATT/ATA/ATG). Additionally, within barnacle mitogenomes, three distinct types of stop codons are identified (T-, TA-, and TAA), with TAA being the most prevalent. Moreover, the T- and TA- stop codons only occur in COI, COIII, and ND4 of N. marisindica and S. stearnsi ( Table 1 ).

The findings show that except for leucine and serine have six to eight codons, the other amino acids of the two examined species have two to four changeable codons ( Figure 3 ). The codon usage analysis demonstrates that the following codon sequences appear with higher frequency: phenylalanine (UUU), isoleucine (AUU), and leucine (UUA). The codons with U or C in the third position occur more frequently than those with A or G. Different from the shallow-sea barnacle, the main amino acids in N. marisindica are asparagine and lysine ( Figure 3 ).

A total of 47,593,764 and 47,752,536 raw reads are generated from S. stearnsi and N. marisindica, respectively. BUSCO assessment of assembly quality are showed in Supplementary Figures 1 , 2 . The complete single-copy and complete duplicated of the two origanisms are over 80%. Through the transcriptome assemblies, there are 98,686 unigenes identified with an N50 length of 916 bp for S. stearnsi, and 53,562 unigenes with an N50 length of 1,721 bp for N. marisindica. Among these unigenes, 38,666 (39.18%) and 20,628 (38.51%) are successfully annotated in S. stearnsi and N. marisindica, respectively ( Table 2 ).

The Gene Ontology functional cluster analysis indicates that 26,718 and 16,746 unigenes are successfully clustered in the three functional categories of GO (48 functional items). The major biological processes are the cellular, metabolic, and single-organism processes, with N. marisindica lacking the biological phase item and S. stearnsi lacking the rhythmic process item. Within the cellular components, the cell and cell part items are identified. The major molecular functions are binding and catalytic activity ( Figure 4 ). In addition, The KEGG functional cluster analysis shows that 6,845 and 12,042 unigenes in N. marisindica and S. stearnsi are clustered into five functional groups, respectively, and more genes are annotated in the endocrine system, translation, and signal transduction pathways ( Figure 5 ).

From the complete mitochondrial genome of the 11 barnacle species, namely N. marisindica, Conchoderma hunteri (Owen), Lepas anserifera Linnaeus, L. australis Darwin, L. anatifera Linnaeus (Lepadidae), S. stearnsi, Arcoscalpellum epeeum Chan (Scalpellidae), Altiverruca navicula Hoek (Verrucidae), Capitulum mitella Linnaeus (Pollicipedidae), Glyptelasma annandalei Pilsbry (Poecilasmatidae), and Ibla cumingi Darwin (Iblidae) ( Table 3 ), 13 PCGs were extracted for positive selection pressure analysis. The positive selection analysis yielded results indicating that the COIII and ND2 genes (ω > 1.0, P<0.05) of deep-sea species experience positive selection pressure. The three-dimensional structure model was utilized to identify and highlight specific sites within the genes that are under positive selection pressure, such as 22I, 29T, and 41E in the COIII gene, as well as 73F, 75S, and 76L in the ND2 gene ( Table 4 ; Figure 6 ).

Genes that are subjected to positive selection tend to respond to natural selection. We identified 148 pairwise orthologs within N. marisindica, S. stearnsi and two reference species in one-to-one correspondence to confirm the PSGs ( Supplementary Table 3 , Supplementary Figure 3 ) and performed an adaptive evolutionary analysis on them. The findings reveal the presence of five orthologs, namely chromosome 3 open reading frame 10 (C3ORF10), E74-like factor (ELF), ubiquitin-like protein ATG12 (ATG12), uracil phosphoribosyltransferase (UPRT), and small subunit ribosomal protein S21e (RPS21), all exhibiting ω values greater than 1.0. Notably, only RPS21 did not fall below the significance threshold of 0.05 ( Table 5 ; Supplementary Figure 3 ). These PSGs are mainly enriched in the RIG-I-like receptor signaling, regulation of autophagy, apoptosis-fly, regulation of the actin cytoskeleton, pyrimidine metabolism, foxO signaling, and ribosome pathways ( Figure 7 ; Table 5 ). It is predicted that the ORFs of the five genes under positive selection pressure are all complete. Except for C3ORF10, the amino acid domain positions of the other genes are predicted, there is no signal peptide, and the molecular weight range from 9.21 to 27.23 kDa. Additionally, the isoelectric point is relatively similar, ranging from 6.05 to 9.88 ( Supplementary Table 4 ).

Previous research has identified that certain barnacles exhibit similar arrangement structures within their mitogenomes (Shen et al., 2015). We also find similar genetic structures in the mitogenomes of the deep-sea and shallow barnacles. The AT content in the complete mitogenomes of both organisms, exceeded the GC content, aligning with the base composition observed in the majority of invertebrates (Liu et al., 2018; Wang et al., 2021; Sato et al., 2023). More negative AT skews are observed in ND1, ND4, ND4L, and ND5 genes of the two barnacles ( Figure 2 ). This phenomenon might be attributed to the deep-sea habitat of these organisms, where DNA replication and transcription processes are subject to the mutations and the natural selection resulting in a bias towards specific nucleotides (Sahyoun et al., 2014; Yu et al., 2019). It is speculated that the different genetic structures of these four genes might be closely related to the pressure of natural selection in the extreme habitats. However, how nucleotide bias affects the initiation and orientation of gene replication and what selective pressures induced nucleotide skews need further study. Furthermore, the amino acid usage patterns observed in N. marisindica are similar to those of the shallow barnacle S. stearnsi ( Figure 3 ), which provides more evidence for the barnacle nucleotide composition and gene skew results. Generally, more than 95% of energy comes from the oxidative phosphorylation of mitochondria (Auger et al., 2021). In the abyssal zone, the dark and absence of autotrophic primary producers pose constraints on the availability of sustenance (Gattuso et al., 2006). The observed nucleotide compositional biases between N. marisindica and S. stearnsi indicated that, in their mitogenomes, comparable alterations occurred during the evolutionary progression which might be a means of adapting to the extreme conditions and nutritional scarcity of the deep-sea environment.

Mitochondria serve as a vital pathway for aerobic respiration within cells. To cope with the extreme environment, selection pressure would have been put on the mitogenome, and non-synonymous mutations are regarded as an adaptive substitution to meet the needs of a new environment (Srivastava et al., 2020). The results of our study indicate that the mitogenomes in N. marisindica experienced significant selective pressure, particularly the COIII and ND2 genes, which have specific functional protein sites ( Table 4 ; Figure 6 ). There is a prevalent belief that functional proteins contain positive selection sites, which serve as a foundation for adaptive responses to alterations in the environment (Rausher and Huang, 2016). COIII is one of the core members of cytochrome C oxidase, a key rate-limiting enzyme at the end of the mitochondrial respiratory chain, which produces ATP and changes oxygen into water for cell metabolism (Timón-Gómez et al., 2018). At the level of mitochondrial transcription and translation, during the respiratory chain embedding in the inner mitochondria, the respiratory supercomplex factor protein family acts in response to biological hypoxia by binding to COIII to adapt to environmental and metabolic changes (Lobo-Jarne and Ugalde, 2018). Through bioinformatics analysis, several adaptive residues potentially presented in COIII and ND2 of alvinocaridid shrimp (Alvinocarididae) inhabiting the deep-sea hydrothermal vent environment also have been identified by researchers. These proteins played an important role in regulating mitogenome intracellular energy metabolism and contributed to the adaptation of organisms to the hydrothermal vent environment (Hui et al., 2018b). Consequently, we speculate that the PCGs help N. marisindica to resist multiple pressures when adapting to the deep-sea environment.

The NADH dehydrogenases are the largest enzyme complex in the mitochondrial electron transport chain; thus, mutations in these subunits may disturb the efficiency of the electron transport process (Da Fonseca et al., 2008). The NADH dehydrogenase genes of deep-sea anemones exhibited novel genetic structures, which were evidence of adaptation to the deep-sea environment (Zhang et al., 2017). Likewise, the hydrothermal vent environment may potentially exert potential selection pressure on NADH dehydrogenase genes from mitogenomes in N. marisindica. This pressure may influence the modulation of mitochondrial complexes and electron transport efficiency and hence, adapting to cold and hypoxic conditions (Luo et al., 2012). Given the crucial role of mitochondria in aerobic respiration and their indispensability for metabolism and energy production, they likely aid the deep-water barnacle in surmounting the constraints imposed by hypoxia, cold and food availability. Additionally, hydrothermal vents are also associated with hydrogen sulfide (H₂S), and the primary cytotoxic effect of H₂S is to interrupt the mitochondrial respiratory chain by directly inhibiting cytochrome C oxidase (Cooper and Brown, 2008). Therefore, we infer that the mitogenomes of N. marisindica are subjected to multiple selection pressures which led to adaptation to the hydrothermal vent environment under extreme deep-sea conditions.

Since they cannot swim in the ocean, barnacles are constrained to inhabit a comparatively restricted narrow region, rendering the process of acclimating to the deep-sea milieu arduous. The Southwest Indian Ridge (SWIR) is a typical ultraslow-spreading oceanic ridge, and the Longqi hydrothermal field is the first confirmed active field on the SWIR. Due to the general lack of deep-sea sediments in this hydrothermal zone, which can provide good channels for hydrothermal circulation and spillage, there are a large number of sulfides (Liang et al., 2023). The extreme environment impacts on population vicariance and dispersal, both significant biological processes. However, few researches have focused on the biodiversity of this region (Zhou et al., 2018). Therefore, our research endeavors to fill this gap by examining the effects of N. marisindica and S. stearnsi in the marine habitat, thereby enhancing our understanding of the adaptations in the harsh deep-sea environment.

ATG12, a core autophagy protein. The interaction between ATG12 and its active complex can form autophagosome. Autophagy is an evolutionarily conserved catabolic process that generally promotes cell survival by eliminating intracellular pathogens and removing damaged organelles, playing a pivotal role in cellular homeostasis and adaptation to adverse environments (Chen et al., 2014). During autophagy, ATG12 binds the E2 ubiquitin conjugating enzyme to form sulfolipid intermediate, which further mediate the binding of phosphatidylethanolamine to LC3 to promote autophagy formation (Yuan et al., 2022). Although the regulation of this process remains incompletely understood, it has been found to induce the formation of the autophagosome membrane and type I interferons (IFNs). The process have been shown to play a crucial role in the antiviral immunity of Larimichthys crocea (Richardson) (Wei et al., 2021), Epinephelus coioides (Hamilton) (Li et al., 2019), and Ctenopharyngodon idella (Valenciennes in Cuvier and Valenciennes) (Chu et al., 2019), shedding light on the potential involvement of ATG12 in antiviral immunity by stimulating both autophagy and type I interferon response. Consequently, it is hypothesized that ATG12 may eliminate intracellular pathogens and damaged celles within N. marisindica by regulating autophagy homeostasis, which have contributed to its adaptation of harsh deep-sea environments. Furthermore, previous studies have demonstrated the involvement of ATG12 in various cell biological processes such as viability, apoptosis, and autophagy (Wei et al., 2019). Meanwhile, ATG12 has been found to directly modulate the mitochondrial apoptotic pathway by binding and inactivating prosurvival Bcl-2 family members (Rubinstein et al., 2011). In contrast to barnacles residing in shallow waters, the survival of N. marisindica exposed to diverse toxic substances, including H₂S. This exposure may have potentially facilitated the evolution of more effective immune responses and antiapoptotic metabolism (Cooper and Brown, 2008; Liang et al., 2023).

In Addition, the ELF (OG15174) and UPRT (OG15292) genes are also positively selected from the enriched pathways of apoptosis-fly and pyrimidine metabolism in N. marisindica, and they are both involved in adaptive immune responses. ELF, a pleiotropic transcription factor that modulates disease susceptibility and mediates several biotic and abiotic stress processes, playing a significant role in the host innate immune defense system (Suico et al., 2017). pH, temperature, and salinity can affect the immune response of crustaceans. ELF was involved in regulatory processes in organisms through synergistic or retrograde interactions, maintaining osmotic pressure homeostasis (Zhang et al., 2019). Previous researches has demonstrated that ELF was significantly expressed in the hepatopancreas and gill tissue of Penaeus monodon Fabricius under salt stress, showing that ELF was involved in the immune stress process of organisms by regulating immune responses and the development of immune-related cells (Si et al., 2022). Simultaneously, UPRT maintains the immunocompetent and regulates the sensitivity of organisms to harmful cells (Koyama et al., 2000; Hasegawa et al., 2013). UPRT is a pyrimidine salvage enzyme. It can convert 5-flurouracil to a toxic compound, which kills the deleterious cells but subsequent conversion to other metabolites enhances its anti-apoptotic effect, by irreversibly inhibiting thymidylate synthase, DNA synthesis and the processing of rRNA and mRNA (Gopinath and Ghosh, 2008). Thus, it is plausible that these PSGs can potentially contribute to maintaining the immune system of N. marisindica, enabling it to withstand the adverse conditions of the deep sea, including low temperatures and toxic substances. Furthermore, high hydrostatic pressure and low temperatures have caused the inhibition of ligand binding and ion channel function which in turn affected signal transduction (Shinozaki and Yamaguchi-Shinozaki, 2000; Siebenaller, 2000). C3ORF10 has been proposed to regulate actin dynamics, and its expression played a role in regulating actin dynamics during the process of neurite outgrowth (Stradal and Scita, 2006). C3ORF10 is a member of the SCAR/WAVE complex, which is the effector of Rac1 to regulate actin dynamics, regulating Arp2/3-dependent actin polymerization and the number of axonal collaterals and dendrites (Wang et al., 2013). Consequently, the genes that have undergone positive selection and involved energy utilization, immune response, and signal transduction metabolism, might enhance the ability of N. marisindica to effectively employ various diverse adaptive strategies within the extreme conditions of the deep-sea environment. N. marisindica also exhibits enrichment of several significant pathways involving genes subject to positive selection, including the apoptosis-fly, regulation of actin cytoskeleton, pyrimidine metabolism, and ribosome pathways. Functionally, the ribosome pathway (K02971) influences almost every aspect of cell and organismal biology including maintaining cell homeostasis in response to changing environmental conditions and acute stress (Ni and Buszczak, 2023). Moreover, impairment of translation can result in the activation of the ribotoxic stress response, implying that ribosomes serve as a signal for the metabolic regulation of stress response (Snieckute et al., 2022), and the inhibition of ribosomal activity hinders the progression of the cell cycle, thereby regulating the adaptive stress response (Stoneley et al., 2022). Consequently, the presence of genes enriched in the ribosomal pathway indicates that the organisms have experienced environmental stress in their habitat.