Integrated RNA-seq and RNAi Analysis of the Roles of the Hsp70 and SP Genes in Red-Shell Meretrix meretrix Tolerance to the Pathogen Vibrio parahaemolyticus.

The "Wanlihong" Meretrix meretrix (WLH-M) clam is a new variety of this species that has a red shell and stronger Vibrio tolerance than ordinary M. meretrix (ORI-M). To investigate the molecular mechanisms responsible for the WLH-M strain's tolerance to Vibrio, we challenged clams with Vibrio parahaemolyticus and then assessed physiological indexes and conducted transcriptome analysis and RNA interference experiments. The mortality, tissue bacterial load, and hemocyte reactive oxygen species level of ORI-M were significantly higher than those of WLH-M, whereas the content and activity of lysozyme were significantly lower. Gene Ontology functional annotation analysis and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis revealed that immune and metabolic pathways were enriched in Vibrio-challenged clams. The expressions of the heat shock protein 70 (Hsp70) and serine protease (SP) genes, which are involved in antibacterial immunity, were significantly upregulated in WLH-M but not in ORI-M, while the expression of the kynurenine 3-monooxygenase gene, a proinflammatory factor, was significantly downregulated in WLH-M. RNA interference experiments confirmed that Hsp70 and SP downregulation could result in increased mortality of WLH-M. Therefore, we speculate that Hsp70 and SP may be involved in the antibacterial immunity of WLH-M in vivo. Our data provided a valuable resource for further studies of the antibacterial mechanism of WLH-M and provided a foundation for the breeding of pathogen-resistant strains.

Introduction

The clam Meretrix meretrix (Mollusca: Bivalvia) is a traditional marine aquaculture species that is widely distributed throughout the Indo-Pacific region (Li et al. 2010; Zhang et al. 2012; Yue et al. 2015). However, repeated mass mortalities in recent years have posed a serious threat to the Chinese mariculture industry and to coastal ecology (Tang et al. 2006; Kim et al. 2011). Although multiple factors are responsible for the high death rates in cultured M. meretrix, previous research has shown that Vibrio is the main pathogen that caused mass mortality (Yue et al. 2010, 2011). Invertebrates lack an adaptive immune system and cannot exert specific immune responses, and their defense consists of a variety of non-specific cellular and humoral components. M. meretrix has only its innate immune system to respond to invasion of pathogens and other stress factors (Medzhitov and Janeway 1997; Arala-Chaves and Sequeira 2000). To date, however, few studies using transcriptome and RNA interference (RNAi) have been conducted to identify the antibacterial mechanisms of M. meretrix.

Vibrio is a Gram-negative bacterium that is widely present in the marine environment, and Vibrio parahaemolyticus is one of the main disease-causing pathogens causing the diseases of marine organisms. Xiao et al. (2017) recently reported that the presence of V. parahaemolyticus can cause acute hepatopancreas necrosis disease in shrimp, and Zhang et al. (2019) found that it induced death of mud crabs by reducing the metabolites associate with energy biosynthesis and the innate immune system. Gong et al. (2021) reported that leukocytes, alkaline phosphatase, and complement were primarily involved in resisting high-dose infection of tongue sole Cynoglossus semilaevis by Vibrio vulnificus. In another study, two small RNA libraries from hemocytes of the pearl oyster Pinctada fucata infected with Vibrio alginolyticus were constructed and sequenced using high-throughput Illumina deep sequencing technology, and the results indicated that ten of the differentially expressed miRNAs might participate in the host immune response to infection (Wang et al. 2017).

The “Wanlihong” M. meretrix strain (named WLH-M, registration No.: GS-01–007-2014) is a new clam variety that is characterized by its red shell, fast growth, and strong stress resistance, which are traits we selected for in our laboratory over five consecutive generations. This strain’s growth rate, nutritional value, and tolerance to Vibrio are significantly higher than those of ordinary M. meretrix (ORI-M). To better understand the antibacterial ability and potential antibacterial mechanism of WLH-M, our specific goals were to (1) determine if WLH-M has stronger antibacterial ability than ORI-M based on immune indexes; (2) identify differentially expressed genes (DEGs) related to WLH-M antibacterial function in the transcriptome; and (3) use RNAi to test whether heat shock protein 70 (Hsp70) and serine protease (SP) genes are involved in the antibacterial reaction of WLH-M. Our results were added to the gene database of clams and provided a basis for further studies of the antimicrobial mechanisms of clams.

Material and Methods

Animal Collection and Acclimation

