Skip to main content
Download PDF

Genome-wide DNA methylation reveals potential epigenetic mechanism of age-dependent viral susceptibility in grass carp

Immunity & Ageing volume 19, Article number: 28 (2022) Cite this article

Abstract

Background

Grass carp are an important farmed fish in China that are infected by many pathogens, especially grass carp reovirus (GCRV). Notably, grass carp showed age-dependent susceptibility to GCRV; that is, grass carp not older than one year were sensitive to GCRV, while those over three years old were resistant to this virus. However, the underlying mechanism remains unclear. Herein, whole genome-wide DNA methylation and gene expression variations between susceptible five-month-old (FMO) and resistant three-year-old (TYO) grass carp were investigated aiming to uncover potential epigenetic mechanisms.

Results

Colorimetric quantification revealed that the global methylation level in TYO fish was higher than that in FMO fish. Whole-genome bisulfite sequencing (WGBS) of the two groups revealed 6214 differentially methylated regions (DMRs) and 4052 differentially methylated genes (DMGs), with most DMRs and DMGs showing hypermethylation patterns in TYO fish. Correlation analysis revealed that DNA hypomethylation in promoter regions and DNA hypermethylation in gene body regions were associated with gene expression. Enrichment analysis revealed that promoter hypo-DMGs in TYO fish were significantly enriched in typical immune response pathways, whereas gene body hyper-DMGs in TYO fish were significantly enriched in terms related to RNA transcription, biosynthesis, and energy production. RNA-seq analysis of the corresponding samples indicated that most of the genes in the above terms were upregulated in TYO fish. Moreover, gene function analysis revealed that the two genes involved in energy metabolism displayed antiviral effects.

Conclusions

Collectively, these results revealed genome-wide variations in DNA methylation between grass carp of different ages. DNA methylation and gene expression variations in genes involved in immune response, biosynthesis, and energy production may contribute to age-dependent susceptibility to GCRV in grass carp. Our results provide important information for disease-resistant breeding programs for grass carp and may also benefit research on age-dependent diseases in humans.

Background

Grass carp (Ctenopharyngodon idellus) are an important farmed fish in China, accounting for more than 18% of the total freshwater aquaculture production in this country. The production of grass carp reached 5.53 million tons in 2019, indicating the great commercial value of this fish [1]. Nevertheless, grass carp reovirus (GCRV)-induced hemorrhagic disease is a major threat to the grass carp aquaculture industry [2]. Consequently, the disease susceptibility of grass carp to GCRV is of major interest to fish-breeding geneticists aiming to identify strategies for disease-resistant breeding [3,4,5,6,7]. Our previous study revealed that age is an important factor that influences the susceptibility of grass carp to GCRV infection. GCRV infection in grass carp not older than one year caused death and hemorrhagic symptoms, while no dead individuals were observed in grass carp over three years of age after virus infection [8]. Moreover, we found that the level of immune response and host biosynthesis and metabolism are involved in the age-dependent viral susceptibility of grass carp [9]. However, the reason for the different immune response levels or host biosynthesis and metabolism levels in grass carp of different ages remains poorly understood.

In mammals, age is also a crucial factor that affects disease outcomes [10]. Pseudorabies virus (PRV) infection causes more severe clinical diseases in newborns than in older individuals [11]. Infants and young children are more sensitive to reovirus-induced encephalitis than adults [12]. Mice display an age-dependent acceleration of mortality due to influenza virus (IAV) infection, and are useful for modeling human aging and the outcomes of IAV infection [13]. It has been suggested that the maturation of the innate immune system and the expression of type I interferon genes may contribute to age-related viral susceptibility [14]; however, the reason for the differential gene expression patterns between young and adult patients is still unclear.

DNA methylation is an important mediator of gene expression regulation, and many genes have been reported to be influenced by DNA methylation at different ages [15,16,17]. Transporter associated with antigen processing 1 (TAP1) gene promoter methylation level was decreased, while mRNA expression was increased in Sutai piglets from birth to weaning age (8, 18, 30, and 35 days old), which may contribute to the resistance of piglets to Escherichia coli F18 at 35 days of age [18]. Genome-wide DNA methylation analysis of CD4+ and CD8+ T cells from younger and older individuals has revealed an increased number of methylation changes and higher methylome variation in CD8+ T cells with age, implying a link between age-related epigenetic changes and impaired T cell function [19]. A genome-wide DNA methylation survey of skeletal muscle between young and old healthy humans showed a dynamic interrelationship between DNA methylation, gene expression, age, and exercise [20]. Therefore, DNA methylation appears to play an important role in age-related gene expression and is thought to underlie many age-related physiological phenomena, such as age-dependent susceptibility to viral infection.

In this study, we investigated genome-wide DNA methylation and gene expression variations between susceptible five-month-old (FMO) grass carp and resistant three-year-old (TYO) grass carp to reveal the potential epigenetic mechanism of age-dependent viral susceptibility in grass carp. Differentially methylated regions (DMRs) and differentially expressed genes (DEGs) were identified. The correlation between DNA methylation and gene expression was determined. Moreover, the roles of the two metabolism-related genes during viral infection were investigated. We believe our results will provide valuable information for disease-resistant breeding of grass carp and may also benefit research on age-dependent diseases in humans.

Results

Global methylation profile and whole-genome bisulfite sequencing

A viral challenge experiment was carried out for approximately 200 FMO and 200 TYO grass carps. For the FMO group, a total mortality of 85% was reached at 14 days, with the first death recorded as early as 6 days post-infection (Fig. 1A). However, none of the individuals in the TYO group died (Fig. 1A). This result was consistent with our previous report [9], further confirming age-dependent susceptibility to GCRV in grass carp. Moreover, we analyzed the global DNA methylation status of spleens from the two groups by colorimetric quantification of 5-methylcytosine. As shown in Fig. 1B, the global methylation level in TYO fish was higher than that in FMO fish, although this difference was not significant.

Fig. 1
figure 1

Age-dependent susceptibility to grass carp reovirus (GCRV) and global methylation profile of the two fish groups. A The cumulative mortality of the two grass carp groups after GCRV infection. B Global DNA methylation level of two grass carp groups. Spleen samples from the two groups were obtained, and the global DNA methylation level (5-mc) was analyzed. Data are shown as the mean ± standard deviation of three replicates. C The percentage of methylcytosines in CG, CHG, and CHH contexts. D Logo plots of the sequences proximal to sites of mCG, mCHG, and mCHH contexts. E Distribution of DNA methylation levels in the two groups. The y axis indicates the fraction of all mCs that display each methylation level (x axis), where the methylation level is the mC/C ratio at each cytosine

Full size image

