Isolation of blue-green eggshell pigmentation-related genes from Putian duck through RNA-seq

The diversity of avian eggshell colour plays important biological roles in ensuring successful reproduction. Eggshell colour is also an important trait in poultry, but the mechanisms underlying it are poorly understood in ducks. This study aimed to provide insights into the mechanism of blue-green eggshell colour generation. Here, white-shelled ducks (HBR) and blue-green-shelled ducks (HQR) were selected from Putian black ducks, and white-shelled ducks (BBR) were selected from Putian white ducks. Transcriptional changes in the shell gland were analysed using RNA-sequencing on the Illumina HiSeq 2500. Twenty-seven individual cDNA libraries were sequenced and generated an average of 7.35 million reads per library; 70.6% were mapped to the duck reference genome, yielding an average of 13,794 genes detected, which accounted for approximately 86.39% of all 15,967 annotated duck genes. A total of 899 differentially expressed genes (DEGs) were detected between the HQR and BBR groups, and 373 DEGs were detected between the HQR and HBR groups. We analysed the DEGs in the HQR-vs-BBR and HQR-vs-HBR comparisons. None of these DEGs were directly involved in the eggshell pigmentation process in HQR-vs-HBR, while UDP-glucuronosyltransferase 2A2 (UGT2A2) and UDP-glucuronosyltransferase 1–1-like (UGT1–1-like), which participate in biliverdin breakdown, were two of the DEGs in HQR-vs-BBR. In the RT-qPCR results, delta-aminolevulinic acid synthase 1 (ALAS1) and EPRS glutamyl-prolyl-tRNA synthetase were significantly upregulated in the HBR group compared with the HQR and BBR groups (P < 0.05). Haem oxygenase (HMOX1) was significantly downregulated in BBR compared with HQR and HBR (P < 0.05). Biliverdin reductase A (BLVRA), GUSB glucuronidase beta, cytochrome c-type haem lyase, protohaem IX farnesyltransferase and UGT2A2 were significantly upregulated in HBR and BBR compared with HQR (P < 0.05). We conducted a comparative transcriptome analysis of the shell glands of Putian white ducks and Putian black ducks. None of the differentially regulated pathways were directly involved in the eggshell pigmentation process in the HQR-vs-HBR comparison, while 2 DEGs related to biliverdin breakdown were found in HQR-vs-BBR. Based on the RT-qPCR results, we can speculate that both HQR and HBR can produce biliverdin, but HBR cannot accumulate it. Compared with HQR, BBR produced less biliverdin and did not accumulate it.


