Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

Parallel analysis of miRNAs and mRNAs suggests distinct regulatory networks in Crassostrea gigas infected by Ostreid herpesvirus 1

BMC Genomics volume 21, Article number: 620 (2020) Cite this article

Abstract

Background

Since 2008, the aquaculture production of Crassostrea gigas was heavily affected by mass mortalities associated to Ostreid herpesvirus 1 (OsHV-1) microvariants worldwide. Transcriptomic studies revealed the major antiviral pathways of the oyster immune response while other findings suggested that also small non-coding RNAs (sncRNA) such as microRNAs might act as key regulators of the oyster response against OsHV-1. To explore the explicit connection between small non-coding and protein-coding transcripts, we performed paired whole transcriptome analysis of sncRNA and messenger RNA (mRNA) in six oysters selected for different intensities of OsHV-1 infection.

Results

The mRNA profiles of the naturally infected oysters were mostly governed by the transcriptional activity of OsHV-1, with several differentially expressed genes mapping to the interferon, toll, apoptosis, and pro-PO pathways. In contrast, miRNA profiles suggested more complex regulatory mechanisms, with 15 differentially expressed miRNAs (DE-miRNA) pointing to a possible modulation of the host response during OsHV-1 infection. We predicted 68 interactions between DE-miRNAs and oyster 3′-UTRs, but only few of them involved antiviral genes. The sncRNA reads assigned to OsHV-1 rather resembled mRNA degradation products, suggesting the absence of genuine viral miRNAs.

Conclusions

We provided data describing the miRNAome during OsHV-1 infection in C. gigas. This information can be used to understand the role of miRNAs in healthy and diseased oysters, to identify new targets for functional studies and, eventually to disentangle cause and effect relationships during viral infections in marine mollusks.

Background

Mollusk aquaculture is regarded as the fastest growing food production sector and plays a vital role in solving the problem of feeding future human generations [1, 2]. However, the recurrence of infectious disease as well as intensive mono-specific farming of bivalve species threaten this production sector [3, 4]. Global oyster production has particularly suffered from Pacific oyster mortality syndrome (POMS), a deadly disease mainly caused by Ostreid herpesvirus 1 (OsHV-1) infection [5]. OsHV-1 is one of the two dsDNA viruses of the Malacoherpesviridae family [6], that can infect oysters and other bivalve species [7,8,9,10], whereas the other, Haliotid herpesvirus 1 (HaHV-1), causes a disease called abalone viral ganglioneuritis in gastropod species [11, 12]. Different genomic variants of these two viruses have been identified in several disease outbreaks and the virus presence in healthy individuals or non-susceptible species, as reservoirs for viral particles, makes virus eradication almost impossible. Hence, the understanding of virus life strategies and host-pathogen interactions in different host species is essential for the prevention and management of mass mortalities [13]. At present, the lack of mollusk cell lines [13] hampers the in vitro propagation of malacoherpesviruses, even if cultured hemocytes have been recently proposed as a tool to propagate the virus and to perform small-scale infection trials [14, 15]. The available knowledge on Malacoherpesviridae biology refers to in vivo infections, either based on laboratory trials or on the rare cases of naturally infected bivalves and, mostly, by high-throughput (HT) sequencing [16]. Notably, recent studies suggest that the weakening of immune defenses consequent to OsHV-1 infection is likely to induce critical changes in the oyster-associated microbiota and, hence, fatal secondary infections [17, 18]. The transcriptional host response has been characterized in Crassostrea gigas and Scapharca broughtonii infected with OsHV-1 [19,20,21] and, more recently, for Haliotis spp. infected with HaHV-1 [22, 23], revealing distinct antiviral responses in the different virus-host combinations. Infected C. gigas mainly activated the interferon pathway to produce antiviral molecules known as interferon stimulated genes (e.g., viperin) [19, 21, 24]. Differently, S. broughtonii infected with the same virus mostly modulated apoptosis-related genes similar to H. diversicolor infected with HaHV-1 [20, 22]. These mechanistic differences show that we still lack a comprehensive understanding of the transcriptional response to viral infections in mollusks.

Besides the expression of coding RNAs, non-coding RNAs (ncRNAs) such as small and long non-coding RNAs (sncRNAs, lncRNAs) can play an important role during viral infections [25]. ncRNAs are known as post-transcriptional regulators [26], and they can participate in virus-host interactions through a variety of mechanisms [27, 28]. Among the ncRNAs, miRNAs are the most frequently studied due to their involvement in various diseases and their potential as diagnostic and prognostic biomarkers [29]. Mature miRNAs are short sequences (20–23 nt), which originate from stem-loop precursor RNAs (pre-miRNAs) and act as post-transcriptional repressors by targeting untranslated regions (UTRs) of messenger RNAs (mRNAs) in biological processes such as immunity, development, cell behavior and host-microorganism interactions [30, 31]. After the first discovery of lin-14 in Caenorhabditis elegans [32], miRNAs have been identified in many animals, plants, and viruses, suggesting independent miRNA origins throughout the major evolutionary lineages [26]. Virus-encoded miRNAs have been reported in vertebrate herpesviruses, for instance in Epstein-Barr virus [33] where they perform immuno-modulatory functions [34]. Regarding invertebrate viruses, miRNAs have only been reported in the White Spot Syndrome Virus (WSSV) [35]. WSSV miRNAs target the JAK-STAT signaling pathway to interfere with the immune system of the shrimp [36], and their sequences are actively edited by host enzymes, like the adenosine deaminase acting on dsRNAs (ADAR-1) [37].

In bivalves, only a few studies investigated the presence and role of miRNAs in relation to biomineralization [38,39,40] and neuro-immunity [41, 42], whereas a single study investigated the roles of miRNAs in scallops (Chlamys farreri) infected by an OsHV-1 variant [43] (the related data are not public). Therefore, we have examined the relationship between coding and non-coding RNAs by parallel HT-sequencing of mRNA and sncRNA in oysters naturally infected by OsHV-1. Based on the coupled transcriptomic landscapes of mRNAs and sncRNAs we add an additional facet to the characterization of the molecular actions and counteractions in OsHV-1 infections of oysters.

Results

OsHV-1 infection in field-exposed oysters

We detected variable amounts of OsHV-1 DNA in 15 oysters collected from Goro lagoon on May 17, 2016 (Fig. 1a). In 8 of the 15 samples, the expression levels of OsHV-1 ORF104, a proxy for total viral transcription, were higher than 1% of the expression level of the host housekeeping gene Elongation factor-1-alpha (Fig. 1a). We further selected 6 oyster gill samples, representative of a range of OsHV-1 DNA and RNA levels, to perform parallel high-throughput sequencing of sncRNAs and (poly(A)-tail selected) mRNAs. In detail, we sequenced two samples (S2 and S5) with high viral RNA: DNA ratios (δ) and low absolute amounts of OsHV-1 DNA, two samples (S1 and S4), with an intermediate δ value and intermediate to high levels of OsHV-1 DNA and two samples (S3 and S6) with the lowest δ values and intermediate to high levels of OsHV1 DNA (Fig. 1b, Table 1). Taking into consideration the presence of OsHV-1 DNA, high and intermediate δ values suggest active OsHV-1 in the early phases of the infection, whereas low δ values suggest a limited transcription of abundantly present OsHV-1. This latter situation could be interpreted as a late infection stage. Overall, the high-throughput sequencing of the six samples yielded 94.2 M sncRNA reads (size range: 18–40 nt) and 380 M mRNA reads (Table 1).

Fig. 1
figure 1

OsHV-1 RNA and DNA levels in the gill tissue of 15 OsHV-1-positive oysters (C. gigas, S1-S15) sampled in the Goro lagoon, Italy. a. Quantification of OsHV-1 transcription by RT-qPCR (bars, left axis) and OsHV-1 DNA loads (black diamonds, right axis). The OsHV-1 ORF104 transcript levels normalized to the expression of the C. gigas housekeeping gene elongation factor 1-apha, was considered as a proxy of the viral transcriptional activity. OsHV-1 DNA content was measured as DNA copy number per μl. Grey bars indicate the six samples selected for sncRNA and mRNA HT-sequencing. In these six samples, red dots represent the number of RNA reads mapping to the OsHV-1 genome. The samples denoted by grey bars were selected for RT-qPCR analysis. b. Subdivision of samples S1-S6 based on the ratio of OsHV-1 RNA over DNA (δ value) grouping pairs of samples into low, mid, or high δ samples

Full size image
Table 1 HT-sequencing results showing the amounts of mRNA and sncRNA reads in the oyster samples. Sample classification was based on δ values. Amount of OsHV-1 DNA (copies per μl), number of clean reads and number of oyster and OsHV-1 reads are reported for each library. For the sample S6, details of the ribo-depleted library are also reported
Full size table

miRNA expression during OsHV-1 infection

The size distributions of the sncRNA reads in the six samples showed a clear 20–22 nt peak, typical of the presence of miRNAs, with a secondary peak around 30 nt, possibly related to the presence of piwi-interacting RNAs (Fig. 2a). The percentage of the sncRNA reads mapping to the oyster genome ranged from 90.7 to 97.1%, whereas only a small fraction of the reads mapped to the OsHV-1 genome (0.03–0.07%, Table 1). Overall, 46% of the sncRNA reads mapped to the 151 C. gigas miRNA precursors retrieved from MirGeneDB v.2.0 [44], showing a clear 22-nt peak (Fig. 2c). The reads not mapping to the 151 oyster miRNA precursors mostly found a match in the C. gigas genome, showing a clear 29-nt peak (Fig. 2d), while the reads matching to the OsHV-1 genome showed shorter size (17–19 nt, Fig. 2e). We verified the presence of the minimal miRNA annotation criteria for the 151 oyster miRNA precursors, including read coverage on both miRNA arms, 5′ read homogeneity and the absence of read mapping in the surroundings of the miRNA arms [45]. Accordingly, we confirmed most of these oyster miRNA predictions with the following exceptions: i) we could not find reads mapping to the star arm for some miRNAs such as Cgi-mir-96-P3, Cgi-mir-87, Cgi-novel-4, Cgi-novel-18 and several Cgi-mir-184 isoforms, ii) we found equal coverage of both mature and star arms in the case of 4 miRNA precursors, and iii) we reverted the mature and star predictions because of differential coverages for other 7 miRNA precursors (Cgi-mir-36, Cgi-mir-1992, Cgi-novel-1, 2, 8, 13 and Cgi-novel-15).

