Transcriptomic Analysis of the Reproductive Effects on Immunity in the Clam Meretrix Petechialis During the Breeding Season


 Background: Immune defense and reproduction are both physiologically demanding processes, therefore the trade-off is expected between them in breeding season. Although this balance strategy has been extensively studied in many species, it is still rarely noticed in mollusc. Moreover, summer mortality of marine bivalve often occurs during breeding season and reproduction is suspected to play a significant role for the mass death. Results: To address if spawning affects immunity to cause clam death, we performed transcriptome and gene expression analyses in the clam Meretrix petechialis pre-/post-spawning. DEGs enrichment analysis revealed important immune signaling pathways and key genes changed after spawning. Further analysis showed females up-regulated genes involved in apoptosis, TLR signal pathway and heat shock, whereas males down-regulated complement-related genes after spawning. Additionally, both genders of clams up-regulated its immune response level to against Vibrio infection after spawning revealed by the changes of four immune-related DEGs. Conclusions: The spawning indeed affected the clam immune resistance. The immune and defense related genes and pathways obtained here provide important insights into the molecular mechanisms of physiological acclimation and immune response under reproduction-influence and offers some clues to explain the possible reason for clam summer mortality.


Background
In natural condition, the resources acquired by an organism are usually not enough to satisfy its needs in growth, reproduction and survival, which will lead to trade-off of resource allocation [1]. It means that organisms exibly invest their energy to different priorities [1][2][3][4][5]. One commonly observed trade-off occurs between reproduction and immune due to they are both physiologically demanding processes [6]. Recent works in this eld can be divided into two broad categories: the immune impacting on reproduction and reproductive efforts affecting immunity, which also allows us to understand this trade-off by a comprehensive perspective. Luong and Polak [7] reported the selected lines of Drosophila nigrospiracula showed a decreased fecundity after increased their resistance against an ectoparasite mite. Similar effects were also founded in other species such as insect and sh [8,9]. On the other hand, investigation on the reproductive cycle suggested that reproductive process can weak the immunocompetence for a seasonally reproducing sh [6]. Furthermore, the trade-off is also known to be affected by factors such as sex, body size, seasonal changes and parasite infection [10].
As we know, mass mortality of marine bivalve often occurs during the summer breeding season [11][12][13]. Spawning is considered to be a key cause of this result, because gametogenesis is suspected to be a period of intensive physiological change and most of the energy acquired is used for the production of gametes [14,15]. Although the trade-off between reproduction and immune has been deeply studied in the invertebrate model organisms, Drosophila and Caenorhabditis, it remains understudied in mollusc.
Fortunately, mollusc shares the conserved immune components and pathways with the above invertebrate model organisms, such as the immune effectors and immune signal transduction cascades, etc [16]. They all display a highly effective immune response that consists of interconnected cellular and humoral components, although lacking an adaptive immune response comparable to vertebrates [17][18][19].
To date, the association of summer mortality and immunity in mollusc has been extensively studied. For example, a review in oyster Crassostrea gigas summarized the omics approaches to investigate hostpathogen interactions in mass mortality outbreaks [20]. Besides, pathogene caused mass mortality were widely reported in scallop Chlamys farreri [21], clam Chamelea gallina [22] and abalone Haliotis gigante [23]. In our laboratory, we have investigated the immune response in the clam Meretrix petechialis to pathogenic bacteria associated with clam mass mortality in recent years [24,25]. However, the relationship study between reproduction and immunity from a global perspective remains rare in mollusc because the weak research background and the lack of feasible method. RNA-seq as an e cient and high-throughput technique is commonly used to investigate the transcription of genes including immunity, reproduction and growth, etc. [26][27][28], thus, it provides a possible way for the correlation research of reproduction and immunity [29].
In the present study, we selected the clam M. petechialis, an important commercial mollusc cultured in China, to investigate the changes of immune characteristics affected by production in breeding season by RNA-seq. Moreover, we further examined the reproductive effects on the clam immune response against pathogen invasion. These results will help us to understand the trade-off of reproduction and immune in breeding season, and provide insights into the causes of clam summer mortality.

Experimental animals
Matured clams were collected from one cultured population in Wenzhou, China. The clams were acclimated in a 2000 L tank with aerated seawater and fed with Isochrysis galbana for ve days at 28 °C. Then we took the pre-spawning individuals, identi ed the males and females by dissecting the gonads and grouped "pre_M" and "pre_F" for the male and female respectively. Parental clams were induced spawning by running seawater while post-spawning individuals were classi ed as "post_M" and "post_F" respectively, according to their germ cell morphology. Then the 9 clams from each group were dissected, and 3 hepatopancreas were mixed as one replicate and total of 3 replicates for RNA-Seq.