Background
Bird eggshell colour displays enormous diversity and has multiple functions, such as avoiding predation [1], strengthening the structure of the shell [2], filtering harmful solar radiation [3] and decreasing trans-shell microbial contamination [4]. Colourful biological pigments are deposited in the shell gland during eggshell formation, resulting in different eggshell colours. Protoporphyrin IX, biliverdin and Zn-biliverdin chelate are three main pigments responsible for eggshell colouration. Protoporphyrin IX is an immediate precursor of haem, which causes reddish or brown background eggshell colour [5]. Biliverdin is a byproduct of the breakdown of haemoglobin and gives eggshells a blue or green colour [6]. Protoporphyrin is a pro-oxidant [7], while biliverdin possesses antioxidant activity [8]. The sexually selected eggshell colouration (SSEC) hypothesis assumes a positive association between female quality and the amount of biliverdin deposited into the eggshell during reproduction [9]. At the same time, darker and brown eggshells reflect poorer maternal condition [10,11]. Identifying the molecular mechanisms involved in eggshell colouration is key to analysing the selection and evolution of this trait.
White and blue-green are the two major eggshell colour phenotypes in ducks, and eggs with a blue-green colour are more acceptable to consumers in Asia [12]. Improving egg colour in domestic poultry has long been of commercial interest, and the mechanism of eggshell colouration in laying hens has been studied previously. Most researchers have suggested that protoporphyrin IX and biliverdin are first synthesized in the shell gland and then secreted and deposited into the eggshell layers [13][14][15]. Some researchers believe that a specific enzymatic system is present in the shell glands of birds laying blue-green eggshell colour eggs that is absent from other birds [16]. The biosynthetic pathway of eggshell colour-related pigments is well established. Briefly, SLC25A38 (solute carrier family 25, member 38) imports glycine into the mitochondrial matrix from the cytosol [17]. ALAS1 catalyses the rate-limiting step in the condensation of succinyl coenzyme and glycine to delta-aminolevulinic acid (ALA) [18]. Following its synthesis, ALA is exported to the cytosol and converted to coproporphyrinogen III [19]. ABCB6 belongs to the ATP-binding cassette family that transports coproporphyrinogen III back to mitochondria [20]. Coproporphyrinogen oxidase (CPOX) then catalyses the conversion of coproporphyrinogen III to protoporphyrinogen IX [21], and protoporphyrinogen IX is oxidized into protoporphyrin IX via protoporphyrinogen oxidase [21], which gives eggshells a brown or pink colour [5]. Ferrochelatase (FECH) catalyses the terminal step of haem biosynthesis by inserting ferrous ions into protoporphyrin IX [22]. Haem oxygenase (HMOX) participates in the haem degradation pathway, which converts protohaem into biliverdin, which gives eggshells a blue-green colour [6,23,24]. Biliverdin can be converted into bilirubin reversibly [25]. UDP-glucuronosyltransferase catalyses the formation of bilirubin β-diglucuronide [26], which is then converted into D-urobilinogen by GUSB glucuronidase beta. Finally, D-urobilinogen is oxidized into urobilin and stercobilin and then excreted out of the body. However, it remains unclear which key genes are associated with blue-green eggshell colouration in ducks.
Putian white duck and Putian black duck were both bred from Tadorna shelducks. Here, we chose a blue-green eggshell line of Putian black ducks, a white eggshell line of Putian black ducks and a white eggshell line of Putian white ducks as a research model to study genetic differences in eggshell colour formation. RNA extracted from the shell gland was analysed using RNA-seq to identify candidate genes that participate in eggshell colouration. These results will further elucidate the transcriptional mechanism of blue-green eggshell colour generation and provide N, unknown base rates higher than 10%; Q20, sequencing error rates lower than 1% a foundation for further studies of the molecular basis of eggshell colour pigmentation in avian species.

Overview of RNA-seq results
To identify blue-green eggshell colour-related genes in Putian black ducks, we conducted a transcriptome analysis of shell glands in HQR, HBR and BBR. The main results of RNA-seq are shown in Tables 1 and 2 and Fig. 1. The number of clean reads generated from each library ranged from 64,378,314 to 84,938,126. The Q20 value was more than 97% (Table 1). After removing low-quality reads, adaptor sequences and rRNA reads, 68.35-72.47% of clean reads were mapped uniquely to the Anas platyrhynchos genome, and only a small proportion of them were mapped to multiple locations in the genome (  1).

DEG functional annotation
The DEGs from HQR, HBR and BBR were annotated using the Gene Ontology (GO) database, and the genes were classified into three categories. In biological process, biological regulation, cellular process and single-organism process were the most frequent terms in both comparisons; in cellular component, cell and cell part were the top terms in both comparisons; binding was observed to occur most frequently in molecular function in both comparisons (Fig. 2). Pathway enrichment analysis was carried out using the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database. As shown in Fig. 3, cell adhesion molecules (CAMs), neuroactive ligand-receptor interaction,  PPAR signalling pathway, glycosaminoglycan biosynthesisheparan sulfate/heparin, leukocyte transendothelial migration, hepatitis C and phagosome were significantly enriched in the HQR-vs-HBR comparison (P < 0.05). A study has found that arachidonic acid could induce oviposition and pigment secretion from the shell gland in Japanese quail [27], and glycine is the synthetic basis for protoporphyrin IX and biliverdin. In the HQR-vs-BBR comparison, arachidonic acid metabolism and glycine, serine and threonine metabolism were significantly enriched (P < 0.05), which suggested that they are related to eggshell pigmentation.

