Transcriptome analysis of Catarina scallop (Argopecten ventricosus) juveniles treated with highly-diluted immunomodulatory compounds reveals activation of non-self-recognition system.

Marine bivalve hatchery productivity is continuously challenged by apparition and propagation of new diseases, mainly those related to vibriosis. Disinfectants and antibiotics are frequently overused to prevent pathogen presence, generating a potential negative impact on the environment. Recently, the use of highly diluted compounds with immunostimulant properties in marine organisms has been trailed successfully to activate the self-protection mechanisms of marine bivalves. Despite their potential as immunostimulants, little is known about their way of action. To understand their effect, a comparative transcriptomic analysis was performed with Argopecten ventricosus juveniles. The experimental design consisted of four treatments formulated from pathogenic Vibrio lysates at two dilutions: [(T1) Vibrio parahaemolyticus and Vibrio alginolyticus 1D; (T2) V. parahaemolyticus and V. alginolyticus 7C]; minerals [(T3) PhA+SiT 7C], scorpion venom [(T4) ViT 31C]; and one control (C1) hydro-alcoholic solution (ethanol 1%). The RNA sequencing (RNAseq) analysis showed a higher modulation of differentially expressed genes (DEG) in mantle tissue compared to gill tissue. The scallops that showed a higher number of DEG related to immune response in mantle tissue corresponded to T1 (V. parahaemolyticus and V. alginolyticus lysate) and T3 (Silicea terra® - Phosphoric acid®). The transcriptome analysis allowed understanding some interactions between A. ventricosus juveniles and highly-diluted treatments.


Introduction
Marine bivalve hatchery production has been challenged by apparition and propagation of new diseases [1], mainly those related with Vibrio spp. [2,3], leading to high mortality and economical losses. Recent studies have suggested that proliferation of pathogenic vibrio species will increase in a near future due to sea surface water warming [4,5], and most of the time, hatchery conditions have been favourable to Vibrio spp. [6]. Traditionally, hatcheries around the world have routinely implemented the prophylactic application of antibiotics and disinfectants to prevent proliferation or kill pathogenic bacteria, as standard prophylactic or therapeutic procedures [7,8].
Chemotherapeutic applications have been restricted by governmental institutions [2] for their undesirable effects, such as selection of antibiotic-resistant bacteria [7,8], their implications in dissolution of mercury compounds when contaminated water is in contact with marine sediment [9] and alterations in the beneficial gastrointestinal microbiota of cultured organisms [10].
Immunostimulants have been proposed as promising and sustainable alternatives to avoid massive bivalve mortalities in hatchery and reduce or replace the use of antibiotics [11]. Some immunomodulatory compounds can prophylactically activate the organisms' self-protection instead of killing the pathogen by chemical compounds, which may compromise the organisms' health and environment. Within the immunostimulant category, highly-diluted immunomodulatory compounds (HDIC) formulated from bacterial lysates, scorpion venom, silica and phosphoric acid have been recently suggested as eco-friendly and sustainable options to activate the immune system without affecting the general condition index of marine organisms [12]. HDIC are inexpensive to produce because serial dilutions at decimal (D, 1:10) or centesimal (C, 1:100) scale [13] are used, and the action mode of these treatments depends on their dilution grade. Low decimal dilution (D) effects are attributed to diluted molecules and nanoparticles while high centesimal dilution (C) effects are attributed to nanoparticles and electromagnetic fields [14], which are sensed by traditional molecular receptors or photoelectrochemical sensing system activated by ultra-weak photo emission signals [15]. These highly-diluted molecules work as low intensity danger signals that generate endogenous amplification by hormesis effect, time-dependent sensitisation and stochastic resonance process [16,17].
Low concentrations of attenuated bacteria [18,19] have already proven to activate immune response in different organisms and highly-diluted concentrations of bacterial lysates, Silicea terra1 (Similia1, CDMX, MX), Phosphoric acid1 (Similia1, CDMX, MX) and Vidatox1 (Labiofam1, Habana, Cuba) have been successfully used in marine organisms to improve response against pathogens [12]. Previous studies have demonstrated that survival of Argopecten ventricosus juveniles increased when the organisms were previously therapeutically treated by Silicea terra1 and Phosphoric acid1 31C, and then challenged against highly pathogenic bacteria, such as Vibrio parahaemolyticus [20]. Furthermore, increase in growth, enzymatic antioxidant activity, and haemocyte proliferation were reported when juvenile scallops were treated prophylactically with diluted and highly-diluted pathogenic bacterial lysates (1D and 7C) for a 21-day period [21]. Although those immunomodulatory compounds are known to allow organisms to activate immune response and increase survival when they face pathogen infections, little is known about how they work on the immune response of marine organisms.
Transcriptomic (RNA-seq) analysis allows discovering potential and novel action mechanisms. Particularly, these analyses have provided information to analyse complex relations among disease, drugs, and organisms [22]. Moreover, RNA-seq has been used in marine bivalves to understand physiological stress [23,24], toxicological effects [25,26], and resistance to disease and immunology [27,28]. In this sense, the RNA-seq analysis helps to understand the molecular response activated by the effect of time-dependent diluted immunomodulatory treatments, which allowed identifying the down-and up-regulated genes of juvenile scallops in response to HDIC.
Because the immune defense mechanisms in response to HDIC formulated by bacterial lysates, Silicea terra1, Phosphoric acid1 and Vidatox1 remain unknown, this study performed a comparative analysis of the mantle and gill transcriptome profile treated with those HDIC. These data will provide important information about the genes and mechanisms that are being regulated in A. ventricosus scallop and contribute to understanding how HDIC acts on marine organisms' response.
This study also provided for the first time whole transcriptome data of the Catarina scallop A. ventricosus juveniles, which was selected as a model organism because of its importance as fishing resource in Baja California Sur, Mexico. Its wild populations have been decreasing throughout the years [29], and spat production at hatchery level faces high mortality events as it is a species highly susceptible to vibrosis [30], one of the main bacteria present in hatchery culture water.

