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)
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.
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.
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.
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 . OsHV-1 is one of the two dsDNA viruses of the Malacoherpesviridae family , 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 . At present, the lack of mollusk cell lines  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 . 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 . ncRNAs are known as post-transcriptional regulators , 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 . 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 , miRNAs have been identified in many animals, plants, and viruses, suggesting independent miRNA origins throughout the major evolutionary lineages . Virus-encoded miRNAs have been reported in vertebrate herpesviruses, for instance in Epstein-Barr virus  where they perform immuno-modulatory functions . Regarding invertebrate viruses, miRNAs have only been reported in the White Spot Syndrome Virus (WSSV) . WSSV miRNAs target the JAK-STAT signaling pathway to interfere with the immune system of the shrimp , and their sequences are actively edited by host enzymes, like the adenosine deaminase acting on dsRNAs (ADAR-1) .
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  (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.
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).
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 , 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 . 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).
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.
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).
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).
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 . 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 . 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.
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 . 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 , 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  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.
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.
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 , 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  and Tegillarca granosa . 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 , 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 . 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 . 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 . Two other miRNAs, mir-375 and mir-750, were also highly responsive to WSSV infection in Panaeus monodon . Additionally, mir-315 regulated the pro-PO system during WSSV infection, thereby inhibiting the spread of the virus . 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 , 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 . 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 . 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 , although the biological meaning of this editing during oyster-OsHV-1 interaction is still unknown.
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 . 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 , 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 , 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 . RNAi was reported during WSSV infection in shrimp  as well as in other arthropod species , but this study rather supports previous findings indicative of a marginal role for RNA interference in the antiviral response of bivalve species [66, 67].
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).
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 . 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) . 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) . 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 . 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  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  and MG561751 (OsHV-1) , respectively) as well as on the predicted oyster miRNA precursors obtained from MirGeneDB (http://mirgenedb.org/) , 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  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 , 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 . 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 . 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) . Putative non-coding transcripts were further screened for the presence of conserved RNA structures, using the Rfam v.13.0 database  with Infernal v1.1 .
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 , 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 . 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 , 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 , 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) .
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).
Reads Per Kilobase Million
Quantitative reverse transcription PCR
short non-coding RNA
Transcript Per Million
FAO. Fao yearbook. fishery and aquaculture statistics 2016. http://www.fao.org/3/i9942t/I9942T.pdf: FOOD & AGRICULTURE ORG; 2018.
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.
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.
Arzul I, Corbeil S, Morga B, Renault T. Viruses infecting marine molluscs. J Invertebr Pathol. 2017;147:118–35.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Rosani U, Venier P. Oyster RNA-seq data support the development of Malacoherpesviridae genomics. Front Microbiol. 2017;8:1515.
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.
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.
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.
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.
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.
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.
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.
Green TJ, Raftos D, Speck P, Montagnani C. Antiviral immunity in marine molluscs. J Gen Virol. 2015;96(9):2471–82.
Bartoszewski R, Sikorski AF. Editorial focus: entering into the non-coding RNA era. Cell Mol Biol Lett. 2018;23(1):45.
Bartel DP. Metazoan MicroRNAs. Cell. 2018;173(1):20–51.
Liu W, Ding C. Roles of LncRNAs in viral infections. Front Cell Infect Microbiol. 2017;7:205.
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.
Vishnoi A, Rani S. MiRNA biogenesis and regulation of diseases: an overview. Methods Mol Biol. 2017;1509:1–10.
Bhaskaran M, Mohan M. MicroRNAs: history, biogenesis, and their evolving role in animal development and disease. Vet Pathol. 2014;51(4):759–74.
Pal AS, Kasinski AL. Animal models to study MicroRNA function. Adv Cancer Res. 2017;135:53–118.
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.
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.
Naqvi AR. Immunomodulatory roles of human herpesvirus-encoded microRNA in host-virus interaction. Rev Med Virol. 2019;20:e2081.
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.
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.
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].
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.
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.
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.
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.
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.
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.
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].
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.
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.
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.
Axtell MJ, Meyers BC. Revisiting criteria for plant MicroRNA annotation in the era of big data. Plant Cell. 2018;30(2):272–84.
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.
Enright AJ, John B, Gaul U, Tuschl T, Sander C, Marks DS. MicroRNA targets in Drosophila. Genome Biol. 2003;5(1):R1.
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.
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.
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.
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.
BASKERVILLE S, BARTEL DP. Microarray profiling of microRNAs reveals frequent coexpression with neighboring miRNAs and host genes. RNA. 2005;11(3):241–7.
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].
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.
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].
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.
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.
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.
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.
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.
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.
Swevers L, Liu J, Smagghe G. Defense mechanisms against viral infection in Drosophila: RNAi and non-RNAi. Viruses. 2018;01:10(5).
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.
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.
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.
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.
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.
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.
Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal. 2011;17(1):10–2.
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.
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.
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.
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.
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.
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.
Nawrocki EP, Eddy SR. Infernal 1.1: 100-fold faster RNA homology searches. Bioinformatics. 2013;29(22):2933–5.
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.
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.
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.
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.
We thank the Province of Ferrara (Italy) for the access to the Goro lagoon water records.
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.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
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.
RT-qPCR analysis of selected miRNAs.
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.
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.
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.
Oyster sampling data.
Comparison of two mRNA enrichment methods for dual RNA-seq analysis of C. gigas infected with OsHV-1.
About this article
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
- C. gigas