Fig. 2
figure 2

sncRNA analysis. a. Size distribution of the sncRNA reads in the six libraries (size range: 18–40 nucleotides). b. Cumulative (S1-S6) expression of C. gigas miRNAs classified according to MirGeneDB v.2.0 (57 miRNAs contributing to less than 0.1% to the global expression were removed; data are available in Additional file 1). c. Cumulative size distribution of the sncRNA reads mapping to the oyster miRNA precursors. D, E. The sncRNA reads not matching to oyster miRNAs were mapped to the oyster (D) or to the OSHV-1 genome (E)

Full size image

We could compute expression values for 132 out of 151 miRNAs (Additional file 1), although 57 of these 132 miRNAs cumulatively contributed less than 0.1% to the global expression and included most of the Cgi-mir-184 isoforms as well as 6 out of 7 novel oyster miRNAs (Cgi-novel-1, 2, 5, 10, 21, 22, Additional file 1). Cgi-mir-10-P2, Cgi-mir-1, Cgi-mir-279 and Cgi-mir-184-P7 resulted to be the most expressed miRNAs, accounting for 18.9, 11.6, 11.2 and 10.9% of the global miRNA expression, respectively (Fig. 2b and Additional file 1). Among them, Cgi-mir-184-P7 and Cgi-mir-10-P2 displayed the most stable expression levels over the six tested samples, with a coefficient of variation of 3 and 9%, respectively. These miRNAs were tested as references for RT-qPCR analysis. Two Cgi-mir-375 isoforms, Cgi-mir-750, Cgi-mir-1175 and Cgi-novel-19 were the most variable miRNAs among the 6 samples, with a CV > 100%. Four out of these (two Cgi-mir-375 isoforms, Cgi-mir-750 and Cgi-mir-1175) mapped in near vicinity to each other on the oyster genome. Next to this cluster, another set of neighboring miRNAs (Cgi-mir-12 and two Cgi-mir-216 isoforms) also showed considerable variation in their expression levels (Additional file 1).

Overall, the miRNA expression changes seemed to be influenced by OsHV-1 activity. Considering the expression of all miRNAs together in a Principal Component Analysis (PCA), we could confirm the three δ groups, with the ‘high’ group separated by axis one, and the ‘low’ group separated by axis 2 (Fig. 3a). Accordingly, we searched for differentially expressed miRNAs (DE-miRNAs) in pairwise comparisons between the groups. The strongest response could be observed in the ‘low’ group, which had 11 DE-miRNAs compared to the ‘mid’ group and 10 DE-miRNAs compared to the ‘high’ group (with 8 DE-miRNAs shared between both contrasts, Fig. 3 and Table 2). The list of the 15 DE-miRNAs included only one novel miRNA (Cgi-novel-19). Cgi-novel-10 was expressed almost exclusively in the “low δ “samples S3 and S6, but high variation in expression levels yielded no statistically significance when applying the same criteria as for the other DE-miRNAs.

Fig. 3
figure 3

a. Principal Component Analysis of the miRNA expression in samples S1-S6. b. Venn diagram depicting the differentially expressed miRNAs in the pairwise comparisons of the samples grouped by δ values

Full size image
Table 2 Differentially expressed miRNAs in pairwise comparisons among samples S1-S6. Table heading: miRNA ID, the comparison in which a specific miRNA resulted to be a differentially expressed and miRNA expression values per samples are reported for the 15 DE-miRNAs
Full size table

Validation of miRNA expression by RT-qPCR

According to overall expression patterns, we selected 8 miRNAs for expression validation by RT-qPCR using the six RNA samples used for sncRNA sequencing plus two samples selected among the 15 and representing a case of low (sample S7) and a high infection (sample S8). We selected Cgi-mir-184-P7 as housekeeping miRNAs because of its low variation of mapped reads among samples and stability compared to the spiked RNA Sp6. Cgi-mir-133, Cgi-mir-315, Cgi-mir-1985 and Cgi-Novel-19 were chosen within the DE-miRNAs to cover different expression ranges. Additionally, we included two miRNAs (Cgi-mir-750 and Cgi-novel-10) because of their contrasted expression patterns in sncRNA-seq data. After data normalization, the correlation between sncRNA-seq and RT-qPCR expression levels for the six samples (S1-S6) ranged from a r2 of 0.85 to 0.99 (Additional file 2). Following RT-qPCR, also the two additional samples showed miRNA expression values that matched the trends obtained by HTS, further supporting the link between expression of the majority of miRNAs consistently depends on the OsHV-1 infection intensity (Additional file 2). When normalizing miRNA expression values to sample S5 (one of the two samples denoted by a high δ value and characterized by the low OsHV-1 DNA load and transcription), Cgi-mir-750 was one of the most expressed and induced miRNAs (Additional file 2 and Fig. 4). The oyster specific miRNA Cgi-novel-10 was also highly induced in samples with comparatively higher viral activity, fitting well to the sample grouping based on δ values. In contrast, the other oyster specific miRNA Cgi-novel-19 did not fit according to such sample grouping. Cgi-mir-133 and Cgi-mir-1985 and Cgi-mir-315, which showed intermediate expression levels (Additional file 2), were mildly induced or downregulated compared to sample S5 (Fig. 4).

Fig. 4
figure 4

RT-qPCR analysis. The expression values of six selected oyster miRNAs over 8 samples are expressed as fold change versus the values measured in the sample S5. Samples are colored following the delta values, green for high value (S5), bluish for intermediate values (S1 and S4), reddish for low values (S3 and S6). The two additional samples (S7 and S8) for which we cannot compute a delta value are reported in white (low infection) or black (high infection)

Full size image

Expression of C. gigas coding genes during OsHV-1 infection

Overall, 87–92.4% of the mRNA reads mapped to the oyster genome under the applied parameters (Table 1). The PCA analysis based on the whole expression profiles outlined two groups along axis 2 (Fig. 5). These groups did not completely match the grouping of relative OsHV-1 expression (i.e. δ value), but rather corresponded to the overall OsHV-1 transcriptional activity in the samples (Fig. 1). Within the samples showing high viral transcription (S1, S4 and S6), S6 was separated along axis 1 from samples S1 and S4, while samples with low OsHV-1 transcription (S2, S3 and S5) clustered together on both axes (Fig. 5). The comparison between these two groups identified 403 DEGs, 224 up-regulated and 179 down-regulated (Additional file 3). Noticeably, a considerable proportion of DEGs corresponded to proteins with unknown function annotated as “hypothetical protein”. Among the upregulated DEGs, we found components of the interferon pathway (IRF8, 6.8x; IRF2, 3.0x and IFRD1, 2.8x), components of the Toll pathway (MyD88, 28.6x; IRAK4, 4.2x; cact, 3.4x; IKBKE, 2.3x) caspases (Casp-8, 2.8x and Casp-7, 2.0x), as well as other components known to be involved in oyster antiviral pathways (viperin, metalloproteinase inhibitors, baculoviral IAP repeat-containing proteins, dual specificity protein phosphatase 3, Bcl-2-like protein 1, SOCS2, Dual specificity mitogen-activated protein kinase kinase 3, Additional file 3). One gene, encoding for the interferon-stimulated enzyme ADAR-1, was strongly upregulated (6.4x) and its expression correlated with the OsHV-1 RNA levels in the six samples (Fig. 5b). Also, upregulated DEGs included several receptors possibly linked to the neuroendocrine system, like FMRF-amide, prostaglandin, dopamine, melatonin and IL-17 receptors as well as genes involved in the Pro-PO system (tyr-3 and Laccase-2, Additional file 3).

Fig. 5
figure 5

a. PCA analysis of samples S1-S6 based on the C. gigas mRNA expression profiles. The samples divided according to the presence of low or high OsHV-1 transcription levels. b. Correlation between the expression values of oyster ADAR-1 and the number of OsHV-1 RNA reads in the six samples

Full size image

Tracing the OsHV-1 transcription during infection

