Genome-Wide DNA Methylation Variations Between Grass Carp with Different Ages

doi:10.21203/rs.3.rs-63484/v1

Download PDF

Research article

Genome-Wide DNA Methylation Variations Between Grass Carp with Different Ages

https://doi.org/10.21203/rs.3.rs-63484/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background: Grass carp is an important farmed fish in China that infected by many pathogens, especially grass carp reovirus (GCRV). Notably, grass carp showed age-dependent susceptibility to GCRV, while the mechanism remains unclear. Herein, we performed a genome-wide survey of differences in DNA methylation and gene expression between five months old grass carp (FMO, sensitive to GCRV) and three years old grass carp (TYO, resistant to GCRV) aim to uncover the mechanism.

Results: Colorimetric quantification revealed global methylation level of TYO fish was higher than that of FMO fish. Whole-genome bisulfite sequencing (WGBS) of two groups revealed 6,214 differentially methylated regions (DMRs) and 4,052 differentially methylated genes (DMGs), with most of DMRs and DMGs showed hypermethylation patterns in TYO fish. Correlation analysis indicated that DNA hypomethylation in promoter negative correlated with gene expression, whereas positive correlation was found between gene-body DNA hypermethylation and gene expression. Enrichment analysis revealed that promoter hypo-DMGs in TYO fish were significant enriched in pathways involved in immune response while gene-body hyper-DMGs in TYO fish were significant enriched in terms related to RNA transcription, biosynthetic, and energy production. RNA-seq indicated these terms or pathways involved in immune response, biosynthetic, and energy production also significant enriched for the up-regulated genes in TYO fish.

Conclusions: Collectively, these results revealed the genome-wide DNA methylation variations between grass carp with different ages. DNA methylation and gene expression variations in genes involved in immune response, biosynthetic, and energy production may contributed to the age-dependent susceptibility to GCRV in grass carp. Our results will provide important information for the disease-resistant breeding programs of grass carp and may also benefit to the research of age-dependent diseases in human.

Epigenetics & Genomics

Grass carp

Age-dependent susceptibility

Grass carp reovirus

Genome-wide DNA methylation

RNA-seq

immune response

biosynthetic

energy production

Grass carp (Ctenopharyngodon idellus) is an important farmed fish in China, which accounts for more than 18% of total freshwater aquaculture production in this country. Production of grass carp reached 6.01 million tons in 2016, implying the great commercial value of this fish [1]. Nevertheless, the grass carp hemorrhagic disease that caused by grass carp reovirus (GCRV) is an important threat for the aquaculture of grass carp [2]. Consequently, grass carp and the virus are of particular interest to fish breeding geneticists aiming to identify strategies for disease-resistant breeding [3–7]. Our previous study revealed that age is an important factor determined the susceptibility to GCRV infection. Grass carp no more than one year old was sensitive to GCRV, while the fish over three years old was resistant to the virus [9]. Understanding the mechanism will be of benefit to grass carp disease-resistant breeding programs.

In mammals, age is also a crucial factor that affects disease outcomes [10]. Pseudorabies virus (PRV) infection caused more severe clinical disease in newborns than in older individual [11]. Infants and young children are more sensitive to the reovirus induced encephalitis than adults [12]. Mice display an age-dependent acceleration of mortality to Influenza virus Influenza virus (IAV) infection and are useful to model human aging and the outcomes to IAV infection [13]. Further research suggested that the maturation of the innate immune system and the expression of type I IFN genes may contribute to age-related susceptibility to virus infection [14], whereas why the these genes existed different expression patterns between young and adults remain unclear.

The age-dependent susceptibility to virus infection was also reported in several fish. Nervous necrosis virus (NNV) infection in barramundi (Lates calcarifer) caused nervous necrosis at 3 and 4 weeks of age, whereas developed subclinical symptom when fish age was more than 5 weeks [15]. Spring viremia of carp virus (SVCV) challenge in north American fish species revealed that fish of younger age classes were more vulnerable to SVCV infection than older fish [16]. Nevertheless, the mechanism in fish remains largely uncovered.

DNA methylation is an important mediator of gene expression regulation and many genes were reported to be influenced by DNA methylation at different ages. TAP1 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 piglets resistant to Escherichia coli F18 at 35 days of age [17]. Genome-wide DNA methylation analysis of CD4 + and CD8 + T cells from younger and older individuals revealed increased number of methylation changes and higher methylome variation in CD8 + T cells with age, implying the link between age-related epigenetic changes and impaired T cell function [18]. Genome-wide DNA methylation survey of skeletal muscle between young and old healthy human showed a dynamic inter-relationship between DNA methylation, gene expression, age, and exercise [19]. Therefore, it seems like that DNA methylation play an important role in age-related gene expression and is thought to underlie many age-related physiologic phenomena, for example, the age-dependent susceptibility to virus infection.