Scallop acquisition and experimental design
The non-governmental association Noroeste Sustentable (NOS) provided 1500 juvenile scallops (average length 1.98 ± 0.1 cm) from Bahia de La Paz, Mexico to perform the experiment. Scallops were acclimated for one week in a nursery upwelling recirculating system with constant food concentration of 150 000 cell mL -1 (Isochrysis galbana; Chaetoceros calcitrans; cellular proportion 1:1). Filtered seawater (1 μm, activated carbon and ultraviolet (UV) irradiation) at 24˚C and 38.5 ± 0.5 UPS salinity were continuously supplied to allow water change totally every day.
After acclimatization, scallops were transferred to 15 experimental units (36-L container) with 52 scallops each (three experimental units per treatment). During the experimental period (21 days), scallops were kept in an open-recirculating flow system that provided filtered and treated seawater (1 μm, activated carbon and UV irradiation) with a blend of I. galbana and C. calcitrans microalgae (199 607 794 cell/organism/day) at 23.5 ± 0.5˚C, and 38.5 ± 0.5 UPS salinity. A 21-day experimental period was sufficient to strengthen the immune system of this species with HDIC-based treatments [20].
At the end of the experiment, organisms for each experimental condition were taken. The mantle and gill tissues of three sampled scallops of each experimental replicate unit were excised, separately fixed in RNAlater solution (#AM7020, ThermoFisher, Scientific, Waltham, MA, USA) and stored at -80˚C for analysis. Mantle and gill tissue were selected as target organs because they have been highly implicated with immune response and have been related with pattern recognition receptors in marine bivalves [31,32].

RNA extraction and cDNA synthesis
To extract RNA, the tissue (mantle and gill) samples from three scallops (~100 mg) were collected and pooled. Each pool was considered as the biological sample. For the Real-Time quantitative polymerase chain reaction (RT-qPCR) analysis, three pools of each tissue were used per treatment ( Fig 1A). For transcriptomic analysis, two pools (~100 mg each) were used per treatment ( Fig 1A). The RNA extraction was assessed following the methodology previously described for the scallop Nodipecten subnodosus RNA extraction [35], using the TriPure reagent (Roche Diagnostics, Indianapolis, IN, USA), followed by ethanol/chloroform purification and DNAse cleaning. Samples were determined for RNA concentration using a NanoDrop 2000/ 2000c1 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), and electrophoresis was performed to evaluate RNA quality. Good quality RNA from a total of 20 Pools (two per treatment per tissue) were sent to the Laboratorio Nacional de Genómica para la Biodiversidad (LANGEBIO, Laboratory of Genomics Services, CINVESTAV, Campus Irapuato, GTO, MX) to create and sequence the RNA-seq libraries using the Illumina Next-Seq1 (San Diego, CA, USA). For RT-qPCR analysis, the cDNA synthesis was performed following the methodology previously reported to quantify gene relative expression levels in N. subnodosus [35].

Illumina sequencing
For RNA-seq library preparation, the two RNA pools were selected from mantle and gill of each experimental condition. A total of 20 libraries were prepared using the RNA extractions. In addition, 12 libraries were prepared with the organisms treated with HDIC formulated by sodium metasilicate-Phosphoric acid 1D plus two controls (C2: diluted ethanol 1:100 and C3: nothing added) to increase assembly and annotation quality. Libraries were validated by 2100 Bioanalyzer Instrument (Agilent, Santa Clara, CA, USA) and sequenced with Illumina NextSeq1 (San Diego, CA, USA) platform on mode paired-end (2 X 150 high). Generation of libraries and sequencing were performed by LANGEBIO (Laboratory of genomics services, CINVESTAV, Campus Irapuato, GTO, MX).