To further reveal the potential epigenetic mechanisms underlying age-dependent viral susceptibility in grass carp, we performed whole-genome bisulfite sequencing (WGBS) on samples collected from FMO and TYO fish at base-pair resolution. Three replicates were processed for each group, yielding six libraries (n = 5 for each library). All libraries were sequenced on an Illumina HiSeq X Ten platform to generate 150 bp paired-end reads. Each library yielded > 33.68 GB clean base and > 23.06 × mean coverage depth. The bisulfite conversion rates for all libraries exceeded 99.88%. Moreover, for all libraries, more than 82.78% of the genomic sites were covered by at least five unique reads, and more than 70.70% of the sites were covered by at least ten unique reads (Table 1). Pearson’s correlation analysis and principal component analysis of the six libraries showed that samples from each group clustered together (Additional file 1: Fig. S1A and S1B). Collectively, these results confirm the high quality and repeatability of the WGBS data and its suitability for further analysis. WGBS data were deposited in the Sequence Read Archive (SRA) at the National Center for Biotechnology Information (NCBI) (accession number: PRJNA638254).

Table 1 Summary of WGBS data
Full size table

In the detected methylation sites, the average methylated cytosine (mC) percentages of whole genomic cytosines were 4.99 and 4.80% for FMO and TYO fish, respectively. In the CG context, approximately half of the cytosines were methylated (51.42% in FMO fish and 49.40% in TYO fish), whereas the methylation rates of cytosines in the CHG and CHH contexts (where H is A, C, or T) were no more than 0.10 and 0.09%, respectively (Additional file 3: Table S1). Among the mC sites that were identified, more than 98% were in the CG context, whereas no more than 0.4 and 1.5% were in the CHG and CHH contexts, respectively (Fig. 1C). The genome sequence preference proximal to the sites of methylated CG (mCG), methylated CHG (mCHG), and methylated CHH (mCHH) contexts was also analyzed. There was no sequence preference in the mCG-flanking regions or upstream of the mCHG and mCHH contexts; however, the base following the mCHG and mCHH contexts was almost always adenine, followed by thymine, while cytosine was observed less often (Fig. 1D). The methylation level (ML), defined as the proportion of reads covering each mC relative to the total reads covering the sites, was calculated. The results showed that more than 80% of the mC sites displayed high ML (ML ≥ 70%) (Fig. 1E).

Gene methylation profile

To characterize the methylation profile of grass carp genes, the relative ML in the context of gene regions and their upstream and downstream regions was calculated. In general, the relative ML of mCGs in the gene body regions was higher than that in the 5′ upstream and 3′ downstream regions in both groups, whereas the ML in mCHG and mCHH was more complex (Fig. 2A and Additional file 1: Fig. S1C). Interestingly, a sharp decrease in ML was observed across the boundaries of gene body regions and upstream or downstream regions (Fig. 2A). We further divided the genome into different elements, including the promoter, 5′ untranslated region (UTR), exon, intron, 3′ UTR, and repeats. The results suggested that the promoter and 5′ UTR had relatively low ML in mCG, followed by the exon, whereas the repeats had the highest ML (Fig. 2B and Additional file 1: Fig. S1D).

Fig. 2
figure 2

Distribution of mCs in different genome regions and the differentially methylated regions (DMRs) identified between the two groups. A Relative methylation level in gene body regions and of their upstream and downstream regions. B Relative methylation level in different gene elements (promoter, 5′ UTR, exon, intron, 3′ UTR, and repeat regions). C The heatmap of DMRs. Red represents hypermethylation, whereas blue represents hypomethylation. D The distribution of DMRs in the genome. Black dots indicate hyper-DMRs, while gray dots represent hypo-DMRs

Full size image

DMRs between FMO and TYO fish

To further investigate the epigenetic differences between the two fish groups, the WGBS data of TYO fish were compared with those of FMO fish, and DMRs were identified using DSS software. As most mC sites occurred in CG contexts (Fig. 1B), only DMRs in CG contexts were considered for further study. A total of 6214 DMRs were identified between the two fish groups, including 5261 hyper-DMRs (hypermethylation in TYO fish compared with FMO fish) and 953 hypo-DMRs (hypomethylation in TYO fish compared with FMO fish) (Additional file 4: Table S2). The number of hyper-DMRs was greater than that of hypo-DMRs. Interestingly, we also found that the mean ML of DMRs in TYO fish was higher than that in FMO fish (Additional File 1: Fig. S1E). The heatmap of the DMRs also showed a similar result (Fig. 2C). The length distribution of the DMRs was analyzed, and the results showed that most of them were shorter than 400 bp (Additional file 1: Fig. S1F). Circos analysis in a circular layout revealed that the DMRs were distributed uniformly in the genome, except for supercontig C101000030 (Fig. 2D), in which more DMRs were identified. We identified 4052 DMGs, including 2855 DMGs that harbored DMRs in gene body regions (from transcription start site (TSS) to transcription end site (TES)) and 1197 DMGs containing DMRs in promoter regions (Additional file 5: Table S3).

Correlation between DNA methylation and gene expression

The gene expression profiles of the two groups were investigated by RNA-seq using the corresponding samples. Three replicates were performed for each group, and six libraries were obtained. All libraries were sequenced on an Illumina X Ten platform to generate 150 bp paired-end reads. Each library generated ≥7.5 GB of clean data and showed Q30 ≥ 94% (Table 2), implying sufficient quality of RNA-seq data. RNA-seq data were deposited in the SRA at NCBI (accession number: PRJNA634937). To analyze the correlation between DNA methylation and gene expression, the DNA methylation levels and gene expression levels of the gene body regions, 5′ upstream and 3′ downstream, in the two fish groups were analyzed. We first divided the genes into four groups according to their mRNA expression levels: none (unexpressed genes, FPKM < 1), low (FPKM < bottom 25% of expressed genes), medium (bottom 25% of expressed genes ≤ FPKM < top 25% expressed genes), and high (FPKM ≥ top 25% expressed genes). The DNA methylation levels in the different groups were compared. In general, the methylation levels in the four groups began to decrease from 2 kb upstream of the TSS of the genes but increased downstream of the TSS. In the 5′ upstream region, the genes with high mRNA expression levels showed the lowest DNA methylation levels, whereas the unexpressed genes displayed the highest methylation levels (Fig. 3A and Additional file 2: Fig. S2A), suggesting a negative correlation between promoter methylation and gene expression. Nevertheless, the correlation between DNA methylation of the gene body or 3′ downstream and mRNA expression is complex. Moreover, we also classified genes into five categories based on their methylation levels, from the bottom 20% to the top 20%, corresponding to the 1st to 5th groups, and the mRNA expression levels of these groups were investigated. The 1st group of genes in the promoter regions showed the highest expression levels (Fig. 3B and Additional File 2: Fig. S2B), which further implies that DNA methylation in promoter regions, especially hypomethylation, is negatively associated with gene expression. Interestingly, the 5th genes in the gene body region showed the highest expression levels (Fig. 3C and Additional file 2: Fig. S2C), indicating that DNA methylation in gene body regions, particularly hypermethylation, is positively correlated with gene expression.