Our study was approved by the Animal Care and Use Committee of Ningbo Wanli University, and the M. meretrix was treated humanely. Adult ORI-M (shell length: 40.08 ± 3.27 mm; wet weight: 16.53 ± 4.25 g) and WLH-M (shell length: 42.06 ± 2.10 mm; wet weight: 20.97 ± 2.56 g) were collected from the same farm in Ningbo Bay, Ningbo, China (29°38′ N and 121°46′ E). Prior to the experiment, the clams were acclimated for 14 days in an 800 L of sand-filtered seawater tank (temperature, 21 ± 0.5 °C; pH, 8.06 ± 0.01; salinity, 23.00 ± 0.01‰). Clams were fed with the microalga Chlorella vulgaris, and the seawater in the tank was replaced once daily during the acclimation period.

Determination of Survival Rate and Bacterial Load

The V. parahaemolyticus strain (strain NO. ATCC17802, Solarbio, Beijing, China) selected for the bacterial infection experiment was cultured in the laboratory. The clams were randomly divided into two 200-L tanks, 300 ORI-M and 300 WLH-M. Each group was divided into four groups, a control group with nothing and three challenge groups with 106, 107, and 108 CFU/mL of V. parahaemolyticus, after which mortality was observed until 120 h. During Vibrio immersion infection, the seawater was changed every 24 h and fresh V. parahaemolyticus was added. The concentration of V. parahaemolyticus was counted with blood cell counting plate (three repetitions) before challenge.

The total RNA of V. parahaemolyticus in gills and hepatopancreas of 10 ORI-M and 10 WLH-M after 96 h of challenge at three concentrations were extracted. Following Chen et al. (2011) and Guo et al. (2014), we used the heat resistant direct hemolysin (TDH) gene as the pathogenic gene for detecting enrichment of V. parahaemolyticus in M. meretrix, and pvuA was used as the internal reference gene (Table 1). The relative expression of TDH in ORI-M and WLH-M was detected by quantitative real-time polymerase chain reaction (qRT-PCR) analysis.

Table 1

Primers used in the qRT-PCR and RNAi

Gene name	Sequence (5′-3′)
Quantitative RT-PCR (qRT-PCR)
TDH-F	GTAAAGGTCTCTGACTTTTGGAC
TDH-R	TGGAATAGAACCTTCATCTTCACC
pvuA-F	CAAACTCACTCAGACTCCA
PvuA-R	CGAACCGAT TCAACACG
KMO-F	CCTGGATACCCATCTACAC
KMO-R	GCATACACGATCAACCCT
Hsp70-F	ACGGGTTGGACTGGAAT
Hsp70-R	TCTCCTGCGTTGGTTTG
SP-F	CCACCCTATGGATTTGTC
SP-R	TACGCTGTCACCTCTGCTTAT
β-actin-F	TTGTCTGGTGGTTCAACTATG
β-actin-R	TCCACATCTGCTGGAAGGTG
RNA interference (RNAi)
Hsp70-F	GGUACGUGUGUAAGGACAUTT
Hsp70-R	AUGUCCUUACACACGUACCTT
SP-F	GGUUGGAUGCAAAGCUCUUTT
SP-R	AAGAGCUUUGCAUCCAACCTT

Estimation of the Intracellular Reactive Oxygen Species Content in Hemocytes

Twenty individuals were taken from the control and experimental groups of ORI-M and WLH-M, respectively, to estimate the intracellular reactive oxygen species (ROS) concentration in hemocytes after challenge to 107 CFU/mL of V. parahaemolyticus for 96 h (the clams were from the “Determination of Survival Rate and Bacterial Load” section experiment). Following the protocol provided by the kit (Nanjing Jiancheng Bio, Inc., Jiangsu, China), we collected 150 μL of hemocytes from each sample and incubated it with fluorochrome (provided by the kit) at 20 °C in the dark for 45 min. Blood smears were washed three times with phosphate-buffered saline (PBS) and subsequently used to determine the ROS-specific fluorescence intensity using a fluorescence spectrophotometer (Spark, Tecan, Männedorf Switzerland) at an excitation/emission wavelength of 488/525 nm. According to the methods of Zha et al. (2019) and Guo et al. (2017), the relative content of intracellular ROS was calculated by dividing the fluorescence intensity of the sample by that of the control.

Measurement of the Concentration and Activity of Lysozyme in the Gill