We analyzed the OsHV-1 expression profiles in the 3 samples showing a sufficient number of viral reads (samples S1, S4 and S6, Table 1). We took advantage of the strand-specificity of our RNA-seq libraries to map the reads with directional constraints on the OsHV-1 genome, either in the expected coding-sense orientation, or as antisense mapping. While the “sense” mapping enabled us to get a more precise quantification of gene expression by counting the reads belonging to a given mRNA, the “antisense” mapping is a measure of antisense transcription levels, possibly produced by unknown viral transcripts located on opposite strands. The OsHV-1 “sense” profiles showed expression of most of the predicted genes, as previously reported for this viral family [18, 19, 22]. The most expressed OsHV-1 genes in our dataset were: ORF107, a protease called assemblin; ORF104 and ORF82, putative capsid proteins; ORF88, a possible membrane protein; ORF99, an apoptosis inhibitor. Other highly expressed viral genes denoted proteins with unknown functions, like ORF33, ORF119, ORF113, ORF41 and ORF127 (Additional file 4). Our samples showed lower averaged expression levels, but averaged expression peaks 2.4 times higher, in comparison to other seven RNA-seq samples resulting from experimental OsHV-1 infection (Fig. 6a). As a result of the heterogenous expression, the 5 most expressed OsHV-1 genes contributed 40–50% of the total expression in samples S1, S4 and S6, while they only contributed 21–26% of the total expression in previous studies [18, 19] (Fig. 6b). We also found that a small proportion of the reads from each sample (6.4–10.6%) mapped opposite to the expected coding directions. Most of these reads were however located at the ORF boundaries, suggesting that this phenomenon is caused by overlapping UTRs. This might reflect the dense distribution of OsHV-1 genes in combination with the lack of knowledge about the extension of OsHV-1 UTRs. In contrast, we also observed longer antisense signals along ORF100 (DNA polymerase) and ORF22, possibly indicating the presence of antisense genes that were not yet annotated (Additional file 5). Taking advantage of the abundant sncRNA reads, we aimed to investigate if OsHV-1 encoded genuine miRNAs. When grouping the putative OsHV-1 sncRNA reads into 19,997 clusters of identical reads, most clusters (97.3%) were represented by less than 10 reads. The distribution of OsHV-1 sncRNA read lengths was skewed towards shorter reads and did not show the distinctive peaks typical of sncRNAs (Fig. 2c-e). This distribution alone does not lend a lot of support for the presence of genuine miRNAs along the OsHV-1 genome. We nevertheless investigated the presence of possible structured RNAs further by using the VMir tool [46]. This resulted in the identification of 236 hairpins, covered by a total of 1456 sncRNA reads. However, the coverage graphs of these hairpins did not fulfill the minimal criteria for the identification of bona fide miRNAs. Therefore, based on our data, we suggested that OsHV-1 does not encode genuine miRNAs and that OsHV-1 sncRNA reads rather originated from mRNA degradation. Further evidence for a mRNA origin of the viral sncRNA reads comes from the matching SNP patterns in sncRNA and mRNA reads at positions edited by ADAR-1. The expression of oyster ADAR-1 correlated with the quantity of OsHV-1 RNA (Fig. 5b) and ADAR-1 exerts its enzymatic activity by editing dsRNA with a mechanism known as A-to-I editing, thus generating an identifiable footprint of ‘G’ mismatches, as we previously demonstrated [47]. In samples S1 and S6 we could identify 79 and 110 ADAR-mediated SNPs occurring with low frequencies (mean of 5.8 and 4%, respectively). Since almost all these edited positions could be traced also in the sncRNA reads, a degradation of full-length mRNA seems more likely than the expression of genuine viral-encoded miRNAs.

Fig. 6
figure 6

OsHV-1 expression analysis. a. The expression values of samples S1, S4 and S6 were compared with the ones obtained from 7 samples obtained from an experimental infection with OsHV-1 [18]. b. Expression profiles of selected OsHV-1 genes in all 10 samples (the samples of the experimental infection are named as time point (hpi, hours post infection)

Full size image

miRNA/mRNA interactions during OsHV-1 infection

Since miRNAs are expected to regulate the expression of coding genes by interacting with 3′-UTRs, we investigated the presence of possible miRNA-mRNA expression correlations among the six samples (Fig. 7). When comparing the correlations of predicted miRNA-mRNA interactions and especially the predicted interactions between differentially expressed DE-miRNAs and mRNAs with the Null distribution of all possible correlations we observed a lack of positive correlations. As expected, we found a comparatively strong increase of negative correlations, suggesting that miRNAs mostly repress gene expression [48]. In detail, the miRNAs with the highest proportions of strong correlations (i.e. the top and bottom 2.5% of all possible correlations) were Cgi-mir-277, Cgi-mir-9, Cgi-mir-315, Cgi-mir-1 and Cgi-novel-10 (Fig. 6b). While Cgi-mir-9 and Cgi-mir-1 are highly expressed miRNAs with little variation among samples, Cgi-mir-315 is included in the list of DE-miRNAs. Since UTRs are not annotated in the available oyster genome [49], we identified the 3′-UTRs by mapping our RNA-seq reads on the oyster genome with a mapper allowing the presence of large gaps due to introns. According to the available gene annotations, we could predict 2074 3′-UTRs longer than 30 nt. These sequences showed an average length of 441 nt. Using this dataset, miranda [50] predicted 1425 possible matches targeting 358 oyster genes. A total of 68 interactions, targeting 50 genes were assigned to the 15 DE-miRNAs (Table 3). However, only a few of these interactions affected differentially expressed oyster genes. The latter included Chloride channel protein 7, Laminin subunit beta-2, Collagen alpha-1(XIV) chain, Endo-1,6-beta-D-glucanase BGN16.3, Kelch-like protein 20, Achaete-scute-like protein 1, Tribbles-like protein 2 and two unknown proteins (Table 3). Next to the DE-miRNAs we also focused on miRNA-mRNA interactions of miRNAs that showed the highest proportions of strong miRNA-mRNA correlations (e.g. Cgi-mir-277 and Cgi-novel-10, Fig. 6b), or showing the highest coefficient of variation among the six samples (Cgi-mir-750 and Cgi-mir-1175, see Additional file 1). These two miRNA groups showed a similar number of matches, with few genes included in the DEG list. Among these, the match between Cgi-mir-277 and a serine/threonine protein kinase stood out, that we found to be preferentially expressed in samples of OsHV-1-infected oysters.

Fig. 7
figure 7

miRNA-mRNA expression correlation. a. Density plot of Pearson correlations for the combination of all oyster genes vs. all miRNA (black line), for all the miRNA-mRNA interactions predicted by miranda (red dotted line) and for the predicted interactions involving DE-miRNAs (blue line). b. Proportion of strong correlations for predicted miRNA-mRNA interactions (green bars correspond to the five miRNAs with the highest number of interactions)

Full size image
Table 3 Putative miRNA-mRNA interactions. The possible miRNA-mRNA interactions are listed for the DE-miRNAs and for other miRNAs of interests. Bolded hits represent interactions involving DEGs. The gene IDs and putative annotations are reported according to the oyster genome annotation (v.9) and the IDs without annotation refer to proteins with unknown function (annotated as “hypothetical protein”). The underlined match to Histone deacetylase 6) refers to a miRNA-mRNA interaction described also in humans
Full size table

Using miranda we also identified 307 matches between oyster miRNAs and the OsHV-1 genome, and 24 of these involved DE-miRNAs. However, the absence of information regarding the extent of viral UTRs as well as the low number of samples with enough mRNA data (n = 3) prevented the possibility to link miRNA matches to a given viral gene.

Discussion

To increase our understanding of the mechanistic role of sncRNAs as gene expression regulators during OsHV-1 infection in oysters, we simultaneously sequenced the sncRNAs and mRNAs from six oysters naturally infected with OsHV-1 in the Goro lagoon (Italy). While we could clearly see the onset of antiviral oyster immunity in response to OsHV-1 infection, expression patterns of oyster miRNAs differed from mRNA transcription in several aspects. Most obvious was the clustering of miRNAs according to virus activity, which was set off against a grouping by viral transcript abundance. Furthermore, only few of the miRNA-mRNA correlations affected antiviral genes responding to OsHV-1 infection. This lack of consistency between miRNA and mRNA profiles related to the antiviral response might be a sign of the inherent variability occurring in natural, uncontrolled, infections. Alternatively, the weak correlation between miRNA and mRNA transcription profiles may indicate a limited regulatory role of miRNAs in oyster antiviral processes. Either way, our exploration of the miRNAome landscape in response to OsHV-1 infection indicated sophisticated miRNA regulatory networks with only loose connections to the oysters’ antiviral immune response.

Expression of miRNA diversity in oysters

By using only the 151 miRNA predictions available for oyster in the MirGeneDB [44], we adopted a conservative approach that can prevent the inclusion of false positive miRNAs, typically found in sncRNA sequencing studies of model organisms [45, 48, 51]. Expression analysis revealed that Cgi-mir-10-P2 (previously reported as mir-100), Cgi-mir-279, Cgi-mir-184-P7 and Cgi-mir-1 accounted for most of the sncRNA reads, making Cgi-mir-184-P7 the most suitable housekeeping miRNA for data normalization. The high expression levels of these miRNAs is in agreement with previous studies in oyster [38, 52], Chlamys farreri [43] and Tegillarca granosa [53]. For other miRNAs (like Cgi-mir-67 or Cgi-let-7), we measured lower expression levels than those reported before [38, 52]. We also observed reversed amounts of mapped reads between the mature and star arms for some miRNAs (Cgi-mir-36, Cgi-mir-1992, Cgi-novel-1, 2, 8, 13 and Cgi-novel-15). Although both arms can be selected to produce mature and functional miRNAs, a preference towards one of the two arms is the rule [26], resulting in a bias of the number of mapped reads among arms. The mechanism allowing such a selection is still unknown and arm-imbalance was recently reported as a way to modify miRNA targets during cancer [54]. However, most of the observed “reversions” came from novel oyster miRNAs, possibly because limited knowledge is available on these new oyster-specific miRNA families.