Table 2 Summary of RNA-seq data
Full size table
Fig. 3
figure 3

DNA methylation and mRNA expression of genes in different categories. A The methylation levels of different gene groups (with different mRNA expression levels) in the five-month-old group. Genes were divided into four groups according to the mRNA expression levels: none, low, medium, and high, and the DNA methylation levels in different groups were compared. B mRNA expression levels of different gene categories (with different methylation levels) in promoter regions of the five-month-old group. Genes were classified into five categories based on the methylation levels, from the bottom 20% to the top 20%, corresponding to the 1st to 5th groups, and the mRNA expression levels were compared. C mRNA expression levels of different gene categories (with different methylation levels) in gene body regions of the five-month-old group. D Correlation between DNA methylation and gene expression of promoter hypo-DMGs. E Correlation between DNA methylation and gene expression of gene body hyper-DMGs

Full size image

Furthermore, the DMRs obtained by WGBS were selected, and the relationship between DNA methylation levels and mRNA expression levels was analyzed using Spearman’s rank correlation method [21]. For the DMRs in the promoter regions, both hyper-DMRs and hypo-DMRs displayed a slightly negative correlation between the DNA methylation levels and mRNA expression levels, especially promoter hypo-DMRs (Fig. 3D and Additional file 2: Fig. S2D). However, for the DMRs in the gene body regions, only the hyper-DMRs showed a slightly positive correlation (r = 0.22, p < 0.05) (Fig. 3E), whereas no correlation was observed in the hypo-DMRs (p > 0.05) (Additional file 2: Fig. S2E). These findings further suggest that hypomethylation in promoter regions and hypermethylation in gene body regions are associated with gene expression.

Functional enrichment analysis of DMGs

Due to hypomethylation in promoter regions, and hypermethylation in gene body regions is associated with gene expression. Therefore, we performed an enrichment analysis for promoter hypo-DMGs (hypomethylated DMGs in promoter regions) and gene body hyper-DMGs (hypermethylated DMGs in gene body regions). GO enrichment analysis revealed that promoter hypo-DMGs were not significantly enriched in any GO terms (corrected P value = 1), whereas the gene body hyper-DMGs were enriched in GO terms related to transcription, biosynthesis, and metabolic processes, such as regulation of transcription (DNA-templated), regulation of nucleic acid-templated transcription, regulation of RNA biosynthetic process, regulation of RNA metabolic process, and regulation of cellular macromolecule biosynthetic process. The top ten significantly enriched GO terms are listed in Table 3. In addition, KEGG enrichment analysis was performed. Interestingly, some typical immune-related pathways, such as herpes simplex infection, cytokine-cytokine receptor interaction, RIG-I-like receptor signaling pathway, and Toll-like receptor signaling pathway, were significantly enriched in promoter hypo-DMGs (Fig. 4A). The gene body hyper-DMGs were enriched in KEGG pathways, such as focal adhesion, ECM-receptor interaction, melanogenesis, dorso-ventral axis formation, adherens junction, and some pathways involved in metabolism and biosynthesis (propanoate metabolism and biosynthesis of amino acids) (Fig. 4B). Detailed information on the KEGG enrichment analysis is shown in Additional file 6: Table S4.

Table 3 GO enrichment analysis of promoter hypo-DMGs and gene body hyper-DMGs (top ten terms)
Full size table
Fig. 4
figure 4

Enrichment analyses of DMGs in different categories and mRNA expression profiles of genes in different categories. A KEGG enrichment analyses of DMGs in promoter hypo-DMGs. B KEGG enrichment analyses of DMGs in gene body hyper-DMGs. C Scatterplots of genes involved in the immune response. D Scatterplots of genes involved in protein synthesis. E Scatterplots of genes involved in energy metabolism. The red dots represent genes upregulated in TYO fish, while the blue dots represent genes downregulated in FMO fish

Full size image

Gene expression patterns of promoter hypo-DMGs and gene body hyper-DMGs

Hypomethylation in the promoter regions was thought to be associated with transcriptional activation (Fig. 3B), and KEGG enrichment analysis showed that promoter hypo-DMGs were enriched in classical innate immune pathways (Fig. 4A). It is well known that the immune response plays an important role in host defense against pathogen infection. Therefore, we analyzed the expression profiles of immune-related genes between two groups, and the results showed that most of them (74.87%) were upregulated in TYO fish (Fig. 4C). Moreover, it was proposed that gene body hypermethylation involved in gene expression (Fig. 3C) and gene body hyper-DMGs were enriched in pathways related to transcription, biosynthesis, and metabolism (Table 3). Protein synthesis and energy metabolism also benefit organisms that are resistant to viral infection. We therefore investigated the mRNA expression patterns of genes involved in protein biosynthesis and energy metabolism to further confirm this hypothesis. As expected, most of the genes involved in biosynthesis (88.31%) and energy metabolism (86.96%) were also upregulated in TYO fish (Fig. 4D, E).

Confirmation of WGBS data by BS-PCR and qPCR

These results indicated that DNA methylation in promoter regions was negatively associated with gene expression, and that promoter hypo-DMGs were enriched in immune-related pathways (Fig. 4A). Therefore, two immune-related genes, signal transducer and activator of transcription 1b (stat1b) and TNF receptor-associated protein 1 (trap1), were selected for further confirmation. The CpG loci of candidate DMRs in the gene promoter were amplified by bisulfite sequencing PCR (BS-PCR) and mRNA expression levels were determined by qPCR. As shown in Fig. 5A and B, BS-PCR showed that the DNA methylation level of the promoter regions of two immune-related genes in FMO fish was higher than that in TYO fish, while the mRNA expression levels presented opposite trends, implying that DNA methylation in promoter regions was negatively associated with gene expression. Moreover, two genes involved in protein biosynthesis and energy metabolism, glyceraldehyde-3-phosphate dehydrogenase glyceraldehyde-3-phosphate dehydrogenase, spermatogenic (gapdhs), and lactate dehydrogenase A (ldha), were selected, and the CpG loci of DMRs in the gene body region were amplified for BS-PCR verification and mRNA expression levels were determined by qPCR. As expected, both of the two genes showed higher DNA methylation levels in TYO fish than in FMO fish, and the mRNA expression also showed similar trends (Fig. 5C and D), further suggesting a positive correlation between gene body DNA methylation and gene expression.

Fig. 5
figure 5

