Transcriptome analysis reveals key information on improving duck yolk lipid contents induced by dietary fish oil or flaxseed oil

ABSTRACT

Genes participating in the changes of yolk lipids by unsaturated oils remain elusive. In this study, ducks were randomly assigned to 3 groups, fed with basal-diet (C), basal-diet-plus-2%-flaxseed-oil (FL), or basal-diet-plus-2%-fish-oil (F), respectively. A digital-gene-expression-profiling analysis of livers were performed after four weeks. Compared to C, 36 and 25 differently-expressed-genes (DEG) were identified in FL and F, respectively. These genes participated in lipid metabolism pathways. In both F and FL, the sterol 12 α-hydroxylase 1 gene was down regulated, as well as the fatty acid synthase gene and stearoyl-CoA desaturase gene. In FL, the fatty acid desaturase 2 gene (FADS2) and elongation of the very long-chain fatty acids protein 2 gene (ELOVL2) were up regulated, and the CYP2U (cytochrome P450, family 2, subfamily U) was down regulated. However, the FADS2 and ELOVL2 were down regulated in F to possibly prevent the redundancy of body n-3 PUFAs. In conclusion, FL and F had different molecular pathways to reduce yolk cholesterol and increase n-3 PUFAs. However, CYP8B1 and FADS2 were the key genes in the regulation of both of them. α-linolenic acid in the FL might be transformed to longer n-3 PUFAs. In addition, this study adds sequence information to duck.

KEYWORDS:

• Duck (laying)
• RNA-Seq
• dietary oils
• cholesterol
• n-3 PUFAs

Table

Abbreviation	Explaination
n-3 PUFAs	ω-3 polyunsaturated fatty acids
DGE	digital-gene-expression
DEG	differently expressed genes
CYP8B1	sterol 12 α-hydroxylase 1
CYP7A1	cholesterol 7 α-hydroxylase
FASN	fatty acid synthase
ACSL	long chain fatty acyl-CoA synthetase
CPT1	carnitine–palmitoyl CoA transferase 1
SCD	stearoyl-CoA desaturase
ELOVL2	elongation of the very long-chain fatty acids enzyme 2
DHA	docosahexaenoic acid, C22:6
EPA	eicosapentaenoic acid, C20:5
FA	fatty acid
NRC	National Research Council
ALA	α-linolenic acid, C18:3
CYP2U	cytochrome P450, family 2, subfamily U gene
CYP2J	cytochrome P450, family 2, subfamily J gene
CYP2C	cytochrome P450, family 2, subfamily C gene
CYP51	14α-demethylase gene
LCAT	lecithin–cholesterol acyltransferase gene
LPIN	lipin-1 gene
ALDH9A1	aldehyde dehydrogenase 9 family, member A1 gene/ aldehyde dehydrogenase 1 family, member A1 gene
LPL	lipoprotein lipase gene

Introduction

Duck egg is one of the most important components of human diet, especially in Asia. It contains important nutrition, such as fatty acids and cholesterol. However, people are advised to limit their intake of egg yolks due to the high fat and cholesterol contents. Among the yolk fatty acids, omega-3 polyunsaturated fatty acids (n-3 PUFAs), especially DHA (docosahexaenoic acid, C22:6) and EPA (eicosapentaenoic acid, C20:5), are ‘good fatty acids’ that essential for vertebrates and play an important role in their growth, development, and reproduction. Recent studies suggest that the limits to dietary cholesterol should be abandoned, as it almost has no influences on cardiovascular diseases and diabetes in the human (Geiker et al. 2018). However, there are good reasons for longstanding recommendations that dietary cholesterol should be limited to below 200 mg/day in persons at risks of vascular diseases (David 2016). Therefore, much attention has been focused on poultry eggs with a low level of cholesterol and a high level of n-3 PUFAs.

Since the first report from Cruickshank (1934) that the fatty acid (FA) compositions of eggs could be significantly influenced by the addition of various dietary FAs. Other studies also demonstrated that dietary PUFAs could not only alter egg FA compositions but also reduce yolk cholesterol amount in avian species (Couch and Saloma 1973; Liu et al. 2011; Deng et al. 2012). Our previous study (Du X et al., 2017) proved that the ducks with optimal egg-production supplemented with 2% fish oil or 2% flaxseed oil on their basal diets increased yolk n-3 PUFA content and reduced yolk cholesterol content. However, the underlying molecular mechanisms remain undefined.

Recently-developed digital gene expression (DGE) profiling technology has been widely applied to identify specific and abundant genes from millions of sequencing tags among experimental groups. The liver plays an important role in the maintenance of cholesterol and fatty acids homeostasis. It is the major site for the biosynthesis of cholesterol and the conversion of other FAs to PUFAs (polyunsaturated fatty acids) in vivo. In this study, the livers of Shan Partridge ducks were collected, and DGE profiling was performed to elucidate the molecular mechanisms of how flaxseed oil and fish oil lower cholesterol and increase n-3 PUFAs contents in the egg-laying ducks. To increase the reliability of the gene function annotation, the de novo assembly of RNA-Seq tags was used. This study provides molecular information on the role of flaxseed oil and fish oil in the production of eggs with high n-3 PUFA and low cholesterol.

