Due to its status as an endangered species, sampling of P. nobilis was carried out with the permission of regional and national authorities for animal welfare (for Spain, Generalitat de Catalunya—Identificació de l'expedient: SF/0003/23; Bank of species of the Mar Menor INF/2020/0017 promoted by the General Directorate of the Mar Menor and the University of Murcia; for Italy, Prot. MATTM 0016478 05/03/2020). A total of 30 animals were included in the study, collected from July 2021 to May 2023. Sampling was performed on natural populations along the Mediterranean coasts of Italy (n = 4) and Spain (n = 18) and on animals maintained in captivity in two facilities in Spain (n = 8). The natural sites in Spain and Italy were selected according to the presence of remaining populations of P. nobilis, despite the occurrence of mortality outbreaks (Figure 1).

In Catalonia (Spain), samples were collected from two estuarine populations in the Ebro Delta, where residual populations still exist: Fangar Bay (northern hemidelta) and Alfacs Bay (southern hemidelta), which are distributed in four different monitored areas (28). Mortality was first observed in 2018, leaving only a few individuals by 2021 (37). The mortality rate progressively increased from 33.5% in September 2021 to 75% in June 2023 (Prado pers. comm). Some of the remaining animals from this area are maintained in the aquarium facilities of the Instituto de Investigación en Medio Ambiente y Ciencia Marina (IMEDMAR) of the Universidad Católica de Valencia (UCV), in Valencia, Spain, and were included in the study. Another remaining population in Spain is the one located at the Mar Menor lagoon (Murcia), in the southwestern part of the Mediterranean. In this area, P. nobilis mortality was first observed in the summer of 2016. Initially, the P. nobilis mortality was correlated with an environmental collapse that became critical in 2016, 2019, and 2021 and reduced the population by >99% (38); for the present work, two surviving localities were sampled in the Baron and La Manga areas (39); some animals (n = 2) from this area were transferred to the seawater tanks of the University of Murcia Aquarium. Animals (n = 6) had been kept in captivity at IMEDMAR and Murcia Aquarium facilities for 6 months to 1 year. Until recently, P. nobilis was widely distributed in the Venice Lagoon (Italy), with the largest and densest colonies over marine tidal deltas and the outer part of central basins (40). Mortality onset was observed at the end of 2020, when the population decreased between 60 and 80% (unpublished data). The animals (n = 4) were sampled near Ottagono Aberoni Island, approximately 1 km away from the Malamocco inlet (Sigovini, personal communication).

For each animal included in the study (from the field and in captivity), 2 ml of hemolymph was collected via a 23-gauge needle from the posterior adductor muscle as a non-destructive sampling. The presence of clinical signs of disease was considered during sampling. Previous reports in the literature indicate that outward signs of disease start with behavioral changes, including shell gaping, slow valve closure speed when touched, and mantle retraction into the valves from the edge of the shell (11). These signs were difficult to recognize in animals in the field and those maintained in captivity since they mostly appear at very late stages of disease (11). In Catalonia, on July 21 and July/August 2022, most of the animals in the study (19) displayed apparent slow valve closure. In Venice Lagoon, in May 2023, animals displayed quick valve closure and visible mantle at the border of the valves instead. Nevertheless, 2 months later, remarkable mortality was observed in the sampling area. Animals maintained in the Murcia Aquarium and IMEDMAR-UCV did not show apparent signs of disease, but in both cases, up to 40–60% of the animals suffered sudden mortality in the following 6 months, primarily attributed to manipulation from artificial spawning attempts.

After the procedure, the animals were immediately relocated to the seabed/tanks, and their health was monitored over the following days. A volume of ~2 ml of hemolymph was collected for each animal. A 1 ml hemolymph aliquot was placed in a 1 ml RNA Later (Sigma-Aldrich, Italy)-labeled tube placed on ice for molecular analyses, and the other half was preserved for ultrastructural analyses. A small hemolymph aliquot (20 μl) was placed on a hemocytometer (Burker chamber) to determine the total cell number ( × 10⁵/ml) and define the total hemocyte count (THC) (Table 1).

For all the animals included in the study, part of the hemolymph was also fixed in 2.5% glutaraldehyde in PBS for 3 h. Cells were postfixed in 1% osmium tetroxide for 60 min and 0.25% uranyl acetate overnight. The cells were then dehydrated through a graded ethanol series and embedded in EPON resin overnight at 37°C, 1 day at 45°C, and 1 day at 60°C. Ultrathin sections (80 nm) were cut parallel to the substrate and placed onto a 200-mesh copper grid. Bright-field TEM images were obtained on the dried sample by using an FEI TECNAI G2 200 kV s-twin microscope operating at 120 kV (Thermo Fisher Scientific, Waltham, USA). Digital images were acquired with an Olympus VELETA camera.

Total RNA was extracted from the hemocytes of six samples of P. nobilis collected from three sites of the wild population in Alfacs Bay, Catalonia (Table 2) using the PureLink RNA Mini Kit (Thermo Fisher Scientific, Waltham, Massachusetts, USA). Following DNase treatment, RNA concentration and quality were measured using a Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific). Strand-specific cDNA libraries were prepared, and Illumina sequencing (NovaSeq 6000, 150 bp paired-end) was carried out by Eurofins Genomics (Ebersberg, Germany). Raw reads underwent trimming and adapter clipping using Trimmomatic 0.4 (41).