We measured the enzymatic activity and concentration of lysozyme (LZM) in the gill tissue samples using bivalves LZM ELISA kits FK-97441 and FK-97442, respectively (FKBIO, Shanghai, China). Twenty individuals were taken from the control group and the experimental group (107 CFU/mL V. parahaemolyticus challenge for 96 h) of ORI-M and WLH-M, respectively. The gill tissue from each sample was homogenized with PBS on ice and centrifuged at 12,000 × g for 10 min at 4 °C, and then the supernatant was collected (the clams were from the “Determination of Survival Rate and Bacterial Load” section experiment). Ten microliters of supernatant was mixed with 40 μL of diluent in a microwell plate, and the mixture was incubated at 37 °C for 30 min following the manufacturer’s protocol. After the microwells were washed five times with a wash buffer, 50 μL of the conjugate reagent was added, followed by incubation at 37 °C for another 30 min. The microwell was washed with the wash buffer again before adding 50 μL each of chromogenic reagents of A and B. After 15 min of mixing in the dark, stop buffer was added to terminate the chromogenic reaction. The values of absorbance were measured at the wavelength of 450 nm using a microplate reader (Multiskan GO, Waltham, MA, Thermo). The concentration and activity of LZM of each sample were subsequently determined by referring to corresponding standard curves.

Sample Collection, RNA Isolation, and Illumina Sequencing

The clams were randomly divided into two 80-L tanks, 50 ORI-M and 50 WLH-M. Each group was divided into two groups, a control group and a challenge group with 107 CFU/mL of V. parahaemolyticus. After 24 h, gills and hepatopancreas were taken from the four groups, immediately frozen in liquid nitrogen. The gills of three replicates and hepatopancreas of three replicates in each group were preserved in dry ice and sent to Novogene Co., Ltd. (Beijing, China) for transcriptome sequencing. The remaining samples were stored at − 80℃ for total RNA extraction.

Total RNA was extracted from gill and hepatopancreas tissues of ORI-M and WLH-M in the experimental and control groups using the TRIzol method (Sangon Biotech, Shanghai, China). The degradation and concentration of the total RNA were checked using gel electrophoresis on 1% agarose gels and the NanoPhotometer® spectrophotometer (Implen, Westlake Village, CA, USA), respectively. RNA integrity was assessed using the RNA Nano 6000 Assay Kit with the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). All RNA of ORI-M and WLH-M was performed on the Illumina hiseq2000 platform of Novogene Co., Ltd. (Beijing, China). According to the high-quality data, the longest transcript unigene obtained by splicing was used as the reference sequence for further expression quantification, functional annotation, and other analysis. Following Sanhueza (2019), we added four fluorescent-labeled dNTPs, DNA polymerase, and linker primers to the sequenced flow cell for amplification. When each sequencing cluster extended the complimentary chain, each fluorescent-labeled dNTP released the corresponding fluorescence. The sequencer captured the fluorescent signal and converted the optical signal into a sequencing peak through computer software to obtain the fragments’ sequence information.

DEGs and Enrichment Analysis

The fragments per kilobase of exon per million fragments mapped (FPKM) was calculated based on the mapped unigene fragments using RNA-seq by Expectation-Maximation, which is the most commonly used method for measuring gene expression (Trapnell et al. 2010; Li and Dewey 2011). DESeq provides statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution, with the variance and the mean linked by local regression (Anders and Huber 2010). In this study, we used the two calculation methods to screen for DEGs for which expression changed significantly in WLH-M but not in ORI-M after V. parahaemolyticus challenge. Based on the DEGs’ function annotation, we then screened for DEGs related to antibacterial immunity.

Gene Ontology (GO) enrichment analysis of the DEGs was implemented using the GOseq R package-based Wallenius non-central hyper-geometric distribution, which can adjust for gene length bias in DEGs (Young et al. 2010). Kyoto Encyclopedia of Genes and Genomes (KEGG) is a database that can be used to understand high level functions and utilities of the biological system, such as the cell, organism, and ecosystem, based on molecular level information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies (Kanehisa et al. 2008).

Gene Expression Analysis

We measured the expression of three DEGs (kynurenine 3-monooxygenase (KMO), Hsp70, and SP) (Table 1) in the gill and hepatopancreas tissues of ORI-M and WLH-M at each time point using the Quantagene q225 Quantitative Real-Time PCR System (Kubo, Beijing, China). All primers (listed in Table 1) were designed by primer premier 5 and synthesized by TsingKe Biotech (Hangzhou, China). The analysis was carried out in a 10 µL reaction volume containing 1.0 µL of cDNA, 5.0 µL of 2 × ChamQ Universal SYBR Master Mix, 0.5 µL of each primer (10 pmol/µL), and 3.0 µL of RNase Free dH2O. The qRT-PCR cycling process was conducted as follows: one cycle of 95 °C for 3 min followed by 40 amplification cycles of 95 °C for 10 s, 60 °C for 30 s, and 95 °C for 15 s. We used the β-actin gene as the internal reference to calculate the relative expression levels of the tested genes (Jiang et al. 2017). Specificity of primers was confirmed by the dissociation curve of the amplification products at the end of each PCR.

Effect of RNAi on WLH-M After V. parahaemolyticus Challenge

