Transcriptome analysis of the livers of ducklings hatched normally and with assistance.

OBJECTIVE: "Hatchability" is an important economic trait in domestic poultry. Studies on poultry hatchability focus mainly on the genetic background, egg quality, and incubation conditions, whereas the molecular mechanisms behind the phenomenon that some ducklings failed to break their eggshells are poorly understood.
METHODS: In this study, the transcriptional differences between the livers of normally hatched and assisted ducklings were systematically analyzed.
RESULTS: The results showed that the clean reads were de novo assembled into 161,804 and 159,083 unigenes (≥200-bp long) by using Trinity, with an average length of 1,206 bp and 882 bp, respectively. The defined criteria of the absolute value of log2 fold-change ≥1 and false discovery rate≤0.05 were differentially expressed and were significant. As a result, 1,629 unigenes were identified, the assisted ducklings showed 510 significantly upregulated and 1,119 significantly down-regulated unigenes. In general, the metabolic rate in the livers of the assisted ducklings was lower than that in the normal ducklings; however, compared to normal ducklings, glucose-6-phosphatase and ATP synthase subunit alpha 1 associated with energy metabolism were significantly upregulated in the assisted group. The genes involved in immune defense such as major histocompatibility complex (MHC) class I antigen alpha chain and MHC class II beta chain 1 were downregulated in the assisted ducklings.
CONCLUSION: These data provide abundant sequence resources for studying the functional genome of the livers in ducks and other poultry. In addition, our study provided insight into the molecular mechanism by which the phenomenon of weak embryos is regulated.

INTRODUCTION

“Hatchability” is directly responsible for the economic benefit of a hatchery. Schaal and Cherian investigated the hatchability of broiler and turkey eggs in the U.S. from 1985 to 2005 and found that advances in genetic selection, nutrition, and management of broiler and turkey flocks during this time period did not lead to an increase in hatchability [1]. The industry economic loss was >$500 million in 2005, because of the lack of improved hatchability [1]. Hatchability of artificially incubated duck eggs is low (65% to 82%) compared with that of domesticated chickens (81% to 85%); therefore, this process has higher positive potential in the duck industry. Several factors that relate to hatchability were explored, which can be loosely divided into 2 groups: breeder factor and hatchery factor. The former mainly includes nutrition, strain and age of the flock, egg characteristics, storage time and conditions [2-5]. Hatchery factors include egg handling, storage, and incubation and conditions such as humidity, temperature, and turning frequency [6,7]. Reasons for low hatchability might be attributed to the impropriety of one or more factors mentioned, or all of the above. Hodgetts [7] reported that the variations in egg sizes, age, and degree of contamination are the main factors influencing on the hatchability of ducks. Harun et al [8] reported that hatching ability was associated with egg characteristics and metabolic rate. Normally hatched ducklings had a higher metabolic rate than that in non-hatched and assisted ducklings (i.e., those that began to pip but did not completely hatch from the egg after more than 48 h were assisted in breaking the eggshell) [8]. Although studies on a variety of factors that influence hatchability and result in the physiological response for embryonic development during the incubation period have been extensively conducted in birds, the molecular basis for the assistance that ducklings require in hatching remains largely unknown.

Recently, the next-generation high-throughput sequencing methods, termed RNA sequencing (RNA-seq), have been widely applied used for transcriptome analysis in eukaryotes such as rice, peanut, humans, yeast, and mice [9-12]. Compared with other high-throughput sequencing technologies, Illumina sequencer can generate large amounts of data in a very cost-efficient manner, and allow for de novo assembly for non-redundant (NR) sequences to provide comprehensive sequence resources and characterize RNA expression patterns in non-model organisms [13,14].

The liver, a central organ of intermediary energy metabolism, plays a major role in processes such as gluconeogenesis, glycogenesis, glycogenolysis, and lipid metabolism. It also provides a variety of functions, including protein synthesis, hormone production, and detoxification [15]. Moreover, the liver has important functions in innate and adaptive immunity [16]; Therefore, it was selected for sampling to explore the differences in metabolism and immune defense between the normally hatched and assisted ducklings using the Illumina RNA-Seq method. In this study the influencing factors of breeder and hatchery were excluded as much as possible. Under the circumstances, we identified a large number of differentially expressed genes in the livers of assisted and normal ducklings. In addition, this de novo transcriptome assembly method can provide the platform by which to generate novel transcriptomic resources for non-model species whose complete genomes are not available.