Confirmation of WGBS data by BS-PCR and qPCR. A DNA methylation and gene expression profiles in the promoter region of stat1b gene. B DNA methylation and gene expression profiles in the promoter region of trap1 gene. C DNA methylation and gene expression profiles in the gene body region of gapdh gene. D DNA methylation and gene expression profiles in the gene body region of idha gene. Top panels: schematic diagram of promoter regions or gene body regions of four genes. The circles indicate the candidate CpG loci in the DMRs. Middle panels: Validation of CpG methylation by BS-PCR. Fifteen subclones from each group were selected for the BSP analysis. The solid circles represent the methylated CpG sites, and the open circles indicate the unmethylated CpG sites. Bottom panels: Comparison of gene expression levels between two groups by qPCR

Full size image

The antiviral effects of metabolism-related genes

The above results showed that most of the genes involved in the immune response and energy metabolism were upregulated in TYO fish when compared with FMO fish. It is well known that immune-related genes play important roles in host defense against pathogen infection; however, the role of metabolism-related genes during viral infection is unclear. Therefore, two metabolism-related genes, gapdhs and ldha, were selected for overexpression or knockdown to investigate their roles in GCRV infection. As shown in Fig. 6, after overexpression of two genes in grass carp ovary (GCO) cells (Fig. 6A), the relative copy number of GCRV nonstructural protein NS80 and structural protein VP5 was significantly lower than that in the negative control (Fig. 6B). In contrast, when the expression of the two genes was knockdown by specific siRNAs (Fig. 6C and E), the relative copy numbers of NS80 and VP5 were significantly higher than those in the negative control (Fig. 6D and F). Collectively, these results indicated the antiviral effects of two metabolism-related genes.

Fig. 6
figure 6

The antiviral effects of two metabolism-related genes. A Western blotting confirmation of the overexpression of two genes in transfected cells. Cells were transfected with the corresponding plasmids for 24 h and harvested for western blotting analysis. Rabbit anti-HA antibody and goat anti-rabbit IgG were used as the primary and secondary antibodies, respectively. B The relative copy number of genes encoding GCRV nonstructural protein NS80 and structural protein VP5 in cells transfected with different plasmids. C RT-qPCR confirm the knockdown of gapdhs by siRNAs. Cells were transfected with the corresponding siRNA for 24 h and harvested for RT-qPCR analysis. D The relative copy number of genes encoding NS80 and VP5 in cells transfected with different siRNAs. E RT-qPCR confirm the knockdown of ldha by siRNAs. F The relative copy number of genes encoding NS80 and VP5 in cells transfected with different siRNAs. Significant differences (P < 0.05) between the control and treated groups are indicated with asterisks (*)

Full size image

Discussion

Age-dependent viral susceptibility has been reported in several fish species. Nervous necrosis virus (NNV) infection in barramundi (Lates calcarifer) causes nervous necrosis at 3 and 4 weeks of age, but develops subclinical symptoms when fish age is more than 5 weeks [22]. Spring viremia of carp virus (SVCV) challenge in North American fish species revealed that younger fish were more vulnerable to SVCV infection than older fish [23]. In our previous study, we found that grass carp showed age-dependent susceptibility to reovirus, and the immune response level or biosynthesis and metabolism levels may account for this phenomenon [9]. Understanding the underlying mechanisms of different immune responses or metabolism patterns in fish of different ages will benefit disease-resistant breeding programs.

DNA methylation, the most widely understood type of epigenetic modification, has been reported to play a crucial role in transcriptional regulation [24, 25]. It is generally accepted that DNA hypomethylation in gene promoters is linked to transcriptional activation, whereas hypermethylation in these regions is correlated with transcriptional repression [26, 27]. Nevertheless, some evidence also shows that promoter hypermethylation appears to be associated with high transcriptional activity [28], indicating that DNA methylation contributes to transcriptional regulation in a more complex and dynamic manner. DNA methylation has been reported to be correlated with the disease resistance traits of Chinese tongue sole and Nile tilapia [29, 30]. Therefore, we hypothesized that DNA methylation may be involved in age-dependent viral susceptibility of grass carp. Here, we report a genome-wide survey of the DNA methylation status of spleens from susceptible FMO grass carp versus resistant TYO grass carp to reveal the potential epigenetic mechanisms of this phenomenon.

Interestingly, we found that the general DNA methylation pattern of grass carp was consistent with that of other species [20, 31,32,33]. For example, most mC sites occurred in the CG context, with relatively low ML in the promoter or 5′ UTR regions and high ML in repeat regions, and a negative correlation between methylation of promoter or 5′ UTR regions and gene expression. These results indicate a conserved role for DNA methylation during evolution. Interestingly, we showed that genes with the highest DNA methylation levels in the gene body regions had the highest mRNA expression levels in both groups (Fig. 3C and Additional file 2: Fig. S2C), suggesting that hypermethylation in gene body regions was positively correlated with mRNA expression. The positive correlation was further confirmed in two biosynthesis- or metabolism-related genes using BS-PCR and qPCR (Fig. 5). Several literatures also revealed the positive correlation between gene body methylation and gene expression [34, 35]. Moreover, it was reported that gene body hypermethylation is frequent occurred in constitutively expressed genes, whereas least frequent in variable expressed genes [36, 37]. Enrichment analysis revealed that gene body hyper-DMGs were involved in genes related to transcription, biosynthesis, and metabolism, all of them are constitutively expressed genes. It was proposed that the function of gene body hypermethylation is most likely homeostatic, which may maintain gene expression in response to developmental and/or environmental stresses [38]. Therefore, the high gene body methylation levels in TYO fish may imply that the TYO fish was more stable than FMO fish after GCRV infection, which may be one of the reasons for age-dependent viral susceptibility in grass carp.

Correlation analysis revealed DNA hypomethylation in promoter regions and hypermethylation in gene body regions associated with gene expression. Therefore, we performed enrichment analysis for promoter hypo-DMGs and gene body hyper-DMGs. The results suggested that promoter hypo-DMGs were significantly enriched in classical immune-related pathways (Fig. 4A), whereas the gene body hyper-DMGs were enriched in terms related to transcription, biosynthesis, and metabolism (Table 3 and Fig. 4B). Coincidently, RNA-seq using the corresponding samples also revealed that most of the genes participating in the immune response, biosynthesis, and metabolism were upregulated in resistant TYO fish (Fig. 4C-E). The role of immune-related genes in viral infection is well documented [39], whereas the functions of metabolism-related genes is still unclear. Gene function analysis by gene overexpression or knockdown indicated the antiviral effects of metabolism-related genes. The upregulation of metabolism-related genes indicated high metabolism levels in TYO fish, which was beneficial for host defense against pathogen infection. Previous studies have also suggested that elevated host metabolism promotes resistance to pathogen infection in fish [9, 40, 41]. Therefore, based on these results, we conclude that different immune response levels or different biosynthesis and metabolism levels between FMO and TYO fish may be mediated by DNA methylation.

Conclusion