Identification of blue-green eggshell colour-associated candidate genes
Biliverdin is the main pigment responsible for blue-green eggshell colour. Here, a total of 4 delta-ALA synthesis-related genes, 4 protoporphyrin IX synthesis-related genes, 3 biliverdin synthesis-related genes, 4 haem breakdown-related genes, 4 biliverdin breakdown-related genes and 4 genes not directly involved in biliverdin metabolism but participating in porphyrin and chlorophyll metabolism were identified among the candidate genes in Fig. 4. ALAS1, EPRS glutamyl-prolyl-tRNA synthetase, aminolevulinate dehydratase (ALAD) and glutamyl-tRNA synthetase 2 (EARS2) participate in the synthesis of ALA [18,28]. Meanwhile, ALAS1 is also the rate-limiting enzyme in ALA synthesis, and this gene was upregulated in the HBR group [18] and downregulated in the HQR and BBR groups. EPRS glutamyl-prolyl-tRNA synthetase and ALAD were upregulated in BBR and downregulated in HQR and HBR. EARS2 was upregulated in the HQR and HBR groups but downregulated in the BBR group. Hydroxymethylbilane synthase (HMBS), uroporphyrinogen III synthase (UROS), uroporphyrinogen decarboxylase (UROD) and coproporphyrinogen oxidase (CPOX) participate in the conversion of ALA into protoporphyrin IX [21,29,30]. HMBS was upregulated in BBR and downregulated in HQR and BBR. UROS was upregulated in HQR and HBR but downregulated in BBR, UROD was upregulated in HQR and BBR but downregulated in HBR, and CPOX was upregulated in HQR but downregulated in HBR and BBR. FECH, a rate-limiting enzyme that is involved in the synthesis of protohaem [31], was upregulated in the HQR and HBR groups but downregulated in the BBR group. Haem oxygenase (HMOX) is the rate-limiting enzyme that participates in the transformation of protohaem into biliverdin [24] and was upregulated in HQR and HBR but downregulated in BBR. Protohaem IX farnesyltransferase and cytochrome c oxidase assembly protein COX15 homologue (COX15) catalyse protohaem into haem A [32]. Protohaem IX farnesyltransferase and COX15 were downregulated in the HQR and HBR groups but upregulated in the BBR group. Protohaem IX farnesyltransferase 2 was upregulated in HQR and downregulated in HBR and BBR. Cytochrome c-type haem lyase participates in the transformation of protohaem into cytochrome c [33], and it was upregulated in the HBR group and downregulated in the HQR and BBR groups. The suppression of the enzymes involved in biliverdin synthesis may be the main reason that Putian white ducks lay white eggshells. In Putian black ducks, the expression levels of biliverdin synthesis-related enzymes show similar patterns in both blue-shelled and white-shelled ducks. In relation to biliverdin breakdown-related genes, BLVRA, UDP-glucuronosyltransferase (UGT) and GUSB glucuronidase beta catalyse biliverdin into D-urobilin [34,35]. BLVRA was upregulated in the HQR and HBR groups but downregulated in the BBR group. UGT1-1-like was upregulated in HQR and downregulated in HBR and BBR, while RichFactor is the ratio of the differentially expressed number of genes in the pathway and the total number of genes in the pathway. The higher the RichFactor, the higher the degree of enrichment. QValue is the P-value after the multiple hypothesis test correction, in the range of 0 to 1; the closer the QValue is to zero, the more significant the enrichment GUSB glucuronidase beta was downregulated in HQR and upregulated in HBR and BBR. UGT2A2 was downregulated in the HQR and HBR groups but upregulated in the BBR group.

Verification of DEGs by qPCR
To confirm the accuracy and reproducibility of the transcriptome analysis results, seven genes were selected for RT-qPCR validation. Primers and candidate genes are shown in Table 3. The expression profiles of seven candidate genes were determined using RT-qPCR, which were consistent with the RNA-seq results (Fig. 5), thus confirming our transcriptome analysis.

Transcript levels of porphyrin and chlorophyll metabolism-related genes by qPCR
In the HQR-vs-BBR comparison, UGT2A2 and UGT1-1-like, which participate in biliverdin breakdown, were two DEGs. Here, 23 porphyrin and chlorophyll metabolism-related genes were analysed by RT-qPCR (Fig. 6). None of these genes were differentially expressed in the RNA-seq results except for UGT2A2 and UGT1-1-like, while some showed significant differences in the RT-qPCR results. Compared with the RT-qPCR results, the expression profiles of 10 genes were not consistent with the RNA-seq data. In the RT-qPCR results, ALAS1 and EPRS glutamyl-prolyl-tRNA synthetase were significantly upregulated in the HBR group compared with the HQR and BBR groups (P < 0.05). HMOX1 was significantly downregulated in BBR compared with HQR and HBR (P < 0.05). BLVRA, GUSB, cytochrome c-type haem lyase, protohaem IX farnesyltransferase and UGT2A2 were significantly upregulated in the HBR and BBR groups compared with the HQR group (P < 0.05). In HBR, 5 biliverdin breakdown-related genes were upregulated. In BBR, 1 biliverdin synthesis-related gene was downregulated, while 5 biliverdin breakdown-related genes were downregulated. The expression patterns of these genes may be the main cause of white eggshell colour.