Materials and methods

Animal care

In this study, the Shan Partridge ducks were randomly selected from one hereditary line at Guowei Poultry Industry Co., Ltd. China. All experimental ducks were grown in individual cages under natural temperature and light conditions in May in Zhejiang Province. Feed and water were provided ad libitum. Animals used in our study were cared in strict accordance with the Institutional Animal Care and Use Committee (IACUC) of Zhejiang Academy of Agricultural Sciences.

Experimental design

A flock of 135 Shan Partridge ducks with optimal egg production were randomly selected and divided into 3 groups. Each group contained 45 ducks and were divided equally to 3 pens. The control group (the C group) was fed on a basal diet. The other two groups were fed on basal diets supplemented with either 2% flaxseed oil (FL) or 2% fish oil (F). The ingredients and calculated nutrient levels of the basal diet were shown in Table S1 and were formulated to meet or exceed NRC. The fatty acid (FA) contents of the supplements were shown in Table S2. The ducks were reared for 4 weeks. The laid eggs were collected for three days for the measurement of yolk FA and cholesterol contents. At the end of the experiment, the ducks were sacrificed by cervical dislocation to obtain livers for the de novo transcriptome assembly and digital gene expression analysis. The livers were isolated, frozen immediately in liquid nitrogen, and stored at −80°C until RNA extraction.

Lipid extraction and methyl esterification were done as our previous study (Chen et al. 2017). Briefly, one gram of the yolk sample was ground and mixed with 2 mL of chloroform/methanol (1:2, v/v in 0.001% Butylated hydroxytoluene), 1 mL of chloroform, and 1 mL of water. The mixture was then centrifuged at 1000 r/min for 15 min, and the lower layer was collected subsequently. The above procedure were repeated twice. All the three lower layers were combined together in a tube and evaporated under a stream of N₂ gas. The lipids were then methylated to fatty acid methyl esters (FAMEs) and measured. Briefly, pentadecanoic acid was added to each sample as an internal standard. Derivatization was performed with 1 mL of 3 M methanolic hydrochloric acid at 78°C for 30 min. Then 2 mL of hexane:chloroform (4:1, v/v) were added to each tube. The upper phase was collected after vortex and centrifuge. The above procedure was repeated twice, and the three upper phases were combined and dried under nitrogen gas. The mixture was redissolved in 200 μL of hexane, and transferred into a GC vial with an insert. Then, the lipid compositions were analyzed using a gas chromatography (6890, Agilent, USA) on a HP-5MS column (Agilent 30 m × 0.25 mm × 0.2 μm, Agilent, USA). Helium was used as the carrier gas. The total cholesterol (TC) was determined by the cholesterol oxidase p-aminophenol (COD-PAP) method. Briefly, 1 mL yolk and 9 mL PBS were ground on ice. The mixture was then centrifuged at 1000 r/min for 15 min. Subsequently, the upper layer was collected and measured in a Microplate Reader.

RNA isolation, library construction and sequencing

A two-step approach was applied to analyze the mechanisms of how the flaxseed oil and fish oil diets lowered cholesterol and increased n-3 PUFAs. Three ducks per group were randomly selected. An integrated reference transcriptome was constructed de novo via RNA-Seq by pooling all the 9 samples. Digital gene expression (DGE) libraries were constructed from each duck for differentially expressed gene analysis.

Total RNAs were isolated from the livers using the TRIzol reagent (Life Technologies, USA) according to the manufacturer’s protocol. A Nanodrop ND-1000 spectrophotometer was used to assess the quantity and quality of the RNA samples. All the samples had an A260:A280 ratio >1.8 and an A260:A230 ratio >2.0. An Agilent 2200 TapeStation (Agilent Technologies, USA) was used to evaluate RNA integrity, and all the samples had a RINe >7.0. The messenger RNAs were isolated from total RNA by magnetic Oligo (dT) beads and fragmented to approximate 200 bp. The fragments were reverse-transcribed into double-stranded cDNA. Subsequently, a TruSeq^® RNA LT/HT Sample Prep Kit (Illumina, USA) was used for adaptor ligation and enrichment with a low-cycle. The purified library products were evaluated by the Agilent 2200 TapeStation and Qubit^® 2.0 (Life Technologies, USA) and then diluted to 10 pM for cluster generation in situ on an Illumina HiSeq™ 2500. The de novo transcriptome library was paired-end sequenced at 2 × 100 bp. The 9 libraries for DGE profiling were constructed from each liver sample and single-end sequenced separately at 1 × 50 bp.

De novo transcriptome assembly and annotation

