Hanyuan Zhang
Peng Xu
Yanliang Jiang
Zixia Zhao1
Chuanju Dong
Jian Xu- 1Key Laboratory of Aquatic Genomics, Ministry of Agriculture, CAFS Key Laboratory of Aquatic Genomics and Beijing Key Laboratory of Fishery Biotechnology, Chinese Academy of Fishery Sciences, Beijing, China
- 2Fujian Collaborative Innovation Center for Exploitation and Utilization of Marine Biological Resources, Xiamen University, Xiamen, China
- 3Henan Academy of Fishery Science, Zhengzhou, China
- 4College of Fishery, Henan Normal University, Xinxiang, China
Polyunsaturated fatty acids (PUFAs) are a set of important nutrients that mainly include arachidonic acid (ARA4), docosahexaenoic acid (DHA), eicosapentaenoic acid (EPA), and α-linolenic acid (ALA). Recently, fish-derived PUFAs have been associated with cardiovascular health, fetal development, and improvement of brain functions. Studies have shown that fish muscular tissues are rich in PUFAs, which are influenced by various factors, including genetic variations, regulatory profiles, and methylation status of desaturase genes during fatty acid desaturation and elongation processes. However, the genetic mechanism and the pathways involved in fatty acid metabolism in fishes remain unclear. The overall aim of this study was to assess differences in gene expression responses among fishes with different fatty acid levels. To achieve this goal, we conducted genome-wide association analysis (GWAS) using a 250K SNP array in a population of 203 samples of common carp (Cyprinus carpio) and identified nine SNPs and 15 genes associated with muscular PUFA content. Then, RNA-Seq and whole genome bisulfite sequencing (WGBS) of different groups with high and low EPA, DHA, ARA4, and ALA contents in muscle, liver and brain tissues were conducted, resulting in 6,750 differentially expressed genes and 5,631 genes with differentially methylated promoters. Gene ontology and KEGG pathway enrichment analyses of RNA-Seq and WGBS results identified enriched pathways for fatty acid metabolism, which included the adipocytokine signaling pathway, ARA4 and linoleic acid metabolism pathway, and insulin signaling pathway. Integrated analysis indicated significant correlations between gene expression and methylation status among groups with high and low PUFA contents in muscular tissues. Taken together, these multi-level results uncovered candidate genes and pathways that are associated with fatty acid metabolism and paved the way for further genomic selection and carp breeding for PUFA traits.
Introduction
Aquatic products contribute to a large part of our daily diet, which not only offer high-quality protein but also contain abundant long-chain polyunsaturated fatty acids (PUFA), vitamins, and mineral substances. The fat contents (FAs) of most fishes range from 1 to 4% (Iverson et al., 2002). In general, fish flesh has lower FA but higher PUFA content compared to livestock and poultry meat (Rymer and Givens, 2005; Wood et al., 2008; Henriques et al., 2014; Geay et al., 2015). Omega-6 and omega-3 PUFAs are essential fatty acids (EFAs) that play critical roles in maintaining cell membrane structure. However, as bioactive lipid mediators, these are not synthesized by the human body (Calder, 2013). Long-chain PUFAs, particularly eicosapentaenoic acid (EPA) and docosahexaenoic acid (DHA), positively influence retinal development, neurological development, and cardiovascular system maintenance (Innis, 2008; Guesnet and Alessandri, 2011; Mozaffarian and Wu, 2011). PUFAs can effectively reduce the risk of heart diseases (e.g., coronary heart disease and heart rate variability) (Vedtofte et al., 2011; Enns et al., 2014), prevent the increase of blood viscosity (Zanetti et al., 2015), and act as an anti-glioma agent that promotes the regression and apoptosis of tumors by acting as an intracellular signaling molecule (Witte and Hardman, 2015). Another two PUFAs, namely, α-linolenic acid (ALA) and arachidonic acid (ARA4), have also been reported to be beneficial for human development and health by preventing cardiovascular disease and promoting embryonic and brain development (Davis-Bruno and Tassinari, 2011; Pan et al., 2012).
Fatty acids in fishes originate from any of two sources, namely, synthesis in vivo from non-lipid carbon sources or uptake from dietary lipids. Therefore, muscular PUFA content can be influenced by dietary intake and endogenous metabolism (Davidson, 2013). Various studies have investigated the mechanism underlying PUFA synthesis and degradation, including associated genetic variants, gene expression, and epigenetics (Berton et al., 2016; Chen et al., 2018). The pathways for the synthesis and modification of PUFAs include a variety of enzyme systems. Genetic variants in fatty acid-synthesizing enzymes substantially influence muscular fatty acid levels (Foster, 2012; Beld et al., 2015). For example, delta-5 desaturase (FADS1) and delta-6 desaturase (FADS2) are two enzymes involved in fatty acid metabolism (Ralston et al., 2015). Two haplotypes including the fads1 and fads2 genes have been associated with differences in the synthesis of long-chain PUFAs (Ameur et al., 2012). During fatty acid synthesis, acetyl-CoA serves as the initial specific substrate for the acetyl-CoA carboxylase. Modifications in the microsomal glycerol-3-phosphate pathway in fishes as well as in mammals involve the incorporation of PUFAs into phospholipids and triacylglycerols (Murphy, 2013). During desaturation and elongation processes in freshwater fishes, omega-9 fatty acids could be synthesized by Δ9 desaturase, whereas omega-3 (n-3) and omega-6 (n-6) fatty acids could not be produced due to the lack of the Δ12 and Δ15 desaturase enzymes.
Although PUFAs have a series of important functions, current understanding of the molecular mechanisms related to fish fatty acid metabolism remains limited. QTL mapping and genome-wide association analysis (GWAS) on muscular FA have been reported in the Atlantic salmon, common carp (Cyprinus carpio), and large yellow croaker (Derayat et al., 2007; Sodeland et al., 2013; Kuang et al., 2015; Xiao et al., 2016; Zheng et al., 2016). A previous genome-wide scan identified multiple QTLs located on LG6 and LG23 that influence n-3 PUFA content in the muscle tissues of Asian seabass (Xia et al., 2014). Fatty acid desaturases and elongases (e.g., Δ5 and Δ6 desaturases) of several freshwater and marine fishes have been cloned (Agaba et al., 2004; Zheng et al., 2016). These enzymes play crucial roles in the biosynthesis of the long-chain C20/22 PUFAs from shorter chain C18 PUFAs. It involves a series of complex biosynthesis reactions, including desaturation, elongation and peroxisomal β-oxidation, to convert ALA to n-3 fatty acids EPA and DHA (Kjaer et al., 2016). A previous study demonstrated that there are two alternative pathways for DHA biosynthesis in teleost, (Oboh et al., 2017). One is Sprecher pathway which desaturates 24:5(n-3) to 24:6(n-3) via Δ6 desaturase Fads2, and the other is Δ4 pathway functioned by Δ4 desaturase Fads2. However, the molecular genetic mechanisms of fatty acid-related traits, particularly those involved in important PUFA (e.g., DHA, EPA) pathways, remain unclear.
The common carp has a long breeding history as an important worldwide cultured fish (Bostock et al., 2010). In 2016, the global production of C. carpio reached 4.56 million tons, while 3.50 million tons were cultured in China1. Dozens of C. carpio strains and populations have been cultured, including Yellow River carp, Songpu mirror carp, and Hebao carp. The diverse populations of C. carpio exhibited phenotypic variations in body color, growth, scales, and muscular fatty acid contents. The muscular PUFA content of common carp have been a research topic of interest based on the need to increase aquaculture quality instead of quantity. Based on its economic and scientific importance, genomic resources of C. carpio have been developed and extensively utilized in recent years. The transcriptome and genome assembly of C. carpio have been reported, paving the way for more in-depth investigations (Ji et al., 2012; Xu J. et al., 2012, Xu P. et al., 2014). A high-throughput carp SNP array with 250K SNPs was developed, thereby offering a powerful tool for genetic association studies (Xu J. et al., 2014). Several QTL and GWAS of C. carpio were conducted using the Carp 250K SNP array, including traits of growth and FA (Peng et al., 2016; Zheng et al., 2016). High-throughput sequencing techniques facilitated the multi-omics studies on genome resequencing, genome methylation, and transcriptome analysis on important traits of C. carpio (Jiang et al., 2014; Wang et al., 2014a,b). Gene expression profiles are largely regulated by DNA methylation through transcriptional regulation and chromatin remodeling. The genomic regions with differential methylation levels, known as differentially methylated regions (DMRs), represent the most active regions that may be related to transcriptional regulation. A handful of integrated studies on the expression and epigenetics of various species have been reported. He et al. (2013) investigated transcriptomic and epigenomic variations in maize hybrids and concluded that similar mechanisms may account for the genome-wide epigenetic regulation of gene activity and transposon stability in different organs. Zhang et al. (2015) conducted integrative analysis of transcriptomic and epigenomic data to reveal regulatory patterns for bone mineral density variations, showing consistent association evidence from both mRNA/miRNA expression and methylation data. A recent study conducted comparative transcriptomic and DNA methylation analyses of color trait in Crucian carp and identified several pigmentation-related pathways (Zhang et al., 2017). The above studies indicated the power of multi-omics data analysis. Here, we report our multi-level research on muscle PUFA content traits in common carp using GWAS, RNA-Seq, and methylation analyses. Dozens of fatty acid metabolism-associated candidate genes and pathways were identified by the integration of associated genes, differentially expressed genes (DEGs), and genes in differential methylated regions (DMRs). This study takes a further step to better understand the mechanism of muscular PUFA content and presents potential applications in breeding.
Materials and Methods
Ethics Statement
This study was conducted in accordance with the recommendations of the Care and Use of Animals for Scientific Purposes established by the Animal Care and Use Committee of the Chinese Academy of Fishery Sciences (ACUC-CAFS). The protocol was approved by the ACUC-CAFS. Before the blood and tissue samples were collected, all fishes were euthanized in MS222 solution.
Sample Collection and Phenotypic Measurements
A total of 203 Yellow River 2-year-old carp were randomly selected from a large population that was cultured at the Henan Fishery Research Institute, Henan, China. From each sample, 1 mL of blood was collected in lysis buffer for DNA extraction and genotyping. More than 100 g of dorsal muscle tissue were collected and cryopreserved in dry ice, and then stored at −80°C until fatty acid content determination. More than 5 g of muscle, liver, and brain tissues were acquired and preserved in RNALater (Qiagen, Hilden, Germany) at −80°C for RNA extraction and sequencing. Approximately 20 g of muscle sample per fish were used in measuring total fat and fatty acid contents (N = 12), which were conducted by the Agricultural Products Safety and Quality Supervision Inspection Center in Zhengzhou, Henan according to national standards for food safety. Correlation analysis among traits was conducted using IBM SPSS Statistics 19.0 software, and the correlation matrix plot was drawn using the R package ggplot22.
DNA Extraction, Genotyping, and Quality Control
Genomic DNA was extracted from whole blood using a DNeasy 96 Blood & Tissue Kit (Qiagen, Shanghai, China) following the manufacturer’s protocol. The extracted DNA was quantified using a NanoDrop-1000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, United States). The integrity of DNA was examined on a 1.5% agarose gel by electrophoresis. The final DNA concentration was diluted to 50 ng/μL for genotyping. The total amount of qualified genomic DNA for whole genome sequencing was 2 μg per sample. The common carp 250K SNP array was developed in a previous study using the Affymetrix Axiom genotyping technology (Xu J. et al., 2014). Genotyping was performed by GeneSeek (Lincoln, Nebraska, United States). After genotyping, PLINK v1.93 was used for quality control (Chang et al., 2015). SNPs with low call rate (<95%) or low minor allele frequency (MAF < 5%) were excluded, and samples with <90% genotyping rates were filtered out.
Genome-Wide Association Analysis
Genome-wide association analysis for screening of fat and fatty acid content trait (total fat, 12 individual fatty acids, nine classified fatty acids)-related genes were performed based on genotyping data using the TASSEL version 5.0 software4 (Bradbury et al., 2007). The model “PCA” was used to create a q-matrix, and then the option “wMLM” was used to perform association analysis. The genome-wide significant P-value threshold was adjusted based on Bonferroni correction, and the suggestive threshold was set to 3.445 × 10−5. Associated SNP loci for four PUFA (ALA, ARA4, DHA, and EPA) content values were selected with a corrected P value <0.05. Gene annotation was performed on 10 kb of the flanking regions of the associated SNP loci. The Manhattan plots and Q-Q plots were generated using qqman package of the Comprehensive R Archive Network5.
RNA Sequencing and Differential Gene Expression Analysis
Three tissues (brain, liver, and muscle) of 20 fishes were included in RNA sequencing based on four traits, namely, ALA, ARA4, DHA, and EPA content. From these fishes, four were selected with high content and four were selected with low content in all four traits. Total RNA was extracted from the brain, liver, and muscle tissues using an RNeasy kit (Qiagen, Shanghai, China) following manufacturer’s instructions. The integrity and size distribution of all samples were assessed using a Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, United States). A complementary DNA (cDNA) library was constructed, and high-throughput sequencing was performed using an Illumina HiSeq2500 Sequencing System with paired-end 2 × 150 nucleotide reads (Illumina, San Diego, CA, United States).
Low-quality reads and residual adapter sequences from FASTQ files were filtered and trimmed using Trimmomatic v.0.32 (Bolger et al., 2014). Reads were trimmed when the average Phred score was <20 across four bases in the sliding window. After filtering, reads shorter than 36 bp or single ends were removed. Bowtie2-build indexer (Bowtie2 v.2.3.4.2) (Langmead and Salzberg, 2012) was used to build a Bowtie index from the common carp genome assembly. The filtered reads were aligned to the genome sequence using Tophat2 (Kim et al., 2013). Samtools (Li et al., 2009) was used to index the Tophat2 output bam files, and Cufflinks (Trapnell et al., 2010, 2014) was used to assemble the reconstructed transcripts from the aligned reads using genome and annotation files. These assembled transcript structures were merged into one single dataset with Cuffmerge. Differential expression and regulation at the gene and transcript levels were identified using Cuffdiff. The volcano plots showing gene expression differences were constructed using the ggplot2 package. The Benjamini-Hochberg FDR corrected P value (q value) <0.05 and the absolute log2 fold-change (FC) value >1 were considered differentially expressed genes. The heatmaps of the four PUFAs in three tissues were generated by the FPKM of differentially expressed genes using the pheatmap package. Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using DAVID Bioinformatics Resources 6.8 Tools6 (Huang et al., 2009). Pathways with at least five differentially expressed genes assigned and q value <0.05 were considered enriched. The GO and KEGG plots were constructed using R package clusterProfiler (Yu et al., 2012).
Whole Genome Bisulfite Sequencing and Differential Methylation Analysis
Genomic DNA was extracted from the muscle tissues of the same fishes used in the RNA-Seq study. Genomic DNA was treated using sodium bisulfite, which converts unmethylated cytosine to uracil, then thymine (Wang et al., 2006). Whole genome bisulfite sequencing (WGBS) was performed using an Illumina HiSeq2500 sequencer with 150-bp paired-end sequencing (Illumina, San Diego, CA, United States). The software swDMR7 was used to comprehensively analyze the DMRs from methylation sequencing profiles by a sliding window approach (Wang et al., 2015). The input files were prepared using a WGBS data aligner Bismark (Krueger and Andrews, 2011). The DMR detection and annotation procedures were performed as described below. First, a sliding window with 1,000-bp window size and 100-bp step size was chosen for scanning methylation rates. Second, The Benjamini-Hochberg FDR corrected P value (q value) <0.1 and absolute log2 FC value >1 were considered as potential DMRs. Third, two potential DMRs were merged when their distance was less than the threshold. The merged DMRs were tested by previous steps to guarantee the significance level. This extension step was repeated until P value >0.1. The new extension of potential DMRs were considered as candidate DMRs. Finally, candidate DMRs in the promoter regions were selected for enrichment and correlation analyses. GO and KEGG enrichment analyses were performed using DAVID Bioinformatics Resources 6.8 Tools6 (Huang et al., 2009). The GO and KEGG plots for candidate genes in the DMPs were constructed using R package clusterProfiler (Yu et al., 2012). The DGE and DMP results were correlated using R package ggplot2 (Ginestet, 2011), and transcriptional factor binding site (TFBS) prediction was conducted using AnimalTFDB 3.0 Tools8 (Hu et al., 2018).
Results
Genotyping and Phenotyping of 203 Accessions of C. carpio
On the basis of our previous work (Xu J. et al., 2014; Xu P. et al., 2014), we randomly collected 203 samples from a cultivated population of C. carpio. After DNA extraction and SNP genotyping using Carp 250K SNP array, a raw genotype data of 250,000 SNPs for 203 samples were generated. A total of 193 samples with 108,684 polymorphic SNPs passed the quality control threshold, and 29,026 tag SNPs were chosen by the pruning method (SNPs with LD R2> 0.9 were filtered out) for further association analysis.
We analyzed total FA and 12 fatty acids, then nine classified groups of fatty acids were also calculated (Table 1). Figure 1A shows that oleic acid and linoleic acid comprise 60% of all fatty acids, and the content of the PUFAs (ALA, ARA4, DHA, and EPA) are relatively low. The 22 phenotypes were also analyzed to uncover potential relationships among phenotypes, and the correlation matrix is shown in Figure 1B and Supplementary Table S1. We found that ALA, ARA4, DHA and EPA contents were significantly correlated with each other, possibly suggesting shared or related genes regulating these traits. The correlations between ALA and other three fatty acids (ARA4, DHA, EPA) were moderate (r = 0.398, 0.434, 0.306, respectively; p < 10−8, p < 10−9, p < 10−4, respectively), whereas the measurements of ARA4, DHA, and EPA were highly correlated (ARA4 vs. DHA: r = 0.916, p < 10−76; ARA4 vs. EPA: r = 0.705, p < 10−29; DHA V.S. EPA: r = 0.684, p < 10−27).
Table 1. Summary of fat, nine classified fatty acids, and 12 single fatty acids.
Figure 1. Phenotypic distribution and correlation heatmap. (A) Box plot of phenotype distribution. The X-axis represents phenotypes, and the Y-axis represents the proportion in the total fatty acids. (B) Correlation of all phenotypes. The colors represent the correlation coefficient R2 values: red means positive correlation and blue means negative correlation.
Genome-Wide Association Analysis and Gene Annotation of Identified SNPs
We identified nine SNPs (7, 8, 2, and 1 SNP for ARA4, DHA, EPA, and ALA, respectively) that achieved the suggestive significance threshold (P < 3.445 × 10−5), in which 4 SNPs (4, 3, 0, and 0 SNPs for ARA4, DHA, EPA, and ALA, respectively) surpassed the significance line (P < 1.72 × 10−6) for the content of four fatty acids. To further identify more potentially associated SNPs, we further investigated the results of the DHA_EPA group, which contained two associated SNPs and six suggestive SNPs (Figure 2 and Supplementary Table S2). The Manhattan plots and Q-Q plots are shown in Figure 2 for traits with associated SNPs (ARA4, DHA, and DHA_EPA), in which the Q-Q plots (genomic inflation factor ≈1) indicated reliability of experimental design and data analysis. Genes within the 10-kb regions of the associated and suggestive SNPs were annotated using the common carp genome (Xu P. et al., 2014), and a total of 15 genes were identified (six genes for associated SNPs, 12 genes for suggestive SNPs).
Figure 2. Manhattan plot and Q-Q plot for traits with significant SNPs. (A) Manhattan plot for the GWAS results of muscular ARA4 content. The red line indicates the significant threshold while the blue line shows the suggestive threshold. (B) Q-Q plot for GWAS results of muscular ARA4 content. (C) Manhattan plot for the GWAS results of muscular DHA content. (D) Q-Q plot for the GWAS results of muscular DHA content. (E) Manhattan plot for the GWAS results of muscular DHA + EPA content. (F) Q-Q plot for GWAS results of muscular DHA + EPA content.
Transcriptomic Profiling of Samples With Divergent PUFA Contents
We assessed the fatty acid content of 20 new individuals and selected eight showing extreme PUFA content (H: high content group; L: low content group). The ARA4, EPA, DHA, and ALA content of the 20 individuals are shown in Supplementary Table S3. Three types of tissues (brain, liver, and muscle) were collected from eight individuals, and we used Illumina high-throughput sequencing to generate mRNA transcriptomes. The overall quality of the extracted RNA from 24 tissue samples is summarized in Supplementary Table S4. After sequencing, approximately 176 Gb of raw data were generated from 24 tissue samples, and 173.1 Gb clean data were used for further differential gene expression (DGE) analysis (Supplementary Table S5). Using the Bowtie software, 76.1% of the clean reads could be mapped to the genome, and the DGE of each trait between the high-content and low-content groups was calculated using the Cufflinks software.
As shown in the volcano plots in Figure 3, a large number of genes were identified to be differentially expressed in the brain, with counts of 6,191 genes, 916 genes, 1,010 genes, and 36 genes for ARA4, DHA, EPA, and ALA, respectively (Figure 3 and Supplementary Table S6). We further conducted GO and KEGG enrichment analyses to focus on the most important pathways (Figure 3 and Supplementary Tables S6–S8). Several important lipid and fatty acid metabolism pathways were enriched, including “adipocytokine signaling pathway,” “ARA4 metabolism,” and “glycerolipid metabolism.” Key genes within these enriched pathways included cd36, npy, and lepr. Interestingly, several development and metabolism-related pathways were also enriched (Figures 3C,D) such as “Wnt signaling pathway,” “Notch signaling pathway,” and the “insulin signaling pathway,” indicating possible relationship between growth and muscular fatty acid. Pathways of amino acid metabolism (especially branched chain amino acids) were also enriched such as “valine, leucine, and isoleucine degradation,” indicating related pathways between amino acid and fatty acid metabolism processes.
Figure 3. Volcano plot, Venn diagram, Gene Ontology, and KEGG enrichment for differentially expressed genes in brain tissues. (A) Volcano plots for DEGs in the brain tissues for four traits. Red dots indicate upregulated genes; blue dots indicate downregulated genes. (B) Venn diagram of DEGs in brain tissues showing the number of shared and unique genes for each trait. (C) Gene Ontology enrichment of DEGs in brain tissues. The size of the circles represents gene numbers in each term; colors represent minus logarithms of adjusted P values. (D) KEGG enrichment of DEGs in brain tissue. Size of circles represents gene numbers in each pathway; colors represent minus logarithms of adjusted P values.
In liver tissues, 113, 299, 228, and 60 genes were identified to be differentially regulated in the four groups, respectively (Supplementary Figure S1 and Supplementary Tables S9–S11). Similar to the results of the brain tissues, pathways related to growth and amino acid metabolism were also enriched in the liver tissues such as the “insulin signaling pathway,” “glycine, serine, and threonine metabolism,” and the “Wnt signaling pathway,” including key genes such as prkaa2 and gsk-3β. In the muscles, there were 23, 15, 38, and 4 genes that were differentially expressed in the four groups, respectively (Supplementary Figure S2 and Supplementary Table S12). Volcano plots and Venn diagrams showed a limited number of genes that were differentially expressed in the muscle tissues, thus no GO terms or KEGG pathways were enriched. Supplementary Figure S3 shows the results of cluster analysis of the four traits in three types of tissues.
Patterns of Methylation Variations Among Groups With Distinct PUFA Contents
To assess the global trends of epigenetic variations in different groups, we performed genome-wide pairwise comparisons of each epigenetic modification between the H and L groups. Muscle tissues were collected from the eight samples used in transcriptome analysis, and DNA was extracted for further WGBS. A total of 541.3 Gb of raw data (3,608,912,206 raw reads) were acquired (around 45× coverage for each sample), and 535.1 Gb of clean data were used for alignment to the carp genome (Figure 4A and Supplementary Table S13). Following standard pipelines (see Methods), DMRs were identified (Figure 4B). The most abundant DMRs were located in intergenic repeat regions, followed by introns, exons, and promoter regions (Supplementary Table S14). Because DMRs in the promoter regions (differentially methylated promoters, DMPs) play vital roles in the regulation of transcription (Haque et al., 2016), we selected DMPs for further functional enrichment analysis (Supplementary Tables S15, S16). GO enrichment analysis identified various terms, including “embryonic organ development,” “embryonic organ morphogenesis,” which were related to growth and development (Figure 4C). More direct evidences have been identified in KEGG enrichment such as shared pathways (Figure 4D), namely, “insulin signaling pathway,” “glycerophospholipid metabolism,” “ether lipid metabolism” for the four traits, and fatty acid metabolism-related pathways, including “ARA4 metabolism,” “linoleic acid metabolism,” “biosynthesis of unsaturated fatty acids,” and “adipocytokine signaling pathway.” Furthermore, “steroid hormone biosynthesis” was enriched, indicating the role of fatty acids in hormone metabolism.
Figure 4. Sequencing depth distribution, DMR classification, and enrichment of genes within DMP regions. (A) Accumulative fraction against sequencing depth. Curve lines with different colors represent eight samples used for WGBS, and the Y-axis represents the cumulative ratio of genomic regions mapped by reads. (B) Counts for each DMR classification in different genomic regions for each trait. (C) Gene Ontology enrichment of genes downstream of DMPs. The size of the circles represents the number of genes in each term; colors represent minus logarithms of adjusted P values. (D) KEGG enrichment of genes downstream of DMPs. The size of the circles represents the number of genes in each pathway; colors represent minus logarithms of adjusted P values.
Integrated Analysis of DGE and DMP Results
Although DGE and DMP analyses identified several genes and pathways, the correlation between the two regulatory levels remains unclear. Methylation levels in the promoter regions significantly affect the transcription of downstream genes; therefore, it is necessary to conduct a correlation analysis of genes shared between the DGE and DMP results. Due to the limited number of genes shared between the liver/muscle DGE and DMP results, we focused on the correlation between brain DGE log2(FC) values and DMP gene methylation differences (Supplementary Table S17). Figure 5 shows the results of Pearson correlation analysis, which indicated a significant linear correlation for the ARA4, DHA, and EPA traits. The genes shown in Figure 5 were classified into two categories, namely, positively correlated (blue dots) and negatively correlated (red dots). For the ARA4, DHA, and EPA traits, 457, 35, and 59 genes were included in the correlation analysis. The R2 values of the 221 positively correlated genes and 236 negatively correlated genes in ARA4 trait analysis were 0.68 and 0.55, respectively. For the DHA trait, the R2 values of 21 positively correlated genes and 14 negatively correlated genes were 0.44 and 0.60, respectively. For the EPA trait, R2 values of 46 positively correlated genes and 13 negatively correlated genes were 0.51 and 0.77, respectively. The DGE and DMP results of the liver and muscle tissues were also integrated (Supplementary Table S17), which identified fewer genes compared to those in the brain. To find more evidence supporting the correlation between gene expression and genome methylation, TFBS prediction of the promoters of selected genes from integrated analyses was performed (Supplementary Table S18).
Figure 5. Pearson correlation of gene expression fold changes and differential methylation rates for DEGs in brain tissues and enrichment analysis. Scatter plots for (A) ARA4, (B) DHA, and (C) EPA. The X-axis represents mean difference in methylation ratio between high and low ARA4 content groups; the Y-axis represents the logarithm of fold changes (log2FC). The red dots represent negatively correlated genes and blue dots represent positively correlated genes. The trend lines and formula in each scatter plot represent the correlation coefficients. (D) GO and KEGG enrichment of correlated genes for three traits. The size of circles represents the number of genes in each pathway; colors represent minus logarithms of adjusted P values.
Discussion
GWAS of 203 Individuals and Gene Annotation
Among the identified SNPs in GWAS, the most significant SNP carp222748 has been associated with ARA4 content in C. carpio, which is located within the coding region of the duox2 gene. The Duox2 (dual oxidase 2) protein, which is encoded by duox2, generates hydrogen peroxide, which is required for the activity of thyroid peroxidase and plays a role in thyroid hormone synthesis. Interestingly, another gene, trh (pro-thyrotropin-releasing hormone), was identified downstream of a suggestive SNP carp215797, and the Trh protein stimulates the release of thyrotropin. Thyroid hormone status affects metabolic pathways of ARA4 in mice and human (Fonseca et al., 2014; Yao et al., 2015), and thus it is possible that a similar genetic mechanism may be utilized in ARA4 metabolism in C. carpio. Another associated SNP carp152877 was found to be associated with the ARA4, DHA, and DHA_EPA traits and is located downstream of the ormdl3 gene. The Ormdl3 (ORM1-like protein 3) protein negatively regulates sphingolipid synthesis and may be indirectly involved in endoplasmic reticulum-mediated calcium ion signaling (Breslow et al., 2010). As the genes discovered from GWAS analysis may not be enough to illustrate the underlying mechanism of muscular PUFA metabolism, we sought to uncover more evidences by investigating gene expression and epigenetic statuses of the collected samples.
GO and KEGG Analysis Based on PUFA RNA-Seq Data
GO and KEGG enrichment analyses of the DGE results identified a handful of genes that are involved in enriched pathways that are related to fatty acid metabolism. CD36 has been reported as a high affinity receptor for long-chain fatty acid (FA) uptake, in addition to the contribution of lipid accumulation and FA-initiated signaling (Pepino et al., 2014). Neuropeptide Y (NPY), an orexigenic hypothalamic neuropeptide that is released by the neurons of arcuate nuclei, influences foraging behavior and food intake in mammals and indirectly affects fat accumulation (Levin et al., 2013). Leptin is an adipocytokine that regulates energy intake and expenditure through interactions with the leptin receptor (LEPR). LEPR has been associated with intramuscular fat and FA accumulation in Duroc pigs (Ros-Freixedes et al., 2016). We presume these genes may have similar functions in C. carpio. Insulin has been reported to promote development and fatty acid biosynthesis through the conversion of glucose into triglyceride in the liver, fat, and muscle cells (Laron, 2008; Xu X. et al., 2012; Pramfalk et al., 2016).
Among the potentially associated genes identified in liver tissue pathways, protein kinase α2 (PRKAA2) participates in lipid metabolism and energy homeostasis. A recent study on mice showed that lipopolysaccharide could significantly inhibit PRKAA2 expression (Tzanavari et al., 2016). In the Wnt and insulin signaling pathways, gsk-3β has been enriched in the liver tissues for DHA and EPA. GSK-3β acts as an upstream regulator of the ACSL family and lipid accumulation in hepatocytes (Chang et al., 2011). Previous studies have shown that high levels of amino acids upregulate hepatic fatty acid biosynthetic gene expression in trout hepatocytes (Dai et al., 2015). Despite intragroup variations in each group due to limited samples, distinct expression patterns were observed between the H (high) and L (low) groups. Interestingly, several genes identified in GWAS were also differentially expressed in three tissues, including ormdl3, trh, and nptn. This provided potential correlations between SNP frequency and relative gene expression. Despite the sample size of GWAS and RNA-Seq were both too small to make a comprehensive exploration of all trait-related genes, the results were still very informative for further larger sample analysis.
Predicting Fatty Acid Metabolism-Related Genes With DNA Methylation Data
The core genes within DMP-enriched pathways included gpd1, acsl1, fads2, elovl6, and the pla2 family (pla2g15, pla2g12b, and pla2g6), some of which have already been reported to be involved in PUFA metabolism in various species. For example, acyl-coenzyme A (CoA) syntheses 1 (ACSL1) is a well-studied obesogenic gene that is involved in fatty acid metabolism. The expression of ACSL1 has been reported to be associated with high caloric food intake in mice and human (Joseph et al., 2015). In the biosynthesis of PUFA, Δ6 desaturase (FADS2) has been demonstrated as an important indicator that catalyzes the first denaturation step and influences PUFA synthesis capability in fish. The expression of fads2 is regulated by dietary fatty acid profiles in Japanese seabass and is significantly negatively correlated with CpG methylation rates in the fads2 gene promoter (Xu H. et al., 2014). The elongation of very long chain fatty acid 6 (ELOVL6) mainly catalyzes the elongation of long chain SFAs and MUFAs (C14, C16, and C18), while it is engaged in balancing the overall fatty acid composition in mammals (Corominas et al., 2015).
Integrated Analysis of DGE and DMP Results
Taken together, the significant linear correlation between expression and methylation indicates the epigenetic effects on transcription in addition to genetic factors such as genome variations. GO and KEGG enrichment further confirmed the importance of these correlated genes, and vital pathways were enriched, including “adipocytokine signaling pathway,” “glycerophospholipid metabolism,” “Notch signaling pathway,” “Wnt signaling pathway,” “ion binding,” and “regulation of growth”. Genes enriching these pathways included agrp, ctbp2, and previously discussed acsl1, gpd1, and pla2g15. AGRP has been reported to play similar roles in food intake and high-fat diet preference as NPY in mice and bats, and the upregulation of the agrp gene enhances diet preference and fat gain (Levin et al., 2013). The ctbp2 gene was downregulated and highly methylated in the high-PUFA content group, which agrees with the findings of previous studies that overexpression of C-terminal-binding protein 2 (CTBP2) suppresses lipid accumulation and hepatic glucose uptake (Liu et al., 2017). Fewer correlated genes were identified from the liver and muscle data, partly because of the limited number of samples. The acsl5 gene showed higher expression levels and lower methylation levels in the promoter region in the H group compared to the L group in relation to muscular ALA content. ACSL5 has been reported to act like other ACSL family members and engages in lipid biosynthesis and fatty acid degradation (Bowman et al., 2016). Besides, Fad and Elovl enzymes encoded by fad and elovl genes have been demonstrated to have all the activities and jointly play an important role for the biosynthesis of DHA from C18 PUFA in rabbitfish, as well as in other teleost (Monroig et al., 2012). In mammal and teleost fish, Elovl2, Elovl4, and Elovl5 have been cloned and functionally characterized separately as crucial elongation enzymes in PUFA biosynthesis (Agaba et al., 2005; Jakobsson et al., 2006; Morais et al., 2009; Monroig et al., 2010). It has been reported that the capability of PUFA biosynthesis varies among teleost fish with alternative pathways during evolution (Castro et al., 2016). Therefore, we assume elovl6 which has been enriched in DMP pathway has similar functions in its encoding Eolvl6 enzyme. Fads2 with Δ4, Δ6, and Δ8 desaturase activities, respectively, composed a crucial part of PUFA biosynthesis as well (Li et al., 2010; Oboh et al., 2017). After integration of three types of analyses, the important genes identified in this study were summarized in Table 2, which provided a full view of genes supported by multiple evidences. We used muscles for WGBS and methylation analysis and considered that it could represent genome methylation of an individual. However, we could not ignore the slight differences in methylation among tissues. Investigations on tissue-specific methylation sites using multi-omics analyses and the correlation between DGE and DMP in the same tissue are thus warranted. Our results indicate that fatty acid metabolism-related genes are associated with the growth, hormone, and even amino acid relevant genes, which regulate muscular PUFA content.
Table 2. Summary of genes identified in three types of analyses.
Conclusion
This study investigated the associated genes and the divergence of transcriptomic and epigenomic variations in C. carpio samples with distinct muscular PUFA contents. The phenotypic correlation and GWAS results indicated dozens of potential genes that are associated with PUFA metabolism. Further investigation of DGE patterns of three tissues and gene profiles with DMPs identified important genes that enriched PUFA synthesis- and degradation-related pathways. Integrated analysis showed significant correlations between DGE and DMP and enrichment of correlated genes involved in vital pathways related to fatty acid desaturation, elongation, and transportation. Confirmation of these results using integrated genomic, transcriptomic, and epigenomic profiling with more extensive sequencing of larger samples is warranted. This study may also facilitate the applications of combined genome selection and molecular-associated breeding for multiple traits.
Data Availability
The sequencing datasets of all samples have been deposited at NCBI (PRJNA493161).
Author Contributions
JX initiated and coordinated the research project. HZ and JX conceived and conducted the analysis and drafted the manuscript. PX engaged in sample collection and genotyping analysis. YJ and ZZ assisted in transcriptome and methylation data analysis. JF, RT, and CD took part in trait measurement, tissue manipulation, and enrichment analysis. All the authors read and approved the final manuscript.
Funding
The Central Public-interest Scientific Institution Basal Research Fund, CAFS (Nos. 2016HY-ZD0302, 2016GH02, and 2016HY-JC0301), the National Natural Science Foundation of China (Nos. 31502151 and 31422057), the National High-Technology Research and Development Program of China (2011AA100401), and the National Infrastructure of Fishery Germplasm Resources of China (No. 2018DKA30470) supported this study.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2019.00217/full#supplementary-material
Footnotes
- ^ http://www.fao.org
- ^ https://cran.r-project.org/web/packages/ggplot2/
- ^ https://www.cog-genomics.org/plink2
- ^ http://www.maizegenetics.net/tassel
- ^ http://cran.r-project.org/package=qqman
- ^ https://david.ncifcrf.gov/
- ^ http://sourceforge.net/projects/swdmr/
- ^ http://bioinfo.life.hust.edu.cn/AnimalTFDB/#!/tfbs_predict
References
Agaba, M., Tocher, D. R., Dickson, C. A., Dick, J. R., and Teale, A. J. (2004). Zebrafish cDNA encoding multifunctional fatty acid elongase involved in production of eicosapentaenoic (20 : 5n-3) and docosahexaenoic (22 : 6n-3) acids. Mar. Biotechnol. 6, 251–261. doi: 10.1007/s10126-003-0029-1
Agaba, M. K., Tocher, D. R., Zheng, X. Z., Dickson, C. A., Dick, J. R., and Teale, A. J. (2005). Cloning and functional characterisation of polyunsaturated fatty acid elongases of marine and freshwater teleost fish. Comp. Biochem. Physiol. B Biochem. Mol. Biol. 142, 342–352. doi: 10.1016/j.cbpb.2005.08.005
Ameur, A., Enroth, S., Johansson, A., Zaboli, G., Igl, W., Johansson, A. C. V., et al. (2012). Genetic adaptation of fatty-acid metabolism: a human-specific haplotype increasing the biosynthesis of long-chain Omega-3 and Omega-6 fatty acids. Am. J. Hum. Genet. 90, 809–820. doi: 10.1016/j.ajhg.2012.03.014
Beld, J., Lee, D. J., and Burkart, M. D. (2015). Fatty acid biosynthesis revisited: structure elucidation and metabolic engineering. Mol. Biosyst. 11, 38–59. doi: 10.1039/c4mb00443d
Berton, M. P., Fonseca, L. F., Gimenez, D. F., Utembergue, B. L., Cesar, A. S., Coutinho, L. L., et al. (2016). Gene expression profile of intramuscular muscle in Nellore cattle with extreme values of fatty acid. BMC Genomics 17:972. doi: 10.1186/s12864-016-3232-y
Bolger, A. M., Lohse, M., and Usadel, B. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120. doi: 10.1093/bioinformatics/btu170
Bostock, J., Mcandrew, B., Richards, R., Jauncey, K., Telfer, T., Lorenzen, K., et al. (2010). Aquaculture: global status and trends. Philos. Trans. R. Soc. Lond. B Biol. Sci. 365, 2897–2912. doi: 10.1098/rstb.2010.0170
Bowman, T. A., O’keeffe, K. R., D’aquila, T., Yan, Q. W., Griffin, J. D., Killion, E. A., et al. (2016). Acyl CoA synthetase 5 (ACSL5) ablation in mice increases energy expenditure and insulin sensitivity and delays fat absorption. Mol. Metab. 5, 210–220. doi: 10.1016/j.molmet.2016.01.001
Bradbury, P. J., Zhang, Z., Kroon, D. E., Casstevens, T. M., Ramdoss, Y., and Buckler, E. S. (2007). TASSEL: software for association mapping of complex traits in diverse samples. Bioinformatics 23, 2633–2635. doi: 10.1093/bioinformatics/btm308
Breslow, D. K., Collins, S. R., Bodenmiller, B., Aebersold, R., Simons, K., Shevchenko, A., et al. (2010). Orm family proteins mediate sphingolipid homeostasis. Nature 463, 1048–1053. doi: 10.1038/nature08787
Calder, P. C. (2013). n-3 fatty acids, inflammation and immunity: new mechanisms to explain old actions. Proc. Nutr. Soc. 72, 326–336. doi: 10.1017/S0029665113001031
Castro, L. F., Tocher, D. R., and Monroig, O. (2016). Long-chain polyunsaturated fatty acid biosynthesis in chordates: insights into the evolution of Fads and Elovl gene repertoire. Prog. Lipid Res. 62, 25–40. doi: 10.1016/j.plipres.2016.01.001
Chang, C. C., Chow, C. C., Tellier, L. C., Vattikuti, S., Purcell, S. M., and Lee, J. J. (2015). Second-generation PLINK: rising to the challenge of larger and richer datasets. Gigascience 4:7. doi: 10.1186/s13742-015-0047-8
Chang, Y. S., Tsai, C. T., Huangfu, C. A., Huang, W. Y., Lei, H. Y., Lin, C. F., et al. (2011). ACSL3 and GSK-3beta are essential for lipid upregulation induced by endoplasmic reticulum stress in liver cells. J. Cell. Biochem. 112, 881–893. doi: 10.1002/jcb.22996
Chen, S., Chai, M., Tian, C., Li, Y., Deng, T., Wu, H., et al. (2018). Genetic variants of fatty acid elongase 6 in Chinese Holstein cow. Gene 670, 123–129. doi: 10.1016/j.gene.2018.05.073
Corominas, J., Marchesi, J. A., Puig-Oliveras, A., Revilla, M., Estelle, J., Alves, E., et al. (2015). Epigenetic regulation of the ELOVL6 gene is associated with a major QTL effect on fatty acid composition in pigs. Genet. Sel. Evol. 47:20. doi: 10.1186/s12711-015-0111-y
Dai, W., Panserat, S., Plagnes-Juan, E., Seiliez, I., and Skiba-Cassy, S. (2015). Amino acids attenuate insulin action on gluconeogenesis and promote fatty acid biosynthesis via mTORC1 signaling pathway in trout hepatocytes. Cell Physiol. Biochem. 36, 1084–1100. doi: 10.1159/000430281
Davidson, M. H. (2013). Omega-3 fatty acids: new insights into the pharmacology and biology of docosahexaenoic acid, docosapentaenoic acid, and eicosapentaenoic acid. Curr. Opin. Lipidol. 24, 467–474. doi: 10.1097/MOL.0000000000000019
Davis-Bruno, K., and Tassinari, M. S. (2011). Essential fatty acid supplementation of DHA and ARA and effects on neurodevelopment across animal species: a review of the literature. Birth Defects Res. B Dev. Reprod. Toxicol. 92, 240–250. doi: 10.1002/bdrb.20311
Derayat, A., Houston, R. D., Guy, D. R., Hamilton, A., Ralph, J., Spreckley, N., et al. (2007). Mapping QTL affecting body lipid percentage in Atlantic salmon (Salmo salar). Aquaculture 272, S250–S251. doi: 10.1016/j.aquaculture.2007.07.046
Enns, J. E., Yeganeh, A., Zarychanski, R., Abou-Setta, A. M., Friesen, C., Zahradka, P., et al. (2014). The impact of omega-3 polyunsaturated fatty acid supplementation on the incidence of cardiovascular events and complications in peripheral arterial disease: a systematic review and meta-analysis. BMC Cardiovasc. Disord. 14:70. doi: 10.1186/1471-2261-14-70
Fonseca, T. L., Werneck-De-Castro, J. P., Castillo, M., Bocco, B. M., Fernandes, G. W., Mcaninch, E. A., et al. (2014). Tissue-specific inactivation of type 2 deiodinase reveals multilevel control of fatty acid oxidation by thyroid hormone in the mouse. Diabetes 63, 1594–1604. doi: 10.2337/db13-1768
Foster, D. W. (2012). Malonyl-CoA: the regulator of fatty acid synthesis and oxidation. J. Clin. Invest. 122, 1958–1959. doi: 10.1172/JCI63967
Geay, F., Wenon, D., Mellery, J., Tinti, E., Mandiki, S. N. M., Tocher, D. R., et al. (2015). Dietary linseed oil reduces growth while differentially impacting LC-PUFA synthesis and accretion into tissues in eurasian perch (Perca fluviatilis). Lipids 50, 1219–1232. doi: 10.1007/s11745-015-4079-8
Ginestet, C. (2011). ggplot2: elegant graphics for data analysis. J. R. Stat. Soc. Ser. A Stat. Soc. 174, 245–245. doi: 10.1007/978-0-387-98141-3
Guesnet, P., and Alessandri, J. M. (2011). Docosahexaenoic acid (DHA) and the developing central nervous system (CNS) - implications for dietary recommendations. Biochimie 93, 7–12. doi: 10.1016/j.biochi.2010.05.005
Haque, M. M., Nilsson, E. E., Holder, L. B., and Skinner, M. K. (2016). Genomic Clustering of differential DNA methylated regions (epimutations) associated with the epigenetic transgenerational inheritance of disease and phenotypic variation. BMC Genomics 17:418. doi: 10.1186/s12864-016-2748-5
He, G., Chen, B., Wang, X., Li, X., Li, J., He, H., et al. (2013). Conservation and divergence of transcriptomic and epigenomic variation in maize hybrids. Genome Biol. 14:R57. doi: 10.1186/gb-2013-14-6-r57
Henriques, J., Dick, J. R., Tocher, D. R., and Bell, J. G. (2014). Nutritional quality of salmon products available from major retailers in the UK: content and composition of n-3 long-chain PUFA. Br. J. Nutr. 112, 964–975. doi: 10.1017/S0007114514001603
Hu, H., Miao, Y. R., Jia, L. H., Yu, Q. Y., Zhang, Q., and Guo, A. Y. (2018). AnimalTFDB 3.0: a comprehensive resource for annotation and prediction of animal transcription factors. Nucleic Acids Res. 47, D33–D38. doi: 10.1093/nar/gky822
Huang, D. W., Sherman, B. T., and Lempicki, R. A. (2009). Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat. Protoc. 4, 44–57. doi: 10.1038/nprot.2008.211
Innis, S. M. (2008). Dietary omega 3 fatty acids and the developing brain. Brain Res. 1237, 35–43. doi: 10.1016/j.brainres.2008.08.078
Iverson, S. J., Frost, K. J., and Lang, S. L. C. (2002). Fat content and fatty acid composition of forage fish and invertebrates in Prince William Sound, Alaska: factors contributing to among and within species variability. Mar. Ecol. Prog. Ser. 241, 161–181. doi: 10.3354/meps241161
Jakobsson, A., Westerberg, R., and Jacobsson, A. (2006). Fatty acid elongases in mammals: their regulation and roles in metabolism. Prog. Lipid Res. 45, 237–249. doi: 10.1016/j.plipres.2006.01.004
Ji, P. F., Liu, G. M., Xu, J., Wang, X. M., Li, J. T., Zhao, Z. X., et al. (2012). Characterization of common carp transcriptome: sequencing, de novo assembly, annotation and comparative genomics. PLoS One 7:e35152. doi: 10.1371/journal.pone.0035152
Jiang, Y. L., Zhang, S. H., Xu, J., Feng, J. X., Mahboob, S., Al-Ghanim, K. A., et al. (2014). Comparative transcriptome analysis reveals the genetic basis of skin color variation in common carp. PLoS One 9:e108200. doi: 10.1371/journal.pone.0108200
Joseph, R., Poschmann, J., Sukarieh, R., Too, P. G., Julien, S. G., Xu, F., et al. (2015). ACSL1 is associated with fetal programming of insulin sensitivity and cellular lipid content. Mol. Endocrinol. 29, 909–920. doi: 10.1210/me.2015-1020
Kim, D., Pertea, G., Trapnell, C., Pimentel, H., Kelley, R., and Salzberg, S. L. (2013). TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 14:R36. doi: 10.1186/gb-2013-14-4-r36
Kjaer, M. A., Ruyter, B., Berge, G. M., Sun, Y., and Ostbye, T. K. (2016). Regulation of the Omega-3 fatty acid biosynthetic pathway in atlantic salmon hepatocytes. PLoS One 11:e0168230. doi: 10.1371/journal.pone.0168230
Krueger, F., and Andrews, S. R. (2011). Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics 27, 1571–1572. doi: 10.1093/bioinformatics/btr167
Kuang, Y. Y., Zheng, X. H., Lv, W. H., Cao, D. C., and Sun, X. W. (2015). Mapping quantitative trait loci for flesh fat content in common carp (Cyprinus carpio). Aquaculture 435, 100–105. doi: 10.1016/j.aquaculture.2014.09.020
Langmead, B., and Salzberg, S. L. (2012). Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359. doi: 10.1038/nmeth.1923
Laron, Z. (2008). Insulin–a growth hormone. Arch. Physiol. Biochem. 114, 11–16. doi: 10.1080/13813450801928356
Levin, E., Yom-Tov, Y., Hefetz, A., and Kronfeld-Schor, N. (2013). Changes in diet, body mass and fatty acid composition during pre-hibernation in a subtropical bat in relation to NPY and AgRP expression. J. Comp. Physiol. B 183, 157–166. doi: 10.1007/s00360-012-0689-0
Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., et al. (2009). The sequence alignment/map format and SAMtools. Bioinformatics 25, 2078–2079. doi: 10.1093/bioinformatics/btp352
Li, Y., Monroig, O., Zhang, L., Wang, S., Zheng, X., Dick, J. R., et al. (2010). Vertebrate fatty acyl desaturase with Delta4 activity. Proc. Natl. Acad. Sci. U.S.A. 107, 16840–16845. doi: 10.1073/pnas.1008429107
Liu, P. L., Shi, L., Cang, X. M., Huang, J. R., Wu, X., Yan, J., et al. (2017). CtBP2 ameliorates palmitate-induced insulin resistance in HepG2 cells through ROS mediated JNK pathway. Gen. Comp. Endocrinol. 247, 66–73. doi: 10.1016/j.ygcen.2017.01.018
Monroig, O., Rotllant, J., Cerda-Reverter, J. M., Dick, J. R., Figueras, A., and Tocher, D. R. (2010). Expression and role of Elovl4 elongases in biosynthesis of very long-chain fatty acids during zebrafish Danio rerio early embryonic development. Biochim. Biophys. Acta 1801, 1145–1154. doi: 10.1016/j.bbalip.2010.06.005
Monroig, O., Wang, S. Q., Zhang, L., You, C. H., Tocher, D. R., and Li, Y. Y. (2012). Elongation of long-chain fatty acids in rabbitfish Siganus canaliculatus: cloning, functional characterisation and tissue distribution of Elovl5-and Elovl4-like elongases. Aquaculture 350, 63–70. doi: 10.1016/j.aquaculture.2012.04.017
Morais, S., Monroig, O., Zheng, X. Z., Leaver, M. J., and Tocher, D. R. (2009). Highly unsaturated fatty acid synthesis in atlantic salmon: characterization of ELOVL5-and ELOVL2-like elongases. Mar. Biotechnol. 11, 627–639. doi: 10.1007/s10126-009-9179-0
Mozaffarian, D., and Wu, J. H. (2011). Omega-3 fatty acids and cardiovascular disease: effects on risk factors, molecular pathways, and clinical events. J. Am. Coll. Cardiol. 58, 2047–2067. doi: 10.1016/j.jacc.2011.06.063
Murphy, R. C. (2013). Lipid Biochemistry, Mass Spectrometry, and Lipidomics: A Journey in Five Decades. Washington, DC: American Chemical Society, 245.
Oboh, A., Kabeya, N., Carmona-Antonanzas, G., Castro, L. F. C., Dick, J. R., Tocher, D. R., et al. (2017). Two alternative pathways for docosahexaenoic acid (DHA, 22:6n-3) biosynthesis are widespread among teleost fish. Sci. Rep. 7:3889. doi: 10.1038/s41598-017-04288-2
Pan, A., Chen, M., Chowdhury, R., Wu, J. H. Y., Sun, Q., Campos, H., et al. (2012). alpha-Linolenic acid and risk of cardiovascular disease: a systematic review and meta-analysis. Am. J. Clin. Nutr. 96, 1262–1273. doi: 10.3945/ajcn.112.044040
Peng, W., Xu, J., Zhang, Y., Feng, J., Dong, C., Jiang, L., et al. (2016). An ultra-high density linkage map and QTL mapping for sex and growth-related traits of common carp (Cyprinus carpio). Sci. Rep. 6:26693. doi: 10.1038/srep26693
Pepino, M. Y., Kuda, O., Samovski, D., and Abumrad, N. A. (2014). Structure-function of CD36 and importance of fatty acid signal transduction in fat metabolism. Annu. Rev. Nutr. 34, 281–303. doi: 10.1146/annurev-nutr-071812-161220
Pramfalk, C., Pavlides, M., Banerjee, R., Mcneil, C. A., Neubauer, S., Karpe, F., et al. (2016). Fasting plasma insulin concentrations are associated with changes in hepatic fatty acid synthesis and partitioning prior to changes in liver fat content in healthy adults. Diabetes 65, 1858–1867. doi: 10.2337/db16-0236
Ralston, J. C., Matravadia, S., Gaudio, N., Holloway, G. P., and Mutch, D. M. (2015). Polyunsaturated fatty acid regulation of adipocyte FADS1 and FADS2 expression and function. Obesity 23, 725–728. doi: 10.1002/oby.21035
Ros-Freixedes, R., Gol, S., Pena, R. N., Tor, M., Ibanez-Escriche, N., Dekkers, J. C., et al. (2016). Genome-wide association study singles out SCD and LEPR as the two main loci influencing intramuscular fat content and fatty acid composition in Duroc Pigs. PLoS One 11:e0152496. doi: 10.1371/journal.pone.0152496
Rymer, C., and Givens, D. I. (2005). n-3 fatty acid enrichment of edible tissue of poultry: a review. Lipids 40, 121–130. doi: 10.1007/s11745-005-1366-4
Sodeland, M., Gaarder, M., Moen, T., Thomassen, M., Kjoglum, S., Kent, M., et al. (2013). Genome-wide association testing reveals quantitative trait loci for fillet texture and fat content in Atlantic salmon. Aquaculture 408, 169–174. doi: 10.1016/j.aquaculture.2013.05.029
Trapnell, C., Roberts, A., Goff, L., Pertea, G., Kim, D., Kelley, D. R., et al. (2014). Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks (vol 7, pg 562, 2012). Nat. Protoc. 9, 2513–2513. doi: 10.1038/nprot.2012.016
Trapnell, C., Williams, B. A., Pertea, G., Mortazavi, A., Kwan, G., Van Baren, M. J., et al. (2010). Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat. Biotechnol. 28, 511–515. doi: 10.1038/nbt.1621
Tzanavari, T., Varela, A., Theocharis, S., Ninou, E., Kapelouzou, A., Cokkinos, D. V., et al. (2016). Metformin protects against infection-induced myocardial dysfunction. Metabolism 65, 1447–1458. doi: 10.1016/j.metabol.2016.06.012
Vedtofte, M. S., Jakobsen, M. U., Lauritzen, L., and Heitmann, B. L. (2011). Dietary alpha-linolenic acid, linoleic acid, and n-3 long-chain PUFA and risk of ischemic heart disease. Am. J. Clin. Nutr. 94, 1097–1103. doi: 10.3945/ajcn.111.018762
Wang, C., Wachholtz, M., Wang, J., Liao, X. L., and Lu, G. Q. (2014a). Analysis of the skin transcriptome in two oujiang color varieties of common carp. PLoS One 9:e90074. doi: 10.1371/journal.pone.0090074
Wang, C., Zhang, Z., Yao, H., Zhao, F., Wang, L., Wang, X., et al. (2014b). Effects of atrazine and chlorpyrifos on DNA methylation in the liver, kidney and gill of the common carp (Cyprinus carpio L.). Ecotoxicol. Environ. Saf. 108, 142–151. doi: 10.1016/j.ecoenv.2014.06.011
Wang, Y., Zheng, W. L., Luo, J. F., Zhang, D. D., and Lu, Z. H. (2006). In situ bisulfite modification of membrane-immobilized DNA for multiple methylation analysis. Anal. Biochem. 359, 183–188. doi: 10.1016/j.ab.2006.09.001
Wang, Z., Li, X. F., Jiang, Y., Shao, Q. Z., Liu, Q., Chen, B. Y., et al. (2015). swDMR: a sliding window approach to identify differentially methylated regions based on whole genome bisulfite sequencing. PLoS One 10:e0132866. doi: 10.1371/journal.pone.0132866
Witte, T. R., and Hardman, W. E. (2015). The effects of omega-3 polyunsaturated fatty acid consumption on mammary carcinogenesis. Lipids 50, 437–446. doi: 10.1007/s11745-015-4011-2
Wood, J. D., Enser, M., Fisher, A. V., Nute, G. R., Sheard, P. R., Richardson, R. I., et al. (2008). Fat deposition, fatty acid composition and meat quality: a review. Meat Sci. 78, 343–358. doi: 10.1016/j.meatsci.2007.07.019
Xia, J. H., Lin, G., He, X. P., Yunping, B., Liu, P., Liu, F., et al. (2014). Mapping quantitative trait loci for Omega-3 fatty acids in Asian Seabass. Mar. Biotechnol. 16, 1–9. doi: 10.1007/s10126-013-9524-1
Xiao, S. J., Wang, P. P., Dong, L. S., Zhang, Y. G., Han, Z. F., Wang, Q. R., et al. (2016). Whole-genome single-nucleotide polymorphism (SNP) marker discovery and association analysis with the eicosapentaenoic acid (EPA) and docosahexaenoic acid (DHA) content in Larimichthys crocea. Peerj 4:e2664. doi: 10.7717/peerj.2664
Xu, H., Dong, X., Ai, Q., Mai, K., Xu, W., Zhang, Y., et al. (2014). Regulation of tissue LC-PUFA contents, Delta6 fatty acyl desaturase (FADS2) gene expression and the methylation of the putative FADS2 gene promoter by different dietary fatty acid profiles in Japanese seabass (Lateolabrax japonicus). PLoS One 9:e87726. doi: 10.1371/journal.pone.0087726
Xu, J., Zhao, Z., Zhang, X., Zheng, X., Li, J., Jiang, Y., et al. (2014). Development and evaluation of the first high-throughput SNP array for common carp (Cyprinus carpio). BMC Genomics 15:307. doi: 10.1186/1471-2164-15-307
Xu, P., Zhang, X. F., Wang, X. M., Li, J. T., Liu, G. M., Kuang, Y. Y., et al. (2014). Genome sequence and genetic diversity of the common carp, Cyprinus carpio. Nat. Genet. 46, 1212–1219. doi: 10.1038/ng.3098
Xu, J., Ji, P. F., Zhao, Z. X., Zhang, Y., Feng, J. X., Wang, J., et al. (2012). Genome-Wide SNP discovery from transcriptome of four common carp strains. PLoS One 7:e48140. doi: 10.1371/journal.pone.0048140
Xu, X., Gopalacharyulu, P., Seppanen-Laakso, T., Ruskeepaa, A. L., Aye, C. C., Carson, B. P., et al. (2012). Insulin signaling regulates fatty acid catabolism at the level of CoA activation. PLoS Genet. 8:e1002478. doi: 10.1371/journal.pgen.1002478
Yao, X., Sa, R., Ye, C., Zhang, D., Zhang, S., Xia, H., et al. (2015). Effects of thyroid hormone status on metabolic pathways of arachidonic acid in mice and humans: a targeted metabolomic approach. Prostaglandins Other Lipid Mediat. 11, 11–18. doi: 10.1016/j.prostaglandins.2015.03.005
Yu, G. C., Wang, L. G., Han, Y. Y., and He, Q. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16, 284–287. doi: 10.1089/omi.2011.0118
Zanetti, M., Grillo, A., Losurdo, P., Panizon, E., Mearelli, F., Cattin, L., et al. (2015). Omega-3 polyunsaturated fatty acids: structural and functional effects on the vascular wall. Biomed. Res. Int. 2015:791978. doi: 10.1155/2015/791978
Zhang, J. G., Tan, L. J., Xu, C., He, H., Tian, Q., Zhou, Y., et al. (2015). Integrative analysis of transcriptomic and epigenomic data to reveal regulation patterns for BMD variation. PLoS One 10:e0138524. doi: 10.1371/journal.pone.0138524
Zhang, Y. Q., Liu, J. H., Fu, W., Xu, W. T., Zhang, H. Q., Chen, S. J., et al. (2017). Comparative transcriptome and DNA methylation analyses of the molecular mechanisms underlying skin color variations in Crucian carp (Carassius carassius L.). BMC Genet. 18:95. doi: 10.1186/s12863-017-0564-9
Keywords: common carp, polyunsaturated fatty acids, GWAS, transcriptome, methylation
Citation: Zhang H, Xu P, Jiang Y, Zhao Z, Feng J, Tai R, Dong C and Xu J (2019) Genomic, Transcriptomic, and Epigenomic Features Differentiate Genes That Are Relevant for Muscular Polyunsaturated Fatty Acids in the Common Carp. Front. Genet. 10:217. doi: 10.3389/fgene.2019.00217
Received: 29 September 2018; Accepted: 27 February 2019;
Published: 15 March 2019.
Edited by:
Hooman Moghadam, SalmoBreed AS, NorwayReviewed by:
Jie Mei, Huazhong Agricultural University, ChinaKieran G. Meade, Teagasc, The Irish Agriculture and Food Development Authority, Ireland
Tereza Manousaki, Hellenic Centre for Marine Research (HCMR), Greece
Copyright © 2019 Zhang, Xu, Jiang, Zhao, Feng, Tai, Dong and Xu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jian Xu, xuj@cafs.ac.cn