Discussion
To help elucidate the molecular basis of eggshell pigmentation, we obtained a global view of transcriptomes from the shell glands of the HQR, HBR and HQR groups. However, none of these DEGs were directly involved in the biosynthesis or transportation of biliverdin in the HBR and HQR groups, while UGT2A2 and UGT1-1-like, which participate in biliverdin breakdown, were two DEGs in the HQR-vs-BBR comparison. In HQR-vs-HBR, CAMs, neuroactive ligand-receptor interaction, the PPAR signalling pathway, glycosaminoglycan biosynthesis-heparan sulfate/heparin, leukocyte transendothelial migration, hepatitis C and phagosome were significantly enriched. None of these pathways or DEGs were related to the synthesis of biliverdin. In the HQR-vs-BBR comparison, arachidonic acid metabolism and glycine, serine and threonine metabolism were significantly enriched (P < 0.05). Arachidonic acid participates in the local synthesis of prostaglandin. A previous study has found that both prostaglandin and arachidonic acid can induce pigment secretion, which is consistent with our results [27]. Next, we analysed the expression level of biliverdinrelated genes by RT-qPCR. ALAS1 and EPRS glutamyl-prolyl-tRNA synthetase participate in the biosynthesis of ALA, which is utilized for the biosynthesis of protoporphyrin IX and biliverdin. It is well established that HMOX is the rate-limiting enzyme in the catabolism of haem into biliverdin, free iron, and carbon monoxide [36]. A study has also found that Dongxiang blue-shelled chickens have higher HMOX1 expression levels and enzymatic activity than do Dongxiang brown-shelled chickens [37]. HMOX1 was significantly downregulated in BBR compared to HQR and HBR,  which is consistent with the findings. However, HMOX1 and HMOX2 were not significantly different between the HQR and HBR groups. This result is consistent with the finding of Wang et al. [37] that HMOX1 and HMOX2 were not significantly different between the blue-green-shelled and white-shelled duck groups. BLVRA, GUSB, cytochrome c-type haem lyase, protohaem IX farnesyltransferase and UGT2A2 were significantly upregulated in HBR and BBR compared with HQR. Thus, we speculate that HQR and HBR can produce biliverdin, but HBR cannot accumulate it. This result is also similar to the findings of Japanese scientists, who found that both red coral and blue coral can produce biliverdin, but red coral cannot accumulate biliverdin in its skeleton [38]. Compared with the HQR group, the BBR group produced less biliverdin and did not accumulate it (Fig. 7).

Conclusions
This study conducted a comparative transcriptome analysis of the shell glands of Putian white ducks and Putian black ducks. In the HQR-vs-HBR comparison, none of the DEGs were directly involved in the eggshell pigmentation process. In the HQR-vs-BBR comparison, UGT2A2 and UGT1-1-like were considered to be two DEGs related to biliverdin breakdown. In the KEGG analysis, glycine, serine and threonine metabolism and arachidonic acid metabolism may play potential roles in the synthesis and secretion of pigments, based on the HQR-vs-BBR comparison. From the RT-qPCR results, 5 biliverdin breakdown-related genes were upregulated in the HBR group. In the BBR group, 1 biliverdin synthesis-related gene was downregulated, and 5 biliverdin breakdown-related genes were downregulated. Thus, we can speculate that HQR and HBR can produce biliverdin, but HBR cannot accumulate it. Compared with HQR, BBR produces less biliverdin and cannot accumulate it. These results will provide new insight into further functional genomic research on the mechanism of blue-green eggshell colouration.