In conclusion, genome-wide variations in DNA methylation and gene expression in grass carp of different ages were investigated. We identified 6214 DMRs and 4052 DMGs, with most DMRs and DMGs showing hypermethylation patterns in TYO fish. Correlation analysis indicated that DNA hypomethylation in promoter regions and hypermethylation in gene body regions are involved in transcriptional activation. DNA methylation and gene expression variations in genes involved in immune response, biosynthesis, and energy metabolism may contribute to age-dependent viral susceptibility in grass carp.

Methods

Experimental fish, GCRV challenge experiment, and sample collection

A total of 500 fish, containing 250 full-sib FMO and 250 full-sib TYO grass carp, representing fish sensitive and resistant to GCRV infection, respectively, were used in the study. All fish were bred at the Guan Qiao Experimental Station, Institute of Hydrobiology, Chinese Academy of Sciences (CAS), and acclimatized in aerated fresh water at 26–28 °C for one week before processing. To ensure that the grass carp were not previously exposed to GCRV, the fish were randomly selected for PCR detection of GCRV using GCRV-specific primers (Additional file 7: Table S5). The GCRV-free fish were fed commercial feed twice daily, and the water was exchanged daily.

After no abnormal symptoms were observed in the two groups, 200 fish from each group were subjected to a viral challenge experiment. Fish were infected with GCRV (GCRV subtype II, 2.97 × 103 RNA copies/μL) at a dose of 20 μL per gram of body weight by intraperitoneal injection. The injected fish were carefully monitored, and the number of dead fish in each group was counted daily. The experiment was concluded, and the total mortality was calculated when no mortality was recorded for seven consecutive days.

The remaining fish, which were not included in the viral challenge experiment, were used for sample collection. The sex ratio (male: female) of the sampled fish was approximately 1:1. In each group, 15 fish containing three biological replicates (n = 5 for each biological duplicate) were collected. Spleen samples were removed rapidly for analysis after the fish were anesthetized and euthanized using MS-222. The obtained samples were stored at − 80 °C until DNA or RNA extraction.

Whole-genome bisulfite sequencing

Genomic DNA was extracted from the spleen using the DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany). DNA purity and concentration were measured using a NanoPhotometer® spectrophotometer (IMPLEN, USA) and a Qubit® 2.0 Flurometer (Life Technologies, USA). DNA of sufficiently high quality was used for library construction. A total of 5.2 μg DNA spiked with 26 ng λDNA was fragmented by sonication to 200–300 bp with Covaris S220, followed by end repair and adenylation. Cytosine-methylated barcodes were ligated to sonicated DNA, according to the manufacturer’s instructions. These DNA fragments were then treated twice with bisulfite using an EZ DNA Methylation-Gold™ Kit (Zymo Research, USA), and the resulting single-stranded DNA fragments were PCR-amplified using KAPA HiFi HotStart Uracil + ReadyMix. The library concentration was quantified using a Qubit® 2.0 Fluorometer and quantitative PCR, and the insert size was assayed using an Agilent Bioanalyzer 2100 system. Libraries were sequenced on an Illumina HiSeq X Ten platform, and 150 bp paired-end reads were generated.

Data analysis

Raw data reads in fastq format were initially processed through a series of quality control (QC) procedures using in-house Perl scripts to obtain clean data. QC standards were as follows: (1) reads with ≥10% unidentified nucleotides (N) were removed; (2) reads with > 50% bases with phred quality < 5 were removed; (3) reads with > 10 nt aligned to the adapter were removed, allowing ≤10% mismatches; and (4) putative PCR duplicates generated by PCR amplification in the library construction process (reads 1 and 2 of two paired-end reads that were completely identical) were removed. The Q20, Q30, and GC content of the clean data were calculated, and all downstream analyses were performed using high-quality clean data. The reference genome of grass carp and the clean reads were transformed into bisulfite-converted sequences (C-to-T and G-to-A converted). Bismark software (version 0.16.3) was used to align the bisulfite-treated reads to the reference genome [42]. Sequence reads that produced the unique best alignment from the two alignment processes (original top and bottom strands) were then compared to the normal genomic sequence, and the methylation state of all cytosine positions in the read was inferred. The same reads aligned to the same regions of the genome in the two alignment processes were regarded as duplicate reads. The sequencing depth and coverage were summarized using duplicate reads. The methylation extractor results were transformed into bigWig format for visualization using the IGV browser. The sodium bisulfite non-conversion rate was calculated as the percentage of cytosines sequenced at cytosine reference positions in the lambda genome.

To calculate the methylation level of the sequence, we divided it into multiple bins with a size of 10 kb. The sum of methylated and unmethylated read counts in each window/bin was calculated. The methylation level (ML) for each window or C site shows the fraction of methylated Cs (mC) and is defined as ML = reads (mC)/reads (mC + umC), where umC is unmethylated Cs. The calculated ML was further corrected using the bisulfite nonconversion rate according to previous studies [43].

Differentially methylated regions analysis

Differentially methylated regions (DMRs) analysis of the two groups/conditions was performed using DSS software [44]. The core of DSS is a new dispersion shrinkage method for estimating the dispersion parameter from the gamma-Poisson or beta-binomial distributions. A DSS possesses three characteristics for detecting DMRs. First, spatial correlation and proper utilization of information from neighboring cytosine sites can help improve the estimation of methylation levels at each cytosine site and hence improve DMRs detection. Second, the read depth of cytosine sites provides information on precision that can be exploited to improve statistical tests for DMRs detection. Finally, the variance among biological replicates provided the information necessary for a valid statistical test to detect DMRs. Differentially methylated genes (DMGs) were identified as genes whose gene body regions (from TSS to TES) or promoter regions (upstream 2 kb from the TSS) overlapped with the DMRs. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of DMGs were performed using the GOseq R package and KOBAS software [45, 46].

RNA -seq

RNA from the spleens was isolated using TRIzol reagent (Invitrogen, USA) according to the manufacturer’s protocol. RNA of sufficiently high quality was used for library construction. Sequencing libraries were generated using the NEBNext Ultra RNA Library Prep Kit for Illumina (New England Biolabs, USA), following the manufacturer’s protocol. Libraries were sequenced on an Illumina X Ten platform, and 150 bp paired-end reads were generated. The output raw data were initially processed with in-house Perl scripts to obtain clean data by removing adapter, poly-N, and poor-quality data, as described above. The clean reads were mapped to the reference genome of grass carp using Hisat2 software [47], and gene expression levels were calculated using the FPKM (expected number of fragments per kilobase of transcript sequence per million base pairs sequenced) method [48]. Differential expression analysis of the two groups/conditions was performed using the DESeq package [49]. Genes with an adjusted P value < 0.05 (q value < 0.05) in DESeq analysis were assigned as DEGs.

Global DNA methylation measurement