In the study, we performed a genome-wide survey of differences in DNA methylation and gene expression in spleen tissues between two age stages in grass carp: five months old grass carp (sensitive to GCRV) and three years old grass carp (resistant to GCRV). Global DNA methylation patterns and gene expression profiles were investigated between two fish with different ages. The differentially methylated regions (DMRs) and differentially expressed genes (DEGs) were identified between two groups. Moreover, the correlation between DNA methylation and gene expression were investigated. We believe our results will provide important information for the disease-resistant breeding of grass carp and may also benefit to the research of age-dependent diseases in human.

Global methylation profile and whole-genome bisulfite sequencing

Epigenetic modifications, such as DNA methylation level, have been shown to progressively accumulate during growth and then decreased in the ageing process [20]. Therefore, we analyzed the global DNA methylation status of spleen samples from two groups by colorimetric quantification of 5-methylcytosine. As shown in Fig. 1a, the global methylation level of TYO fish was higher than that of FMO fish, although this difference is not significantly. The gradually increase of global DNA methylation level in TYO fish was in accordance with the reports in human, implying the potential role of DNA methylation in age-related gene regulation.

To further reveal the mechanism underlying age-dependent susceptibility to GCRV in grass carp, we performed whole-genome bisulfite sequencing (WGBS) on spleen samples collected from FMO and TYO fish at base-pair resolution. Three duplicates were processed for each group, yielding a total of 6 libraries. All libraries were sequenced on an Illumina Hiseq X Ten platform to generate 150 bp pair-end reads. Each library yielded > 33.68 GB clean base and > 23.06 × mean coverage depth. The bisulfite conversion rates for all libraries were over 99.88%. Moreover, for all libraries, more than 82.78% genomic sites were covered by at least five unique reads and more than 70.70% sites were covered by at least ten unique reads (Table 1). Pearson correlation analysis and principal component analysis of the 6 libraries showed that samples from each group clustered together (additional file 1: Figure S1). Collectively, these results confirmed the high quality and repeatability of the WGBS data and suitability for further analysis. The sequencing data have been deposited in the Sequence Read Archive (SRA) at the National Center for Biotechnology Information (NCBI) (accession number: PRJNA638254).

Table 1

Summary of WGBS data
Sample	duplicates	clean reads	clean bases(G)	Mapped reads	Unique Mapping rate (%)	BS conversion rate (%)	Mean coverage (×)	1× coverage (%)	5× coverage (%)	10× coverage (%)
FMO	FMO-a	131443366	37.55	91721180	69.78	99.88	27	91.6	85.56	76.16
	FMO-b	129686319	36.75	78836313	60.79	99.90	23.06	91.07	82.78	70.7
	FMO-c	117516406	33.68	82684543	70.36	99.91	23.86	91.85	86.2	75.31
TYO	TYO-a	126812227	36.45	95438882	75.26	99.89	28.5	91.57	84.77	74.1
	TYO-b	133532050	38.15	90027308	67.42	99.88	26.83	91.63	83.57	70.99
	TYO-c	121522621	34.73	84701266	69.7	99.89	25.31	91.02	83.57	71.52

In the detected methylation sites, the average methylated cytosines (mCs) percentage of whole genomic cytosines was 4.99% and 4.80% for FMO and TYO fish, respectively. In the CG contexts, about half of cytosines were methylated (51.42% in FMO fish and 49.40% in TYO fish), while the methylation rates of the cytosines in CHG and CHH contexts (where H is A, C, or T) were on more than 0.10% and 0.09%, respectively (additional file 2: Table S1). Among the mC sites that identified, more than 98% were in the CG contexts, while no more than 0.4% and 1.5% in CHG and CHH contexts (Fig. 1b). The genome sequence preference proximal to the sites of methylated CG (mCG), methylated CHG (mCHG), and methylated CHH (mCHH) contexts were also analyzed. It was found that no sequence preference in the mCG-flanking regions or upstream of mCHG and mCHH contexts; however, the base following mCHG and mCHH contexts was almost always an adenine, following by thymine, while cytosine was observed less often (Fig. 1c). Moreover, the methylation level (ML), defined as the proportion of reads covered each mCs relative to the total reads covered the sites, was calculated. Results showed that more than 80% of the mC sites displayed high ML (ML ≥ 70%) (Fig. 1d).

Gene methylation profile