Ethics statement
All animal experiments were reviewed and approved by the Institutional Animal Care and Use Committee of the College of Animal Science, Fujian Agriculture and Forestry University. All of the following procedures were performed strictly according to the regulations and guidelines established by the committee.  Table 3, and GAPDH was used as the reference gene. The error bars on each column indicate the SD based on three replicates. Different lowercase letters above the bars indicate statistically significant differences at P < 0.05 (one-way ANOVA, Duncan's tests). Different letters indicate significant differences between groups, while the same letter represents no significant difference. The genes and gene_ids are as follows: carbonic anhydrase 4 (CA4, ncbi_101789643); prostatic acid phosphatase-like (ncbi_101805018); carbonic anhydrase 2 (ncbi_101802744); ovostatin (ncbi_101801176); CUNH2orf40 (ncbi_101802643); phospholipase A2 group XIIA (PLA2G12A, ncbi_101795972); polymeric immunoglobulin receptor-like (PIGR, ncbi_101802412)

Animals and tissue collection
In this study, the Putian black ducks and Putian white ducks were all bred by Shishi Conservation and Research Centre of Waterfowl Genetic Resources. We randomly collected nine Putian white ducks that produced white eggshells, nine Putian black ducks that produced white eggshells and nine Putian black ducks that produced blue-green eggshells. Three shell gland samples of the same phenotype were mixed in equal amounts for further mRNA sequencing. According to their individual oviposition routines, ducks were stunned and killed 3-5 h before estimated oviposition. Then, the shell gland samples were washed with 0.9% isotonic saline and preserved in liquid nitrogen for further RNA extraction. All the ducks were slaughtered at the same time, and their intact shell glands were harvested.

Sequence read quality control, mapping, and annotation
Quality of the raw reads was assessed using FastQC [39]. Raw reads that contained adapter, more than 10% unknown bases or low-quality reads were removed to obtain high-quality clean reads. Next, rRNA was removed from the high-quality clean reads using Bowtie by mapping to a ribosome database [40]. The remaining high-quality clean reads were aligned to the reference genome (Anas platyrhynchos) using TopHat2 [41]. Subsequently, the Cufflinks reference annotation based transcript (RABT) assembly method was used to assemble these mapped reads into possible transcripts and then generate a final transcriptome assembly [42,43]. The GO and KEGG enrichment analyses were performed with a Q-value cut-off of 0.05.

RNA-seq analysis
In this study, the read numbers mapped to each gene were counted using edgeR. The FPKM (fragments per kilobase of transcript per million mapped reads) value of each gene was calculated based on the length of the gene and the read counts mapped to this gene. The resulting P-values were adjusted using the Benjamini and Hochberg approach for controlling the false discovery rate (FDR) [44]. Genes with an adjusted P-value of < 0.05 and | Log2 (fold change) | > 1 were identified as differentially expressed genes (DEGs) by edgeR.

Real-time PCR assay
cDNA synthesis was performed using the GoScript™ Reverse Transcription System (Promega, USA) following the manufacturer's instructions. Quantitative gene expression was assessed by real-time PCR with a CFX96 Touch Deep Well Real-Time PCR Detection System (Bio-Rad, USA). Target genes were amplified using the primers in Table 3: GAPDH was employed as an endogenous control. The reaction mixture contained 6.25 μL 2× GoTaq® qPCR Master Mix (Promega, USA) following the manufacturer's instructions. Each reaction was run in triplicate. The qRT-PCR conditions were an initial predenaturation step at 95°C for 3 min, followed by 39 cycles of 95°C for 20 s and 60°C for 30 s. Fluorescence data were collected and analysed with the Bio-Rad CFX Manager software and normalized to those of GAPDH. Melting curves were generated after the qPCR was completed.

Real-time PCR analysis
Differences between groups were analysed using ANOVA. Experimental data are expressed as the mean ± SE from three independent experiments. P-values less than 0.05 were considered statistically significant. Statistical analyses were performed by using SPSS version 10.0. Fig. 7 A hypothesis explaining that differential expression of biliverdin synthesis-related genes can cause blue-green eggshell colour. Gene names shown in red represent rate-limiting enzymes. The blue upward arrows represent upregulation in HBR compared to HQR. The blue downward arrow represents downregulation in HBR compared to HQR. The yellow upward arrows represent upregulation in BBR compared to HQR. The yellow downward arrow represents downregulation in HBR compared to HQR