Overall, we identified 15 oyster miRNAs that were differentially expressed among samples. The resulting miRNA expression profiles suggested a sample grouping according to the ratio of OsHV-1 RNA over DNA (i.e. δ values, the relative transcriptional activity), whereas expression of C. gigas coding genes rather grouped according to absolute viral transcription levels. This would indicate the existence of more subtle regulative mechanisms depending upon OsHV-1 stage that control miRNA expression, compared to the regulation of coding genes. Notably, some of the DE-miRNAs were located in clusters in the oyster genome (Cgi-mir-750 with Cgi-mir-1175 and Cgi-mir-12 with Cgi-mir-216) and showed similar expression values between δ groups, suggesting co-expression and potential synergistic functions, previously shown for humans [55]. Several of the DE-miRNAs belong to miRNA families that were previously described to modulate immunity during viral infections in other organisms. This included mir-12 that was found upregulated in WSSV-infected shrimps, where it modulates phagocytosis, apoptosis and antiviral immunity [56]. Two other miRNAs, mir-375 and mir-750, were also highly responsive to WSSV infection in Panaeus monodon [57]. Additionally, mir-315 regulated the pro-PO system during WSSV infection, thereby inhibiting the spread of the virus [58]. Both in sncRNA-seq and RT-qPCR data, Cgi-mir-750 was one of the most induced miRNAs in samples with high OsHV-1 infection and it was the miRNA with the highest number of matched among DEGs.

As previously reported by using SSH libraries [59], an upregulation of pro-PO related genes could be detected during OsHV-1 infection and an increase of the PO activity was also evident in C. farreri infected with a OsHV-1 congener [60]. Accordingly, we found a few pro-PO related DEGs upregulated in samples showing OsHV-1 activity, including a laccase and a tyrosinase. A similar laccase (lac-2) was recently reported as strongly upregulated during WSSV infection in Litopenaeus vannamei [61]. Among the oyster’s countermeasures against the OsHV-1 infection we can include the interferon-like pathway, the toll-pathway as well as apoptosis and pro-PO activity in a conceptual model of the oyster response to OsHV-1 infection (Fig. 8). Moreover, the trace of ADAR-1-mediated editing of viral dsRNAs in samples S1 and S6, together with a good correlation between ADAR-1 expression and OsHV-1 transcription confirm our previous findings suggesting the main role of this enzyme in editing exogenous dsRNAs [47], although the biological meaning of this editing during oyster-OsHV-1 interaction is still unknown.

Fig. 8
figure 8

Conceptual model of possible molecular events describing OsHV-1 and Crassostrea gigas interactions in the infected oyster cells. Replicating OsHV-1 DNA amplifies the production of viral RNAs, dsRNAs and proteins necessary for virion assembly and responsible of some host-pathogen interactions (e.g. OsHV-1 IAPs). The binding of viral DNAs and dsRNAs to specific oyster receptors, namely endogenous TLRs (eTLRs) and RIG-I/MDA5 proteins, activates the Toll and Interferon pathway, respectively (red boxes) and leads to the transcription of antiviral effectors (green boxes). Pro-apoptotic genes, like caspases, Pro-PO elements, like tyrosinases and laccases, and interferon stimulated genes, like viperin and ADAR-1 are upregulated during OsHV-1 infection. These antiviral effectors control the virus, which counteracts by expressing anti-apoptotic viral genes (IAPs, like ORF99). Oyster ADAR-1 edits dsRNAs with a mechanism known as A-to-I editing, producing G mismatches that impair dsRNAs, and possibly making the edited dsRNAs less effective in activating dsRNA receptors, while the impact on OsHV-1 replication is unknown. A few oyster miRNAs are regulated during OsHV-1 infection, but their function in controlling host and viral genes remain unclear

Full size image

Conceivably, all these well-known antiviral pathways are poorly interconnected with the predicted miRNA matches. In fact, for some predicted miRNA-mRNA interactions we found strong correlations in expression values indicative of miRNA-mediated co-regulation processes. However, we could only identify a limited number of possible interactions between miRNAs and 3′-UTRs of mRNAs. As regards the matches of DE-miRNAs reacting to OsHV-1 infection, we found that these matches mostly involved genes that were neither directly related to antiviral responses, nor included in the list of DEGs. The possible interaction between Cgi-mir-277 and a serine/threonine protein kinase likely represents an exception here, since this transcript was upregulated in our and other OsHV-1-infected samples, and serine/threonine protein kinases are known to be involved in differential resistance to viral infection [62]. Therefore, cgi-mir-277 could target a functionally conserved mechanism in viral immune responses. The other DE-miRNAs may be involved in indirect regulative mechanisms or target gene transcripts not yet known to be involved in antiviral oyster processes, but overall, the regulatory role of miRNAs in the antiviral response of the oyster seems to be limited.

Interaction of miRNAs with viral genes

Oyster miRNAs could also target viral infections by binding to viral genes and thereby limiting or promoting viral gene expression. Out of the prediction of 307 positions the OsHV-1 genome, 40 matches were exclusive to oyster miRNA families, which might be a result of co-evolution between C. gigas and OsHV-1. And while we found several possible miRNA-OsHV-1 matches including highly expressed viral genes (e.g. ORF88) or genes in the vicinity of highly expressed genes (e.g. ORF108 near ORF107, i.e. the most expressed OsHV-1 gene), we could not find compelling evidence for a functional relationship between oyster miRNAs and viral genes due to the low number of samples with sufficient viral reads. Furthermore, the extent of OsHV-1 UTRs is presently unknown and the high rate of antisense transcription that we found is a clear indicator for the poorly understood complexity of the OsHV-1 transcriptome. Based on this and on the fact that most of the OsHV-1 genes have unknown function, any assumption of a functional meaning of these matches would be very speculative. Yet, the enrichment of negative correlations in predicted miRNA-mRNA interactions along with the previously reported matches from other organisms still suggests, that the miRNA machinery can play an important role in regulating transcription in oysters, but this role might be limited for the interaction with OsHV-1.

For the virus itself, it was not yet clear whether it utilizes miRNA to regulate its own transcription or to modulate host immunity. Our search for OsHV-1-encoded miRNAs suggested that OsHV-1 does not utilize miRNAs. This was not self-evident since other viruses related to OsHV-1 [6], as well as WSSV [16, 63] can encode miRNAs. However, the OsHV-1 sncRNA reads showed a size distribution lacking abundant size classes typical of miRNAs. While genome-wide prediction tools suggested the presence of RNA hairpins along the OsHV-1 genome [16], our data did not support evidence for increased expression of these loci, nor coverage profiles consistent with the presence of genuine miRNAs. Therefore, we rather suggest that OsHV-1 sncRNA reads originated from degraded mRNAs. Also, the size profile of OsHV-1 sncRNA reads suggested that RNAi is not active against OsHV-1 in oyster, since a size-bias would be expected in case of RNAi-mediated RNA degradation [64]. RNAi was reported during WSSV infection in shrimp [35] as well as in other arthropod species [65], but this study rather supports previous findings indicative of a marginal role for RNA interference in the antiviral response of bivalve species [66, 67].

Conclusions

The parallel characterization of sncRNA and mRNA data from six oysters naturally infected with OsHV-1 revealed the revealed the regulatory network of the oyster miRNActome. Generally, miRNAs could serve molecular targets for future studies trying to control disease outbreaks. While the overall miRNA expression profiles and some specific miRNAs (for example, Cgi-mir-750, Cgi-mir-277, Cgi-mir-375) correlated with the transcriptional activity of the virus, we found only few predicted interactions with genes known to be involved in viral immunity. In contrast, we observed large differences in the mRNA transcriptional profiles, characterized by the upregulation of several antiviral and immune pathways. The miRNAs transcription profiles, however, showed differential expression of few miRNAs possibly linked to their regulative roles in the onset or in the control of the OsHV-1 infection. The coupling of mRNA and sncRNA sequencing showed that while several major immune pathways are activated during an infection, with miRNAs only playing a limited direct role in the viral infection of oysters (Fig. 8). Nevertheless, the altered expression profiles of miRNAs in response to viral transcriptional activity could indicate that more indirect roles exist for some of these miRNAs possibly leading to functional consequences in terms of host-pathogen interactions (e.g. Cgi-novel-10, Fig. 8).

Methods

C. gigas deployment in open field and sampling

Triploid Crassostrea gigas seed of French origin (T6 size, 0.15 g on average) was placed in lantern-like baskets at 0.5–1 m depth in the Goro lagoon (44°48.728′N 12°17.905′E, North Adriatic Sea, Italy) on March 19th, 2016. The growing oysters were sampled from Apr 26th to Jun 21st (26/4, 03–10–17-24-31/05 and 07–14-21/06), a time span which previously included OsHV-1 infections and mortality [68]. Within a month the oyster spat reached a shell length of 2.46 ± 0.12 cm (N = 30, 17th May) and continued to grow during the whole monitoring period (Additional file 6). From the oysters collected on May 17th (2.3–2.7 cm shell length, N = 45), a fragment of mantle and the whole gills were individually dissected on ice, immerged in RNA Stabilization Reagent (RNAlater, Qiagen, Milano, Italy) and stored at − 80 °C. Fragments of mantle and gill tissues were used for OsHV-1 DNA quantification and for individual RNA purification.

Analysis of OsHV-1 DNA