Transcriptome, annotation and statistical analyses
The raw reads obtained by the Illumina NextSeq1 platform (San Diego, CA, USA) were first analysed with the FastQC software [36] to determine quality sequences; then, the reads were filtered using the Trimmomatic program [37] by removing low-quality reads (Q < 25), ambiguity ("N"), reads under 50 bp and trimming the adapters. The high-quality filtered reads were used to perform the de novo transcriptome assembly with Trinity software [38] with the "non-normalized reads" option. This process was performed in a WorkStation with 10 cores and 512 GB RAM. Once the assembly process was completed, the quality and integrity of the transcriptome were analysed. The high-quality filtered reads were mapped to assembled transcripts with RSEM software which used Bowtie2 for the alignment [39], and the expression values were normalised to transcripts per million (TPM) for each transcript. Those with low expression values (TPM < 0.01) were removed. The assembled transcripts were compared against the non-redundant nucleotide and protein databases of the National Center for Biotechnology Information (NCBI) using BLAST with an e-value < 1e -5 , sequences with hit to viruses, bacteria, fungus, and so on (considered contaminant sequences) were removed using the SeqClean tool [ To determine the functional annotation based on homology, the filtered transcriptome was compared against the public protein non-redundant nucleotides (NCBI), Swiss-Prot and RefSeq databases using BLASTX (e-value < 1e -5 ). In addition, a BLASTX (e-value < 1e-5 ) against predicted proteins of the Crassosrea gigas, Crassostrea virginica, Mytilus galloprovincialis and Mizuhopecten yessoensis genomes was performed. The functional annotation was assigned according to the BLAST hits. Finally, the functional categorisation was designated based on the comparison with the Gene Ontology (GO), InterPro and Kyoto Encyclopaedia of Genes (KEGG) databases using the BLAST2GO software [46]. The raw reads were deposited in the Gene Expression Omnibus of NCBI under accession number PRJNA596225.

Differential expression and enrichment analyses
Only were 20 libraries corresponding to the treatments (T1, T2, T3, T4) and control (C1) considered for the posterior analysis. High-quality alignment of the filtered reads to the filtered transcriptome and quantification of reads were carried out using the RSEM software [39]. Raw counts of each transcript of each library were used to generate a matrix with all the experimental conditions and replicates. These raw counts were normalised to fragments per kilobase per million mapped reads (FPKM). To determine dispersion between samples and biological replicates, the Pearson correlation value was calculated between all experimental libraries. Only were libraries with Pearson's correlation > 0.88 between replicates kept for subsequent analyses.
Once every biological replicate was validated, transcripts lowly expressed were removed and only transcripts with a FPKM > 0.3 in at least two replicates of at least one treatment were kept, as edgeR reduce analysis quality when very low expressed genes are use [47]. Differentially expressed genes in Catarina scallop juveniles treated with the HDIC (Fig 1B) versus control without treatment were identified with the edgeR package [47], using a statistical method based on the generalised linear model (GLM). To estimate the variance between samples, tagwise dispersion was calculated. Transcripts with FDR < 0.01 and Log 2 FoldChange >1 (Log 2 FC) were considered differentially expressed. DEG transcripts were grouped by tissue as down-and up-regulated for each treatment, and then Venn diagrams were generated including all the comparisons between DEG of treatments.
Gene Ontology and KEGG category enrichment analyses were performed using BLAS-T2GO [46] and KOBAS [48], respectively, with FDR < 0.05. The function to the most specific terms was used to reduce the result-set of over-represented GO terms. Heatmaps and hierarchical clustering were performed using the gplots and hclust in R [49].

RT-qPCR validation
To validate differential expression from RNAseq data, 10 target genes (Log 2 FC > 1.5 and FDR < 0.01) and four constitutive genes (CV% < 0.3, Log 2 FC between -1.5 and +1.5 and FPKM > 0.3) from the transcriptome were selected (S1 Table). Primers were designed using the software Primer3 [50] and then analysed with OligoAnalyzer Tool [51]. The most stable of the constitutive genes for each tissue was selected using RefFinder software [52] to express relative expression. Primers were evaluated and target sequences validated by sequencing analysis (Macrogen, Seoul, Korea).
For the RT-qPCR analysis, three biological replicates were assessed per treatment (each biological replicate was performed by triplicate) using previously synthesised cDNA from mantle and gill. RT-qPCR analysis was performed as previously reported to quantify gene relative expression levels in N. subnodosus [35] following Minimum Information for Publication of Quantitative Real-Time PCR Experiments (MIQE) guidelines [53]. Percentage of efficiency of each primer pair was obtained from efficiency curves from six serial dilutions (1:5) [53]. Reactions of RT-qPCR were performed with the equipment using 2X Eva-Green1 (ThermoFisher, Scientific, Waltham, MA, USA) PCR mix. Reaction mixtures were made in 15 μL, including 10 μL mix (0.45 U of GoTaq Flexi DNA polymerase (Promega, Madison, WI, USA), 2.5 mM MgCl2, 1x Go Taq Flexi Buffer, 0.2 mM dNTP Mix (Promega, Madison, WI, USA), 2x EvaGreen fluorescent dye (Biotium, Inc. Fremont, CA, USA), 0.15-0.45 μM each primer) and 5μL cDNA (diluted to 50 ng μL -1 ). The amplification parameters were 4 min at 94˚C, followed by 39 cycles of 9 s at 95˚C, 30 s at 60˚C and 30 s at 72˚C; finally, a melt curve was added to confirm amplification of specific products. Data were analysed using Ct values and the 2 −ΔΔCt method [54], and relative expression was normalised to the abundance of sodium/potassium-transporting ATPase subunit alpha-like gene for mantle tissue and cAMP-dependent protein kinase catalytic subunit-like isoform X3 gene for gill tissue.
According with the 2 −ΔΔCt values, data for each tissue and condition selected were converted to natural logarithm and analysed using a student T-test (independent). Statistical comparison between the treatment and its respective control were performed, and statistical significances were set at p < 0.05. The log 2 FC was assessed using treatment vs control data to validate the expression of selected genes compared to the transcriptome results and by a linear regression. The software Statistica1 10.0 (StatSoft, Tulsa, OK, USA) was used to perform statistical analysis.

Ethics statements
All the scallops used in this work were handled in accordance with the Official Mexican Standard protocols (NOM-062-ZOO-1999). Argopecten ventricosus juveniles were provided from the non-governmental association Noroeste Sustentable (NOS) and transported to aquarium tanks at CIBNOR with all the required permits from the federal agency CONA-PESCA. The animals used in this study were produced by Acuacultura Robles hatchery in captivity for experimental purposes. Organisms were kept in optimal culture conditions to avoid stressful conditions and no harmful effect has been detected using HDIC in marine organisms.
According to the Internal Committee for the Care and Use of Laboratory Animals (CIC-UAL) recommendations, the number of sampled organisms contemplated "the rule of maximizing information published and minimizing unnecessary studies". In this sense, for this work, we considered the minimum number of organisms needed to obtain a high-quality transcriptome.

De novo assembly and functional annotation of the juvenile scallop Argopecten ventricosus transcriptome
To evaluate the effect of the immunomodulatory compounds on gene regulation, 20 libraries from juvenile Catarina scallops were considered and treated with the four HDIC (Fig 1B) versus control C1 (two libraries for experimental condition of each tissue) plus 12 libraries using NexSeq illumina in paired-end mode (2 X 150). A total of 423 million raw pair reads were generated. After the process of filtering reads, 372 million high-quality pair reads were used for the de novo transcriptome assembly. The final assembly generated 191 432 transcripts with a maximum length of 36 697 bp, minimum length of 184 bp, an average length of 863 bp and N50 value of 1123 bp (Fig 2A). The transcript size distribution showed a higher abundance of transcripts from 100 to 600 bp and lower abundance of transcripts from 2000 to 2200 bp ( Fig  2B). The analysis of transcriptome completeness using the BUSCO tool against the database of Metazoa odb9, which consisted of benchmarking universal single-copy of conserved orthologous genes, showed that the percentage of complete and fragmented genes was 79.4% and 17%, respectively ( Fig 2D).
Functional annotation results showed that 62 548 (32%) of the transcripts were annotated with at least one database. The Mizuphecten yesseniosis genome was the database that allowed the highest number of annotated transcripts (60157 hits). M. galloprovincialis genome and KEGG database showed the lowest number of annotated transcripts (15 425 and 17 487, respectively) (Fig 2A). The distribution of the BLAST top hits against the non-redundant protein (NCBI) database showed that the species with the highest number of top hits were the oyster Crassostrea gigas (18 339 hits) and the pectinid Mizuphecten yesseniosis (12 104 hits) ( Fig 2C).
As a result, 47 313 transcripts were annotated from the functional annotation of the A. ventricosus juvenile transcriptome in Blast2go. A total of 29 116 transcripts were assigned to some biological processes, and the most represented were related to biosynthesis (2749), signal transduction (2 711) and cellular nitrogen compound metabolic process (2485) (Fig 3A). For molecular functions 30 274 transcripts were annotated. The most represented molecular function categories were related to ion binding (6662), oxidoreductase activity (2233) and transmembrane transporter activity (2076) (Fig 3A). In the category of cellular components 19 886 transcripts were recorded. The cellular component with higher number of transcripts corresponded to protein containing complex (3305), cytoplasm (1329) and intracellular (1246) components ( Fig 3A). The KEGG annotation mapped 17 487 transcripts and showed the metabolic pathways mostly related to biosynthesis of secondary metabolites (222), thermogenesis (139) and endocytosis (128) (Fig 3B).
The Venn diagram showed that juvenile scallops had a specific response to each treatment as most of the DEG were not shared between treatments. In most of the cases, the mantle recorded a higher DEG number by the treatment effect compared to gill. The scallops that recorded the highest number of up-and down-regulated DEG in mantle were from T1/C1 (1841 up, 491 down) while the lowest number of up-regulated DEG were recorded in scallops from T2/C1 (257 up, 275 down) (Fig 4A). The DEG profiles in mantle were clustered, showing that T3 and T4 were clustered in the same group and T1 and T2 in other groups separately ( Fig 4B). In gill, treatments with the highest number of up-and down-regulated DEG were recorded in scallops from T2/C1 (248 up, 309 down) and T3/C1 (235 up, 483 down). The lowest number of up-and down-regulated DEG in gills were recorded in T4/C1 (94 up, 245 down) (Fig 4C). T1 and T2, formulated by V. parahaemolyticus and V. alginolyticus lysate in different concentration, were clustered in the same group. T3 (Silicea terra1 and Phosphoric acid1) and T4 (Vidatox1, formulated by commercial brands, were clustered in a different group (Fig 4D).

Co-regulation of gene expression and functional enrichment analysis
Hierarchical clustering analysis of the differentially expressed transcripts allowed identifying groups of co-expressed transcripts for each tissue and treatment. Furthermore, the GO biological processes and metabolic pathways enriched for each group were analysed (Figs 5 and 6). In mantle, the clustering analyses showed 20 clusters; from these clusters seven were enriched by GO terms (p-value < 0.01). Interestingly, three groups were related to oxidative phosphorylation (2,12,14), and four groups (cluster 13, 15, 16 and 20) were enriched in biological processes related to non-self-recognition and immune response processes. Cluster 13 comprised five transcripts related to focal adhesion and five related to endocytosis (Fig 5). Most of the DEG from cluster 13 in mantle were down-regulated mainly by the effect of T4/C1 (mean Log 2 FC -3.7) (S1 Fig). A total of five transcripts from cluster 15 were related to apoptosis process. DEG from cluster 15 were mainly up-regulated in all treatments with a higher number of up-regulated DEG in T3/C1 (mean Log 2 FC 5.8). In cluster 16, 44 transcripts were related to pathogenic Escherichia coli infection pathway, 31 to cytoskeleton organization and eight to immune response. All treatments from cluster 16 in mantle were mainly up-regulated with the highest number of up-regulated DEG in T1/C1 (mean Log 2 FC 2.5) (S1 Fig). Cluster 20 comprised 65 transcripts related to signal transduction and 18 to pathogenic E. coli infection pathway. Cluster 20 showed up-and down-regulated DEG with the highest number of up-regulated DEG in T4/C1 (mean Log 2 FC 0.7) and the highest number of down-regulated DEG in T1/C1 (mean Log 2 FC -0.5) (Figs 5 and S1). In mantle the enrichment of clustered co-expressed genes were in accordance to the enrichment of all DEG as shown in (S2 and S3 Figs).
On the other hand, 17 groups were generated in gill, of which six showed enriched GO terms (p-value < 0.01, Fig 6). The enrichment analysis showed that three groups (6, 7 and 12) were related to oxidative phosphorylation and three (clusters 5, 8 and 16) were related to nonself-recognition and immune response (Fig 6). In cluster 5 from gill tissue, 11 transcripts were related to cellular response to stimulus (Fig 6). Scallops from cluster 5 showed a DEG downregulation mainly by effect of T1/C1 (mean Log 2 FC -7) (S1 Fig). Cluster 8 comprised 10 transcripts related to NF-kappa B signalling pathway, 44 to signal transduction and 17 to stress response. In cluster 8 the DEG were mainly down-regulated in all treatments with a higher number of down-regulated DEG in T2/C1 (mean Log 2 FC -1). A total of 13 transcripts from cluster 16 were related to transmembrane transport in gill tissue. Transcripts for all the treatments in gill were mainly up-regulated in cluster 16 with the highest number of up-regulated DEG in T1/C1 (mean Log 2 FC 2.5) (Figs 6 and S1). As recorded in mantle, the enrichment of clustered co-expressed genes from gill tissue were in accordance to the enrichment of all DEG shown in (S2 and S3 Figs).

Identification of genes associated to immune response regulated by effect of highly-diluted immunomodulatory compounds (HDIC)
Since HDIC have shown to improve the ability of these organisms to cope with pathogen infections [20], probably through immune system modulation [12,21], the transcripts related to the immune system were identified and their expression levels analysed in response to the treatments. Based on the functional annotation and GO and KEGG categorisation, a total of 193 transcripts were identified related to the non-self-recognition receptor, internalisation and immune system. The expression analyses showed that the tissue with higher DEG number related to immune response was mantle (175 transcripts) compared to gill (28 transcripts) (S3 and S4 Tables). The representative metabolic pathways with higher DEG numbers were phagosome, endocytosis and NF-kappa B signal, but non-self-recognition receptors were also found.
In gill, most of the pathways were down-regulated; the treatment with the highest DEG was T2/C1 (up: 1, down: 15), and the one with the lowest DEG was T4/C1 (up: 1, down: 3). As in mantle, scallops from T2/C1 showed a higher number of down-regulated DEG (15) in gill compared with the up-regulated DEG (1). The metabolic pathways with higher numbers of downregulated DEG in T2/C1 were recorded in the NF-kappa B signal pathway (4) and endocytosis pathway (up: 1, down: 2). The gill of the organisms treated with T2/C1 up-regulated only one gene, and it was related to endocytosis. The metabolic pathways with higher numbers of up-regulated DEG in T4/C1 were observed in NOD-like receptor signalling (1) (Fig 7 and S4 Table).

Validation of RNA-seq data by qPCR analysis
Ten relevant DEG related to the immune response and regulated by the four treatments in this study were selected to validate the transcriptome data by RT-qPCR analysis. All contrasts were made using the control group as a reference; relative expression was normalised to the abundance of sodium/potassium-transporting ATPase subunit alpha-like gene in mantle tissue and cAMP-dependent protein kinase catalytic subunit-like isoform X3 gene in gill tissue because they were the most stable genes among the experimental conditions. The up-and down-regulated genes selected from the transcriptome had the same pattern in RT-qPCR analysis. IAP, HSP90; ERIS genes significantly (p < 0.001) up-regulated in organism mantle by effect of T1. In T2 gene HSP90 (p < 0.001) from the mantle significantly up-regulated and FILA (p < 0.01) down-regulated. The genes ERIS, TRAF3 (p < 0.01) and SRC (p < 0.001) up-regulated in the organism mantle from T3, and in T4 they promoted up-regulation of HEXA gene (p < 0.001) (Fig 8A). In gill T1 up-regulated (p < 0.05) the expression of TLR2 gene; T2 promoted the up- https://doi.org/10.1371/journal.pone.0233064.g008 regulation (p < 0.001) of TRAF3 gene; T3 regulated the expression of TLR3 gene (p < 0.01) (Fig 8A). Finally, the correlation between the expression of the genes in mantle and gill using Log 2 FC data from RT-qPCR and RNA-seq analyses was 0.80 and 0.94, respectively (Fig 8B and  8C), which were considered good.

Discussion
Most of the transcriptomes assembled in marine bivalves did not report fragmentation percentage in their results; however, fragmentation is known to be typical in de novo transcriptome assemblies [55]. The assembly statistics in this study was according to most published marine bivalve transcriptomes referred by the authors as high quality assembly [56,57,58] and showed a higher quality assembly compared to others [27,59,60].
The effectiveness of the use of HDIC in bivalves to improve their survival against pathogen infection challenges has been widely demonstrated, our work team has proven that the application of HDIC formulated with scorpion venom, Phosphoric acid1, Silicea terra1, sodium metasilicate and Vibrio lysate (V. alginoliticus-V. parahemolityics) strengthened marine organisms (bivalves, fish and shrimp) self-defense, allowing them to survive infection challenges, reduce parasite proliferation, increase energy reserves and growth; even some of them outperformed antibiotics when they were used prophylactically in A. ventricosus culture [20,21,61,62].
The efficiency of these HDIC to protect marine organisms against infections has been attributed to increases of the antioxidant system response, haemocyte count, energetic reserves and enzymes related to food assimilation [20,21,63]. Furthermore, attenuated V. splendidus cells activated the immune response in the scallop Argopecten purpuratus [18] and micro-sized silica [64], phosphorus [65] and toxins [66] activating antioxidant response, which allowed organisms to face pathogen infections by promoting free radical eradication. This information considered that the evaluated HDICs may have the potential to modulate immune response at high and low concentrations; however, the mechanisms regulated at low dilutions have not been elucidated yet. Information about how genes and mechanisms are being regulated in marine organisms treated with HDIC is scarce. In this sense, interestingly, the results from the transcriptomic analysis in this study suggested that the immune response was activated by most of the HDIC but not like conventional immunostimulants, such as polysaccharides, nutrients, oligosaccharides, herbs, antibacterial peptides and microorganisms [11] used at higher concentrations, which increased production of bactericidal and cytotoxicity activities, lysozyme and antimicrobial peptides [11,67,68]. In this study, HDIC regulated transcripts associated to recognition of non-self-molecules and internalisation of particles suggested that HDIC were not able to trigger all the immune response mechanisms due their high dilution. This result may explain those [61] in organisms of S. rivoliana treated with HDIC where only the genes IL-1β and MyD88 up-regulate before the challenge with V. parahaemolyticus. The results in this study suggested that the non-self-recognition system activation by effect of HDIC improved the response against infections. When the defense mechanisms were activated before the infection process, the organisms had more probabilities to successfully overcome a disease, which was the principle of the immunostimulants used in aquaculture [11].
It is worth to highlight that most of the DEG related to non-self-recognition, internalisation and immune response were detected in mantle compared to gill tissue. This result can be explained because mantle and mucosal epithelial cells are one of the first barriers that are in contact with the environment, and epithelial cells have been recognised as the first line of defense against organic, inorganic and pathogen intruders [69]. Epithelial cells are able to endocytose biotic and abiotic particles and activate signals to proliferate haemocytes, the principal effector of the immune response [70], which in turn can also eliminate intruders by endocytosis and enhance immune response [70]. Both epithelial cells and haemocytes can detect intruders by many pattern recognition receptors (PRR) that allow organisms to detect potential danger signals [69]. Thus, the up-regulated transcripts related to non-self-recognition, internalisation and immune response in mantle may be due to the participation of epithelial mantle cell, migrating haemocytes or both. Future research should clarify and corroborate which cells are interacting and recognising the HDIC.
The recognition mechanisms of marine bivalves that detect non-self-molecules are integrated by the PRR, which are activated by pathogen-associated molecular patterns (PAMPs), endogenous ligands (e.g. HSP70) and damage associated molecular patterns (DAMPs). The activation of the PRR allow organisms to activate the immune response depending on the stimuli received [59]. In agreement with the previous information, the results in this study showed that the up-regulated transcripts are related to PRR from the toll-like receptors family, which are highly conserved in the animal kingdom and have proven to activate the immune response in marine bivalves [59]. Specifically, this study detected the induction of TLR2 by the effect of T1 and T2, and TLR4 in T2 and T3. Moreover, T3 up-regulated TLR3. The activation of TLR2 and TLR4 in T1 and T2 was attributed to their formulation (V. alginoliticus and V. parahemolityics) because TLR2 and TLR4 are Pattern recognition receptors (PRR) located in plasma membrane and mainly recognise bacterial PAMPs, such as lipopolysaccharides [71] that can be found in gram negative bacteria like V. alginoliticus and V. parahemolityics [72]. In bibliography no data was found between the interaction of TLR3 and Silicea terra1 or Phosphoric acid1, which were the components in T3, but we only know that TLR3 is a cytosolic receptor that detects non-self-nucleic acids [73]. In this sense, this study has reported for the first time that TLR3 is associated to the exposure of highly-diluted Silicea terra1 or Phosphoric acid1. Additionally, these results showed that T3 formulated by Silicea terra1 and Phosphoric acid1 promoted the induction of the TLR4 gene. Accordingly, evidence that SiO 2 nanoparticles are able to increase the TLR4a gene expression has been observed in Dino zebrafish embryos [74], and S. rivoliana juveniles treated with Silicea terra1 -Phosphoric acid1 up-regulated MyD88 transcript only when organism were challenged. It should be noted that MyD88 regulation is highly related to TLR4 [75], which suggested that the up-regulation of TLR4 in T3 was associated to Silicea terra1 action. In T4, formulated with Vidatox1, no PRR regulated were observed, which may explain the low regulation of transcripts related to immune response. The information above suggested that HDIC had been detected mainly through toll-like receptors in treatments T1, T2 and T3, while Vidatox1 may have detected them by an alternative route in T4.
When small intruders (e.g. Bacteria, LPS) are detected via PRR, many mechanisms can be activated to respond against the potential danger signals, such as PAMPs and DAMPs [70]. The mechanisms activated in this study were dependent on each treatment and their dilutions. Treatments T1 and T2 mainly modulated the expression of transcripts related to endocytosis, lysosome, phagosome, proteasome and NF-kappa B pathways. Interestingly, T1, which had a higher concentration of bacterial compounds than T2, up-regulated transcripts associated to those processes while T2 regulated them negatively. Other authors have reported that the presence of PAMPs activated endocytosis, lysosome, and phagosome mechanisms [69,70] to destroy potentially harmful molecules or disease-causing microorganisms. The activation of those mechanisms can also activate the signalling cascade to express immune related genes [69,70] via NF-kappa B pathway signal [76]. This result suggested that the organisms treated with T1 recognised and responded to the presence of PAMPs or DAMPs in the HDIC formulated by V. parahaemolyticus and V. alginolyticus lysate.
In this study, most of the transcripts that up-regulated in NF-kappa B pathway by T1 were modulators at the first level of the immune response and effectors of the anti-apoptosis process. In this sense, T1 allowed activating a conservative response. The results showed that T1 also up-regulated genes related to sense microbial viability (STING), which allowed organisms to detect living bacteria [77], bind Vibrio protein (CALR) that can be released from the cell and neutralise vibrio [78,79] and heat shock proteins (HSP70 and HSP90), related to showing antigen (adaptive immune system) [80], protein fold [81], activating Toll-like receptors [82] and mediating immune response in bacterial challenges [83]. The above suggested that T1 allowed the activation of pathways related to recognition, neutralisation, internalisation and destruction of Vibrio. Although T2 did not activate transcripts associated to mechanisms of endocytosis, lysosome, phagosome, proteasome and NF-kappa B pathways, it has been successfully used in strengthening marine organisms against pathogens [12,13,61]. These results showed that the induction of TLRs in scallops by T2 may have been enough to protect organisms against pathogens. Because T2 was highly diluted compared to T1, the concentration of molecules from the V. parahaemolyticus and V. alginolyticus lysate might have not been enough to activate the internalisation and destruction of non-self-molecule mechanisms; this situation allowed organisms to tolerate the presence of these molecules and increase the expression of the transcripts related to detection mechanisms of potential danger signals, such as TLRs, which may imply a lower energetic cost for the organisms. Perhaps, as previously reported, PAMPs used at high-doses induced tolerance while at lower ones they allowed trained sensitisation, and in very low doses, they did not activate immune response [84]. Higher dilutions from the HDIC used in this study were in accordance with literature [84], which proposed that very low doses of PAMPs did not activate innate immune response. However, for the first time in Argopecten ventricosus this study detected the up-regulation of transcripts related to TLRs without activating the immune response. Contrary to very low doses, the low doses of HDIC in this study showed an activation of detection, neutralisation and internalisation responses against intruders without activating all the immune response mechanisms. The above suggests that a prolonged exposure time to low PAMPs doses may also generate some tolerance because the activation of immune response implies high ATP amounts. The reason why innate cells depend on PAMPs dose and exposure time may be explained by the energetic cost; if organisms sense a very low or low concentrations of PAMPs, the energetic balance is maintained; the organisms can invest energy in processes related to self-protection depending on the dose and time of exposure to the stressor (e.g. LPS).
Interestingly, organisms treated with T3 (Silicea terra1 -Phosphoric acid1) induced transcripts related to bacterial challenge (STING, CALR), antioxidant system (CAT, PRDX6) [59,85], cell migration (CTTN) [86], amplification of immune response (CD40LG, found in adaptive response) [87], and the receptor associated to actin-assembly machinery on the cytoplasmic side of the phagosome (P2X7) [88]. In previous publications, Silicea terra1 and Phosphoric acid1 had recorded an improvement of A. ventricosus juvenile growth, antioxidant activity [21] and survival when organisms were challenged against pathogenic bacteria [20]. Moreover, silicic acid, which has been reported in silica solutions [89], may activate catalase (CAT) and superoxide dismutase (SOD) activities in mouse brain from organisms intoxicated with aluminium [90]; in S. rivoliana, the use of Silicea terra1 and Phosphoric acid1 increased survival when organisms were challenged against pathogenic bacteria [61]. These results suggested that the activation of non-self-recognition, endocytosis and antioxidant system by T3 was the reason for survival when organisms were challenged.
As mentioned before, T4 behaved differently compared to the rest of the treatments and mainly activated the melanogenesis and haematopoietic processes and a gene associated to calcium transport (TRPM2), which have been linked to immune response when organisms are challenged against bacteria [70,71]; however, the immunomodulatory effect of Vidatox1 was not clear. Vidatox1 has demonstrated to increase SOD activity and survival of L. vannamei when it has been challenged against pathogenic bacteria [91], but it has also been related to proliferation of hepatocellular carcinoma in cultured mouse cells [92]. As in the other treatments, oxidative phosphorylation in mitochondria was modulated by the effect of T4. The activation of the immune response may be modulated by alternative routes as inflammasomes that are activated by DAMPs [71] and highly related to the production of mitochondria ROS (mROS), which have been implied in the activation of the immune response [71,88]. This result could explain the hepatocellular carcinoma proliferation in cultured mouse cell when Vidatox1 was used since proliferation of mROS may cause oxidation damage in cancer cells [93]; the positive effect recorded in challenged shrimp, as mROS, may also have been related to the activation of immune response and killing pathogens [71]. Thus, T4 may be more related to the activation of immune response by mitochondria signalling. Interestingly, T4 was formulated by the most diluted compound, suggesting that molecules highly diluted mainly act by their magnetic charge [14], sensed by traditional molecular receptors or photoelectrochemical sensing system activated by ultra-weak photo emission signals [15]. In Arabidopsis thaliana their defense mechanisms may be activated by the sound of caterpillars (Pieris rapae) eating their leaves [94] via mitochondria signalling; their mROS production has been implied in vibration recognition [95]. Since mitochondria have been related to the activation of the immune response by vibration recognition mechanisms, future research should analyse the effect of Vidatox1, taking into account the importance of mitochondria and oxidative phosphorylation process.
It is worth to notice that all the treatments modulated the oxidative phosphorylation process in mitochondria, an organelle that is not only related to vibration recognition but also plays a crucial role in supplying ATP since it is necessary for cell internalisation process, acidification or the phagosome, production of mROS for eradication of bacteria and activation of immune signalling cascades [71,88]. Nonetheless, the up-or down-regulated process could not be detected due to transcriptome fragmentation (17%), which diminished the robustness of expression estimates in fragmented genes [55].
This study used different HDIC, most of them showing up-regulation in the PAMPs recognition system, which may explain why marine organisms increased survival challenges against pathogens in previous studies [12].

Conclusion
This study opens new insight into the activation of Argopecten ventricosus self-protection mechanism using HDIC, which allowed answering questions and opening the possibility of generating others about interactions between organisms and their exposure to low toxin concentration, contaminants, minerals, pathogens and other stress factors during a prolonged time. The results of this study implied that lower dilutions activated mechanisms of immune response in a higher level in a shorter period while higher dilutions activated only the first defense mechanisms in A. ventricosus, which may be dependent on the exposure time between the organisms and the HDIC. These results showed the complex dynamics between the non-self-recognition mechanisms of A. ventricosus juveniles and HDIC in a long time exposure. Additionally, mitochondria, besides its role in energy production, may play an important role in defense response activation and signalling process when organisms are exposed to HDIC. Furthermore, this is the first assembled transcriptome of the scallop A. ventricosus juveniles using RNA-seq technology that allows us to generate biological markers for future investigation of this important resource, whose populations are declining in Baja California Sur, México.