Comparative transcriptome analysis provides clues to molecular mechanisms underlying blue-green eggshell color in the Jinding duck (Anas platyrhynchos)
© The Author(s). 2017
Received: 15 February 2017
Accepted: 7 September 2017
Published: 12 September 2017
In birds, blue-green eggshell color (BGEC) is caused by biliverdin, a bile pigment derived from the degradation of heme and secreted in the eggshell by the shell gland. Functionally, BGEC might promote the paternal investment of males in the nest and eggs. However, little is known about its formation mechanisms. Jinding ducks (Anas platyrhynchos) are an ideal breed for research into the mechanisms, in which major birds lay BGEC eggs with minor individuals laying white eggs. Using this breed, this study aimed to provide insight into the mechanisms via comparative transcriptome analysis.
Blue-shelled ducks (BSD) and white-shelled ducks (WSD) were selected from two populations, forming 4 groups (3 ducks/group): BSD1 and WSD1 from population 1 and BSD2 and WSD2 from population 2. Twelve libraries from shell glands were sequenced using the Illumina RNA-seq platform, generating an average of 41 million clean reads per library, of which 55.9% were mapped to the duck reference genome and assembled into 31,542 transcripts. Expression levels of 11,698 genes were successfully compared between all pairs of 4 groups. Of these, 464 candidate genes were differentially expressed between cross-phenotype groups, but not for between same-phenotype groups. Gene Ontology (GO) annotation showed that 390 candidate genes were annotated with 2234 GO terms. No candidate genes were directly involved in biosynthesis or transport of biliverdin. However, the integral components of membrane, metal ion transport, cholesterol biosynthesis, signal transduction, skeletal system development, and chemotaxis were significantly (P < 0.05) overrepresented by candidate genes.
This study identified 464 candidate genes associated with duck BGEC, providing valuable information for a better understanding of the mechanisms underlying this trait. Given the involvement of membrane cholesterol contents, ions and ATP levels in modulating the transport activity of bile pigment transporters, the data suggest a potential association between duck BGEC and the transport activity of the related transporters.
From the blue eggs of the whinchat (Saxicola rubetra) to the red eggs of the black-capped donacobius (Donacobius atricapilla), from the speckled eggs of the Japanese quail (Coturnix japonica) to the eggs of the common murre (Uria aalge) with the calligraphic lines, bird eggs display an enormous diversity in eggshell color and maculation pattern, fascinating biological, evolutionary, genetic, and ecological researchers. As a feature exclusively appearing in birds among vertebrates , eggshell color can function in camouflage , reinforcement of eggshell structure , filtering of solar radiation , thermal protection , or indication of female quality , playing important biological roles in ensuring birds successful reproduction. However, the molecular mechanisms underlying these diverse eggshell colors are a poorly understood aspect in avian biology.
Blue-green eggshell color (BGEC) was confirmed to have existed in the extinct upland moa (Megalapteryx didinus)  and in an oviraptorosaur species , suggesting that the color is evolutionarily conserved throughout diverse lineages of Aves . Based on the antioxidant properties of biliverdin, the pigment responsible for BGEC, sexual selection theory indicates that its deposition is related to female antioxidant capacity and egg fertility, inducing males to allocate more parental care in their offspring . There are multiple avian species that lay BGEC eggs, including the whinchat (Saxicola rubetra) , eastern bluebird (Sialia sialis) , spotless starling (Sturnus unicolor) , and some domestic fowls [12–14]. Among these avian species, chickens are unique in that the mutation controlling BGEC was identified as a retrovirus insertion in the SLCO1B3 gene . Ducks are another type of domestic fowl that lay BGEC eggs. In an earlier study, we have examined expression of SLCO1B3 in shell glands of blue-shelled ducks and sequence variants. But, we detected neither expression of SLCO1B3 nor the retrovirus insertion, suggesting that duck BGEC is controlled by different mechanisms . Further identifying the molecular basis of duck BGEC will advance our understanding to the evolution of the trait in the avian phylogeny, and will also improve breeding for the eggshell color in the ducks.
Eggshell color is caused by eggshell pigments. In the final phase of egg formation, the egg enters the shell gland (the uterine part of the oviduct) and stays there for 18–20 h to form eggshell, during which eggshell pigments are secreted into the uterine fluid and progressively deposited in the eggshell, with the deposition rate accelerating in the last 3–5 h before oviposition [16, 17]. Protoporphyrin and biliverdin are two main eggshell pigments and are responsible for red-brown and blue-green eggshell colors, respectively . BGEC is caused predominately by biliverdin-IX and its zinc chelate . Although biliverdin-IX and protoporphyrin simultaneously appear in eggshells of various color, they have different distribution patterns. Protoporphyrin-IX is largely distributed in the cuticle of the eggshell; with heavy scrubbing, the pigment can be removed exposing a white eggshell . In contrast, biliverdin-IX is distributed throughout the eggshell, with the inside of the eggshell being as blue as the outside . This distribution pattern implies that biliverdin deposition should take place together with the calcification of the eggshell.
Biliverdin is a bile pigment that is derived from the oxidative degradation of heme in the reticuloendothelial cells of spleen, liver, bone marrow, and, to a smaller extent, kidney [19, 20]. In mammals, biliverdin is subsequently reduced to bilirubin by biliverdin reductase; the latter is taken up from the blood by hepatocytes, further secreted into bile, and finally excreted in vitro via urine and feces . However, in birds, biliverdin is the predominant end product of heme metabolism due to low biliverdin reductase activity . In terms of the origin of the biliverdin appearing in the eggshell, some researchers believe that it derives from the shell gland, where the pigment is synthesized and secreted [23, 24].
Kennedy and Vevers proposed that BGEC was caused by a specific enzymatic system present in the shell glands of birds laying BGEC eggs that is absent from other birds . At present, it is well known that heme oxygenase (HO) is the rate-limiting enzyme catalyzing the degradation of heme into biliverdin, suggesting that HO is an important candidate for BGEC . In humans and other mammals, HO consists of two active isozymes termed HO-1 and HO-2, which are the products of two different genes, HMOX1 and HMOX2, respectively . HO-1 is an inducible form of HO, and its expression can be induced by all kinds of stimuli and agents that cause oxidative stress and pathological conditions . In contrast, HMOX2 is constitutively expressed at high levels in the testis and brain . Wang et al. found that HMOX1 is highly expressed in the shell glands of blue-shelled chickens compared to that in brown-shelled chickens, which suggests that upregulation of HMOX1 was associated with chicken BGEC . However, it remains unclear whether there is a similar association beween HO expression and duck BGEC.
Apart from this synthesis mechanism, Liu et al. speculated that duck BGEC should be associated with a transport mechanism in the shell gland that allows biliverdin to be more efficiently transported into the uterine fluid, causing BGEC . A similar transport mechanism has been suggested to be present in the chicken . Wang et al. found that chicken BGEC was caused by a retrovirus insertion which activated expression of SLCO1B3 exclusively in the shell glands of blue-shelled chickens . Functionally, OATP1B3 (gene product) is a membrane transporter that mediates the uptake of bile pigments into the hepatocytes . In birds, biliverdin transport mechanisms are poorly understood to date. In the bile pigment transport system of the liver, it is well known that bilirubin is taken up by hepatocytes via organic anion transporting polypeptide (OATP) and excreted into the bile or blood plasma by multidrug resistance protein (MRP) . OATP1A2, OATP1B3, OATP2B1, MRP1 (ABCC1), MRP2 (ABCC2), and MRP3 (ABCC3) are several known members mediating bilirubin transport [28, 29]. Given their highly similar molecular structures and physicochemical properties , biliverdin and bilirubin may share these transporters.
RNA-seq that can sequence the complete set of transcripts in a cell provides an effective method for identifying numerous phenotype-associated differentially expressed genes (DEGs). Using an Illumina RNA-seq platform, a total of 12 libraries (n = 3 per group) from the shell glands of two blue-shelled duck (BSD) and two white-shelled duck (WSD) groups were sequenced. Transcriptomic profiles were compared between all pairs of four groups. This study enables identification of BGEC-associated DEGs to establish a sound foundation for analyses into the mechanisms underlying duck BGEC. In addition, massive high quality reads and transcript assemblies reported here provide valuable resources for further annotation of the duck draft genome .
RNA-seq data and assembly of the transcriptome
Summary of RNA-seq quality, read counts, mapping rates and transcript assemblies for 12 libraries
Clean read number
Number of clean bases (Gb)
Percentage of Q20 bases (%)
Read mapping rate (%)
Unique hits to reference genome (%)
Number of transcripts
Percentage of annotated transcripts (%)
Mean length (bp) of transcripts
Mean ± SD
40,979,283 ± 10,406,949
5.73 ± 1.93
96.3 ± 0.6
55.9 ± 3.3
96.9 ± 1.5
31,542 ± 5308
71.3 ± 3.1
1839 ± 86
2849 ± 142
Identification of BGEC-associated candidate genes
Functional annotation and overrepresentation analysis of candidate genes
Expression analysis of ion transporters related to eggshell mineralization
During eggshell formation, a series of ions, including Ca2+, HCO3 −, Na+, K+, H+, and Cl− participate in the mineralization process . Here, a total of 13 Ca2+, 10 Na+, 8 K+, 5 H+, and 1 Cl− transporter and two HCO3 − synthesis genes were identified among the candidate genes (Additional files 2, and 4). Among the 13 Ca2+ transport genes, ATP2C2, HTP1B, PKD2, RASA3, CACNB2, TRPC1, TRPC4, and TRPV2 were further annotated as calcium ion import (GO:0070509) and positive regulation of calcium ion import (GO:0090280) (Additional file 5). These Ca2+ import genes were almost all downregulated in the BSD groups, with the exception of HTR1B and ATP2C2 (Fig. 5). In contrast, three Ca2+ pumps (ATP2A2, ATP2B2, and ATP2C2) were all upregulated (Fig. 5). Genes mediating Na+, K+, H+, and Cl− transport were almost all upregulated, with the exception of one Na+ (PKD2) and two K+ (PKD2 and KCNF1) transport genes (Fig. 5). In addition, two HCO3 − synthesis genes (CA2 and CA8) were upregulated in the BSD groups (Fig. 5).
We identified a considerable number of candidate genes involved in ion transport coupled ATP hydrolysis (i.e. Ca2+ pumps and Na+/K+ exchangers) and synthesis (i.e. H+ transporters). We found five ion transporters (ATP12, ATP1B1, ATP2A2, ATP2B2, and ATP2C2) with ATPase coupled ion transmembrane transporter activity (GO:0042625), one H+ transporter (ATP5J) with ATP synthesis coupled proton transport (GO:0015986), and two candidate genes (ENSAPLG00000000921 and CYC) with mitochondrial ATP synthesis coupled electron transport (GO:0042775). The above ATP metabolism genes were all upregulated in the BSD groups (Fig. 5).
Verification of DEGs by qPCR
RNA-seq data were verified by qPCR. A total of nine genes were selected for verification, representing three expression types: (1) SLCO1A2, SLCO1B3, and ABCC2, which were only barely detected by RNA-seq (Fig. 3); (2) HMOX1 and ABCC1, which were expressed, but did not show a significant expression difference or a consistent trend among four cross-phenotype comparisons (Fig. 3); (3) CA2, ATP1B1, ATP2B2 and SLC14A1, four candidate genes that were all upregulated in the BSD groups (Additional file 2). Expression of SLCO1A2, SLCO1B3, and ABCC2 remained undetected by qPCR (Additional file 6). Expression of the other seven genes was consistent with the RNA-seq results (Additional file 6). The correlation coefficient between qPCR and RNA-seq was 0.914 (P < 0.01), indicating good reproducibility of the differential expression results obtained by RNA-seq (Additional file 6).
Eggshell color, a principle component of avian reproductive system, is related to reproduction in wild birds , and to customer preference [41, 42] and egg nutrition [43, 44] in domestic fowls. Although considerable research efforts have focused on the biological roles of eggshell color and the biochemistry of eggshell pigments [17, 24, 45], molecular mechanisms underlying eggshell colors are poorly understood so far. Using the Illumina RNA-seq, we obtained a global view of transcriptomes from shell glands of BSD and WSD. Further, expression levels of 11,698 genes were successfully compared between all pairs of two BSD and two WSD groups, resulting in the discoveries of 464 BGEC-associated candidate genes. This study is among the first to employ RNA-seq to robustly discover BGEC-associated DEGs, further to provide insight into the molecular mechanisms.
BGEC is caused by biliverdin . Thus, any mechanisms that enhance the biosynthesis or transport efficiency of the pigment in the shell gland can cause BGEC. Here, the expression of two biliverdin biosynthesis genes, HMOX1 and HMOX2, was not significantly different between BSD and WSD groups. This result is compatible with the finding of Liu et al. that BSD and WSD have similar enzymatic activities of HO in the shell glands , suggesting that duck BGEC is not caused by a biliverdin synthesis mechanism.
In all kingdoms of life, anomalous transport of pigments themselves or pigment-related substrates is another important pigmentation mechanism, which has been confirmed to be responsible for eye color in D.melanogaster  and silkworm ; shell color in pacific oyster ; and coat color in horse , tiger , and chicken . A pigment transport mechanism responsible for bird BGEC has been proposed by Liu et al.  and Wang et al. . OATP1A2, OATP1B3, OATP2B1, MRP1, MRP2 and MRP3 are well-known transporters mediating bile pigment transport [28, 29]. It is therefore tempting to speculate that upregulation of these transporters increase transport efficiency of biliverdin, further causing duck BGEC. But, current results did not obviously support the speculation (Fig. 3). However, there is another possibility that it is due to transport activity rather than to upregulation of transporters that duck BGEC is caused.
Cholesterol contents in the membrane have a great effect on the activities of bile salt transporters. A well-characterized example is derived from the association of intrahepatic cholestasis with the mutations in ATP8B1 . ATP8B1 is a member of P4 ATPase family and exerts an important role in resisting bile salt-mediated cholesterol extraction from the canalicular membrane . ATP8B1 deficiency reduces cholesterol contents in the canalicular membrane, which impairs the activities of the bile salt export pumps, further causing the cholestasis . In this study, we found that the steroid biosynthsis process was highly enriched by 9 candidate genes which were almost all upregulated in the BSD except for SRD5A2 (Fig. 5). ABCA1 that mediates cholesterol export in peripheral tissues was downregulated in the BSD (Fig. 5) . Although the exact interaction mechanisms of cholesterol and biliverdin transporters have not been fully understood, our results and the known relationship between membrane cholesterol contents and the transport activity of bile salt transporters urge us to speculate that upregulation of cholesterol biosynthesis genes and downregulation of a cholesterol export gene (ABCA1) may exert a positive action for membrane cholesterol contents, which promotes the activity of biliverdin transporters logically.
Apart from the potential cholesterol-related mechanism, ion transport can be implicated into the formation of BGEC. Eggshell mineralization refers to transmembrane transfers of multiple ions including Ca2+, HCO3 −, Na+, K+, H+, and Cl− . During the process, concentrations of Ca2+, K+ and H+ in the uterine fluid continuously increase from 8 h to 18 h post ovulation; ones of Na+ and Cl− generally decline . A similar increasing process has been observed for biliverdin from 12 h to 23.5 h post ovulation . The synchronism implies that an ion/biliverdin exchange or symport mechanism can exist in the shell glands of the birds laying BGEC eggs, by which biliverdin deposition can occur together with eggshell mineralization.
Firstly, ion transporters have been suggested to be the promising candidates for the eggshell mechanical properties . In this study, although a considerable number of ion transporters were differentially expressed between BSD and WSD, there is no evidence supporting that blue and white duck eggs have a significant difference in eggshell quality. The absence of difference between blue and white eggs has been confirmed in the chicken .
Secondly, CA2 and CA8 were upregulated in the BSD (Fig. 5), which can increase the yield of intracellular HCO3 −. Redundant HCO3 − can be transported into the blood plasma via OATP2B1 (Fig. 7). The consequence of balancing the concentration of intracellular HCO3 − is that more biliverdin is exchanged into glandular cells (Fig. 7). In addition, given that two Na+/H+ exchangers were upregulated, an acidic extracellular pH may further enhance the transport activity of OATP2B1 (Fig. 7).
Thirdly, there were three Ca2+ (ATP2A2/2B2/2C2) and two Na+ (ATP1B1, ATP12A) pump, three ATP synthetic (ATP5J, ENSAPLG00000000921 and CYC) and one ATP transport (SLC25A24) genes which were upregulated in the BSD among 464 candidate genes (Fig. 5). To maintain intracellular ion and ATP homeostasis, potential effect caused by upregulation of these genes can be balanced by a competition for ATP between ion pumps and MRPs. The consequence of the competition may be more biliverdin transferred into the uterine fluid, finally causing the BGEC (Fig. 7).
As we know, administration of coccidiostat Nicarbazin leads to eggshell depigmentation , and acetazolamide suppresses the activity of carbonic anhydrases . Therefore, it will be helpful to accurately characterize the association of these candidate genes with duck BGEC, and to establish their relative contribution by incorporating the inhibitory models in the future study.
This study found 464 candidate genes potentially associated with duck BGEC whereas none of them were directly implicated into biliverdin biosynthesis or transport. In contrast, cholesterol biosynthsis and ion transport processes were significantly overrepresented by candidate genes. Given potential roles of membrane cholesterol contents, ions and intracellular ATP levels in modulating the transport activity of bile pigment transporters, the data suggests that duck BGEC might be associated with the change of the transport activity of the related transporters. Our present and earlier  results support the possibility that the abnormal transport of biliverdin may be shared by avian species laying BGEC eggs despite of the causative mutations probably differing from each other. Sexual selection theory interpreted BGEC as a signal reflecting the genetic quality of the females . Our balance hypothesis speculates that BGEC is a consequence of balancing potential effect caused by differential expression of ion transport and synthesis genes, supporting that a structural significance can exist for BGEC.
Even though the currently available duck genome annotation has been shown to be quite incomplete , we are thus sure to have implemented the best approach we could. A new, much improved duck genome reference with annotation is underway by international efforts but will only become available in the future. Our data of 492 million high quality reads and assembled transcripts will aid significantly in the annotation efforts by enriching the databases during future annotations.
Ducks and collection of shell glands
A Chinese indigenous breed, Jinding duck (Anas platyrhynchos), was used in this study. Jinding duck is originated from the Fujiang province of China, and characterized by an excellent laying performance . In the Jinding ducks, the eggshell color is not fixed, where major individuals lay BGEC eggs, but there are minor ducks laying white eggs. Two Jinding duck populations were independently raised in two duck farms at the Baiyang Lake of HeBei province of China. There was no gene flow between the two duck populations. Three blue-shelled and three white-shelled ducks were randomly selected from population 1, and the other three blue-shelled and three white-shelled ducks were collected from population 2. In the text, blue- and white-shelled ducks from population 1 were respectively defined as BSD1 and WSD1, and ducks from population 2 were named as BSD2 and WSD2. All birds were 35 weeks old and had access to food and water ad libitum.
Before sampling, a part of the ducks were caged individually. The oviposition was monitored every 30 min. According to the oviposition routines, ducks were stunned and killed at 3–5 h before estimated oviposition. Shell glands (the uterine part of the oviduct) were excised, rinsed briefly with 0.9% isotonic saline and immediately immersed in RNAstore RNA stabilization reagent (CWBIO) overnight at 4 °C and finally stored at −80 °C. The care and use of ducks and collection of shell glands were performed according to the guidelines of the Animal Care and Use Committee of Northwest A&F University.
RNA-Seq library preparation and Illumina-Solexa sequencing
Total RNA was extracted from shell glands with TRIzol reagent (Invitrogen) according to the manufacture’s protocol. The sequencing libraries were prepared for each duck using NEBNext® Ultra™ RNA Library Prep Kit (NEB). Briefly, mRNA was isolated from total RNA with poly-T oligo-attached magnetic beads and then fragmented. First-strand cDNA was synthesized using M-MuLV reverse transcriptase (RNase H-) and random primers. The cDNA was further converted into double stranded cDNA. Following the reverse transcription, the cDNA was end-repaired, adenylated, and ligated to NEBNext Adaptors. The cDNA fragments of 150–200 bp in length were selected with AMPure XP system (Beckman Coulter). The size-selected, adaptor-ligated cDNA was treated with 3 μL USER enzyme (NEB). PCR was performed with Phusion High-Fidelity DNA polymerase (NEB), universal PCR primers and Index (X) Primer. Clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer’s instructions. After the clustering was completed, these libraries were sequenced using a paired-end 2 × 125 bp lane (BSD1 and WSD1) and a 2 × 150 bp lane (BSD2 and WSD2) on an Illumina HiSeq 2500 platform.
Acquisition and analysis of RNA-Seq data
A total of 12 libraries (n = 3 per group) from shell glands of two BSD and two WSD groups were sequenced. Raw reads of fastq format were firstly filtered to generate clean reads by removing reads containing adapters or ambiguous nucleotides (the proportion of “N” base exceeding 10%), and reads of low quality (the proportion of bases with Phred quality score of less than 20 exceeding 50%). The filtering steps were completed by Novogene Co. (Beijing, China) using their in-house Perl scripts. The duck reference genome (BGI duck 1.0 Assembly, release-82) and gene model annotation files (GTF file) were downloaded from Ensembl database [59, 60]. Indexes of the reference genome were built using Bowtie v2.2.7 . The paired-end clean reads were aligned to the reference genome by Tophat v2.0.14 using the parameters described by Tropnell C et al. with little modification made in the parameter of N that was set to four .
Following alignment, mapped reads were handled by Cufflinks v2.2.0 to generate a transcriptome assembly for each library. All assemblies from 12 libraries and reference transcripts were merged by Cuffmerge to generate a uniform basis for calculating transcript and gene expression . Finally, the mapped reads, the merged assembly, and the duck reference genome were submitted to Cuffdiff to calculate expression levels, run bias detection and correction algorithm, and test the statistical significance of changes in abundance of each gene between all pairs of 4 groups .
Cuffdiff estimates abundance at the transcript level using a linear statistical model . The abundance of each gene is calculated by adding up the expression levels of a group of transcripts from the gene . Abundance is reported as fragments per kilobase of transcript per million mapped fragments (FPKM), which is calculated according to the number of reads uniquely mapped to the gene . In this method, read counts are normalized by gene length, and then by sequencing depth, which ensures that abundance estimates are proportional to the number of reads produced by each gene and are less influenced by these variables [62, 63].
When more than two groups are handled, Cuffdiff performs all possible pair-wise comparisons between these groups . To test for significant differences, the square root of the Jensen-Shannon divergence of the relative abundances in each of two groups is calculated. The variance of this Jensen-Shannon distance under the null hypothesis of no change in relative abundance is estimated. Using the estimated variance, one-sided t-test is performed to evaluate the significance in the abundance changes of each gene between two groups . The P-value was adjusted using the Benjamini and Hochberg’s approach for controlling type I errors due to multiple testing . Genes with a q value <0.05 were determined to be differentially expressed.
Screening of candidate genes and functional annotation
Candidate genes were obtained according to three conditions as follows: (1) an assembled gene was annotated to unique one reference gene; (2) gene expression was detected among the four duck groups, and differential expression analysis was successfully performed across all six pair-wise comparisons; (3) expression levels of genes were consistently significantly different between all cross-phenotype groups (BSD1 vs. WSD1, BSD1 vs. WSD2, BSD2 vs. WSD1 and BSD2 vs. WSD2), but not for between same-phenotype groups (BSD1 vs. BSD2 and WSD1 vs. WSD2).
Functional annotation and overrepresentation tests of candidate genes were completed using the PANTHER Classification System . Complete GO molecular function, biological process, and cellular component categories were used as function annotation and overrepresentation analysis data sets. Because the duck genes have not been deposited in the PANTHER Classification System to date, chicken orthologues of these candidate genes (Additional file 7) were submitted to the PANTHER for functional analysis. Chicken orthologues were found by BLAT alignment between duck genes and chicken genome (Galgal4, Ensembl release 82), with a cut-off E-value of 0.00001 . The genes with the highest alignment score were identified as the chicken orthologues. Overrepresentation tests against the chicken genome background were performed using the statistical overrepresentation test option. Statistical significance was assessed using the binomial test. Bonferroni correction was performed to control type I errors caused by multiple tests. GO terms with an adjusted P value less than 0.05 were considered to be significantly overrepresented. If there was an ancestor-child hierarchy relationship among significant GO terms, only the most specific child terms were considered.
Real-time quantitative RT-PCR (qPCR) verification of RNA-seq data
The RNA-seq data were verified by qPCR using the same RNA samples used for RNA-seq. The concentration of total RNA was measured using a NanoDrop 2000c spectrophotometer (Thermo Fisher Scientific). One μg of RNA was used as a template from which first-strand cDNA was synthesized using a SuperRT cDNA Synthesis Kit (CWBIO). QPCR was performed using an UltraSYBR Mixture Kit (CWBIO). The total volume of the qPCR mixture was 20 μL containing 10 μL of 2 × UltraSYBR Mixture (High ROX), 0.4 μL of forward primer (10 μM), 0.4 μL of reverse primer (10 μM), 1 μL of cDNA, and 8.2 μL of ddH2O. Reactions were carried out on a Bio-Rad iQ™5 Thermal Cycler (Bio-Rad) with the following reaction conditions: denaturation at 95 °C for 10 min; 40 cycles of denaturation at 95 °C for 15 s, annealing and elongation for 1 min at 60 °C, and a melting curve process from 55 to 95 °C. Nine genes were selected as qPCR-verified targets. The primer sequences for these genes were listed in Additional file 8. Three biological replicates per groups and two technological replicates per sample were included in the qPCR. GAPDH was used as the endogenous reference gene, and the WSD2 group was set as the criterion. Expression levels of genes in each group were presented as expression folds relative to the criterion using the 2-ΔΔCT method .
This study was supported by the National Natural Science Foundation of China (Grant No. 31401051) and China Postdoctoral Science Foundation (2014M550510 and 2015T81060).
Availability of data and materials
Clean data from transcriptome sequencing has been deposited in the NCBI Sequence Read Archive database (https://trace.ncbi.nlm.nih.gov/Traces/sra/) with accession number: SRP097820.
ZPW designed the research framework, carried out the experiments, performed the bioinformatics analyses, and drafted the manuscript. GHM participated in data analysis. YB and GHM performed the qPCR verification experiment. GHM and YD collected and prepared the experimental materials. RFL contributed to the research design and reviewed the manuscript. LHS participated in manuscript formatting and editing. All authors read and approved the final manuscript.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interest.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Cassey P, Portugal SJ, Maurer G, Ewen JG, Boulton RL, Hauber ME, Blackburn TM. Variability in avian eggshell colour: a comparative study of museum eggshells. PLoS One. 2010;5:e12054.View ArticlePubMedPubMed CentralGoogle Scholar
- Skrade PDB, Dinsmore SJ. Egg crypsis in a ground-nesting shorebird influences nest survival. Ecosphere. 2013;4:151.View ArticleGoogle Scholar
- Gosler AG, Higham JP, Reynolds SJ. Why are birds’ eggs speckled? Ecol Lett. 2005;8:1105–13.View ArticleGoogle Scholar
- Lahti DC. Population differentiation and rapid evolution of egg color in accordance with solar radiation. Auk. 2008;125:796–802.View ArticleGoogle Scholar
- Bakken GS, Vanderbilt VC, Buttemer WA, Dawson WR. Avian eggs: thermoregulatory value of very high near infra-red reflectance. Science. 1978;200:321–3.View ArticlePubMedGoogle Scholar
- Moreno J, Osorno JL. Avian egg colour and sexual selection: does eggshell pigmentation reflect female condition and genetic quality? Ecol Lett. 2003;6:803–6.View ArticleGoogle Scholar
- Igic B, Greenwood DR, Palmer DJ, Cassey P, Gill BJ, Grim T, Brennan PLR, Bassett SM, Battley PF, Hauber ME. Detecting pigments from colourful eggshells of extinct birds. Chemoecology. 2010;20:43–8.View ArticleGoogle Scholar
- Wiemann J, Yang T, Sander PNN, Schneider M, Engeser M, Kath-Schorr S, Müller CE, Sander M. The blue-green eggs of dinosaurs: how fossil metabolites provide insights into the evolution of bird reproduction. Peer J Pre Prints. 2015;3:e1080v1. https://doi.org/10.7287/peerj.preprints.1080v1 Google Scholar
- Brulez K, Mikšík I, Cooney CR, Hauber ME, Lovell PG, Maurer G, Portugal SJ, Russell D, Reynolds SJ, Cassey P. Eggshell pigment composition covaries with phylogeny but not with life history or with nesting ecology traits of British passerines. Ecol Evol. 2016;6:1637–45.View ArticlePubMedPubMed CentralGoogle Scholar
- Siefferman L, Navara KJ, Hill GE. Egg coloration is correlated with female condition in eastern bluebirds (Sialia sialis). Behav Ecol Sociobiol. 2006;59:651–6.View ArticleGoogle Scholar
- López-Rull I, Miksik I, Gil D. Egg pigmentation reflects female and egg quality in the spotless starling Sturnus unicolor. Behav Ecol Sociobiol. 2008;62:1877–84.View ArticleGoogle Scholar
- Punnett RC. Genetic study in poultry: IX. The blue egg. J Genet. 1933;27:465–70.View ArticleGoogle Scholar
- Liu HC, Hsiao MC, Hu YH, Lee SR, Cheng WTK. Eggshell pigmentation study in blue-shelled and white-shelled ducks. Asian-Australas J Anim Sci. 2010;23:162–8.View ArticleGoogle Scholar
- Ito S, Tsudzuki M, Komori M, Mizutani M. Celadon: an eggshell color mutation in Japanese quail. J Hered. 1993;84:145–7.View ArticlePubMedGoogle Scholar
- Wang ZP, Qu LJ, Yao JF, Yang XL, Li GQ, Zhang YY, Li JY, Wang XT, Bai JR, Xu GY, et al. An EAV-HP insertion in 5′ flanking region of SLCO1B3 causes blue eggshell in the chicken. PLoS Genet. 2013;9:e1003183.View ArticlePubMedPubMed CentralGoogle Scholar
- Jonchère V, Réhault-Godbert S, Hennequet-Antier C, Cabau C, Sibut V, Cogburn LA, Nys Y, Gautron J. Gene expression profiling to identify eggshell proteins involved in physical defense of the chicken egg. BMC Genomics. 2010;11:57.View ArticlePubMedPubMed CentralGoogle Scholar
- Lang MR, Wells JW. A review of eggshell pigmentation. Worlds Poult Sci J. 1987;43:238–45.View ArticleGoogle Scholar
- Steggerda M, Hollander WF. Observations on certain shell variations of hens’ eggs. Poult Sci. 1944;23:459–61.View ArticleGoogle Scholar
- Maines MD. The heme oxygenase system: a regulator of second messenger gases. Annu Rev Pharmacol Toxical. 1997;37:517–54.View ArticleGoogle Scholar
- Hudson MF, Smith KM. Bile pigments. Chem Soc Rev. 1975;4:363–99.View ArticleGoogle Scholar
- Wang X, Chowdhury JR, Chowdhury NR. Bilirubin metabolism: applied physiology. Curr Paediatr. 2006;16:70–4.View ArticleGoogle Scholar
- Lin GL, Himes JA, Cornelius CE. Bilirubin and biliverdin excretion by the chicken. Am J Physiol. 1974;226:881–5.PubMedGoogle Scholar
- Zhao R, Xu GY, Liu ZZ, Li XY, Yang N. A study on eggshell pigmentation: biliverdin in blue-shelled chickens. Poult Sci. 2006;85:546–9.View ArticlePubMedGoogle Scholar
- Kennedy GY, Vevers HG. Eggshell pigments of the araucano fowl. Comp Biochem Physiol B. 1973;44B:11–25.View ArticleGoogle Scholar
- Schuller DJ, Wilks A, Ortiz de Montellano PR, Poulos TL. Crystal structure of human heme oxygenase-1. Nat Struct Biol. 1999;6:860–7.View ArticlePubMedGoogle Scholar
- Wang ZP, Liu RF, Wang AR, Li JY, Deng XM. Expression and activity analysis reveal that heme oxygenase (decycling) 1 is associated with blue egg formation. Poult Sci. 2011;90:836–41.View ArticlePubMedGoogle Scholar
- Stieger B, Hagenbuch B. Organic anion transporting polypeptides. Curr Top Membr. 2014;73:205–32.View ArticlePubMedPubMed CentralGoogle Scholar
- Keppler D. The roles of MRP2, MRP3, OATP1B1, and OATP1B3 in conjugated hyperbilirubinemia. Drug Metab Dispos. 2014;42:561–5.View ArticlePubMedGoogle Scholar
- Pascolo L, Fernetti C, Garcia-Mediavilla MV, Ostrow JD, Tiribelli C. Mechanisms for the transport of unconjugated bilirubin in human trophoblastic BeWo cells. FEBS Lett. 2001;495:94–9.View ArticlePubMedGoogle Scholar
- Kaur H, Hughes MN, Green CJ, Naughton P, Foresti R, Motterlini R. Interaction of bilirubin and biliverdin with reactive nitrogen species. FEBS Lett. 2003;543:113–9.View ArticlePubMedGoogle Scholar
- Huang Y, Li Y, Burt DW, Chen H, Zhang Y, Qian W, Kim H, Gan S, Zhao Y, Li J, et al. The duck genome and transcriptome provide insight into an avian influenza virus reservoir species. Nat Genet. 2013;45:776–83.View ArticlePubMedPubMed CentralGoogle Scholar
- Truong AD, Hong YH, Lillehoj HS. RNA-seq profiles of immune related genes in the spleen of necrotic enteritis-afflicted chicken lines. Asian Australas J Anim Sci. 2015;28:1496–511.View ArticlePubMedPubMed CentralGoogle Scholar
- Rowley JW, Oler AJ, Tolley ND, Hunter BN, Low EN, Nix DA, Yost CC, Zimmerman GA, Weyrich AS. Genome-wide RNA-seq analysis of human and mouse platelet transcriptomes. Blood. 2011;118:e101–11.View ArticlePubMedPubMed CentralGoogle Scholar
- Zhu F, Yuan JM, Zhang ZH, Hao JP, Yang YZ, Hu SQ, Yang FX, Qu LJ, Hou ZC. De novo transcriptome assembly and identification of genes associated with feed conversion ratio and breast muscle yield in domestic ducks. Anim Genet. 2015;46:636–45.View ArticlePubMedGoogle Scholar
- Martin JA, Wang Z. Next-generation transcriptome assembly. Nat Rev Genet. 2011;12:671–82.View ArticlePubMedGoogle Scholar
- Conesa A, Madrigal P, Tarazona S, Gomez-Cabrero D, Cervera A, McPherson A, Szcześniak MW, Gaffney DJ, Elo LL, Zhang X, et al. A survey of best practices for RNA-seq data analysis. Genome Biol. 2016;17:13.View ArticlePubMedPubMed CentralGoogle Scholar
- Nys Y, Gautron J, Garcia-Ruiz JM, Hincke MT. Avian eggshell mineralization: biochemical and functional characterization of matrix proteins. CR Palevol. 2004;3:549–62.View ArticleGoogle Scholar
- Mi HY, Muruganujan A, Casagrande JT, Tomas PD. Large-scale gene function analysis with the PANTHER classification system. Nat Protoc. 2013;8:1551–66.View ArticlePubMedGoogle Scholar
- Vit O, Petrak J. Integral membrane proteins in proteomics. How to break open the black box? J Proteome. 2017;153:8–20.View ArticleGoogle Scholar
- Jonchère V, Brionne A, Gautron J, Nys Y. Identification of uterine ion transporters for mineralization precursors of the avian eggshell. BMC Physiol. 2012;12:10.View ArticlePubMedPubMed CentralGoogle Scholar
- Liu HC, Cheng WT. Eggshell pigmentation: a review. J Chin Soc Anim Sci. 2010;39:75–89.Google Scholar
- Odabaşi AZ, Miles RD, Balaban MO, Portier KM. Changes in brown eggshell color as the hen ages. Poult Sci. 2007;86:356–63.PubMedGoogle Scholar
- Butler MW, McGraw KJ. Eggshell coloration reflects both yolk characteristics and dietary carotenoid history of female mallards. Funct Ecol. 2013;27:1176–85.View ArticleGoogle Scholar
- Wang XL, Zheng JX, Ning ZH, Qu LJ, Xu GY, Yang N. Laying performance and egg quality of blue-shelled layers as affected by different housing systems. Poult Sci. 2009;88:1485–92.View ArticlePubMedGoogle Scholar
- Reynolds SJ, Martin GR, Cassey P. Is sexual selection blurring the functional significance of eggshell coloration hypotheses? Anim Behav. 2009;78:209–15.View ArticleGoogle Scholar
- Dermauw W, Van Leeuwen T. The ABC gene family in arthropods: comparative genomics and role in insecticide transport and resistance. Insect Biochem Mol Biol. 2014;45:89–110.View ArticlePubMedGoogle Scholar
- Kômoto N, Quan GX, Sezutsu H, Tamura T. A single-base deletion in an ABC transporter gene causes white eyes, white eggs, and translucent larval skin in the silkworm w-3(oe) mutant. Insect Biochem Mol Biol. 2009;39:152–6.View ArticlePubMedGoogle Scholar
- Feng D, Li Q, Yu H, Zhao X, Kong L. Comparative transcriptome analysis of the pacific oyster crassostrea gigas characterized by shell colors: identification of genetic bases potentially involved in pigmentation. PLoS One. 2015;10:e0145257.View ArticlePubMedPubMed CentralGoogle Scholar
- Mariat D, Taourit S, Guérin G. A mutation in the MATP gene causes the cream coat colour in the horse. Genet Sel Evol. 2003;35:119–33.View ArticlePubMedPubMed CentralGoogle Scholar
- Xu X, Dong GX, Hu XS, Miao L, Zhang XL, Zhang DL, Yang HD, Zhang TY, Zou ZT, Zhang TT, et al. The genetic basis of white tigers. Curr Biol. 2013;23:1031–5.View ArticlePubMedGoogle Scholar
- Gunnarsson U, Hellström AR, Tixier-Boichard M, Minvielle F, Bed'hom B, Ito S, Jensen P, Rattink A, Vereijken A, Andersson L. Mutations in SLC45A2 cause plumage color variation in chicken and Japanese quail. Genetics. 2007;175:867–77.View ArticlePubMedPubMed CentralGoogle Scholar
- Folmer DE, Elferink RP, Paulusma CC. P4 ATPases - lipid flippases and their role in disease. Biochim Biophys Acta. 1791;2009:628–35.Google Scholar
- Gelissen IC, Harris M, Rye KA, Quinn C, Brown AJ, Kockx M, Cartland S, Packianathan M, Kritharides L, Jessup W. ABCA1 And ABCG1 synergize to mediate cholesterol export to apoA-I. Arterioscler Thromb Vasc Biol. 2006;26:534–40.View ArticlePubMedGoogle Scholar
- Nys Y, Hincke MT, Arias JL, Garcia-Ruiz JM, Solomon SE. Avian eggshell mineralization. Avian Biol Res. 1999;10:143–66.Google Scholar
- Duan Z, Chen S, Sun C, Shi F, Wu G, Liu A, Xu G, Yang N. Polymorphisms in ion transport genes are associated with eggshell mechanical property. PLoS One. 2015;10:e0130160.View ArticlePubMedPubMed CentralGoogle Scholar
- Sadjadi M, Renden JA, Benoff FH, Harper JA. Effects of the blue egg shell allele (O) on egg quality and other economic traits in the chicken. Poult Sci. 1983;62:1717–20.View ArticleGoogle Scholar
- Eastin WC, Spaziani E. On the mechanism of calcium secretion in the avian shell gland (uterus). Biol Reprod. 1978;19:505–18.View ArticlePubMedGoogle Scholar
- Hoque MA, Skerratt LF, Rahman MA, Alim MA, Grace D, Gummow B, Rabiul Alam Beg AB, Debnath NC. Monitoring the health and production of household Jinding ducks on Hatia Island of Bangladesh. Trop Anim Health Prod. 2011;43:431–40.View ArticlePubMedGoogle Scholar
- Duck reference genome in the Ensembl database. ftp://ftp.ensembl.org/pub/release-82/fasta/anas_platyrhynchos/dna/. Accessed 20 Apr 2016.
- GTF file for duck reference genome. ftp://ftp.ensembl.org/pub/release-82/gtf/anas_platyrhynchos. Accessed 20 Apr 2016.
- Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10:R25.View ArticlePubMedPubMed CentralGoogle Scholar
- Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and cufflinks. Nat Protoc. 2012;7:562–78.View ArticlePubMedPubMed CentralGoogle Scholar
- Sims D, Sudbery I, Ilott NE, Heger A, Ponting CP. Sequencing depth and coverage: key considerations in genomic analyses. Nat Rev Genet. 2014;15:121–32.View ArticlePubMedGoogle Scholar
- Cufflinks. http://cole-trapnell-lab.github.io/cufflinks/. Accessed 9 Apr 2017.
- Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, Salzberg SL, Wold BJ, Pachter L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28:511–5.View ArticlePubMedPubMed CentralGoogle Scholar
- BLAST/BLAT search. http://asia.ensembl.org/Multi/Tools/Blast?db=core. Accessed 10 Jan 2017.
- Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitivative PCR and the 2-ΔΔCT method. Methods. 2001;25:402–8.View ArticlePubMedGoogle Scholar