To characterize methylation profile of grass carp genes, the relative ML in the contexts of gene regions and of their upstream and downstream regions was calculated. In general, the relative ML of mCGs in the gene-body regions was higher than in the 5 ’ upstream and 3 ’ downstream regions in both group, while the ML in mCHG and mCHH was more complex (Fig. 2a). Interestingly, a sharp decrease of 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 that including promoter, 5’ untranslational regions (UTR), exon, intron, 3’ UTR, and repeat. Results suggested that promoter and 5’ UTR had relatively low ML in mCG, following by the exon, while the repeat had the highest ML (Fig. 2b, 2c). Moreover, the heat map of genomic elements suggested that genomic elements from the same group had the same methylation trends and a clear separation was observed between two groups, further implying the repeatability of duplicates and the distinction between fish with different ages (Fig. 2c).

DMRs between FMO and TYO fish

To further investigate the difference between two fish, the data of TYO fish group were compared with FMO fish group, and the DMRs were identified by using the DSS software. Due to most of mC sites were occurred in CG contexts, thereby only DMRs in CG contexts was considered for further study (Fig. 1b). A total of 6,214 DMRs were identified between the two fish groups, including 5,261 hyper-DMRs (hypermethylation in TYO fish compared with FMO fish) and 953 hypo-DMRs (hypomethylation in TYO fish compared with FMO fish) (additional file 3: Table S2). Apparently, the number of hyper-DMRs was more than that of hypo-DMRs. Interestingly, we also found that the mean ML of DMRs in TYO fish was higher than FMO fish (additional file 4: Figure S2a).The heat map of DMRs also showed similar result (Fig. 3a). The length distribution of DMRs was analyzed, and results showed that most of them were shorter than 400 bp (additional file 4: Figure S2b). Circos analysis in a circular layout revealed that the DMRs distributed uniformly in the genome, except the supercontigs C101000030 (Fig. 3b), in which more DMRs were identified. We identified a total of 4,052 DMGs, including 2,855 DMGs that harbored DMRs in gene-body regions (from TSS to TES) and 1,197 DMGs contained DMRs in promoter regions (additional file 5: Table S3).

Functional enrichment analysis of DMGs

To investigate the possible role of genes that showed differential methylation status, GO enrichment analysis was performed. For the gene-body hyper-DMGs, the significant enriched GO terms were related to “regulation”, such as regulation of transcription (DNA-templated), regulation of nucleic acid-templated transcription, regulation of RNA biosynthetic process, and regulation of RNA metabolic process, implying the DMGs involved in RNA transcription and biosynthesis. Whereas the gene-body hypo-DMGs were enriched in GO terms related to modulation of cell apoptosis or cell death, but the enrichment are not significantly. The promoter hyper-DMGs were also enriched in “regulation” related terms, such as regulation of RNA biosynthetic process, regulation of RNA metabolic process, and regulation of nucleobase-containing compound metabolic process, while the promoter hypo-DMGs were not significant enriched in any GO terms. The top five significant enriched GO terms were listed in Table 2 and the detail information was shown in additional file 6: Table S4.

Table 2

GO enrichment analysis of DMGs (Top five terms)
Categories	GO terms	Corrected P Value	DMR genes
gene-body hyper-DMGs	regulation of transcription, DNA-templated	1.44E-07	262
	regulation of nucleic acid-templated transcription	1.44E-07	262
	regulation of RNA biosynthetic process	1.44E-07	262
	regulation of RNA metabolic process	1.44E-07	262
	regulation of nucleobase-containing compound metabolic process	1.44E-07	263
gene-body hypo-DMGs	modulation by virus of host apoptotic process	0.21255	3
	modulation of programmed cell death in other organism	0.21255	3
	modulation of apoptotic process in other organism	0.21255	3
	modulation by symbiont of host programmed cell death	0.21255	3
	modulation by symbiont of host apoptotic process	0.21255	3
promoter hyper-DMGs	regulation of RNA biosynthetic process	0.00084731	118
	regulation of RNA metabolic process	0.00084731	118
	cellular macromolecule metabolic process	0.00084731	280
	macromolecule metabolic process	0.00084731	308
	regulation of nucleobase-containing compound metabolic process	0.00084731	118
promoter hypo-DMGs	regulation of cell shape	1	2
	regulation of cell morphogenesis	1	2
	electron transporter, transferring electrons within the cyclic electron transport pathway of photosynthesis activity	1	2
	tumor necrosis factor receptor binding	1	2
	tumor necrosis factor receptor superfamily binding	1	2