MATERIALS AND METHODS

Ethics statement

The animal care and use protocols were approved by the Institutional Animal Care and Use Committee of Zhejiang A&F University and conducted in accordance with the Regulations for the Administration of Affairs Concerning Experimental Animals (China, 1988) and the Measures for the administration of experimental animals (Zhejiang, China, 2009).

Liver collection

We collected 120 fertile eggs in 1 d from a batch of Shaoxing ducks at age of 300 d that had a similar genetic background and the same breeding management, and then measured several indices including egg weight, eggshell thickness, and egg shape, the mean values of which were 63.96±5.10 g, 0.34±0.03 mm, 1.36±0.05, respectively. The 60 eggs within an average value for the 3 characteristics were chosen to be incubated at 37°C and 60% relative humidity. After 27 d of incubation, ducklings that began to pip and either could or could not escape from the eggs after >48 h were defined as either “normally hatched” or “assisted,” respectively (hereinafter referred to as normal ducklings and assisted ducklings, respectively). From each group, 6 randomly selected ducklings were killed by cervical dislocation, and all liver tissues were separated and immediately snap frozen in liquid nitrogen.

RNA isolation and cDNA library construction

Total RNA was isolated from the livers of the 6 ducklings in each group using Trizol reagent (Takara, Osaka, Japan) following the manufacturer’s protocol. The RNA samples were treated with RNase-free DNase I (Takara, Japan) to remove residual genomic DNA. The RNA quality and quantity were checked using a spectrophotometer (NanoVue, GE Healthcare, Piscataway, NJ, USA), and the integrities were examined using 1.2% formaldehyde degeneration gels containing 0.1% ethidium bromide. Equal quantities of RNA from each tissue of a group (n = 6) were pooled for cDNA library construction. In brief, 10 μg of total RNA from 2 pools were used to isolate poly (A) mRNA by using oligo (dT) magnetic beads. Then, the mRNA was digested by a fragmentation buffer to produce short pieces, which were used for first-strand cDNA synthesis by using random hexamer primers, followed by the synthesis of second-strand cDNA by using a buffer, dNTPs, RNase H, and DNA polymerase I. After the short fragments were end repaired and ligated with two different adapters, the fragments of approximately 200 bp were separated by agarose gel electrophoresis and followed by polymerase chain reaction (PCR) amplification. Finally, the cDNA library was built and sequenced in two Hiseq 2000 lanes.

De novo assembly and functional annotation

Trinity [17] was used to de novo assemble the clean reads, which were obtained by removing the adaptor and low-quality sequences from the raw reads. The overlapped reads were first connected to form contigs, and then mapped back to the contigs so that the different contigs and their distances from the same transcript were estimated and assembled to form the longer sequence, which were unigenes. The contigs were again mapped back to the unigenes to obtain the final unigenes that could not be extended on either end. The other sample was assembled according to the same pipeline as described above. For transcriptomic analysis at the macro level, the unigenes from the two samples were assembled and redundant sequences removed to obtain all unigenes that were as long as possible.