To assess the function of the Hsp70 and SP genes in the antibacterial process of WLH-M, these genes were knocked down in vivo via double-stranded RNA (dsRNA)-mediated RNAi. Healthy WLH-M averaging 19.53 ± 1.12 g in weight were obtained from our laboratory and acclimatized in a lattice culture box with aeration for 7 days. The seawater temperature was 21 ± 0.5 °C and the salinity was 23‰ during the experiments. The clams were fed with C. vulgaris powder, and the seawater and the V. parahemolyticus were changed daily. For the RNAi experiment, 120 clams were divided randomly into five groups (Hsp70 interference, SP interference, PBS, no injection, and no injection/no challenge control). For each individual in the two interference groups, 20 µL (20 µM/L) of the gene-specific dsRNA was injected into the adductor muscle after anesthesia by eugenol, and the treated clams were cultured with 107 CFU/mL V. parahaemolyticus (Table 1). Un-injected clams and clams injected with 20 µL of PBS prior to V. parahaemolyticus challenge were used as control groups, and the no injection/no challenge group served as the time 0 control.

After 48 h and 96 h of exposure, five individuals were randomly chosen from each experimental group for gene expression analysis. Total RNA of each individual was extracted from the hepatopancreas tissue with TRIzol. The RNA quality and quantity were detected using NanoReady micro-volume spectrophotometers (LifeReal, Hangzhou, China). We used the HIScript®III RT SuperMix for qRT-PCR (+ gDNA wiper) kit (Vazyme, Nanjing, China) to purify total RNA and synthesize first-strand cDNA. Expression levels of Hsp70 and SP in the hepatopancreas at each time point were measured with the Quantagene q225 Quantitative Real-Time PCR System. Table 1 lists the specific primers used for qRT-PCR.

Statistical Analysis

All data are expressed as the means ± standard deviation from at least three sets of independent experiments. Data were processed using GraphPad Prism 8.0.1 and Origin 2021 statistical software, and significance among groups were determined by t-test. The relative expression level of genes was calculated by the 2–ΔΔCT method (Livak and Schmittgen 2001). p-values < 0.05 were considered to be statistically significant.

Results

Effect of V. parahaemolyticus Exposure on the Mortality and Bacterial Load of Clams

The mortality of both ORI-M and WLH-M increased as bacterial concentration increased, but the survival rate of ORI-M was lower than that of WLH-M at all three pathogen concentrations (Fig. 1A, B). The qRT-PCR detection showed that the bacterial load in the gills and hepatopancreas of ORI-M and WLH-M increased as the concentration of V. parahaemolyticus increased and that the bacterial load in ORI-M was significantly higher than that of WLH-M (Fig. 1C, D).

Fig. 1

The change of ORI-M and WLH-M with different concentrations of V. parahaemolyticus challenge. A, B Survival rate. C, D Bacterial load. *Significance is p < 0.05

Impact of V. parahaemolyticus Exposure on the Intracellular ROS Content of Clam Hemocytes

We measured ROS level to assess the effect of V. parahaemolyticus on free radical production in the hemocyte of ORI-M and WLH-M. Under V. parahaemolyticus stress, the ROS level in ORI-M hemocytes increased rapidly to a value that was 2.12 times higher than that of the control group (p < 0.05). However, the ROS content in WLH-M hemocytes did not change significantly after V. parahaemolyticus challenge (Fig. 2).

Fig. 2

Comparison of ROS content in hemocyte of ORI-M and WHL-M. VP indicates V. parahaemolyticus challenge. *Significance is p < 0.05

Effect of V. parahaemolyticus Challenge on LZM Content and Activity of Clams

The content and enzymatic activity of LZM in the gills of ORI-M and WLH-M after V. parahaemolyticus challenge are shown in Fig. 3A, B, respectively. When the clams were challenged by V. parahaemolyticus, the LZM content and activity of WLH-M were significantly higher than those of ORI-M, which were 1.7 times and 1.8 times higher than those of ORI-M, respectively (p < 0.05).

Fig. 3

Comparison of lysozyme content and activity of the ORI-M and WLH-M. A Lysozyme content. B Lysozyme activity. Note that the VP indicates V. parahaemolyticus challenge. *Significance is p < 0.05

Transcriptome Sequencing and Assembly