Genomic DNA from the spleens of FMO and TYO fish was extracted using the traditional phenol–chloroform protocol with RNase treatment [50]. The global DNA methylation level was measured using the MethylFlash™ Global DNA Methylation (5-methylcytosine, 5-mC) ELISA Easy Kit (Epigentek, USA), according to the manufacturer’s protocol. The amount of input DNA for each assay was 100 ng, to ensure optimal quantification. Three replicates were performed for each group. Data are presented as the mean ± standard deviation of three replicates.

Bisulfite sequencing PCR

Bisulfite sequencing PCR (BS-PCR) was performed to confirm the DNA methylation levels of the DMGs between the two groups. DNA samples were treated with bisulfite using the EZ DNA Methylation Gold™ Kit (Zymo Research, USA), according to the manufacturer’s instructions. DMGs related to immune response or biosynthesis were selected as candidates for further verification. Genomic sequences of the candidate DMGs were extracted from the grass carp draft genome [3] and submitted to the online software MethPrimer [51] to obtain the CpG island region and CpG loci. Specific methylation primers (Additional file 7: Table S5) for BS-PCR were designed using the MethPrimer software to amplify the CpG island region or CpG loci that showed different methylation levels. The lengths of the amplified regions ranged from 150 to 350 bp. Moreover, an outer forward or reverse primer was designed for each amplified region to improve the specificity of the amplification. The obtained DNA fragments were cloned into the pMD18-T vector (TaKaRa, Japan). Fifteen clones from each group were randomly selected for sequencing to evaluate methylation status.

RT-qPCR

RT-qPCR was used to investigate the mRNA expression levels of DMGs in the two groups. Total RNA was isolated using TRIzol reagent (Thermo Fisher, USA) and first-strand cDNA was obtained using a ReverTra Ace kit (Toyobo, Japan). RT-qPCR was performed using a fluorescence quantitative PCR instrument (Bio-Rad, USA). Each reaction mixture contained 0.8 μL forward and reverse primers (for each primer), 1 μL cDNA template, 10 μL 2× SYBR green master mix (Toyobo, Japan), and 7.4 μL ddH2O. DMGs related to immune response or biosynthesis were selected as candidates for RT-qPCR verification. The cDNA sequences of candidate DMGs were extracted from the grass carp genome, and specific primers were designed using the Primer Premier 5software (Additional file 7: Table S5). Three replicates were included for each sample, and β-actin was used as an internal control to normalize the gene expression. The program was as follows: 95 °C for 10 s; 40 cycles of 95 °C for 15 s, 56 °C for 30 s, and 72 °C for 30 s; and melt curve construction. Relative expression levels were calculated using the 2-△△Ct method [52]. The data represent the mean ± standard deviation of three replicates.

Gene function analysis

Two metabolism-related genes, gapdhs and ldha, were selected for further study to investigate their roles during GCRV infection. The complete ORF sequences of the two genes were amplified, fused with an HA tag, and then inserted into the pcDNA3.1 vector (Additional file 7: Table S5). The resulting plasmids were transfected into the GCO cells. The empty vector pcDNA3.1 was also transfected at the same time as a negative control. After transfection for 24 h, cells were infected with GCRV at a multiplicity of infection (MOI) of 1. Cells were collected at 24 and 48 h post-infection, and the relative copy number of genes encoding the GCRV nonstructural protein NS80 or structural protein VP5 was determined by RT-qPCR. Moreover, siRNAs specifically targeting these two genes were synthesized to generate the corresponding siRNAs (Additional file 7: Table S5). The siRNA sequences were randomly scrambled and used as negative controls. All siRNAs used in this study were synthesized by GenePharma (Shanghai, China). GCO cells were transfected with the siRNAs and then infected with GCRV at a MOI of 1 24 h after transfection. Cells were harvested 24 and 48 h post-infection, and the relative copy number of genes encoding NS80 or VP5 was determined by RT-qPCR.

Statistical analysis

Statistical significance between different groups was determined by one-way analysis of variance (ANOVA) and Fisher’s least significant difference (LSD) post-test. Differences were considered statistically significant at P < 0.05. P < 0.05 was denoted by *.

Availability of data and materials

The WGBS and RNA-seq data generated in this study were deposited in the Sequence Read Archive (SRA) at the National Center for Biotechnology Information (NCBI) under accession numbers PRJNA638254 and PRJNA634937. Other data generated or analyzed during this study are included in this published article.

Abbreviations

CAS:

Chinese Academy of Sciences

DEGs:

Differentially expressed genes

DMGs:

Differentially methylated genes

DMRs:

Differentially methylated regions

FMO:

Five months old

GCRV:

Grass carp reovirus

IAV:

Influenza virus

mCs:

Methylated cytosines

ML:

Methylation level

NCBI:

National Center for Biotechnology Information

NNV:

Nervous necrosis virus

PRV:

Pseudorabies virus

SRA:

Sequence Read Archive

SVCV:

Spring viremia of carp virus

TES:

Transcription end site

TSS:

Transcription start site

TYO:

Three year old

UTR:

Untranslated regions

WGBS:

Whole-genome bisulfite sequencing