RNA-Seq and data analysis
For RNA-Seq, the hepatopancreas were sent to BGI (Shenzhen, China). After total RNA isolation, cDNA libraries preparation, RNA-seq was performed by BGISEQ-500 platform. Raw data (raw reads) were ltered by SOAPnuke, removing reads containing adapters, more than 5% unknown nucleotides and low-quality reads (more than 50% bases with Q-value ≤ 15) [30]. Then the assembled de novo transcriptome [24] was used as the reference database and the clean reads were aligned with the reference database using bowtie2 (v2.2.5) [31]. The expression level of genes and transcripts was calculated using RSEM [32]. For gene annotation, assembled unigenes were aligned with the TFs (AnimalTFDB), Nr (NCBI non-redundant protein sequences), GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) database.
Differentially expressed genes (DEGs) analysis Differential expression analysis between two groups, including pre_F vs. post_F, pre_M vs. post_M, pre_F vs. pre_M and post_F vs. post_M was performed by DESeq R package software (1.10.1). Adjusted P value (Q-values) of 0.01 and log2 (fold_change) of 2 were set as the threshold for signi cantly differential expression. DEGs were classi ed to signaling pathways by KEGG pathway analysis, while GO enrichment analysis was used to describe DEGs functions by mapping to functional categories. Based on the GO and KEGG annotation results and the o cial classi cation, the functions of DEGs were performed by the phyper R package, and enrichment was considered for FDR corrected p-values (Qvalue) less than 0.05.

Vibrio challenge and sample collection
The Vibrio challenge experiment was conducted with immersion method as described in Yue et al. [33] and Vibrio parahaemolyticus strain (MM21) isolated from M. petechialis was used in the challenge. Clams were divided into four groups as the groups in RNA-seq (pre_F, post_F, pre_M and post_M) and immersed in seawater with 1×10 7 CFU ml −1 V. parahaemolyticus. The seawater through sand lter and V. parahaemolyticus is renewed every day. The hepatopancreas tissues from 5 clams in each group were sampled at 0, 12, 24 hours post-immersed and then frozen in liquid nitrogen for RNA extraction.

Quantitative real-time PCR analysis
In order to verify the expression pattern obtained by RNA-seq, ten immune-related DEGs were randomly selected and detected by qPCR. The template used in the qPCR was the cDNA applied to RNA-seq. Primer premier 5 was used to design the primer pairs to amplify the selected genes. The sequences of primer pairs were listed in Additional le 1 Table S1. The mean value of β-actin and EF1α were employed as the internal references to normalize the relative expression levels among samples [34]. Bio Rad CFX 96 Real time PCR system of SYBR premix (TaKaRa, Japan) was used in qPCR and the PCR process were 95 °C for 30 s and 40 cycles of 95 °C for 5 s, 60 °C for 30 s, followed by the melting curve determination. Similarly, the expression patterns of four immune related DEGs were examined by qRT-PCR in the clams under Vibrio immersion challenge. The relative gene expression was analyzed using the 2 -△△CT method [35]. All data were expressed as the mean ± standard deviation and applied to F text. Through analysis by t-test using software SPSS version 22.0, the difference was considered statistically signi cant when P< 0.05.