In addition, KEGG enrichment analysis was also carried out. Focal adhesion is the most significant enriched pathway of gene-body hyper-DMGs, following by ECM-receptor interaction, melanogenesis, dorso-ventral axis formation, and adherens junction. For the gene-body hypo-DMGs, Melanogenesis, GnRH signaling pathway, ubiquitin mediated proteolysis, vascular smooth muscle contraction, ubiquinone and other terpenoid-quinone biosynthesis were the top five enriched pathways. Interestingly, immune related pathways, such as herpes simplex infection, cytokine-cytokine receptor interaction, and RIG-I-like receptor signaling pathway, were significant enriched in the promoter hypo-DMGs, whereas the biosynthesis and metabolism related pathways were the predominant enriched pathways in the promoter hyper-DMGs (Fig. 4). The detail information of the KEGG enrichment analysis was shown in additional file 7: Table S5.

Gene expression profiles between two groups

The gene expression profiles between two groups were investigated by RNA-seq using the same samples that used for WGBS. Three duplicates were performed for each group and then obtained a total of 6 libraries. All libraries were sequenced on an Illumina Hiseq X Ten platform to generate 150 bp pair-end reads. Each library generated ≥ 7.5 GB clean data and showed Q30 ≥ 94% (additional file 8: Table S6), implying the sufficient quality of RNA-seq data and suitability for further analysis. The sequencing data have been deposited in the SRA at the NCBI (accession number: PRJNA634937). A total of 9,319 DEGs (additional file 9: Table S7) were identified when data of two groups were compared, including 4,639 up-regulated and 4,680 down-regulated genes. The DEGs were selected for annotation in order to assign the potential function. For the up-regulated DEGs, GO enrichment analysis revealed that catalytic activity was the most significant enriched GO terms and also had the most number of enriched genes. In addition, many terms belonged to the cellular component category, such as intracellular part, cell part, and cytoplasmic part, were significant enriched. The down-regulated DEGs were enriched in GO terms related to signal transduction, enzyme activity, and cell response to stimulation (additional file 10: Table S8). KEGG enrichment analysis showed that metabolism, biosynthesis, and energy production related pathways, such as oxidative phosphorylation, proteasome, ribosome biogenesis in eukaryotes, and metabolic pathway, were significant enriched in the up-regulated DEGs. Moreover, some pathways participated in host defense (lysosome and phagosome), were also enriched. Focal adhesion, ECM-receptor interaction, phosphatidylinositol signaling system, vascular smooth muscle contraction, and adherens junction were the top five enriched pathways in the down-regulated DEGs, following by several growth and development related signaling pathways, such as calcium signaling pathway, mTOR signaling pathway, and VEGF signaling pathway (additional file 11: Table S9 and additional file 12: Figure S3).

Correlation between DNA methylation and gene expression

To analyze the correlation between DNA methylation and gene expression, the DNA methylation level and gene expression level of the gene-body regions, 5’ upstream, and 3’ downstream in two fish groups were analyzed. As shown in additional file 13: Figure S4, the DNA methylation and gene expression of gene-body regions, 5’ upstream, and 3’ downstream in both groups were presented. Furthermore, we divided genes into four groups according to the mRNA expression level: none, low, medium, and high, and the DNA methylation level in different groups were compared. In general, the methylation level in the four groups began to decrease from 2 kb upstream of the TSS of the genes while was increased downstream the TSS. In the 5’ upstream, the genes with high mRNA expression level showed lowest DNA methylation level, whereas the unexpressed genes displayed the highest methylation level (Fig. 5a and additional file 14: Figure S5a). In contrast, the correlation between gene-body or 3’ downstream methylation and mRNA expression was more complex. In addition, we also classified genes into five categories based the methylation level, from the bottom 20% to the top 20%, corresponding to the 1st to 5th groups, respectively, and mRNA expression level of these groups were investigated. It could be seen that 1st group genes in the promoter regions showed the highest expression level, whereas mRNA expression of other groups was complex (Fig. 5b and additional file 14: Figure S5b). Meanwhile, the 5th genes in the gene-body regions presented the highest expression level, while genes in other groups shared similar expression patterns (Fig. 5c and additional file 14: Figure S5c). Therefore, based on the above findings, it could be speculated that DNA hypomethylation in promoter regions and DNA hypermethylation in gene-body regions promoted gene expression.

DNA methylation and mRNA expression of DEGs and DMRs