To identify RNA viruses eventually infecting the P. nobilis hemocytes, the trimmed, good-quality reads were aligned to the P. nobilis genome (ASM1616189v1) using Bowtie2 (42), and the unmapped reads were extracted using SAMtools (43).

A total of 33,097,808 paired unmapped reads were assembled using Trinity 2.12.0 (44), and the abundance estimation of contigs was performed with RSEM (45).

The Trinity-assembled contigs (140,922) were used as queries in a BLASTX search against the viral protein database (https://ftp.ncbi.nlm.nih.gov/refseq/release/viral/, downloaded on 6 June 2023). To focus the annotation on RNA viruses, the assembled contigs were analyzed using Virsorter2 (46) by selecting only the classification of RNA viruses. Additionally, VirBot, an RNA viral contig detector for metagenomic data (47), was utilized with the sensitive option. The contigs identified as positive by both VirSorter2 and VirBot were subjected to analysis using ORFfinder (https://www.ncbi.nlm.nih.gov/orffinder/). The resulting ORFs were then subjected to a BLASTP search against the UniProtKB/Swissprot, RefSeq protein, and non-redundant protein sequence databases on 4 July 2023.

The predicted sequences of the RNA-dependent RNA polymerase (RdRp) protein encoded by the identified RNA viruses were used for BLASTP searches against the viral protein database, and the highest scoring homologous sequences were downloaded. Multiple amino acid alignments of the selected RdRp sequences were performed using the Constraint-based Multiple Alignment Tool (COBALT, https://www.ncbi.nlm.nih.gov/tools/cobalt/re_cobalt.cgi). The best amino acid substitution model was searched, and the maximum likelihood tree was generated using the LG + G + I model and 500 bootstrap replicates using MEGAX (48), resulting in a final dataset of 317 characters.

Total RNA was extracted from the hemocytes of other 5 animals P. nobilis samples collected from different areas in Alfacs Bay, Catalonia (Table 3), as described before. RNA (500 ng) was reverse-transcribed using Maxima First Strand cDNA Synthesis (Thermo Scientific). Specific primer pairs were designed to amplify the two most abundant RNA viruses identified by the in silico analysis (Table 4), and quantitative real-time PCR experiments were conducted. Amplification reactions were performed in a total volume of 10 μl using 5 μl of PowerUp SYBR Green Master Mix (Applied Biosystems), 0.2 μM of each specific primer, and adjusted to 10 μl with distilled water. The reactions were conducted in technical duplicates using P. nobilis elongation factor-1 (EF-1) as a reporter gene (49, 50) (Table 4). The amplification conditions were as follows: 1 cycle for 2 min at 50°C; 1 cycle for 2 min at 95°C; 40 cycles of amplification at 95°C for 15 s; 60°C for 18 s; and 72°C for 1 min. The reactions were conducted in the QuantStudio 3 Real-Time PCR System (Thermo Scientific). The cDNA of the samples identified as Trab8a and XIR2a (the same used in the RNA-seq experiment) was subjected to real-time PCR amplification to confirm the results of the in silico analysis, and two negative controls were run without cDNA.

The THC values varied among individuals and areas over the years. The mean hemocyte count showed the lowest values from animals maintained in captivity in the IMEDMAR-UCV and Murcia Aquarium (between 3.5 × 10⁴ and 1.60 × 10⁵ ml⁻¹ cells) compared with those from the natural population in Catalonia (1.90–2.42 × 10⁵ ml⁻¹ cells) (Table 1). The most represented cells were granulocytes (~10 μm), displaying a small peripheral nucleus, and a small quantity of hyalinocytes (~7 μm), presenting a large central nucleus and small cytoplasm. All the hemocytes collected from animals in all areas of Italy and Spain from 2021 to 2023 (100%) displayed ultrastructural features of viral infection with associated cell death. These hemocytes all exhibited complex and unique membrane relocations, displaying numerous cytoplasmic vesicles associated with damaged mitochondria, a total lack of granules, the presence of protein aggregates, and, in some cases, cellular buddings. The formation of invaginations occurred at the membrane of various organelles, including the ER, endolysosomes, and mitochondria (Figures 2, 3). Features of viral infection at the hemocyte level are represented in Figure 2.

In all immune cells, aggregates of modified membranes occupied large areas of the perinuclear cytoplasm. Numerous membranous vesicles forming virions at various stages of self-assembly were visible. Vesicles in the cytoplasm were clearly associated with a dilated ER and nuclear membrane (Figures 3C, D). These vesicles had different dimensions and were constituted by double membranes, reported in the literature as DMVs or as autophagosome-like vesicles. The small vesicles measured 80–100 nm, whereas larger DMVs measured ~800–1,000 nm (Figures 3C–E). The presence of viral particles related to cytoplasmic vesicles, or mitochondrial membranes, was visible in all immune cells of P. nobilis (Figure 3F). Smaller vesicles were mostly used for virus genome replication at the cytoplasmic level, whereas larger DMVs supported a later step of virus production, specifically provirion maturation and the formation of complex factories. The DMVs displayed evident double or more complex membranes (Figure 4A). In earlier phases, the DMVs were filled by an unassembled viral amorphous and granular material; vesicles then showed one or two ring-like structures captured in the lumen vesicle, as also reported in other PVs (34) (Figure 4B). Viral replication factories are typically constituted by complex membrane remodeling emerging from DMVs formed by viral particles and proteins assembling mature virions, in some cases packed together. The virion is an icosahedral, non-enveloped, small (20 nm) particle with no discernible projections (Figures 4C, D).

Once the virus is assembled, viral release can occur through different mechanisms. Viruses can be released without lysis through so-called secretory autophagy, which occurs after the fusion of DMVs with the plasma membrane and are released as exosomes. Exosome typical dimensions ranged between 800 and 1,000 nm. Exosomes were secreted after their fusion with cell membrane vesicles and released as single-membraned virus-filled vesicles (Figures 4E, F). The four samples from Venice Lagoon, defined in the field as apparently healthy, displayed mature viral particles spreading into the cytoplasm (Figure 4D), small DMVs, and active factories.

The predominant infective strategy observed in the PV of P. nobilis, typical for non-enveloped viruses, was a non-lytic “unconventional secretion” using extracellular vesicles (EVs) (Figure 5). PV-infected cells secreted EVs carrying viral RNA particles. Vesicles ranged in diameter from 100 to 700 nm. Virions may be released from cells by a budding process from both granulocytes and hyalinocytes (Figures 5A–C). Infective vesicles were observed carrying either mature viral particles (Figure 5D) or vesicle factories containing still immature virions and DMVs (Figures 5E, F). EVs may be taken up by recipient cells by endocytosis or fusion with the plasma membrane. The vesicle is internalized, via endocytosis, in an endocytic vacuole surrounded by DMVs to form mature viral particles (Figures 6A, B). Infective vesicles transporting mature virions and injecting their contents were also visible (Figures 6C, D).

Dead hemocytes were observed during infections. The highest number of damaged cells was observed in animals maintained in captivity at IMEDMAR-UCV and Murcia Aquarium, surrounded by an elevated number of EVs or containing endocytic vesicles. Damaged cells appeared empty, shrunk, contained few vesicles, and presented membrane rupture and apoptotic bodies (Figure 7).

Sequencing of six stranded cDNA libraries retrieved from hemocyte samples of six P. nobilis using Illumina technology generated a total of 162,139,300 paired-end raw reads (Table 2), subsequently deposited in the NCBI Sequence Read Archive (SRA) under the BioProject accession PRJNA999583. After trimming and mapping against the P. nobilis genome, the number of unmapped paired reads was 33,097,808 (20.4%). These unmapped reads were used to assemble the viral RNA genomes present in the P. nobilis hemocytes.

BLASTX analysis of the 140,922 contigs assembled using Trinity resulted in 61,634 positive matches (43.7%, e-value < 10–20) against the viral protein database, encompassing DNA and RNA viruses as well as bacteriophages. Since this study was performed to identify RNA viruses in P. nobilis hemocytes, we focused specifically on this viral group using VirSorter2 and VirBot software tools. Both analyses identified the same set of five contigs, ranging from 5,210 to 8,876 nucleotides (nt) in length, all classified, with a score of 1 (the highest possible score with Virsorter2), into the realm Riboviria.

These five sequences exhibited the highest BLASTX matches with different Riboviria polyproteins 1 and 2 (Table 5), mainly of Picornavirales. Notably, two of these contigs (TRINITY_DN35348_c1_g1 and TRINITY_DN9545_c0_g1) displayed the highest number of mapped reads, particularly in the Trab samples (Site 1, Trabucador), indicating their abundance among the five identified viruses. Additionally, during the global BLASTX analysis, we discovered the contig TRINITY_DN58105_c0_g1_i1 (841 nt), which matched the predicted replication-associated protein of Chaetoceros tenuissimus RNA virus type II (BAP99820.1). Manual overlap with TRINITY_DN35348_c1_g1 (5,210 nt, matched with the predicted structural protein of Chaetoceros tenuissimus RNA virus type II, BAP99821.1) resulted in a 6,023 nt consensus sequence.

Figure 8 presents the genomic organization of the two most prevalent viruses identified herein. The sequence consensus_DN35348_c1_g1_DN58105_c0_g1 (A) is partial, lacking the complete 5′-end, while the DN9545_c0_g1 sequence (B) is complete. These two sequences were deposited in GenBank under accessions OR448788 and OR448789. The ORFfinder analysis identified two large ORFs encoding two viral polyproteins in both sequences. Within polyprotein 1, there were conserved domains of RdRp and RNA helicase (the latter was detected only in sequence DN9545). In contrast, structural polyproteins 2 have two rhv-like domains (picornavirus capsid protein domain-like), a CRPV capsid domain (a family of capsid proteins found in positive stranded ssRNA viruses), and a short VP4 domain (a family of minor capsid proteins located within the viral capsid at the interface between the external protein shell and packaged RNA genome). The presence of two large ORFs encoding these specific domains suggests that these viral sequences belong to the Picornavirales order.

The genomic organization of the remaining three viruses identified in P. nobilis hemocytes found herein (DN33054, DN37302, and DN33329) was deposited in GenBank under accession numbers OR448790, OR448791, and OR448792, respectively (Supplementary Figure 1). Similar to the two most abundant viruses, the genomic organization of DN33054 and DN37302 presented two large ORFs, with ORF1 encoding RdRp and RNA helicase and ORF2 encoding two rhv-like domains, a CRPV capsid domain and a short VP4 domain. In contrast, the genomic organization of DN33329 was quite different, with three partially overlapping main ORFs encoding a peptidase and RdRp, similar to the structure of the Sobelivirales order.

The maximum likelihood tree was generated based on the alignment of the conserved region of the RdRp proteins of the viruses identified in the P. nobilis hemocytes with homologous sequences downloaded from GenBank (Figure 9). Almost all of the RdRp sequences downloaded from GenBank belonged to Picornavirales from marine or freshwater environmental samples. A restricted number of sequences belonged to Sobelivirales, and two of them (Penaeus vannamei picornavirus and Wenzhou shrimp virus 8) to Dicistoviridae, which are associated with mollusk and crustacean diseases (27, 51). Among the viral sequences identified in P. nobilis hemocytes, DN33329 clustered with those in the order Sobelivirales, along with viruses from rivers, soil, and plants, agreeing with the genomic organization previously described (Supplementary Figure 1). Additionally, the remaining four sequences from P. nobilis belonged to the order Picornavirales, as expected from their genomic organization; among them, the two most abundant were clustered within the family Marnaviridae.

In terms of ultrastructure, both samples from the Trabucador area (Trab8a and Trab8b) that had the highest number of mapped viral reads displayed numerous cytoplasmic viral factories and extracellular infective vesicles spreading into the hemolymph, suggesting an active infective phase.

The quantitative real-time PCR experiments conducted on the cDNA of the P. nobilis hemocytes collected at different sampling points in Catalonia were positive for the primer pair PNChetoF/R designed on the reconstructed DN35348_DN58105 sequence, similar to a Chaetoceros tenuissimus RNA virus type II. These results confirm the in silico analysis of viral read abundance (Table 5), with Trab8a and XIR2a showing the highest infection level and XIR5 showing a low infection level (Figure 10). With the same primer pair, viral RNA was detected in all other examined samples, presenting variable infection levels. The primer pair PNPicorF/R designed for the sequence DN9545 did not provide a detectable amplification signal.