Western blot analysis
The protein expression of MpDOUX and MpMITF were detected using Western blots. Protein Analysis Kit (BestBio, China) was used to extract the total protein from the hepatopancreases of clams. The polyclonal antibody used in this study included DOUX polyclonal antibody (catalog no. A8583, ABclonal, USA) and MITF polyclonal antibody (catalog no. A1255, ABclonal, USA). The detailed description of the western blot method was referred to Yue et al. [36]. In detail, the samples at each time point (12 h, 24 h) after the Vibrio challenge were mixed by protein from 3 individuals. 12% SDS-PAGE was used to separate each equal amounts of protein and then protein was transferred from gel onto a PVDF membrane. After blocked with 5% non-fat milk for 2 h, this membrane was incubated with the primary rabbit anti-DOUX or anti-MITF antibody or rabbit anti-β-actin antibody (ABclonal, USA) at 4℃ overnight. Then the membrane was washed 3 times for 10 min with TBST and incubated with the secondary goat anti-rabbit IgG-HRP antibody (ZSGB-BIO, China) for 1 h. Finally, after washing 3 more times, the protein on membrane was a visualization using Sparkjade ECL super (Cat # ED0015-B, Sparkjade).

Sequencing assembly and annotation
A total of 263.21M raw reads for reference database was generated based on BGISEQ-500 sequencer and 262.05M clean reads were remained after ltering out the low quality reads. For the clean reads, more than 97.89% of the bases had a phred value > 20, and more than 89.35% of them had a phred value > 30. An overview of the assembly results was provided in Additional le 2 Table S2. The raw data were uploaded to the NCBI sequence read archive database with the accession number PRJNA632458.
In total, 48,339 unigenes were generated and annotated by searching the sequences against the NR, KEGG, GO and TFs databases, which produced 17,737 (36.69%), 13,336 (27.59%), 6,128 (12.68%) and 1,756 (3.63%) hits respectively (Table 1). There were only 1,752 (3.62%) were annotated in all of the databases, and 17,847 (36.92%) were annotated in at least one database. Among these, GO classi cation analysis showed that 48,339 unigenes were categorized into 21 biological processes, 16 cellular components and 9 molecular functions. The KEGG classi cation analysis showed that total unigenes were assigned to 6 KEGG categories in the top level: 12 metabolisms, 11 human diseases, 10 organismal systems, 4 cellular processes, 4 genetic information processing and 3 environmental information processing.

DEGs analysis and immune related pathway identi cation
DEGs analysis among the four groups Pre_F, Post_F, Pre_M and Post _M was applied to identify gene changes pre/post-spawning. All DEGs with the absolute value of fold change > 2, adjust P-value < 0.01 are showed in Fig.1. A total of 1,805 DEGs were detected between the Pre_F vs. Post_F, of which 722 and 715 unigenes were up-regulated and down-regulated respectively. For the Pre_M vs. Post_M comparison, 2,488 DEGs were detected, of which 772 and 1716 unigenes were up-regulated and down-regulated respectively. As shown in Fig. 1, the number of up-regulated genes were basically equal between male groups (Pre_M vs. Post_M) and female groups (Pre_F vs. Post_F), but the number of down-regulated genes in male was more than double in female after spawning.
To further explore the function of the identi ed DEGs, GO classi cation analyses were conducted and 7 categories attracted our attention, which were cellular process, metabolic process, biological regulation, cell, organelle, binding and catalytic activity. In these categories, while the number of DEGs in male was at least 40% more than that of DEGs in female, the proportion of down-regulated genes in male was at least 50% more than that of down-regulated genes in female. These categories were marked in Figure 2, and a list of DEGs numbers was provided in Table 2. Additionally, in the majority of KEGG classi cation categories, the number of DEGs in both genders were approximately in the same, except "human diseases" that the number of DEGs in male was 80% 256/141 more than that of DEGs in female, and the proportion of down-regulated genes in male was 171% (192/39) more than that of down-regulated genes in female. The function of the more down-regulated DEGs in male is to participate in metabolism, catalytic activity, cell function and disease, which implied male reduced more these competences after spawning.
GO term and KEGG pathway enrichment analyses were performed to explore the DEGs related to immune regulation. We screened differential genes on GO term by keywords "immune" in Pre_F vs. Post_F and Pre_M vs. Post_M comparisons, and enriched 11 and 4 differential genes, respectively. Here, both quantity of DEGs and the type of GO enrichment were more in females than in males (Additional le 3, Fig.S1). To better understand the immune-related pathways changes pre/post-spawning, we evaluated signi cant enriched pathways about "Immune system" of secondary level of KEGG pathway enrichment, and enriched 19 and 20 differential genes in female and male respectively (Fig. 3). Interestingly, common pathways with high value enriched in different genders included IL-17 signaling pathway, NOD-like receptor signaling pathway, TLR signal pathway, RIG-I-like receptor signaling pathway, which played an important role in innate immunity. The results here implied that immune-related cytokines and in ammation were indeed affected by spawning.

Differences expression of immune-related genes between reproductive states
Based on the DEGs enriched by KEGG pathway in the Section 3.2 (19 and 20 genes in female and male respectively), 22 genes were further screened by Nr annotation. Based on their function in pathway, we classi ed them into seven categories: complement and coagulation cascades, apoptosis, TLR signal pathway, NFκB signal pathway, leukocyte transendothelial migration, MAPK signal pathway and others (Fig. 4). The expressions of genes in complement and coagulation cascades are almost 2/3 signi cantly down-regulated after spawning in males, whereas the analysis showed the expressions of genes in apoptosis, TLR signal pathway, NFκB signal pathway and HSP family are signi cantly up-regulated after spawning in females. It means that male decreases immune-competent of complement, while female improves immune function in apoptosis, antimicrobial peptides and heat shock.
To validate the reliability of DEGs identi ed by RNA-Seq, 10 out of these 22 genes were selected for further relative expression quanti ed by qRT-PCR. These genes belong to NFκB signal pathway (IKK, DTHD1), TLR signal pathway (TLR4), apoptosis (BIRC7, IAP and CREM), MAPK signal pathway (ARAFlike, MAPKK3), complement and coagulation cascades (CTRP2) and others immune-related genes (PSME3). The expressions of these genes pre/post-spawning showed in Fig.5. Compared their expression in RNA-seq (Fig. 4), 8/10 genes had a common expression pattern, only BIRC7 and IAP are not consistent in male.

Gene expression changes pre/post-spawning under Vibrio challenge
The transcriptome analysis revealed that immune-related genes changes pre/post-spawning in clam M. petechialis in natural condition. Here, we want to know how these genes expression pre/post-pawning under the pathogen challenge. Four immune-related DEGs (ARAF-like, CREM, TLR4 and TBK1) were selected for further examination. As shown in Fig. 6, without Vibrio challenge (0h), the expression patterns of these DEGs were consistent with the changes revealed by RNA-seq under natural condition ( Fig. 4), that was, after spawning their expressions were signi cantly increased in females but there was no signi cant difference in males. Under Vibrio challenge, the four genes also showed up-regulated expression in females at 12h or 24h; in male, the non-differentially expressed genes showed a signi cant up-regulated expression after spawning at 24h. These results indicated that the response of these immune-related genes became more sensitive and their expression changes more drastically after spawning under the pathogen challenge. It means pathogen reinforced the expression differences of immune-related genes between pre-spawning groups and post-spawning groups.

Assessment of immune status of clam under Vibrio challenge
To further verify the effects of spawning, we chose two previously reported genes, DOUX and MITF both involving the immune defense of M. petechialis [24,37], to detect the immune changes of the clams pre/post-spawning. As shown in Fig. 7, under bacterial immersion challenge, the mRNA expression of MpDOUX and MpMITF were signi cantly up-regulated after spawning in both genders. We also tested the protein changes of MpDOUX and MpMITF in female clams by western blot because we perceived that many immune genes had more signi cant differences after spawning in female based on transcriptome data (Fig. 4). As shown on Fig. 8 and Additional le 4, the protein expressions increased at 24 h for MpDOUX and at 12 h for MpMITF in the post-spawning individuals under Vibrio challenge. The upregulation at the transcription and protein levels of these maker genes further con rmed our conclusion that pathogen reinforced the expression differences of immune-related genes between pre-spawning and post-spawning groups.

Discussion
The existence of temporal changes in immunocompetence associated with reproduction is widely accepted [38]. The purpose of present study was to examine the effect of spawning on molluscan immunity. Through transcriptome analysis, we compared and analyzed the globle changes of genes in clams pre/post-spawning. According to the annotation and enrichment analysis, we obtained immunerelated pathways and key genes and they showed signi cant changes after spawning. Further Vibrio challenge experiment clari ed that immune-related genes expression changed after spawning during pathogen invasion.
Previous studies demonstrated the trade-off between reproduction and immunity can in uenced by sex [39,40], and here the transcriptome analysis showed gene expression differences between male and female in the clam M. petechialis. In detail, the number of down-regulated genes in male is more than double in female after spawning. Through GO and KEGG enrichment, we found these decreased genes were distributed in metabolism, catalytic activity, disease and so on. The effect of reproduction on metabolism caused researchers attention very early, Williams [5] reported high production e ciency resulted in some undesirable effects that mostly related to metabolic imbalance. In our results, the decrease in metabolic capacity occurred primarily in males, possibly due to the differences in hormone levels between males and females, as Zera and Zhao [41] suggested that hormonal tight controlled of metabolism.
Immunological performance may be subject to rapid temporal changes due to possible resource reallocation between the immune system and reproduction, or through immunomodulation by reproductive hormones [42]. In our study, major immune-related pathways and genes were explored by enrichment analysis. The pathways changed pre/post-spawning in both genders were IL-17 signaling pathway and three classic PRRs-related pathways (NOD-like receptor, TLR and RIG-I-like receptor pathways), which involved in cytokine secretion, in ammation and apoptosis [43][44][45][46]. By further analysis of up/downregulated genes we found the up-regulated DEGs in apoptosis, TLR signal pathway and heat shock protein after spawning only occurred in females. However, DEGs in complement and coagulation cascades were down-regulated after spawning in males. Similarly, a positive relationship between antimicrobial capacity and reproductive investment in female was found in insect during breeding season [10], moreover, the latest research in sea cucumber revealed that the complement activation were weakened after spawning [47]. All these results suggested that females improve their immunity (like apoptosis, antimicrobial peptides and heat shock), whereas males weaken part of their immunity (such as the complement system) after spawning. One explanation in vertebrates for the down-regulation of male immune defense is that testosterone weakens immunity [48], but it was unclear what mechanisms in invertebrate.
Mollusc usually reproduce in summer and they are threatened by more enriched pathogens in the water simultaneously, resulting in mass mortality in breeding season. Therefore, we further focused on the response of immune pathways and genes to pathogen stress before and after spawning. In this study, Vibrio was selected as pathogen to challenge this clam, and four immune DEGs from different pathways were used to test their expression differences between pre/post-spawning groups under Vibrio challenge. Compared to pre-spawning individuals, the four genes were all up-regulated expression in the postspawning clams under Vibrio challenge. This result indicated pathogen could reinforce the differences between different reproductive state. The expression changes of these genes that involved in apoptosis, TLR signaling pathway and MAPK signaling pathway [49][50][51] can re ect the immune status, however, in this clam these genes have not been well characterized before, and we need more marker genes to verify our nding in M. petechialis. In our previous works, several immune genes have been identi ed from M. petechialis and their expression could indicate the immune response of this clam to Vibrio infection [52,53]. Among them, microphthalmia-associated transcription factor (MITF) is a basic helix-loop-helixleucine zipper protein that plays a key role in cell proliferation and immune defense, and Dual oxidase (DUOX) is a main ROS source and plays role in host resistance against infection by diverse pathogens [37,54]. In present study, we used MpMITF and MpDOUX as makers to evaluate the immune response changes of this clam. Their up-regulation expression at the transcription and protein levels after spawning under the Vibrio challenge further veri ed that spawning can strengthen immune response to pathogen invasion.
Although part of our results differ from the theory of trade-off, they coincide with the ndings in sh that the changes between pre-spawning and post-spawning showed different trends under differential measure parameters of the immune system [38]. These results further support the idea that the trade-off could be affected by factors such as sex, individual condition and parasite infection [10]. Anyhow, we investigated the trade-off between reproduction and immune in natural condition and pathogenic stress environment, which can help us to understand the interaction between these two complex physiological processes. The combination of the current results and previous research may give us some clues to explain the possible reason for clam summer mortality. As an important nding of this study, the downregulate genes related to the complement system in males will reduce their ability to withstand adversity stress. For females, although immune function are enhanced after spawning, the large investment of immunization and reproduction consume huge energy and direct somatic damage by the immune response can also be responsible for the impact on survival [55].

Declarations Ethics approval and consent to participate
This study does not involve endangered invertebrates. According to the national regulation (Fisheries Law of the People's Republic of China), no permission is required to collect the animals and no formal ethics approval is required for this study.

Consent for publication
Not applicable.

Availability of data and material
The datasets analysed during the current study are available from the corresponding author on reasonable request. BL designed the project, interpreted the data and revised the manuscript. DW performed the experiment and contributed to data interpretation rafted the manuscript. All authors read and approved the nal manuscript.   DEGs associated with immune function in the hepatopancreas of M. petechialis. The color of each cell represents log2fold change in the genes expression level of each comparison. Red is up-regulated and blue is down-regulated after spawning compared to pre-spawning, but only genes with a |log2FoldChange| > 2 and a p-value threshold after FDR correction of 0.05 were considered signi cantly differentially-expressed genes marked by ▲.

Figure 6
Relative mRNA expression of DEGs in the hepatopancreases of M. petechialis at 0-24 hour after immersion in V. parahaemolyticus by qRT-PCR. Error bars represent the SD. The asterisk (*) in postspawning clam represents signi cant differences found when compared to pre-spawning one at each time point (P < 0.05). "F" and "M" means "female" and "male" respectinely. Figure 7