All unigenes from the de novo assembly of the two samples were compared to the NCBI NR protein and the UniProt database (http://uniprot.org) with BLASTX (E-value<10-5) to ascertain putative annotation. All unigenes were annotated using gene ontology (GO) terms, and placed into the three categories and functional groups using Blast2GO (http://www.blast2go.org/). The Kyoto Encyclopedia of Genes and Genomes (KEGG) database [18] was used to annotate the pathway of all unigenes.

Analysis of differential expression

Unigene expression level was estimated by counting the fragments per kilobase exon per million mapped fragments (FPKM) by using edgeR [19]. The different expression of unigenes was identified between two samples based on the assembled unigenes as a reference, and unigenes were considered significant when the false discovery rate (FDR) <0.05 and the absolute value of log2 ratio >1.

Quantitative reverse transcription polymerase chain reaction validation

Quantitative reverse transcription polymerase chain reaction (qRT-PCR) was performed to validate four differentially expressed genes in the cDNA pools composed of six individual embryo livers from each group. Gene-specific primers were designed using Primer 5.0 and shown in Table 1. β-actin was used as a reference gene. qRT-PCR was run in triplicate using ABI 7300 (Applied Biosystems, Foster City, CA, USA) at 94°C for 3 min, 94°C for 10 s, and 40 cycles at 60°C for 30 s. RT-PCR was performed in a 25-μL reaction mixture containing 1 μL cDNA template, 1× THUNDERBIRD SYBR qPCR Mix, 1× ROX reference dye (TOYOBO, Tokyo, Japan), and 0.4 μm of each primer. The relative expression levels of the tested genes were calculated as 2−ΔCt method (ΔCt = Ct target gene − Ct β-actin). Data are presented as means±standard error of the mean. The differences between groups were examined using a t test and considered significant when p<0.05.

Table 1

Primer sequences used for qRT-PCR validation of the differentially expressed genes in the livers of normal and assisted ducklings

Gene name	Binding sites	Nucleotide sequences (forward and reverse) (5′→3′)	Product size (bp)
ATP5A1 F	22–42	GGTATCAGCTCCCAAAATACG	246
ATP5A1 R	249–267	CAATCCCACTCTATCTCCAC
FEZ2 F	545–564	CAGCAGAACAGGTGATTGAA	105
FEZ2 R	630–649	GGGAAAGTATGGAGAGACGA
MHC I F	512–530	GAGGAAGGGACTGATGCTGA	196
MHC I R	688–707	GAGATGGGTCGTGGGTAGAA
MHC II F	482–501	GTAGGTCCAGTCCCCGTTCT	184
MHC II R	647–665	AGTTGAGCCCAAGGTGAGG
β-actin F	685–701	ATGTCGCCCTGGATTTCG	164
β-actin R	828–848	CACAGGACTCCATACCCAAGAA

qRT-PCR, quantitative reverse transcription polymerase chain reaction.

RESULTS

RNA-seq and de novo assembly of liver transcriptome

RNA-Seq yielded 77,572,048 and 79,921,784 raw reads with 101-bp paired-ends from the livers in normal and assisted ducklings, respectively. Low-quality reads and adapter sequences were removed by screening reads of >25 bp, resulting in 73,975,798 and 77,977,436 clean reads, corresponding to 7.87 and 7.47 Gb from normal and assisted ducklings, respectively. The percent guanine–cytosine (GC) of each sample was 47.08 and 46.85 in normal and assisted ducklings, respectively. Then, the clean reads were de novo assembled into 159,083 and 161,804 unigenes (≥200 bp long) using Trinity, with an average length of 882 bp and 1,206 bp, and an N50 length of 1,769 bp and 3,059 bp in normal and assisted ducklings, respectively. The length distribution of assembled unigenes for each sample is presented in Figure 1. Table 2 shows the statistical information of RNA-seq reads and de novo assembly of duck-liver transcriptome. The unigenes from 2 samples including overlapped sequences were further assembled and the redundancy removed to acquire 74,046 NR unigenes (all unigenes). The merged transcriptome sequences were deposited in NCBI’s Transcriptome Shotgun Assembly (TSA) database (http://www.ncbi.nlm.nih.gov/, accession number: SRP026134).

Figure 1

Size distributions of the unigenes. (A) and (B) indicate the size distribution of unigenes (>200-bp long) generated by de novo assembly in normal and assisted ducklings, respectively.

Table 2

Summary of transcriptome sequences and de novo assembly from normal and assisted ducklings, respectively

	Normal eggs	Assisted eggs
Raw reads number	77,572,048	79,921,784
Data product (G)	7.83	7.99
Clean reads number	73,975,798	77,977,436
Read length (bp)	101	101
Unigene number (≥ 200 bp)	159,083	161,804
Unigene N50 length (bp)	1,769	3,059
Average unigene length (bp)	862	1,206
GC percent (%)	47.08	46.85

GC, guanine–cytosine.

Annotation and classification of the unigenes

All unigenes aligned to NCBI NR and UniProt protein databases were merged according to optimal matching (BLASTX-hit) using BLASTX with an E-value cutoff of 1e-5. As a result, 74,045 unigenes were matched to 20,804 known proteins. In addition, a majority of BLASTX-hit unigenes (51,982; 73.39%) showed >70% sequence identity, in which 2,238 sequences revealed the identity of 100% (Figure 2). We further determined the distribution of top hit species based on nr annotation and found that 43% of total annotated unigenes were similar to Gallus gallus sequences, followed by Meleagris gallopavo (15%), Anolis carolinensis (6%), and Monode lphisdomestica (4%) (Figure 3A). Only 0.19% of the top matched hits were for Anas platyrhynchos, which suggested that the duck protein sequences available in the NCBI database were very limited.

Figure 2

Sequence identity distributions. Vertical histogram of green and red represents the number of all BLAST-hit transcripts and annotated proteins.

Figure 3

(A) Species distribution of the best hit unigenes. (B) Gene ontology (GO) annotation results. All assembled transcripts were annotated in three major functional categories: cellular component, molecular function, and biological process.

To gain further insight into the distribution of transcript functions in duck liver at the macro level, GO functional annotation of all assembled transcripts was performed using BLAST2GO [20] (http://www.blast2go.org/), after which GO functional categorization of all unigenes was conducted using R. In total, 26,332 unigenes were assigned 34 GO terms, which are distributed under three main classifications (cellular component, molecular function, and biological process), and “cell part”, “binding” and “cellular process” were the most representative categories, respectively (Figure 3B).

In addition, 8,965 unigenes were assigned to 136 KEGG pathways. Purine metabolism was the most predominant, with 2,129 members, followed by the cysteine, methionine, and pyrimidine metabolism, and the phosphatidylinositol signal system, with 1,018 members, 1,004 members, and 752 members, respectively.

Analysis of differentially expressed genes in the livers of assisted and normal ducklings

FPKM was used to quantify the expression level of the transcriptome data from the livers of assisted and normal ducklings, and the criteria of log2 fold-change ≥1 and FDR≤0.05 were defined as significantly differentially expressed. As a result, 1,629 unigenes were identified, from which 510 were significantly upregulated and 1,119 significantly down-regulated in the assisted ducklings compared to those in the normal ducklings. A volcano plot was used to portray the magnitude distribution of the differentially expressed unigenes (Figure 4A). Among these, 886 were annotated to 413 unique protein accessions using BLASTX against NR and UniProt. The functions of the 413 genes were determined by consulting literature materials but not GO annotation, because, on the one hand, the number of annotated transcripts in ducklings could be confined, and on the other hand, the majority of annotated unigenes were matched to known genes of Gallus gallus, which have not been extensively studied.

Figure 4

(A) Volcano plot of differentially expressed unigenes. The horizontal axis represents the −log error value, and the vertical axis indicates the fold change values of expression levels for significantly different unigenes in the assisted and normal ducklings. (B) Quantitative polymerase chain reaction (qPCR) analysis of the selected genes in the livers of assisted and normal ducklings. The mRNA level of each gene was normalized to that of β-actin. Each column represents the mean of gene expression±standard error in 4 livers of normal and assisted ducklings, respectively. ** Indicates a significant difference at p<0.01.

In this study, we focused on the genes involved in metabolic and immune system processes. The results showed that 61 and 24 genes involved in metabolic and immune-defense processes, respectively, were identified. Among the genes in the assisted group, 61 genes performed metabolic functions, of which 14 and 47 genes were upregulated and downregulated, respectively. Of the 24 genes associated with immune defense, 5 were upregulated and 19 were downregulated.

qRT-PCR validation of differentially expressed genes in the livers of assisted and normal ducklings

To further confirm the differential gene expression obtained from the sequencing data, we analyzed by quantitative RT-PCR the expression levels of ATP synthase subunit alpha 1 (ATP5A1), fasciculation and elongation protein zeta-2 (FEZ2), MHC class I antigen alpha chain (HLA-A), and MHC class II beta chain 1 (BLB1) genes, which were mainly involved in metabolism and immune defense. The results verified that the expression level of ATP5A1 was significantly upregulated, and those of FEZ2, MHC I, and MHC II genes were significantly downregulated in the livers of assisted ducklings compared to those in the livers of the normal ducklings (p<0.01), which was consistent with the RNA-seq results (Figure 4B).

DISCUSSION

In June 2013, the entire duck genome sequence was presented and published in Nature Genetics [21]. The reads by RNA-seq in our study were conducted by de novo assembly using Trinity because of the lack of genomic sequence information at that time. As a result, 74,045 unigenes were matched to 20,804 known proteins after all unigenes were annotated, which provides abundant sequence resources for studying the functional genome and identifying genes of interest in the livers of ducks and other poultry. The differential analysis showed that 1,629 unigenes were significantly changed, of which 886 were annotated to known proteins by BLASTX against NR and Uniprot databases; however, some unigenes could match with the same protein and the reasons could be that the unigenes were presented in a different isomer or the different sequences of the same gene.

According to Harun et al [8] reports that hatching ability was associated with egg characteristics and metabolic rate, and liver is an important organ that is directly involved in the metabolic and immune-response processes; therefore, we paid close attention to the different genes involved in metabolism and the immune process between assisted and normal ducklings. We found that the mRNA expressions of 14 genes involved in the metabolic process were upregulated and 47 were downregulated in assisted ducklings compared with those in normal ducklings. The results indicated that the metabolic rate in the livers of assisted ducklings was lower than that in the livers of normal ducklings. Moreover, it was worth noting that compared to the normal ducklings, in the assisted ducklings, the mRNA level of glucose-6-phosphatase (G6PC) and ATP5A1 genes were significantly upregulated.

It has been estimated that >90% of the energy requirement of the developing embryo is derived from the oxidation of yolk lipids [22]. Embryos require sufficient oxygen to oxidize yolk nutrients. As the embryo develops, oxygen gradually becomes limited and the embryos switch to anaerobic metabolism [23], through which they gradually decrease fat metabolism and increase glycogen metabolism. G6PC, a critical enzyme for providing glucose during starvation and diabetes [24], is involved in the terminal step in gluconeogenic and glycogenolytic pathways and catalyzes the hydrolysis of glucose 6-phosphate (G6P) to glucose and inorganic phosphate. G6PC, highly expressed in the liver and kidneys and less expressed in the intestines and pancreas [25], maintains the blood glucose levels during the fasting state. It is reported that long-chain fatty acids promote the G6Pase mRNA expression in cultured fetal hepatocytes [26] and adult rat liver [27]. It is of interest that the mRNA level of elongation of very long-chain fatty acid protein 2 genes (ELOVL2) was significantly upregulated in assisted ducklings compared to that in normal ducklings. ATP5A1 is an important protein in the electron transport chain during oxidative phosphorylation that releases highly efficient energy and produces ATP [28,29]. According to the above analysis, we can conclude that the overall metabolism of assisted ducklings was lower than that of normal ducklings, but that some genes involved in energy metabolism were significantly upregulated in the assisted ducklings. These changes reflect that the assisted ducklings must accelerate their gluconeogenic and glycogenolytic pathways to provide sufficient energy to resolve the challenges that they face when hatching from the eggs; therefore, the regulation mechanism of these genes involved in energy metabolism in assisted ducklings needs to be further researched.

In our study, the expression level of genes mainly related to the immune system from assisted ducklings was significantly downregulated compared to those in normal ducklings, except for alpha-2-macroglobulin-like protein 1 (A2ML1), C-C motif chemokine 4 homolog (CCL4), complement receptor type 2 (CR2), tripartite motif protein 39 (TRIM39), and immunoresponsive 1 homolog (IRG1). The reasons might be that i) the development of the liver was imperfect or the presence of liver disorder in assisted ducklings or ii) assisted ducklings were infected by some kind of virus or bacteria, wherein although the immune system is defensive under normal circumstances, it failed to perform its normal function because of the immune disorder. One of the important functions of the immune system is to protect against infections by invading pathogens. The immune system is able to attack pathogens through both innate and adaptive immunity [30], and the liver has a number of major functions in these 2 fundamental pathways. The innate immune system in the liver is made up of Kupffer cells, dendritic cells, natural killer (NK) and NK T (NKT) cells, inflammatory cytokines, chemokines, acute phase proteins, and complement cells that coordinate to eliminate invading pathogens. The acute-phase response is launched when macrophages release cytokines in response to infection or injury, and upregulated cytokines increase the expression of intercellular adhesion molecule 1 (ICAM-1) [31]. Then, these adhesion molecules trap activated CD8+ T cells in the mouse liver [32]. In acute-phase proteins (APPs), A2ML1, CR2, and CCL4 mRNA levels were upregulated in assisted ducklings compared with those in normal ducklings. The mRNA levels of interferon-induced guanylate-binding protein 1 (GBP1), interferon alpha-inducible protein 6 (IFI6), beta-defensin 9 (GAL9), ICAM-1, HLAA, and BLB1 genes were significantly downregulated in the liver of assisted ducklings compared to that in normal ducklings. Central to maintaining immunity is that MHC class I and class II molecules present antigenic peptides at the surface of CD8+ and CD4+ T cells, respectively. This recognition process can eliminate the infected cells and broadly influence other immune responses. The MHC class is a critically important gene cluster that maintains immunity against viruses. This may indicate that the embryonic liver of assisted ducklings was infected or inflamed, and the body needed to initiate an immune response to restore homeostasis; however, the body failed to recognize and induce the appropriate defense signals. In conclusion, the low rate of embryonic metabolism and immunity responses are very important reasons for the hatching failure in ducklings.