The presence of OsHV-1 DNA was measured by quantitative PCR. Total DNA was extracted from ~ 25 mg w.w. of oyster gill and mantle using the QIAamp® DNA Mini Kit following the manufacturer’s instructions (Qiagen). The purified DNA samples were quantified with a NanoDrop spectrophotometer (ThermoFisher Scientific, Waltham, MA, USA) and diluted to 5 ng/μl for the amplification reaction targeting the catalytic subunit of the viral DNA polymerase (AY509253 ORF100) [69]. Five μl of DNA were added to a qRT-PCR reaction mix composed of 12.5 μl SsoFast™ EvaGreenR Supermix (Bio-Rad Laboratories, Segrate, Milano), 1.25 μl of each primer (HVDP-F: 5′ ATTGATGATGTGGATAATCTGTG 3′; HVDP-R: 5′ GGTAAATACCATTGGTCTTGTTCC 3′) diluted at the concentration of 0.5 μM, and 5 μl of water. The amplification reactions were performed in a Rotor-Gene Q thermocycler (Qiagen) as follows: 1 cycle of polymerase activation at 95 °C for 5 min; 40 cycles of amplification at 95 °C for 30 s, 60 °C for 1 min, and 72 °C for 45 s and a final step for melting temperature curve analysis from 65 to 95 °C (10 s/step, ramp rate 0.5 °C/sec) [70]. The absolute number of OsHV-1 DNA copies/μl was determined by comparing the Ct values resulting from a standard curve. Plasmidic DNA including the OsHV-1 target region was serially diluted 1:10 in the range 10–106 DNA copies/μl and used to compose the standard curve.

RNA extraction, library preparation and sequencing

Total RNA was extracted from individual oyster gills according to the Trizol manufacturer’s instructions (Thermofisher Scientific), quantified with a Qubit 2.0 fluorometer (Thermofisher Scientific), and the RNA quality was tested with an Agilent bioanalyzer using the RNA6000 pico kit by automated electrophoresis and fluorescence signal detection (Thermofisher Scientific). A Real Time PCR (qRT-PCR) approach was used to estimate the OsHV-1 transcription levels in individual oysters. Briefly, the expression values of OsHV-1 ORF104 and oyster housekeeping gene Elongation factor 1-alpha (El1α) were compared using the delta Ct method [71]. One μg of total RNA per sample was used to prepare sncRNA and mRNA libraries. In particular, sncRNA libraries were prepared according to TruSeq Small RNA Library Prep Kit (Illumina), whereas mRNA libraries were based on poly(A) RNA-selection or rRNA-depletion. Both mRNA library types were prepared according to the Illumina paired-end technology to retain the strand information of the reads (stranded libraries). The poly(A) library was prepared using the mRNA-Seq library Prep Kit V2 (Lexogen, Austria) following the manufacturer instructions, while the Ribo-0 library was prepared according to the TruSeq Stranded Total RNA kit (using a Ribo-zero Gold kit, Illumina, US). HT-sequencing was outsourced and carried out on a HiSeq 2500 instrument (2 × 150, Admera Health, USA) or on a HiSeq High-Output v4 instrument (2 × 125, DNA Sequencing Center at Brigham Young University, USA).

Analysis of small non-coding RNA-seq data

High-throughput sncRNA data were trimmed using cutadapt v.1.18 [72] implemented into Trimgalore! v.0.6.0 (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/), setting the minimal quality threshold to PHRED25 after removing adaptor sequences. FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was used to verify the absence of sequencing adaptors as well as to verify the quality statistics of each dataset. The clean reads were uploaded in CLC Genomic Workbench v.12.0 (Qiagen, US), they were size selected (18–40 nt) and mapped on the oyster and OsHV-1 reference genomes (GCA_000297895 [49] and MG561751 (OsHV-1) [73], respectively) as well as on the predicted oyster miRNA precursors obtained from MirGeneDB (http://mirgenedb.org/) [44], applying 0.9 and 1 for similarity and length fraction parameters, respectively. The mapping graphs of the 151 oyster miRNA precursors were manually inspected to verify the presence of the minimal annotation criteria. These included the presence of coverage for both miRNA arms, absence of reads mapped in the surroundings of the annotated miRNAs and 5′ read homogeneity. Subsequently, the sncRNA reads were clustered per sample and only perfect matches and were counted. Only clusters with a minimal representation of 100 reads were further considered and annotated according to the MirGeneDB oyster entries. miRNA expression values were calculated per sample as number of mapped read and compared using a proportions-based test [74] with FDR-corrected p-values. miRNAs with p-value lower than 0.01 and absolute fold-change higher than 2 were considered as differentially expressed (DE-miRNAs).

Annotation of oyster UTRs and miRNA/mRNA expression analysis

To identify the 3′ untranslated regions (3′-UTRs) of C. gigas mRNAs, we mapped the mRNA reads on the oyster genome (oyster_v9) using the large gap mapper implemented in the CLC, that allowed the presence of gaps due to introns. The resulting mapping was used to extend the predicted genes (limited to CDS regions) to untranslated regions. 3′-UTRs longer than 30 nt were counted and used, together with the OsHV-1 genome, to verify the presence of putative miRNA targets using miranda [50], applying a conservative minimal score of 155 and a minimal energy of − 20. To investigate the interaction of miRNAs and mRNA we computed Pearson’s correlation coefficient for all possible combination of the expression values of all genes and all miRNAs with expression values to obtain a Null distribution of all possible correlations. We then took subsets of a) all the interactions predicted by miranda and b) the interactions involving only differentially expressed miRNAs among these predicted interactions. The distributions of both predicted sets were then evaluated against the null distributions to identify regions where predicted interactions deviated from the null distribution.

Selection of the best OsHV-1 reference genome

Trimgalore was used to trim RNA-seq reads, applying a minimal quality of PHRED30, a minimal read length of 80 bp, and only validated paired reads were considered. To determine the bivalve Malacoherpesviridae most suitable as a reference genome, the whole read dataset was mapped on three available OsHV-1-μvar genomes, applying 0.95 and 0.95 of similarity and length fraction parameters, respectively. Read mapping resulted in similar matches to the three genomes, while SNP calling identified 59, 83 and 80 variants, for MG561751 (Italy), KY242785 and KY271630 (France and Ireland) genomes, respectively. Based on few synonymous and non-synonymous SNPs and considering the geographical vicinity of the OsHV-1-PT isolate (Porto Tolle, North Adriatic Sea), we selected MG561751 genome as reference for the purposes of this study.

Analysis of poly(a) tail selection and rRNA ribosomal RNA depletion performances

Preliminary to proceed with high-throughput sequencing, we tested the performance of an rRNA depletion approach in comparison with poly(A)-tail selection on one sample (S6). To calculate the amount of C. gigas and OsHV-1 reads (defined as on-target reads), all clean reads were mapped to the oyster and OsHV-1 reference genomes using the CLC mapping tool (Qiagen, Denmark), setting 0.8 and 0.5 for the similarity and length parameters, respectively. Unmapped reads were collected and remapped on the reference genomes while allowing the presence of large gaps (introns or large structural variations), using the CLC large gap read mapping tool, and applying 0.9 and 0.9 for the similarity and length parameters, respectively. The remaining unmapped reads were de-novo assembled using the CLC assembler tool (with a minimal contig length of 200 bp, bubble and word sizes set to automatic) and the obtained contigs were subjected to ORF prediction with the transdecoder tool by applying default parameters [75]. According to ORF predictions, the contigs were preliminary classified into coding or putative non-coding transcripts and were blasted (blastx) against the NCBI nr-protein database (downloaded the 10th of Sept. 2018). Blast results were used for species assignment, sequence annotation and identification of possible lncRNAs using the Blast2GO suite [76]. The identification of conserved domains on the predicted proteins was carried out with HMMer v.3.1 based on the whole Pfam-A domain collection (cut-off E-value of 10− 5) [77]. Putative non-coding transcripts were further screened for the presence of conserved RNA structures, using the Rfam v.13.0 database [78] with Infernal v1.1 [79].

To estimate the level of rRNA depletion, the trimmed reads of poly(A) and Ribo-0 libraries were mapped on a reference oyster rRNA sequence, obtained concatenating the 13 rRNA sequences annotated in the C. gigas genome as well as other oyster rRNAs annotated as ‘hypothetical protein’, identified by blastn against the nr NCBI database. The level of strand specificity of the libraries was determined by mapping the trimmed reads of each library on 7 oyster housekeeping genes [80], selected for the absence of known antisense transcription (EKC19952, Ubiquitin-conjugating enzyme E2D2; EKC42233, S-phase kinase-associated protein 1; EKC41722, Heterogeneous nuclear ribonucleoprotein A2/B1; EKC37135, Heterogeneous nuclear ribonucleoprotein Q; EKC32788, Eukaryotic translation elongation factor 2; EKC23295, Glyceraldehyde-3 phosphate-dehydrogenase and EKC33063, Elongation factor-1α). The coverage graphs were manually inspected to exclude the possible presence of non-annotated antisense transcripts. Subsequently, the number of reads mapped in the sense direction over the total number of mapped reads was taken as an estimation of the strand-specificity of the library, according to [81]. Limited to viral reads, and to better compare viral expression profiles between the two libraries and between sense and antisense directions, the expression values were computed both as total number of mapped reads and as Reads Per Kilobase Million (RPKM) either mapping the reads on the whole genome or, separately, on each viral ORF. In line with other reports [82], we showed that poly(A) selection underestimates the expression levels of non-polyadenylated transcripts like ncRNAs and histones, but, on the contrary, it produces a higher amount of on-target reads. This latter point was the reason for using poly(A) selection for this study. However, we still considered ribosomal rRNA depletion as a valid alternative for future investigations, probably useful to reveal OsHV-1 non-polyadenylated transcription not easily detectable with poly(A) reads. Further details are reported in Additional file 7.

RNA-seq expression analysis

Expression profiles were computed by mapping all clean reads on the virus and host genomes, applying 0.8 for both the length and similarity parameters. Owing to the stranded libraries, reads were mapped using a strand constraint (either sense or antisense mapping). For C. gigas, the expression values were computed as Transcript Per Million (TPM) to normalize for the different sequencing yield [83], whereas for OsHV-1 we used Reads Per Kilobase Million (RPKM) values because of the high difference in the number of total mapped reads among samples. Genes were regarded as differentially expressed (DEG) if presenting an absolute fold change higher than 2 with an FDR corrected p-value lower than 0.01 (Baggerley’s test).