References

  1. Fisheries Bureau of Ministry of Agriculture in China. China fishery statistical yearbook of 2019. Beijing: China Agriculture Press; 2020. p. 21–37.

    Google Scholar 

  2. Rao Y, Su J. Insights into the antiviral immunity against grass carp (Ctenopharyngodon idella) reovirus (GCRV) in grass carp. J Immunol Res. 2015;2015:670437.

    Article  PubMed  PubMed Central  Google Scholar 

  3. Wang Y, Lu Y, Zhang Y, Ning Z, Li Y, Zhao Q, et al. The draft genome of the grass carp (Ctenopharyngodon idellus) provides insights into its evolution and vegetarian adaptation. Nat Genet. 2015;47:625–31.

    Article  CAS  PubMed  Google Scholar 

  4. He L, Zhang A, Pei Y, Chu P, Li Y, Huang R, et al. Differences in responses of grass carp to different types of grass carp reovirus (GCRV) and the mechanism of hemorrhage revealed by transcriptome sequencing. BMC Genomics. 2017;18(1):452.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  5. Ji J, Rao Y, Wan Q, Liao Z, Su J. Teleost-Specific TLR19 Localizes to Endosome, Recognizes dsRNA, Recruits TRIF, Triggers both IFN and NF-κB Pathways, and Protects Cells from Grass Carp Reovirus Infection. J Immunol. 2018;200(2):573–85.

    Article  CAS  PubMed  Google Scholar 

  6. Xue Y, Jiang X, Gao J, Li X, Xu J, Wang J, et al. Functional characterisation of interleukin 34 in grass carp Ctenopharyngodon idella. Fish Shellfish Immunol. 2019;92:91–100.

    Article  CAS  PubMed  Google Scholar 

  7. Shen Y, Wang L, Fu J, Xu X, Yue GH, Li J. Population structure, demographic history and local adaptation of the grass carp. BMC Genomics. 2019;20(1):467.

    Article  PubMed  PubMed Central  Google Scholar 

  8. Chen G, Xiong L, Wang Y, He L, Huang R, Liao L, et al. Different responses in one-year-old and three-year-old grass carp reveal the mechanism of age restriction of GCRV infection. Fish Shellfish Immunol. 2019;86:702–12.

    Article  CAS  PubMed  Google Scholar 

  9. He L, Zhu D, Liang X, Li Y, Liao L, Yan C, et al. Multi-omics sequencing provides insights into age-dependent susceptibility of grass carp (Ctenopharyngodon idellus) to Reovirus. Front Immunol. 2021;12:694965.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Chan YH, Ng LFP. Age has a role in driving host immunopathological response to alphavirus infection. Immunology. 2017;152(4):545–55.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Verpoest S, Cay B, Favoreel H, De Regge N. Age-dependent differences in pseudorabies virus Neuropathogenesis and associated cytokine expression. J Virol. 2017;91(2):e02058–16.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Wu AG, Pruijssers AJ, Brown JJ, Stencel-Baerenwald JE, Sutherland DM, Iskarpatyoti JA, et al. Age-dependent susceptibility to reovirus encephalitis in mice is influenced by maturation of the type-I interferon response. Pediatr Res. 2018;83(5):1057–66.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Smith CA, Kulkarni U, Chen J, Goldstein DR. Influenza virus inoculum volume is critical to elucidate age-dependent mortality in mice. Aging Cell. 2019;18(2):e12893.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  14. Giraldo D, Wilcox DR, Longnecker R. The type I interferon response and age-dependent susceptibility to herpes simplex virus infection. DNA Cell Biol. 2017;36(5):329–34.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Anastasiadi D, Piferrer F. A clockwork fish: age prediction using DNA methylation-based biomarkers in the European seabass. Mol Ecol Resour. 2020;20:387–97.

    Article  CAS  PubMed  Google Scholar 

  16. Mayne B, Korbie D, Kenchington L, Ezzy B, Berry O, Jarman S. A DNA methylation age predictor for zebrafish. Aging (Albany NY). 2020;12:24817–35.

    Article  CAS  Google Scholar 

  17. Unnikrishnan A, Freeman WM, Jackson J, Wren JD, Porter H, Richardson A. The role of DNA methylation in epigenetics of aging. Pharmacol Ther. 2019;195:172–85.

    Article  CAS  PubMed  Google Scholar 

  18. Dong W, Yin X, Sun L, Wang J, Sun S, Zhu G, et al. Age-associated methylation change of TAP1 promoter in piglet. Gene. 2015;573(1):70–4.

    Article  CAS  PubMed  Google Scholar 

  19. Tserel L, Kolde R, Limbach M, Tretyakov K, Kasela S, Kisand K, et al. Age-related profiling of DNA methylation in CD8+ T cells reveals changes inimmune response and transcriptional regulator genes. Sci Rep. 2015;5:13107.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Zykovich A, Hubbard A, Flynn JM, Tarnopolsky M, Fraga MF, Kerksick C, et al. Genome-wide DNA methylation changes with age in disease-free human skeletal muscle. Aging Cell. 2014;13(2):360–6.

    Article  CAS  PubMed  Google Scholar 

  21. Kretzmer H, Bernhart SH, Wang W, Haake A, Weniger MA, Bergmann AK, et al. DNA methylome analysis in Burkitt and follicular lymphomas identifies differentially methylated regions linked to somatic mutation and transcriptional control. Nat Genet. 2015;2015(47):1316–25.

    Article  CAS  Google Scholar 

  22. Jaramillo D, Hick P, Whittington RJ. Age dependency of nervous necrosis virus infection in barramundi Lates calcarifer (Bloch). J Fish Dis. 2017;40(8):1089–101.

    Article  CAS  PubMed  Google Scholar 

  23. Emmenegger EJ, Sanders GE, Conway CM, Binkowski FP, Winton JR, Kurath G. Experimental infection of six north American fish species with the North Carolina strain of spring viremia of carp virus. Aquaculture. 2016;450:273–82.

    Article  Google Scholar 

  24. Weber M, Hellmann I, Stadler MB, Ramos L, Pääbo S, Rebhan M, et al. Distribution, silencing potential and evolutionary impact of promoter DNAmethylation in the human genome. Nat Genet. 2007;39(4):457–66.

    Article  CAS  PubMed  Google Scholar 

  25. Lister R, Pelizzola M, Dowen RH, Hawkins RD, Hon G, Tonti-Filippini J, et al. Human DNA methylomes at base resolution show widespread epigenomic differences. Nature. 2009;462(7271):315–22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Meissner A, Mikkelsen TS, Gu H, Wernig M, Hanna J, Sivachenko A, et al. Genome-scale DNAmethylation maps of pluripotent and differentiated cells. Nature. 2008;454(7205):766–70.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Li J, Li R, Wang Y, Hu X, Zhao Y, Li L, et al. Genome-wide DNA methylome variation in two genetically distinct chicken lines using MethylC-seq. BMC Genomics. 2015;16:851.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Smith J, Sen S, Weeks RJ, Eccles MR, Chatterjee A. Promoter DNA Hypermethylation and paradoxical gene activation. Trends Cancer. 2020;6:392–406.

    Article  CAS  PubMed  Google Scholar 

  29. Xiu Y, Shao C, Zhu Y, Li Y, Gan T, Xu W, et al. Differences in DNA methylation between disease-resistant and disease-susceptible Chinese tongue sole (Cynoglossus semilaevis) families. Front Genet. 2019;10:847.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Hu Q, Ao Q, Tan Y, Gan X, Luo Y, Zhu J. Genome-wide DNA methylation and RNA analysis reveal potential mechanism of resistance to Streptococcus agalactiae in GIFT strain of Nile Tilapia (Oreochromis niloticus). J Immunol. 2020;204:3182–90.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Jin L, Jiang Z, Xia Y, Lou P, Chen L, Wang H, et al. Genome-wide DNA methylation changes in skeletal muscle between young and middle-aged pigs. BMC Genomics. 2014;15:653.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  32. Murphy PJ, Cairns BR. Genome-wide DNA methylation profiling in zebrafish. Methods Cell Biol. 2016;135:345–59.

    Article  CAS  PubMed  Google Scholar 

  33. Zhang B, Ban D, Gou X, Zhang Y, Yang L, Chamba Y, et al. Genome-wide DNA methylation profiles in Tibetan and Yorkshire pigs under high-altitude hypoxia. J Anim Sci Biotechnol. 2019;10:25.

    Article  PubMed  PubMed Central  Google Scholar 

  34. Yang X, Han H, De Carvalho DD, Lay FD, Jones PA, Liang G. Gene body methylation can alter gene expression and is a therapeutic target in cancer. Cancer Cell. 2014;26(4):577–90.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Arechederra M, Daian F, Yim A, Bazai SK, Richelme S, Dono R, et al. Hypermethylation of gene body CpG islands predicts high dosage of functional oncogenes in liver cancer. Nat Commun. 2018;9(1):3164.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  36. Kim MY, Zilberman D. DNA methylation as a system of plant genomic immunity. Trends Plant Sci. 2014;19(5):320–6.

    Article  CAS  PubMed  Google Scholar 

  37. Bewick AJ, Schmitz RJ. Gene body DNA methylation in plants. Curr Opin Plant Biol. 2017;36:103–10.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Huang X, Zhang X, Zong L, Gao Q, Zhang C, Wei R, et al. Gene body methylation safeguards ribosomal DNA transcription by preventing PHF6-mediated enrichment of repressive histone mark H4K20me3. J Biol Chem. 2021;297(4):101195.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Dai Z, Li J, Hu C, Wang F, Wang B, Shi X, et al. Transcriptome data analysis of grass carp (Ctenopharyngodon idella) infected by reovirus provides insights into two immune-related genes. Fish Shellfish Immunol. 2017;64:68–77.

    Article  CAS  PubMed  Google Scholar 

  40. Zeng ZH, Du CC, Liu SR, Li H, Peng XX, Peng B. Glucose enhances tilapia against Edwardsiella tarda infection through metabolome reprogramming. Fish Shellfish Immunol. 2017;61:34–43.

    Article  CAS  PubMed  Google Scholar 

  41. Jiang M, Chen ZG, Zheng J, Peng B. Metabolites-enabled survival of crucian carps infected by Edwardsiella tardain high water temperature. Front Immunol. 2019;2019(10):1991.

    Article  CAS  Google Scholar 

  42. Krueger F, Andrews SR. Bismark: a flexible aligner and methylation caller for bisulfite-Seq applications. Bioinformatics. 2011;27(11):1571–2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Lister R, Mukamel EA, Nery JR, Urich M, Puddifoot CA, Johnson ND, et al. Global epigenomic reconfiguration during mammalian brain development. Science. 2013;341(6146):1237905.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  44. Park Y, Wu H. Differential methylation analysis for BS-seq data under general experimental design. Bioinformatics. 2016;32(10):1446–53.

    Article  CAS  PubMed  Google Scholar 

  45. Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11(2):R14.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  46. Xie C, Mao X, Huang J, et al. KOBAS 2.0: web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 2011;39:W316–22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28(5):511–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Bio. 2010;11(10):R106.

    Article  CAS  Google Scholar 

  50. Clements DN, Wood S, Carter SD, Ollier WE. Assessment of the quality and quantity of genomic DNA recovered from canine blood samples by three different extraction methods. ResVet Sci. 2008;85:74–9.

    CAS  Google Scholar 

  51. Li LC, Dahiya R. MethPrimer: designing primers for methylation PCRs. Bioinformatics. 2002;18:1427–31.

    Article  CAS  PubMed  Google Scholar 

  52. 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:402–8.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