To identify the genes that were differentially expressed between ORI-M and WLH-M after V. parahaemolyticus challenge, we conducted RNA-seq and compared the overall transcriptome expression in the gills and hepatopancreas of ORI-M and WLH-M before and after challenge. We produced eight RNA-seq libraries using gill and hepatopancreas tissues from ORI-M and WLH-M from the treatment and control groups to generate a reference transcriptome. All raw reads were deposited to the Sequence Read Archive of the NCBI (accession number PRJNA821346). Library sequencing yielded a total of 564,859,660 raw reads from the eight libraries. After filtering out low-quality reads, 554,965,199 high-quality reads were retained and de novo assembled. We used Trinity software to select the longest transcript of each gene as the unigene of the eight libraries. The average GC content was 38.77%, the average values of Q20 and Q30 were 98.55% and 95.37%, and the sequencing error rate was only 0.02%. After assembly, 535,748,066 transcripts and 181,536,600 unigenes were obtained (Supplementary Table S1).

Functional Annotation and Classification of Unigenes and Identification and Analysis of DEGs

Because M. meretrix is not a model animal, we used seven databases to annotate the unigenes, thus providing more effective information for our research. Out of the 181,536,600 unigenes, only 193,622 were annotated. The numbers of unigenes that could be annotated in the different databases were 41,667 in NCBI non-redundant protein sequences (NR, 21.51%), 12,493 in NCBI nucleotide sequences (NT, 6.45%), 16,545 in Kyoto Encyclopedia of Genes and Genomes (KO, 8.54%), 25,565 in a manually annotated and reviewed protein sequence database (SwissProt, 13.2%), 44,293 in Protein family (PFAM, 22.87%), 44,286 in GO (22.87%), and 11,864 in euKaryotic Ortholog Groups (6.12%) (Supplementary Table S2).

In the gills of ORI-M, we detected 4560 DEGs after challenge with V. parahaemolyticus (1,484 up-regulated and 3,076 down-regulated genes), whereas we detected 792 DEGs (240 upregulated and 552 downregulated genes) in the gills of WLH-M after exposure to the pathogen. For the hepatopancreas, we found 312 DEGs (142 upregulated and 170 downregulated genes) in the ORI-M group and 212 DEGs (125 upregulated and 87 downregulated genes) in the WLH-M group after challenge with V. parahaemolyticus (Supplementary Table S3). The Venn diagram (Supplementary Fig. 1) and volcano plots (Supplementary Fig. 2) can intuitively show the DEGs.

Functional Analysis of DEGs

To further understand the changes in the patterns of transcriptional expression in ORI-M and WLH-M after V. parahaemolyticus challenge, we subjected the DEGs to GO enrichment analysis, which included placing the DEGs into three major functional categories: biological processes, cellular components, and molecular functions.

Based on the comparison between the gills of ORI-M and WLH-M, the most DEGs found in the biological processes category were related to biosynthetic process, organic substance biosynthetic process, and transmembrane transport; those in the cellular component category were integral component of membrane, intrinsic component of membrane, and cell projection; and those in the molecular functions category were ligase activity, lyase activity, and carboxy-lyase activity (Fig. 4A, B). The comparison between the hepatopancreas of ORI-M and WLH-M showed that the most DEGs found in the biological processes category were related to cell adhesion, biological adhesion, cellular component organization, and cellular component assembly; those in the cellular component category were cation binding, metal ion binding, and extracellular region; and those in the molecular functions category were chitin binding, protein transporter activity, and P-P-bond-hydrolysis-driven protein transmembrane transporter activity (Fig. 4C, D).

Fig. 4

GO annotation of differential expression genes between the ORI-M and WLH-M. A, B The comparison between the gills and hepatopancreas of ORI-M in the control group and experimental group. C, D The comparison between the gills and hepatopancreas of WLH-M in the control group and the experimental group. Note that the BP indicates biological process, CC indicates cellular component, MF indicates molecular function

We then classified the function of annotated genes in the gills of pathogen-exposed ORI-M and WLH-M using the KEGG database. The 20 most significantly enriched KEGG pathways were analyzed to identify those that play important roles in immune- and metabolic-related signaling pathways. The four most enriched pathways were related to cancer, toxoplasmosis, apoptosis-multiple species, and ubiquitin-mediated proteolysis (Fig. 5A, B). The pathways most enriched in the hepatopancreas of ORI-M and WLH-M were involved in in hippo signaling pathway-multiple species, hippo signaling pathway-fly, notch signaling pathway, and amino sugar and nucleotide sugar metabolism (Fig. 5C, D).

Fig. 5

KEGG pathway of differential expression genes between the ORI-M and WLH-M. A, B The comparison between the gills and hepatopancreas of ORI-M in the control group and the experimental group. C, D The comparison between the gills and hepatopancreas of WLH-M in the control group and the experimental group. Note that the color depth of the circle indicates the size of the q-value. The smaller the q-value is, the closer the color is to red. The size of the dots indicates how many genes are differentially expressed in each pathway

Screening of Candidate DEGs and Verification by qRT-PCR