The DEGs that obtained by RNA-seq were chose for further analysis to reveal the DNA methylation level of DEGs between two fish group. As shown in Fig. 5d, the DNA methylation level of gene-body regions of all DEGs in TYO fish group was slightly higher than that in FMO fish group, while DNA methylation level of promoter of DEGs in two groups was similar (Fig. 5e). Furthermore, the DMRs that obtained by WGBS were selected and the relationship between DNA methylation level and mRNA expression level were analyzed. For the DMRs in the promoter regions, both hyper-DMRs and hypo-DMRs displayed slightly negative correlation between DNA methylation level and mRNA expression level, and this correlation was more evidently in the promoter hypo-DMRs (Fig. 6a and additional file 15: Figure S6a). However, for the DMRs in the gene-body regions, only the hyper-DMRs showed slightly positive correlation (r = 0.22, p < 0.05) (Fig. 6b), while no correlation was observed in the hypo-DMRs (p > 0.05) (additional file 15: Figure S6b). These findings further suggested that hypomethylation in promoter regions and hypermethylation in gene-body regions involved in transcriptional activation.

Genes involved in age-dependent susceptibility to GCRV

Enrichment analysis revealed gene-body hyper-DMGs were enriched in terms or pathways related to RNA transcription, biosynthesis, and energy metabolism. Moreover, it was proposed that gene-body hypermethylation promoted gene expression (Fig. 5c). We then investigated the mRNA expression pattern of genes involved in these terms or pathways to further confirmed the standpoint. As expected, for genes involved in protein biosynthesis and energy production, most of them were up-regulated in TYO fish (Fig. 6c, 6d). The promoter hypo-DMGs were enriched in the classical innate immune pathways and hypomethylation in promoter regions was thought to involve in transcriptional activation (Fig. 5b). We therefore analyzed the expression profiles of immune-related genes, and results showed that most of them presented up-regulation patterns in TYO fish (Fig. 6e). Genes involved in growth and development was also investigated, while most of them showed the opposite trend, down-regulated in TYO fish group (Fig. 6f), implying the FMO fish was concentrated in cell growth and development. It is well known that immune response play an important role in host defense against pathogen. Protein synthesis and energy production are also benefit to organism resistant to virus infection. Therefore, the DNA methylation and gene expression variations of genes in these pathways may responsible for the age-dependent susceptibility to GCRV in grass carp with different ages.

DNA methylation, the most widely and well-understood type of epigenetic modification, was reported to play a crucial role in transcriptional regulation [21, 22]. In general, DNA hypomethylation in gene promoters is usually linked to transcriptional activation while hypermethylation in these regions is considered to correlate with transcriptional repression [23, 24]. In fish, DNA methylation was reported to participate in skin color variations in crucian carp [25], phenotypic differences, environmental adaptation, and sex chromosome evolution in threespine stickleback [26–28], behavioral effects in zebrafish after bisphenol A exposure [29], sex determination in half-smooth tongue sole [30], and sex-specific phenotypes in hybrid tilapia [31]. Moreover, it was also considered that DNA methylation was correlated with the disease resistance trait in grass carp and Chinese tongue sole [32, 33]. Grass carp showed age-dependent susceptibility to GCRV. Therefore, we hypothesized that DNA methylation may involve in the disease resistance ability of three years old fish. Here, we reported for the first time a genome-wide survey of DNA methylation status of spleen DNA from five months old grass carp versus three years old grass carp and revealed the correlation between DNA methylation and mRNA expression. We believed that our results will benefit to the disease-resistant breeding of grass carp and may also provide useful information the research of age-dependent diseases in human.

We performed WGBS on spleen DNA from FMO and TYO grass carp. To confirm the reliability and repeatability of sequencing data, three duplicates were performed for each group. Therefore, the size and coverage depth of sequencing data obtained in our study was more than that reported in chicken, pigs, and Chinese tongue sole [24, 33, 34]. Interestingly, we found that general DNA methylation pattern of grass carp is consistent with other species [19, 35, 36]. For example, most of mC sites were occurred in the CG contexts, the relative low ML in the promoter or 5’ UTR regions and high ML in repeat regions, negative correlation between promoter or 5’ UTR regions methylation and gene expression. These results implied the conserved role of DNA methylation during evolution. The specific roles of DNA methylation in gene-body regions remain unclear [37, 38]. Interestingly, we showed genes with highest methylation level in the gene-body regions presented the highest expression level in both groups (Fig. 5c and Figure S5C), suggesting that hypermethylation in gene-body regions positive correlated with mRNA expression. Similar results of positive correlation were also observed in arabidopsis, chicken, and human [24, 39, 40]. However, further study of DMRs revealed that negative correlation in promoter regions and positive correlation in gene-body regions are not obviously, probably due to the DNA methylation is only one of the factors that affect gene expression [34, 35].