Not applicable

Funding

This work was supported by the National Natural Science Foundation of China (32073017 and U20A2063), the Youth Innovation Promotion Association CAS (2021338), the Strategic Priority Program (A) of CAS (XDA24030203), and State Key Laboratory of Freshwater Ecology and Biotechnology (2019FBZ05). The funding bodies had no role in the design of the study; collection, analysis, and interpretation of data; or in writing the manuscript.

Author information

Authors and Affiliations

  1. State Key Laboratory of Freshwater Ecology and Biotechnology, Institute of Hydrobiology, Chinese Academy of Sciences, Wuhan, 430072, China

    Libo He, Xinyu Liang, Qian Wang, Cheng Yang, Yongming Li, Lanjie Liao, Zuoyan Zhu & Yaping Wang

  2. University of Chinese Academy of Sciences, Beijing, 100049, China

    Libo He, Xinyu Liang & Qian Wang

  3. Innovative Academy of Seed Design, Chinese Academy of Sciences, Beijing, 100101, China

    Yaping Wang

Authors
  1. Libo He
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Xinyu Liang
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Qian Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Cheng Yang
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Yongming Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Lanjie Liao
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Zuoyan Zhu
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Yaping Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

YPW, ZYZ, and LBH conceived the study and supervised the research. LBH conducted the sample preparation, data analysis, bioinformatics analysis, and wrote the manuscript. YML and LJL contributed to the fish culture. XYL, QW, and CY performed the DNA/RNA isolation and data collection. YPW and ZYZ edited the manuscript. All the authors approved the final version of the manuscript and the decision to submit the manuscript for publication. The author(s) read and approved the final manuscript.

Corresponding authors

Correspondence to Libo He or Yaping Wang.

Ethics declarations

Ethics approval and consent to participate

Animal welfare and experimental procedures were performed in accordance with the Guide for the Care and Use of Laboratory Animals (Ministry of Science and Technology of China, 2006), and the protocol was approved by the committee of the Institute of Hydrobiology, CAS. All surgeries were performed under MS-222 anesthesia, and efforts were made to minimize the suffering of the fish.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

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

Supplementary Information

Additional file 1: Fig. S1.

Whole-genome bisulfite sequencing of the FMO and TYO fish groups. A and B Pearson’s correlation analysis (A) and principal component analysis (B) of the six samples from the two groups. C and D Distribution of mCs in different genomic regions (C) and in different gene elements (D). E and F Methylation levels (E) and length distribution (F) of DMRs.

Additional File 2: Fig. S2.

Correlation between DNA methylation and gene expression (A) The methylation levels of different gene groups (with different mRNA expression levels) in the TYO group. B The mRNA expression levels of different gene categories (with different methylation levels) in the promoter region of the TYO group. C mRNA expression levels of different gene categories (with different methylation levels) in the gene body regions of the TYO group. D Correlation between DNA methylation levels and gene expression levels of promoter hyper-DMGs. E Correlation between DNA methylation levels and gene expression levels of gene body hypo-DMGs.

Additional file 3: Table S1.

Methylated cytosines percentage of different context

Additional file 4: Table S2.

Detailed information on DMRs identified between the two fish groups.

Additional file 5: Table S3.

Detailed information on the DMGs between the two fish groups.

Additional file 6: Table S4.

Detailed information on KEGG enrichment analysis of DMGs.

Additional file 7: Table S5.

Oligonucleotide sequences used in the study.

Rights and permissions

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

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

He, L., Liang, X., Wang, Q. et al. Genome-wide DNA methylation reveals potential epigenetic mechanism of age-dependent viral susceptibility in grass carp. Immun Ageing 19, 28 (2022). https://doi.org/10.1186/s12979-022-00285-w

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12979-022-00285-w

Keywords

Immunity & Ageing

ISSN: 1742-4933

Contact us