According to the FPKM expression level, we identified 16 DEGs with significant differential expression in WLH-M but not in ORI-M after V. parahaemolyticus challenge (Supplementary Table S4). We then selected the DEGs KMO, Hsp70, and SP, which are related to antimicrobial immunity based on GO analysis, KEGG enrichment analysis, and database annotations, for qRT-PCR analysis to confirm that they were enriched in WLH-M after pathogen challenge. The expression levels of these three DEGs changed significantly at different time points after V. parahaemolyticus challenge (p < 0.05).

As shown in Fig. 6, KMO was specifically expressed in the gills of ORI-M. At 48, 72, and 96 h post-challenge, the relative expression level of KMO in ORI-M was significantly higher than that of WLH-M. Hsp70 and SP were also specifically expressed in the hepatopancreas of ORI-M; the relative expression level of Hsp70 in WLH-M was significantly higher than that in ORI-M at 24, 48, and 72 h, and the relative expression level of SP in WLH-M was significantly higher than that of ORI-M at 24, 48, 72, and 96 h post-challenge. Thus, the qRT-PCR analysis confirmed that these three DEGs were involved in the host immune response during the stress of V. parahaemolyticus exposure, indicating that these genes might play a role in the antibacterial response of WLH-M.

Fig. 6

The qRT-PCR validation of KMO, Hsp70, and SP. A The expression of KMO gene in gills. B, C The expression of Hsp70 and SP genes in hepatopancreas. *Significance is p < 0.05

Mortality of WLH-M After Hsp70 an SP RNAi

After RNAi, the expression levels of Hsp70 and SP changed over time. Compared with the time 0 control group, expression of Hsp70 and SP in the no injection group showed a significant upward trend after infection with V. parahaemolyticus, whereas no significant change in expression was detected in the PBS injection group after challenge. In both RNAi experimental groups, the expressions of Hsp70 and SP genes were downregulated significantly at 48 h and 96 h (p < 0.05) (Fig. 7A, B). The survival rates of WLH-M in the above five groups were 100%, 95%, 95%, 65%, and 80%, respectively (Fig. 7C).

Fig. 7

The change of WLH-M after RNAi. A, B Relative expression of Hsp70 and SP. C Survival rate of WLH-M. *Significance is p < 0.05

Discussion

WLH-M is a new clam variety with a red shell that was generated by selective breeding in our lab. Compared with ordinary clams, WLH-M exhibits stronger Vibrio tolerance and a high rate of survival in clam aquaculture. Its stronger tolerance to V. parahaemolyticus was confirmed by bacterial challenge in the laboratory. However, the molecular mechanism by which WLH-M responds to Vibrio infection is poorly understood, although some immune response characteristics have been reported in M. meretrix and other clams after bacterial challenge (Liang et al. 2017; Ren et al. 2017).

In our study, we exposed ORI-M and WLH-M to V. parahaemolyticus via immersion challenge, which mimicked natural infection conditions in seawater. The death rate and tissue bacterial load (gill and hepatopancreas) of ORI-M were significantly higher than those of WLH-M under different concentrations of V. parahaemolyticus (p < 0.05). With the increase of V. parahaemolyticus concentration in water, the bacterial load and mortality of ORI-M and WLH-M increase, which means that the group with higher Vibrio concentration has more Vibrio adhesion, invasion, and proliferation in the host. Some researchers found in the study of pathogen infection of M. meretrix and Litopenaeus vannamei that different concentrations of Vibrio attack will have a great impact on the mortality of the host, and high concentration of Vibrio stimulation will advance the time of large-scale death of the host, and lead to higher mortality (Wang et al. 2019; Phuoc et al. 2009).

We found that the ROS level of ORI-M was significantly higher than that of WLH-M (p < 0.05). Some previous studies reported that the excessive production of oxygen radicals in organisms will lead to lipid peroxidation, DNA damage, biofilm damage, and even damage to the integrity of tissues (Storey 2006; Yang et al. 2015). Under normal physiological conditions, antioxidant enzymes such as superoxide dismutase and catalase in the antioxidant defense system of mollusks can help remove ROS, enhance their antioxidant and immune capacity, and protect the body from oxidative stress (Nagalakshmi and Prasad 1998). When ORI-M were challenged with V. parahaemolyticus, the bacterial toxin caused the clams to begin the ROS stress response and produce a large amount of ROS. In contrast, the antioxidant enzyme system in WLH-M reduced the level of ROS to a value lower than that of ORI-M. This observation is consistent with previous results showing that ROS content of Chlamys farreri increased significantly at 24, 48, and 72 h after Vibrio anguillarum challenge (Fan et al. 2007). Ingestion of toxin-producing dinoflagellates by the scallop Patinopecten yessoensis and the mussel Mytilus edulis can lead to ROS oxidative damage and an increase in ROS concentration (Ma 2012).