Before the assembly, high quality clean reads were generated by removing adapter sequences and low quality reads in the employ of fastX toolkit. The Trinity were used for the production of contigs by extension through links of overlapped reads. Subsequently, the contigs were further clustered with TGICL software (Pertea et al. 2003) to extract scaffolds. Unigenes with the length equal to or longer than 200 bp were retained to compare with Kyoto Encyclopedia of Genes and Genomes (Kanehisa and Goto 2000) (KEGG, https://www.genome.jp/), Gene Ontology (Ashburner et al. 2000) (GO), Clusters of Orthologous Group of proteins/ Eukaryotic Orthologous Groups of proteins (Tatusov et al. 2003) (COG/KOG), NCBI non-redundant protein database (NR), and SwissProt protein database by BLASTX (version 2.2.14) with an E-value threshold of 1e⁻⁵. The results were noted and gathered to obtain the homologous protein sequences and related pathways. The de novo assembly results served as a reference transcriptome for DGE profiling.

The top BLASTX hits of NCBI NR annotations were used for identifying putative homologous proteins, and each unigene was annotated with GO terms using Blast2 GO (Conesa et al. 2005). The GO classification was performed using WEGO (Ye et al. 2006) according to molecular function, biological process, and cellular component ontologies to provide the background of gene functional category tags and gene functional study. A GO enrichment analysis for the selected gene set with RPKM values >100 was performed and visualized in Cytoscape (v2.8.3). The cytoscape plugin BiNGO (v2.44) (Maere et al. 2005) was customized with the GO annotation results of the Shan Partridge duck de novo transcriptome and then used to perform the hypergeometric statistical test of significance (corrected p-value <0.05) to assess GO term enrichment. To adjust for multiple hypothesis testing, a Benjamini and Hochberg false discovery rate (FDR) correction was performed.

In order to investigate the metabolic pathway annotations of unigenes, the unigenes were aligned using the KEGG database. Enzyme codes were acquired for each sequence, and Enzyme Commission numbers were used to color and retrieve the KEGG pathway maps by KOBAS (Xie et al. 2011) with a corrected p-value <0.05.

Analysis of differentially expressed genes between treatments

The clean reads of three DGE profiles were gathered and mapped against the reference transcriptome with Tophat (v1.4.0) (Trapnell et al. 2009) to get a synthetic evaluation of sequencing coverage regions and depths. Two mismatches per read and 2 gaps per read were permitted in the comparisons. The RPKM of each gene was also calculated and used to estimate the levels of gene expression with p-values adjusted into q-values using the FDR method. Q-values ≤ 0.001 and |log2(Fold change)| > 1 served as the threshold of significance for DEG (Audic and Claverie 1997).

Enrichment analyses of DEG using KEGG were applied. Hypergeometric distribution of each DEG was constructed to analyze the most significantly (P < 0.01) enriched KEGG pathways, through which genes with specific functions were identified.

Quantitative real-time PCR validation

RNAs from 6 individuals per group were prepared as detailed above and reverse-transcribed into cDNA by Reverse Transcriptase M-MLV (RNase H-) (TaKaRa, Code No.: D2639A). Quantitative real-time PCR (qRT-PCR) was implemented for the DEG that had shown a relationship with FA or cholesterol metabolism. The duck beta-actin gene (Actin-B) served as an internal control gene, as it stable expressed in liver and ovary. Primer sequences were designed by Primer Premier 5 and shown in Table 1. QRT-PCR with SYBR Real-time PCR Premix Ex TaqTM (Tli RNaseH Plus) (TaKaRa, Code No.:RR420A) was performed on a 7300 System (Eppendorf, Germany). The results were normalized to the expression level of the constitutive Actin-B, and a relative quantitative method (2^−ΔΔCt) was used to evaluate the quantitative variation.

Table 1. Primers used in qRT-PCR.

Gene	Forward primer (5′-3′)	Reverse primer (5′-3′)	Product size (bp)
Actin-B	ATGTCGCCCTGGATTTCG	CACAGGACTCCATACCCAAGAA	165
FASN	AGCCTGCTCCTCCAGTCTCA	GGTCAAAAGAACGACGACAACA	226
HSD17B7	CCTCTTTGGGCACTTTATCCTG	CGTCAGGTCAGTGGCGTATTTG	184
FADS2	CCTCCTTGCATTCATCCTT	CAGCTTCTTCTTGCCATACTC	277
ELOVL2	TTTCCATTACCATACCCCTCTCC	GTCTTTGCCTTCTGCTACATCCTT	128
ELVL6	CGTCTTCAGAATACTTGGTG	GTTCGGGTGCTTTGCTTAG	158
ACSL	GTGGCAAGAAGTGTGGGGTAG	CTGGCTTTGGAGGTAGAGGTTTT	173
LPL	GAGGGAACCTGATTCAAACG	CATCCAGTCAATAAACATAGCG	124
CYP8B1	ATGCGGGATCGTTTCAGGTT	TGATGTCAGTTGGTGGGCTC	93

Results

Yolk cholesterol and fatty acids

The effect of different oils on yolk cholesterol is shown in Table 2. The yolk cholesterol level in the C group was significantly (P < 0.05) higher than those in the F and FL groups, whereas there was no significant difference between the F and FL groups.

Table 2. Influences of dietary oils on yolk cholesterol content (mmol/L yolk) and fatty acid proportions (g/100 g total fatty acids).

Group	Cholesterol level	SFAs	MUFAs	PUFAs	n-3 PUFAs	ALA	EPA	DHA
C	16.12^a	36.25^a	47.00	16.55^c	0.66^c	0.33^c	0.04^c	0.30^b
FL	14.08^b	30.00^c	44.75	24.96^a	3.34^a	3.32^a	0.55^b	1.58^a
F	14.05^b	32.75^b	39.45	21.73^b	3.04^b	0.75^b	0.77^a	1.40^a

Notes: C: ducks supplied with a basal diet; FL: ducks supplied with a basal diet plus 2% flaxseed oil; F: ducks supplied with a basal diet plus 2% fish oil; SFAs: saturated fatty acids; MUFAs: mono-unsaturated fatty acids; PUFAs: poly-unsaturated fatty acids; n-3 PUFAs: ω-3 polyunsaturated fatty acids; ALA: α-linolenic acid; EPA: eicosapentaenoic acid; the superscripts represent significant differences in the same column.

The effect of different oils on the proportions of yolk FAs is shown in Table 2. There was no significant difference in mono-unsaturated fatty acids (MUFAs) between any 2 groups. Saturated fatty acids (SFAs) in the C group were significantly higher than the FL and F groups, whereas the PUFA content in the C group was significantly lower than the FL and F groups. These results suggested that the diets with 2% flaxseed oil or 2% fish oil could reduce the proportions of SFA and increase the content of PUFAs (both n-6 PUFA and n-3 PUFA) in the yolks. This was in agreement with the findings that supplementation of UFA can increase the concentration of the same UFA in the yolks (Cherian et al. 1996; Antruejo et al. 2011). The concentrations of EPA and DHA in the F group were significantly higher than those in other groups. Significantly increased ALA (α-linolenic acid, C18:3) was observed in the FL group in comparison to the C and F groups.

Illumina sequencing and de novo assembly

To obtain an overall liver transcriptome for Shan Partridge ducks, a mixture of equal amounts of total RNAs isolated from different groups was used to construct the sequencing library. Illumina sequencing produced 61.6 million raw reads, 53.2 million of which were retained following the removal of adapter sequences and low quality reads. The deep sequencing data has been uploaded to the NCBI Sequence Read Archive, the SRA accession of de novo Transcriptome is SRR7146572. 55,426 unigenes were assembled with an N50 of 1943 bp, which was the length of the shortest contig at 50% of the total transcriptome length.

For the annotation and classification of the assembled unigenes, sequence similarity searches were performed against public protein databases: NCBI NR, SwissProt, GO, KOG, and KEGG. There were 27,376 unigenes with a significant similarity (E-value < 1e⁻⁵) to the existing records in the databases (Xls. S1). The transcript abundance and BLASTX Hit ratio of the unigenes were evaluated relative to the sequence length (Figure 1a). Sorting the BLASTX hits by species revealed the top two species: Anas platyrhynchos (36.9%) and Gallus gallus (20.12%) (Figure 1b). The reference transcriptome was further utilized for the analysis of DEG between the dietary treatments of Shan Partridge ducks.

Figure 1. The length and species distribution of unigenes from the de novo assembled liver transcriptome. (A) Species distribution of the unigenes in the Shan Partridge duck de novo liver transcriptome with significant BLASTX Hits. Top scoring BLASTX Hits against the NCBI-NR protein database are depicted. The number of BLAST hits per species is shown on the x-axis. The 8 most represented species with proportions of more than 1% were shown in the histogram. (B) The distribution of sequence lengths was shown. Most of the sequences longer than 1500 had BLASTX hits.

Functional annotation of Shan Partridge duck de novo transcriptome

To understand the putative functional categories of the annotated unigenes, GO was employed to classify the annotated transcripts (Figure 2 and Xls. S2). A total of 12,549 unigenes were annotated with one or more terms. The terms were categorized into GO functional groups distributed across three categories: biological processes (10,027), cell components (8913), and molecular functions (9947). KEGG provides a basic database for the systematic analysis of gene functions in terms of the metabolic networks of gene products (Kanehisa et al. 2012). The Unigenes annotated by KEGG orthologous proteins were mapped to the metabolic pathways. In total, 25,702 unigenes were mapped to 278 reference pathways. In addition to the GO and KEGG analyses, a KOG analysis was performed to further evaluate the function of the assembled unigenes (Table 3). Of 20,577 unigenes assigned to KOG, ‘Signal transduction mechanisms’ (20.6%) was the most represented category, followed by ‘General function prediction’ (19.7%), ‘Post-translational modification, protein turnover, chaperones’ (9.7%), ‘Cytoskeleton’ (9.4%), ‘Extracellular structures’ (7.0%), and ‘Transcription’ (6.8%). These results provide a reference for the investigation of specific functions, processes and pathways, as well as for the identification of novel genes involved in the liver lipid-related metabolism in the Shan Partridge duck.

Figure 2. Gene Ontology classification of the de novo assembled liver transcriptome of the Shan Partridge duck. Gene Ontology annotation results of the genes from the de novo assembled Shan Partridge Duck liver transcriptome were mapped to the categories of GO terms. The y-axis and x-axis indicate the numbers of genes in each cluster and the names of clusters, respectively. Molecular functions, biological processes and cellular components are shown in different colors.

Table 3. KOG functional classification of Shan Partridge duck de novo transcripts.

Description	Code	Unigenes
INFORMATION STORAGE AND PROCESSING
Translation, ribosomal structure and biogenesis	J	4855
RNA processing and modification	A	9219
Transcription	K	13,751
Replication, recombination and repair	L	3326
Chromatin structure and dynamics	B	4140
METABOLISM
Energy production and conversion	C	3712
Carbohydrate transport and metabolism	G	4863
Amino acid transport and metabolism	E	4741
Nucleotide transport and metabolism	F	2750
Coenzyme transport and metabolism	H	1320
Lipid transport and metabolism	I	6383
Inorganic ion transport and metabolism	P	5505
Secondary metabolites biosynthesis, transport and catabolism	Q	2732
CELLULAR PROCESSES AND SIGNALING
Cell cycle control, cell division, chromosome partitioning	D	4649
Nuclear structure	Y	807
Defense mechanisms	V	4476
Signal transduction mechanisms	T	54,730
Cell wall/membrane/envelope biogenesis	M	5609
Cell motility	N	1090
Cytoskeleton	Z	18,981
Extracellular structures	W	13,982
Intracellular trafficking, secretion, and vesicular transport	U	9785
Posttranslational modification, protein turnover, chaperones	O	19,461
POORLY CHARACTERIZED
General function prediction only	R	52,350
Function unknown	S	12,451

Notes: C: ducks supplied with a basal diet; FL: ducks supplied with a basal diet plus 2% flaxseed oil; F: ducks supplied with a basal diet plus 2% fish oil; Description: the description of classifications and subclasses of KOG; Code: the code of KOG; Unigenes: the number of unigenes (unigenes are the splicing result of scaffolds).

Expression analysis and DEG

Nine DGE libraries (C1, C2, C3, FL1, FL2, FL3, F1, F2 and F3, see Materials and Methods) were constructed and sequenced. The deep sequencing data have been uploaded to the NCBI Sequence Read Archive. The SRA accessions are SRR7141761 (C1), SRR7141760 (C2), SRR7141763 (C3), SRR7141762 (FL1), SRR7141765 (FL2), SRR7141764 (FL3), SRR7141767 (F1), SRR7141766 (F2), and SRR7141759 (F3). The high quality of each library was confirmed by the Q30 and GC contents of the reads (Table 4). More than 10 million clean reads were obtained from each library, and 88% to 92.05% of them were mapped to the reference transcriptome, and over 73% of the reads were uniquely mapped. The numbers of reads aligned to a unigene ranged from 1 to 252,361, with a mean of 155.3 for the C group, 157.8 for the FL group, and 160.0 for the F treatment, respectively. Based on the number of uniquely aligned reads, each unigene was assigned an RPKM value, which was normalized to the total reads obtained in each individual library. The mean RPKM values of transcripts for the C, FL, and F groups were 12.35, 12.16, and 12.48, respectively. The numbers of upregulated DEG were 156, 216, and 231, and the downregulated DEG were 740, 667, and 171 in the comparisons of C vs. FL, C vs. F, and F vs. FL, respectively. Interestingly, the C vs. FL comparison had the least upregulated DEG and the most downregulated DEG. Furthermore, the F vs. FL comparison had the most upregulated DEG and the least downregulated DEG. These results suggested that the regulation mechanism of flaxseed oil and fish oil were different. Compared to the F, the FL treatment induced a larger change in gene expression, which was in consistent with the effects on yolk fatty acid content.

Table 4. Summary for DGE sequencing datasets.

Summary (n)	Sequencing libraries
	C1	C2	C3	FL1	FL2	FL3	FO1	FO2	F3
Raw reads	11,334,129	10,806,004	11,701,202	11,346,626	12,182,722	10,447,954	11,691,944	10,752,454	12,080,842
Clean reads	11,061,345	10,551,958	11,474,637	11,124,302	11,798,152	10,111,839	11,296,512	10,415,348	11,702,087
Q30	95.66%	95.68%	95.71%	95.78%	94.94%	94.59%	94.41%	94.62%	94.95%
GC content	54%	53%	53%	52%	52%	52%	53%	53%	52%
Total mapped	9,734,358 (88.00%)	9,665,526 (91.60%)	10,562,944 (92.05%)	10,232,062 (91.98%)	10,598,500 (89.83%)	9,076,976 (89.77%)	9,954,215 (88.12%)	9,469,041 (90.91%)	10,611,861 (90.68%)
Multiple mapped	1,649,979 (14.92%)	1,225,806 (11.62%)	1,260,867 (10.99%)	1,232,397 (11.08%)	1,295,224 (10.98%)	1,136,794 (11.24%)	1,330,855 (11.78%)	1,232,444 (11.83%)	1,372,684 (11.73%)
Uniquely mapped	8,084,379 (73.09%)	8,439,720 (79.98%)	9,302,077 (81.07%)	8,999,665 (80.90%)	9,303,276 (78.85%)	7,940,182 (78.52%)	8,623,360 (76.34%)	8,236,597 (79.08%)	9,239,177 (78.95%)
Reads map to ‘+'	4,759,056 (43.02%)	4,817,970 (45.66%)	5,275,688 (45.98%)	5,104,924 (45.89%)	5,274,943 (44.71%)	4,523,784 (44.74%)	4,945,692 (43.78%)	4,726,670 (45.38%)	5,288,811 (45.20%)
Reads map to ‘−'	4,975,302 (44.98%)	4,847,556 (45.94%)	5,287,256 (46.08%)	5,127,138 (46.09%)	5,323,557 (45.12%)	4,553,192 (45.03%)	5,008,523 (44.34%)	4,726,670 (45.38%)	5,323,050 (45.49%)

Notes: C: ducks supplied with a basal diet; FL: ducks supplied with a basal diet plus 2% flaxseed oil; F: ducks supplied with a basal diet plus 2% fish oil; each group has three ducks; Raw reads: raw reads count in each sample; Clean reads: filtered raw reads count; Q30: the ratio of Phred >30 bases to total bases; GC content: the ratio of G and C numbers to total bases; Total mapped: total reads mapped to reference genome: Multiple mapped: reads that mapped two or more sites to reference genome; Uniquely mapped: reads that mapped one site to reference genome; Reads map to ‘+': reads mapped to sense strands, Reads map to ‘−': reads mapped to anti-sense strands;

Functional annotation of DEG

The KEGG pathway analysis of DEG was also performed. As expected, the analysis indicated that some DEG resulted from the C vs. FL, C vs. F, and F vs. FL comparisons were related to lipid metabolism pathways (Table 5).

Table 5. DEG related to lipid metabolism and the basic information.

Gene name	Abbreviation of genes	Function of genes	Regulation of FL	Regulation of F
cyclooxygenase 1/prostaglandin H2 synthase 1 gene	COX-1/PTGS1	It is the key enzyme converts ARA to prostaglandin.	Up^a	–^b
cytochrome P450, family 2, subfamily U genes	CYP2U	Enzymes in this subfamily catalyze omega/omega-1 hydroxylation of ARA and the hydroxylation of DHA (Chuang et al. 2004).	Down^c	–
cytochrome P450, family 2, subfamily J genes	CYP2J	Enzymes in this subfamily are active in the metabolism of fatty acids to generate bioactive compounds (Graves et al. 2015).	–	up
cytochrome P450, family 2, subfamily C genes	CYP2C	Enzymes in this subfamily oxidize ARA to epoxyeicosatrienoic acids (EETs) (Lundblad et al. 2005).	dwon	up
14α-demethylase gene	CYP51	It serves as a key player in the transformation of steroids, lipids and vitamins.	up	–
cytochrome P450, family 8, subfamily B, polypeptide 1 gene	CYP8B1	It is known as sterol 12-alpha-hydroxylase, which is important for the synthesis of cholesterol.	down	down
7-alpha-monooxygenase gene	CYP7A1	It is the key enzyme that convert cholesterol to bile acid.	up	–
fatty acid synthase gene	FASN	The main function is synthesizing palmitate.	down	down
long chain fatty acyl-CoA synthetase	ACSL	It catalyzes the formation of acetyl-CoA to breakdown complex fatty acids.	down	down
elongation of very long chain fatty acids protein 2 gene	ELOVL2	Its main function is synthesis of PUFAs, it is the key enzyme that catalyzed ALA to DHA.	–	down
elongation of very long chain fatty acids protein 6 gene	ELOVL6	Its main function is synthesis of SFAs and MUFAs.	down	–
fatty acid desaturase 2 gene	FADS2	It is important in the synthesis of PUFAs.	up	down
lecithin–cholesterol acyltransferase gene	LCAT	It converts cholesterol to cholesterol ester for the synthesis of HDL, LDL and VLDL.	down	down
aldehyde dehydrogenase 9 family, member A1 gene/ aldehyde dehydrogenase 1 family, member A1 gene	ALDH9A1/ ALDH1A1	It catalyzes the dehydrogenation of gamma-aminobutyraldehyde to gamma-aminobutyric acid.	down	down
lipoprotein lipase gene	LPL	It hydrolyzes triglyceride in lipoprotein.	up	–
stearoyl-CoA desaturase gene	SCD	It catalyzes the forming of double bond in stearoyl-CoA, it is responsible for the synthesis of MUFAs from SA.	down	down
carnitine palmitoyltransferase 1 gene	CPT1	It is carnitine acyltransferase that responsible for the formation of acyl carnitines by catalyzing the transfer of the acyl group of a long-chain fatty acyl-CoA from coenzyme A to l-carnitine.	up	up
acyl-CoA Synthetase, Bubblegum Family, member 1 gene	ACSBG	It is thought to play a central role in brain very long-chain fatty acids metabolism and myelinogenesis.	–	down
diacylglycerol O-acyltransferase 1 gene	DGAT1	It catalyzes the conversion of diacylglycerol and fatty acyl CoA to triacylglycerol.	up	up
7-Dehydrocholesterol reductase gene	DHCR7	It catalyzes the production of cholesterol from 7-Dehydrocholesterol using NADPH.	–	up
ectonucleotide pyrophosphatase/phosphodiesterase family member 1 gene	ENPP2	It is important for generating the lipid signaling molecule lysophosphatidic acid (LPA).	up	up
farnesyl-diphosphate farnesyltransferase 1 gene	FDFT1	It catalyzes the dimerization of two molecules of farnesyl diphosphate in a two-step reaction to form squalene.	up	–
galactocerebrosidase gene	GALC	It removes galactose from ceramide derivatives.	down	down
17β-Hydroxysteroid dehydrogenase	HSD17B7	It oxidizes or reduces estrogens and androgens in mammals and regulates the biologic potency of these steroids.	up	down
lipin-1 gene	LPIN	It has phosphatidate phosphatase activity.	down	down
membrane bound O-acyl transferase gene	MBOAT1	It catalyzes the transfer of an acyl group from acyl-coenzyme to one of several different substrates.	down	down

Notes: FL: ducks supplied with a basal diet plus 2% flaxseed oil; F: ducks supplied with a basal diet plus 2% fish oil.

^aUp regulated by FL/F compared to C; ^bNo significant effects induced by FL/F; ^cDown regulated by FL/F.

Compared to the C diet, the FL diet significantly reduced the expression of cholesterol synthesis related gene, CYP8B1 (cytochrome P450, family 8, subfamily B, polypeptide 1), and significantly increased cholesterol decomposition related gene, CYP7A1 (7-alpha-monooxygenase). On the other side, the FL diet significantly reduced the expression of DHA decomposition related gene, CYP2U (cytochrome P450, family 2, subfamily U), and significantly increased the expression of n-3 PUFAs synthesis related gene, FADS2 (fatty acid desaturase 2). In addition, the FL diet regulated SCD (stearoyl-CoA desaturase) and FASN (FA synthese), these two genes are related to short-chain fatty acid synthesis. Some other differently expressed genes were also detected by comparing C vs FL, including glyceride metabolism related genes, such as LCAT (lecithin–cholesterol acyltransferase), LPIN (lipin-1), ALDH9A1 (aldehyde dehydrogenase 9 family, member A1), and LPL (lipoprotein lipase).

Compared to the C diet, the F diet decreased the expression of CYP8B1. On the other hand, the F diet significantly increased long chain fatty acids decomposition related gene expression, ACSL (long chain fatty acyl-CoA synthetase), and reduced n-3 PUFAs synthesis related gene expression, ELOVL2 (elongation of very long chain fatty acids protein 2) and FADS2. In addition, the F diet regulated the expression of CYP2C (cytochrome P450, family 2, subfamily C), CPT1 (carnitine palmitoyltransferase 1), SCD, and FASN, which are possibly related to the reduction of short-chain fatty acid content.

Quantitative real-time PCR validation

In order to validate the differential gene expression, eight genes were selected for qRT-PCR analysis. The results are shown in Figure 3. The gene expression results obtained through qRT-PCR mostly corroborated DEG data, suggesting that the DGE profiling analysis provided reliable data for the Shan Partridge ducks fed with fish oil or flaxseed oil diets.

Figure 3. Verification of the expression of 8 DEG by qRT-PCR in the C, FL, and F groups. Actin-B was used as a reference gene for normalization of gene-expression data. The left Y-axis represents relative gene expression levels measured by qRT-PCR, whereas the right Y-axis represents FPKM of the genes measured by RNA-Seq for the stages using the same RNA samples. C: control group; FL: flaxseed oil group; F: fish oil group. All qRT-PCR data are shown as mean ± standard error.

Discussion

In the study, we confirmed that the 2% flaxseed oil or 2% fish oil diet supplements lowered cholesterol and SFA contents, but increased PUFA contents in the yolks of Shan Partridge ducks. This finding was consistent with previous studies on egg-laying hens that were fed dietary flaxseed oil or fish oil (Basmacioglu et al. 2004; Liu et al. 2011; Du et al. 2017). The results of the present study confirmed that the diets supplemented with 2% flaxseed oil or 2% fish oil for 28 days had the same cholesterol-lowering effect on the yolks of Shan Partridge ducks. It is worth noting that the diet plus 2% flaxseed oil, which has the same EPA and DHA content as our C group diet, significantly increased EPA and DHA in the yolks, suggesting that ALA in the FL group diet might be transformed into EPA and DHA in the liver.

To analyze the DEG in the liver induced by the F and FL diet supplements, the duck liver transcriptome was de novo assembled with an N50 of 1943 bp, and 27,376 unigenes were identified and annotated. The top two species that had the most BLASTX hits with the transcriptome were Anas platyrhynchos (36.9%) and Gallus gallus (20.12%). It is important to note that Anas platyrhynchos is a closely-related species of the Shan Partridge duck. However, due to the incomplete sequence information of Anas platyrhynchos in the public databases, the hit ratio was not very high. The results indicated that a large number of specific genes were present in our de novo assembly, which will provide a resource for the development of functional molecular markers of Shan Partridge ducks.

Differentially expressed genes were further identified between the diet groups, which were annotated functionally via GO functional enrichment and KEGG pathway analyses. GO functional enrichment analysis results showed that the most significantly enriched molecular function and cellular component induced by FL and F supplements were ATP binding and cytoplasm, respectively. The enriched function and cellular component of the F vs. FL comparison were metal ion binding and integration to the membrane, respectively. Metal ion binding is important for the structure and function of proteins and other functional molecules. Compared to flaxseed oil, fish oil contains more DHA and EPA, which are important for the integrity of the membranes (Giusto et al. 2000). The FL diet led to the enrichment of DEG in the oxidation–reduction process compared to the C or F diets. This is possibly because FL diet is ALA rich, which is the precursor to synthesize DHA and EPA in an oxidation–reduction process.

KEGG pathway analysis of the DEG suggested three key points. First, both fish oil and flaxseed oil reduced yolk cholesterol possibly by decreasing the expression of CYP8B1, which is an important enzyme that catalyzes the monooxygenation at the 12α position of 7α-hydroxy-4-cholesten-3-one in the synthesis of cholesterol. Previous studies have shown that CYP8B1 feedback regulates cholesterol 7α-hydroxylase (CYP7A1), which is the rate-limiting enzyme of cholesterol degradation (Li-Hawkins et al. 2002). In addition to CYP8B1, the flaxseed oil diet directly increased the expression of CYP7A1, whereas the fish oil diet had no significant effects on this gene.

Second, both F and FL inhibited FA biosynthesis by reducing the expression of FASN (Munoz et al. 2003). FASN is a key enzyme of lipogenesis. In this study, both fish oil and flaxseed oil supplied sufficient basic FAs for the ducks; thus, the body automatically reduced de novo FA synthesis.

Third, both the fish oil and flaxseed oil diets upregulated CPT1 and downregulated SCD in the biosynthesis of UFAs. Additionally, the F diet significantly downregulated ELOVL2 and upregulated ACSL, whereas flaxseed oil did not. CPT1 is an enzyme responsible for the transport of acyl groups of long-chain fatty acid-CoA into mitochondria for oxidation (Guo et al. 2012). SCD catalyzes the rate-limiting step in the formation of MUFAs. ELOVL2 is essential for the conversion of dietary ALA to DHA as it can elongate DPA (22:5n-3) to 24:5n-3, which is required for DHA synthesis (Gregory and James 2014). ACSL is a ligase that catalyzes the breakdown of complex fatty acids (Soupene and Kuypers 2008). It is suggested that when supplied with oils rich in PUFAs, ducks changed their metabolic pathways to promote the decompostion of fatty acids, enhance the synthesis of DHA and EPA, and inhibit the de novo synthesis of FAs. As fish oil is also rich in DHA and EPA, which can satisfy the body demands, the synthesis key enzyme, ELOVL2, was downregulated to prevent energy dissipation.

In conclusion, both flaxseed oil and fish oil supplements reduced yolk cholesterol and increase yolk n-3 PUFAS. However, the flaxseed oil supplements induced a larger change in gene expression than fish oil supplements. CYP8B1 and FADS2 were the key genes in the regulation of both flaxseed oil and fish oil supplements. ALA in the flaxseed oil might be transformed to n-3 PUFAs, such as DHA and EPA. For another, we assembled de novo duck liver transcriptome as a resource for the further studies.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This project was supported by the Twelfth Five Year Plan of National Science and Technology Project in Rural Areas [grant number 2015BAD03B06]; Zhejiang Major Science and Technology Project of Agricultural (Livestock’s) Breeding [grant number 2016C02054-12]; National Waterfowl Industry Technology System [grant numbers CARS-43-6, CARS-43-39]; and the Natural Science Foundation of Zhejiang Province [grant number LQ17C170003].