Colorimetric quantification of 5-mc showed global methylation level of TYO fish was higher than that of FMO fish (Fig. 1a). Moreover, we identified a total of 6,214 DMRs, in which most of them were hypermethylated in TYO fish. Meanwhile, the ML of DMRs in TYO fish was higher than FMO fish (Figure S2A). These results may suggest that global DNA hypermethylation in TYO fish compared with FMO fish. Previous study showed that DNA methylation variations occur throughout the lifetime. Global DNA methylation levels increased over the first few years of life and then decreased beginning in late adulthood [41, 42]. Grass carp have a lifespan more 20 years and sexual maturity period of grass carp is about 4 ~ 5 years. The TYO grass carp is like the adolescent of human while the FMO fish is similar the infant or childhood of human beings. Therefore, it is understandable that DNA methylation level was higher in TYO fish. Similar result also reported in mouse, which blood showed a general pattern of epigenome-wide hypermethylation with age [43]. The DNA hypermethylation in TYO fish may be involved in the epigenetic reprogramming during development, which is important in regulate gene expression [44, 45].

To reveal the potential role of DMGs, GO and KEGG analysis were performed. Results showed gene-body hyper-DMGs enriched in terms related to RNA transcription, biosynthesis, and metabolism, while gene-body hypo-DMGs were enriched in pathways involved in cell growth, proliferation, and differentiation (such as melanogenesis, GnRH signaling pathway, and calcium signaling pathway). Considering the DNA methylation in gene-body regions, especially DNA hypermethylation, positive correlated gene expression level, therefore, it could be speculated that TYO fish had superiority in biosynthesis and energy production, while the FMO fish were concentrated in cell growth and development. Similar results of enrichment were also presented in the DEGs that obtained by RNA-sEq. The gene-body hyper-DMGs or up-regulated DEGs involved in RNA transcription, biosynthesis, and energy production suggested that the TYO fish is vigorousness, which may be benefit to resistant to virus infection [46, 47]. The promoter hypo-DMGs were significant enriched in classical immune related pathways, such as herpes simplex infection, cytokine-cytokine receptor interaction, apoptosis, and RIG-I-like receptor signaling pathway. Regarding the DNA methylation in promoter regions, particularly DNA hypomethylation, negative correlated with mRNA expression level, we then supposed that the immunologic function of TYO fish was more perfect than FMO fish, which is important in defense against pathogen infection. Coincidently, RNA-seq also revealed that some immune related terms were also significant enriched in the up-regulated DEGs but not in the down-regulated DEGs.

In conclusion, global DNA methylation patterns and gene expression profiles between two fish groups with different ages were investigated. We identified 6,214 DMRs and 4,052 DMGs, with most of DMRs and DMGs showed hypermethylation patterns in TYO fish. DNA hypomethylation in promoter regions negative correlated with gene expression while DNA hypermethylation in gene-body regions positive correlated with gene expression. DNA methylation and gene expression variations in genes involved in immune response, biosynthetic, and energy production may contributed to the age-dependent susceptibility to GCRV in grass carp.

Experimental fish and sample collection

A total number of 100 fish, which contained 50 five months old (FMO) and 50 three years old (TYO) grass carp, on behalf of the sensitive and resistant fish to GCRV infection, respectively, were used in the study. All the fish were obtained from 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. Fish were fed with commercial feed twice a day and water was exchanged daily. If no abnormal symptoms were observed, the fish were selected for further study. In each group, 15 fish that represented three biological duplicates (n = 5 for each duplicate) were random collected, the fish were anesthesia by MS-222 at the concentration of 100 mg/L, and then fish were dissected and spleens were removed rapidly for analysis. The obtained samples were stored at − 80 °C until DNA or RNA extraction. After sample collection, fish were euthanized by immediately cutting off the spinal cord adjacent to the head. All of the experiments involving animals were carried out in accordance with the guide for the care and use of laboratory animals (Ministry of Science and Technology of China), and the protocol was approved by the committee of the Institute of Hydrobiology, CAS. All surgery was performed under MS-222 anesthesia, and all efforts were made to minimize suffering of fish.

Whole-genome bisulfite sequencing

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

Data analysis

Raw data reads in fastq format were initially processed using in-house perl scripts. In this step, clean data (clean reads) were obtained by removing adapter, poly-N and poor quality data. The Q20, Q30, and GC contents of the clean data were calculated, and all downstream analysis was performed using the clean high quality data. The reference genome of grass carp and clean reads were transformed into bisulfite-converted sequences (C-to-T and G-to-A converted). Then, bismark software (version 0.16.3) was used to perform alignments of bisulfite-treated reads to the reference genome [48]. Sequence reads that produce a unique best alignment from the two alignment processes (original top and bottom strand) are then compared to the normal genomic sequence and the methylation state of all cytosine positions in the read was inferred. The same reads that aligned to the same regions of genome were regarded as duplicated ones. The sequencing depth and coverage were summarized using deduplicated 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.