LZM is a widely distributed hydrolase that has lytic activity against bacterial peptidoglycan, which enables it to protect the host against pathogenic infection (Li et al. 2009). The activity of enzymes is an important index of mollusk physiology and immune status. The invertebrate-type LZM plays an important role in the immunity of invertebrates, including bivalves (Itoh et al. 2007; Xue et al. 2007; Kuwano et al. 2013). Gagnaire et al. (2006) reported that the activity of LZM in the blood cells of the oyster Crassostrea gigas decreased significantly under the stress of exposure to polycyclic aromatic hydrocarbons. Zhao et al. (2007) found that the LZM identified in the scallop Chlamys farreri was an important component of the scallop antibacterial system and that its sequence polymorphism was associated with resistance to diseases caused by the bacterium Listonella anguillarum (Li et al. 2021). In the current study, we found that the content and activity of LZM in the gills of WLH-M after V. parahaemolyticus challenge were significantly higher than those of ORI-M (p < 0.05). This indicated that WLH-M had stronger antibacterial and immune ability than ORI-M.

Carotenoids also may be involved in the enhanced immunity of WLH-M compared with ORI-M. WLH-M contain more carotenoids than M. meretrix with other shell colors, and the higher concentration of this pigment is one of the causes of its red shell. We previously extracted carotenoids from M. meretrix strains by organic solvent extraction and found that their content was 1.5 times higher in WLH-M compared to ORI-M. Researchers also have reported that carotenoids could regulate the immune system to resist bacterial invasion because the antioxidants and antioxidant enzymes in carotenoids could inhibit the gene expression of inflammatory factors (Han 1999; Zhang and Diao 2021). Therefore, we speculate that high carotenoid contents might play a role in the stronger bacteria tolerance of WLH-M compared to ORI-M.

As an important application of second-generation sequencing technology, RNA-seq has become an essential approach for transcriptome analysis and quantitative analysis of gene expression in organisms (Wang et al. 2009; Yu et al. 2019). In this study, we constructed transcriptome libraries from gill and hepatopancreas tissues of WLH-M and ORI-M before and after infection by V. parahaemolyticus. Both GO annotations and KEGG pathways were enriched mainly in immune metabolism-related signaling pathways in Vibrio-challenged clams. We found that the expression of 16 genes significantly changed in WLH-M but not in ORI-M after challenge with V. parahemolyticus, and three of these DEGs (Hsp70, SP, and KMO) were related to antimicrobial immunity. The extensive transcriptional data obtained from the gills and hepatopancreas of ORI-M and WLH-M provided us with the necessary references for further studies of the antibacterial immune mechanism of clams.

The Hsp family constitutes the oldest defense system in all organisms ranging from bacteria to humans (Anirudhan et al. 2021). Hsp70 is known to have physiological and ecological importance in the response to pathogen infection and environmental stress. For example, Rungrassamee et al. (2010) reported that the expression levels of Hsp70 and Hsp90 in the prawn Penaeus monodon significantly increased after 3 h of exposure to Vibrio harveyi, and Wang et al. (2013) reported upregulation of Hsp70 after V. parahaemolyticus infection in adult bay scallops Argopecten irradians. Additionally, the expression of Hsp70 in the zebra mussel Dreissena polymorpha showed a time­dependent increase after lipopolysaccharide stimulation (Xu and Faisal 2009). Under stress, the expression level of Hsp70 family genes in cells increases, which can prevent cell atrophy and death by inhibiting the expression of stress kinase and apoptotic genes (Ming et al. 2010; Cellura et al. 2006; Campisi et al. 2003). Other studies have found that Hsp70 can also promote the production of endogenous antioxidant enzymes and reduce the production of ROS by inhibiting key enzymes of oxygen free radicals, such as NADPH reductase (Miao et al. 2007). Therefore, the significant upregulation of HSP70 gene expression in WLH-M also indicates that it inhibits oxygen free radicals in M. meretrix and reduces the production of ROS, thus preventing the damage of WLH-M under the stress of V. parahaemolyticus.

