From head to rootlet: comparative transcriptomic analysis of a rhizocephalan barnacle Peltogaster reticulata (Crustacea: Rhizocephala)

Sampling

Hermit crabs Pagurus minutus infected with Peltogaster reticulata were collected in the Sea of Japan (Marine Biological Station “Vostok” of the Institute of Marine Biology of the Russian Academy of Sciences) (N: 42.893720, E: 132.732755). All the parasites were adults with fully developed externae.

The infected crabs were dissected in filtered sea water. The parasite was removed from the host’s body cavity and the remains of the host tissues were isolated from the interna. The body of each parasite was divided into four parts: 1) the externa (it was separated at the level of the stalk), 2) distal part of the main trunk (the growing part) (approximately the last three millimeters of the main where muscular system is not developed yet²⁷), 3) the main part of the interna located in the abdomen of the host (the main trunk), and 4) the part of interna from the thorax of the host (the thoracic part) (Figure 1). For each body part the pooled sample was prepared, containing material from five parasitic individuals in two biological replicates. The samples were collected into centrifuge tubes and frozen at -80 °C in IntactRNA (Evrogen, Moscow, Russia) reagent according to the manufacturer`s protocol.

Figure 1. Molecular signatures of female Peltogaster reticulata.

(a) Generalized scheme of female P. reticulata in the host. Colour sectors indicate the body parts examined in our study: externa (red), growing part of interna (green), main trunk of interna (purple), and thoracic part (blue). (b) The number of common OMA groups. The colour key on the heatmap shows the number of shared OMA groups between species. (c, e) Venn diagram for a set of genes either included in molecular signatures of the body parts (c) or over-expressed in the body parts (e). (d, f) Multidimensional scaling (MDS) plots for molecular signatures (d) and sets of over-expressed genes (f). Different clusters in MDS plots are marked with colours. Abbreviations: MS/overexp – the molecular signature or set of over-expressed genes for the body part, respectively; rep1/2 – biological replication identifier.

Before RNA isolation, the IntactRNA-fixed samples were rinsed in 0.1M phosphate-buffered saline (PBS). The total RNA was isolated using Quick-RNA MiniPrep (R1054, Zymo Research, Irvine, California, USA) according to the manufacturer’s protocol. The libraries were synthesized using NEBNext Ultra Directional RNA Library Prep Kit for Illumina (E7760, New England BioLabs, Ipswich, Massachusetts, USA) according to the manufacturer protocol. Paired-end sequencing was carried out using Illumina HiSeq 2500 instrument (Illumina, San-Diego, California, USA).

Sampling was conducted in accordance with the European Community Council Directive of November 24, 1986 (86/609/EEC). All possible efforts were made to minimize the number of animals used.

Preparation of reads libraries and de novo transcriptome assembly

The primary quality control of paired-end reads libraries was manually assessed using FastQC (v0.11.5) [https://www.bioinformatics.babraham.ac.uk/projects/fastqc/]. The potential sequencing error identification and correction was performed by Karect²⁸ (v1.0) [https://github.com/aminallam/karect] with the following parameters: --celltype=diploid –matchtype=hamming. Trimmomatic (v0.39) [http://www.usadellab.org/cms/?page=trimmomatic] was used for removing the sequencing adaptors (ILLUMINACLIP:$ADAPTERS:2:30:10:2:TRUE), low-quality read regions (SLIDINGWINDOW:4:20 MAXINFO:50:0.8) as well as reads with a length less than 25 nucleotides (MINLEN:25). Since library preparation and sequencing were performed in the laboratory where researchers also work with human medical samples, the parasite data were checked for the presence of the read pairs with a high identity to the Homo sapiens reference transcriptome (GENCODE v.31) using BBTools (v37.02).

The prepared libraries were pooled and used for de novo reference transcriptome assembly using Trinity²⁹ (v2.5.1) [https://github.com/trinityrnaseq/trinityrnaseq] with k-mer size and required minimal contig length equal to 25 and 200 nucleotides, respectively. The assembled contigs were renamed by adding the four-digit tag of the species, “Pret”, at the beginning of IDs. Isoforms were clustered on all assembled contigs using CDHIT-est³⁰ (v4.7) [http://weizhong-lab.ucsd.edu/cd-hit/] and a sequence identity threshold equal to 95% (-c 0.95), accurate mode (-g 1), and both +/+ and +/- strands alignments (-r 1). TransRate³¹ (v1.0.1) [https://hibberdlab.com/transrate/] was used for the quality assessment of clustered sequences. Only the contigs classified as “good” by TransRate³¹ were included in further analysis.

Removal of potential contamination

The 18S and 28S ribosomal RNA sequences were searched using RNAmmer³² (v1.2) [https://services.healthtech.dtu.dk/service.php?RNAmmer-1.2] from the Trinotate pipeline (v3.1.1) [https://anaconda.org/bioconda/trinotate]. The sequences obtained this way were compared with the NCBI nucleotide database with BLASTn³³ (v2.6.0+) [https://blast.ncbi.nlm.nih.gov/Blast.cgi?PAGE_TYPE=BlastDocs&DOC_TYPE=Download] to identify their possible sources.

We used MCSC (Model-based Categorical Sequence Clustering) Decontamination method³⁴ [https://github.com/Lafond-LapalmeJ/MCSC_Decontamination] for removing potential contamination, with “Arthropoda” as the target taxon and the clustering level equal to 5 (32 clusters). The parsing of BAM-files with reads alignment results created by Bowtie2³⁵ [http://bowtie-bio.sourceforge.net/bowtie2/index.shtml] was performed to extract only reads pairs that mapped to the decontaminated set of contigs.

Quantification of gene expression levels and identification of encoded amino acid sequences

Salmon³⁶ (v1.0.1) was used for expression level quantification (-l ISF –discardOrphans –seqBias –gcBias –validateMappings). The expression quantification results and the transcripts-to-genes map from the Trinity output were provided to the “tximport” package for R to obtain expression levels of genes. The tables with both unaveraged transcripts-per-million (TPM)³⁷ values and TPM values averaged between biological replicates were prepared. Only the sequences with expression levels ≥ 1 TPM in at least one sample were included in further analysis.

TransDecoder (v5.5.0) [https://github.com/TransDecoder/TransDecoder] was used for the identification of the amino acid sequences encoded by assembled contigs. Firstly, the long open reading frames with a length ≥ 100 amino acids (aa) and products of its translation were found. Secondly, identified proteins were compared with the NCBI non-redundant (DIAMOND BLASTp³⁸ (v0.9.22.123) [https://github.com/bbuchfink/diamond/], e-value = 1e-3) and the PfamA³⁹ [http://pfam.xfam.org] (HMMscan (v3.1b2)) databases. Thirdly, the comparison results were provided to the TransDecoder to identify the likely coding regions and to obtain the probable set of proteins.

Reference gene set preparation

In our analysis, the focus was only on genes that successfully passed two filters: 1) noticeable expression level (i.e., the expression is ≥ 1 TPM in at least one sample) and 2) encoding of the proteins with a length greater than or equal to 100 amino acids. Only the longest protein and its coding transcript were selected as representatives for each gene and referred to as “reference sets”. The completeness of the protein reference set was evaluated by comparison with the database of single-copy orthologues of Metazoa (odb-9) using BUSCO^40,41 (v3.0.1) [https://gitlab.com/ezlab/busco] (e-values = 1e-3, mode = proteins).

Sequence annotation

For the annotation of the genes, their nucleotide and amino acid sequences were compared with publicly available databases: NCBI nucleotide collection (nt), NCBI non-redundant (nr), and SwissProt^42,43. The similarity search was carried out with BLASTn megablast³³ (nt) and the sensitive mode of DIAMOND BLASTp³⁸ (amino acid databases), with an expected value (e-value) threshold equal to 1e-3 and a limit up to 10000 for the number of description and alignments⁴⁴. The best BLAST hits (BBH) were selected with a custom script.

The potential domain architecture of the proteins was identified using the PfamA database (HMMScan) and custom script. The proteins were also analysed using the eggNOG-mapper web-resource⁴⁵ (v2) [http://eggnog-mapper.embl.de] with default parameters.

Identification of orthogroups

Orthogroups were identified with the use of OMA standalone program (v2.5.0)⁴⁶ [https://omabrowser.org/standalone/] in three steps. Firstly, the reference proteomes of Amphibalanus amphitrite Darwin, 1854 (UP000440578), Armadillidium nasatum Budde-Lund, 1885 (UP000326759), Armadillidium vulgare Latreille, 1804 (UP000288706), Daphnia magna Straus, 1820 (Strain: Xinb3) (UP000076858), Daphnia pulex Leydig, 1860 (UP000000305), Penaeus vannamei Boone, 1931 (UP000283509), Portunus trituberculatus Miers, 1876 (UP000324222), and Tigriopus californicus Baker, 1912 (UP000318571) were downloaded from the UniProt^42,43 [https://www.uniprot.org] database [accessed 21 October 2021]. Only the sequences with a length equal to or more than 100 amino acids were analysed. The OMA standalone was run with default parameters with the “bottom-up” algorithm for inference of hierarchical orthologous groups (HOGs), without a phylogenetic tree, but with the identification of two Daphnia species as an out-group. Secondly, we reconstructed the phylogenetic tree following the protocol by Dylus et al.⁴⁷. Briefly, using the filter_groups.py provided, we selected OMA groups that included at least eight of the nine crustacean species involved in the analysis. Then, using MAFFT⁴⁸ (v7.487) [https://mafft.cbrc.jp/alignment/software/], multiple protein alignment in each orthogroup was performed (--maxiterate 1000 -localpair). The alignments were concatenated into a supermatrix using the concat_alignments.py. The selection of suitable sites in the supermatrix was carried out using tramAl⁴⁹ (-automated1) [http://trimal.cgenomics.org]. We used the ProtTest program^50,51 (v3.4.2) [https://github.com/ddarriba/prottest3] to determine the most appropriate sequence evolution model. The phylogenetic tree was reconstructed using the IQ-TREE^52,53 (v2.1.4-beta) [http://www.iqtree.org] with the following parameters: -m LG+I+G+F --seed 12345 -B 1000 --nmax 1000. The consensus tree was rooted by the out-group using the “ape”⁵⁴ [https://cran.r-project.org/web/packages/ape/index.html] library for R. Thirdly, the phylogenetic tree was used when the OMA standalone⁴⁶ was re-run with default settings. The construction of a heatmap with the number of common OMA groups between the studied species was performed in RStudio using the “ggplot2” (v3.3.5), “pheatmap” (v1.0.12), and “RColorBrewer” (v1.1-2) libraries.

Reference gene set expression analysis

The “molecular signature” of a body part was defined as a set of genes with an expression level ≥ 2 TPM in the body part. The expression threshold value was chosen in accordance with the results of studies by Wagner, Kin, and Lynch, according to which “genes with more than two transcripts per million transcripts (TPM) are highly likely from actively transcribed genes”⁵⁵. Genes that had an expression ≥ 2 TPM in all the body parts, were classified as “commonly expressed”.

Significant variation of gene expression between samples was detected using “RNentropy” library⁵⁶ (v1.2.2) [https://cran.r-project.org/web/packages/RNentropy/index.html] for R. The analysis was carried out using a table with unaveraged TPM values between replicates. The corrected global sample specificity test P < 0.01 was used according to the Benjamini-Hochberg method, and local sample specificity test P < 0.01.

The overlaps between the molecular signatures and the sets of over-expressed genes were visualized with InteractiVenn⁵⁷ [http://www.interactivenn.net].

Multidimensional scaling

A multidimensional scaling (MDS) analysis of the molecular signatures and the sets of over-expressed genes was performed. The presence / absence matrices were used as input. It was indicated in the matrix for each gene (row) whether the gene was included in the molecular signature of the body part or had an increased expression in it (“1”) or not (“0”). The optimal number of clusters was determined using the “silhouette” method implemented in the “factoextra” library (v1.0.7) for R. The metaMDS function from the “vegan” library (v2.5-7) was used with the following parameters: distance = "manhattan", try = 100, trymax = 100000, autotransform = FALSE, binary = TRUE, k = the optimal number of clusters. The seed was set to 1234 both when the optimal number of clusters was determined and in MDS. To visualize the results, the ggscatter function from “ggpubr” (v0.4.0) library for R was used.

Potential excretory/secretory proteins (ESP) identification and analysis

The in silico identification of potential ESP was performed according to the pipelines described by Garg and Ranganathan⁵⁸. Firstly, all proteins from the reference set were analysed with SignalP⁵⁹ (v5.0b) [https://services.healthtech.dtu.dk/service.php?SignalP-5.0]. Based on the analysis results, the proteins were divided into potential “classical” (SP ≥ 0.5) and “non-classical” (SP < 0.5) ESP. Secondly, the “non-classical” ESP were analysed using SecretomeP⁶⁰ (v1.0) [https://services.healthtech.dtu.dk/service.php?SecretomeP-2.0]. Only proteins that had NN-scores ≥ 0.9 and were simultaneously predicted not to contain a signal peptide were selected. Thirdly, all potential ESP were scanned for the presence of the mitochondrial transit peptide with TargetP⁶¹ (v2.0). The proteins with this signal were excluded. Fourthly, the transmembrane hidden Markov model (TMHMM)⁶² (v2.0c) [https://services.healthtech.dtu.dk/service.php?TMHMM-2.0] was used to model and predict the location and orientation of transmembrane domains in proteins, and only proteins without them were considered as potential ESP. Out of these potential ESP, only those were selected that were included in at least one molecular signature. The overlap analyses between ESP sets were performed using InteractiVenn⁵⁷.

We ran all against all BLASTp searching using DIAMOND³⁸ (--evalue 1e-3) for both classical and non-classical ESP. Similarity search results were used in SiLix⁶³ (v1.2.11) [https://lbbe-web.univ-lyon1.fr/fr/SiLix] (-r 0.9), which assigns proteins to putative gene families. For annotation, all ESP were compared against the NeuroPep⁶⁴ [http://isyslab.info/NeuroPep/] and MetazSecKB⁶⁵ [http://proteomics.ysu.edu/secretomes/animal/] databases using DIAMOND BLASTp³⁸ with the following parameters: --sensitive --max-target-seqs 10000 –evalue 1e-3. The best BLAST hits were selected using a custom script.

Gene set enrichment analysis (GSEA)

The GSEA using “topGO” library (v2.40.0) [https://bioconductor.org/packages/release/bioc/html/topGO.html] for R was performed for 1) whole molecular signatures, 2) sets of over-expressed genes, and 3) sets of potential “classical” and “non-classical” ESP. Only the Gene Ontology (GO) terms describing biological processes were considered. Fisher’s exact test was used and extracted only the terms including at least 10 significant genes (GO terms with p-value <0.01) from the results. Redundancy was reduced with the “rrvgo” library (v1.0.2) [http://bioconductor.org/packages/release/bioc/html/rrvgo.html] for R. The minus log10-transformed p-values were used as scores, org.Dm.eg.db (Genome wide annotation for Fly) as database, relevance as similarity measures methods, and 0.7 as the threshold for reduceSimMatrix. The “wordcloud” (v2.6) library [https://cran.r-project.org/web/packages/wordcloud/index.html] for R was used to build word clouds based on the redundancy reduction results. The more often the parental bioprocess was found in the list of enriched bioprocesses, the larger the word size. Each bioprocess was assigned a colour from the “viridis” (v0.6.1) palette package.

Phylostratigraphy and Transcriptome Age Index (TAI) measuring

The phylostratigraphic analysis of the P. reticulata reference protein set was performed using the “phylostratr” package⁶⁶ (v0.2.1) [https://github.com/arendsee/phylostratr] for R. We used reference proteomes of Crustacea prepared beforehand as well as the prebuilt dataset of prokaryotes, human, and yeast. In total, in addition to the P. reticulata, 139 species were included in the analysis, of which 44 were representatives of the Metazoa. The complete phylogenetic tree of the studied species is presented in the Figure S1 (Extended data⁶⁷). Similarity search between proteins was carried out with BLASTp³³ (v2.6.0+). The tables with BLAST results in “6” output format were used as input for “phylostratr” for protein distribution between phylostrata: 1) “Cellular organisms”, 2) “Eukaryota”, 3) “Opisthokonta”, 4) “Metazoa”, 5) “Eumetazoa”, 6) “Bilateria”, 7) “Protostomia”, 8) “Ecdysozoa”, 9) “Panarthropoda”, 10) “Arthropoda”, 11) “Mandibulata”, 12) “Pancrustacea”, 13) “Crustacea”, 14) “Multicrustacea”, 15) “Hexanauplia”, 16) “Cirripedia”, and 17) “Peltogaster reticulata”.

The phylostratigraphic composition was analysed for the P. reticulata reference gene set, the set of genes with noticeable (≥ 2 TPM) expression in all the female body parts considered, the sets of overexpressed genes as well as the sets of genes encoding potential “classical” and “non-classical” ESP. The results were visualized using “ggplot2”, “viridis” (v0.6.1), and “reshape” (v0.8.8) libraries for R.

Transcriptome Age Index (TAI) definition was performed for P. reticulata body parts using phylostratigraphic results and tables with averaged TPM-values. The analysis was carried out using the “myTAI”⁶⁸ (v0.9.3) [https://github.com/drostlab/myTAI] package for R. Genes with an expression level < 2 TPM at all body parts were excluded. The analysis was carried out on log2(TPM + 1) transformed values. The FlatLineTest function was used to quantify the statistical significance of the global TAI pattern. For analysis with the use PlotRE and PlotBarRE functions, the phylostrata were divided into two groups: “before” (phylostrata 1–13), and “after” (phylostrata 14–17) the division of Crustacea.

Using the pMatrix function from “myTAI”⁶⁸, the contributions of genes to the TAI of body parts were determined. For each body part, 500 genes with the largest contribution were selected out of the genes with the GO annotation. Further, GSEA for the selected gene sets was performed similarly to GSEA for the molecular signature.

Histology

For histological and light-microscopic examination, the dissected internae were fixed with Bouin solution (picric acid (trinitrophenol) 71.5%, paraformaldehyde 24% and acetic acid 4.5%). Paraffin sections (5 μm thick) were made using standard histological methods with the help of a Leica RM-2265 microtome and stained with hematoxylin-eosin. The sections were examined under a Leica DM2500 microscope. The photos were taken with a Nikon DS-Fi1 camera and processed with ImageJ software (FiJi⁶⁹).

Confocal laser scanning microscopy (CLSM)

Samples of interna for immune labelling were fixed with 4% paraformaldehyde (PFA; Sigma-Aldrich) in PBS (Fluka) at 4 °C for four hours, and then rinsed three times with PBS. Prior to immunocytochemical staining, the fixed material was incubated with PBST (PBS + 0.1 % Triton-X100; Sigma-Aldrich) during 24 hours at 4 °C. Then, the samples were incubated in primary antibodies and anti-acetylated α-tubulin (Sigma Aldrich, Germany, T6793, produced in mice) and anti-serotonin (Sigma Aldrich, Germany, S5545, produced in rabbit) for three days. After incubation the specimens were rinsed in PBS three times and incubated in secondary antibodies anti-mouse IgG CFTM 633 (Sigma Aldrich, Germany, SAB4600138) and anti-rabbit IgG CFTM 488A (Sigma Aldrich, Germany, SAB4600030).

The specimens were rinsed with PBS three times and stained with the DAPI nuclei stain (1 μg/ml; Sigma-Aldrich) for 30 min, rinsed in PBS and mounted in DABCO-glycerol. The samples were examined using a Leica TCS SP5 confocal laser scanning microscope in the Resource Center “Microscopy and Microanalysis” of the Research Park of Saint Petersburg State University. The images were processed with ImageJ software (FiJi⁶⁹).

Scanning electron microscope (SEM)

Specimens for SEM were fixed at 4 °C in 2.5% glutaraldehyde, dehydrated in a gradient ethanol series and acetone, critical point-dried in a Hitachi critical point dryer HCP- 2, mounted on stubs, coated with platinum with the use of a Giko IB-5 Ion sputter coater, and viewed under a FEI Quanta 250 scanning electron microscope in the “Taxon” Research Resource Center of the Zoological Institute of the Russian Academy of Sciences.

The de novo assembled transcriptome was characterized by high quality and completeness of assembly

The transcriptomes of the whole body, the thoracic part of the interna, the growing part of the interna and the main trunk of the interna as well as that of the externa of adult P. reticulata were collected and sequenced in two biological replicates (Figure 1a). More than 87% of read pairs remained in all the paired-end libraries after the removal of adapters, poor-quality regions, and short sequences (Table S1, Underlying data⁷⁰). The potential contamination with human-derived sequences did not exceed 4.1% from the total number of trimmed read pairs in each library (Table S1, Underlying data⁷⁰).

All prepared libraries were merged, and the resulting libraries were used as input for Trinity software. In general, 353,130 contigs with lengths greater than or equal to 200 nucleotides were assembled de novo. After the clusterization of similar sequences, the P. reticulata transcriptome included 267,188 contigs. TransRate assembly and optimal scores made up 0.3331 and 0.3835, respectively. More than 95% (256047/267188) of the contigs were well-assembled (“good”) according to the TransRate quality control results. The completeness analysis for “good” contigs using a database of the metazoan single-copy orthologues revealed that 96.3% (Single: 57.8%, duplicated: 38.5%) of the orthologues were assembled completely.

Peltogaster reticulata reference transcriptome included transcripts of 12,620 protein-encoding genes

Given the parasitic lifestyle of P. reticulata, the assembled transcriptome was checked for the presence of potential contamination. According to the RNAmmer analysis results, eight and nine sequences could be classified as 18S and 28S ribosomal RNAs, respectively. The comparison with the NCBI nucleotide database revealed that the contigs aligned successfully with ribosomal sequences from Alveolata (HQ891115.2), Fungi (AY382649.1; CP033152.1; CP030254.1; GQ336996.1; MF611880.1), and Metazoa (AY265359.1; EU082415.1; KY454201.1; EU370441.1; KU052603.1). Among the metazoan hits, only sequences from Crustacea were identified. The identity percentage with the database-derived 18S host (KY454201.1) and parasite (EU082415.1) sequences made up 88% (990/1,123) and 99% (1,200 / 1,201), respectively. In this study, the MCSC hierarchical clustering algorithm was used to remove potential contamination from the assembled transcriptome. The decontaminated transcriptome included 80,779 contigs and contained 81.6% (Single: 56.9%; Duplicated: 24.7%) of completely assembled and 6.2% fragmented single-copy metazoan orthologues. The number of paired-end reads successfully aligned to selected contigs varied from 10.74 (externa, first replicate) to 30.64 million (the thoracic part of the interna, second replicate) (Table S1, Underlying data⁷⁰).

The sequence expression level quantification was performed with Salmon by mapping selected read pairs to the decontaminated transcriptome. The mapping rate ranged from 85% to 91%. The transcript-to-gene map was used to obtain gene expression levels in TPM values. After excluding genes with a low activity in all analysed samples (expression level < 1 TPM), the dataset contained TPM values for 20,980 contigs. According to the TransDecoder results, 32,990 contigs encoded proteins with lengths ≥ 100 amino acids.

Only the protein-coding genes with a noticeable expression were involved in further analysis. After filtering by expression level (≥ 1 TPM in at least one sample) and encoded protein lengths (≥ 100 aa), the reference gene set obtained for P. reticulata contained 12,620 sequences. For each gene, the longest protein encoded by its splice variants was selected as a representative sequence. The comparison with a single-copy metazoan orthologues database revealed that in the reference protein set of P. reticulata 75.6% of the orthologues were assembled completely (Single: 70.4%; Duplicated: 5.2%), whereas 3.4% and 21% of the sequences were fragmented or absent, respectively.

Most of the sequences from the reference sets were annotated successfully

The prepared reference sets were compared with publicly available databases. According to the results, 3,399, 8,502, and 9,690 genes had hits with NCBI nucleotide, SwissProt, and NCBI non-redundant databases, respectively. The overlap analysis revealed that 3,240 genes were successfully annotated using each database. Moreover, 6,090 genes belonged to at least one GO terms. The domain architecture of the encoded protein was identified for 8,803 genes based on the comparison with the PfamA database. Annotation results are presented in Table S2 (Underlying data⁷¹).

The proteins set of P. reticulata was similar to the reference proteome of another cirripede barnacle, Amphibalanus amphitrite.

The identification of orthogroups in P. reticulata and other Crustacea involved in this study was carried out in three stages. First, the proteomes of P. reticulata and eight reference crustacean species were analysed using OMA standalone. As a result, 24,840 OMA groups were discovered.

Secondly, a phylogenetic tree of the studied crustacean species was constructed based on the results (Figure S2, Extended data⁷²). A total of 609 OMA groups were selected, containing at least eight out of nine species. The multiple protein alignment results in each of the OMA groups were concatenated into a supermatrix. After site selection, the final supermatrix contained 282,907 sites. In the resulting phylogenetic tree, P. reticulata was united into the same taxon with another cirripede barnacle, Amphibalanus amphitrite, with full support (Figure S2, Extended data⁷²).

Thirdly, the resulting tree was used to refine the search results for orthogroups. More than 20,000 orthogroups were found: 24,840 and 20,874 OMA and HOGs were identified, respectively. Figure 1b shows the number of common OMA groups for pairs of species. P. reticulata had the largest number of "common" OMA Groups (4,354) with A. amphitrite. For comparison, P. reticulata had no more than 3,490 “common” OMA groups with other crustacean species, the smallest overlap being found with Portunus trituberculatus (1,262 OMA groups).

The externa was clustered separately from the interna based on the gene expression analysis results

The molecular signature of a body part was defined as a set of genes with an expression ≥ 2 TPM in the body part considered. Each molecular signature included at least 8,000 genes: the main trunk of interna (8,070 genes), the growing part (8,148), the thoracic part (8,223), and externa (9,233) (Table S3, Underlying data⁷³). Approximately 54% (6,829 / 12,620) of the genes from the reference set were included into molecular signatures of all body parts (Figure 1c). Figure 1c shows significant overlaps between the parts of the interna and an almost 10-fold difference in the number of “specific” genes between the interna parts and the externa. Based on the gene expression, the body parts were divided into two clusters (Figure 1d). The first cluster contained all parts of the interna, while the second cluster contained only two replicates of the externa.

According to the differential expression analysis results, the number of over-expressed genes varied from 204 (the main trunk of the interna) to 2,224 (the externa) (Table S3, Underlying data⁷³). The number of genes with an increased expression in the interna did not exceed 283. Figure 1e shows that (i) only one gene was over-expressed both in the externa and in the interna, (ii) the number of “shared” over-expressed genes between different interna parts was ≤ 35, (iii) only six genes were over-expressed in the whole interna. The externa clustered separately from the interna, and all interna parts were remote from each other (Figure 1f).

The results of identification of molecular signatures and differentially expressed genes are presented in Table S3 (Underlying data⁷³).

GSEA results and histological studies indicate the presence of germ-like cells in the interna

Figure 2a and b show Venn diagrams for lists of bioprocesses enriched with genes included in the molecular signatures (Figure 2a) and genes with increased expression (Figure 2b) in the female body parts considered. In contrast to the externa, where many active bioprocesses were associated with development, bioprocesses enriched in the interna were mainly connected with metabolism. Comparative analysis revealed that 50 bioprocesses were active in at least two parts of the female body. Among them were “mitotic cell cycle process” (main trunk and thoracic part of interna), “cell division” (growing and thoracic part of interna), “homeostasis of number of cells” (same), “apoptotic signalling pathway” (same), “symbiotic process” (growing and main trunk of interna), “determination of adult lifespan” (same), “immune system process” (growing part, main trunk, and thoracic part of interna), and “autophagy” (same).

Figure 2. Gene set enrichment analysis (GSEA) results for the gene sets.

(a, b) Venn diagram for sets of parental bioprocesses enriched with either genes included in the molecular signature of the body part (a) or genes over-expressed in body part considered (b). (c) Clouds of enriched parental bioprocesses for the body parts. Most often, the bioprocess was found in the list of enriched bioprocesses, the larger the word size. Abbreviations: MS/Over – bioprocesses enriched with genes included either in molecular signature or in sets of over-expressed genes; PB – parental bioprocesses. (d– f) Histological sections of the main trunk of P. reticulata 1 - the wall of the main trunk; 2 - central lumen; 3 - groups of the floating cells; 4 - Nuage body. Scale bars: d - 200µm, e - 20µm, f - 50µm.

The “germ cell development” bioprocess was enriched in the main trunk and the growing part of the interna, whereas “female gamete generation” was only found among the lists of enriched bioprocesses in the thoracic interna part (Figure 2c). At the same time, according to our results, meiosis probably only occurred in the externa (Figure 2c). The results of histological studies also confirmed the presence of germ-like cells in the central lumen of the rootlets of the interna. Groups of floating small round cells with a high nuclear cytoplasmic ratio were found in the central lumen of the main trunk and peripheral rootlets (Figure 2d–f). A Nuage body (Figure 2d–f), which is a marker of germ cells, was present in each cell next to the nucleus.

No common bioprocess was found in the different parts of the female body, enriched with over-expressed genes. The genes with over-expression in the externa were involved in bioprocesses associated with cuticle transformation and development of the nervous system. In the growing interna part, over-expressed genes were involved in “gland development,” “response to nutrient levels”, “organic acid transport”, as well as lipid and fatty acid metabolic processes. In addition to various metabolic processes, “determination of adult lifespan” and “intrinsic apoptotic signaling pathway” were also found among a set of enriched bioprocesses for the main trunk. Bioprocesses associated with responses to various stimuli (bacterium/oxygen-containing compound/wounding), as well as “interspecies interaction between organisms”, “ion transmembrane transport”, “cellular homeostasis”, and “aging” were classified as “enriched” in the thoracic part of the interna.

All GSEA results are presented in Table S4 (Extended data⁷⁴).

Hundreds of genes encoding potential excretory/secretory proteins (ESP) were identified

The identification of potential ESP was performed in silico. A total of 852 "classical" and 282 "non-classical" ESP were found, which, respectively, had or did not have classical N-terminal signal peptides. Figure 3a, b shows Venn diagrams for the sets of genes encoding “classical” (Figure 3a) and “non-classical” (Figure 3b) ESP, respectively, which were included in the molecular signatures of body parts. Approximately 35% (297/852) of the genes encoding “classical” ESP had a noticeable expression level in all body parts considered. At the same time, more than half (143/282) of genes encoding “non-classical” ESP presented this expression pattern. In both cases, the externa had the largest number of “specific” ESP: 325 and 56 “classical” and “non-classical” ESP, respectively. Dozens of ESP (101 “classical” and 35 “non-classical”) were common for the three interna parts considered.

Figure 3. Gene set enrichment analysis (GSEA) results for the identified sets of potential excretory/secretory proteins (ESP).

(a, b) Venn diagram for sets of potential “classical” (a) and “non-classical” (b) ESP encoded by genes from molecular signatures of the body parts. (c) Clouds of parental bioprocesses enriched by genes encoding “classical” ESP. Most often the bioprocess was found in the list of enriched bioprocesses, the larger the word size. (d– e) Scanning electron microscope (SEM) photos of the interna of P. reticulata, 1 - rootlets; 2 – host tissues enlacing rootlets. (f) Confocal laser scanning microscopy (CLSM) photo of the interna of P. reticulata, scale bar 100µm, 3 - host tissues stained with antibodies against α-tubulin. Abbreviations: class/nonclassES – “classical” and “non-classical” ESP, respectively.

Both “classical” and “non-classical” ESP were divided into families based on their sequence similarity. For “classical” ESP, 35 families contained two to four proteins, while only two “non-classical” ESP were combined into one family. Most (579/852) of the “classical” ESP matched the MetazSecKB database, the hits being, e.g., mannose-binding proteins, serine proteinase, and cuticle proteins (Table S5, Extended data⁷⁵). Only 58 “non-classical” ESP matched MetazSecKB, of which 35 were “uncharacterized proteins” (Table S5, Extended data⁷⁵). All ESP were also compared to the NeuroPep database, with which only 14 “classical” ESP had hits (Table S5, Extended data⁷⁵). The latter included cerebellin-1, kininogen-1, insulin-like growth factors I and II, neuroparsin-A, nucleobindin-2, and five representatives of the serpin family.

Figure 3c shows bioprocesses enriched with genes encoding “classical” ESP. Among the body-part-specific bioprocesses were the “molting cycle” and “chitin-based cuticle development” in the externa, “positive regulation of cell communication” and “muscle organ development” in the growing interna part, and “immune response” and “regulation of apoptotic signaling pathway” in the main trunk and the thoracic part of interna, respectively. The majority (21/34) of enriched bioprocesses were common for two or more body parts considered. For example, “regulated exocytosis” (externa and growing part), “mesoderm development” (the growing part and the main trunk), “cell recognition” (same), “proteolysis” (all body parts considered), “motor neuron axon guidance” (same), “cell adhesion” (same), and “neuron recognition” (externa and thoracic part of interna). The activity of the bioprocesses associated with the involvement of the nervous system is consistent with the fact that trophic rootlets of P. reticulata were enlaced by a network of host's neurons marked by a presence of a-tubulin and serotonin (Figure 3d–f).

All ESP analysis results are presented in Table S5 (Extended data⁷⁵).

Significant differences between the TAI of body parts were revealed

Almost all (12,618 / 12,620) P. reticulata genes were distributed across 17 phylostrata, i.e., sets of genes that coalesce to founder genes having a common phylogenetic origin⁷⁶ (Table S6, Extended data⁷⁷) The three largest phylostrata were “Cellular organisms” (32.34%, 4,082 genes), “Eukaryota” (22.25%, 2,809 genes), and species-specific ones (13.7%, 1,729 genes) (Figure 4a). Less than 100 genes were assigned to “Ecdysozoa” (50 genes), “Panarthropoda” (44 genes), “Mandibulata” (63 genes), “Crustacea” (28 genes), and “Hexanauplia” (6 genes). The phylostratum “Cirripedia”, which included two barnacles, A. amphitrite and P. reticulata, consisted of 363 genes.

Figure 4. Phylostratigraphy and evolutionary transcriptomics results.

(a) Phylostratigraphic composition analysis results for P. reticulata reference gene set (“reference”), set of genes with noticeable (≥ 2 transcripts-per-million) expression in all female body parts considered (“common expr”), sets of genes over-expressed (externa/growing/main trunk/thoracic_overexpr) and sets of genes encoding potential “classical” (“classical_ESP”) and “non-classical” (“nonclassical_ESP”) excretory/secretory proteins (ESP). (b, c) Relative mean expression levels of phylostrata which occurred “before” (b) or “after” (c) division of Crustacea. (d) Transcriptome Age Indices (TAI) variation for female body parts considered. A lower TAI value describes an “older” transcriptome, whereas a higher TAI denotes a “younger” one. (e) The cumulative phylostrata contribution to the final (global) TAI profile.

The phylostratigraphic results were also used for composition analysis of various gene sets (Figure 4a). More than half of the genes with expression ≥ 2 TPM in all body parts considered belonged to “Cellular organisms” and “Eukaryota” phylostrata. The proportion of species-specific genes in this set was approximately 7% (490/6,829). Noticeable differences in the contributions of different phylostrata to the sets of over-expressed genes were found. For example, approximately 31% (88/283) of such genes in the thoracic interna part belonged to the species-specific phylostratum. In contrast, in other parts of the body, the proportion of species-specific genes from the total number of over-expressed genes did not exceed 16%. A complex phylostratigraphic composition was also revealed for genes encoding both "classical" and "non-classical" ESP. “Cellular organisms” phylostratum made the greatest contribution to the “classical” ESP (32.16%), while the species-specific phylostratum contributed most to the “non-classical” ESP (38.3%).

The phylostrata were divided into two groups: prior to the divergence of Crustacea (from “Cellular organisms” to “Crustacea”) and after this event (from “Multicrustacea” to “P. reticulata”). The relative expression patterns of the phylostrata are shown in Figure 4b, c. All phylostrata except the species-specific one had the highest relative expression in the externa, while the highest expression of the species-specific phylostratum was recorded in the thoracic part of the interna. At the same time, 14 out of 17 phylostrata had the least expression in the main trunk of the interna, the remaining three being “Bilateria”, “Pancrustacea”, and “Multicrustacea”.

One metric to quantify transcriptome conservation on a global scale is the TAI⁷⁸, which denotes the average transcriptome age throughout the biological process of interest⁶⁸. The TAI was measured for each part of the female P. reticulata body. The higher the value of the TAI, the greater the contribution of the “young” phylostrata. Significant differences between the TAI of different body parts were revealed. The TAI of the thoracic part of the interna was the highest (4.33), whereas the TAI of the growing part of the interna was the lowest (4.21) (Figure 4d). The partial contribution of the “P. reticulata” phylostratum to the TAI of the body part was about approximately three times greater than the contribution of any other phylostratum (Figure 4e).

The top-500 with the largest contribution to the TAI of body parts were extracted from the genes with the GO annotation and performed GSEA for these gene sets (Figure 5). Among the “specific” bioprocesses were “embryonic organ development” and “male sex differentiation” in the externa, “cell fate commitment involved in the formation of primary germ layer” and “positive regulation of chemotaxis” in the growing part of the interna, “formation of primary germ layer” and “gland morphogenesis” in the main trunk of the interna, and “response to external stimulus” and “cell population proliferation” in its thoracic part. Only four bioprocesses were common to all female body parts considered: “stem cell population maintenance”, “regulation of anatomical structure morphogenesis”, “chitin-based cuticle development”, and “animal organ morphogenesis”. Developmental processes were registered in each of the female body parts studied (Figure 5).

Figure 5. Gene set enrichment analysis (GSEA) results for top-500 annotated genes with the largest contribution to Transcriptome Age Index (TAI).

Clouds of parental bioprocesses enriched with top-500 genes with both Gene Ontology (GO) annotation and large contribution to TAI of the body part under consideration. The more often the bioprocess was found in the list of enriched bioprocesses, the larger the word size.

Phylostratigraphy and evolutionary transcriptomics results are presented in Tables S6 and S7 (Extended data^77,79), respectively.