In order to calculate the methylation level of the sequence, we divided the sequence into multiple bins, with size is 10 kb. The sum of methylated and unmethylated read counts in each window/bin was calculated. 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 the unmethylated Cs. Calculated ML was further corrected with the bisulfite non-conversion rate according to previous studies [45].

Differentially methylated regions analysis

Differentially methylated regions (DMRs) analysis of two groups/conditions was identified using the DSS software [49], in which a sliding-window method was used. The window was set as 1,000 bp and step length to 100 bp, and at least 10 informative CpGs in windows were considered. Fisher exact test was used to detect the DMRs. Differentially methylated genes (DMGs) were identified as genes whose gene-body regions (from transcription start site (TSS) to transcription end site (TES)) or promoter regions (upstream 2 kb from the TSS) have an overlap with the DMRs. GO and KEGG enrichment analysis of the DMGs were implemented by the GOseq R package and KOBAS software [50, 51].

RNA -seq

RNA of spleens that collected above was isolated using TRizol reagent (Invitrogen, USA) according the manufacturer’s protocol. RNA of sufficiently high quality was used in 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 Hiseq X Ten platform and 150 bp pair-end reads were generated. The output raw data reads were processed as described previously in order to get the clean data [4]. The clean reads were mapped to the reference genome of grass carp by Hisat2 software [52] and gene expression levels were calculated by FPKM (expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced) method [53]. Differential expression analysis of two groups/conditions was performed using the DESeq package [54]. Genes with an adjusted P-value < 0.05 (q value < 0.05) in DESeq analysis were assigned as differentially expressed genes (DEGs). GO and KEGG enrichment analysis of DEGs were performed by GOseq R package and KOBAS software [50, 51].

Global DNA methylation measurement and bisulfite sequencing PCR

Genomic DNA from the spleens of two groups was extracted using the traditional phenol–chloroform protocol with RNase treatment. Global DNA methylation level was measured with 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 duplicates were performed for each group. Data were shown as the mean ± standard deviation of three replicates.

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 years old; UTR: Untranslational regions; WGBS: Whole-genome bisulfite sequencing.

Ethics approval and consent to participate

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

Consent for publication

Not applicable.

Availability of data and materials

The WGBS and RNA-seq data reported in the study have been deposited in the Sequence Read Archive (SRA) at the National Center for Biotechnology Information (NCBI) (accession number: PRJNA638254 and PRJNA634937). Other data generated or analysed during this study are included in additional files.

Competing interests

The authors declare that they have no competing interests

Funding

This work was supported by the National Natural Science Foundation of China (Grant NO. 31702332 and 31721005) and the Strategic Priority Program (A) of Chinese Academy of Science (Grant NO. XDA24030203). The funding bodies had no role in study design, data collection, analysis and interpretation, decision to public, or writing of the manuscript.

Authors' contributions

YW, ZZ, and LH designed, supervised, and coordinated the study. LH performed the study, analyzed and interpreted the data, and wrote the manuscript. YL contributed to fish culture. DZ, PC, RH, and LL performed DNA/RNA isolation and data collection. YW and ZZ revised the manuscript. All authors gave approval to the final version of the manuscript and the decision to submit the work for publication.

Acknowledgements

Not applicable.

[1] Food and Agriculture Oranization of the United Nations. Fishery and aquaculture statistics yearbook. Rome; 2016.