SPs are digestive enzymes found in marine bacteria (Li et al. 2016), marine fungi (Jasmin et al. 2010; Benmrad et al. 2019), and marine invertebrates (Fatma et al. 2017; Harish and Uppuluri 2018). When Drosophila are infected by fungi or Gram-positive bacteria, activation of SP makes Spaetzle protein become an active ligand, which binds to Toll and activates the Toll/dif pathway, leading to the synthesis and release of antimicrobial peptides (Ligoxygakis et al. 2002). There is a type of SP in the digestive tract of the silkworm Bombyx mori that has strong antiviral activity against Bombyx mori nuclear polyhedrosis virus and can play an antiviral role in the process of virus infection (Nakazawa et al. 2004). SPs play an important role in the non-specific immune defense system of shellfish. When pathogenic bacteria invade the body, the serine protease in the cell wall is released outside the cell in the form of the zymogen, which is activated after binding with a specific protease. The activated serine protease mediates phenol oxidase (Piao et al. 2005), which inhibits and kills pathogenic bacteria by identifying pathogenic surface molecules and producing enzyme cascade reactions such as agglutination, and antibacterial substances (lysozyme), and promoting blood cell phagocytosis (Rao et al. 2010). When Vibrio anguillarum infection can induce the expression of serine protease gene in scallops, it is found that the expression level of serine protease gene is significantly increased 16 h after infection, indicating that they participate in the immune response of scallops and play a role in immune defense (Zhu et al. 2006). In this study, after the stress of V. parahaemolyticus, the expression of serine protease in the hepatopancreas of WLH-M was significantly higher than that of ORI-M, indicating that serine protease participated in the immune response process of WLH-M and played a certain bactericidal role.

KMO is a key rate-limiting enzyme in canine urea metabolism to 3-hydroxylamine, and its increased expression can amplify the inflammatory response and lead to cell oxidative damage (Gao 2017). Inhibition of KMO has been recognized as a potential approach to reducing the production of neurotoxic kynurenine pathway metabolites and moderating the inflammatory response during infection and inflammation. Elevated plasma quinolinate levels in rhesus monkeys infected with simian immunodeficiency virus can be prevented by treatment with the KMO inhibitor CHDI-340246 (Swainson et al. 2019). Li et al. (2018) reported that BcKMO affected the virulence of the fungus Botrytis cinerea by regulating the activity of cell wall degradation enzymes, toxins, acid production, and expression of pathogenicity-related genes. Potential use of KMO inhibitors to modulate inflammation in patients with COVID-19 has also been reported (Collier et al. 2021).

These studies showed that all three genes were involved in inflammatory responses, so we selected them for qRT-PCR to further confirm their role in the immune response in WLH-M after bacterial challenge. We found that Hsp70 and SP were specifically expressed in the hepatopancreas and that their expression level was significantly upregulated in WLH-M after V. parahaemolyticus challenge compared to ORI-M. In contrast, KMO was specifically expressed in the gills and was significantly downregulated in pathogen-challenged WLH-M compared to ORI-M.

RNAi is a highly evolutionally conserved process of post-transcriptional gene silencing by which dsRNA, when introduced into a cell, causes sequence-specific degradation of homologous mRNA sequences (Almeida and Allshire 2005). Silencing effects could be used to control infectious disease by interfering with pathogen transmission, development, and proliferation within the host (Kang and Hong 2008). In recent years, the use of RNAi has emerged as an effective tool to study gene function in aquatic animals (Tang et al. 2019). Although the emergence and spread of Vibrio in the aquaculture setting has substantial negative economic and environmental impacts, the potential use of RNAi to develop treatment strategies against these diseases has received little attention. Wang et al. (2011) studied numerous aspects of the molluscan immune response against pathogens, and they used RNAi to confirm the existence of the C. farreri Toll-like receptor (TLR) signaling pathway and its involvement in bivalve immune defense. Additionally, they found a lower apoptosis level and higher cumulative mortality after bacterial challenge with Listonella anguillara. These results collectively indicate that several immune responses, including apoptosis, antibacterial activity, and antioxidant activity, can be activated by the TLR signaling pathway in scallops. Shen et al. (2018) injected mussels with melanocyte inducing transcription factor dsRNA and detected a decrease in expression of tyrosinase, which suggested that this transcription factor could positively regulate the expression of tyrosinase, which is an important enzyme involved in melanin synthesis. In our experiment, the expression levels of Hsp70 and SP were significantly downregulated after RNAi in WLH-M challenged with V. parahaemolyticus challenge, and their survival rate also decreased. The significant increase in mortality after RNAi and bacterial challenge indicated that these two genes might be involved in the antibacterial process of WLH-M.

Conclusion

In summary, the physiological indexes, and the transcriptome sequencing of gill and hepatopancreas of ORI-M and WLH-M after V. parahaemolyticus challenge were conducted in our study. The DEGs detected in this study, especially those related to the inflammatory reaction in M. meretrix, provided useful information for elucidating the molecular mechanisms responsible for the clam’s antibacterial response to Vibrio infection. RNAi results showed that two genes, Hsp70 and SP, were involved in the antibacterial immunity of WLH-M. Our findings enriched the M. meretrix transcriptome database and provided insight into the antibacterial immunity mechanism of WLH-M against V. parahaemolyticus.

Below is the link to the electronic supplementary material.

Supplementary file1 (DOCX 1014 KB)