28 ISJ 19: 28-36, 2022 ISSN 1824-307X RESEARCH REPORT Identification of five picorna-like viruses associated with the endangered cave- dwelling bivalve Congeria kusceri (Bole, 1962) A Scapolatiello1, U Rosani2, C Manfrin1, S Puljas3, A Pallavicini1, M Gerdol1* 1University of Trieste, Department of Life Sciences, Trieste, Italy 2University of Padova, Department of Biology, Padova, Italy 3University of Split, Faculty of Science, Split, Croatia This is an open access article published under the CC BY license Accepted November 16, 2021 Abstract Congeria kusceri is a bivalve mollusk species endemic to the Dinaric Karst, which displays unique adaptations that have allowed its survival in the subterranean environment with small morphological changes compared with its fossil relatives. Anthropic activities have recently impacted the surface flow of the Neretva river, impairing the seasonal water cycle that has characterized the habitat of this species for hundreds of thousands of years. The lack of an adequate water supply, together with pollution from agricultural and farm water runoff, are posing a serious threat to C. kusceri, as evidenced by the sharp population decline observed in several locations during the past few decades. Due to the limited knowledge available about the basic biology of this filter-feeding species, the precise factors that may affect its health status and reproduction and therefore represent a hazard for its conservation are still unclear. Here, through a transcriptomic approach, we describe the nearly- complete genomes of five C. kusceri-associated RNA viruses belonging to the Picornaviridae family and phylogenetically related with picorna-like viruses previously described in other Mollusca. Although it is presently unknown whether these viruses may have a detrimental effect on bivalve health, we observed a significant increase of viral load during the summer season. Key Words: virome; Karst; subterranean; Picornaviridae; transcriptome; RNA-sequencing Introduction Congeria kusceri (Bole, 1962) is a freshwater bivalve endemic to the Dinaric Karst, a vast region geologically characterized by the presence of carbonate rocks and located in the Balkan Peninsula, between the Adriatic Sea and the Pannonian Basin (Zupan Hajna, 2019). C. kusceri is currently listed as a vulnerable (VU) species in the European Red List of non-marine mollusks (Cuttelod et al., 2011) and as a critically endangered (CR) species in the Croatian Red List of subterranean fauna (Bilandžija and Jalžić, 2009). A member of a once widespread genus of freshwater bivalves that inhabited the Paratethys Sea, C. kusceri is a living fossil which maintained an apparent morphological stasis since the end of Miocene (approximately 5.3 Mya). The development of several unique adaptations, such as a K-selected reproductive strategy, larval brooding, and a long _________________________________________ Corresponding author: Marco Gerdol University of Trieste Department of Life Sciences Via Giorgieri 5, 34127 Trieste, Italy E-mail: mgerdol@units.it lifespan (Morton and Puljas, 2013), allowed the survival of this species in the subterranean environment to the present day (Stepien et al., 2001; Bilandžija et al., 2013). C. kusceri, whose current range of distribution includes just eight known sites in caves in to the Neretva river basin in Croatia and Bosnia and Herzegovina, was long thought to be the only member of a single pan- Dinaric genus. However, recent molecular studies have established that two other reproductively isolated congeneric species survive in the Dinaric karst, i.e. Congeria mulaomerovici and Congeria jalzici (Bilandžija et al., 2013). The adaptation of C. kusceri to a highly stable cave environment poses a serious threat for its survival, since this species is not expected to be able to withstand significant long-term environmental alterations. In particular, the annual cycles which typically characterize the waters influx in Karst caves play an important role in regulating shell growth, reproduction and other important aspects of Congeria’s life cycle (Morton and Puljas, 2013; Puljas et al., 2014; Jovanović Glavaš et al., 2016). The building of dams, hydropower plants and channels for agricultural land use over the past few 29 decades have undoubtedly modified the water regime of the Neretva river basin, with a severe impact on the amount of water that seasonally floods the cave systems where C. kusceri lives. These alterations are now threatening the survival of this species, which has already disappeared from some of the caves where large populations had been historically documented and is quickly declining in others (Bilandžija et al., 2021). Moreover, agriculture and farm runoff can lead to organic and inorganic pollution in the subterranean environment (Sasakova et al., 2018). This, together with the progressive salinization of the lower parts of the Neretva basin due to seawater infiltration, poses a further threat for the unique habitat of Congeria, which is thought to have remained unchanged for hundreds of thousands of years (Vranješ et al., 2013). Such environmental changes could have a significant impact on the cave microbial communities, as well on the Congeria- associated microbiota, primarily consisting of chemolithoautotrophic bacteria (Rigonni et al., 2021). It is presently unknown whether alterations in the microbiome and in the virome found in Karst subterranean waters can pose a threat for C. kusceri survival and no scientific literature is available on this subject. Nevertheless, alterations of the bacterial and viral communities associated with bivalve tissues have been clearly linked with mass mortality events and with a general decline of the health of marine bivalve populations (Barbosa Solomieu et al., 2015; Lasa et al., 2019). Picornaviridae are among the most frequently described viruses associated with bivalve mollusks. Even though their detection is often the result of bioaccumulation from ingested food particles or pathogens of bivalve associated parasites (Crespo- González et al., 2008), some of them can infect mollusks and replicate in their tissues (Renault, 2016). In particular, ultrastructural studies have identified on several occasions Picorna-like viral-like particles (VLPs) in bivalves affected with granulocytomas ( Rasmussen, 1986; Jones et al., 1996) and numerous reports have suggested that Picorna-like viruses may be the main responsible or contributors to mass mortality events observed in different species of marine bivalves (Carballal et al., 2003; Renault and Novoa, 2004). It is however still unclear whether these RNA viruses should be just considered as bivalve pathogens or rather as a component of the virome associated with these filter-feeding organisms under normal physiological conditions, considering their frequent detection though metatranscriptomic approaches in absence of detectable disease (Moreira et al., 2012; Rosani and Gerdol, 2017; Rosani et al., 2019). The information available about Picornaviridae infections in freshwater bivalves is extremely limited. High viral loads of different Picorna-like viruses have been previously observed in association with mass mortality events of Actinonaias pectorosa in the Clinch river (USA), even though they were unlikely to be the cause of the pathology (Richard et al., 2020). Although other picornaviruses have been also tentatively indicated as potential pathogens in North American freshwater bivalves (Goldberg et al., 2019), on some occasions they have been linked with concurrent infections by parasitic trematodes (Ip and Desser, 1984). Considering the role potentially played by Picornaviridae as bivalve pathogens and the high susceptibility of cave-dwelling animals to environmental alterations, the detection of these viruses in C. kusceri may represent a potential threat for the conservation of this highly endangered species. Materials and methods Viral genome identification and characterization The transcriptome of Congeria kusceri was generated through the de novo assembly of RNA- sequencing data obtained from multiple tissues (gills, mantle, gonads, digestive gland and adductor muscle). All tissues were carefully dissected with the aid of a scalpel, following the severing of the adductor muscle which allowed valve opening, from a total of ten individuals sampled between July and September 2018. The full methodology is described in detail elsewhere (Scapolatiello et al., 2021). This sequence dataset was screened for the presence of contigs encoding proteins containing a RNA- dependent RNA polymerase (RdRp) domain. In brief, the viral RdRp HMM profiles belonging to the clan RdRP (CL0027) were obtained from Pfam (Finn et al., 2009) and used as queries for HMMER v.3.3.2 (Finn et al., 2011) against the virtually- translated proteome of the species (obtained with TransDecoder, https://github.com/TransDecoder/). Positive matches were identified with an e-value threshold of 1x10-5. The reliability of the selected contigs was evaluated by the visual inspection of the uniformity of mapped Illumina reads along the full length of the contigs and double checked through the comparison with an alternative de novo assembly obtained with metaSPAdes v.3.10 (Nurk et al., 2017), which fully confirmed the sequences obtained. The presence of Open Reading Frames (ORFs) was evaluated with the CLC Genomics Workbench 20 (Qiagen, Hilden, Germany) and the encoded proteins were: (i) used as queries for a BLASTp search (Altschul et al., 1990) against the NCBI nr protein database, looking for significant homologies with sequences from previously described viruses; (ii) screened with InterProScan v.5 to detect conserved protein domains (Jones et al., 2014); (iii) analyzed with HHPRED, looking for structural similarities not detectable using the two methods described above (Söding et al., 2005). This approach allowed the identification of five viral genomes, whose completeness was assessed based on the pairwise comparison with the sequences of the closest available relatives detected in public databases. Phylogenetic analyses Since the five picorna-like viruses identified in C. kusceri displayed a different genome architecture and were all incomplete, the most conserved protein domain in Picornaviridae, i.e. the RdRp, was extracted to create a multiple sequence alignment (MSA) with MUSCLE v.3.8.31 (Edgar, 2004). The regions corresponding to the RdRp domain were similarly extracted from known viruses displaying a 30 significant homology with the five genomes under scrutiny, based on the detection of a significant BLASTp and pairwise identity at the amino acid level, higher than 40%. The previously published transcriptome of the closely related species Dreissena rostriformis bugensis (Péden et al., 2019), deposited in the NCBI Transcript Shotgun Assembly database under the accession ID GHIX00000000.1, was similarly screened to detect other uncharacterized viruses. Redundant sequences (i.e. those displaying a pairwise similarity higher than 90%) were removed with CD-HIT (Li and Godzik, 2006). The MSA was trimmed with Gblocks v.0.91b in order to remove highly divergent and phylogenetically uninformative regions (Talavera and Castresana, 2007). Then, this sequence set was analyzed with modeltest-ng v.0.1.3 (Darriba et al., 2019) to identify the most appropriate model of molecular evolution according with the corrected Akaike Information Criterion (Cavanaugh, 1997), which was found to be a LG model (Le and Gascuel, 2008), with a proportion of invariable sites and a gamma-distributed rate of variation among sites. Phylogenetic inference was carried out with MrBayes v.3.2.7a (Huelsenbeck and Ronquist, 2001), using two independent Monte Carlo Markov Chain (MCMC) analyses in parallel running for 250,000 generations each, sampling trees each 1,000 generations. The first 25% the obtained trees were discarded during the burn-in process and the convergence of all the estimated parameters of the model was evaluated with Tracer v.1.7.1 (Rambaut et al., 2018). Quantification of viral abundance The trimmed Illumina reads obtained from the 10 samples available were mapped against the five reference viral genomes with the CLC Genomics Workbench v.20, using stringent parameters (length fraction = 0.75, similarity fraction = 0.98) to avoid non-specific matches. Based on the number of reads mapped on each viral genome, the average read coverage per base was calculated for each virus and each sample, obtaining a normalized measure of abundance that could allow a direct comparison among samples. In detail, the five sampled tissues were gills, mantle, gonad, digestive gland and adductor muscle, which were taken from individuals collected in the Pit in Predolac in the Neretva river basin (Croatia) in July and September 2018, i.e. at the beginning and the end of summer, two critical time points for the life cycle of C. kusceri, as well as for the annual water cycle of the Karst subterranean environment. Each sample derived from a pool of five different individuals. The sequencing data analyzed in this study have been deposited in the Sequenced Read Archive (SRA) under the BioProject PRJNA704250 (SRR13767952-SRR13767961). Results and discussion Viral genome organization We identified five RNA viruses characterized by the presence of a complete ORF encoding a protein with a RNA-directed RNA polymerase (RdRp) domain. Upon further scrutiny, and based on the homology detected with other known Picornaviridae, all five genomes were found to be partial, as they were missing a region of variable length at the 5’- end. Since the library preparation protocol used a poly-A selection step, this result was expected, due to the enrichment of Illumina reads towards the 3’- end of poly-adenylated mRNAs (Love et al., 2016), including those of viral origin (Rosani et al., 2019). The five viruses were named “Congeria picorna-like- virus”, followed by a progressive number (from 1 to 5), thereafter abbreviated as Cplv1, Cplv2, Cplv3, Cplv4 and Cplv5. Overall, while the genomes of all five viruses encoded, as expected, the presence of a RNA- directed RNA polymerase (RdRp), a picornain 3C peptidase and a number of structural proteins, their organization largely varied. In line with previous reports about other Picornaviridae isolated from invertebrates (Shi et al., 2016), structural and non- structural proteins were organized in polyproteins. These could be either encoded by a single or by two distinct ORFs and they could be interchangeably encoded either at the 5’ or the 3’ ends of the viral genome. However, the regions encoding the picornain and RdRp domains were always detected in pair (with the peptidase located in a 3’ position compared with RdRp). In detail, the sequence of Congeria picorna-like virus 1 (Cplv1) was 6,133 bp long and harbored two ORFs. Based on the comparison with its known close relatives, the assembled viral genome was expected to lack approximately 3 Kb at its 5’ end, corresponding to the N-terminal region of the replicative polyprotein, which contained the picornain and RdRp domains. The second complete ORF encoded the structural polyprotein (Fig. 1). Unlike Cplv1, the 5,614 bp genome of Congeria picorna-like virus 2 (Cplv2) included a single partial ORF, which was predicted to lack about 3 Kb of coding sequence. The architecture of the encoded polyprotein included a picornain/RdRp domain pair in the N-terminal part and the structural proteins of the capsid at the C-terminal end (Fig. 1). The genome of Congeria picorna-like virus 3 (Cplv3) was 7,585 bp long and was predicted to lack approximately 1 Kb of sequence at its 5’ end, likely making it the most complete out of the five viral genomes identified in the present work. The first incomplete ORF encoded the structural polyprotein, whereas the second complete ORF encoded a polyprotein which included the picornain/RdRp domain pair at the C-terminal end (Fig. 1). In addition, a portion of this polyprotein showed structural homology with the 2C non-structural protein of other Picornaviridae, which bears ATPase and RNA helicase activities (Cheng et al., 2013) and is thought to be involved in the inhibition of the TNF-α-mediated activation of NF-κB in the host (Li et al., 2016). The assembled genome of Congeria picorna- like virus 4 (Cplv4) was the shortest out of the five genomes identified in this study, with a length of 3,817 bp. The only detectable ORF encoded the picornain/RdRp domain pair in its C-terminal part. Similar to Cplv3, Cplv4 had a region with structural homology with the 2C non-structural protein (Fig. 1). It is presently unknown whether the capsid proteins 31 Fig. 1 Schematic representation of the five partial viral genomes identified in the present study. The serrated margins at the 5’ end indicate missing sequence. The main detectable viral domains (i.e. the capsid proteins, the RNA-directed RNA polymerase, the 3C peptidase and the 2C ATPase/helicase) are indicated. White regions indicate the spacer regions located between ORF1 and ORF2 and the 3’UTR in this virus are encoded by the N-terminal part of the same polyprotein or by a second ORF located in the 5’ missing end of the viral genome. The genome of Congeria picorna-like virus 5 (Cplv5) was 5,230 bp long and it was predicted to lack ~3 Kb of sequence at its 5’ end. The first, incomplete ORF displayed the usual picornain/RdRp domain pair at the C-terminal end, whereas the second complete ORF encoded the structural polyprotein (Fig. 1). Phylogenetic placement of the Congeria picorna-like viruses Although all the five Congeria picorna-like viruses were novel, Cplv3 was placed in a long branch of the phylogenetic tree, distantly related with most known Picornaviridae, but displayed a very high sequence identity (75-90% at the amino acid level) with a small group of viruses previously described in freshwater environments and occasionally found associated to Mollusca (Fig. 2). In detail, Cplv3 was phylogenetically related with a few nearly-identical viruses, placed in a single cluster by CD-HIT and indicated by “freshwater picorna-like virus cluster 1” in Figure 2. This included the Trichosanthes kirilowii picorna-like virus, the Forsythia suspensa picorna-like virus (Yang et al., 2021), the snail Beihai picorna-like virus 105/niflavirus (Ng et al., 2012; Shi et al., 2016) and with the freshwater mussel Clinch dicistro-like virus 1 (Richard et al., 2020). On the other hand, Cplv1 belonged to a large monophyletic and heterogeneous group of Picornaviridae (Fig. 2). The closest know relative to Cplv1 included the Beihai picorna-like virus 70, 71 and 72, as well as the Wenzhou picorna-like virus 26 and Atrpec virus 1, which displayed sequence homology at the amino acid in the range of ~40% for the structural polyprotein and ~50% for the replicative polyprotein. These viruses have been reported in different invertebrates, including Cnidaria and Crustacea, but also in Mollusca (i.e. the Beihai picorna-like virus 70 and Wenzhou picorna-like virus 26) (Shi et al., 2016). Intriguingly, albeit quite distantly related with Cplv1, we found that another virus detected in D. rostriformis bugensis (GenBank accession ID: GHIX01007600.1), one of the closest extant relatives of Congeria spp. (Stepien et al., 2001; Bilandžija et al., 2013;), was placed in the same large clade of Picornaviridae (Fig. 2), which may therefore show some degree of host specificity to Mollusca. Cplv4 was placed in a different, well distinct group of Picornaviridae (Fig. 2) and only showed weak sequence homology with known picorna-like viruses isolated from several different invertebrates, which most notably included Mollusca, such as in the case of the Beihai mollusks virus 1 and 2 (Shi et al., 2016). Interestingly, albeit distantly related, two viruses detected in D. rostriformis bugensis (GenBank accession IDs: GHIX01044252.1 and GHIX01014380.1) also belonged to this clade, which 32 Fig. 2 Unrooted phylogenetic tree displaying the evolutionary relationship between the five Congeria picorna-like viruses (indicated by red dots), other picorna-like viruses identified in the D. rostriformis bugensis transcriptome (indicated by blue dots) and other known Picornaviridae. The tree is represented as a 50% majority consensus tree and posterior probability support values are only shown for major nodes. The scale bar indicates the number of substitutions per site suggests that these viruses may be part of the regular virome of Dreissenidae and possibly of other bivalves (Fig. 2). Cplv2 and Cplv5 did not show any significant match with known Picornaviridae, only showing sequence homologies with less than 30% identity at the amino acid level. They were therefore both placed as orphan taxa in long branches of the phylogenetic tree, which nevertheless created a highly supported monophyletic clade with the Cplv3 and the freshwater picorna-like virus cluster 1 (Fig. 2). Tissue and temporal distribution of the five C. kusceri picorna-like viruses Overall, the five viruses displayed a very similar pattern of relative abundance, both across tissues and between the two sampling time points, even though it was not possible, due to the pooling strategy used prior to RNA-sequencing, to evaluate whether such abundances were influenced by significant inter-individual differences. All viruses could be detected in C. kusceri in July 2018, albeit at relatively low levels. Higher viral loads were identified in the mantle tissue, followed by gonads and digestive gland, with little or no detectable signal in the gills and adductor muscle (Fig. 3), suggesting that the viruses were poorly replicating in bivalve tissues. In September 2018, higher viral loads were recorded in all tissues, with an increase by several orders of magnitude in the digestive gland (by ~50X for Cplv5 and by > 100X for the other 4 viruses). Cplv3 was the most abundant in this tissue, reaching an average read coverage close to 800X. The relative abundance of Cplv3 was ~3.7 higher than Cplv1, ~5.5X higher than Cplv2, ~13.1X higher than Cplv4 and ~17.4X higher than Cplv5. Four out of the five viral sequences were identified as differentially expressed transcripts in the statistical analysis for differential gene expression (a Kal’s Z-test on proportions (Kal et al., 1999) carried out between the September 2018 and July 2018 digestive gland samples (Scapolatiello et al., 2021). In detail, differential expression was 33 Fig. 3 Relative abundance of the 5 Congeria kusceri picorna-like viruses in five different tissues in June and September 2018, calculated based on the number of reads mapped on the reference genomes supported by false Discovery Rate-corrected p- values equal to 0 (Cplv3), 1.82 x 10-12 (Cplv1), 4.26 x 10-9 (Cplv2) and 8.31 x 10-3 (Cplv5). The significant viral loads recorded in the digestive gland may have different interpretations. In absence of ultrastructural observations, the precise location of the viral particles cannot be defined, leaving three alternative explanations open for the detection of viral RNA: (i) bioaccumulation of the virus from exogenous sources by filter-feeding; (ii) association of the virus with parasites or symbionts of C. kusceri; (iii) active viral replication in the digestive gland of C. kusceri. The bioaccumulation of Picornaviridae by filter- feeding has been documented in bivalves on several occasions, in particular for what concerns wastewater-borne human enteric viruses, even though such viruses are likely unable to infect bivalves and replicate in their tissues (Hansman et al., 2008). It is therefore possible that the five picorna-like viruses described in this work are the result of influx of surface waters containing virus- associated food particles in the caves, with no direct impact on the health of C. kusceri. Nevertheless, this hypothesis is challenged by relatively poor influx of water from the external environment which occurs during the summer season at the Jama u Predolcu (Croatia). The possibility that the detection of these five viruses was linked with the presence of other parasitic infections, as previously reported in some marine species (Crespo-González et al., 2008), would certainly deserve further study. Although several freshwater bivalve species have been previously identified as viable trematode hosts (Brian and Aldridge, 2019; Taskinen, 1998; Marszewska and Cichy, 2015; Müller et al., 2015), no scientific literature is available on this subject for C. kusceri. As of note, no 18S/28S rRNA or COI markers suggestive of the presence of eukaryotic parasites was detected in the assembled transcriptome. Future studies should also aim at investigating possible host-virus links, such as the presence of Endogenous Viral Elements (EVEs) (Blair et al., 2020), or at mapping CRISPR spaceromes (Shmakov et al., 2020) in the genome of C. kusceri, once this resource will become available. Concerning the possibility of active viral replication in the digestive gland, it needs to be remarked that histological lesions linked with picorna-like VLPs have been previously observed in this tissue in other bivalve species, supporting the idea that the digestive gland may serve as a preferential viral replication site also in C. kusceri (Hine and Wesney, 1997). Conclusions The five different picorna-like viruses identified in C. kusceri represent novel additions to the poorly known group of bivalve-associated RNA viruses and provides some insights in the virtually unexplored viral landscape associated with subterranean waters. Moreover, our bioinformatics-based detection approach confirms the usefulness of analyzing RNA-sequencing datasets for the identification of uncharacterized RNA viruses in bivalve hosts (Rosani and Gerdol, 2017; Rosani et al., 2019). 34 It is worth noting that Cplv2, Cplv3 and Cplv5 belonged to a well-supported monophyletic clade that only comprised a few known viruses previously isolated from freshwater sources, which may therefore represent an understudied group of picorna-like viruses associated with freshwater bivalves, and possibly with other freshwater invertebrates. All viruses displayed a marked increase in abundance throughout the summer season, which clearly indicates that this threatened bivalve species can be exposed to high viral loads in its natural environment. While it is presently unknown whether these viruses may have a detrimental effect on the health of this species, their phylogenetic relatedness with other Picornaviridae found in different molluscan species suggest they may be not the product of bioaccumulation by filter- feeding of food particles. Further ultrastructural investigations could clarify whether the presence of VLPs in the tissues of Congeria are linked with tissue lesions, helping to define whether alterations of the virome associated with this species should be considered as an additional risk factor to be monitored for the conservation of this endangered cave-dwelling bivalve. In light of the previous experience built in the conservation of freshwater bivalves living in surface waters (Brian et al., 2021), we believe that an improved characterization of the potentially pathogenic microorganisms associated with Congeria and its unique habitat, should be considered as a priority for a better planning of future translocation-based conservation efforts. Acknowledgements This work was supported by the Microgrants funding program of the University of Trieste - Transcriptomic investigation of the endemic cave- dwelling “living fossil” bivalve Congeria kusceri. References Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 215: 403-410, 1990. Barbosa Solomieu V, Renault T, Travers MA. Mass mortality in bivalves and the intricate case of the Pacific oyster, Crassostrea gigas. J Invertebr Pathol. 131: 2-10, 2015. Bilandžija H, Jalžić B. Dinaric cave clam, Congeria kusceri Bole, in: Ozimec R (ed), Red Book of Croatian Cave Dwelling Fauna. Ministarstvo kulture, Državni Zavod za zaštitu prirode, Zagreb, Croatia, pp. 67-68, 2009. Bilandžija H, Morton B, Podnar M, Cetković H. Evolutionary history of relict Congeria (Bivalvia: Dreissenidae): unearthing the subterranean biodiversity of the Dinaric Karst. Front Zool. 10: 5, 2013. Bilandžija H, Puljas S, Gerdol M. Hidden from Our Sight, but Not from Our Impact: The Conservation Issues of Cave Bivalves. Preprints. 2021050023, 2021. Blair CD, Olson KE, Bonizzoni M. The Widespread Occurrence and Potential Biological Roles of Endogenous Viral Elements in Insect Genomes. Curr. Issues Mol. Biol. 34: 13-30, 2020. Brian JI, Aldridge DC. Endosymbionts: An overlooked threat in the conservation of freshwater mussels? Biol. Conserv. 237: 155- 165, 2019. Brian JI, Ollard IS, Aldridge DC. Don't move a mussel? Parasite and disease risk in conservation action. Conserv. Lett. 14: e12799, 2021. Carballal MJ, Villalba A, Iglesias D, Hine PM. Virus- like particles associated with large foci of heavy hemocytic infiltration in cockles Cerastoderma edule from Galicia (NW Spain). J Invertebr Pathol. 84: 234-237, 2003. Cavanaugh JE. Unifying the derivations for the Akaike and corrected Akaike information criteria. Stat. Probab. Lett. 33: 201-208, 1997. Cheng Z, Yang J, Xia H, Qiu Y, Wang Z, Han Y, et al. The nonstructural protein 2C of a Picorna- like virus displays nucleic acid helix destabilizing activity that can be functionally separated from its ATPase activity. J Virol. 87: 5205-5218, 2013. Crespo-González C, Rodríguez-Domínguez H, Soto-Búa M, Segade P, Iglesias R, Arias- Fernández C, et al. Virus-like particles in Urastoma cyprinae, a turbellarian parasite of Mytilus galloprovincialis. Dis Aquat Organ. 79: 83-86, 2008. Cuttelod A, Seddon M, Neubert E. European Red List of Non-marine Molluscs. European Union, 2011. Darriba D, Posada D, Kozlov AM, Stamatakis A, Morel B, Flouri T. ModelTest-NG: a new and scalable tool for the selection of DNA and protein evolutionary models. bioRxiv 612903, 2019. Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 32: 1792-1797, 2004. Finn RD, Clements J, Eddy SR. HMMER web server: interactive sequence similarity searching. Nucleic Acids Res. 39: W29-W37, 2011. Finn RD, Mistry J, Tate J, Coggill P, Heger A, Pollington JE, et al. The Pfam protein families database. Nucleic Acids Res. 38: D211-D222, 2009. Goldberg TL, Dunn CD, Leis E, Waller DL. A Novel Picorna-Like Virus in a Wabash Pigtoe (Fusconaia flava) from the Upper Mississippi River, USA. Freshw. Mollusk Biol. Conserv. 22: 81-84, 2019. Hansman GS, Oka T, Li TC, Nishio O, Noda M, Takeda N. Detection of human enteric viruses in Japanese clams. J. Food Prot. 71: 1689- 1695, 2008. Hine PM, Wesney B. Virus-like particles associated with cytopathology in the digestive gland epithelium of scallops Pecten novaezelandiae and toheroa Paphies ventricosum. Dis. Aquat. Org. 29: 197-204, 1997. Huelsenbeck JP, Ronquist F. MRBAYES: Bayesian inference of phylogenetic trees. Bioinforma. Oxf. Engl. 17: 754-755, 2001. Ip HS, Desser SS. A picornavirus-like pathogen of Cotylogaster occidentalis (Trematoda: Aspidogastrea), an intestinal parasite of freshwater mollusks. J. Invertebr. Pathol. 43: 197-206, 1984. 35 Jones JB, Scotti PD, Bearing SC, Wesney B. Virus- like particles associated with marine mussel mortalities in New Zealand. Dis. Aquat. Organ. 25: 143-149, 1996. Jones P, Binns D, Chang HY, Fraser M, Li W, McAnulla C, et al. InterProScan 5: genome- scale protein function classification. Bioinforma. Oxf. Engl. 30: 1236-1240, 2014. Jovanović Glavaš O, Jalžić B, Bilandžija H. Population density, habitat dynamic and aerial survival of relict cave bivalves from genus Congeria in the Dinaric karst. Int. J. Speleol. 46: 13-22, 2016. Kal AJ, van Zonneveld AJ, Benes V, van den Berg M, Koerkamp MG, Albermann K, et al. Dynamics of gene expression revealed by comparison of serial analysis of gene expression transcript profiles from yeast grown on two different carbon sources. Mol. Biol. Cell. 10: 1859-1872, 1999. Lasa A, Cesare A, Tassistro G, Borello A, Gualdi S, Furones D, et al. Dynamics of the Pacific oyster pathobiota during mortality episodes in Europe assessed by 16S rRNA gene profiling and a new target enrichment next-generation sequencing strategy. Environ. Microbiol. 21: 4548-4562, 2019. Le SQ, Gascuel O. An improved general amino acid replacement matrix. Mol. Biol. Evol. 25: 1307- 1320, 2008. Li Q, Zheng Z, Liu Y, Zhang Z, Liu Q, Meng J, et al. 2C Proteins of Enteroviruses Suppress IKKβ Phosphorylation by Recruiting Protein Phosphatase 1. J. Virol. 90: 5141-5151, 2016. Li W, Godzik A. Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics 22: 1658- 1659, 2006. Love MI, Hogenesch JB, Irizarry RA. Modeling of RNA-seq fragment sequence bias reduces systematic errors in transcript abundance estimation. Nat. Biotechnol. 34: 1287-1291, 2016. Marszewska A, Cichy A. Unionid clams and the zebra mussels on their shells (Bivalvia: Unionidae, Dreissenidae) as hosts for trematodes in lakes of the Polish lowland. Folia Malacol. 23: 149-154, 2015. Moreira R, Balseiro P, Planas JV, Fuste B, Beltran S, Novoa B, et al. Transcriptomics of in vitro immune-stimulated hemocytes from the Manila clam Ruditapes philippinarum using high- throughput sequencing. PloS One 7: e35009, 2012. Morton B, Puljas S. Life-history strategy, with ctenidial and pallial larval brooding, of the troglodytic ‘living fossil’ Congeria kusceri (Bivalvia: Dreissenidae) from the subterranean Dinaric Alpine karst of Croatia. Biol. J. Linn. Soc. 108: 294-314, 2013. Müller T, Czarnoleski M, Labecka AM, Cichy A, Zając K, Dragosz-Kluska D. Factors affecting trematode infection rates in freshwater mussels. Hydrobiologia 742: 59-70, 2015. Ng TFF, Marine R, Wang C, Simmonds P, Kapusinszky B, Bodhidatta L, et al. High variety of known and new RNA and DNA viruses of diverse origins in untreated sewage. J. Virol. 86: 12161-12175, 2012. Nurk S, Meleshko D, Korobeynikov A, Pevzner PA. metaSPAdes: a new versatile metagenomic assembler. Genome Res. 27: 824-834, 2017. Péden R, Poupin P, Sohm B, Flayac J, Giambérini L, Klopp C, et al. Environmental transcriptomes of invasive dreissena, a model species in ecotoxicology and invasion biology. Sci. Data 6: 234, 2019. Puljas S, Peharda M, Morton B, Giljanović NŠ, Jurić I. Growth and Longevity of the “Living fossil” Congeria kusceri (Bivalvia: Dreissenidae) from the Subterranean Dinaric Karst of Croatia. Malacologia 57: 353-364, 2014. Rambaut A, Drummond AJ, Xie D, Baele G, Suchard MA. Posterior Summarization in Bayesian Phylogenetics Using Tracer 1.7. Syst. Biol. 67: 901-904, 2018. Rasmussen LPD. Virus-associated granulocytomas in the marine mussel, Mytilus edulis, from three sites in Denmark. J. Invertebr. Pathol. 48: 117- 123, 1986. Renault T. Picornalike Viruses of Mollusks. In: Kibenge FSB and Godoy MG (eds), Aquaculture Virology, Elsevier, Amsterdam, Netherlands, pp 529-531, 2016. Renault T, Novoa B. Viruses infecting bivalve molluscs. Aquat. Living Resour. 17: 397-409, 2004. Richard JC, Leis E, Dunn CD, Agbalog R, Waller D, Knowles S, et al. Mass mortality in freshwater mussels (Actinonaias pectorosa) in the Clinch River, USA, linked to a novel densovirus. Sci. Rep. 10: 14498, 2020. Rigonni H, Bilandžija H, Engel AS. Food web analysis for a stygobitic community in Croatian karst using stable isotopes. Proceedings of the 18th International Congress of Speleology, July, Savoie Mont Blanc, France, 2021. Rosani U, Gerdol M. A bioinformatics approach reveals seven nearly-complete RNA-virus genomes in bivalve RNA-seq data. Virus Res. 239: 33-42, 2017. Rosani U, Shapiro M, Venier P, Allam B. A Needle in A Haystack: Tracing Bivalve-Associated Viruses in High-Throughput Transcriptomic Data. Viruses 11: 205, 2019. Sasakova N, Gregova G, Takacova D, Mojzisova J, Papajova I, Venglovsky J, et al. Pollution of Surface and Ground Water by Sources Related to Agricultural Activities. Front. Sustain. Food Syst. 2: 42, 2018. Scapolatiello A. Transcriptome assembly and analysis of the endangered cave-dwelling freshwater bivalve Congeria kusceri. Master thesis, University of Trieste, pp. 97, 2021. Shi M, Lin XD, Tian JH, Chen LJ, Chen X, Li CX, et al. Redefining the invertebrate RNA virosphere. Nature 540: 539-543, 2016. 36 Shmakov SA, Wolf YI, Savitskaya E, Severinov KV, Koonin EV. Mapping CRISPR spaceromes reveals vast host-specific viromes of prokaryotes. Commun. Biol. 3: 1-9, 2020. Söding J, Biegert A, Lupas AN. The HHpred interactive server for protein homology detection and structure prediction. Nucleic Acids Res. 33: W244-248, 2005. Stepien CA, Morton B, Dabrowska KA, Guarnera RA, Radja T, Radja B. Genetic diversity and evolutionary relationships of the troglodytic ‘living fossil’ Congeria kusceri (Bivalvia: Dreissenidae). Mol. Ecol. 10: 1873-1879, 2001. Talavera G, Castresana J. Improvement of phylogenies after removing divergent and ambiguously aligned blocks from protein sequence alignments. Syst. Biol. 56: 564-577, 2007. Taskinen J. Cercarial production of the trematode Rhipidocotyle fennica in clams kept in the field. J. Parasitol. 84: 345-349, 1998. Vranješ M, Prskalo M, Džeba T. Hidrologija i hidrogeologija sliva Neretve i Trebišnjice, osvrt na izgradnju dijela he sustava-gornji horizonti. E-Zbonik Electron. Collect. Pap. Fac. Civ. Eng. 5: 1-23, 2013. Yang S, Shan T, Wang Y, Yang J, Chen X, Xiao Y, et al. Virome of riverside phytocommunity ecosystem of an ancient canal, Research Square 25620, 2021. Zupan Hajna N. Dinaric karst-Geography and geology, in: White, W.B., Culver, D.C., Pipan, T. (Eds.), Encyclopedia of Caves (Third Edition). Academic Press, pp. 353-362, 2019.