[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.

[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(6): 625-31.

[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.

[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-585.

[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.

[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.

[9] 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-712.

[10] Chan YH, Ng LFP. Age has a role in driving host immunopathological response to alphavirus infection. Immunology. 2017, 152(4):545-555.

[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.

[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-1066.

[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.

[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-334.

[15] 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-1101.

[16] 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-282.

[17] 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.

[18] 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.

[19] 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.

[20] Jones MJ, Goodman SJ, Kobor MS. DNA methylation and healthy human aging. Aging Cell. 2015; 14(6): 924-32.

[21] 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.

[22] 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.

[23] 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.

[24] 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.

[25] Zhang Y, Liu J, Fu W, Xu W, Zhang H, Chen S, et al. Comparative Transcriptome and DNA methylation analyses of the molecular mechanisms underlying skin color variations in Crucian carp (Carassius carassius L.). BMC Genet. 2017; 18(1):95.

[26] Smith G, Smith C, Kenny JG, Chaudhuri RR, Ritchie MG. Genome-wide DNA methylation patterns in wild samples of two morphotypes of threespine stickleback (Gasterosteus aculeatus). Mol Biol Evol. 2015; 32(4): 888-95.

[27] Artemov AV, Mugue NS, Rastorguev SM, Zhenilo S, Mazur AM, Tsygankova SV, et al. Genome-Wide DNA Methylation Profiling Reveals Epigenetic Adaptation of Stickleback to Marine and Freshwater Conditions. Mol Biol Evol. 2017; 34(9): 2203-2213.

[28] Metzger DCH, Schulte PM. The DNA Methylation Landscape of Stickleback Reveals Patterns of Sex Chromosome Evolution and Effects of Environmental Salinity. Genome Biol Evol. 2018; 10(3):775-785.

[29] Olsvik PA, Whatmore P, Penglase SJ, Skjærven KH, Anglès d'Auriac M, Ellingsen S. Associations Between Behavioral Effects of Bisphenol A and DNA Methylation in Zebrafish Embryos. Front Genet. 2019; 10: 184.

[30] Shao C, Li Q, Chen S, Zhang P, Lian J, Hu Q, et al. Epigenetic modification and inheritance in sexual reversal of fish. Genome Res. 2014; 24(4): 604-15.

[31] Wan ZY, Xia JH, Lin G, Wang L, Lin VC, Yue GH. Genome-wide methylation analysis identified sexually dimorphic methylated regions in hybrid tilapia. Sci Rep. 2016; 6: 35903.

[32] Shang X, Yang C, Wan Q, Rao Y, Su J. The destiny of the resistance/susceptibility against GCRV is controlled by epigenetic mechanisms in CIK cells. Sci Rep. 2017; 7(1): 4551.

[33] 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.

[34] 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.

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

[36] 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.

[37] Holliday R, Pugh JE. DNA modification mechanisms and gene activity during development. Science. 1975; 187(4173): 226-32.

[38] Zemach A, McDaniel IE, Silva P, Zilberman D. Genome-wide evolutionary analysis of eukaryotic DNA methylation. Science. 2010; 328(5980): 916-9.

[39] Cokus SJ, Feng S, Zhang X, Chen Z, Merriman B, Haudenschild CD, et al. Shotgun bisulphite sequencing of the Arabidopsis genome reveals DNA methylation patterning. Nature. 2008; 452(7184): 215-9.

[40] Laurent L, Wong E, Li G, Huynh T, Tsirigos A, Ong CT, et al. Dynamic changes in the human methylome during differentiation. Genome Res. 2010; 20(3): 320-31.

[41] Gensous N, Bacalini MG, Franceschi C, Meskers CGM, Maier AB, Garagnani P. Age-Related DNA Methylation Changes: Potential Impact on Skeletal Muscle Aging in Humans. Front Physiol. 2019; 10: 996.

[42] 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-185.

[43] Kochmanski J, Marchlewicz EH, Cavalcante RG, Sartor MA, Dolinoy DC. Age-related epigenome-wide DNA methylation and hydroxymethylation in longitudinal mouse blood. Epigenetics. 2018; 13(7): 779-792.

[44] Wang X, Bhandari RK. DNA Methylation Dynamics During Epigenetic Reprogramming of Medaka Embryo. Epigenetics. 2019; 14(6): 611-622.

[45] 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.

[46] Li Y, Jia A, Wang Y, Dong L, Wang Y, He Y, et al. Immune effects of glycolysis or oxidative phosphorylation metabolic pathway in protecting against bacterial infection. J Cell Physiol. 2019; 234(11): 20298-20309.

[47] Frederick M, Skinner HD, Kazi SA, Sikora AG, Sandulache VC. High expression of oxidative phosphorylation genes predicts improved survival in squamous cell carcinomas of the head and neck and lung. Sci Rep. 2020; 10(1): 6380.

[48] Krueger F, Andrews SR. Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics. 2011; 27(11): 1571-2.

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

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

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

[51] 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-322.

[53] 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.

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

Download PDF

Version 1

posted

You are reading this latest preprint version

Genome-Wide DNA Methylation Variations Between Grass Carp with Different Ages

Status:

Version 1

Abstract

Figures

Background

Results

Global methylation profile and whole-genome bisulfite sequencing

Gene methylation profile

DMRs between FMO and TYO fish

Functional enrichment analysis of DMGs

Gene expression profiles between two groups

Correlation between DNA methylation and gene expression

DNA methylation and mRNA expression of DEGs and DMRs

Genes involved in age-dependent susceptibility to GCRV

Discussion

Conclusion

Methods

Experimental fish and sample collection

Whole-genome bisulfite sequencing

Data analysis

Differentially methylated regions analysis

RNA -seq

Global DNA methylation measurement and bisulfite sequencing PCR

Abbreviations

Declarations

References

Supplementary Files

Status:

Version 1