Validation of miRNA expression by RT-qPCR

The expression levels of 8 selected miRNAs were tested by RT-qPCR using the miRCURY LNA miRNA SYBR Green PCR kit (Qiagen). The LNA primers were purchased from catalogue products, if available, or designed using the GeneGlobe platform (https://geneglobe.qiagen.com/) (Additional file 2). First-strand cDNAs were synthesized by starting from 50 ng of total RNA of samples S1-S8 and using the miRCURY LNA RT Kit (Qiagen), according to manufacturer’s instructions with the following cycle: 60 min at 42 °C, for 5 min at 95 °C and immediately cooled down to 4 °C. The obtained cDNAs were mixed in a unique pool and 5 dilutions were prepared, from 1:10 to 1:200 to test primer efficiency in a preliminary RT-qPCR plate. All the designed primers showed a high efficiency (r2 > 0.95) when tested over serial cDNA dilutions, with the single exception of Cgi-Novel-19 primer (0.88 of efficiency). Final RT-qPCR reactions were carried out using 3 μl of 1:100 cDNAs in a 10 μl of final reaction mixture (5 μl of 2X Master Mix, 0.5 μl of Rox passive reference dye, 1 μl of primer, 0.5 μl of water). Amplification cycles were performed on an Applied Biosystems 7900HT Fast Real-Time PCR System in a MicroAmp Fast Optical 384-Well Reaction Plate (Life Technologies) as follows: 95 °C for 2 min and 40 cycles of 95 °C for 10 s and 56 °C for 1 min. At the end of the reaction, a dissociation curve analysis was performed to ascertain the primer specificity. Each qPCR assay was carried out in triplicate on the same plate for each primer. Two stable housekeeping miRNAs were chosen as reference (Cgi-mir-10-P2 and Cgi-mir-184-P7). The relative expression ratio of the selected target gene was based on the delta–delta Ct method (2 − ΔΔCt) [71].

Availability of data and materials

The raw reads are deposited in the NCBI Short Reads Archive under project accession IDs PRJNA484109. For comparative purposes, other RNA-seq datasets were retrieved from the NCBI SRA database (PRJNA423079).

Abbreviations

cDNA:

complimentary DNA

CDS:

Coding sequence

DE:

Differentially expressed

HaHV-1:

Haliotid herpesvirus-1

HTS:

High-throughput sequencing

miRNA:

microRNA

mRNA:

messenger RNA

ncRNA:

non-coding RNA

OsHV-1:

Ostreid herpesvirus-1

RPKM:

Reads Per Kilobase Million

RT-qPCR:

Quantitative reverse transcription PCR

sncRNA:

short non-coding RNA

TPM:

Transcript Per Million

UTR:

Untranslated region

References

  1. FAO. Fao yearbook. fishery and aquaculture statistics 2016. http://www.fao.org/3/i9942t/I9942T.pdf: FOOD & AGRICULTURE ORG; 2018.

    Google Scholar 

  2. Kroodsma DA, Mayorga J, Hochberg T, Miller NA, Boerder K, Ferretti F, et al. Tracking the global footprint of fisheries. Science. 2018;359(6378):904–8.

    PubMed  CAS  Google Scholar 

  3. Pernet F, Lupo C, Bacher C, Whittington RJ. Infectious diseases in oyster aquaculture require a new integrated approach. Phil. Trans. R. Soc. B. 2016;371:20150213. https://doi.org/10.1098/rstb.2015.0213.

  4. Arzul I, Corbeil S, Morga B, Renault T. Viruses infecting marine molluscs. J Invertebr Pathol. 2017;147:118–35.

    PubMed  Google Scholar 

  5. Barbosa Solomieu V, Renault T, Travers M-A. Mass mortality in bivalves and the intricate case of the Pacific oyster, Crassostrea gigas J Invertebr Pathol 2015;131:2–10.

  6. Davison AJ, Trus BL, Cheng N, Steven AC, Watson MS, Cunningham C, et al. A novel class of herpesvirus with bivalve hosts. J Gen Virol. 2005;86(Pt 1):41–53.

    PubMed  CAS  Google Scholar 

  7. Bai C, Gao W, Wang C, Yu T, Zhang T, Qiu Z, et al. Identification and characterization of ostreid herpesvirus 1 associated with massive mortalities of Scapharca broughtonii broodstocks in China. Dis Aquat Org. 2016;118(1):65–75.

    PubMed  CAS  Google Scholar 

  8. López Sanmartín M, Power DM, de la Herrán R, Navas JI, Batista FM. Experimental infection of European flat oyster Ostrea edulis with ostreid herpesvirus 1 microvar (OsHV-1μvar): mortality, viral load and detection of viral transcripts by in situ hybridization. Virus Res. 2016;217:55–62.

    PubMed  Google Scholar 

  9. Ren W, Chen H, Renault T, Cai Y, Bai C, Wang C, et al. Complete genome sequence of acute viral necrosis virus associated with massive mortality outbreaks in the Chinese scallop, Chlamys farreri Virol J 2013;10:110.

  10. Batista FM, Arzul I, Pepin J-F, Ruano F, Friedman CS, Boudry P, et al. Detection of ostreid herpesvirus 1 DNA by PCR in bivalve molluscs: a critical review. J Virol Methods. 2007;139(1):1–11.

    PubMed  CAS  Google Scholar 

  11. Chang PH, Kuo ST, Lai SH, Yang HS, Ting YY, Hsu CL, et al. Herpes-like virus infection causing mortality of cultured abalone Haliotis diversicolor supertexta in Taiwan. Dis Aquat Org. 2005;65(1):23–7.

    PubMed  Google Scholar 

  12. Corbeil S, Williams LM, McColl KA, Crane MSJ. Australian abalone (Haliotis laevigata, H. rubra and H. conicopora) are susceptible to infection by multiple abalone herpesvirus genotypes. Dis Aquat Org. 2016;119(2):101–6.

    PubMed  CAS  Google Scholar 

  13. Yoshino TP, Bickham U, Bayne CJ. Molluscan cells in culture: primary cell cultures and cell lines. Can J Zool 2013 1;91(6):391–404.

  14. Morga B, Faury N, Guesdon S, Chollet B, Renault T. Haemocytes from Crassostrea gigas and OsHV-1: a promising in vitro system to study host/virus interactions. J Invertebr Pathol. 2017;150:45–53.

    PubMed  CAS  Google Scholar 

  15. Ji A, Li X, Fang S, Qin Z, Bai C, Wang C, et al. Primary culture of Zhikong scallop Chlamys farreri hemocytes as an in vitro model for studying host-pathogen interactions. Dis Aquat Org. 2017;125(3):217–26.

    PubMed  CAS  Google Scholar 

  16. Rosani U, Venier P. Oyster RNA-seq data support the development of Malacoherpesviridae genomics. Front Microbiol. 2017;8:1515.

    PubMed  PubMed Central  Google Scholar 

  17. King WL, Siboni N, Williams NLR, Kahlke T, Nguyen KV, Jenkins C, et al. Variability in the composition of Pacific oyster microbiomes across oyster families exhibiting different levels of susceptibility to OsHV-1 μvar disease. Front Microbiol. 2019;10:473.

    PubMed  PubMed Central  Google Scholar 

  18. de Lorgeril J, Lucasson A, Petton B, Toulza E, Montagnani C, Clerissi C, et al. Immune-suppression by OsHV-1 viral infection causes fatal bacteraemia in Pacific oysters. Nat Commun. 2018 Oct 11;9(1):4215.

    PubMed  PubMed Central  Google Scholar 

  19. Rosani U, Varotto L, Domeneghetti S, Arcangeli G, Pallavicini A, Venier P. Dual analysis of host and pathogen transcriptomes in ostreid herpesvirus 1-positive Crassostrea gigas. Environ Microbiol. 2015 Nov;17(11):4200–12.

    PubMed  CAS  Google Scholar 

  20. Bai C-M, Rosani U, Xin L-S, Li G-Y, Li C, Wang Q-C, et al. Dual transcriptomic analysis of Ostreid herpesvirus 1 infected Scapharca broughtonii with an emphasis on viral anti-apoptosis activities and host oxidative bursts. Fish Shellfish Immunol. 2018;82:554–64.

  21. He Y, Jouaux A, Ford SE, Lelong C, Sourdaine P, Mathieu M, et al. Transcriptome analysis reveals strong and complex antiviral response in a mollusc. Fish Shellfish Immunol. 2015;46(1):131–44.

    PubMed  CAS  Google Scholar 

  22. Bai CM, Zhang SM, Li YN, Xin LS, Rosani U, Wang CM. Dual transcriptomic analysis reveals a delayed antiviral response of haliotis diversicolor supertexta against haliotid herpesvirus-1. Viruses. 2019;1(4):383. https://doi.org/10.3390/v11040383.

  23. Neave MJ, Corbeil S, McColl KA, Crane MSJ. Investigating the natural resistance of Blackfoot p-a%%KERN_ERR%%ua Haliotis iris to abalone viral ganglioneuritis using whole transcriptome analysis. Dis Aquat Org. 2019;135(2):107–19.

    PubMed  CAS  Google Scholar 

  24. Green TJ, Raftos D, Speck P, Montagnani C. Antiviral immunity in marine molluscs. J Gen Virol. 2015;96(9):2471–82.

    PubMed  CAS  Google Scholar 

  25. Bartoszewski R, Sikorski AF. Editorial focus: entering into the non-coding RNA era. Cell Mol Biol Lett. 2018;23(1):45.

    PubMed  PubMed Central  Google Scholar 

  26. Bartel DP. Metazoan MicroRNAs. Cell. 2018;173(1):20–51.

    PubMed  PubMed Central  CAS  Google Scholar 

  27. Liu W, Ding C. Roles of LncRNAs in viral infections. Front Cell Infect Microbiol. 2017;7:205.

    PubMed  PubMed Central  Google Scholar 

  28. Liu T-Y, Zhang Y-C, Lin Y-Q, Hu Y-F, Zhang Y, Wang D, et al. Exploration of invasive mechanisms via global ncRNA-associated virus-host crosstalk. Genomics. 2019;112(2):1643–50.

  29. Vishnoi A, Rani S. MiRNA biogenesis and regulation of diseases: an overview. Methods Mol Biol. 2017;1509:1–10.

    PubMed  CAS  Google Scholar 

  30. Bhaskaran M, Mohan M. MicroRNAs: history, biogenesis, and their evolving role in animal development and disease. Vet Pathol. 2014;51(4):759–74.

    PubMed  CAS  Google Scholar 

  31. Pal AS, Kasinski AL. Animal models to study MicroRNA function. Adv Cancer Res. 2017;135:53–118.

    PubMed  PubMed Central  CAS  Google Scholar 

  32. Lee RC, Feinbaum RL, Ambros V. The C. Elegans heterochronic gene lin-4 encodes small RNAs with antisense complementarity to lin-14. Cell. 1993;75(5):843–854.

  33. Pfeffer S, Zavolan M, Grässer FA, Chien M, Russo JJ, Ju J, et al. Identification of virus-encoded microRNAs. Science. 2004;304(5671):734–6.

    PubMed  CAS  Google Scholar 

  34. Naqvi AR. Immunomodulatory roles of human herpesvirus-encoded microRNA in host-virus interaction. Rev Med Virol. 2019;20:e2081.

    Google Scholar 

  35. He Y, Zhang X. Comprehensive characterization of viral miRNAs involved in white spot syndrome virus (WSSV) infection. RNA Biol. 2012 Jul;9(7):1019–29.

    PubMed  CAS  Google Scholar 

  36. Ren Q, Huang Y, He Y, Wang W, Zhang X. A white spot syndrome virus microRNA promotes the virus infection by targeting the host STAT. Sci Rep. 2015;5:18384.

    PubMed  PubMed Central  CAS  Google Scholar 

  37. Cui Y, Huang T, Zhang X. RNA editing of microRNA prevents RNA-induced silencing complex recognition of target mRNA. Open Biol 2015 ;5(12). Available from: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4703055/. [cited 2018 Jun 12].

  38. Jiao Y, Zheng Z, Du X, Wang Q, Huang R, Deng Y, et al. Identification and characterization of microRNAs in pearl oyster Pinctada martensii by Solexa deep sequencing. Mar Biotechnol. 2014;16(1):54–62.

    PubMed  CAS  Google Scholar 

  39. Tian R, Zheng Z, Huang R, Jiao Y, Du X. miR-29a participated in nacre formation and immune response by targeting Y2R in Pinctada martensii. Int J Mol Sci. 2015;16(12):29436–45.

    PubMed  PubMed Central  CAS  Google Scholar 

  40. Zheng Z, Du X, Xiong X, Jiao Y, Deng Y, Wang Q, et al. PmRunt regulated by pm-miR-183 participates in nacre formation possibly through promoting the expression of collagen VI-like and Nacrein in pearl oyster Pinctada martensii. PLoS One. 2017;12(6):e0178561.

    PubMed  PubMed Central  Google Scholar 

  41. Chen H, Zhou Z, Wang L, Wang H, Liu R, Zhang H, Song L. An invertebrate-specific miRNA targeted the ancient cholinergic neuroendocrine system of oyster. Open Biol. 2016;6:160059. https://doi.org/10.1098/rsob.160059.

  42. Chen H, Xin L, Song X, Wang L, Wang W, Liu Z, et al. A norepinephrine-responsive miRNA directly promotes CgHSP90AA1 expression in oyster haemocytes during desiccation. Fish Shellfish Immunol. 2017;64:297–307.

    PubMed  CAS  Google Scholar 

  43. Chen G, Zhang C, Jiang F, Wang Y, Xu Z, Wang C. Bioinformatics analysis of hemocyte miRNAs of scallop Chlamys farreri against acute viral necrobiotic virus (AVNV). Fish Shellfish Immunol. 2014;37(1):75–86.

    PubMed  CAS  Google Scholar 

  44. Fromm B, Domanska D, Hoye E, Ovchinnikov V, Kang W, Aparicio-Puerta E, et al. MirGeneDB 2.0: The metazoan microRNA complement Genomics; 2018. Available from: http://biorxiv.org/lookup/doi/10.1101/258749. [cited 2019 Aug 14].

  45. Fromm B, Billipp T, Peck LE, Johansen M, Tarver JE, King BL, et al. A uniform system for the annotation of vertebrate microRNA genes and the evolution of the human microRNAome. Annu Rev Genet. 2015;49:213–42.

    PubMed  PubMed Central  CAS  Google Scholar 

  46. Grundhoff A, Sullivan CS, Ganem D. A combined computational and microarray-based approach identifies novel microRNAs encoded by human gamma-herpesviruses. RNA. 2006;12(5):733–50.

    PubMed  PubMed Central  CAS  Google Scholar 

  47. Rosani U, Bai C-M, Maso L, Shapiro M, Abbadi M, Domeneghetti S, et al. A-to-I editing of Malacoherpesviridae RNAs supports the antiviral role of ADAR1 in mollusks. BMC Evol Biol. 2019;19(1):149.

    PubMed  PubMed Central  Google Scholar 

  48. Axtell MJ, Meyers BC. Revisiting criteria for plant MicroRNA annotation in the era of big data. Plant Cell. 2018;30(2):272–84.

    PubMed  PubMed Central  CAS  Google Scholar 

  49. Zhang G, Fang X, Guo X, Li L, Luo R, Xu F, et al. The oyster genome reveals stress adaptation and complexity of shell formation. Nature. 2012;490(7418):49–54.

    PubMed  CAS  Google Scholar 

  50. Enright AJ, John B, Gaul U, Tuschl T, Sander C, Marks DS. MicroRNA targets in Drosophila. Genome Biol. 2003;5(1):R1.

    PubMed  PubMed Central  Google Scholar 

  51. Tarver JE, Taylor RS, Puttick MN, Lloyd GT, Pett W, Fromm B, et al. Well-annotated microRNAomes do not evidence pervasive miRNA loss. Genome Biol Evol. 2018;10(6):1457–70.

    PubMed  PubMed Central  CAS  Google Scholar 

  52. Zhou Z, Wang L, Song L, Liu R, Zhang H, Huang M, et al. The identification and characteristics of immune-related microRNAs in haemocytes of oyster Crassostrea gigas. PLoS One. 2014;9(2):e88397.

    PubMed  PubMed Central  Google Scholar 

  53. Bao Y, Zhang L, Dong Y, Lin Z. Identification and comparative analysis of the Tegillarca granosa haemocytes microRNA transcriptome in response to cd using a deep sequencing approach. PLoS One. 2014;9(4):e93619.

    PubMed  PubMed Central  Google Scholar 

  54. Zhang Z, Pi J, Zou D, Wang X, Xu J, Yu S, et al. microRNA arm-imbalance in part from complementary targets mediated decay promotes gastric cancer progression. Nat Commun. 2019;10(1):1–16.

    Google Scholar 

  55. BASKERVILLE S, BARTEL DP. Microarray profiling of microRNAs reveals frequent coexpression with neighboring miRNAs and host genes. RNA. 2005;11(3):241–7.

    PubMed  PubMed Central  CAS  Google Scholar 

  56. Shu L, Zhang X. Shrimp miR-12 Suppresses White Spot Syndrome Virus Infection by Synchronously Triggering Antiviral Phagocytosis and Apoptosis Pathways. Front Immunol. 2017;8. Available from: https://www.frontiersin.org/articles/10.3389/fimmu.2017.00855/full. [cited 2019 Oct 23].

  57. Kaewkascholkul N, Somboonviwat K, Asakawa S, Hirono I, Tassanakajon A, Somboonwiwat K. Shrimp miRNAs regulate innate immune response against white spot syndrome virus infection. Dev Comp Immunol. 2016;60:191–201.

    PubMed  CAS  Google Scholar 

  58. Jaree P, Wongdontri C, Somboonwiwat K. White Spot Syndrome Virus-Induced Shrimp miR-315 Attenuates Prophenoloxidase Activation via PPAE3 Gene Suppression. Front Immunol. 2018;9. Available from: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6178132/. [cited 2019 Oct 23].

  59. Renault T, Faury N, Barbosa-Solomieu V, Moreau K. Suppression substractive hybridisation (SSH) and real time PCR reveal differential gene expression in the Pacific cupped oyster, Crassostrea gigas, challenged with Ostreid herpesvirus 1. Dev Comp Immunol. 2011;35(7):725–35.

    PubMed  CAS  Google Scholar 

  60. Xing J, Lin T, Zhan W. Variations of enzyme activities in the haemocytes of scallop Chlamys farreri after infection with the acute virus necrobiotic virus (AVNV). Fish Shellfish Immunology. 2008;25(6):847–52.

    PubMed  CAS  Google Scholar 

  61. Chen Y-H, Song F, Miao Y-T, He H-H, Lian Y-Y, Li X-C, et al. A novel Laccase gene from Litopenaeus vannamei is involved in the immune responses to pathogen infection and oxidative stress. Dev Comp Immunol. 2020;105:103582.

    PubMed  CAS  Google Scholar 

  62. Chandan RK, Singh AK, Patel S, Swain DM, Tuteja N, Jha G. Silencing of tomato CTR1 provides enhanced tolerance against tomato leaf curl virus infection. Plant Signal Behav. 2019 Mar 4;14(3):e1565595.

    PubMed  PubMed Central  Google Scholar 

  63. Mushegian A, Karin EL, Pupko T. Sequence analysis of malacoherpesvirus proteins: Pan-herpesvirus capsid module and replication enzymes with an ancient connection to “Megavirales”. Virology. 2018;513:114–28.

    PubMed  CAS  Google Scholar 

  64. Golyaev V, Candresse T, Rabenstein F, Pooggin MM. Plant virome reconstruction and antiviral RNAi characterization by deep sequencing of small RNAs from dried leaves. Sci Rep. 2019;9(1):1–10.

    Google Scholar 

  65. Swevers L, Liu J, Smagghe G. Defense mechanisms against viral infection in Drosophila: RNAi and non-RNAi. Viruses. 2018;01:10(5).

    Google Scholar 

  66. Waldron FM, Stone GN, Obbard DJ. Metagenomic sequencing suggests a diversity of RNA interference-like responses to viruses across multicellular eukaryotes. PLoS Genet. 2018 Jul;14(7):e1007533.

    PubMed  PubMed Central  Google Scholar 

  67. Rosani U, Shapiro M, Venier P, Allam B. A needle in a haystack: tracing bivalve-associated viruses in high-throughput transcriptomic data. Viruses. 2019;11(3):205. https://doi.org/10.3390/v11030205.

  68. Domeneghetti S, Varotto L, Civettini M, Rosani U, Stauder M, Pretto T, et al. Mortality occurrence and pathogen detection in Crassostrea gigas and Mytilus galloprovincialis close-growing in shallow waters (Goro lagoon, Italy). Fish Shellfish Immunol. 2014;41(1):37–44.

    PubMed  Google Scholar 

  69. Webb SC, Fidler A, Renault T. Primers for PCR-based detection of ostreid herpes virus-1 (OsHV-1): application in a survey of New Zealand molluscs. Aquaculture. 2007;272(1):126–39.

    CAS  Google Scholar 

  70. Segarra A, Baillon L, Faury N, Tourbiez D, Renault T. Detection and distribution of ostreid herpesvirus 1 in experimentally infected Pacific oyster spat. J Invertebr Pathol. 2016;133:59–65.

    PubMed  Google Scholar 

  71. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(−Delta Delta C(T)) method. Methods. 2001;25(4):402–8.

    CAS  Google Scholar 

  72. Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal. 2011;17(1):10–2.

    Google Scholar 

  73. Abbadi M, Zamperin G, Gastaldelli M, Pascoli F, Rosani U, Milani A, et al. Identification of a newly described OsHV-1 μvar from the North Adriatic Sea (Italy). J Gen Virol. 2018;99(5):693–703.

  74. Baggerly KA, Deng L, Morris JS, Aldaz CM. Differential expression in SAGE: accounting for normal between-library variation. Bioinformatics. 2003 Aug 12;19(12):1477–83.

    PubMed  CAS  Google Scholar 

  75. Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, et al. De novo transcript sequence reconstruction from RNA-seq using the trinity platform for reference generation and analysis. Nat Protoc. 2013;8(8):1494–512.

    CAS  Google Scholar 

  76. Götz S, García-Gómez JM, Terol J, Williams TD, Nagaraj SH, Nueda MJ, et al. High-throughput functional annotation and data mining with the Blast2GO suite. Nucleic Acids Res. 2008;36(10):3420–35.

    PubMed  PubMed Central  Google Scholar 

  77. Finn RD, Coggill P, Eberhardt RY, Eddy SR, Mistry J, Mitchell AL, et al. The Pfam protein families database: towards a more sustainable future. Nucleic Acids Res. 2016;44(D1):D279–85.

    CAS  Google Scholar 

  78. Kalvari I, Argasinska J, Quinones-Olvera N, Nawrocki EP, Rivas E, Eddy SR, et al. Rfam 13.0: shifting to a genome-centric resource for non-coding RNA families. Nucleic Acids Res. 2018;46(D1):D335–42.

    PubMed  CAS  Google Scholar 

  79. Nawrocki EP, Eddy SR. Infernal 1.1: 100-fold faster RNA homology searches. Bioinformatics. 2013;29(22):2933–5.

    PubMed  PubMed Central  CAS  Google Scholar 

  80. Du Y, Zhang L, Xu F, Huang B, Zhang G, Li L. Validation of housekeeping genes as internal controls for studying gene expression during Pacific oyster (Crassostrea gigas) development by quantitative real-time PCR. Fish Shellfish Immunol. 2013;34(3):939–45.

    PubMed  CAS  Google Scholar 

  81. O’Grady T, Cao S, Strong MJ, Concha M, Wang X, Splinter Bondurant S, et al. Global bidirectional transcription of the Epstein-Barr virus genome during reactivation. J Virol. 2014;88(3):1604–16.

    PubMed  PubMed Central  Google Scholar 

  82. Boone M, De Koker A, Callewaert N. Capturing the ‘ome’: the expanding molecular toolbox for RNA and DNA library construction. Nucleic Acids Res. 2018 Apr 6;46(6):2701–21.

    PubMed  PubMed Central  CAS  Google Scholar 

  83. Wagner GP, Kin K, Lynch VJ. A model based criterion for gene expression calls using RNA-seq data. Theory Biosci. 2013 Sep;132(3):159–64.

    PubMed  CAS  Google Scholar 

Download references

Acknowledgments

We thank the Province of Ferrara (Italy) for the access to the Goro lagoon water records.

Funding

This work is partially granted by the H2020 project VIVALDI (Scientific basis and tools for preventing and mitigating farmed mollusk diseases; GA 678589) of the European Commission. The funding body had no role in designing the experiments nor in the interpretation of results. Open Access funding provided by Projekt DEAL.

Author information

Authors and Affiliations

  1. Department of Biology, University of Padova, 35121, Padova, Italy

    Umberto Rosani & Paola Venier

  2. Coastal Ecology Section, AWI - Alfred Wegener Institute Helmholtz Centre for Polar and Marine Research, Wadden Sea Station Sylt, 25992, List, Germany

    Umberto Rosani & K. Mathias Wegner

  3. Istituto Zooprofilattico delle Venezie, Legnaro, Italy

    Miriam Abbadi & Giuseppe Arcangeli

  4. Centre for Shellfish Research & Department of Fisheries and Aquaculture, Vancouver Island University, Nanaimo, BC, V9R 5S5, Canada

    Timothy Green

  5. Yellow Sea Fisheries Research Institute, Chinese Academy of Fishery Sciences, Qingdao, 266071, China

    Chang-Ming Bai

  6. CRIM Laboratory - Delta Institute, Goro, Ferrara, Italy

    Edoardo Turolla

Authors
  1. Umberto Rosani
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Miriam Abbadi
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Timothy Green
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Chang-Ming Bai
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Edoardo Turolla
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Giuseppe Arcangeli
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. K. Mathias Wegner
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Paola Venier
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

UR and PV designed the experiments and drafted the manuscript, ET curated the field trial and samplings, GA provided laboratory organization and equipment for bivalve manipulation and virus diagnosis. MA curated the first processing of the sampled bivalves and analyzed viral DNA loads, UR quantified viral RNA loads, UR and MW analyzed miRNA data, UR and TG comparatively analyzed library data, PV analyzed environmental data and provided research funding, UR performed all other analyses, UR, CMB, TG, MW and PV discussed the results. All the authors read and approved the final manuscript.

Corresponding authors

Correspondence to Umberto Rosani or Paola Venier.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1.

miRNA expression values. miRNA ID according to MirGeneDB, average expression value, percentage of total miRNA expression, expression levels in the six samples and coefficient of variation of the expression values are reported for 132 oyster miRNAs. miRNAs reported in bold resulted differentially expressed (Table 2), while underlined miRNAs contributed less to 0.01% to the total expression.

Additional file 2.

RT-qPCR analysis of selected miRNAs.

Additional file 3

mRNA expression values. Gene ID, fold change, FDR-corrected p-value, gene description and the six expression values as TPM are reported for the 403 differentially expressed genes in the comparison between samples S1, S4 and S6 versus S2, S3 and S5. Gene directly involved in oyster antiviral pathway are highlighted in yellow, whereas other genes discussed in the text are highlighted in light blue.

Additional file 4.

OsHV-1 expression values. Gene ID, expression values in samples S1, S4, S6 plus in other 7 datasets are reported for the OsHV-1 genes.

Additional file 5.

RNA-seq coverage of the OsHV-1 genome obtained with Ribo-depleted and polyA-selected libraries generated from sample S6. The coverage graph along the 204 kb OsHV-1 genome was reported in a 0-3000x scale for the Ribo-0 library mapped in sense direction (A), in antisense direction (C) and for the poly(A) library in sense direction (D) and antisense direction (F). The red arrows (B and E) depicted the OsHV-1 ORF annotations. The blue rectangle highlighted the viral DNA polymerase, ORF100.

Additional file 6.

Oyster sampling data.

Additional file 7

Comparison of two mRNA enrichment methods for dual RNA-seq analysis of C. gigas infected with OsHV-1.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Rosani, U., Abbadi, M., Green, T. et al. Parallel analysis of miRNAs and mRNAs suggests distinct regulatory networks in Crassostrea gigas infected by Ostreid herpesvirus 1. BMC Genomics 21, 620 (2020). https://doi.org/10.1186/s12864-020-07026-7

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-020-07026-7

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us